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

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

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

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