Diagnostics for Structured Covariance Misspecification

Why Diagnose a Structured Covariance Fit?

A structured covariance model replaces an unrestricted covariance matrix with a more interpretable or estimable pattern: for example, zero covariance between certain blocks, zero covariance beyond a band, equal covariances at equal lags, or a Kronecker product. Those restrictions can stabilize estimation, but they can also hide dependence the data contain. Diagnostics ask whether the restrictions leave systematic patterns in the residuals.

Throughout, suppose $Y_1,\ldots,Y_n$ are independent observations of a $p$-dimensional vector, with $Y_r\sim N_p(\mu,\Sigma)$. A fitted model supplies a mean estimate $\hat\mu$ and a positive-definite covariance estimate $\hat\Sigma$. These diagnostics assess the covariance pattern conditional on the fitted mean model. If observations are not independent replicates, the appropriate likelihood and residual calculations must reflect their sampling design.

Definition: The residual for replicate $r$ is $e_r=Y_r-\hat\mu$. The empirical residual covariance is $S_e=n^{-1}\sum_{r=1}^n e_re_r^\mathsf{T}$. A standardized residual covariance discrepancy for coordinates $i,j$ is $$ D_{ij}=\frac{(S_e)_{ij}-\hat\Sigma_{ij}}{\sqrt{\hat\Sigma_{ii}\hat\Sigma_{jj}}}. $$ The denominator puts the discrepancy on a scale relative to the fitted marginal standard deviations. It is not, by itself, a standard error or a test statistic.

It is also useful to inspect the standardized residual products $e_{ri}e_{rj}/\sqrt{\hat\Sigma_{ii}\hat\Sigma_{jj}}$ across replicates, rather than only their average. A few unusual replicates can make a covariance summary look like a persistent pattern. Because the same data are used to estimate $\hat\mu$ and $\hat\Sigma$, residual products are generally dependent and their averages need not have the same distribution as if the parameters were known.

What Residual Covariances Should Look Like

A first reference point is the distinction between a known and an estimated mean. This matters even before allowing for the additional effects of fitting a covariance model.

Result: If $\mu$ is known and $S_\mu=n^{-1}\sum_{r=1}^n(Y_r-\mu)(Y_r-\mu)^\mathsf{T}$, then $E[S_\mu]=\Sigma$. If instead $\bar Y=n^{-1}\sum_rY_r$ is used and $S_{\bar Y}=n^{-1}\sum_r(Y_r-\bar Y)(Y_r-\bar Y)^\mathsf{T}$, then $E[S_{\bar Y}]=\frac{n-1}{n}\Sigma$.

To verify the second statement, write $Y_r-\bar Y=(Y_r-\mu)-(\bar Y-\mu)$. Expanding the outer products and summing gives $$ \sum_{r=1}^n(Y_r-\bar Y)(Y_r-\bar Y)^\mathsf{T} =\sum_{r=1}^n(Y_r-\mu)(Y_r-\mu)^\mathsf{T} -n(\bar Y-\mu)(\bar Y-\mu)^\mathsf{T}. $$ The expectation of the first sum is $n\Sigma$, while $\operatorname{Var}(\bar Y)=\Sigma/n$, so the expectation of the subtracted term is $\Sigma$. Dividing by $n$ proves the result. In particular, the factor $(n-1)/n$ is a finite-sample effect of estimating the common mean; it is not evidence of covariance misspecification. $\square$

Under a fitted structure, inspect whether the residual covariance discrepancies align with the restrictions. Under a block-diagonal model, persistent residual covariance between blocks is a warning. Under a banded model, covariance remaining systematically beyond the fitted bandwidth suggests the band may be too narrow. Under a Toeplitz model, covariances at the same lag should be similar, after accounting for sampling variability. These are pattern checks, not claims that each empirical entry must equal its model value.

Worked Example: A Residual Pattern Beyond a Banded Model

Suppose four ordered measurements have fitted marginal variances equal to 1, and the fitted covariance is banded with bandwidth 1. Thus, the fitted covariances at lags 2 and 3 are zero. The empirical residual covariance matrix from many replicates is $$ S_e= \begin{pmatrix} 1.00&0.42&0.24&0.18\\ 0.42&1.00&0.39&0.22\\ 0.24&0.39&1.00&0.41\\ 0.18&0.22&0.41&1.00 \end{pmatrix}. $$

Since the fitted marginal variances are all 1, standardizing a covariance discrepancy leaves its numerical value unchanged. At lag 2, the three empirical covariances are $0.24,0.39,0.22$, so the mean lag-2 discrepancy is $$ (0.24+0.39+0.22)/3=0.85/3\approx 0.283. $$ At lag 3, the one available covariance is $0.18$. Both are positive despite the model's zero-covariance restrictions. This is a pattern worth investigating: it suggests dependence outside the fitted band, though the values alone do not establish statistical significance.

