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

Douglas Bates bates at stat.wisc.edu
Mon Apr 12 16:05:58 CEST 2010 

lme won't return a correlation of exactly -1 but that is a deficiency
of lme, not lmer.

Mixed-effects software should be able to handle the case of a singular
variance-covariance matrix for the random effects.  Because most
numerical methods, including those in lme, iterate with respect to the
precision matrix (inverse of the variance-covariance matrix) they
usually declare convergence at a non-singular matrix.  That is not an
advantage.  If the MLEs or REML estimates are at a singular covariance
the user should learn this.

On Mon, Apr 12, 2010 at 8:51 AM, Viechtbauer Wolfgang (STAT)
<Wolfgang.Viechtbauer at stat.unimaas.nl> wrote:
> Dear Kevin,
>
> I am curious to see what happens if you fit the same model with lme and if you play around with the optimizer used. Also, try changing the definition of the intercept. In particular:
>
> library(nlme)
> res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin)
> summary(res)
>
> res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin,
>           control=list(msVerbose=T, opt="nlm"))
> summary(res)
>
> insulin$Dose <- insulin$Dose - 8
>
> res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin)
> summary(res)
>
> res <- lme(iAUC~Treatment+Dose, random=~Dose|Subject, data=insulin,
>           control=list(msVerbose=T, opt="nlm"))
> summary(res)
>
> Best,
>
> --
> Wolfgang Viechtbauer                        http://www.wvbauer.com/
> Department of Methodology and Statistics    Tel: +31 (43) 388-2277
> School for Public Health and Primary Care   Office Location:
> Maastricht University, P.O. Box 616         Room B2.01 (second floor)
> 6200 MD Maastricht, The Netherlands         Debyeplein 1 (Randwyck)
>
>
> ----Original Message----
> From: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Kevin E.
> Thorpe Sent: Monday, April 12, 2010 15:23 To: Douglas Bates
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Another case of -1.0 correlation of random
> effects
>
>> 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=insu
>>>> lin)
>>>> 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
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

More information about the R-sig-mixed-models mailing list