Deriving Variance Components for Nested Hierarchies

Why a Nested Observation Has Several Sources of Variation

An observation in a hierarchy can differ from the overall mean for several reasons at once. A measurement from a student, for example, may reflect a school-wide deviation, an additional deviation specific to that student’s class, and individual-level variation. The marginal variance describes the total variation among observations after the random effects have been averaged over. Deriving it shows exactly how the levels contribute.

The result below is exact for a Gaussian nested random-intercept model. It also gives the variance of the random part of the linear predictor in a generalized linear mixed model with the same independent random intercepts. For a non-Gaussian outcome, however, that linear-predictor variance is generally not the marginal variance of the observed response.

Model and Assumptions

Consider three nested levels: level-2 units $i$, level-1 units $j$ within unit $i$, and observations $k$ within unit $j$ and unit $i$. For example, $i$ could index schools, $j$ classes within schools, and $k$ students within classes. Let $x_{ijk}^{\mathsf T}\beta$ be the fixed-effects part of the model.

Definition: A three-level nested random-intercept model is $$ Y_{ijk}=x_{ijk}^{\mathsf T}\beta+a_i+b_{ij}+\varepsilon_{ijk}, $$ where $a_i$ is the level-2 random intercept, $b_{ij}$ is the level-1 random intercept nested within level 2, and $\varepsilon_{ijk}$ is the observation-level residual. The standard variance decomposition assumes that these terms have mean zero, finite variances, and are mutually independent across terms and units, with $$ \operatorname{Var}(a_i)=\tau_2^2,\qquad \operatorname{Var}(b_{ij})=\tau_1^2,\qquad \operatorname{Var}(\varepsilon_{ijk})=\sigma^2. $$

The fixed-effects values $x_{ijk}^{\mathsf T}\beta$ are treated as fixed when taking the model-based variance. The random intercepts and residual are random. The normality assumptions $a_i\sim N(0,\tau_2^2)$, $b_{ij}\sim N(0,\tau_1^2)$, and $\varepsilon_{ijk}\sim N(0,\sigma^2)$ are needed if we also want the marginal distribution of $Y_{ijk}$ to be normal; they are not needed just to add variances. Independence and the stated variances are sufficient for the variance calculation.

This model is a special case of an $L$-level hierarchy. Each additional nested random intercept contributes one term to the observation’s random part. The lowest-level residual contributes another term, unless the model specifies a different observation-level error structure.

Deriving the Marginal Variance

“Marginal” means that the random effects are not held at particular values: their variation is included. Subtract the fixed part from the observation. Its random deviation from the fixed-effects mean is $a_i+b_{ij}+\varepsilon_{ijk}$. The variance of a sum equals the sum of variances plus twice each pairwise covariance. Under the independence assumption, those covariances are zero.

\begin{aligned} \operatorname{Var}(Y_{ijk}\mid x_{ijk}) &=\operatorname{Var}(a_i+b_{ij}+\varepsilon_{ijk})\\ &=\operatorname{Var}(a_i)+\operatorname{Var}(b_{ij}) +\operatorname{Var}(\varepsilon_{ijk})\\ &\quad+2\operatorname{Cov}(a_i,b_{ij}) +2\operatorname{Cov}(a_i,\varepsilon_{ijk}) +2\operatorname{Cov}(b_{ij},\varepsilon_{ijk})\\ &=\tau_2^2+\tau_1^2+\sigma^2. \end{aligned}

The conditioning on $x_{ijk}$ indicates that the fixed covariate values are held fixed; it does not condition on the random effects. Since each random term has mean zero, the marginal mean is $E(Y_{ijk}\mid x_{ijk})=x_{ijk}^{\mathsf T}\beta$. Thus the marginal variance is the total spread around that fixed-effects mean.

For $L$ nested random-intercept levels, let $\tau_\ell^2$ denote the variance at level $\ell$, and let $\sigma^2$ denote the observation-level residual variance. If the random intercepts and residual are independent, then

\operatorname{Var}(Y\mid x)=\sum_{\ell=1}^{L}\tau_\ell^2+\sigma^2.

