Correlated Random Coefficients Across Levels

Why Random-Coefficient Covariances Matter

In a generalized linear mixed model, observations can be dependent because they share one or more grouping units. Random intercepts allow the baseline conditional linear predictor to vary across units; random slopes allow the association between a predictor and the conditional linear predictor to vary. When those random coefficients are correlated, the direction and strength of the induced dependence can change with the predictor value. With more than one grouping level, it is also important to distinguish correlations among coefficients at the same level from covariance between coefficients at different levels.

Write $Y_i$ for response $i$, $x_i$ for a scalar predictor, and $\eta_i$ for its conditional linear predictor. A useful random-intercept-and-slope form is $ \eta_i = \mathbf{x}_i^\mathsf{T}\boldsymbol{\beta} + \sum_{\ell=1}^L \mathbf{z}_i^\mathsf{T}\mathbf{b}_{\ell,g_\ell(i)} $, where $\mathbf{z}_i=(1,x_i)^\mathsf{T}$, $g_\ell(i)$ identifies the level-$\ell$ unit associated with observation $i$, and $\mathbf{b}_{\ell,j}$ is that unit's random intercept-and-slope vector. The response distribution and link complete the conditional model: $E(Y_i\mid\mathbf{b})=\mu_i$ and $g(\mu_i)=\eta_i$.

The fixed effects describe the conditional linear predictor when the random effects are zero. They are not generally population-average effects: with a nonlinear link, averaging the conditional response over the random-effects distribution does not usually preserve the linear predictor $\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}$.

Definition: For a random intercept and slope, the coefficient vector at a grouping unit is $\mathbf{b}=(b_0,b_1)^\mathsf{T}$, where $b_0$ is the unit-specific intercept deviation and $b_1$ is its slope deviation. Its covariance matrix is $$ \mathbf{G}= \begin{pmatrix} \tau_0^2 & \tau_{01}\\ \tau_{01} & \tau_1^2 \end{pmatrix}. $$ Here $\tau_0^2$ and $\tau_1^2$ are the intercept and slope variances, and $\tau_{01}$ is their covariance. Their correlation is $\rho=\tau_{01}/(\tau_0\tau_1)$, provided both standard deviations are positive. A valid covariance matrix must be positive semidefinite; in this two-coefficient case, $\tau_0^2\geq 0$, $\tau_1^2\geq 0$, and $\tau_{01}^2\leq\tau_0^2\tau_1^2$.

From Coefficient Covariance to Predictor Dependence

For a single grouping unit, the random contribution to the predictor at $x$ is $b_0+xb_1=\mathbf{z}^\mathsf{T}\mathbf{b}$. Thus, for two observations with predictor values $x_i$ and $x_j$ sharing that unit, the covariance of their random contributions is $\mathbf{z}_i^\mathsf{T}\mathbf{G}\mathbf{z}_j$. Expanding gives

\operatorname{Cov}(b_0+x_i b_1,\ b_0+x_j b_1) =\tau_0^2+(x_i+x_j)\tau_{01}+x_ix_j\tau_1^2.

This expression shows why coefficient correlation matters. The covariance term $(x_i+x_j)\tau_{01}$ can increase or decrease dependence depending on its sign and on the predictor values. For example, a positive intercept-slope covariance tends to increase the random contribution together at higher values of $x$; a negative covariance can offset some of the intercept-driven covariance. The full expression, rather than the sign of $\rho$ alone, determines the result.

At several levels, add the covariance contributions from the random effects that the two observations share, if the levels are independent. In a nested model, two observations in the same class share both their class and school effects; observations in different classes at the same school share only the school effect. In a crossed model, such as subjects crossed with items, two observations can share a subject, an item, both, or neither. Each shared grouping unit contributes its own design-vector covariance.

Worked Example: Independent Random Slopes in a Nested Model

Suppose a hypothetical Gaussian identity-link mixed model has school and class random intercepts and slopes. The school-level covariance matrix is $\mathbf{G}_S=\begin{pmatrix}0.36&0.12\\0.12&0.09\end{pmatrix}$, and the class-level matrix is $\mathbf{G}_C=\begin{pmatrix}0.16&-0.04\\-0.04&0.04\end{pmatrix}$. Assume vectors are independent across units and between levels. Calculate the covariance of the random predictor contributions for observations at $x_i=0$ and $x_j=1$.

