[R] discrepancy between paired t test and glht on lme models
Greg Snow
538280 at gmail.com
Fri Mar 30 16:41:22 CEST 2012
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 <pdalgd at gmail.com> 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]]
>>
>> ______________________________________________
>> 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
>
> ______________________________________________
> 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.
--
Gregory (Greg) L. Snow Ph.D.
538280 at gmail.com
More information about the R-help
mailing list