Covariance in Cross-Classified Models

Why Shared Classifications Determine Covariance

In a cross-classified design, an observation can belong to a level of each of several classifications, and the classifications are not simply nested within one another. For example, a response might be associated with both a school and a clinician. Two responses from the same school can come from different clinicians; two responses from the same clinician can come from different schools. The random intercepts induce dependence when observations share a classification level.

The key calculation is not to add every variance component in the model. For a pair of distinct observations, include only the components for classifications whose levels they share. In a Gaussian model with an identity link, this gives the marginal response covariance directly. In a GLMM with a nonlinear link, the shared random effects still determine dependence, but the response-scale covariance generally is not just a sum of variance components.

The notation \(q_k(r)\) identifies the level of classification \(k\) assigned to response \(r\). The random-intercept assumptions mean that two observations sharing a level have the same random term for that classification, while random terms for different levels are uncorrelated. Conditional independence means that any marginal association between distinct responses comes from the random effects, not from an additional residual dependence after those effects are known.

The Covariance Formula for an Identity-Link Model

First consider a Gaussian mixed model with identity link and an observation-level error: $$ Y_r=\mathbf{x}_r^\top\boldsymbol{\beta}+\sum_{k=1}^{K}u_{k,q_k(r)}+\varepsilon_r, $$ where the errors are independent of the random intercepts and of one another, with mean zero and variance \(\sigma^2\). For distinct observations \(r\) and \(s\), the fixed parts do not vary, and independent random terms contribute zero covariance. A shared random intercept contributes its variance.

The formula follows by expanding the covariance of the two random-effect sums. For classification \(k\), if the observations share a level, both contain the same random intercept, whose covariance with itself is \(\tau_k^2\). If they do not share a level, they contain independent intercepts and that covariance is zero. Covariances between different classifications are zero by the assumed independence.

The observation-level errors add \(\sigma^2\) to a response’s variance, but not to the covariance of two distinct responses, because their errors are independent. On the diagonal of the covariance matrix, each response is paired with itself and all its random intercepts are shared; the residual variance is also present. Thus the diagonal is not obtained by applying the distinct-observation formula while omitting \(\sigma^2\).

What Changes With a Nonlinear Link?

For a general GLMM, write \(\mu_r=E(Y_r\mid\mathbf{u})=g^{-1}(\eta_r)\), where \(\mathbf{u}\) denotes all random effects. Conditional independence for distinct responses gives the law-of-total-covariance result $$ \operatorname{Cov}(Y_r,Y_s)=\operatorname{Cov}(\mu_r,\mu_s),\qquad r\ne s. $$ Indeed, the law of total covariance states that the marginal covariance is the expected conditional covariance plus the covariance of the conditional means. The first term is zero here for distinct observations, leaving the second. Unlike the identity-link case, the conditional means are nonlinear functions of the random effects. Their covariance therefore depends on the link, the fixed predictors, and the random-effects distributions, not merely on a sum of variances.

A useful explicit case is a Poisson GLMM with a log link and independent normal random intercepts. For responses \(r\) and \(s\), let \(S\) be the sum of the random intercepts they share, and let \(V=\operatorname{Var}(S)\). Write their conditional means as \(\lambda_r=\exp(a_r+S+U_r)\) and \(\lambda_s=\exp(a_s+S+U_s)\), where \(a_r,a_s\) contain fixed terms and \(U_r,U_s\) contain the unshared random terms. Under the assumptions that \(S,U_r,U_s\) are jointly independent, centered normal variables, the shared term gives $$ \operatorname{Cov}(Y_r,Y_s)=E(\lambda_r\lambda_s)-E(\lambda_r)E(\lambda_s) =E(\lambda_r)E(\lambda_s)\{\exp(V)-1\}. $$ This follows from the normal moment-generating function \(E[\exp(tZ)]=\exp(t^2\operatorname{Var}(Z)/2)\) for a mean-zero normal variable \(Z\). In the product \(\lambda_r\lambda_s\), the shared term appears as \(2S\), so its exponential moment is \(\exp(2V)\); in the product of the separate marginal means, the corresponding factor is \(\exp(V)\). The unshared terms cancel from the ratio, yielding the stated factor \(\exp(V)-1\). This result requires the normal random-effects assumption and the log link; it is not a general GLMM covariance formula.

Reading the Shared-Classification Pattern

A reliable way to apply the formula is to make a short membership ledger for the pair. Record, for each classification, whether the two observations have the same level. In the identity-link model, multiply each “same” entry by its variance component and add; entries marked “different” contribute zero. In a nonlinear model, this ledger identifies the shared random effects, but a response-scale calculation must then use the appropriate link and distribution assumptions.

Assumptions, Diagnostics, and Common Pitfalls

The covariance formula is only as appropriate as its assumptions. The identity-link sum requires independent random intercepts across classifications and levels, the stated variance components, and independent observation-level errors for the diagonal variance formula. The nonlinear example additionally requires normal random effects and a log link. Conditional independence is important in both settings: if responses remain dependent after conditioning on all random effects, their conditional covariance must also be included.

Useful checks begin with the design. Confirm that the pair really shares the levels indicated by the membership records, and that the factors are represented as crossed rather than incorrectly treated as nested. Then inspect model fit: convergence warnings, variance estimates near zero or implausibly large, and residual patterns can indicate a misspecified random-effects structure or distribution. For a Gaussian model, residual dependence remaining among pairs after accounting for their shared levels challenges the conditional-independence assumption. For a GLMM, examine residual and dispersion diagnostics appropriate to the response distribution; excess variation can make a simple conditional-variance specification inadequate. Raw response covariances also need cautious interpretation in GLMMs, because the mean and link affect covariance on the response scale.

A common mistake is to add every random-intercept variance when calculating the covariance of a pair. That gives the random-intercept contribution to a single response’s variance, not the pair’s covariance. The correction is to include a component only when the pair shares a level of that classification. A second mistake is to use the Gaussian identity-link sum on the observed counts in a Poisson log-link model. There, shared effects determine the dependence, but the response-scale covariance must be derived from the conditional means and the distributional assumptions.

This calculation complements, rather than replaces, other covariance descriptions. A nested random-intercept model has a structured sharing pattern: observations in different lower-level units may still share a higher-level intercept. A cross-classified model uses the same principle, but a pair can share any subset of the crossed classifications. One could also describe the covariance through a full covariance matrix, but the shared-level formula makes each contribution interpretable. These random-effects representations are not a substitute for checking whether the design and assumptions support the covariance structure; sparse overlap, misspecified random effects, or dependence not represented in the model can limit the interpretation.