The lag-1 entries, $0.42,0.39,0.41$, are comparatively similar. That observation is consistent with a common lag-1 covariance but does not validate the full model. In particular, the diagnostics should be interpreted alongside uncertainty, possible outliers, and whether the mean model has left systematic residuals.

Lag and Block Summaries

When coordinates have an order, individual covariance entries can be noisy. A lag summary pools entries that share the same separation. For standardized coordinates, define the empirical residual correlation at pair $i,j$ by $$ \hat r_{ij}=\frac{(S_e)_{ij}}{\sqrt{(S_e)_{ii}(S_e)_{jj}}}, $$ when both empirical variances are positive. For lag $h$, let $\mathcal P_h=\{(i,j):|i-j|=h\}$, and use the average of $\hat r_{ij}$ over this set. Compare that profile with the fitted covariance-implied correlations. A Toeplitz covariance has equal covariances at a given lag, and because its diagonal entries are all equal, it also has equal correlations at a given lag. For a stationary correlation pattern with constant marginal variance, the correlation profile itself is also constant within each lag.

Result: If $\Sigma$ is Toeplitz, so that $\Sigma_{ij}=\gamma_{|i-j|}$, then for any nonempty lag set $\mathcal P_h$, the average $ |\mathcal P_h|^{-1}\sum_{(i,j)\in\mathcal P_h}\Sigma_{ij}$ equals $\gamma_h$. Thus, systematic differences among empirical entries at the same lag are relevant to assessing the Toeplitz restriction.

This identity follows directly from the defining equality: every term in the sum is $\gamma_h$, so the sum is $|\mathcal P_h|\gamma_h$, and division by the number of terms gives $\gamma_h$. In data, entries fluctuate around their population values; the identity does not imply that observed lag averages are identical or independent across lags. $\square$

For a block structure, summarize residual covariances by block pair, not just by coordinate. For each pair of blocks, inspect the corresponding residual covariance submatrix or a scale-adjusted summary such as the average of its standardized entries. A block-diagonal model predicts zero population covariances across distinct blocks. If the fitted model allows a shared cross-block component, the relevant check is instead whether the observed submatrix departs from that specified pattern. Pooling can help reveal a weak repeated departure, but a single large entry may signal an influential replicate rather than a broad structural failure.

Worked Example: Checking a Block Restriction

Suppose measurements form two blocks, each with two coordinates. The fitted model sets every cross-block covariance to zero, and all fitted marginal variances are 4. The cross-block entries of the empirical residual covariance are $0.8, 0.4, 0.6,$ and $0.2$.

Each standardized discrepancy is the covariance divided by $\sqrt{4\cdot4}=4$. The four values are therefore $0.20,0.10,0.15,$ and $0.05$. Their average is $$ (0.20+0.10+0.15+0.05)/4=0.50/4=0.125. $$

All four estimates are positive, so the pattern is more suggestive of a small positive cross-block association than one isolated large entry would be. Yet the average is not a p-value: its variability depends on the number of independent replicates, the fitted mean and covariance, and the dependence among residual products. A parametric simulation under the fitted model can provide a reference distribution for such a chosen summary.

Separable Covariance and Residual Structure

A separable covariance over two indexed dimensions has the form $\Sigma=A\otimes B$, where $A$ and $B$ are covariance factors and $\otimes$ denotes the Kronecker product. Such a model imposes joint restrictions: for indices $a,a'$ in the first dimension and $b,b'$ in the second, the covariance is $A_{aa'}B_{bb'}$. A diagnostic can therefore examine whether covariance patterns across one dimension are proportional across levels of the other dimension, after accounting for scale.

One simple implication, when the relevant entries are nonzero, is that ratios across the first dimension do not depend on the second-dimension pair: $$ \frac{\Sigma_{(a,b),(a',b)}}{\Sigma_{(c,b),(c',b)}}=\frac{A_{aa'}}{A_{cc'}}, $$ for pairs where the denominator is nonzero and the indexing uses the same $B_{bb}$ factor. In practice, inspect standardized residual covariance displays organized by both indices, or compare observed covariance subpatterns with those implied by the fitted factors. A departure may indicate nonseparability, but can also arise from misspecified marginal variances or mean structure. Factor scale normalization handles representational nonidentifiability; it does not remove the model's substantive covariance restrictions.

Worked Example: Comparing Two Lag Profiles Across Groups

Consider two groups of repeated measurements. For each group, suppose the empirical residual covariances at lag 1 are $0.30$ and $0.15$, while the corresponding covariances at lag 2 are $0.20$ and $0.10$. Assume the coordinates are ordered and these entries are already on a common covariance scale.

The group-1 to group-2 ratio is $0.30/0.15=2$ at lag 1 and $0.20/0.10=2$ at lag 2. This is consistent with proportional covariance profiles across these two lags, one pattern a separable model can produce. It is not proof of separability: two lags are only a small part of the full covariance pattern, and empirical ratios are unstable when denominators are small.

