[R-sig-ME] Another case of -1.0 correlation of random effects
bolker at ufl.edu
Fri Apr 9 14:26:15 CEST 2010
Kevin E. Thorpe wrote:
> I know this has come up a couple times recently, but I'm still not sure
> what to do about it in my data. Note that my sessionInfo() will be at
> the bottom.
> My data come from a crossover trial and are balanced.
> > str(gluc)
> 'data.frame': 96 obs. of 4 variables:
> $ Subject : int 1 2 3 5 6 7 10 11 12 13 ...
> $ Treatment: Factor w/ 2 levels "Barley","Oat": 1 1 1 1 1 1 1 1 1 1 ...
> $ Dose : int 8 8 8 8 8 8 8 8 8 8 ...
> $ iAUC : num 110 256 129 207 244 ...
> > xtabs(~Treatment+Dose,data=gluc)
> Treatment 0 2 4 8
> Barley 12 12 12 12
> Oat 12 12 12 12
> I plot the data (attached as gluc.pdf, if it comes through).
> From the plot, I think I want to fit the model as:
> It could possibly be argued that the (Treatment|Subject) part is not
> needed. When I fit this, I got -1.0 correlation within the Dose random
> effects. To simplify, I will fit a simpler model, since the issue persists.
> > lmer(iAUC~Dose+(Dose|Subject),data=gluc,subset=Treatment=="Oat")
> Linear mixed model fit by REML
> Formula: iAUC ~ Dose + (Dose | Subject)
> Data: gluc
> Subset: Treatment == "Oat"
> AIC BIC logLik deviance REMLdev
> 562.6 573.9 -275.3 563.1 550.6
> Random effects:
> Groups Name Variance Std.Dev. Corr
> Subject (Intercept) 8274.324 90.9633
> Dose 16.214 4.0266 -1.000
> Residual 4862.319 69.7303
> Number of obs: 48, groups: Subject, 12
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 309.352 30.539 10.130
> Dose -14.424 3.596 -4.012
> Correlation of Fixed Effects:
> Dose -0.647
> Now, a plot created by (and attached as lmlist.pdf):
> shows (I think) a strong negative correlation between the intercept and
> slope random effects for Dose.
> So, I would appreciate some advice on how I might specify these random
> effects correctly.
It's not entirely clear whether Subjects are unique within Treatments
(OK, maybe we should say "nested") or whether the same Subject can have
both Treatments -- what would xtabs(~Subject+Dose+Treatment) look like?
*If* each Subject appears in only one Treatment (i.e., 5 doses per
Subject, each Subject is either Oat or Barley, 12 Subjects per
Treatment), which is what I initially took for granted, then it makes
perfect sense to me that you can't specify Treatment|Subject, because
Treatment does not vary within Subject -- there doesn't seem to be any
way you can recover information about how the effect of Treatment varies
Since your plots do look reasonably linear within Subject,
seems best to me.
For what it's worth, this example closely parallels the Orthodont
example in the nlme package ...
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / people.biology.ufl.edu/bolker
GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc
More information about the R-sig-mixed-models