how to compute Confident intervals for BLUPS for lme function in nlme library
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.
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.
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