Lagged term in time series with stationary errors: too good to be true?

Question

I often have datasets where there are many animals, in several treatment groups, and each animal's body weight is measured at regular intervals over the course of its lifetime. The response of body weight is non-monotonic over time: when the animals are young, they grow. Then there is a varying steady-state period. Then, as the animals get old and sick their weight starts to decline. Treatments can alter the rate of weight-gain, onset and duration of the steady-state, the maximum weight attained, and the onset and rate of the terminal decline. Currently I use mixed-effect models and attempt to make the weight monotonic in various ways, all of which either don't have a satisfactory interpretation and/or don't produce well-fitting models.

I tried simply including a lagged body weight term as a random slope, and got beautiful residuals. Here is the R code.

library(nlme); library(gridExtra);
## get anonymized data, an data.frame named lagbwtest will be created
source('http://pastebin.com/download.php?i=MmSDVBEr');
## id = grouping variable, treatment1 & treatment2 = fixed effects
## time = time, y = response variable (body weight), lagy = lagged response
lagbwtest<-na.omit(lagbwtest);
## model with a lagged response + time random slope
fit01<-lme(y ~ (treatment1+treatment2)*time, data=lagbwtest, random=~time+lagy|id, method="ML", na.action=na.exclude);
## model without lagged random slope
fit02<-lme(y ~ (treatment1+treatment2)*time, data=lagbwtest, random=~time|id, method="ML", na.action=na.exclude);
## model with polynomial time term
fit03<-update(fit02,.~(treatment1+treatment2)*poly(time),random=~poly(time)|id);
## model with spline time term
fit04<-update(fit02,.~(treatment1+treatment2)*ns(time),random=~ns(time)|id);
## make a list to automate doing the same thing to each fit
myfits<-list(lag=fit01,nolag=fit02,poly=fit03,spline=fit04);

If we compare the goodness of fit, the lag model seems far better than the others:

with(myfits,anova(lag,nolag,poly,spline));

The lag model gives different predictions than the others (please forgive the baroque lapply idioms-- just trying to keep the code as brief as I can):

do.call(grid.arrange,lapply(names(myfits),function(ii) plot(myfits[[ii]],fitted(.,level=0)~time|treatment1:treatment2,pch='.',abline=0,main=ii)))

The lagged model gives a different result than all the others. My questions are:

1. Have I really solved the problem of finding the right functional form for such data, or have I just swept it under the rug by making the time-dependent response part of the error term?
2. Given that the error is stationary (confirmed in the ACF and resid vs time plots at the end of this post) is this a valid way to introduce an autocorrelated variable into a linear model? If no, why not?
3. Can the lag model be used to predict and test hypotheses about the maximum y value attained in each treatment combination?
4. The level=1 (animal-level) plots for the lagged model generally look better than the other modesl but level=0 (group-level) plots all look horrible and it occurs to me that I don't know how the level=0 residuals from a well-fitted mixed-effect model should look. Do they even matter?
5. Why do the non-lagged models all come out so similar?

Here is the code for the diagnostic plots. For two of them the actual images are included:

## autocorrelation
do.call(grid.arrange,lapply(names(myfits),function(ii) plot(ACF(myfits[[ii]],form=~time|id),main=ii)));
## spatial correlation (variograms)
do.call(grid.arrange,lapply(names(myfits),function(ii) plot(Variogram(myfits[[ii]],form=~time|id),main=ii)))
## individual animal residuals vs fitted
do.call(grid.arrange,lapply(names(myfits),function(ii) plot(myfits[[ii]],resid(.,type='pearson',level=1)~fitted(.,level=1),pch='.',abline=0,main=ii)))

## individual animal residuals vs time
do.call(grid.arrange,lapply(names(myfits),function(ii) plot(myfits[[ii]],resid(.,type='pearson',level=1)~time,pch='.',abline=0,main=ii)))

## treatment-group residuals vs fitted
do.call(grid.arrange,lapply(names(myfits),function(ii) plot(myfits[[ii]],resid(.,type='pearson',level=0)~fitted(.,level=0),pch='.',abline=0,main=ii)))
## treatment-group residuals vs time
do.call(grid.arrange,lapply(names(myfits),function(ii) plot(myfits[[ii]],resid(.,type='pearson',level=0)~time,pch='.',abline=0,main=ii)))

By the way, I also tried updating the nolag model with correlation=corAR1() (not shown) and it did not make the shapes of any of the plots more similar to those of the lag model, not even the ACF plots.