Hello All. Thanks to Frederic Gosselin, Don MacQueen, and Thomas Hotz for
suggested solutions to this problem, namely that multicomp doesn't have a
specific method for lme models. The answer is to use the functions fixef
and vcov to extracted the fixed effects parameter estimates, and
variance-covariance matrix, from the lme object, and then submit these
quantities to multicomp.default. I did this two different ways, using the
Machines data set that comes with S-Plus and lme. This is a two-way mixed
model design with replication, with Machines the fixed and Worker the random
effect. In the first attempt, I had lme produce estimated means for each
machine by specifying no intercept for the model. See below.
> machmixed.lme <- lme(score ~ Machine - 1, data = Machines, random = ~ 1 |
Worker/Machine)
> x <- fixef(machmixed.lme)
> df <- 10
> vmat <- vcov(machmixed.lme)
> multicomp.default(x, vmat, df.residual = df)
95 % simultaneous confidence intervals for specified
linear combinations, by the Tukey method
critical point: 2.7415
response variable:
intervals excluding 0 are flagged by '****'
Estimate Std.Error Lower Bound Upper Bound
MachineA-MachineB -7.97 2.18 -13.9 -2.0000 ****
MachineA-MachineC -13.90 2.18 -19.9 -7.9500 ****
MachineB-MachineC -5.95 2.18 -11.9 0.0181
The resulting intervals and tests were identical to those produced by JMP,
so we must be on the right track. I also attempted a solution using a model
including an intercept. The model parameters are not the means in this
case, and must be manipulated to generate the means for each machine, and an
associated variance-covariance matrix.
> options(contrasts = c(factor = "contr.treatment", ordered = "contr.poly"))
> machmixed.lme <- lme(score ~ Machine, data = Machines, random = ~ 1 |
Worker/Machine)
> params <- fixef(machmixed.lme)
> nparams <- length(params)
> muest <- cbind(rep(1, nparams), contr.treatment(nparams))
> muvec <- as.numeric(muest %*% params)
> names(muvec) <- c("A", "B", "C")
> vmat <- muest %*% vcov(machmixed.lme) %*% t(muest)
> df <- 10
> multicomp.default(muvec, vmat, df.residual = df)
95 % simultaneous confidence intervals for specified
linear combinations, by the Tukey method
critical point: 2.7415
response variable:
intervals excluding 0 are flagged by '****'
Estimate Std.Error Lower Bound Upper Bound
A-B -7.97 2.18 -13.9 -2.0000 ****
A-C -13.90 2.18 -19.9 -7.9500 ****
B-C -5.95 2.18 -11.9 0.0181
This gives the same answer as before, but I think this approach is more
general. Thanks again for the help.
Best wishes,
John Reeve
************************************
John D. Reeve
Dept. of Zoology
1125 Lincoln Drive
Southern Illinois University
Carbondale, IL 62901-6501
(618) 453-6670 (office)
(618) 453-2806 (fax)
jreeve@zoology.siu.edu
www.science.siu.edu/zoology/reeve/
************************************
|