Covariance in Three-Level Nested Models

Why the nesting structure determines covariance

In a three-level hierarchy, observations belong to classrooms, and classrooms belong to schools. Two observations from the same classroom share both the classroom and school effects. Observations in different classrooms at the same school share the school effect only. Observations from different schools share neither. This structure determines which random effects contribute to the covariance of a pair.

The main point requires one qualification. Adding random-effect variances gives the covariance of the linear predictor, and it gives the marginal response covariance directly in a Gaussian identity-link model. In a non-Gaussian GLMM, the response mean is a nonlinear function of the random effects, so the response-scale covariance generally is not just a sum of variance components.

Definition: Index schools by $i$, classrooms within school $i$ by $j$, and observations within classroom $ij$ by $k$. A three-level nested random-intercept model has linear predictor $$ \eta_{ijk}=\mathbf{x}_{ijk}^{\mathsf T}\boldsymbol{\beta}+b_i+c_{ij}, $$ where $b_i$ is the school random intercept and $c_{ij}$ is the classroom random intercept. Use the standard assumptions $b_i\sim N(0,\sigma_S^2)$, $c_{ij}\sim N(0,\sigma_C^2)$, with all school effects independent of one another, all classroom effects independent of one another, and the two sets mutually independent. The classroom effect $c_{ij}$ is specific to its school and classroom.

For a general GLMM, the conditional distribution of $Y_{ijk}$, given the random effects, belongs to the chosen exponential-family model, with conditional mean $\mu_{ijk}$ and link $g(\mu_{ijk})=\eta_{ijk}$. Conditional on all random effects and covariates, distinct observations are assumed independent. A Gaussian identity-link model additionally has an independent residual $\epsilon_{ijk}\sim N(0,\sigma_e^2)$, so $Y_{ijk}=\eta_{ijk}+\epsilon_{ijk}$. The residual contributes to the variance of one observation, but not to the covariance between distinct observations.

Deriving the pairwise covariance

First consider the linear predictors of two distinct observations, indexed $a=(i,j,k)$ and $a'=(i',j',k')$. Fixed-effect terms are nonrandom, so they do not contribute to covariance. The random-effect terms contribute only when the two observations share the corresponding grouping unit.

\operatorname{Cov}(\eta_{ijk},\eta_{i'j'k'}) = \operatorname{Cov}(b_i+c_{ij},b_{i'}+c_{i'j'}) = \begin{cases} \sigma_S^2+\sigma_C^2, & i=i',\ j=j',\\ \sigma_S^2, & i=i',\ j\ne j',\\ 0, & i\ne i'. \end{cases}

For example, if $i=i'$ and $j\ne j'$, the school effect is the same random variable in both predictors, giving covariance $\operatorname{Var}(b_i)=\sigma_S^2$. The classroom effects are different and independent, and they are independent of the school effect, so their covariance contributions are zero. If both indices match, both random intercepts are shared. If the school indices differ, every random effect in one predictor is independent of every random effect in the other.

In the Gaussian identity-link model, the same result gives the marginal covariance of distinct responses: their independent residuals add no cross-covariance. The marginal variance of one response is $\sigma_S^2+\sigma_C^2+\sigma_e^2$. Thus the diagonal entries of the response covariance matrix include the residual variance, while off-diagonal entries follow the shared-group pattern.

Theorem: Under the stated Gaussian identity-link assumptions, two distinct observations in the same classroom have marginal covariance $\sigma_S^2+\sigma_C^2$; two observations in different classrooms at the same school have covariance $\sigma_S^2$; and observations in different schools have covariance $0$. A single observation has marginal variance $\sigma_S^2+\sigma_C^2+\sigma_e^2$.

Worked covariance calculations

Worked Example: Gaussian Covariances for Three Pair Types

Suppose a hypothetical Gaussian model for a continuous outcome has estimated variance components $\sigma_S^2=0.16$, $\sigma_C^2=0.09$, and $\sigma_e^2=0.75$, in squared outcome units. Consider distinct observations in each of the three possible relationships.

