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

[Torsten Hothorn]

> Message: 64
> Date: Wed, 14 May 2003 11:39:14 +0200
> From: "Dieter Menne" <dieter.menne at menne-biomed.de>
> Subject: [R] Multiple comparison and lme (again, sorry)
> To: "R-Help" <r-help at stat.math.ethz.ch>
> Message-ID:
> 	<JLEPLGAANFCEAEDCAGJNEEBNCFAA.dieter.menne at menne-biomed.de>
> Content-Type: text/plain;   charset="iso-8859-1"
>
> 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))

`asympt = TRUE' must be there for a logistic regression model and
it works with mvtnorm_0.5-12 (on CRAN).

>
> -------
> 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")))

you need to specify the `how.many' argument:

 ## need to specify how.many
R> contrasts(treat,how.many = 3) <- zapsmall(
                                    mginv(contrMat(table(treat),
                                                   type="Tukey")))

R> attr(treat, "contrasts")
         T1-P       T2-P      T2-T1
P  -0.3333333 -0.3333333  0.0000000
T1  0.3333333  0.0000000 -0.3333333
T2  0.0000000  0.3333333  0.3333333

 ## however, this will not work!
R> example.lme<-lme(response~treat, random=~1|subj)
Error in MEEM(object, conLin, control$niterEM) :
        Singularity in backsolve at level 0, block 1

The problem is a basic one: AFAIK, none of the `model.fit' functions deals
with singular design matrices without additional work (have a look at
TukeyHSD for what that means) as induced for example by Tukey contrasts
(I would be happy to be wrong here).

And that is the reason why we can't use `lm.fit' for parameter estimation
for arbitrary contrasts within `multcomp'. The `multcomp' package, as a
workaround, estimates the parameters of a linear model (!) via the
Moore-Penrose inverse.

For all other models, we currently only offer a low-level interface
(`csim{int,test}') which assumes that the
parameters / contrasts and their covariance matrix are there. As long as
those values can be computed via `coef' and `vcov', everything is fine
(and we are currently working on methods for `lm' and `glm' to the
`simint' and `simtest' generics). If not, as it is the case here, there
is no general solution (and I can't answer this particular question).


Best,

Torsten

> example<-data.frame(response, treat, subj)
> example.lme<-lme(response~treat, data=example, random=~1|subj)
> summary(example.lme)
>
> #summary(csimtest.....)
>
>
> ------------------------------
>
> _______________________________________________
> R-help at stat.math.ethz.ch mailing list  DIGESTED
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
>
> End of R-help Digest, Vol 3, Issue 14
> *************************************
>
>