This sum is a variance, in squared outcome units. Its square root is the marginal standard deviation, in the original outcome units. A component’s share of the total variance is its variance divided by the total variance. These shares sum to one when all terms in the decomposition are included.

Worked Example: Three Levels of Variation in Scores

Suppose a hypothetical Gaussian model describes a test score, measured in points, for students nested in classes within schools. Its fitted variance components are $\tau_2^2=0.16$ points squared for schools, $\tau_1^2=0.36$ points squared for classes, and $\sigma^2=1.00$ points squared for students. Assume the random intercepts and residual satisfy the independence assumptions above.

The marginal variance for one student’s score is the sum of the three components:

\operatorname{Var}(Y_{ijk}\mid x_{ijk}) =0.16+0.36+1.00=1.52\text{ points}^2.

The marginal standard deviation is $\sqrt{1.52}=1.233$ points. The school share is $0.16/1.52=0.1053$, or 10.53%; the class share is $0.36/1.52=0.2368$, or 23.68%; and the student-level residual share is $1.00/1.52=0.6579$, or 65.79%. The percentages sum to 100.0% after rounding. They describe the model’s decomposition of variation among individual scores, not causal effects of schools or classes.

How the Hierarchy Depth Changes the Total

Every included nested level contributes to the variance of an observation belonging to that level. The same rule applies whether the hierarchy has two levels or many; what changes is the number of terms in the random part. A level that is absent from the model contributes no separate component, and a random intercept at a higher level is shared by every observation in its group.

Worked Example: A Two-Level Nested Model

In a hypothetical Gaussian model, measurements are taken from clients $j$ within clinics $i$. Suppose the clinic-intercept variance is $0.49$ units squared and the residual variance is $1.21$ units squared. There is no additional random intercept for a level between clinic and observation. Assume clinic effects are independent of the residuals, and residuals for distinct observations are independent.

For one measurement, $Y_{ij}=x_{ij}^{\mathsf T}\beta+a_i+\varepsilon_{ij}$, so

\operatorname{Var}(Y_{ij}\mid x_{ij}) =0.49+1.21=1.70\text{ units}^2.

The marginal standard deviation is $\sqrt{1.70}=1.304$ units. The clinic component accounts for $0.49/1.70=0.2882$, or 28.82%, of this variance. The remaining 71.18% is observation-level residual variance. The calculation is for a single measurement; the shared clinic effect also means two measurements in the same clinic are dependent.

The decomposition extends directly to a hierarchy with more levels. Each component must be tied to the level it describes, and the nesting indices must match the sampling structure. The following calculation uses four levels to make the bookkeeping explicit.

Worked Example: Adding a Fourth Nested Level

Suppose a hypothetical Gaussian outcome in milligrams is recorded from samples nested within batches, batches within production runs, and runs within facilities. A model assigns variance components of $0.09$, $0.16$, and $0.25$ milligrams squared to facilities, runs, and batches, respectively, and residual variance $0.50$ milligrams squared to samples. Assume all four terms are mutually independent and have mean zero.

A sample’s random deviation contains one intercept from each of its four levels plus its residual. Therefore:

\begin{aligned} \operatorname{Var}(Y\mid x) &=0.09+0.16+0.25+0.50\\ &=1.00\text{ milligrams}^2. \end{aligned}

The marginal standard deviation is $\sqrt{1.00}=1.00$ milligram. The facility, run, batch, and residual shares are 9%, 16%, 25%, and 50%, respectively. Their sum is 100%. This decomposition does not mean that every sample has a distinct facility, run, and batch effect: observations sharing a level also share that level’s random intercept.

What the Formula Does—and Does Not—Say in a GLMM

In a generalized linear mixed model, the random intercepts enter the linear predictor. If the levels have independent normal random intercepts with variances $\tau_1^2,\ldots,\tau_L^2$, the random-intercept contribution to the linear-predictor variance is $\sum_{\ell=1}^{L}\tau_\ell^2$. This is a variance on the linear-predictor scale. With a nonidentity link, it is generally not the variance of the observed response.

