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, [EMAIL PROTECTED] 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@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>>
>
>
--
Brian D. Ripley, [EMAIL PROTECTED]
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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
