[Rd] (PR#8905) Recommended package nlme: bug in predict.lme when an independent variable is a polynomial

[Renaud Lancelot]

Many thanks for your very useful comments and suggestions.

Renaud

2006/5/30, Prof Brian Ripley <ripley at stats.ox.ac.uk>:
> On Tue, 30 May 2006, Prof Brian Ripley wrote:
>
> > This is not really a bug.  See
> >
> > http://developer.r-project.org/model-fitting-functions.txt
> >
> > for how this is handled in other packages. All model-fitting in R used to
> > do this (and it is described in the White Book and MASS1-3).
> >
> > predict.lme does not use model.frame as described in that URL.  Dr Bates'
> > recent response to another query applies here: lmer is more standard and I
> > suggest you try it instead.   (I don't think anyone is going to be
> > rewriting lme to use model.frame: it is essentially in maintainence mode.)
>
> Another workaround is to use poly(..., raw=TRUE).  I don't actually see
> anything in this report using predict.lme, but compare
>
> > predict(fm, Newdata)
>       M01      M01      M01      M01      M02      M02
> 21.01353 25.63852 32.30504 38.59640 17.24007 21.86507
> attr(,"label")
> [1] "Predicted values (mm)"
> > fm <- lme(distance ~ poly(age, 3, raw=TRUE) + Sex, data = Orthodont,
>              random = ~1)
> > predict(fm, Newdata)
>       M01      M01      M01      M01      M02      M02
> 25.52963 26.51111 27.99259 29.43703 21.75617 22.73765
> attr(,"label")
> [1] "Predicted values (mm)"
>
>
> > On Sat, 27 May 2006, renaud.lancelot at cirad.fr wrote:
> >
> >> Full_Name: Renaud Lancelot
> >> Version: Version 2.3.0 (2006-04-24)
> >> OS: MS Windows XP Pro SP2
> >> Submission from: (NULL) (82.239.219.108)
> >>
> >>
> >> I think there is a bug in predict.lme, when a polynomial generated by poly() is
> >> used as an explanatory variable, and a new data.frame is used for predictions. I
> >> guess this is related to * not * using, for predictions, the coefs used in
> >> constructing the orthogonal polynomials before fitting the model:
> >>
> >>> fm <- lme(distance ~ poly(age, 3) + Sex, data = Orthodont, random = ~ 1)
> >>>
> >>> # data for predictions
> >>> Newdata <- head(Orthodont)
> >>> Newdata$Sex <- factor(Newdata$Sex, levels = levels(Orthodont$Sex))
> >>>
> >>> # "naive" model matrix for predictions
> >>> mm1 <- model.matrix(~ poly(age, 3) + Sex, data = Newdata)
> >>>
> >>> # "correct" model matrix for predictions
> >>> p <- poly(Orthodont$age, 3)
> >>> mm2 <- model.matrix(~ poly(age, 3, coefs = attr(p, "coefs")) + Sex, data =
> >> Newdata)
> >>>
> >>> data.frame(pred1 = predict(fm, level = 0, newdata = Newdata),
> >> +            pred2 = mm1 %*% fixef(fm),
> >> +            pred3 = head(predict(fm, level = 0)),
> >> +            pred4 = mm2 %*% fixef(fm))
> >>     pred1    pred2    pred3    pred4
> >> 1 18.61469 18.61469 23.13079 23.13079
> >> 2 23.23968 23.23968 24.11227 24.11227
> >> 3 29.90620 29.90620 25.59375 25.59375
> >> 4 36.19756 36.19756 27.03819 27.03819
> >> 5 18.61469 18.61469 23.13079 23.13079
> >> 6 23.23968 23.23968 24.11227 24.11227
> >>
> >> Best regards,
> >>
> >> Renaud
> >>
> >> ______________________________________________
> >> R-devel at r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-devel
> >>
> >>
> >
> >
>
> --
> Brian D. Ripley,                  ripley at stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>


-- 
Renaud LANCELOT
Département Elevage et Médecine Vétérinaire (EMVT) du CIRAD
Directeur adjoint chargé des affaires scientifiques

CIRAD, Animal Production and Veterinary Medicine Department
Deputy director for scientific affairs

Campus international de Baillarguet
TA 30 / B (Bât. B, Bur. 214)
34398 Montpellier Cedex 5 - France
Tél   +33 (0)4 67 59 37 17
Secr. +33 (0)4 67 59 39 04
Fax   +33 (0)4 67 59 37 95