# Estimates of the variance of the variance component of a mixed effects model

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.

• lmer() is the function, the package is called lme4. Jul 15 '15 at 14:14
• Variance of the variance seems rather awkward in a frequentist MLE framework -- this might feel more natural expressed as hyperparameters in a Bayesian framework. Jul 15 '15 at 15:48
• 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). Jul 26 '15 at 8:51
• @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. Jan 15 '19 at 21:54
• @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. Jan 15 '19 at 22:00

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")
 "gb"


(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.

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")
 "gb"


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.

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)
> 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  • But where does it provide the estimate of$Var(\widehat{\sigma_{a}^2})\$? Jul 15 '15 at 15:05
• 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). Jul 15 '15 at 15:09
• @Livius You are misunderstanding the question. The OP is asking about the variance of the variance components. Jul 15 '15 at 15:17
• @Wolfang Ah, yes it seems I have. Jul 15 '15 at 15:46