The design vectors are $\mathbf{z}_i=(1,0)^\mathsf{T}$ and $\mathbf{z}_j=(1,1)^\mathsf{T}$. If the observations share a school, its contribution is $\mathbf{z}_i^\mathsf{T}\mathbf{G}_S\mathbf{z}_j=0.36+0.12=0.48$. If they also share a class, the class contribution is $\mathbf{z}_i^\mathsf{T}\mathbf{G}_C\mathbf{z}_j=0.16-0.04=0.12$.

Therefore, observations in the same class have random-predictor covariance $0.48+0.12=0.60$. Observations in different classes at the same school have covariance $0.48$, and observations in different schools have covariance $0$ under the stated independence assumptions.

The matrices are admissible: their diagonal entries are nonnegative, and their determinants are $0.36(0.09)-0.12^2=0.018$ and $0.16(0.04)-(-0.04)^2=0.0048$, both nonnegative. The negative class intercept-slope covariance reduces the class contribution at this pair of predictor values; it does not make the total covariance negative.

Correlations Within and Between Levels

A covariance matrix such as $\mathbf{G}_S$ describes within-level coefficient dependence: for example, whether schools with higher-than-average intercept deviations tend to have higher or lower slope deviations. It does not, by itself, describe dependence between a school coefficient and a class coefficient. A common multilevel specification assumes random-effect vectors are independent across levels. A model that instead allows cross-level covariance must state that structure explicitly and ensure the resulting joint covariance is valid.

