[R-sig-ME] Modelling growth data: time polynomials as random effects

[Reinhold Kliegl]

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.

Reinhold

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
>