Covariance From Crossed Random Intercepts

Why sharing classification levels matters

In a crossed design, each observation belongs to a level of one classification and a level of another. For example, a rating might be associated with both a rater and a target; a response might be associated with both a person and an item. Two observations can share the first level, the second, both, or neither. Those relationships determine which random intercepts they have in common, and therefore how their latent predictors are correlated.

The introductory material on covariance from a nested random-intercept model treats dependence through shared levels in a hierarchy. With crossed intercepts, the key distinction is that the two classifications vary independently: sharing one does not imply sharing the other. The derivation below first gives the exact covariance on the latent linear-predictor scale, then clarifies what changes when the response is modeled through a nonlinear link.

Model and assumptions

Index observations by \(r=1,\ldots,n\). Let \(a(r)\) identify the level of classification A for observation \(r\), and let \(b(r)\) identify its level of classification B. A two-crossed-random-intercept GLMM has conditional mean \(\mu_r\) and linear predictor

Here \(g\) is the link function, \(\mathbf{x}_r\) is the fixed-effect covariate vector, \(\boldsymbol{\beta}\) is its coefficient vector, and \(u_a\) and \(v_b\) are the random intercepts for levels \(a\) and \(b\). Assume

Within each classification, random intercepts for distinct levels are independent; the two collections of random intercepts are mutually independent. The usual Gaussian specification takes \(u_a\sim N(0,\sigma_A^2)\) and \(v_b\sim N(0,\sigma_B^2)\), independently. For the response model, observations are assumed conditionally independent given all the random effects and covariates, with the chosen GLM response distribution and mean \(\mu_r=g^{-1}(\eta_r)\). The covariates are treated as fixed when calculating the covariance.

Deriving the four latent covariance cases

The fixed-effect parts of \(\eta_r\) and \(\eta_s\) are constants in this calculation, so they contribute no covariance. Expanding the remaining terms gives

The cross-covariances between an A intercept and a B intercept are zero by independence. Within classification A, two intercepts have covariance \(\sigma_A^2\) if they are the same random variable, and zero if they belong to distinct levels. The corresponding statement holds for classification B. Thus

where \(I\{\cdot\}\) equals 1 when its condition holds and 0 otherwise. This immediately yields the four cases:

The last case has a sum because the observations share two independent random intercepts. If they are distinct observations in the same A-by-B combination, they still share both intercepts. The formula is for covariance of their latent predictors, not automatically for covariance of the observed responses.

When the covariance is also a response covariance

For a Gaussian model with an identity link, the response can be written \(Y_r=\eta_r+\epsilon_r\), where the residuals are independent of the random intercepts and of each other, with mean zero and variance \(\sigma_e^2\). For distinct observations, residual covariance is zero, so the four latent covariance cases also give the response covariance cases. For a single observation,

In a general GLMM, however, the response mean is \(\mu_r=g^{-1}(\eta_r)\), a possibly nonlinear function of the random effects. For distinct observations, the law of total covariance gives

because conditional independence makes the first term zero. This second covariance is generally not equal to \(\sigma_A^2\), \(\sigma_B^2\), or their sum. Its value depends on the link, fixed predictors, random-effect distributions, and variances. If the observations share neither level, their random effects are disjoint and independent under the stated assumptions, so their responses are marginally independent and their covariance is zero. When they share a level, the response covariance must account for the inverse link.

A nonlinear-link example

The following calculation illustrates why response-scale covariance needs a separate derivation. It uses a specified, symmetric two-point distribution for an intercept rather than the conventional Gaussian assumption; the conditional Bernoulli model and the covariance identity remain valid. The example isolates the effect of a nonlinear link.

Assumptions, diagnostics, and alternatives

The four-case latent covariance formula relies on independent random-effect collections, independence across distinct levels within each collection, and a model containing the two stated random intercepts. Conditional independence is additionally needed for the simple response-scale law-of-total-covariance expression for distinct observations. If there is further dependence—such as serial dependence, another shared random effect, or a residual correlation structure—the corresponding covariance contribution must also be considered.

A practical diagnostic begins with the design, not with a covariance estimate. Verify the A and B membership labels for each row, then determine which pairs actually share each level. A data-coding error can turn a pair classified as “neither” into one that shares an intercept, or the reverse. Next, check whether residual dependence remains after accounting for the modeled structure. For a GLMM, simulated-response checks or residual summaries organized by shared A and B levels can reveal patterns the random intercepts fail to capture. A near-zero estimated variance component means that the corresponding classification may contribute little estimated dependence, but does not prove that the design labels are irrelevant.

A common mistake is to report \(\sigma_A^2\) as the response covariance for every pair sharing A. That is exact for the latent predictor and for distinct responses under the Gaussian identity-link model with independent residuals. In a nonlinear GLMM, the response covariance is instead \(\operatorname{Cov}(g^{-1}(\eta_r),g^{-1}(\eta_s))\), which must be evaluated using the model’s link and random-effect distribution. Another mistake is to add a residual variance to the covariance of two distinct observations: independent residuals add to each observation’s variance, not to their pairwise covariance.

The crossed pattern differs from a nested model because sharing a level in one classification does not determine whether a level in the other is shared. A model with a single nested hierarchy can be a simpler alternative when the sampling design genuinely has that structure. Omitting a classification is another simpler option, but it cannot represent dependence induced by that source of shared membership. Neither alternative repairs a mismatch between the assumed conditional response distribution and the data.