[R-sig-ME] Another case of -1.0 correlation of random effects
Kevin E. Thorpe
kevin.thorpe at utoronto.ca
Mon Apr 12 15:00:58 CEST 2010
Kevin E. Thorpe wrote:
> Ben Bolker wrote:
>> Ken Knoblauch wrote:
>>> Kevin E. Thorpe <kevin.thorpe at ...> writes:
>>>> 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 ...
>>> Shouldn't you make Subject into a factor?
>> It would make the plot a little bit prettier but I don't think it
>> matters in this case because variable that appears as a grouping
>> variable (i.e. on the right of the | ) is automatically treated as a
>> factor? I think?
>> Since it is really a crossover trial, it would seem reasonable in
>> principle to have the (Treatment|Subject) random effect in there as
>> well. I'm not sure what to do about the -1 correlation: it seems the
>> choices (not necessarily in order) are (1) throw up your hands and say
>> there's not enough data to estimate independently; (2) try WinBUGS,
>> possibly with slightly informative priors; (3) try using lme4a to create
>> profiles of the parameters and see if you can figure out what's
> Let's see. I wish (1) was an option. (2) would be promising if my
> knowledge of BUGS and Bayesian methods filled more than a thimble.
> Thanks to Jarrod for his suggestion in response to this. I'll take a
> look at that too. Option (3) is probably worth a go too.
> Aside from the fact that the Dose variable are the actual doses and not
> categories, and we all know not to categorize continuous variables, what
> are your thoughts on treating Dose as a factor (since it seems to behave)?
> Thanks all for taking the time to provide your suggestions.
Regarding lme4a: how do I obtain it? I guess that comes down to, what
is the repository to give to install.packages()? Does it require a
different Matrix package than the one I have, which is, 0.999375-33 and
if so, how do I not break my current lme4/Matrix combination?
By the way, the problems with these data get stranger. In a different
outcome from the same trial the following results from a model fitting
Linear mixed model fit by REML
Formula: iAUC ~ Treatment + Dose + (Treatment | Subject) + (Dose | Subject)
AIC BIC logLik deviance REMLdev
1956 1982 -968 1983 1936
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 4.2678e-02 2.0659e-01
TreatmentOat 1.6115e+07 4.0144e+03 0.000
Subject (Intercept) 3.0430e+08 1.7444e+04
Dose 1.5173e+06 1.2318e+03 -1.000
Residual 3.1907e+07 5.6486e+03
Number of obs: 96, groups: Subject, 12
Estimate Std. Error t value
(Intercept) 40142.4 5146.7 7.800
TreatmentOat 1340.3 1634.7 0.820
Dose -2675.1 405.5 -6.597
Correlation of Fixed Effects:
Dose -0.922 0.000
As you can see, I get a 0 correlation within one set of random effects
and -1.0 in the other. Also, the fact the fixed effects estimates are
huge makes me suspicious.
Note that if I drop the treatment portions and fit the Dose model to
only one treatment, the correlation is again -1.0 and the fixed effects
Kevin E. Thorpe
Biostatistician/Trialist, Knowledge Translation Program
Assistant Professor, Dalla Lana School of Public Health
University of Toronto
email: kevin.thorpe at utoronto.ca Tel: 416.864.5776 Fax: 416.864.3016
More information about the R-sig-mixed-models