[R-sig-ME] Another case of -1.0 correlation of random effects

[Ben Bolker]

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