discrepancy between paired t test and glht on lme models

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

discrepancy between paired t test and glht on lme models

Rajasimhan Rajagovindan
Hi folks,



I am working with repeated measures data and I ran into issues where the
paired t-test results did not match those obtained by employing glht()
contrasts on a lme model. While the lme model itself appears to be fine,
there seems to be some discrepancy with using glht() on the lme model
(unless I am missing something here).  I was wondering if someone could
help identify the issue. On my actual dataset the  differences between
glht() and paired t test is more severe than the example provided here.



I am using glht() for my data since I need to perform pairwise comparisons
across multiple levels, any alternate approach to performing posthoc
comparisons on lme object is also welcome.

I have included the code and the results from a mocked up data (one that I
found online) here.





require(nlme)

require(multcomp)



dv <- c(1,3,2,2,2,5,3,4,3,5)

subject <- factor(c("s1","s1","s2","s2","s3","s3","s4","s4","s5","s5"))

myfactor <- factor(c("f1","f2","f1","f2","f1","f2","f1","f2","f1","f2"))

mydata <- data.frame(dv, subject, myfactor)

rm(subject,myfactor,dv)

attach(mydata)



# paired t test (H0: f2-f1 = 0)

t.test(mydata[myfactor=='f2',1],mydata[myfactor=='f1',1],paired=TRUE)

# yields :  t = 3.1379, df = 4, p-value = 0.03492, mean of the differences=
1.6





# lme (f1 as reference level)

fit.lme <- lme(dv ~ myfactor, random =
~1|subject,method="REML",correlation=corCompSymm(),data=mydata)

summary(fit.lme) # yields identical results as paired t test

# f2-f1:   t = 3.1379, df = 4, p-value = 0.0349



summary(glht(fit.lme,linfct=mcp(myfactor="Tukey")))

# while test statistic is comparable, p value is different

# have noticed cases where the differences between glht() and paired t test
is more severe



########### sample outputs from the script ###############


######### things appear ok here and match paired t test results
#############

summary(fit.lme)

Linear mixed-effects model fit by REML

 Data: mydata

       AIC      BIC    logLik

  36.43722 36.83443 -13.21861



Random effects:

 Formula: ~1 | subject

        (Intercept)  Residual

StdDev:   0.7420274 0.8058504



Correlation Structure: Compound symmetry

 Formula: ~1 | subject

 Parameter estimate(s):

          Rho

-0.0009325763

Fixed effects: dv ~ myfactor

            Value Std.Error DF  t-value p-value

(Intercept)   2.2 0.4898979  4 4.490732  0.0109

myfactorf2    1.6 0.5099022  4 3.137857  0.0349

 Correlation:

           (Intr)

myfactorf2 -0.52



Standardized Within-Group Residuals:

        Min          Q1         Med          Q3         Max

-1.45279696 -0.53193228  0.03481143  0.58490026  1.09867599



Number of Observations: 10

Number of Groups: 5



######### result differs from paired t test  !!!!!

 summary(glht(fit.lme,linfct=mcp(myfactor="Tukey")),test=adjusted("none"))



         Simultaneous Tests for General Linear Hypotheses



Multiple Comparisons of Means: Tukey Contrasts





Fit: lme.formula(fixed = dv ~ myfactor, data = mydata, random = ~1 |

    subject, correlation = corCompSymm(), method = "REML")



Linear Hypotheses:

             Estimate Std. Error z value Pr(>|z|)

f2 - f1 == 0   1.6000     0.5099   3.138   0.0017 **     <<<<------

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- none method)



####################

platform       i386-pc-mingw32

arch           i386

os             mingw32

system         i386, mingw32

status

major          2

minor          13.1

year           2011

month          07

day            08

svn rev        56322

language       R

version.string R version 2.13.1 (2011-07-08)

        [[alternative HTML version deleted]]


______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: discrepancy between paired t test and glht on lme models

Peter Dalgaard-2

On Mar 28, 2012, at 20:23 , Rajasimhan Rajagovindan wrote:

