Plug-In Inference With Estimated Covariance Parameters

Why Plug-In Inference Needs Care

Generalized least squares (GLS) uses the error covariance structure to weight observations and account for their dependence. In applications, that structure is rarely known: its parameters are estimated from the data, and those estimates are inserted into the GLS calculation. This raises two questions. How should coefficient uncertainty be calculated after this substitution, and does estimating the covariance parameters itself add a first-order uncertainty term?

Under a correctly specified covariance model and suitable regularity conditions, the usual plug-in covariance matrix gives the first-order asymptotic uncertainty of feasible GLS coefficients. This does not make the formula an exact finite-sample result. The distinction matters when interpreting confidence intervals, especially in samples where the asymptotic approximation may be rough.

The Plug-In GLS Covariance Matrix

Consider the linear model \(y=X\beta+u\), where \(y\) is an \(n\)-vector of outcomes, \(X\) is a full-column-rank \(n\times p\) regressor matrix, and \(\beta\) is the \(p\)-vector of coefficients. Suppose, conditional on \(X\), the errors have mean zero and covariance \(\operatorname{Var}(u\mid X)=\Omega(\theta_0)\). Here \(\Omega(\theta_0)\) is positive definite and belongs to a specified covariance family indexed by a finite-dimensional parameter \(\theta_0\).

If \(\theta_0\) were known, GLS would weight using \(\Omega(\theta_0)^{-1}\). In feasible GLS, an estimate \(\hat{\theta}\) replaces \(\theta_0\). Write \(\hat{\Omega}=\Omega(\hat{\theta})\). The feasible GLS coefficient estimate and its conventional plug-in covariance estimate are

