General Covariance Derivation for Crossed Classifications

Why Derive the Covariance?

In a crossed design, an observation can share a level of one classification with another observation while differing on other classifications. For example, two measurements might come from the same person, use different items, and be collected at the same site. Their marginal covariance depends on every random effect they share, not just on whether they are “repeated” observations. A design-matrix derivation keeps these contributions distinct and also accommodates correlated random coefficients, such as a random intercept and random slope.

The key operation is to express each observation’s random contribution as a design vector multiplied by a random-effects vector. The covariance between two observations then follows from the covariance matrix of those random effects. This extends the variance-component calculations for simpler nested and crossed random-intercept models.

Model and Assumptions

Definition: Let $Y$ be an $n$-vector of Gaussian responses and $X$ a fixed-effects design matrix. For each of $F$ crossed classifications, let $Z_f$ be the random-effects design matrix and $b_f$ the vector of random coefficients. The model is $$ Y=X\beta+\sum_{f=1}^{F}Z_f b_f+\varepsilon. $$ For classification $f$, suppose there are $m_f$ levels and $q_f$ random coefficients per level. The level-specific coefficient vectors are independent and identically distributed as $N(0,G_f)$, so $b_f\sim N(0,I_{m_f}\otimes G_f)$. Assume the random-effect vectors for different classifications are mutually independent, and independent of $\varepsilon$, with $\varepsilon\sim N(0,R)$. Each $G_f$ and $R$ must be symmetric positive semidefinite.

These are assumptions about the model, not consequences of a crossed design. In particular, the model assumes zero means for the random effects, independence across levels within a classification, and independence between classifications. The covariance matrix $G_f$ need not be diagonal: its off-diagonal entries allow coefficients at the same level to be correlated. The residual covariance $R$ may be a general positive semidefinite matrix; independent residual errors with common variance $\sigma^2$ are the special case $R=\sigma^2 I_n$.

For observation $i$, its row of $Z_f$ has the vector $z_{if}^{T}$ in the block for the level of classification $f$ assigned to that observation, and zeros in all other level blocks. Here $z_{if}$ is a $q_f$-vector of coefficient covariates. A random intercept has a component equal to 1; a random slope has a component equal to the corresponding covariate value.

Marginal Covariance: Matrix Derivation

Conditionally on the random effects, the mean of $Y$ is $X\beta+\sum_f Z_f b_f$. Marginally, the random effects have mean zero, so $E(Y)=X\beta$. To obtain the covariance, subtract this mean and expand the covariance of the sum. Independence makes every covariance between distinct random-effect classifications zero, as does independence between the random effects and residuals. Therefore,

\begin{aligned} \operatorname{Var}(Y) &=\operatorname{Var}\left(\sum_{f=1}^{F}Z_f b_f+\varepsilon\right)\\ &=\sum_{f=1}^{F}Z_f\operatorname{Var}(b_f)Z_f^{T}+R\\ &=\sum_{f=1}^{F}Z_f(I_{m_f}\otimes G_f)Z_f^{T}+R =:V. \end{aligned}

The second equality uses $\operatorname{Var}(A b)=A\operatorname{Var}(b)A^T$ for a fixed matrix $A$. Under the stated Gaussian assumptions, this also gives the marginal distribution $Y\sim N(X\beta,V)$. The covariance result itself follows from the second-moment assumptions even without Gaussianity.

The corresponding pairwise expression makes the role of shared levels explicit. Let $a_f(i)$ denote the level of classification $f$ assigned to observation $i$. For distinct observations $i$ and $j$,

\operatorname{Cov}(Y_i,Y_j) = \sum_{f=1}^{F} \mathbf{1}\{a_f(i)=a_f(j)\} \,z_{if}^{T}G_f z_{jf} +R_{ij}.

If the observations do not share a level of classification $f$, their corresponding random-effect contribution is zero because different level-specific vectors are independent. If they do share a level, the contribution is the covariance of two linear combinations of that level’s coefficients. The residual term $R_{ij}$ is added separately. With independent residual errors, it is zero for distinct observations. For a single observation, its variance is the sum of all its random-effect contributions and $R_{ii}$.

Correlated Random Coefficients

For a random intercept and random slope on a covariate $x$, the design vector is $z(x)=(1,x)^T$. Write the coefficient covariance matrix as

G= \begin{pmatrix} \tau_{00} & \tau_{01}\\ \tau_{01} & \tau_{11} \end{pmatrix}.