> Hi folks,
>
>
>
> I am working with repeated measures data and I ran into issues where the
> paired t-test results did not match those obtained by employing glht()
> contrasts on a lme model. While the lme model itself appears to be fine,
> there seems to be some discrepancy with using glht() on the lme model
> (unless I am missing something here).  I was wondering if someone could
> help identify the issue. On my actual dataset the  differences between
> glht() and paired t test is more severe than the example provided here.


You might want to move to the R-sig-ME (mixed effects) mailing list for up to date advice.

The basic issue appears to be that glht is not smart enough to deal with degrees of freedom so it uses an asymptotic z-test instead of a t-test. Infinite df, basically, and since 4 is a pretty poor approximation of infinity, you get your discrepancy.

It's not that surprising, given that lme() itself is pretty poor at figuring out df in some cases. Especially if you have to deal with cross-stratum effects, the calculation of appropriate degrees of freedom is nontrivial. Some recent developments allow the calculation of Kenward-Roger for the lmer() models, but I wouldn't know to what extend this carries to glht-style testing.

>
>
>
> I am using glht() for my data since I need to perform pairwise comparisons
> across multiple levels, any alternate approach to performing posthoc
> comparisons on lme object is also welcome.
>
> I have included the code and the results from a mocked up data (one that I
> found online) here.
>
>
>
>
>
> require(nlme)
>
> require(multcomp)
>
>
>
> dv <- c(1,3,2,2,2,5,3,4,3,5)
>
> subject <- factor(c("s1","s1","s2","s2","s3","s3","s4","s4","s5","s5"))
>
> myfactor <- factor(c("f1","f2","f1","f2","f1","f2","f1","f2","f1","f2"))
>
> mydata <- data.frame(dv, subject, myfactor)
>
> rm(subject,myfactor,dv)
>
> attach(mydata)
>
>
>
> # paired t test (H0: f2-f1 = 0)
>
> t.test(mydata[myfactor=='f2',1],mydata[myfactor=='f1',1],paired=TRUE)
>
> # yields :  t = 3.1379, df = 4, p-value = 0.03492, mean of the differences=
> 1.6
>
>
>
>
>
> # lme (f1 as reference level)
>
> fit.lme <- lme(dv ~ myfactor, random =
> ~1|subject,method="REML",correlation=corCompSymm(),data=mydata)
>
> summary(fit.lme) # yields identical results as paired t test
>
> # f2-f1:   t = 3.1379, df = 4, p-value = 0.0349
>
>
>
> summary(glht(fit.lme,linfct=mcp(myfactor="Tukey")))
>
> # while test statistic is comparable, p value is different
>
> # have noticed cases where the differences between glht() and paired t test
> is more severe
>
>
>
> ########### sample outputs from the script ###############
>
>
> ######### things appear ok here and match paired t test results
> #############
>
> summary(fit.lme)
>
> Linear mixed-effects model fit by REML
>
> Data: mydata
>
>       AIC      BIC    logLik
>
>  36.43722 36.83443 -13.21861
>
>
>
> Random effects:
>
> Formula: ~1 | subject
>
>        (Intercept)  Residual
>
> StdDev:   0.7420274 0.8058504
>
>
>
> Correlation Structure: Compound symmetry
>
> Formula: ~1 | subject
>
> Parameter estimate(s):
>
>          Rho
>
> -0.0009325763
>
> Fixed effects: dv ~ myfactor
>
>            Value Std.Error DF  t-value p-value
>
> (Intercept)   2.2 0.4898979  4 4.490732  0.0109
>
> myfactorf2    1.6 0.5099022  4 3.137857  0.0349
>
> Correlation:
>
>           (Intr)
>
> myfactorf2 -0.52
>
>
>
> Standardized Within-Group Residuals:
>
>        Min          Q1         Med          Q3         Max
>
> -1.45279696 -0.53193228  0.03481143  0.58490026  1.09867599
>
>
>
> Number of Observations: 10
>
> Number of Groups: 5
>
>
>
> ######### result differs from paired t test  !!!!!
>
> summary(glht(fit.lme,linfct=mcp(myfactor="Tukey")),test=adjusted("none"))
>
>
>
>         Simultaneous Tests for General Linear Hypotheses
>
>
>
> Multiple Comparisons of Means: Tukey Contrasts
>
>
>
>
>
> Fit: lme.formula(fixed = dv ~ myfactor, data = mydata, random = ~1 |
>
>    subject, correlation = corCompSymm(), method = "REML")
>
>
>
> Linear Hypotheses:
>
>             Estimate Std. Error z value Pr(>|z|)
>
> f2 - f1 == 0   1.6000     0.5099   3.138   0.0017 **     <<<<------
>
> ---
>
> Signif. codes:  0 Œ***‚ 0.001 Œ**‚ 0.01 Œ*‚ 0.05 Œ.‚ 0.1 Œ ‚ 1
>
> (Adjusted p values reported -- none method)
>
>
>
> ####################
>
> platform       i386-pc-mingw32
>
> arch           i386
>
> os             mingw32
>
> system         i386, mingw32
>
> status
>
> major          2
>
> minor          13.1
>
> year           2011
>
> month          07
>
> day            08
>
> svn rev        56322
>
> language       R
>
> version.string R version 2.13.1 (2011-07-08)
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