\hat{\beta}_{F} = (X'\hat{\Omega}^{-1}X)^{-1}X'\hat{\Omega}^{-1}y, \qquad \widehat{\operatorname{Var}}_{\text{plug-in}}(\hat{\beta}_{F}) = (X'\hat{\Omega}^{-1}X)^{-1}.

The diagonal entries of the second matrix estimate the coefficient variances; their square roots are the estimated standard errors. Off-diagonal entries estimate covariances between coefficient estimates. The formula assumes the covariance model is correctly specified and does not include a separate allowance for covariance misspecification.

Definition: A plug-in GLS covariance estimate is obtained by replacing the unknown covariance matrix in the known-covariance GLS variance formula with its fitted value: \((X'\hat{\Omega}^{-1}X)^{-1}\). It estimates the coefficient covariance under the fitted GLS model; it is generally an asymptotic approximation when covariance parameters are estimated.

Worked Examples of Plug-In Inference

Worked Example: An Intercept With Correlated Pairs

Suppose 20 hypothetical measurements are observed in 10 independent pairs, and the model has an intercept only. The fitted covariance matrix for each pair is \(\hat{\Omega}_{b}=\begin{pmatrix}2&1\\1&2\end{pmatrix}\), with zero covariance between different pairs. The 10 observed pairs are (11,13), (12,12), (10,14), (13,15), (12,14), (11,12), (14,13), (12,13), (13,14), and (10,10). Find the plug-in standard error and a 95% asymptotic confidence interval for the mean coefficient.

The pair covariance is positive definite: its eigenvalues are 1 and 3. Its inverse is \(\hat{\Omega}_{b}^{-1}=\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\). For an intercept-only model, each pair contributes \(1'\hat{\Omega}_{b}^{-1}1=2/3\) to \(X'\hat{\Omega}^{-1}X\). Thus the total information is \(10(2/3)=20/3\), and the plug-in variance is its reciprocal: \[ \widehat{\operatorname{Var}}_{\text{plug-in}}(\hat{\beta}_{F}) = (20/3)^{-1}=3/20=0.15. \] The inverse covariance gives equal weight to the two measurements within each pair, so the GLS intercept is their overall average. The pair sums total 248, giving \(\hat{\beta}_{F}=248/20=12.4\).

The standard error is \(\sqrt{0.15}=0.3873\). Using the 1.96 standard-normal critical value, the interval is \(12.4\pm1.96(0.3873)=12.4\pm0.7591\), or approximately \([11.641,13.159]\). Under the stated model and asymptotic approximation, this procedure gives a 95% confidence interval for the population intercept, in the same measurement units as the observations.

Worked Example: Estimating a Common Error Scale

In a hypothetical intercept-only regression with 25 independent observations, suppose the observed values have sample mean 10.0 and residual sum of squares 96. Assume normally distributed errors with covariance \(\Omega(\sigma^2)=\sigma^2 I\). Using the Gaussian maximum-likelihood estimate of the variance, calculate the plug-in standard error and 95% asymptotic interval for the intercept.

The maximum-likelihood estimate is \(\hat{\sigma}^2=96/25=3.84\). The fitted covariance is \(3.84I\), so the plug-in coefficient variance is \[ (X'\hat{\Omega}^{-1}X)^{-1} = \left(25/3.84\right)^{-1} = 3.84/25 = 0.1536. \] Therefore the standard error is \(\sqrt{0.1536}=0.3919\). The 95% asymptotic interval is \(10.0\pm1.96(0.3919)=10.0\pm0.7681\), or \([9.232,10.768]\) after rounding. It is an asymptotic interval based on the estimated variance, not an exact finite-sample interval.

Worked Example: A Contrast Uses the Full Covariance Matrix

Suppose a feasible GLS fit has two coefficients for standardized predictors, both measured in outcome points per predictor standard deviation. The estimates are \(\hat{\beta}_1=1.2\) and \(\hat{\beta}_2=-0.4\), and the plug-in covariance matrix is \(\begin{pmatrix}0.040&0.012\\0.012&0.090\end{pmatrix}\). Construct a 95% asymptotic interval for the contrast \(\beta_1+2\beta_2\).

The estimated contrast is \(1.2+2(-0.4)=0.4\). For \(c=(1,2)'\), its estimated variance is \(c'\widehat{\operatorname{Var}}_{\text{plug-in}}(\hat{\beta}_{F})c\), including the covariance between coefficients: \[ 0.040+2(1)(2)(0.012)+4(0.090)=0.448. \] The standard error is \(\sqrt{0.448}=0.6693\). The interval is \(0.4\pm1.96(0.6693)=0.4\pm1.312\), or approximately \([-0.912,1.712]\) outcome points per standardized-predictor contrast. This interval concerns the stated linear combination, not either coefficient separately.

Why Covariance-Parameter Estimation Is Often First-Order Negligible

The key result is that, under regularity conditions, feasible GLS has the same first-order limiting distribution as GLS that uses the true covariance matrix. In particular, assume the covariance family is correctly specified, \(\hat{\theta}\) is consistent, the inverse covariance matrix varies smoothly near \(\theta_0\), the regressors provide information growing proportionally to \(n\), and suitable moment and dependence conditions give the relevant laws of large numbers and central limit theorem. The number of coefficients and covariance parameters is fixed.

Theorem: Under these conditions, with \(A=\operatorname{plim}_{n\to\infty} n^{-1}X'\Omega(\theta_0)^{-1}X\) positive definite, \[ \sqrt{n}(\hat{\beta}_{F}-\beta)\ \xrightarrow{d}\ N(0,A^{-1}), \qquad n(X'\hat{\Omega}^{-1}X)^{-1}\ \xrightarrow{p}\ A^{-1}. \] Thus the plug-in covariance matrix consistently estimates the first-order covariance of feasible GLS. The result does not assert exact finite-sample coverage.

A short expansion explains why estimating \(\theta_0\) does not usually change the first-order coefficient distribution. For any covariance parameter value \(\theta\) near \(\theta_0\), define \(W(\theta)=\Omega(\theta)^{-1}\) and \(Q_n(\theta)=n^{-1}X'W(\theta)X\). Since \(y=X\beta+u\), the GLS estimate using \(\theta\) satisfies

\hat{\beta}(\theta)-\beta = Q_n(\theta)^{-1}\frac{X'W(\theta)u}{n}.

With \(Q_n(\theta)\) converging to a positive-definite limit and the stated moment and dependence conditions, the residual-weighted term \(X'W(\theta)u/n\) is \(O_p(n^{-1/2})\). Smoothness then makes the derivative of \(\hat{\beta}(\theta)\) with respect to \(\theta\) also \(O_p(n^{-1/2})\) locally: differentiating changes both \(Q_n(\theta)^{-1}\) and the weighted residual term, each multiplying a term of that order. A local expansion therefore gives

\hat{\beta}(\hat{\theta})-\hat{\beta}(\theta_0) = O_p(n^{-1/2})\,(\hat{\theta}-\theta_0) = o_p(n^{-1/2}),

because consistency means \(\hat{\theta}-\theta_0=o_p(1)\). The difference is smaller than the \(n^{-1/2}\) sampling scale of the coefficient estimate. Consequently, it disappears after multiplying by \(\sqrt{n}\), leaving the same first-order limit as GLS with known covariance. If \(\hat{\theta}\) is root-\(n\) consistent, this argument gives the stronger order \(O_p(n^{-1})\) for the difference under the stated smoothness conditions.

Confidence Intervals and Interpretation

For coefficient \(j\), let \(\hat{V}_{jj}\) be the \(j\)th diagonal entry of the plug-in covariance matrix. A two-sided \(100(1-\alpha)\%\) asymptotic confidence interval is

\hat{\beta}_{F,j} \ \pm\ z_{1-\alpha/2}\sqrt{\hat{V}_{jj}}, \qquad \hat{V}=(X'\hat{\Omega}^{-1}X)^{-1},

where \(z_{1-\alpha/2}\) is the corresponding standard-normal quantile. For example, \(z_{0.975}\approx1.96\) for a 95% interval. A confidence interval describes a procedure with its stated long-run coverage under the model and regularity conditions; it does not mean there is a 95% probability that a fixed coefficient lies in the particular computed interval.

The same principle applies to a fixed linear contrast \(c'\beta\): estimate it by \(c'\hat{\beta}_{F}\) and estimate its variance by \(c'\hat{V}c\). Keeping the off-diagonal entries is essential whenever coefficients are correlated. In applications, report the estimated coefficient, its standard error, interval, and units, and identify the covariance model on which the plug-in calculation depends.

Common Mistakes

Treating plug-in inference as exact. Substituting \(\hat{\Omega}\) does not make the known-covariance variance formula exact in finite samples. It provides a first-order asymptotic approximation under the stated assumptions.
Adding a nuisance-parameter variance term automatically. In correctly specified feasible GLS, consistent estimation of the covariance parameters is often first-order negligible for the coefficient estimator. That conclusion relies on the regularity conditions; it is not a universal rule for every two-step estimator.
Dropping covariances in a contrast. The variance of \(c'\hat{\beta}_{F}\) is \(c'\hat{V}c\), not just a sum of the coefficient variances. Cross-products can increase or decrease the result.
Using the interval as a probability statement about the parameter. The 95% refers to the long-run coverage of the interval procedure, not a 95% probability that a fixed parameter is in the observed interval.

Key takeaway: For correctly specified feasible GLS, use \((X'\hat{\Omega}^{-1}X)^{-1}\) for first-order plug-in coefficient inference. Its justification is asymptotic equivalence to GLS with known covariance parameters, not an assertion of exact finite-sample coverage.

Check Your Understanding

Answer each question using the assumptions and formulas stated above.

In the plug-in GLS covariance formula, which unknown object is replaced by its fitted value?
Why does consistency of \(\hat{\theta}\), together with the stated smoothness and regularity conditions, make its effect on \(\hat{\beta}_{F}\) smaller than the \(n^{-1/2}\) scale?
What is the estimated variance of a linear contrast \(c'\hat{\beta}_{F}\), and why might the off-diagonal covariance entries matter?
Does a 95% plug-in asymptotic interval have exact 95% coverage in every finite sample? Explain.
If the covariance model is correctly specified and all stated conditions hold, what is the first-order limiting covariance of \(\sqrt{n}(\hat{\beta}_{F}-\beta)\)?