[R-sig-ME] how to compute Confident intervals for BLUPS for lme function in nlme library
Douglas Bates
bates at stat.wisc.edu
Wed Jan 14 15:46:52 CET 2009
On Wed, Jan 14, 2009 at 12:50 AM, Reinhold Kliegl
<reinhold.kliegl at gmail.com> wrote:
> Well, intervals() is a pretty good nlme response, I think, even if it
> was not exactly answering the question. The lme4 equivalent is
> doplot() which produces a so-called "caterpillar plot" of conditional
> means.
There's a typo there. The function is dotplot(), not doplot().
> For example:
> ?ranef
>
> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
> rr1 <- ranef(fm1, postVar = TRUE)
> dotplot(rr1,scales = list(x = list(relation = 'free')))[["Subject"]]
> It is important to be clear about the distinction between random
> effect estimates of variances and correlations, provided as part of
> the model output and the conditional means for groups/units based on
> ("predicted with") them. At least, correlations computed from
> conditional means can differ strongly from the correlation estimates.
> The short story: You should only trust the model estimates. As Greg
> Lee said, conditional means are to be handled with care; you can not
> apply any of the usual inferential statistics to them. Nevertheless,
> the caterpillar plots are highly informative and diagnostic about, for
> example, whether you need a random effect to account for unit-by-unit
> variability. They may also lead you to discover distinct subgroups
> suggestive of a fixed effect (for a future model). There is a
> manuscript at the top of my publications page for download (and
> constructive feedback) walking through an example.
I agree with all of the above. As Reinhold says, the BLUPs are the
conditional mean vector of the distribution of the random effects
given Y = y, the observed data and evaluated at the parameter
estimates. Sometimes I use the term "conditional modes". For a
linear mixed model the conditional distribution is multivariate
Gaussian so the mean and the mode coincide. In generalized linear
mixed models or nonlinear mixed models they do not necessarily
coincide and what is calculated is the conditional mode.
The intervals shown in the dotplot come from the evaluation of the
conditional means and the conditional variances of the random effects
given the observed data. You can obtain the conditional variances by
including the argument postVar = TRUE in a call to ranef. (The
argument is "postVar" and not "condVar" because the person who
requested it used the term "posterior variance", which is what this is
if you take a Bayesian or empirical Bayes approach.)
> Reinhold Kliegl
>
> On Wed, Jan 14, 2009 at 2:03 AM, Greg Lee <sp8ial at gmail.com> wrote:
>> Hello again Maria,
>>
>> I was on autopilot when I answered the first time. The intervals()
>> function in nlme provides confidence intervals for the estimated model
>> components (fixed effects, random effect variances and correlation
>> coefficients, if present), which is not what you were asking.
>>
>> Actually, the more I think about your question the less certain I am. The
>> coefficients in a mixed effects model are (potentially) comprised of both
>> fixed- and random-effects. But the random effects are not "estimated" so
>> much as "predicted". The only (assumed) requirement on them is that they be
>> normally distributed, with standard deviation estimates as reported by
>> intervals().
>>
>> Does it even make sense to seek confidence intervals in this context? I am
>> not sure, and leave it for wiser heads to respond.
>>
>> Regards,
>> Greg
>>
>>
>>
>> 2009/1/14 Greg Lee <sp8ial at gmail.com>
>>
>>> Hi Maria,
>>>
>>> Try:
>>>
>>> intervals(ssb_model2.1)
>>>
>>> Cheers,
>>> Greg
>>>
>>>
>>>
>>> 2009/1/14 Maria Jose Juan Jorda <mjuanjorda at gmail.com>
>>>
>>>> Dear all,
>>>>
>>>>
>>>> How do you compute the confident intervals for the estimations of the
>>>> random
>>>> coefficients? I have not seem a lot of discussion about this in the forum.
>>>> I
>>>> wonder why. Is there any theory background about this?
>>>>
>>>> Description of my data:
>>>>
>>>> I have 29 time series of data, biomass over time for 29 populations. I
>>>> want
>>>> a hierarchical linear model with varying intercepts and slopes. So it
>>>> computes 29 intercepts, 29 slopes and the overall mean of the intercept
>>>> and
>>>> slope and the variation among intercepts and slopes.
>>>>
>>>> *Here you have my linear mixed model:
>>>> *
>>>> ssb_model2.1<-lme(log(SSB_tonnes)~year2,random=~year2|id, data=ssb)
>>>>
>>>> *Data:*
>>>> SSB_tonnes: spawning stock biomass (ssb), continuous variable
>>>> year2: time in years, continuous variable
>>>> id: there are 29 id, in other words 29 populations of fish
>>>>
>>>> *Here, in the summary of the ouput* I find the oveall mean of the
>>>> interceps
>>>> and slopes and their variation.
>>>> *> summary(ssb_model2.1)*
>>>> Linear mixed-effects model fit by maximum likelihood
>>>> Data: ssb
>>>> AIC BIC logLik
>>>> 1575.647 1605.880 -781.8235
>>>>
>>>> Random effects:
>>>> Formula: ~year2 | id
>>>> Structure: General positive-definite, Log-Cholesky parametrization
>>>> StdDev Corr
>>>> (Intercept) 3.30025584 (Intr)
>>>> year2 0.02960239 -0.674
>>>> Residual 0.42259598
>>>>
>>>> Fixed effects: log(SSB_tonnes) ~ year2
>>>> Value Std.Error DF t-value p-value
>>>> (Intercept) 13.190739 0.6271420 1110 21.033098 0.0000
>>>> year2 -0.019053 0.0058418 1110 -3.261529 0.0011
>>>> Correlation:
>>>> (Intr)
>>>> year2 -0.69
>>>>
>>>> Standardized Within-Group Residuals:
>>>> Min Q1 Med Q3 Max
>>>> -5.34384192 -0.39116381 0.04824214 0.41871509 4.60300635
>>>>
>>>> Number of Observations: 1140
>>>> Number of Groups: 29
>>>>
>>>> *Here I have the estimations for the random coefficient. 29 interceps and
>>>> 29
>>>> slopes. How do you compute the confidence intervals for these estimations?
>>>> please any advise.
>>>> *
>>>> *> coef(ssb_model2.1)*
>>>> (Intercept) year2
>>>> 1 11.490207 -0.017008685
>>>> 2 11.169138 0.008690365
>>>> 3 13.373769 -0.030298582
>>>> 4 13.680485 -0.024891089
>>>> 5 15.498426 -0.033080606
>>>> 6 14.869951 -0.012568273
>>>> 7 13.016336 -0.025133918
>>>> 8 13.968938 -0.030511346
>>>> 9 13.079165 -0.017713727
>>>> 10 12.935920 -0.054212415
>>>> 11 9.124597 -0.026398183
>>>> 12 11.493282 -0.005128911
>>>> 13 10.371121 0.012636602
>>>> 14 18.273836 -0.081414152
>>>> 15 15.679493 -0.015946371
>>>> 16 5.397496 0.027787751
>>>> 17 14.479558 -0.024691465
>>>> 18 16.388451 -0.030843356
>>>> 19 18.326202 -0.098614391
>>>> 20 12.031710 -0.010605874
>>>> 21 10.352634 -0.008145436
>>>> 22 14.501495 -0.038909443
>>>> 23 14.415786 0.009461663
>>>> 24 6.915476 0.025257749
>>>> 25 7.247051 0.016988203
>>>> 26 14.519202 -0.026469286
>>>> 27 16.403880 -0.035931279
>>>> 28 18.120565 0.014417736
>>>> 29 15.407256 -0.019264309
>>>>
>>>> Thanks
>>>>
>>>> --
>>>> >))):) >))):) >))):) >))):) >))):) >))):) >))):) >))):)
>>>>
>>>> Maria Jose Juan Jorda
>>>>
>>>> AZTI - Tecnalia / Unidad de Investigación Marina
>>>> Herrera Kaia Portualde z/g
>>>> 20110 Pasaia, Gipuzkoa, Spain
>>>>
>>>> Recursos Marinos y Pesquerias
>>>> Depart. Biologia Animal, Vegetal y Ecologia
>>>> Universidade A Coruña
>>>> Campus A Zapateira s/n
>>>> 15071, A Coruña, Spain
>>>>
>>>> Tel. Oficina +34-981167000 ext. 2204
>>>> Tel. Mobil + 34-671072900
>>>> Fax. +34981167065
>>>> mjuan at pas.azti.es
>>>> mjuanjorda at gmail.com
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>>>
>>>
>>>
>>>
>>
>>
>> --
>> --------------
>> Greg Lee
>> Biometrician
>> Tasmanian Institute of Agricultural Research
>> New Town Research Laboratories
>> University of Tasmania
>> 13 St Johns Avenue,
>> New Town, 7008
>> Australia
>> Ph: +613 6233 6858
>> Fax: +613 6233 6145
>>
>> [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> 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