For a school vector $\mathbf{u}_s$ and a class vector $\mathbf{v}_c$ nested within school $s$, one coherent construction is $\mathbf{u}_s\sim N(\mathbf{0},\mathbf{G}_S)$ and $\mathbf{v}_c=\mathbf{A}\mathbf{u}_s+\mathbf{e}_c$, with $\mathbf{e}_c\sim N(\mathbf{0},\mathbf{R})$ independent across classes and schools and independent of the school effects. It follows that $\operatorname{Cov}(\mathbf{u}_s,\mathbf{v}_c)=\mathbf{G}_S\mathbf{A}^\mathsf{T}$ and $\operatorname{Var}(\mathbf{v}_c)=\mathbf{A}\mathbf{G}_S\mathbf{A}^\mathsf{T}+\mathbf{R}$. This construction also determines the covariance between two class vectors in the same school: $\operatorname{Cov}(\mathbf{v}_c,\mathbf{v}_{c'})=\mathbf{A}\mathbf{G}_S\mathbf{A}^\mathsf{T}$ for distinct classes. The shared school component makes the specification coherent across all classes, rather than merely specifying one pairwise covariance in isolation.

Worked Example: A Structured School-Class Correlation

Consider hypothetical school effects with $\mathbf{G}_S=\begin{pmatrix}0.25&0.05\\0.05&0.09\end{pmatrix}$. Let $\mathbf{A}=0.5\mathbf{I}$, where $\mathbf{I}$ is the two-by-two identity matrix, and let $\mathbf{R}=\begin{pmatrix}0.09&0\\0&0.04\end{pmatrix}$. Thus, each class coefficient vector is half its school's vector plus an independent class deviation. Compare the random-predictor covariance at $x_i=0$ and $x_j=1$ for two different classes in the same school.

The total random coefficient vector in the predictor for a class-level observation is $\mathbf{u}_s+\mathbf{v}_c=1.5\mathbf{u}_s+\mathbf{e}_c$. For different classes, the $\mathbf{e}$ vectors are independent, so only the shared school component contributes. The covariance is $1.5^2\mathbf{z}_i^\mathsf{T}\mathbf{G}_S\mathbf{z}_j =2.25(0.25+0.05)=2.25(0.30)=0.675$.

For observations in the same class, the class deviation adds $\mathbf{z}_i^\mathsf{T}\mathbf{R}\mathbf{z}_j=0.09$, giving total covariance $0.675+0.09=0.765$. As a check, the cross-level covariance matrix is $\operatorname{Cov}(\mathbf{u}_s,\mathbf{v}_c)=0.5\mathbf{G}_S =\begin{pmatrix}0.125&0.025\\0.025&0.045\end{pmatrix}$. This is not the same quantity as the intercept-slope covariance within either level. The construction guarantees a valid joint distribution because it is generated from $\mathbf{u}_s$ and independent $\mathbf{e}_c$ vectors with valid covariance matrices.

For crossed groupings, the standard independent-level construction is often easier to interpret: give subjects and items separate covariance matrices, and assume the subject and item effects are independent. Two records sharing a subject then receive the subject-level covariance contribution; two records sharing an item receive the item-level contribution. A cross-covariance between subject and item effects is a different, stronger assumption. It must be defined for the whole crossed structure, not inferred from the within-subject or within-item matrices. An arbitrary collection of pairwise cross-level covariances may fail to produce a positive-semidefinite joint covariance matrix.

Why Latent-Predictor Covariance Is Not Always Response Covariance

For a Gaussian response with identity link and conditionally independent errors, covariance of the random predictor contributions translates directly into covariance of the responses. In a GLMM with a nonlinear link, that shortcut is generally wrong. Given random effects, suppose distinct observations are conditionally independent. The law of total covariance gives

\operatorname{Cov}(Y_i,Y_j) =E\{\operatorname{Cov}(Y_i,Y_j\mid\mathbf{b})\} +\operatorname{Cov}\{E(Y_i\mid\mathbf{b}),E(Y_j\mid\mathbf{b})\} =\operatorname{Cov}\{\mu_i(\mathbf{b}),\mu_j(\mathbf{b})\},\quad i\ne j.

The first term is zero by conditional independence, but the second is the covariance of the conditional means after applying the inverse link. It is not generally equal to $\operatorname{Cov}(\eta_i,\eta_j)$. If the conditional distribution has additional dependence after conditioning on the random effects, the first term need not be zero and must also be included.

Worked Example: Response Covariance in a Poisson Log-Link Model

Suppose two hypothetical counts are conditionally independent Poisson responses with log link and share a random intercept-and-slope vector $\mathbf{b}\sim N(\mathbf{0},\mathbf{G})$, where $\mathbf{G}=\begin{pmatrix}0.25&0.10\\0.10&0.16\end{pmatrix}$. Their fixed linear predictors are $\log(2)$ at $x_i=0$ and $\log(3)$ at $x_j=1$. Find their marginal means and covariance.

The random parts of the linear predictors have variances $v_i=0.25$ and $v_j=0.25+2(0.10)+0.16=0.61$. Their covariance is $c_{ij}=0.25+0.10=0.35$. For a normal random variable $W$ with mean $m$ and variance $v$, $E\{\exp(W)\}=\exp(m+v/2)$. Hence the marginal response means are $m_i=2\exp(0.25/2)=2\exp(0.125)\approx2.266$ and $m_j=3\exp(0.61/2)=3\exp(0.305)\approx4.070$.

Because the two random predictors are jointly normal, their exponentials have covariance $m_i m_j\{\exp(c_{ij})-1\}$. Therefore, $\operatorname{Cov}(Y_i,Y_j)\approx(2.266)(4.070)\{\exp(0.35)-1\} \approx9.223(0.4191)\approx3.865$. The counts have positive marginal covariance. This exact result uses the Poisson log-link model, joint normal random effects, and conditional independence; it should not be transferred unchanged to a different response family or link.

Derivation and Practical Checks

The covariance expression for a shared random-coefficient vector follows directly from bilinearity. Let $U_i=\mathbf{z}_i^\mathsf{T}\mathbf{b}$ and $U_j=\mathbf{z}_j^\mathsf{T}\mathbf{b}$, with $E(\mathbf{b})=\mathbf{0}$ and $\operatorname{Var}(\mathbf{b})=\mathbf{G}$. Then

\begin{aligned} \operatorname{Cov}(U_i,U_j) &=E\left[(\mathbf{z}_i^\mathsf{T}\mathbf{b})(\mathbf{z}_j^\mathsf{T}\mathbf{b})\right]\\ &=\mathbf{z}_i^\mathsf{T}E(\mathbf{b}\mathbf{b}^\mathsf{T})\mathbf{z}_j\\ &=\mathbf{z}_i^\mathsf{T}\mathbf{G}\mathbf{z}_j. \end{aligned}

If several independent random-effect terms are present, the covariance of their sum is the sum of the covariance contributions. If terms are correlated, cross-covariance terms must be added as well. This is why an independence assumption is not a cosmetic detail: it determines which terms belong in the calculation.

Before interpreting an estimated matrix, check that its diagonal entries are variances, its off-diagonal entries are covariances, and all predictor units are understood. Changing the origin or scale of a predictor changes the numerical intercept-slope covariance and correlation, even when it represents the same underlying model expressed in different coordinates. Also inspect whether the data contain enough within-unit predictor variation to support random slopes at each level. A covariance matrix that is numerically estimated does not guarantee that its components are well determined.

Near-zero variance: a slope variance near zero suggests little information for that random slope in the fitted model, but can also reflect weak within-unit predictor variation.
Correlation near $1$ or $-1$: this can indicate a boundary or nearly singular covariance estimate, not necessarily a strong scientific relationship that is precisely known.
Unexpected dependence: residual patterns may remain if the grouping structure or conditional response model is inadequate. Random-coefficient covariance is not a substitute for checking model fit.
Implausible cross-level structure: a cross-level covariance should have a clear structural interpretation and a valid joint covariance model; it should not be added merely because separate within-level matrices are available.

Specification	What it allows	Interpretive caution
Independent random effects across levels	Separate coefficient covariance matrices at each level; dependence arises from shared units.	It rules out cross-level covariance, rather than estimating it as zero from the data in every possible design.
Correlated intercept and slope within a level	Intercept deviations and slope deviations at the same grouping level may covary.	The covariance is predictor-scale dependent; interpret it alongside the predictor coding.
Structured cross-level covariance	Coefficients at related levels can covary according to an explicit joint model.	The full covariance structure must be coherent and supported by the design.

Common Mistake and Key Takeaway

A common mistake is to read a positive intercept-slope correlation as meaning that any two observations in the same group must have more positively correlated responses. That conclusion does not follow. The covariance of their random predictor contributions also depends on $x_i$, $x_j$, the variances, the grouping units shared, and any cross-level covariance. In a nonlinear GLMM, response covariance further depends on the inverse link and response distribution.

The practical alternative is to state the random-effects assumptions first, calculate covariance on the latent-predictor scale using the relevant design vectors and shared units, and then translate to response-scale dependence using the conditional mean model. Keep within-level coefficient covariance distinct from covariance between levels.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For $\mathbf{G}=\begin{pmatrix}0.25&-0.05\\-0.05&0.09\end{pmatrix}$, find the covariance of random contributions at $x_i=0$ and $x_j=2$ for observations sharing that unit.
In a nested model with independent school and class effects, what random-predictor covariance contributions apply to observations in different classes at the same school?
Why does a positive covariance between random intercept and slope not guarantee a positive covariance between the random contributions at every pair of predictor values?
In a conditionally independent Poisson GLMM with a log link, why is response covariance not generally equal to random-predictor covariance?

Solutions. (1) $\mathbf{z}_i=(1,0)^\mathsf{T}$ and $\mathbf{z}_j=(1,2)^\mathsf{T}$, so the covariance is $0.25+2(-0.05)=0.15$. (2) They share the school effect, so include $\mathbf{z}_i^\mathsf{T}\mathbf{G}_S\mathbf{z}_j$; they do not share a class effect, so include no class-level contribution. (3) The full expression is $\tau_0^2+(x_i+x_j)\tau_{01}+x_ix_j\tau_1^2$, whose terms depend on both predictor values; a negative slope contribution can outweigh the other terms. (4) The inverse log link transforms the random predictors before forming conditional means. The marginal response covariance is the covariance of those conditional means (under conditional independence), not simply the covariance of the linear predictors.