Marginal Covariance Decomposition

Why Decompose the Marginal Covariance?

Repeated biomarker measurements on the same subject tend to vary together. Some of that dependence can arise because subjects have different underlying trajectories; some can arise from residual variation around each subject’s trajectory. The marginal covariance decomposition makes these sources explicit. It is useful both for interpreting a fitted mixed model and for checking whether the covariance it implies is plausible for the measurement schedule.

The focus here is the covariance of the observed measurement vector after the subject-specific random effects have been averaged over. The earlier material on random intercepts, random slopes, and residual error distinguishes the roles of those terms; the matrix decomposition puts their contributions together.

Model and Assumptions

For subject $i$, collect the $m_i$ repeated biomarker measurements in the vector $Y_i$. A linear mixed model can be written as

Y_i = X_i\beta + Z_i b_i + \varepsilon_i .

Here $X_i$ is the $m_i$-by-$p$ fixed-effects design matrix, $\beta$ is the vector of fixed-effects coefficients, $Z_i$ is the $m_i$-by-$q$ random-effects design matrix, and $b_i$ is the $q$-vector of subject-specific random effects. The vector $\varepsilon_i$ contains the residual errors for the $m_i$ observations.

Definition: Suppose $E(b_i)=0$, $\operatorname{Var}(b_i)=D$, $E(\varepsilon_i)=0$, and $\operatorname{Var}(\varepsilon_i)=R_i$. Assume $b_i$ and $\varepsilon_i$ are independent. Then the marginal mean and covariance of $Y_i$, after averaging over the random effects and residual errors, are $$ E(Y_i)=X_i\beta,\qquad \operatorname{Var}(Y_i)=V_i=Z_i D Z_i^{\mathsf T}+R_i . $$

The matrix $D$ describes the covariance among a subject’s random effects. The matrix $R_i$ describes the covariance among that subject’s residual errors. Both are covariance matrices and therefore must be symmetric and positive semidefinite. Their entries have units determined by the effects they represent; for example, a biomarker variance has units of biomarker-units squared.

The independence assumption matters to the stated sum. If $b_i$ and $\varepsilon_i$ are dependent, the variance also includes their cross-covariance terms. The decomposition as written does not require normality: it follows from the stated means, variances, and independence. If $b_i$ and $\varepsilon_i$ are jointly Gaussian, then $Y_i$ is Gaussian with mean $X_i\beta$ and covariance $V_i$. In many basic random-intercept and random-slope examples, residuals are assumed independent with common variance, so $R_i=\sigma^2 I_{m_i}$; the general formula allows other specified residual covariance matrices.

Why the Covariance Is a Sum

The fixed-effects term $X_i\beta$ is constant when computing the covariance for a given design, so it contributes no variance. Set $A=Z_i b_i$. The centered random part of $Y_i$ is $A+\varepsilon_i$. Expanding its covariance gives

\begin{aligned} \operatorname{Var}(Y_i) &=E\left[(A+\varepsilon_i)(A+\varepsilon_i)^{\mathsf T}\right]\\ &=E(AA^{\mathsf T})+E(\varepsilon_i\varepsilon_i^{\mathsf T}) +E(A\varepsilon_i^{\mathsf T})+E(\varepsilon_i A^{\mathsf T}). \end{aligned}

Because $E(b_i)=0$ and $E(\varepsilon_i)=0$, the first two terms are $Z_iD Z_i^{\mathsf T}$ and $R_i$, respectively. Independence gives $E(A\varepsilon_i^{\mathsf T})=0$ and $E(\varepsilon_i A^{\mathsf T})=0$. Therefore the two variance contributions add:

\operatorname{Var}(Y_i)=Z_iD Z_i^{\mathsf T}+R_i .

This is a matrix identity, not a claim that the two sources can be separated by simply inspecting one observed subject’s measurements. It describes the covariance implied by the model across repetitions of the subject-generating process.

Worked Example: Random Intercept With Independent Residuals

Suppose three measurements are taken on each hypothetical subject. The model has a random intercept $b_{0i}$ with variance $4$ biomarker-units squared, and independent residual errors with variance $9$ biomarker-units squared. There are no random slopes.

