Thank you to all for all the background theoretical information and
practical information.
I followed your advise and example below to get the intervals for the random
effects coefficients.
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
rr1 <- ranef(fm1, postVar = TRUE)
dotplot(rr1,scales = list(x = list(relation = 'free')))[["Subject"]]
This works well in the lmer fucntion.
For my models I am using the nlme package with the lme function. The reason
I am using lme function from the nlme package ( and not the lmer fucntion
from the lme4 package) it is because i am using the option weights with
varIdent to correct for the different variances in each of my groups and the
option of corr structure with AR1 to correct for the temporal
autocorrelation. To my understanding the lme4 package can not do this yet.
is this true
So being stuck with lme (but happy!), anybody has come around the solution
to get/calculate the confidence intervals for the random coefficients? or
this is not possible for some reason?
Any advise?
Again thanks a lot for all the information before. It was usefull.
Maria
On Wed, Jan 14, 2009 at 3:46 PM, Douglas Bates wrote:
> On Wed, Jan 14, 2009 at 12:50 AM, Reinhold Kliegl
> 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 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
> >>
> >>> Hi Maria,
> >>>
> >>> Try:
> >>>
> >>> intervals(ssb_model2.1)
> >>>
> >>> Cheers,
> >>> Greg
> >>>
> >>>
> >>>
> >>> 2009/1/14 Maria Jose Juan Jorda
> >>>
> >>>> 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@pas.azti.es
> >>>> mjuanjorda@gmail.com
> >>>>
> >>>> [[alternative HTML version deleted]]
> >>>>
> >>>>
> >>>> _______________________________________________
> >>>> R-sig-mixed-models@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@r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>
> >>
> >
> > _______________________________________________
> > R-sig-mixed-models@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
--
>))):) >))):) >))):) >))):) >))):) >))):) >))):) >))):)
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@pas.azti.es
mjuanjorda@gmail.com
[[alternative HTML version deleted]]