[R] polynomial predict with lme

[Prof Brian Ripley]

I don't believe predict.lme implements makepredictcall (which is the magic 
used in normal model-fitting), nor does it use model.frame.default.

So the answer would appear to be that there is no reason to expect it to 
work with poly().

See http://developer.r-project.org/model-fitting-functions.txt
for the standard paradigm.

On Thu, 6 Apr 2006, i.m.s.white wrote:

> Does lme prediction work correctly with poly() terms?
> In the following simulated example, the predictions
> are wildly off.
>
> Or am I doing something daft?
>
> Milk yield for five cows is measured weekly for 45 weeks.
> Yield is simulated as cubic function of weekno + random
> cow effect (on intercept) + residual error.
> I want to recover an estimate of the fixed curve.
> ###############
> library(nlme)
> set.seed(1)
> ncows <- 5; nweeks <- 45; week <- 1:nweeks
> mcurve <- 25 + 0.819*week - 0.0588*week^2 + 0.000686*week^3
> cow.eff <- rnorm(ncows)
> week <- rep(week, ncows)
> cow <- gl(ncows,nweeks)
> yield <- mcurve + cow.eff[cow] + rnorm(ncows*nweeks)
>
> lmefit <- lme(yield ~ poly(week,3), random = ~1|cow)
> summary(lmefit) # seems OK
> someweeks <- seq(5,45,by=5)
> new <- data.frame(week=someweeks)
> predicts <- predict(lmefit, new, level=0)
> print(predicts) # not even close
>
> #plot(week, yield, las=1)
> #lines(someweeks, predicts)
> ###############
>
>

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