Here $\tau_{00}$ is the random-intercept variance, $\tau_{11}$ is the random-slope variance, and $\tau_{01}$ is their covariance. For two observations sharing that level, the contribution to their covariance is

z(x_i)^T G z(x_j) = \tau_{00}+\tau_{01}(x_i+x_j)+\tau_{11}x_i x_j.

Thus the shared-level covariance can vary with both covariate values. Setting $\tau_{01}=0$ is an additional restriction, not a general property of random slopes. More generally, $G_f$ represents covariance among all $q_f$ random coefficients at a level, and the product $z_{if}^{T}G_fz_{jf}$ handles them together, including their cross-covariances.

Worked Example: Two Crossed Random Intercepts

Suppose hypothetical measurements are indexed by participant and task. The model has an independent participant random intercept with variance $4.00$, a task random intercept with variance $1.50$, and independent residual errors with variance $2.00$. Assume all random effects have mean zero and are mutually independent. Consider two distinct measurements from the same participant on different tasks.

Each random-intercept design entry is 1. The pair shares a participant, but not a task. The residual covariance is zero because the records are distinct and the residuals are independent. The pairwise formula gives

\operatorname{Cov}(Y_i,Y_j) =(1)(4.00)(1)+(0)(1.50)(1)+0 =4.00.

For either observation, the participant, task, and residual contributions all enter its variance:

\operatorname{Var}(Y_i)=4.00+1.50+2.00=7.50.

The marginal correlation is therefore $4.00/\sqrt{7.50(7.50)}=0.5333$. These calculations follow from the specified model: the shared participant contributes to covariance; the different tasks and independent residuals do not. If instead the observations came from different participants but the same task, their covariance would be $1.50$.

Worked Example: A Correlated Random Intercept and Slope

Suppose a hypothetical clinic-level random intercept and slope have covariance matrix

G= \begin{pmatrix} 9.00 & -1.20\\ -1.20 & 0.64 \end{pmatrix},

and the residual variance is $1.50$, independently across records. Two distinct patients at the same clinic have covariate values $x_i=-1$ and $x_j=2$. The design vectors are $z_i=(1,-1)^T$ and $z_j=(1,2)^T$. First, the proposed $G$ is a valid covariance matrix: its diagonal entries are positive and its determinant is $9.00(0.64)-(-1.20)^2=4.32>0$, so it is positive definite.

The clinic contribution to the covariance is

\begin{aligned} z_i^T G z_j &=(1,-1) \begin{pmatrix} 9.00 & -1.20\\ -1.20 & 0.64 \end{pmatrix} \begin{pmatrix}1\\2\end{pmatrix}\\ &=(1,-1)\begin{pmatrix}6.60\\0.08\end{pmatrix} =6.60-0.08=6.52. \end{aligned}

There is no residual covariance between the distinct patients, so their marginal covariance is $6.52$. Their marginal variances include the residual variance. For $x_i=-1$, the random-effect variance is $9.00+2(-1)(-1.20)+(-1)^2(0.64)=12.04$, giving $\operatorname{Var}(Y_i)=13.54$. For $x_j=2$, it is $9.00+2(2)(-1.20)+2^2(0.64)=6.76$, giving $\operatorname{Var}(Y_j)=8.26$. Their marginal correlation is $6.52/\sqrt{13.54(8.26)}\approx0.6166$. The negative intercept-slope covariance is included in both variance calculations and in the pairwise covariance; omitting it would change the result.

Worked Example: Adding Contributions From Three Crossed Classifications

Consider hypothetical observations classified by person, item, and site. Each classification has an independent random intercept, with variances $2.25$, $0.81$, and $1.44$, respectively. The independent residual variance is $1.00$. Suppose records A and B are from the same person and site but different items. Their random-intercept design entry is 1 for each classification, and their residual covariance is zero. Thus

\operatorname{Cov}(Y_A,Y_B) =1(2.25)1+0(0.81)1+1(1.44)1+0 =3.69.

Each record has marginal variance $2.25+0.81+1.44+1.00=5.50$, so their marginal correlation is $3.69/5.50\approx0.6709$. A third record from a different person and site but the same item as A would share only the item effect with A; their covariance would be $0.81$. The calculation adds the contributions classification by classification rather than treating “same study” or “same setting” as a single covariance category.

Constructing and Checking the Covariance

A reliable calculation follows from the model specification and the membership design, not from visual guesses about which observations seem related. The following sequence makes the assumptions and arithmetic visible.

Specify coefficient vectors.
For each classification, state which random coefficients occur at each level and their covariance matrix $G_f$.

