[R] Multiple comparison and lme (again, sorry)

[Dieter Menne]

Dear list,

As a reply to my recent mail:

> simint and TukeyHSD work for aov objects.
> Can someone point me to similar functions for lme objects?

Douglas Bates wrote

There aren't multiple comparison methods for lme objects because it is
not clear how to do multiple comparisons for these.  I don't think the
theory of multiple comparisons extends easily to lme models.  One
could use approximations but it is not clear how good the
approximations are.
--------

This should serve as a warning, but let's assume I can live with "only
approximations".

In another thread, Thorsten Hothorn suggested for glm (slightly edited)

library(multcomp)
set.seed(290875)
# a factor at three levels
group <- factor(c(rep(1,10), rep(2, 10), rep(3,10)))
# Williams contrasts
contrasts(group)<-zapsmall(mginv(contrMat(table(group), type="Will")))
# a binary response
z <- factor(rbinom(30, 1, 0.5))
# estimate the model
gmod <- glm( z ~ group, family=binomial(link = "logit"))
summary(gmod)
# exclude the intercept

# Should be the following, but does not work due to a confirmed
# bug in the CRAN-binary version 5.10
#summary(csimtest(coef(gmod)[2:3], vcov(gmod)[2:3,2:3],
#                 cmatrix=diag(2), df=27,asympt=T))
summary(csimtest(coef(gmod)[2:3], vcov(gmod)[2:3,2:3],
                 cmatrix=diag(2), df=27))

-------
This works and can be extended to to lme, but only gives TWO of the three
possible Tukey contrasts. Setting type="Tukey" does not work, as by
assigning to
contrasts(group) only two independent columns are copied (which makes
sense).

# Can someone help me out with the example below: I need P-T1, P-T2, T1-T2

library(nlme)
library(multcomp)
set.seed(701)
sub<-rnorm(50,0,2.0) # subjects random
vr<-0.2
response<-c(rnorm(50,0,vr)+sub,rnorm(50,2,vr)+sub,rnorm(50,3,vr)+sub)
treat<-factor(c(rep("P",50),rep("T1",50),rep("T2",50)))
subj<-factor(rep(1:50,3))
# The following only gives me 2 (instead of 3) contrasts
#contrasts(treat)<-zapsmall(mginv(contrMat(c(50,50,50), type="Tukey")))
example<-data.frame(response, treat, subj)
example.lme<-lme(response~treat, data=example, random=~1|subj)
summary(example.lme)

#summary(csimtest.....)