4
$\begingroup$

Say $y=X\beta+ Zu +\epsilon$ is our mixed effects model where $u=(u_1,..,u_r)$ and $u_{j} \stackrel{i.i.d.}{\sim} N(0, \sigma^2_{a})$ for $j=1,...,r$ and $\epsilon=(\epsilon_1,...,\epsilon_n)$ are i.i.d. $N(0, \sigma^2_{b})$, furthermore $\epsilon_j$ and $u_i$ are also assumed to be independent for all $j$'s and all $i$'s.

I am interested in Var($\widehat{\sigma_{a}^2}$) i.e., the variance of the estimate of $\sigma^2_{a}$. In the R package ${\bf lme4}$ they do provide the M.L.Es but they do not provide the estimate of the Var($\widehat{\sigma_{a}^2}$). I don't think there is any closed form expression for the estimate of Var($\widehat{\sigma_{a}^2}$) and I was wondering if anyone knew of a R implementation (or any easily implementable algorithm) of how to calculate this quantity.

Improve this question
$\endgroup$
6
  • $\begingroup$ lmer() is the function, the package is called lme4. $\endgroup$
    – Livius
    Jul 15, 2015 at 14:14
  • $\begingroup$ Variance of the variance seems rather awkward in a frequentist MLE framework -- this might feel more natural expressed as hyperparameters in a Bayesian framework. $\endgroup$
    – Livius
    Jul 15, 2015 at 15:48
  • $\begingroup$ Not the same thing, but still useful if your underlying concern is the stability/error of the variance estimate: confint() gives confidence intervals for the variance components (as well as the fixed effects). $\endgroup$
    – Livius
    Jul 26, 2015 at 8:51
  • $\begingroup$ @Livius We are substantively interested in variances, including variances of variances. For example if two populations have the same mean systolic blood pressure, but different variances in SBP, the more variable population has worse health (in terms of hypertension/hypotension). Estimating variances of variances—for example state-level variance of variance of state-level SBP gives information about how states compare to one another, and insights in the range of what is possible for states to aim at. $\endgroup$
    – Alexis
    Jan 15, 2019 at 21:54
  • $\begingroup$ @Alexis for your purposes, a distributional model, such as those possible in Stan and made easy to implement via brms, might be the better tool. $\endgroup$
    – Livius
    Jan 15, 2019 at 22:00

5 Answers 5

Reset to default
4
$\begingroup$

Here is the analysis with R-package VCA V1.2:

> library(VCA)
> data(sleepstudy)
> fit <- anovaMM(Reaction~Days*(Subject), sleepstudy)
> inf <- VCAinference(fit, VarVC=TRUE)
> print(inf, what="VCA")



Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

        [Fixed Effects]

     int     Days 
251.4051  10.4673 


        [Variance Components]

  Name         DF    SS          MS         VC        %Total  SD      CV[%]   Var(VC)   
1 total        11.21                        1388.5416 100     37.2631 12.4831           
2 Subject      17    250618.1083 14742.2417 698.5289  50.3067 26.4297 8.8539  94751.0064
3 Days:Subject 17    60322.0013  3548.353   35.0717   2.5258  5.9221  1.9839  204.4845  
4 error        144   94311.5079  654.941    654.941   47.1675 25.5918 8.5732  5914.196  

Mean: 298.5079 (N = 180) 

Experimental Design: unbalanced

Fixed effects are equal and variance components of the ANOVA Type1-estimators are, except for Subject which is a bit larger (conservatively estimated), also equal to REML-estimators. Column "Var(VC)" contains variances of variance components according to Giebrecht and Burns (1985). The complete covariance matrix for variance components can also be extracted:

> vcovVC(fit)
               Subject Days:Subject       error
Subject      94751.006   -128.55799 -1523.85985
Days:Subject  -128.558    204.48451   -47.53872
error        -1523.860    -47.53872  5914.19600
attr(,"method")
[1] "gb"
Improve this answer
$\endgroup$
3
$\begingroup$

In package VCA V1.3 it is possible to use REML-estimation of linear mixed models besides ANOVA-type estimation.

> library(VCA)
> data(sleepstudy)
> fit <- remlMM(Reaction~Days+(Subject)+Days:(Subject), sleepstudy, cov=TRUE)
> fit


REML-Estimation of Mixed Model:
-------------------------------

        [Fixed Effects]

int      Days 
251.40510  10.46729 


        [Variance Components]

Name           DF         VC         %Total    SD        CV[%]     Var(VC)     
1 total        41.025787  1302.10245 100       36.084657 12.088343 82653.906666
2 Subject      9.357189   612.089747 47.007802 24.740448 8.288038  80078.294606
3 Days:Subject 11.714078  35.071663  2.693464  5.922133  1.983912  210.007398  
4 error        145.181043 654.941041 50.298733 25.591816 8.573246  5909.142918 

