R devel - May 2006 - (PR#8905) Recommended package nlme: bug in predict.lme when an independent variable is a polynomial

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=3DTRUE).  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=3DTRUE) + Sex, data =3D
Orthodont,
>              random =3D ~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 =3D
Orthodont, random =3D ~ 1)
> >>>
> >>> # data for predictions
> >>> Newdata <- head(Orthodont)
> >>> Newdata$Sex <- factor(Newdata$Sex, levels =3D
levels(Orthodont$Sex))
> >>>
> >>> # "naive" model matrix for predictions
> >>> mm1 <- model.matrix(~ poly(age, 3) + Sex, data =3D Newdata)
> >>>
> >>> # "correct" model matrix for predictions
> >>> p <- poly(Orthodont$age, 3)
> >>> mm2 <- model.matrix(~ poly(age, 3, coefs =3D attr(p,
"coefs")) + Sex, data =3D
> >> Newdata)
> >>>
> >>> data.frame(pred1 =3D predict(fm, level =3D 0, newdata =3D
Newdata),
> >> +            pred2 =3D mm1 %*% fixef(fm),
> >> +            pred3 =3D head(predict(fm, level =3D 0)),
> >> +            pred4 =3D 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
>

--=20
Renaud LANCELOT
D=E9partement Elevage et M=E9decine V=E9t=E9rinaire (EMVT) du CIRAD
Directeur adjoint charg=E9 des affaires scientifiques

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

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