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

[Julie Marsh]

Dear Prof Kliegl,

Many thanks for your excellent comments.  I centered the TIME data and  
managed to squeeze in poly(I(TIME-100),2) as a random effect. (The  
majority of the ultrasounds occur after Day 100).  However including  
poly(I(TIME-100),3) in the model produced the familiar error message:

> ac.model3 <- lmer(y.ac ~ poly(I(TIME-100),3) +  
> (poly(I(TIME-100),3)|STUDYNO), data=s.workdat, method="ML",  
> na.action=na.omit)

Warning message:
In .local(x, ..., value) : nlminb returned message false convergence (8)

I suspect that you are correct in suggesting that there may be too  
little variation between subjects for any higher order effects.  The  
within subject variation is far greater (across time) than the  
variability between subjects (at each time point).

Thanks again for all your help.

kindest regards,  julie.



Quoting Reinhold Kliegl <reinhold.kliegl at gmail.com>:

> 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
>>
>