[R] polynomial predict with lme

[i.m.s.white]

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

-- 
************************************************
*    I.White                                   *
*    University of Edinburgh                   *
*    Ashworth Laboratories, West Mains Road    *
*    Edinburgh EH9 3JT                         *
*    Fax: 0131 650 6564   Tel: 0131 650 5490   *
*    E-mail: i.m.s.white at ed.ac.uk              *