From Membership to Covariance
In a hierarchical or cross-classified model, two observations can be dependent because they share one or more random-effect units. The membership structure says which units they share; the random-effect variances say how much each shared unit contributes to their covariance. An incidence matrix turns that structure into a calculation that can be checked row by row.
The central construction is simple for a Gaussian model: each random-effect term contributes a matrix of the form \(Z_kG_kZ_k^{\mathsf T}\), and the contributions are added to the residual covariance. The same sum describes covariance in the linear predictor when random effects enter it linearly. For non-Gaussian GLMM responses, however, the marginal response covariance generally does not equal that sum; the distinction is essential.
For independent random-effect terms, write the linear predictor as \(\eta=X\beta+\sum_k Z_kb_k\), with \(E(b_k)=0\) and \(\operatorname{Cov}(b_k)=G_k\). Independence gives \(\operatorname{Cov}(\eta)=\sum_k Z_kG_kZ_k^{\mathsf T}\). In a Gaussian identity-link model with residual covariance \(R\), the marginal covariance of the response is \(V=\sum_k Z_kG_kZ_k^{\mathsf T}+R\). A common special case has \(R=\sigma^2 I\), where \(I\) is the \(n\)-by-\(n\) identity matrix. Fixed effects affect the mean, not this covariance, when their values are treated as fixed.
With random intercepts, \(G_k=\tau_k^2 I\), so the \((i,j)\) entry of a term's contribution is \(\tau_k^2\) if observations \(i\) and \(j\) share a unit and 0 otherwise. Its diagonal entries are also \(\tau_k^2\), because every observation belongs to a unit. The full covariance is therefore obtained by adding the contributions of every shared membership level, as well as the residual variance on the diagonal.
Nested Membership Produces Blocks
For nested random intercepts, a lower-level unit belongs to exactly one higher-level unit. For example, each class belongs to one school. Use one incidence matrix for schools and another for classes; class columns identify distinct classes across the whole data set. Each observation then has one membership entry in each matrix. The school contribution links observations in the same school, while the class contribution links observations in the same class. With observations ordered by school, these patterns appear as blocks, and the smaller class blocks sit within the school blocks.
Worked Example: Schools and Classes
Consider five hypothetical observations ordered as follows: two observations in school A, class A1; one in school A, class A2; and two in school B, class B1. Suppose the model has independent school random intercepts with variance 4, independent class random intercepts with variance 2, and independent residual errors with variance 1. Assume the Gaussian identity-link model, so the marginal response covariance is the sum of these components.
The school incidence matrix has columns A and B; the class incidence matrix has columns A1, A2, and B1. In the stated observation order, their rows are respectively \((1,0),(1,0),(1,0),(0,1),(0,1)\) and \((1,0,0),(1,0,0),(0,1,0),(0,0,1),(0,0,1)\). For any two observations, add 4 if they share a school, add 2 if they share a class, and add 1 to the diagonal for residual variance.
For example, observations 1 and 2 share both school A and class A1, so their covariance is \(4+2=6\). Observations 1 and 3 share school A but not a class, giving covariance 4. Observations 1 and 4 share neither, giving covariance 0. Each diagonal entry is \(4+2+1=7\). Thus
The first three observations form one school block and the last two another. Within each block, the class membership adds covariance only for observations in the same class. This matrix follows from the stated model assumptions; the pattern would change if the class or school membership were recorded incorrectly, or if random effects at different levels were correlated.
Crossed Membership Produces Overlapping Patterns
In a crossed structure, an observation belongs to one unit from each of two classifications, but neither classification is nested within the other. Subjects and items are a common example: a subject can respond to multiple items, and an item can be answered by multiple subjects. A subject incidence matrix and an item incidence matrix produce separate covariance patterns. Two observations can share a subject, share an item, share both, or share neither; each random-intercept term contributes according to its own membership, regardless of the other term.
Worked Example: Subjects Crossed With Items
Suppose four hypothetical responses are observed in this order: subject 1 with item 1, subject 1 with item 2, subject 2 with item 1, and subject 2 with item 2. Assume independent subject random intercepts with variance 1.5, independent item random intercepts with variance 2.5, and independent Gaussian residual errors with variance 0.5.
The diagonal variance for each response is \(1.5+2.5+0.5=4.5\). Responses 1 and 2 share subject 1, so their covariance is 1.5; responses 1 and 3 share item 1, so their covariance is 2.5. Responses 1 and 4 share neither and have covariance 0. Applying this rule to every pair gives
The off-diagonal pattern is not a single set of disjoint blocks: the subject and item patterns overlap in different positions. For instance, observations 2 and 3 share neither classification, even though each shares one classification with observation 1. The result depends on the specified complete set of four subject-item combinations; the same construction works when some combinations are absent, using only the membership rows for observed responses.
Correlated Coefficients and General Incidence Rows
An incidence row need not contain only a single 1. For random slopes, a row contains the values of the predictors associated with that observation. If each group has a random intercept and slope, the row for observation \(i\) in group \(g\) contains \((1,x_i)\) in the columns for that group's coefficients. Let the coefficient covariance matrix be \(\Sigma\). The contribution to the covariance of observations \(i\) and \(j\) in the same group is \(z_i\Sigma z_j^{\mathsf T}\), where \(z_i=(1,x_i)\); observations in different groups have zero contribution from this term when group coefficients are independent. This allows the contribution to vary with predictor values and can be negative for some pairs if the covariance structure and predictor values imply it.
Worked Example: A Correlated Random Intercept and Slope
In one hypothetical group, consider two observations with predictor values \(x_1=0\) and \(x_2=2\). Assume the group's random intercept and slope have mean zero and covariance matrix \(\Sigma=\begin{pmatrix}1&0.2\\0.2&0.25\end{pmatrix}\). Assume the residual errors are independent across observations with variance 0.5. The random-effect rows are \(z_1=(1,0)\) and \(z_2=(1,2)\).
The random-effect variance for observation 1 is \(z_1\Sigma z_1^{\mathsf T}=1\). For observation 2 it is \(z_2\Sigma z_2^{\mathsf T}=1+2(2)(0.2)+2^2(0.25)=2.8\). Their random-effect covariance is \(z_1\Sigma z_2^{\mathsf T}=1+2(0.2)=1.4\). The residual variance adds 0.5 to each diagonal and adds nothing off the diagonal. The resulting marginal covariance matrix is
The larger variance for the second observation arises because its predictor value makes the random slope contribute to its variability. Its covariance with the first observation includes the intercept-slope covariance as well as the shared random intercept. Treating the term as a random intercept alone would miss both features.
Justification and Checks
The matrix sum follows from covariance rules. For independent terms, covariance of a sum is the sum of the covariances, and \(\operatorname{Cov}(Z_kb_k)=Z_k\operatorname{Cov}(b_k)Z_k^{\mathsf T}=Z_kG_kZ_k^{\mathsf T}\). For any vector \(a\) of length \(n\),
because \(G_k\), as a covariance matrix, is positive semidefinite. Therefore every contribution and their sum are positive semidefinite; adding a valid residual covariance \(R\) preserves this property. This is a useful check, not a guarantee that the membership structure or variance assumptions are correct.
Before interpreting a constructed matrix, check that its rows and columns use exactly the same observation order and that every \(Z_k\) has one row per observation. Check symmetry and confirm that each diagonal entry includes every random-effect variance applying to that observation and its residual variance. For random intercepts, verify each off-diagonal entry by listing the shared units. A negative eigenvalue beyond numerical rounding suggests an invalid covariance component, a calculation error, or an invalid parameter value.
The independence assumption between random-effect terms matters. If \(b_k\) and \(b_\ell\) are correlated, the covariance of the linear predictor also includes \(Z_k\operatorname{Cov}(b_k,b_\ell)Z_\ell^{\mathsf T}\) and its transpose. Omitting those terms is justified only when the cross-covariances are zero. Likewise, a residual structure that includes dependence must be represented in \(R\), rather than silently treated as independent residual errors.
GLMM Responses and the Limits of the Matrix Sum
For a GLMM, conditional on the random effects, responses are commonly assumed independent with a specified conditional distribution and mean \(\mu_i(b)\), linked to a linear predictor \(\eta_i=x_i^{\mathsf T}\beta+\sum_k z_{ik}^{\mathsf T}b_k\). The linear predictor still has covariance \(\sum_k Z_kG_kZ_k^{\mathsf T}\). But with a nonlinear link or a non-Gaussian response distribution, this is not generally the marginal covariance matrix of the observed responses.
The exact identity for two responses is \(\operatorname{Cov}(Y_i,Y_j)=E\{\operatorname{Cov}(Y_i,Y_j\mid b)\}+\operatorname{Cov}\{E(Y_i\mid b),E(Y_j\mid b)\}\). Under conditional independence, the first term is zero for \(i\ne j\), so dependence between distinct responses is the covariance of their conditional means across the random effects. On the diagonal, the conditional response variance also contributes. In general, this is a nonlinear calculation involving the response distribution and link, rather than simply adding \(Z_kG_kZ_k^{\mathsf T}\).
Worked Example: Shared Random Intercept in a Probit GLMM
Consider two hypothetical binary responses with conditional success probability \(\Phi(b)\), where \(b\sim N(0,1)\) is a shared random intercept and the responses are conditionally independent given \(b\). Here \(\Phi\) is the standard normal cumulative distribution function. Each marginal success probability is \(1/2\), by symmetry.
To find the joint success probability, introduce independent standard normal errors \(\epsilon_1,\epsilon_2\). A success can be represented by \(b+\epsilon_i>0\). Each sum has variance 2, and their correlation is \(\operatorname{Cov}(b+\epsilon_1,b+\epsilon_2)/2=1/2\). For two zero-mean standard normal variables with correlation \(\rho\), the probability both are positive is \(1/4+\arcsin(\rho)/(2\pi)\). Substituting \(\rho=1/2\) gives \(P(Y_1=1,Y_2=1)=1/4+(\pi/6)/(2\pi)=1/3\). Thus
The random intercept variance on the linear-predictor scale is 1, not \(1/12\). The latter is the response-scale covariance for this particular probit model and conditional mean. This example illustrates why a covariance matrix derived from incidence matrices must be identified with its scale: it may describe the linear predictor exactly without being the marginal covariance of the observed GLMM responses.
Choosing and Checking the Construction
A direct pairwise calculation is often the clearest way to verify a small covariance matrix: for each pair of observations, add the variances of the random-effect terms whose units they share. Incidence matrices are more systematic, especially with several terms or random slopes, and expose dimensions and ordering. A direct unstructured covariance model instead estimates covariance entries without tying them to membership; it may be useful when no meaningful structure is specified, but it does not encode the sharing assumptions of a nested or crossed random-effects model.
Membership structure determines which covariance patterns are allowed, not whether their variance components can be reliably estimated. If the design provides little information to distinguish components, a mathematically valid covariance matrix does not resolve that limitation. Also, zero covariance for a pair under the model means that no modeled random-effect term or residual structure links that pair; it is not proof that the observations are unrelated in every scientific sense.
A frequent mistake is to give every pair in the same school the same covariance, then add class covariance on top without checking which pairs share a class. The correction is to evaluate each term separately: school membership contributes to all within-school pairs, class membership contributes only to within-class pairs, and residual variance contributes to diagonal entries (under independent residuals). Another mistake is to report \(ZGZ^{\mathsf T}+R\) as the exact observed-response covariance for a binary or count GLMM. That expression is exact for the Gaussian identity-link model under its stated assumptions, and it gives the linear-predictor covariance; for non-Gaussian responses, use the law of total covariance on the response scale.
Check Your Understanding
Answer each question, then compare with the solutions below.
- Three observations are in one school: observations 1 and 2 share a class, while observation 3 is in another class. School, class, and residual variances are 3, 1, and 2. What are the diagonal entries and the distinct off-diagonal covariances of the Gaussian marginal covariance matrix?
- In a crossed subject-item model, what random-intercept covariance contribution applies to two observations that share both subject and item, if the subject and item variances are 0.8 and 1.2 and the terms are independent?
- Why does a random-intercept variance of 1 in a probit GLMM not imply that the observed binary responses have covariance 1?
- What does a negative eigenvalue in a proposed covariance matrix indicate, and what should be checked first?
Solutions. (1) Each diagonal is \(3+1+2=6\). Observations 1 and 2 share school and class, so their covariance is \(3+1=4\); observations 1 and 3, and 2 and 3, share only the school, so each covariance is 3. (2) Both independent terms contribute, so the covariance is \(0.8+1.2=2.0\); residual variance contributes only if the same response is being considered. (3) The variance is on the linear-predictor scale. The binary response covariance depends on the link and the distribution of conditional means; it is found from the law of total covariance, not by copying the random-effect variance. (4) A valid covariance matrix must be positive semidefinite. Check the observation ordering, incidence matrices, component covariance matrices, and arithmetic; allow for tiny negative values attributable only to numerical rounding.