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

[Douglas Bates]

On Sun, Apr 6, 2008 at 5:51 AM, Reinhold Kliegl
<reinhold.kliegl at gmail.com> wrote:
> 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.

I agree.  The Pixel data analyzed in "Mixed-effects Models in S and
S-PLUS" (and available in the MEMSS package) is such an example.  The
overal shape of the density versus time curve is quadratic but the
variation between dogs can only be modeled up to the linear term and
the variation in side within dog as an additive shift only.

>  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
>  >
>
>  _______________________________________________
>  R-sig-mixed-models at r-project.org mailing list
>  https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>