For example, a Bernoulli response has conditional variance $p(1-p)$ given its conditional probability $p$. When random effects make $p$ vary, the marginal response variance follows the law of total variance:

\operatorname{Var}(Y) =E\{\operatorname{Var}(Y\mid B)\} +\operatorname{Var}\{E(Y\mid B)\},

where $B$ represents the random effects. Both terms depend on how the link maps the random effects to the response probability. Simply adding the random-effect variances to a Bernoulli variance is not a valid response-scale calculation. In a Gaussian identity-link model with an independent residual, by contrast, the additive decomposition above is the marginal response variance.

Assumptions, Diagnostics, and Alternatives

The additive formula depends on the model’s structure, not just on arithmetic. If random effects at different levels are correlated, the covariance terms in the derivation need not be zero. For two terms $A$ and $B$, for instance, $\operatorname{Var}(A+B)=\operatorname{Var}(A)+\operatorname{Var}(B)+2\operatorname{Cov}(A,B)$. Likewise, correlated residuals or residual variance that changes across groups require a different variance expression.

Useful checks include confirming that each observation has the intended hierarchy membership, inspecting residual patterns for unequal spread or remaining structure, and checking whether estimated variance components are near zero or unusually uncertain. A component at the boundary near zero may indicate little estimated variation at that level, but does not by itself prove the level is irrelevant. A systematic pattern in residuals can indicate that the random-intercept structure or residual assumptions are inadequate. The variance sum remains algebraically correct for the specified independent model; it cannot rescue a poor model for the data.

The main alternatives answer different modeling questions. A one-level model omits between-group variation and may therefore understate dependence among observations sharing a group. Treating group indicators as fixed effects estimates specific group contrasts rather than a distribution of group deviations. A model with varying residual variances or correlated errors changes the observation-level part of the decomposition. The appropriate choice depends on the design and scientific target; the additive random-intercept formula is not a substitute for that decision.

Finally, a decomposition is an interpretation of fitted model components, not a guarantee that the data contain enough information to estimate every component precisely. Sparse replication, few groups, misspecified nesting, or omitted dependence can make the estimated contributions unstable or misleading. Report the levels, units, assumptions, and whether the variance is on the response or linear-predictor scale.

Common Mistake

Pitfall: “The class variance is $0.36$, so a student’s total variance is $0.36$.” This confuses one component with the variance of the full observation. Under the independent three-level Gaussian model, the correct total is $0.16+0.36+1.00=1.52$ points squared. Another error is to add the random-effect variances to an observed binary outcome variance as if the response were Gaussian. In a non-Gaussian GLMM, distinguish linear-predictor variance from marginal response variance.

Check Your Understanding

For each question, give the calculation and interpret the result in context.

A Gaussian outcome has two nested random-intercept variances of $0.25$ and $0.40$, and residual variance $0.85$, all in squared units. Find the marginal variance and standard deviation of one observation.
In a three-level model, the total variance is $2.00$ units squared, and the highest-level component is $0.30$ units squared. What share of total variance is attributed to that level?
Why can’t the independent-components formula be used unchanged if a higher-level random intercept and a lower-level random intercept have nonzero covariance?
In a logistic mixed model, what scale does the sum of independent random-intercept variances describe, and why is it not generally the marginal variance of a binary response?

Solutions. (1) The variance is $0.25+0.40+0.85=1.50$ squared units, so the standard deviation is $\sqrt{1.50}=1.225$ units. (2) The share is $0.30/2.00=0.15$, or 15% of the model’s total variance. (3) The variance of the sum includes twice the covariance: $\operatorname{Var}(A+B)=\operatorname{Var}(A)+\operatorname{Var}(B)+2\operatorname{Cov}(A,B)$. A nonzero covariance adds a term that the independent-components formula omits. (4) The sum describes the random-intercept variance on the linear-predictor scale. The link transforms the random effects into response probabilities, so marginal response variance must also account for variation in those probabilities and the conditional response variance.