Same classroom. The pair shares both random intercepts: $$ \operatorname{Cov}(Y_{ijk},Y_{ijk'})=\sigma_S^2+\sigma_C^2=0.16+0.09=0.25. $$ Same school, different classrooms. Only the school intercept is shared: $$ \operatorname{Cov}(Y_{ijk},Y_{ij'k'})=\sigma_S^2=0.16. $$ Different schools. No random intercept is shared, so the covariance is $0$.

For an individual observation, the marginal variance is $0.16+0.09+0.75=1.00$. The corresponding correlations with distinct observations are $0.25/1.00=0.250$ for the same classroom and $0.16/1.00=0.160$ for the same school but different classrooms. The residual variance appears in this denominator because it contributes to each observation’s variance, but it does not contribute to the pair’s covariance.

Worked Example: Constructing a Covariance Matrix

Use the same hypothetical Gaussian variance components as above. Suppose four observations are recorded in this order: two from classroom 1 at school A, one from classroom 2 at school A, and one from classroom 1 at school B. Each diagonal entry is the marginal variance $1.00$. The first two observations share a classroom, the first three share school A, and the fourth is from a different school.

The first and second observations have covariance $0.25$. Each of the first three pairs involving different classrooms at school A has covariance $0.16$, except the first-second pair already identified as being in the same classroom. Every pair involving the observation from school B and an observation from school A has covariance $0$. Therefore the covariance matrix is $$ \begin{pmatrix} 1.00 & 0.25 & 0.16 & 0\\ 0.25 & 1.00 & 0.16 & 0\\ 0.16 & 0.16 & 1.00 & 0\\ 0 & 0 & 0 & 1.00 \end{pmatrix}. $$

The matrix is symmetric, as covariance requires. Its entries also reflect the nesting: the first two observations share more levels than either does with the third, and no random-effect level is shared with the fourth. A zero between-school covariance follows from the model’s independent school effects and absence of other shared random terms; it is not a universal claim about all observations from different schools.

Why a non-Gaussian response needs a different calculation

For a general GLMM, define $m_a=E(Y_a\mid\mathbf{b})$, where $\mathbf{b}$ denotes all random effects. For two distinct observations that are conditionally independent given those effects, the law of total covariance gives

\operatorname{Cov}(Y_a,Y_{a'}) = E\{\operatorname{Cov}(Y_a,Y_{a'}\mid\mathbf{b})\} + \operatorname{Cov}\{E(Y_a\mid\mathbf{b}),E(Y_{a'}\mid\mathbf{b})\} = \operatorname{Cov}(m_a,m_{a'}).

The first term is zero by conditional independence. The remaining covariance is between conditional means, which are functions of the shared random effects. With a nonlinear inverse link, their covariance is not generally equal to the covariance of the linear predictors. It can be found by integrating the product of conditional means over the joint random-effect distribution and subtracting the product of their marginal means. In contrast, under the Gaussian identity link, each conditional mean is the linear predictor, so the earlier additive formula applies exactly.

Worked Example: Response-Scale Covariance in a Probit GLMM

Consider a hypothetical binary-outcome probit model with no covariates and intercept $0$. The random intercepts have $\sigma_S^2=0.16$ and $\sigma_C^2=0.09$. For this probit model, represent each binary response as $Y=1$ when a latent normal variable $Z$ is positive, with $Z=b_i+c_{ij}+\epsilon$ and independent $\epsilon\sim N(0,1)$. This is a standard latent-variable representation of the probit link.

Each latent variable has variance $0.16+0.09+1=1.25$. For two distinct observations in the same classroom, the latent covariance is $0.16+0.09=0.25$, so their latent correlation is $0.25/1.25=0.200$. For observations in different classrooms at the same school, the latent covariance is $0.16$ and the latent correlation is $0.16/1.25=0.128$. For different schools, that covariance and correlation are zero.

Both latent variables have mean zero. For two standard, mean-zero jointly normal variables with correlation $r$, the probability that both are positive is $1/4+\arcsin(r)/(2\pi)$. Since each binary response has marginal probability $1/2$, their response covariance is $$ \operatorname{Cov}(Y_a,Y_{a'}) = P(Y_a=1,Y_{a'}=1)-\frac14 = \frac{\arcsin(r)}{2\pi}. $$

For the same classroom, this gives $\arcsin(0.200)/(2\pi)\approx0.0321$. For the same school but different classrooms, it gives $\arcsin(0.128)/(2\pi)\approx0.0204$. For different schools, it gives $\arcsin(0)/(2\pi)=0$. These are covariances of binary responses, not the latent covariances $0.25$, $0.16$, and $0$. Because each response has variance $0.5(1-0.5)=0.25$, the corresponding response correlations are approximately $0.128$ and $0.0817$ for the first two pair types.

Assumptions, checks, and alternatives

The pairwise formulas depend on the model’s structure, not just on the fact that the data are grouped. They assume the stated nesting, independent random effects at each level, constant variance components within each level, random intercepts only, and conditional independence of distinct responses. The Gaussian response formulas additionally require a Gaussian identity-link model with independent residual errors of constant variance. If covariates differ, the fixed-effect means may differ, but fixed effects do not alter the covariance formulas when treated as nonrandom.

A useful way to check the assumed pattern is to compare residual or model-based associations across pair types: same classroom, same school but different classrooms, and different schools. Strong unexplained dependence between schools, or substantial within-classroom pattern not captured by the model, can signal omitted clustering, additional dependence, or a misspecified response model. Such checks do not prove the random-effects assumptions. In particular, there may be too few higher-level groups to assess school-level variation reliably.

A random-intercept model implies equal covariance for all pairs of the same type, regardless of covariates or distance within a group. It also implies zero covariance across schools unless another shared term is included. Random slopes, serial dependence, or additional cross-classified effects can change these patterns; they are not represented by the model derived here. For non-Gaussian outcomes, covariance also depends on the inverse link and the distribution of the random effects, so latent-scale variance components alone do not fully describe response-scale association.

Common pitfall: “Two observations share a school, so their covariance is the sum of the school, classroom, and residual variances.” This is wrong for observations in different classrooms: they share only the school intercept, so the Gaussian covariance is $\sigma_S^2$. The classroom effects are independent, and independent residual errors contribute only to the individual response variances. For a non-Gaussian GLMM, even this shared-effect sum describes latent-predictor covariance rather than, in general, response covariance.

A covariance structure is a compact account of which observations share sources of unexplained variation. Compared with treating observations as independent, a nested random-intercept model represents within-school and within-classroom dependence. Compared with a model that gives every pair the same correlation, it distinguishes the stronger association expected within a classroom from the weaker association expected across classrooms in one school. The benefit depends on whether the assumed hierarchy and shared-effect structure are credible.

Check Your Understanding

Answer each question using the nested random-intercept assumptions. Solutions follow the questions.

In a Gaussian identity-link model, what is the covariance of distinct observations from different classrooms at the same school?
Suppose $\sigma_S^2=0.20$, $\sigma_C^2=0.10$, and $\sigma_e^2=0.70$. Find the marginal variance of one observation and the covariance of two observations from the same classroom.
Under the same variance components, what is the covariance of two observations from different schools, and which assumption makes it zero?
Why is adding random-effect variances generally insufficient to obtain the response covariance in a non-Gaussian GLMM?

Solutions. (1) The pair shares only the school intercept, so its covariance is $\sigma_S^2$. (2) The marginal variance is $0.20+0.10+0.70=1.00$; the same-classroom covariance is $0.20+0.10=0.30$. (3) The covariance is $0$, because school effects are independent across schools and no other random effect is shared. (4) Conditional means are transformed by a nonlinear inverse link. Response covariance is the covariance of those conditional means over the random-effect distribution, not generally the covariance of the linear predictors.