[R] discrepancy between paired t test and glht on lme models

[peter dalgaard]

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]]
> 
> ______________________________________________
> R-help at r-project.org 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: pd.mes at cbs.dk  Priv: PDalgd at gmail.com