If instead the second group's lag-2 covariance were $0.16$, its ratio would be $0.20/0.16=1.25$, unlike the lag-1 ratio of 2. This discrepancy would be a candidate nonseparability signal. Before interpreting it, check scale differences and the uncertainty of all four estimates.

Likelihood-Based Checks

Residual summaries are interpretable but may not combine all restrictions efficiently. A likelihood comparison provides a global check when the structured model is nested in an unrestricted covariance model. For independent Gaussian replicates with an unrestricted mean, let $S=n^{-1}\sum_r(Y_r-\bar Y)(Y_r-\bar Y)^\mathsf{T}$, and assume $S$ is positive definite. Its covariance maximum-likelihood estimate in the unrestricted model is $S$. If $\hat\Sigma_0$ is the structured model's fitted covariance, the likelihood-ratio statistic is $T=2\{\ell_{\mathrm{unres}}-\ell_{\mathrm{struct}}\}$.

Result: For the Gaussian model just specified, the likelihood-ratio statistic can be written $$ T=n\left\{\log|\hat\Sigma_0|-\log|S|+\operatorname{tr}(\hat\Sigma_0^{-1}S)-p\right\}. $$ If the structured covariance family is a regular, identifiable, smooth submodel of dimension $q$ inside the positive-definite covariance matrices, and the true covariance is an interior point, Wilks' theorem gives $T\xrightarrow{d}\chi^2_{p(p+1)/2-q}$ under the null as $n\to\infty$.

For the formula, the Gaussian log-likelihood, after maximizing over the common mean, is $$ \ell(\Sigma)=-\frac{np}{2}\log(2\pi)-\frac n2\log|\Sigma|-\frac n2\operatorname{tr}(\Sigma^{-1}S). $$ Subtract the structured log-likelihood from the unrestricted one and multiply by 2. The unrestricted fit substitutes $S$; its trace term is $\operatorname{tr}(S^{-1}S)=p$. The resulting expression is exactly the displayed statistic. Wilks' theorem applies under the stated regularity conditions; it is not a finite-sample identity. $\square$

Worked Example: Testing a Diagonal Covariance

For $n=100$ independent bivariate Gaussian observations, suppose the sample covariance is $$ S=\begin{pmatrix}1&0.40\\0.40&1\end{pmatrix}. $$ Compare a diagonal covariance model with the unrestricted model. The diagonal fit has $\hat\Sigma_0=\operatorname{diag}(1,1)$, while $|S|=1-(0.40)^2=0.84$, $|\hat\Sigma_0|=1$, and $\operatorname{tr}(\hat\Sigma_0^{-1}S)=2$.

Therefore $$ T=100\{\log(1)-\log(0.84)+2-2\}=-100\log(0.84)\approx17.44. $$ The unrestricted covariance has three free entries and the diagonal model has two, giving one degree of freedom under the regularity conditions. The asymptotic reference is $\chi^2_1$; the statistic is large relative to that reference and indicates evidence against zero covariance. This calculation is a likelihood check, not a substitute for inspecting where a more complicated structured model fails.

Interpreting Diagnostics Carefully

Likelihood calibration can fail at boundaries, under nonidentifiability, or when the model is not a regular smooth submodel. Overlapping restrictions require the dimension of the identifiable covariance family, not a naive count of named parameters. In small samples, even regular chi-square approximations can be poor. A parametric bootstrap offers a practical alternative: simulate independent data under the fitted structured model, refit both models to every simulated data set, and compare the observed likelihood statistic or residual summary with its simulated distribution. This procedure must repeat the fitting steps, since treating estimated parameters as known understates uncertainty.

A common mistake is to conclude that a model is adequate because each standardized residual covariance is small, or to reject it because one entry is large. The entries are jointly dependent, their uncertainty varies, and many inspected entries create opportunities for chance extremes. Use summaries chosen to match the restrictions, inspect the full pattern, and account for estimation. Also avoid diagnosing covariance structure before checking the mean: a systematic mean pattern can create apparent residual dependence.

Key takeaway: Match each diagnostic to the restriction being assessed: cross-block entries for block assumptions, off-band entries for banded assumptions, within-lag comparisons for Toeplitz assumptions, and cross-dimension proportionality for separability. Treat residual summaries as evidence to investigate, and use a properly calibrated likelihood or simulation check when formal inference is needed.

Check Your Understanding

Use these questions to test whether you can connect a covariance restriction to an appropriate diagnostic.

Why does the expected residual covariance using the sample mean and divisor $n$ differ from $\Sigma$?
What residual pattern would be concerning for a block-diagonal model?
What does a Toeplitz covariance assumption imply about entries at the same lag?
Why is a standardized residual covariance discrepancy not automatically a z-statistic?
In the likelihood-ratio result, what conditions are needed for the stated chi-square reference distribution?