Rajasimhan Rajagovindan
2012-Mar-28 18:23 UTC
[R] discrepancy between paired t test and glht on lme models
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 differences1.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]]
peter dalgaard
2012-Mar-29 07:36 UTC
[R] discrepancy between paired t test and glht on lme models
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