Mean: 298.5079 (N = 180) 

Experimental Design: unbalanced  |  Method: REML

You find the variance of variance components in column "Var(VC)". The VCA-package uses the lme4-package for REML-estimation, so the fitted model is identical to one using lmer(). Here, the variance of variance components is approximated via the method given in Giesbrecht & Burns (1985).

> vcovVC(fit)
                 Subject Days:Subject       error
Subject      80078.29461    -62.72396 -1657.13070
Days:Subject   -62.72396    210.00740   -51.91447
error        -1657.13070    -51.91447  5909.14292
attr(,"method")
[1] "gb"
Improve this answer
$\endgroup$
2
$\begingroup$

(Leaving my previous answer to the wrong question in tact for posterity, hopefully this time I'm answering the question actually being asked...)

A question about the variance of the variance estimates was recently posted on R-SIG-MIXED-MODELS. Ben Bolker, one of the lme4 authors, has already worked out how to do this for ML estimates, for REML the problem is apparently a bit harder due to the internal parameterization (links below).

The full answer is a bit long, but the basic idea is to use confidence intervals, as I suggested in my comment. Modern lme4 provides only profile and bootstrap confidence intervals for the random-effect components, which aren't as straightforwardly related to the variance/standard error of those estimates as the Wald confidence intervals are, but perhaps provide the better measure of the estimate's variability. If you do want to go the Wald confidence-interval route, from which you can rapidly compute the standard error and hence the variance on those estimates, then check out Ben Bolker's longer explication (with code). There is also an older version that not completely identically in methodology and focus (much in the same way that nlme differs from lme4, that might be worth taking a look at.

Improve this answer
$\endgroup$
1
$\begingroup$

If you are willing to fit the mixed model using ANOVA Type-1 estimation you can use R-package VCA which has two approaches to estimation of the variance of variance components implemented following Searle et al. (1992) "Variance Components" and alternatively an approximation of Giesbrecht and Burns (1985) Two-Stage Analysis Based on a Mixed Model: Large-Sample Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486.

Improve this answer
$\endgroup$
-1
$\begingroup$

The lmer function in lme4 does provide estimates of the variance of the varying slopes/intercepts, both on the variance and the standard deviation scales. 

> library(lme4)
Loading required package: Matrix
Loading required package: Rcpp
> m <- lmer(Reaction ~ Days + (Days|Subject),sleepstudy)
> m
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1743.628
Random effects:
 Groups   Name        Std.Dev. Corr
 Subject  (Intercept) 24.740       
          Days         5.922   0.07
 Residual             25.592       
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.41        10.47  
> summary(m)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
   Data: sleepstudy

REML criterion at convergence: 1743.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9536 -0.4634  0.0231  0.4634  5.1793 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  251.405      6.825   36.84
Days          10.467      1.546    6.77

Correlation of Fixed Effects:
     (Intr)
Days -0.138

As part of the REML or ML calculations, the BLUPs (more generally the conditional modes) are also computed. You can extract them with ranef():

> ranef(m)
$Subject
    (Intercept)        Days
308   2.2585637   9.1989722
309 -40.3985802  -8.6197026
310 -38.9602496  -5.4488792
330  23.6905025  -4.8143320
331  22.2602062  -3.0698952
332   9.0395271  -0.2721709
333  16.8404333  -0.2236248
334  -7.2325803   1.0745763
335  -0.3336936 -10.7521594
337  34.8903534   8.6282835
349 -25.2101138   1.1734148
350 -13.0699598   6.6142055
351   4.5778364  -3.0152574
352  20.8635944   3.5360130
369   3.2754532   0.8722166
370 -25.6128737   4.8224653
371   0.8070401  -0.9881551
372  12.3145406   1.2840295
Improve this answer
$\endgroup$
4
  • $\begingroup$ But where does it provide the estimate of $Var(\widehat{\sigma_{a}^2})$? $\endgroup$
    – user22546
    Jul 15, 2015 at 15:05
  • $\begingroup$ That's the Variance column in the random effects block. In the example I did, there are multiple variance components (intercepts and slopes for Days, grouped by subjects, as well as the residual term, which is strictly speaking also a variance component). $\endgroup$
    – Livius
    Jul 15, 2015 at 15:09
  • 2
    $\begingroup$ @Livius You are misunderstanding the question. The OP is asking about the variance of the variance components. $\endgroup$
    – Wolfgang
    Jul 15, 2015 at 15:17
  • $\begingroup$ @Wolfang Ah, yes it seems I have. $\endgroup$
    – Livius
    Jul 15, 2015 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.