The random-effects design matrix is a column of ones, $Z_i=(1,1,1)^{\mathsf T}$, and $D=(4)$. The residual covariance is $R_i=9I_3$. Thus

Z_i D Z_i^{\mathsf T} = 4\begin{pmatrix} 1&1&1\\ 1&1&1\\ 1&1&1 \end{pmatrix}, \qquad R_i= \begin{pmatrix} 9&0&0\\ 0&9&0\\ 0&0&9 \end{pmatrix}.

Adding the matrices gives

V_i= \begin{pmatrix} 13&4&4\\ 4&13&4\\ 4&4&13 \end{pmatrix}.

Each marginal measurement variance is $4+9=13$ biomarker-units squared. Each covariance between two different measurements is $4$ biomarker-units squared: the shared random intercept contributes $4$, while the independent residual errors contribute zero to off-diagonal entries. The corresponding correlation is $4/\sqrt{13\cdot13}=0.3077$. The model therefore implies positive, equal pairwise correlations, not independent marginal measurements.

Reading Contributions Entry by Entry

For measurements $j$ and $k$ on the same subject, let $z_{ij}^{\mathsf T}$ and $z_{ik}^{\mathsf T}$ be the corresponding rows of $Z_i$. The marginal covariance is

\operatorname{Cov}(Y_{ij},Y_{ik}) =z_{ij}^{\mathsf T}D z_{ik}+(R_i)_{jk}.

The first term is the covariance induced by shared random effects at the two observation times. The second is the residual covariance at those measurements. In particular, the marginal variance at time $j$ is $z_{ij}^{\mathsf T}D z_{ij}+(R_i)_{jj}$. Even if residuals are independent, so that the off-diagonal entries of $R_i$ are zero, the off-diagonal entries of $V_i$ need not be zero.

Worked Example: Random Intercept and Random Slope

Consider a hypothetical biomarker measured at times $0$, $1$, and $2$ years. Let the random intercept and random slope be $b_{0i}$ and $b_{1i}$, respectively. Use

D= \begin{pmatrix} 4&1\\ 1&1 \end{pmatrix}, \qquad R_i=4I_3, \qquad Z_i= \begin{pmatrix} 1&0\\ 1&1\\ 1&2 \end{pmatrix}.

The entries of $D$ represent an intercept variance of $4$, an intercept-slope covariance of $1$, and a slope variance of $1$, with time measured in years. This $D$ is a valid covariance matrix: its determinant is $4(1)-1^2=3>0$, and its diagonal entries are positive. The residual errors have variance $4$ biomarker-units squared and are independent across the three measurements.

First calculate the random-effects contribution. Multiplying $Z_iD Z_i^{\mathsf T}$, or using $4+t+s+ts$ for the random-effects covariance at times $t$ and $s$, gives

Z_iD Z_i^{\mathsf T} = \begin{pmatrix} 4&5&6\\ 5&7&9\\ 6&9&12 \end{pmatrix}.

For example, at times $1$ and $2$, the contribution is $4+1+2+(1)(2)=9$. Now add $R_i$, which adds $4$ to each diagonal entry and nothing to the off-diagonal entries:

V_i= \begin{pmatrix} 8&5&6\\ 5&11&9\\ 6&9&16 \end{pmatrix}.

Thus the marginal variances at times $0$, $1$, and $2$ are $8$, $11$, and $16$ biomarker-units squared. The covariance between the measurements at times $0$ and $2$ is $6$ biomarker-units squared. Their marginal correlation is $6/\sqrt{8\cdot16}=0.5303$. The increasing variances here arise from the random-slope contribution; the residual variance is constant.

Worked Example: Unequal Residual Variances

Suppose, in a separate hypothetical example, a subject has three measurements with a random-intercept variance of $2.25$ biomarker-units squared. Assume independent residual errors with standard deviations $1$, $2$, and $3$ biomarker units at the respective measurements. This setup has $Z_i=(1,1,1)^{\mathsf T}$, $D=(2.25)$, and

