Choosing Among Nested Covariance Structures

Why Compare Covariance Structures?

In GLS, the covariance matrix is part of the model: it determines how observations are weighted and how their dependence is represented. A diagonal covariance may allow unequal variances, while a structured matrix may also represent serial correlation or shared random effects. Choosing among such structures is not simply a matter of finding the most flexible one. Extra covariance parameters can improve fit, but may be weakly identified, unstable, or unnecessary for the estimation and inference at hand.

A useful comparison starts with a specified mean model and a set of plausible covariance models. This tutorial considers nested covariance structures, where one model is obtained from another by imposing restrictions on its covariance parameters. Likelihood-ratio tests and information criteria summarize different aspects of the comparison; neither removes the need to assess whether the candidate structures make sense for the data.

Definition: Let \(y\) have mean \(X\beta\) and covariance matrix \(\Sigma(\theta)\), where \(X\) is the specified design matrix and \(\theta\) denotes covariance parameters. A covariance model \(M_0\) is nested in \(M_1\) if every covariance matrix permitted by \(M_0\) is also permitted by \(M_1\). Equivalently, \(M_0\) is obtained from \(M_1\) by imposing restrictions on its parameters.

For example, a random-intercept model with nonnegative random-intercept variance includes the model with no random-intercept variation when that variance equals zero. Likewise, a covariance model with a freely estimated correlation may contain an independence model at correlation zero, provided zero is in the parameter space. Nesting concerns the set of permitted covariance matrices, not whether the two models have the same number of parameters or are equally plausible.

Likelihood-Ratio Comparisons

Suppose both models use the same response vector, mean specification, and likelihood basis. Let \(\ell_0\) and \(\ell_1\) be their maximized log likelihoods, with \(M_0\) the restricted model and \(M_1\) the larger model. Since \(M_1\) includes all parameter values allowed by \(M_0\), its maximized likelihood cannot be smaller. The likelihood-ratio statistic is

LR=2(\ell_1-\ell_0).

The null hypothesis says that the restrictions defining \(M_0\) hold; the alternative allows the additional covariance parameters in \(M_1\). If the null parameter value lies in the interior of the larger model’s parameter space and standard regularity conditions hold, Wilks’ theorem gives an approximate \(\chi^2_q\) reference distribution, where \(q\) is the number of independent restrictions. The test is an asymptotic procedure, so its accuracy depends on the sample information available for estimating the covariance parameters.

A key exception arises when the null places a variance component at zero. Variances cannot be negative, so zero is on the boundary of the allowed parameter space, not in its interior. The ordinary \(\chi^2_q\) rule need not apply. In the common case of testing one nonnegative variance component equal to zero, with other relevant parameters identifiable under the null and suitable regularity conditions, the asymptotic reference distribution is a mixture: half its mass is at zero and half follows \(\chi^2_1\). More complicated boundary restrictions can produce other mixtures or nonstandard distributions; the simple half-and-half rule should not be applied automatically.

Theorem: For a single nonnegative variance component tested at zero under the conditions above, the asymptotic null distribution of \(LR\) is \(\tfrac12\chi^2_0+\tfrac12\chi^2_1\), where \(\chi^2_0\) denotes a point mass at zero. For an observed \(LR>0\), the upper-tail p-value is \(\tfrac12P(\chi^2_1\geq LR)\). This boundary result replaces, rather than supplements, the ordinary \(\chi^2_1\) reference.

Worked Example: Testing a Variance Component at Zero

Suppose a Gaussian GLS analysis compares a model with no random-intercept variance to the nested model that adds one nonnegative random-intercept variance. The fitted log likelihoods are \(\ell_0=-50.0\) and \(\ell_1=-48.6\). Assume the mean structure is identical in both models, the larger model is correctly maximized, and the conditions for the single-component boundary mixture are appropriate.

The statistic is

$$ LR=2[-48.6-(-50.0)]=2(1.4)=2.80. $$

Because the null variance is on the boundary, use the mixture distribution. For \(LR=2.80\), its upper-tail probability is

p=\tfrac12P(\chi^2_1\geq 2.80) =\tfrac12\{2P(Z\geq\sqrt{2.80})\} =P(Z\geq 1.6733) \approx 0.04713,

where \(Z\sim N(0,1)\). At significance level \(\alpha=0.05\), this is evidence against the zero-variance restriction under the stated model and asymptotic approximation. It is not the probability that the null is true. The conclusion is also specific to this comparison: it does not establish that every feature of the larger covariance model is correct.