Build membership design rows.
Place each observation’s coefficient vector in the block for its assigned level; use zeros in the other level blocks.

Form each covariance contribution.
Use $Z_f(I_{m_f}\otimes G_f)Z_f^T$, or for a pair use the shared-level term $z_{if}^TG_fz_{jf}$.

Add the residual covariance.
Include $R$ once, then check that the resulting $V$ is symmetric positive semidefinite.

The matrix expression also justifies the final check: each term $Z_f(I_{m_f}\otimes G_f)Z_f^T$ is positive semidefinite when $G_f$ is, because for any vector $c$, $c^TZ_f(I_{m_f}\otimes G_f)Z_f^Tc=(Z_f^Tc)^T(I_{m_f}\otimes G_f)(Z_f^Tc)\geq0$. Since $R$ is positive semidefinite, their sum $V$ is too. A computed covariance matrix that is asymmetric or has a negative variance on its diagonal signals an algebra, coding, or model-parameter problem.

Alternatives, Diagnostics, and Limitations

The pairwise formula is usually the quickest hand calculation, especially for random intercepts. It is a special case of the matrix formula, not a substitute for it when coefficients are correlated or the design is complex. The matrix approach is systematic and readily extends to many classifications, while a conditional description focuses on the response given the random effects and does not by itself display the marginal dependence between observations. All three views are compatible when they use the same assumptions.

Several checks can reveal when a proposed derivation or model is going wrong:

Check membership level by level. A nonzero shared-level contribution is included only when the two observations actually share that level. Confusing classification labels can create spurious covariance.
Check design-vector order. The order of components in $z_{if}$ must match the row and column order in $G_f$. Otherwise the intercept-slope covariance may be applied to the wrong coefficients.
Check covariance matrices. Each $G_f$ and $R$ must be symmetric positive semidefinite. A negative variance or invalid covariance matrix cannot describe the assumed random vector.
Compare implied patterns with residual dependence. If observations with no shared random-effect levels still show systematic dependence, the assumed $R$ or random-effects structure may be inadequate.

This derivation depends on the stated independence between classifications. If random effects from different classifications are correlated, cross-covariance terms must also be included; they are not part of the independent-classification formula above. Likewise, this page derives the marginal covariance conditional on the model parameters; it does not establish that the design can identify or estimate every component separately. Even a mathematically valid $V$ can represent an unsuitable dependence structure if the assumptions do not match the data-generating process.

Common Pitfall

Common mistake: “Both observations belong to the same clinic, so their covariance is the random-intercept variance.” This is incomplete when the clinic also has random slopes, when the observations have different covariate values, or when other crossed classifications are shared.

Correction: For every classification, check whether the observations share a level and calculate $z_{if}^{T}G_fz_{jf}$ using the full coefficient covariance matrix. Then add the residual covariance $R_{ij}$. The random-intercept-only shortcut is correct only when it matches the specified model.

Key takeaway: Under independent crossed classifications, the marginal covariance is the sum of the design-matrix contributions $Z_f(I_{m_f}\otimes G_f)Z_f^T$ and the residual covariance $R$. Within a shared level, correlated random coefficients contribute through the full bilinear form $z_{if}^TG_fz_{jf}$.

Check Your Understanding

Use the model assumptions in each question. Solutions follow the questions.

Two distinct observations share a level of a random-intercept classification with variance $3$, but share no other random-effect levels. Their residuals are independent. What is their marginal covariance?
A shared-level random intercept and slope have $G=\begin{pmatrix}4&0.5\\0.5&1\end{pmatrix}$. Find the random-effect covariance contribution for covariate values $x_i=0$ and $x_j=2$.
Why does the marginal covariance matrix remain positive semidefinite under the stated assumptions?
Two observations share no levels of any random-effect classification. Their residual covariance is $0.30$. What is their marginal covariance?

Solutions. (1) The shared random intercept contributes $1(3)1=3$; all other contributions are zero, so the covariance is $3$. (2) The design vectors are $(1,0)^T$ and $(1,2)^T$. Thus $(1,0)G(1,2)^T=(1,0)(5,2)^T=5$, so the contribution is $5$. (3) Each random-effect term $Z_f(I_{m_f}\otimes G_f)Z_f^T$ is positive semidefinite because $G_f$ is; $R$ is also positive semidefinite. Their sum is therefore positive semidefinite. (4) No random-effect contribution is shared, so only $R_{ij}$ remains and the covariance is $0.30$.