--
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: discrepancy between paired t test and glht on lme models

glsnow
I nominate the following paragraph for the fortunes package:

"The basic issue appears to be that glht is not smart enough to deal
with degrees of freedom so it uses an asymptotic z-test instead of a
t-test. Infinite df, basically, and since 4 is a pretty poor
approximation of infinity, you get your discrepancy."



On Thu, Mar 29, 2012 at 1:36 AM, peter dalgaard <[hidden email]> wrote:

>
> On Mar 28, 2012, at 20:23 , Rajasimhan Rajagovindan wrote:
>
>> Hi folks,
>>
>>
>>
>> I am working with repeated measures data and I ran into issues where the
>> paired t-test results did not match those obtained by employing glht()
>> contrasts on a lme model. While the lme model itself appears to be fine,
>> there seems to be some discrepancy with using glht() on the lme model
>> (unless I am missing something here).  I was wondering if someone could
>> help identify the issue. On my actual dataset the  differences between
>> glht() and paired t test is more severe than the example provided here.
>
>
> You might want to move to the R-sig-ME (mixed effects) mailing list for up to date advice.
>
> The basic issue appears to be that glht is not smart enough to deal with degrees of freedom so it uses an asymptotic z-test instead of a t-test. Infinite df, basically, and since 4 is a pretty poor approximation of infinity, you get your discrepancy.
>
> It's not that surprising, given that lme() itself is pretty poor at figuring out df in some cases. Especially if you have to deal with cross-stratum effects, the calculation of appropriate degrees of freedom is nontrivial. Some recent developments allow the calculation of Kenward-Roger for the lmer() models, but I wouldn't know to what extend this carries to glht-style testing.
>
>>
>>
>>
>> I am using glht() for my data since I need to perform pairwise comparisons
>> across multiple levels, any alternate approach to performing posthoc
>> comparisons on lme object is also welcome.
>>
>> I have included the code and the results from a mocked up data (one that I
>> found online) here.
>>
>>
>>
>>
>>
>> require(nlme)
>>
>> require(multcomp)
>>
>>
>>
>> dv <- c(1,3,2,2,2,5,3,4,3,5)
>>
>> subject <- factor(c("s1","s1","s2","s2","s3","s3","s4","s4","s5","s5"))
>>
>> myfactor <- factor(c("f1","f2","f1","f2","f1","f2","f1","f2","f1","f2"))
>>
>> mydata <- data.frame(dv, subject, myfactor)
>>
>> rm(subject,myfactor,dv)
>>
>> attach(mydata)
>>
>>
>>
>> # paired t test (H0: f2-f1 = 0)
>>
>> t.test(mydata[myfactor=='f2',1],mydata[myfactor=='f1',1],paired=TRUE)
>>
>> # yields :  t = 3.1379, df = 4, p-value = 0.03492, mean of the differences=
>> 1.6
>>
>>
>>
>>
>>
>> # lme (f1 as reference level)
>>
>> fit.lme <- lme(dv ~ myfactor, random =
>> ~1|subject,method="REML",correlation=corCompSymm(),data=mydata)
>>
>> summary(fit.lme) # yields identical results as paired t test
>>
>> # f2-f1:   t = 3.1379, df = 4, p-value = 0.0349
>>
>>
>>
>> summary(glht(fit.lme,linfct=mcp(myfactor="Tukey")))
>>
>> # while test statistic is comparable, p value is different
>>
>> # have noticed cases where the differences between glht() and paired t test
>> is more severe
>>
>>
>>
>> ########### sample outputs from the script ###############
>>
>>
>> ######### things appear ok here and match paired t test results
>> #############
>>
>> summary(fit.lme)
>>
>> Linear mixed-effects model fit by REML
>>
>> Data: mydata
>>
>>       AIC      BIC    logLik
>>
>>  36.43722 36.83443 -13.21861
>>
>>
>>
>> Random effects:
>>
>> Formula: ~1 | subject
>>
>>        (Intercept)  Residual
>>
>> StdDev:   0.7420274 0.8058504
>>
>>
>>
>> Correlation Structure: Compound symmetry
>>
>> Formula: ~1 | subject
>>
>> Parameter estimate(s):
>>
>>          Rho
>>
>> -0.0009325763
>>
>> Fixed effects: dv ~ myfactor
>>
>>            Value Std.Error DF  t-value p-value
>>
>> (Intercept)   2.2 0.4898979  4 4.490732  0.0109
>>
>> myfactorf2    1.6 0.5099022  4 3.137857  0.0349
>>
>> Correlation:
>>
>>           (Intr)
>>
>> myfactorf2 -0.52
>>
>>
>>
>> Standardized Within-Group Residuals:
>>
>>        Min          Q1         Med          Q3         Max
>>
>> -1.45279696 -0.53193228  0.03481143  0.58490026  1.09867599
>>
>>
>>
>> Number of Observations: 10
>>
>> Number of Groups: 5
>>
>>
>>
>> ######### result differs from paired t test  !!!!!
>>
>> summary(glht(fit.lme,linfct=mcp(myfactor="Tukey")),test=adjusted("none"))
>>
>>
>>
>>         Simultaneous Tests for General Linear Hypotheses
>>
>>
>>
>> Multiple Comparisons of Means: Tukey Contrasts
>>
>>
>>
>>
>>
>> Fit: lme.formula(fixed = dv ~ myfactor, data = mydata, random = ~1 |
>>
>>    subject, correlation = corCompSymm(), method = "REML")
>>
>>
>>
>> Linear Hypotheses:
>>
>>             Estimate Std. Error z value Pr(>|z|)
>>
>> f2 - f1 == 0   1.6000     0.5099   3.138   0.0017 **     <<<<------
>>
>> ---
>>
>> Signif. codes:  0 Œ***‚ 0.001 Œ**‚ 0.01 Œ*‚ 0.05 Œ.‚ 0.1 Œ ‚ 1
>>
>> (Adjusted p values reported -- none method)
>>
>>
>>
>> ####################
>>
>> platform       i386-pc-mingw32
>>
>> arch           i386
>>
>> os             mingw32
>>
>> system         i386, mingw32
>>
>> status
>>
>> major          2
>>
>> minor          13.1
>>
>> year           2011
>>
>> month          07
>>
>> day            08
>>
>> svn rev        56322
>>
>> language       R
>>
>> version.string R version 2.13.1 (2011-07-08)
>>
>>       [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> --
> Peter Dalgaard, Professor
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: [hidden email]  Priv: [hidden email]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



--
Gregory (Greg) L. Snow Ph.D.
[hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.