Many thanks for your very useful comments and suggestions. Renaud
2006/5/30, Prof Brian Ripley <[EMAIL PROTECTED]>: > 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 > -- 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 ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