Information Criteria: Fit Balanced Against Complexity

Likelihood-ratio tests ask whether a specified restriction is inconsistent with the data under a reference distribution. Information criteria instead rank candidate models by combining fit with a penalty for estimated parameters. For a maximized log likelihood \(\ell\), \(k\) estimated parameters, and sample size \(n\), the usual forms are

AIC=-2\ell+2k, \qquad BIC=-2\ell+k\log(n).

For both criteria, smaller values are preferred among the models being compared. Count all estimated parameters included in the likelihood, including mean-model coefficients and covariance parameters. Use a consistent definition of the likelihood, parameter count, and sample size across candidates. AIC and BIC are not p-values, and a small difference does not prove that one covariance structure is true. AIC is commonly used to compare relative expected predictive information loss; BIC uses a stronger, sample-size-dependent penalty and has a large-sample model-selection interpretation under its assumptions.

With Gaussian data, maximum likelihood (ML) and restricted maximum likelihood (REML) are distinct likelihoods. For covariance comparisons with the same fixed-effect design \(X\), restricted likelihoods can be compared when calculated consistently. Do not compare REML likelihoods across models with different fixed-effect designs: the restricted likelihood depends on that design. ML is the standard basis when comparing models that differ in their fixed effects. Do not mix ML and REML values in a single comparison.

Worked Example: Ranking Three Covariance Models

Suppose three Gaussian models use the same response observations and fixed-effect design. Their maximized ML log likelihoods and total parameter counts are: \(M_0\), \(\ell=-120.4\), \(k=4\); \(M_1\), \(\ell=-116.9\), \(k=5\); and \(M_2\), \(\ell=-115.8\), \(k=7\). Let \(n=80\), and suppose the models are fitted on exactly the same observations.

For \(M_0\), \(AIC=-2(-120.4)+2(4)=240.8+8=248.8\). For \(M_1\), \(AIC=233.8+10=243.8\). For \(M_2\), \(AIC=231.6+14=245.6\). Thus \(M_1\) has the smallest AIC.

For BIC, \(\log(80)\approx4.3820\). The values are

\begin{aligned} BIC_{M_0}&=240.8+4(4.3820)\approx258.33,\\ BIC_{M_1}&=233.8+5(4.3820)\approx255.71,\\ BIC_{M_2}&=231.6+7(4.3820)\approx262.27. \end{aligned}

Here \(M_1\) also has the smallest BIC. The result is a ranking within these candidates using this sample and criterion; it is not a test that \(M_1\) is correct. \(M_2\) fits better in the likelihood sense than \(M_1\), but its additional parameters outweigh that improvement under both penalties in this example.

How the Covariance Choice Changes GLS

For a positive-definite covariance matrix \(\Sigma\), the GLS estimator is

\hat{\beta}_{GLS}=(X^\mathsf{T}\Sigma^{-1}X)^{-1}X^\mathsf{T}\Sigma^{-1}y.

The covariance matrix therefore affects both the coefficient estimate and its model-based covariance. When \(\Sigma\) is treated as known, \(\operatorname{Var}(\hat{\beta}_{GLS}\mid X)= (X^\mathsf{T}\Sigma^{-1}X)^{-1}\). In practice covariance parameters are generally estimated, so plug-in GLS uses an estimated matrix and inference must account for that estimation as appropriate. The formula makes clear why covariance selection is consequential: changing \(\Sigma\) changes the weighting of observations and the estimated precision.

Worked Example: Different Covariance Weights, Different GLS Estimate

Consider an intercept-only model for three hypothetical observations \(y=(1,2,5)^\mathsf{T}\), so \(X=(1,1,1)^\mathsf{T}\). To isolate the weighting effect, suppose the two covariance matrices below are treated as fixed. Under \(M_0\), let \(\Sigma_0=I_3\). Under \(M_1\), let \(\Sigma_1=\operatorname{diag}(1,1,4)\). This numerical illustration compares GLS calculations; it does not claim that either matrix was selected from these three observations.

Under \(M_0\), \(\Sigma_0^{-1}=I_3\), so the GLS estimate is the ordinary mean:

\hat{\beta}_0= \frac{1+2+5}{3} =\frac{8}{3}\approx2.667, \qquad \operatorname{Var}(\hat{\beta}_0\mid X)=\frac{1}{3}.

