Modelling growth data: time polynomials as random effects

You may want to center "TIME" prior to running lmer. Probably, it
 probably better to use "poly(TIME, 3)".

 For example, it works here:

 > (fm1 <- lmer(Reaction ~ poly(Days,2) + (poly(Days,2)|Subject), sleepstudy))

 Linear mixed model fit by REML
 Formula: Reaction ~ poly(Days, 2) + (poly(Days, 2) | Subject)
   Data: sleepstudy
  AIC  BIC logLik deviance REMLdev
  1733 1765 -856.4     1738    1713

Random effects:
  Groups   Name           Variance Std.Dev. Corr
          (Intercept)     1422.41  37.715
          poly(Days, 2)1 54543.55 233.546   0.732
          poly(Days, 2)2 20260.08 142.338  -0.140 -0.025
  Residual                  518.16  22.763
 Number of obs: 180, groups: Subject, 18


 Fixed effects:
               Estimate Std. Error t value
 (Intercept)      298.51       9.05   32.98
 poly(Days, 2)1   403.36      59.57    6.77
 poly(Days, 2)2    32.86      40.54    0.81

 Correlation of Fixed Effects:
            (Intr) p(D,2)1
 ply(Dys,2)1  0.665
 ply(Dys,2)2 -0.114 -0.019

 On a practical side, it could also be that there simple is too little
 variance between subjects for higher-order TIME effects.

 On Sun, Apr 6, 2008 at 12:18 PM, Julie Marsh
 <marshj02 at student.uwa.edu.au> wrote:

 >
 >  Dear LMER Experts,
 >
 >  I am trying to model some growth data based on multiple ultrasound
 >  measurements taken across the time course (in days) of each pregnancy.
 >  The vast majority of subjects have 4+ measurements.  Everything seems
 >  good up to:
 >

 >  > model.ac1 <- lmer(y.ac ~ TIME + I(TIME^2) + I(TIME^3) +
 >  > (TIME|STUDYNO), data=s.workdat, method="ML", na.action=na.omit)

 >

 >  > summary(model.ac1)

 >  Linear mixed-effects model fit by maximum likelihood
 >  Formula: y.ac ~ TIME + I(TIME^2) + I(TIME^3) + (TIME | STUDYNO)
 >     Data: s.workdat
 >     AIC   BIC logLik MLdeviance REMLdeviance
 >   -5051 -5004   2532      -5065        -4978
 >  Random effects:
 >   Groups   Name        Variance   Std.Dev.  Corr
 >   STUDYNO  (Intercept) 6.3577e-02 0.2521444
 >            TIME        1.7435e-06 0.0013204 -0.828
 >   Residual             1.3597e-02 0.1166055
 >  number of obs: 6066, groups: STUDYNO, 1289
 >
 >  Fixed effects:
 >                Estimate Std. Error t value
 >  (Intercept)  1.309e+00  1.313e-01   9.970
 >  TIME         5.825e-02  2.160e-03  26.963
 >  I(TIME^2)   -1.113e-04  1.144e-05  -9.723
 >  I(TIME^3)    7.197e-08  1.957e-08   3.678
 >
 >  Correlation of Fixed Effects:
 >            (Intr) TIME   I(TIME^2
 >  TIME      -0.996
 >  I(TIME^2)  0.988 -0.997
 >  I(TIME^3) -0.976  0.991 -0.998
 >
 >
 >  However as the slope consists of the polynomial combination (TIME +
 >  TIME^2 + TIME^3) intuitively it would seem correct to include these
 >  polynomial terms as  random effects. However, everything goes
 >  pair-shaped when I try:
 >

 >  > test.ac3 <- lmer(y.ac ~ TIME + I(TIME^2) + I(TIME^3) +
 >  > (TIME+I(TIME^2)+I(TIME^3)|STUDYNO), data=s.workdat, method="ML",
 >  > na.action=na.omit)

 >
 >  Warning messages:
 >  1: In .local(x, ..., value) :
 >    Estimated variance-covariance for factor 'STUDYNO' is singular
 >
 >  2: In .local(x, ..., value) :
 >    nlminb returned message false convergence (8)
 >
 >
 >  I am assuming that this is due to the strong correlations between the
 >  polynomials of TIME. I also performed this analysis on the subset of
 >  subjects with 5 or greater ultrasound measurements (in desperation!
 >  n=1024) and obtained the same error message.
 >
 >  My dilemma is explaining why TIME has both a fixed and random
 >  component, whereas TIME^2 and TIME^3 only have a fixed component.  I
 >  suspect I am missing something fundamental !!!
 >
 >  I have also had some fun fitting different AR structures but decided
 >  to strip back the model to the basic correlation structure for this
 >  question. Any help would be very much appreciated.
 >
 >
 >  kindest regards,  julie marsh
 >
 >  PhD Student
 >  Centre for Genetic Epidemiology
 >  University of Western Autsralia
 >  email: marshj02 at student.uwa.edu.au
 >
 >  _______________________________________________
 >  R-sig-mixed-models at r-project.org mailing list
 >  https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
 >