Another case of -1.0 correlation of random effects
Kevin E. Thorpe wrote:
Hello. 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)
Dose
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:
lmer(iAUC~Treatment+Dose+(Treatment|Subject)+(Dose|Subject),data=gluc)
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:
(Intr)
Dose -0.647
Now, a plot created by (and attached as lmlist.pdf):
plot(confint(lmList(iAUC~Dose|Subject,data=gluc,subset=Treatment=="Oat"),pooled=TRUE),order=1)
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 among Subject. Since your plots do look reasonably linear within Subject, iAUC~Treatment*Dose+(Dose|Subject) seems best to me. For what it's worth, this example closely parallels the Orthodont example in the nlme package ...
Ben Bolker 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