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

[Kevin E. Thorpe]

Douglas Bates wrote:
> On Mon, Apr 12, 2010 at 8:00 AM, Kevin E. Thorpe
> <kevin.thorpe at utoronto.ca> wrote:
>> 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 ...
>>>>>>
>>>>>> clip>
>>>>> Shouldn't you make Subject into a factor?
>>>>>
>>>>> Ken
>>>>>
>>>>  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
>>>> happening.
>>> 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.
>>>
>>> Kevin
>>>
>> 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?
> 
> The repository is http://R-forge.R-project.org but be aware that lme4a
> is under active development and not guaranteed against breakage.  It
> would be inadvisable to rely on functions and classes in that package
> to persist.
> 
>> 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
>> attempt.
>>
>> lmer(iAUC~Treatment+Dose+(Treatment|Subject)+(Dose|Subject),data=insulin)
>> Linear mixed model fit by REML
>> Formula: iAUC ~ Treatment + Dose + (Treatment | Subject) + (Dose | Subject)
>>   Data: insulin
>>  AIC  BIC logLik deviance REMLdev
>>  1956 1982   -968     1983    1936
>> Random effects:
>>  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
> 
> And did you notice that you have an Intercept term by Subject in there
> twice?  It is not surprising that the parameter estimates are
> unstable.  You will need to rethink the model.

Thanks Doug.  Any suggestions?  Note that the instability is present 
without the random effects for treatment in the model too.

> 
>>  Residual              3.1907e+07 5.6486e+03
>> Number of obs: 96, groups: Subject, 12
>>
>> Fixed effects:
>>             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:
>>            (Intr) TrtmnO
>> TreatmentOt -0.079
>> 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 are
>> similar.
>>
>>


-- 
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