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]]
______________________________________________
[email protected] 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.
