[R] nlme and multiple comparisons
Pedro J. Aphalo
pedro.aphalo at cc.jyu.fi
Fri Feb 27 13:07:54 CET 2004
I received some suggestions about how to do multiple comparisons,
between levels of a factor used as explanatory variable in the
fixed part of a model in an nlme fit. Many thanks!
I also received a request to summarize. So, here is a summary of my
attempts at following the suggestions, and further questions...
FIRST SUGGESTION
Use the L argument to anova to compute the contrasts.
Then use p.adjust.
Works just fine with nlme objects.
(And results are consistent when changing the
order of levels in the factor used as explanatory variable.)
SECOND SUGGESTION
Use fit.contrast from package gregmisc.
> library(nlme)
> library(gregmisc)
> fit.contrast(fm5bPrun.nlme, "factor(soil.temp)", c(0,-1,1))
Error in nlme.formula(model = log.area ~ fPrunty(day, Asym, Slope,
Curve, :
unused argument(s) (contrasts ...)
>
nlme does not have a contrasts argument as lme has.
fit.contrast has a method for lme, not for nlme.
Does not work with nlme objects.
THIRD SUGGESTION (actually pointing to an old thread from May 2003)
http://maths.newcastle.edu.au/~rking/R/help/03a/5046.html
Based on the example given by Torsten Hothorn for glm,
and suggested to be adaptable to lme. I tried the following
example, based on an example in chapter 6 of Pinheiro and
Bates (2000): (please note that the model used here is not the
final model described in the book and can be improved, but
it is good enough, I think, for this example.)
##begin R code
library(nlme)
library(multcomp)
data(Soybean)
fm1Soy.lis <- nlsList( weight ~ SSlogis(Time, Asym, xmid, scal),
data = Soybean )
#fm1Soy.lis
fm1Soy.nlme <- nlme( fm1Soy.lis )
#fm1Soy.nlme
fm2Soy.nlme <- update( fm1Soy.nlme, weights = varPower() )
soyFix <- fixef( fm2Soy.nlme )
options( contrasts = c("contr.treatment", "contr.poly") )
fm3Soy.nlme <- update( fm2Soy.nlme,
fixed = Asym + xmid + scal ~ Year,
start = c(soyFix[1], 0, 0, soyFix[2], 0, 0, soyFix[3], 0, 0) )
Soybean$YearTuk <- Soybean$Year
contrasts(Soybean$YearTuk) <-
mginv(contrMat(table(Soybean$YearTuk),type="Tukey"))
fm3SoyTuk.nlme <- update( fm2Soy.nlme,
fixed = Asym + xmid + scal ~ YearTuk,
start = c(soyFix[1], 0, 0, soyFix[2], 0, 0, soyFix[3], 0, 0) )
# and the following for the first parameter Asym:
summary(csimtest(fixef(fm3SoyTuk.nlme)[2:3],fm3SoyTuk.nlme$varFix[2:3,2:3],cmatrix=diag(2),df=477))
#This does not give all the pairwise contrasts, so I
#rearange the order of the levels in the factor, and redo.
Soybean$YearTukRev <- factor(Soybean$Year, levels=c(1990,1989,1988))
contrasts(Soybean$YearTukRev) <-
mginv(contrMat(table(Soybean$YearTukRev),type="Tukey"))
fm3SoyTukRev.nlme <- update( fm2Soy.nlme,
fixed = Asym + xmid + scal ~ YearTukRev,
start = c(soyFix[1], 0, 0, soyFix[2], 0, 0, soyFix[3], 0, 0) )
#and the following for the first parameter Asym:
summary(csimtest(fixef(fm3SoyTukRev.nlme)[2:3],fm3SoyTukRev.nlme$varFix[2:3,2:3],cmatrix=diag(2),df=477))
##end R
But now we have two pairs of simultaneous tests, and the test that
is common gives a totally different P-value for the 1988-1990, versus
1990-1988 comparisons. This difference can be clearly seen
with summary directly on the fitted model objects, where it is easier
to find the coefficient estimates for the tests than in the output of
csimtest. Also note that the estimate of the intercept does not change
when reordering the levels, so what does this represent? And the second
coefficient has not changed either, although at least the name of the
coefficient indicates that it represents a different comparison...
> summary(fm3SoyTuk.nlme)
(edited)
[...]
Fixed effects: Asym + xmid + scal ~ YearTuk
Value Std.Error DF t-value p-value
Asym.(Intercept) 16.95143 0.5407434 356 31.34839 0.0000
Asym.YearTuk1989-1988 -9.14093 2.1797271 356 -4.19361 0.0000
Asym.YearTuk1990-1988 -0.62722 2.2891767 356 -0.27400 0.7842
xmid.(Intercept) 51.65615 0.3881693 356 133.07636 0.0000
xmid.YearTuk1989-1988 -0.11277 1.6095886 356 -0.07006 0.9442
xmid.YearTuk1990-1988 -7.21578 1.6281871 356 -4.43179 0.0000
scal.(Intercept) 7.52005 0.0844056 356 89.09412 0.0000
scal.YearTuk1989-1988 -1.20229 0.3470370 356 -3.46444 0.0006
scal.YearTuk1990-1988 -0.39143 0.3667542 356 -1.06727 0.2866
[...]
> summary(fm3SoyTukRev.nlme)
(edited)
[...]
Fixed effects: Asym + xmid + scal ~ YearTukRev
Value Std.Error DF t-value p-value
Asym.(Intercept) 16.95152 0.5407411 356 31.34867 0.0000
Asym.YearTukRev1989-1990 -9.14090 2.1797162 356 -4.19362 0.0000
Asym.YearTukRev1988-1990 9.76819 2.4079472 356 4.05665 0.0001
xmid.(Intercept) 51.65625 0.3881701 356 133.07634 0.0000
xmid.YearTukRev1989-1990 -0.11268 1.6095926 356 -0.07001 0.9442
xmid.YearTukRev1988-1990 7.32858 1.7013868 356 4.30742 0.0000
scal.(Intercept) 7.52006 0.0844060 356 89.09388 0.0000
scal.YearTukRev1989-1990 -1.20228 0.3470388 356 -3.46439 0.0006
scal.YearTukRev1988-1990 1.59373 0.3602370 356 4.42411 0.0000
[...]
So I think that I am doing something wrong with the contrasts, but
I do not know what...
So I ask again for help. (I think I can follow the first suggestion
above with my own data, but I am curious about what is going on here...)
Once again, many thanks in advance.
Pedro.
Pedro J. Aphalo wrote:
> This is only partly a question about R, as I am not quite sure about the
> underlying statistical theory either.
>
> I have fitted a non-linear mixed-effects model with nlme. In the fixed
> part of the model I have a factor with three levels as explanatory
> variable. I would like to use Tukey HSD or a similar test to test for
> differences between these three levels.
>
> I have two grouping factors: 'plant' to which the treatments were
> assigned at random, and 'leaf' which are subsamples.
>
> With summary, with the default setting for contrasts I get two of the
> possible three comparisons. One possibility that I can think of is to
> change the order of the levels in the factor, repeat the fit, use
> summary again, and finally use p.adjust. Would this be valid?
>
> Is there a more elegant solution?
>
> If it matters, the data are slightly unbalanced (missing observations).
>
> Sorry that this message became so long... Thanks in advance for any
> suggestions, and I hope I haven't missed the answer when looking at the
> help pages, FAQ, MASS(3ed) and Pinheiro and Bates' book.
>
I did miss the L argument to anova...
> Pedro.
>
--
==================================================================
Pedro J. Aphalo
Department of Biological and Environmental Science
University of Jyväskylä
P.O. Box 35, 40351 JYVÄSKYLÄ, Finland
Phone +358 14 260 2339
Mobile +358 50 3721504
Fax +358 14 260 2321
mailto:pedro.aphalo at cc.jyu.fi
http://www.jyu.fi/~aphalo/ ,,,^..^,,,
More information about the R-help
mailing list