I have a mixed effects growth model I call SITAR, modified from Ken Beath's (Stat Med 2007;26:2547-64), which consists of a spline curve in transformed age with random intercept (a).  Transformed age TA is a linear function of age, TA = (age – b) * exp(c), where b and c are also random effects.

Using nlme, the model provides an excellent fit to pubertal height data, explaining ~99% of the variance of the corresponding fixed effects model. However there is appreciable autocorrelation as shown by the ACF function. If I adjust for it by adding corr=corAR1() the likelihood increases by 985 units  but the residual SD also increases from 0.006 to 0.039, and the residuals are correspondingly larger. Using corr=corCAR1(~x) gives similar results.

I don't understand what's going on, and wonder if I'm missing something?

The nlme code for a curve with 5 df is as follows:


fitnlme <- function(x, s1, s2, s3, s4, s5, a, b, c) {

    splinecoefs <- as.matrix(cbind(s1, s2, s3, s4, s5))

    as.vector(a + (splinecoefs * as.matrix(ns((x-b)*exp(c),

    knots=knots, Boundary.knots=bounds))) %*% matrix(rep(1, df),  ncol=1))

}

nlme1 <- nlme(y ~ fitnlme(x, s1, s2, s3, s4, s5, a, b, c),

fixed = s1 + s2 + s3 + s4 + s5 + a + b + c ~ 1,

random = a + b + c ~ 1 | id,

…)

nlme2 <- update(nlme1, corr=corAR1())

Thanks,
Tim Cole
--
Tim.Cole@ucl.ac.uk<mailto:Tim.Cole@ich.ucl.ac.uk> Phone +44(0)20 7905 2666 Fax +44(0)20 7905 2381
MRC Centre of Epidemiology for Child Health<http://www.ucl.ac.uk/ich/research-ich/mrc-cech>
UCL Institute of Child Health, London WC1N 1EH, UK
** PLEASE- NOTE NEW EMAIL ADDRESS **

	[[alternative HTML version deleted]]

