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]]