R_i= \begin{pmatrix} 1&0&0\\ 0&4&0\\ 0&0&9 \end{pmatrix}.

The random-intercept contribution is $2.25$ in every matrix entry. Adding the residual contribution gives

V_i= \begin{pmatrix} 3.25&2.25&2.25\\ 2.25&6.25&2.25\\ 2.25&2.25&11.25 \end{pmatrix}.

The covariance of the first and second measurements remains $2.25$ biomarker-units squared, because their residual errors are independent. Their correlation is $2.25/\sqrt{3.25\cdot6.25}=0.4992$. The random-intercept covariance is unchanged across pairs, but the correlations differ because the marginal variances differ. This illustrates why equal random-effects covariance contributions do not guarantee equal marginal correlations.

Assumptions, Checks, and Alternatives

Before interpreting $V_i$, check that the model’s design matrices match the actual observation times and that $D$ and $R_i$ describe the intended sources of variation. In the examples, random-effects and residual terms were centered and independent; the residual covariance was explicitly specified; and the random-effects covariance was valid. A fitted covariance matrix should also be symmetric and positive semidefinite. A materially negative implied variance or an impossible correlation signals an algebra, implementation, or model-specification problem.

Useful diagnostics compare the covariance pattern implied by the model with the pattern in the data, for example by examining residuals and empirical within-subject associations across time. Such comparisons are uncertain with few subjects or sparse repeated measurements. They can still reveal mismatches, such as correlations that appear to decline with time separation when the fitted random-intercept model implies the same covariance for every pair. Diagnostics indicate possible lack of fit; they do not by themselves identify which component should be changed.

A random-effects representation is one way to construct a marginal covariance. An alternative is to specify a covariance matrix directly for the repeated measurements. A direct specification can represent patterns not captured by a chosen random-effects design, but may involve many covariance parameters and may not offer the same subject-level interpretation. Conversely, a random-effects structure can be parsimonious and interpretable, but constrains the covariance to the form $Z_iD Z_i^{\mathsf T}+R_i$. The choice is substantive and should be guided by the measurement schedule and scientific purpose, not by the decomposition alone.

The decomposition also has limits. It gives the covariance implied by the chosen model; it does not establish that the model is correct, or that either component has a causal interpretation. With too few observations per subject, different choices of $D$ and $R_i$ can be difficult to distinguish from the data. Nor does this formula alone address dependence between subjects or missingness mechanisms. Those require additional modeling assumptions.

Common Pitfall

Common mistake: “The residual variance is the marginal variance of each measurement, so the random-effects variance only explains correlation.” This is wrong because the random effects contribute to the diagonal of $V_i$ as well as its off-diagonal entries. For a random intercept with variance $4$ and residual variance $9$, the marginal variance is $4+9=13$, not $9$. The correction is to calculate the full matrix $Z_iD Z_i^{\mathsf T}+R_i$, then read both its diagonal variances and off-diagonal covariances.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

In the random-intercept example with random-intercept variance $4$ and residual variance $9$, what is the covariance between two distinct measurements?
If $D$ is a valid covariance matrix and $R_i$ is positive semidefinite, why is $Z_iD Z_i^{\mathsf T}+R_i$ also positive semidefinite?
For independent residuals, does a zero off-diagonal entry in $R_i$ imply a zero off-diagonal entry in $V_i$? Explain.
In the unequal-residual-variance example, what is the marginal variance of the third measurement, and what contributes to it?

Solutions. (1) It is $4$ biomarker-units squared: the shared random intercept contributes $4$, and independent residuals contribute zero to this covariance. (2) For any vector $a$, $a^{\mathsf T}Z_iD Z_i^{\mathsf T}a=(Z_i^{\mathsf T}a)^{\mathsf T}D(Z_i^{\mathsf T}a)\geq0$, and $a^{\mathsf T}R_i a\geq0$; their sum is nonnegative. (3) No. The random-effects contribution $Z_iD Z_i^{\mathsf T}$ can create a nonzero covariance between measurements. (4) It is $2.25+9=11.25$ biomarker-units squared: $2.25$ from the random intercept and $9$ from the third residual variance.