[R-sig-ME] Modelling growth data: time polynomials as random effects
Reinhold Kliegl
reinhold.kliegl at gmail.com
Sun Apr 6 12:51:24 CEST 2008
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
>
More information about the R-sig-mixed-models
mailing list