Under \(M_1\), the third observation has variance four times as large and receives one-quarter the weight of either other observation. Thus

\hat{\beta}_1 =\frac{1(1)+1(2)+\tfrac14(5)}{1+1+\tfrac14} =\frac{17/4}{9/4} =\frac{17}{9}\approx1.889, \qquad \operatorname{Var}(\hat{\beta}_1\mid X) =\frac{1}{1+1+\tfrac14} =\frac{4}{9}.

The estimate shifts because the value \(5\) is downweighted; the model-based standard error changes from \(\sqrt{1/3}\approx0.5774\) to \(\sqrt{4/9}=0.6667\), in the units of the response. This illustrates the mechanics, not a recommendation to choose a model solely because its coefficients or standard errors are more attractive. In fitted GLS, covariance estimation and uncertainty about the covariance parameters also matter.

A Short Justification of the GLS Formula

The GLS estimator minimizes the weighted residual sum of squares \(Q(\beta)=(y-X\beta)^\mathsf{T}\Sigma^{-1}(y-X\beta)\). Assuming \(X^\mathsf{T}\Sigma^{-1}X\) is invertible, differentiate with respect to \(\beta\), set the gradient to zero, and solve:

\begin{aligned} \frac{\partial Q(\beta)}{\partial\beta} &=-2X^\mathsf{T}\Sigma^{-1}(y-X\beta)=0,\\ X^\mathsf{T}\Sigma^{-1}X\hat{\beta} &=X^\mathsf{T}\Sigma^{-1}y,\\ \hat{\beta} &=(X^\mathsf{T}\Sigma^{-1}X)^{-1}X^\mathsf{T}\Sigma^{-1}y. \end{aligned}

If the mean is correctly specified and \(\operatorname{Var}(y\mid X)=\Sigma\), then \(\hat{\beta}-\beta=(X^\mathsf{T}\Sigma^{-1}X)^{-1}X^\mathsf{T}\Sigma^{-1}(y-X\beta)\). Taking its covariance gives \((X^\mathsf{T}\Sigma^{-1}X)^{-1}\). This result explains why a covariance model affects model-based inference as well as point estimates; under covariance misspecification, the model-based variance formula need not be valid.

Practical Judgment and Common Pitfalls

A structured comparison should begin with a defensible mean model and candidate covariance structures motivated by the observation process. Residual diagnostics can reveal patterns that a candidate model fails to capture, but a favorable likelihood or information criterion does not guarantee adequate diagnostics or robust conclusions. When models are close by an information criterion, stability of the GLS coefficients and their uncertainty across plausible structures is useful context.

Using an ordinary chi-square test at a boundary. A variance component equal to zero lies on the boundary. Check the parameter restrictions and use a justified mixture or other appropriate reference distribution.
Comparing incompatible likelihoods. The likelihoods must be calculated on the same observations and a consistent ML or REML basis. REML comparisons require the same fixed-effect design.
Counting only covariance parameters. For the stated AIC and BIC formulas, \(k\) includes all estimated parameters in the likelihood, not just the covariance parameters.
Treating the smallest criterion as proof. AIC and BIC rank the candidates under a penalty rule; neither establishes that the chosen model is true or that inference is insensitive to misspecification.
Assuming covariance choice only changes standard errors. GLS weights can change the coefficient estimates as well as their estimated covariance, as the worked calculation showed.

Key takeaway: For nested covariance models, compare maximized likelihoods only on a compatible basis, use a reference distribution appropriate to any boundary restrictions, and interpret information criteria as relative rankings. Then examine how plausible covariance choices affect the GLS estimates and their uncertainty.

Check Your Understanding

Answer each question using the definitions and comparison principles developed above.

In a model that adds one nonnegative variance component, why may an ordinary \(\chi^2_1\) likelihood-ratio reference be inappropriate when testing whether that component is zero?
Two nested models have maximized log likelihoods \(-42.0\) and \(-40.5\), with the larger model second. Calculate the likelihood-ratio statistic. What additional information is needed before assigning its p-value?
Why should REML likelihoods not be compared across models with different fixed-effect designs?
For \(n=100\), two models have \((\ell,k)=(-60,4)\) and \((-57,6)\). Calculate both AIC values and identify the smaller one.
In one sentence, explain how a covariance structure can change a GLS coefficient estimate even when the response values and design matrix are unchanged.