[R-sig-ME] Extracting variances of the estimated variance components in lme4

[ONKELINX, Thierry]

Dear Freedom,

Please note that the Std. Dev. In the model output is the standard deviation of the random effects (sigma) and thus exactly the square root of the variance (sigma ^ 2). It is *NOT* the standard error of the variance. That information is not available in the output of lme() or lmer().

If you want confidence intervals on the variance, then you have to do some math. Since (n - 1) s ^2 / sigma ^ 2  ~ chisq(n - 1). You could approxiate it by

library(lme4)
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
s2 <- VarCorr(model)$Subject[1]
n <- nrow(ranef(model)$Subject)
CI <- (n - 1) * s2 / qchisq(c(0.975, 0.025), df = n - 1)
CI

Note that confidence intervals will be very wide unless you have a very large number of random effect levels.

Best regards,

Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces op r-project.org [mailto:r-sig-mixed-models-bounces op r-project.org] Namens Freedom Gumedze
Verzonden: donderdag 3 mei 2012 14:46
Aan: r-sig-mixed-models op r-project.org
Onderwerp: [R-sig-ME] Extracting variances of the estimated variance components in lme4

Dear all,

How does one extract extracting variances of the variance components in lme4?
vcov(model) only gives the covariance matrix of fixed part of the fitted model, while VarCorr(model) only gives the estimated variance components without their corresponding standard errors.
Yes the standard errors are asymptotic but how does one extract them from the fit?

many thanks,
Freedom
>>> Ben Bolker <bbolker op gmail.com> 2012/05/03 02:31 PM >>>
Angelina Mukherjee <angelina.mukherjee88 op ...> writes:

> I have response measures corresponding to 2 patients. The structure
is as
> follows:
>
> Patient 1:  Region 1             Region 2              Region 3
>                S1 S2 S3 S4       S1 S2 S3 S4        S1 S2 S3 S4
>
> Patient 2:  Region 1               Region 2            Region 3
>                 S1 S2 S3 S4       S1 S2 S3 S4        S1 S2 S3 S4


  Hmm.  Do you really have only two patients, i.e. a total of
24 response values?  I understand that you're trying to do a variance decomposition here (no fixed effects, only random effects), but your estimates of variance will be extremely inaccurate based on only two patients (you might want to consider making patient a fixed effect, then you would at least have
6 data points (5 df) for the patient:region variance ...

> Each patient has 3 regions and each region has 4 sub-regions.
(nested
> design)
>
> Fitting   *lme( Response ~ 1, random=~ Patient + Region +Subregion |
> Patient/Region/Subregion )*
> allows me to specify covariance structure for the 'sub-region' term.
>
> But I'm trying to fit a random effects model of the form as I have
only 1
> observation per 'sub-region':
> *lme( Response ~ 1, random=~ Patient + Region | Patient/Region )*
>
> Is there a way I can specify a covariance structure like the
> auto-regressive (to specify that correlation decreases with
distance as
> one moves from Subregion 1 to Subregion 4) for the 'sub-region' term
only
> as it is not included in my random effects model but I'd like to
account
> for the correlation in it?

   I would think that something like

lme(Response~1, random = ~1|Patient/Region,
   correlation=corAR1(~Subregion))

  But I also think you're fitting a more complicated model than can really be supported by 24 data points ...

  Ben Bolker

_______________________________________________
R-sig-mixed-models op r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models





###

UNIVERSITY OF CAPE TOWN

This e-mail is subject to the UCT ICT policies and e-mai...{{dropped:19}}