Coefficient Covariance and Standard Errors in Weighted Least Squares

Why Coefficient Uncertainty Depends on the Weights

Weighted least squares gives more influence to observations with greater precision. That changes not only the fitted coefficients but also their uncertainty: the precision-weighted design determines how much information the data provide about each coefficient. This tutorial derives that covariance and shows how to turn it into standard errors, both when the Gaussian dispersion is known and when it must be estimated.

Throughout, take the design matrix as fixed. Let $Y$ be the $n$-vector of responses, $X$ an $n\times p$ design matrix of full column rank, and $W$ a known diagonal matrix of positive relative precision weights. The Gaussian model is

Y\sim N_n(X\beta,\ \sigma^2W^{-1}).

Thus the $i$th response has variance $\sigma^2/w_i$, where $w_i$ is the $i$th diagonal entry of $W$. The weights specify relative precision; $\sigma^2$ is the common unknown or known scale multiplying all variances. The weighted least-squares estimator is

\hat{\beta}=(X^\mathsf{T}WX)^{-1}X^\mathsf{T}WY.

The inverse exists because $W$ is positive definite and $X$ has full column rank. In particular, $X^\mathsf{T}WX$ is positive definite.

Definition: The covariance matrix of $\hat{\beta}$ is the $p\times p$ matrix $\operatorname{Var}(\hat{\beta})$. Its $j$th diagonal entry is the variance of $\hat{\beta}_j$; its $(j,k)$ entry is the covariance between $\hat{\beta}_j$ and $\hat{\beta}_k$. A coefficient’s standard error is the square root of its variance, or of an estimate of that variance when $\sigma^2$ is unknown.

Deriving the Coefficient Covariance

The coefficient estimator is a linear transformation of $Y$. For a fixed matrix $A$, $\operatorname{Var}(AY)=A\operatorname{Var}(Y)A^\mathsf{T}$. Applying this identity gives the covariance result directly.

Theorem: Under the fixed-design Gaussian model $Y\sim N_n(X\beta,\sigma^2W^{-1})$, with $W$ known positive diagonal and $X$ full column rank, $$ \operatorname{Var}(\hat{\beta})=\sigma^2(X^\mathsf{T}WX)^{-1}. $$

Proof. Set $A=(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W$, so that $\hat{\beta}=AY$. Since $X$, $W$, and $A$ are fixed, the linear-transformation rule for covariance applies:

\begin{aligned} \operatorname{Var}(\hat{\beta}) &=A\operatorname{Var}(Y)A^\mathsf{T}\\ &=(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W (\sigma^2W^{-1}) W X(X^\mathsf{T}WX)^{-1}\\ &=\sigma^2(X^\mathsf{T}WX)^{-1} X^\mathsf{T}WX (X^\mathsf{T}WX)^{-1}\\ &=\sigma^2(X^\mathsf{T}WX)^{-1}. \end{aligned}

Here $W W^{-1}W=W$, since $W$ is invertible and symmetric. This proves the stated covariance matrix. $\square$

The matrix contains more than one uncertainty measure. Its diagonal gives coefficient variances, while its off-diagonal entries describe how the coefficient estimates vary together. If dispersion $\sigma^2$ is known, the standard error for coefficient $j$ is therefore

\operatorname{SE}(\hat{\beta}_j) =\sqrt{\sigma^2\left[(X^\mathsf{T}WX)^{-1}\right]_{jj}}.

Worked Example: Covariance in a Weighted Straight-Line Model

Suppose $x=(0,1,2)$, the model has an intercept and slope, and the relative weights are $1,2,1$. Thus

X= \begin{pmatrix} 1&0\\ 1&1\\ 1&2 \end{pmatrix}, \qquad W= \begin{pmatrix} 1&0&0\\ 0&2&0\\ 0&0&1 \end{pmatrix}.

First calculate the weighted cross-product matrix and its inverse:

X^\mathsf{T}WX= \begin{pmatrix} 4&4\\ 4&6 \end{pmatrix}, \qquad (X^\mathsf{T}WX)^{-1} =\frac{1}{8} \begin{pmatrix} 6&-4\\ -4&4 \end{pmatrix} = \begin{pmatrix} 0.75&-0.50\\ -0.50&0.50 \end{pmatrix}.

The theorem gives $\operatorname{Var}(\hat{\beta})=\sigma^2 \begin{pmatrix}0.75&-0.50\\-0.50&0.50\end{pmatrix}$. For example, if $\sigma^2=4$, then the intercept variance is $4(0.75)=3$, and its standard error is $\sqrt{3}\approx1.732$. The slope variance is $4(0.50)=2$, and its standard error is $\sqrt{2}\approx1.414$. The covariance between the intercept and slope estimates is $4(-0.50)=-2$.

When Dispersion Is Unknown

In many Gaussian regressions, $\sigma^2$ is unknown. The usual unbiased residual estimate, when $n>p$, is the weighted residual sum of squares divided by the residual degrees of freedom:

\hat{\sigma}^2 =\frac{(Y-X\hat{\beta})^\mathsf{T}W(Y-X\hat{\beta})}{n-p}.

This is the degrees-of-freedom-corrected residual variance estimate. Substituting it for $\sigma^2$ in the covariance formula gives the estimated coefficient covariance matrix and estimated standard errors:

\widehat{\operatorname{Var}}(\hat{\beta}) =\hat{\sigma}^2(X^\mathsf{T}WX)^{-1}, \qquad \widehat{\operatorname{SE}}(\hat{\beta}_j) =\sqrt{\hat{\sigma}^2\left[(X^\mathsf{T}WX)^{-1}\right]_{jj}}.

The residual estimate has the stated scale under this model. To see why, define $Z=W^{1/2}Y$ and $B=W^{1/2}X$. Then $Z$ has covariance $\sigma^2 I_n$, and the weighted residual sum of squares is the ordinary squared residual length after projecting $Z$ onto the column space of $B$. The projection matrix is $H=B(B^\mathsf{T}B)^{-1}B^\mathsf{T}$. It is symmetric, idempotent, and has rank $p$. Consequently,

\begin{aligned} E\!\left[(Y-X\hat{\beta})^\mathsf{T}W(Y-X\hat{\beta})\right] &=E[Z^\mathsf{T}(I_n-H)Z]\\ &=\sigma^2\operatorname{tr}(I_n-H)\\ &=\sigma^2(n-p). \end{aligned}

The second equality follows by expanding the quadratic form around the mean of $Z$: its mean contribution vanishes because the mean lies in the column space of $B$, and the covariance contribution is $\sigma^2\operatorname{tr}(I_n-H)$. Since $H$ has rank $p$, its trace is $p$. Therefore $E[\hat{\sigma}^2]=\sigma^2$. This establishes the scale of the plug-in covariance estimate; it does not make $\hat{\sigma}^2$ known or remove the sampling uncertainty in that estimate.

Worked Example: Standard Errors With Unknown Dispersion

Consider an ordinary straight-line model with four observations at $x=0,1,2,3$, unit weights, and responses $y=(1,2,2,5)^\mathsf{T}$. Here $n=4$, $p=2$, and

X^\mathsf{T}X= \begin{pmatrix} 4&6\\ 6&14 \end{pmatrix}, \qquad (X^\mathsf{T}X)^{-1} =\frac{1}{20} \begin{pmatrix} 14&-6\\ -6&4 \end{pmatrix}.

The coefficient estimate is

\hat{\beta} =(X^\mathsf{T}X)^{-1}X^\mathsf{T}y =\frac{1}{20} \begin{pmatrix} 14&-6\\ -6&4 \end{pmatrix} \begin{pmatrix} 10\\21 \end{pmatrix} = \begin{pmatrix} 0.7\\1.2 \end{pmatrix}.

The fitted values are $0.7,1.9,3.1,4.3$; the residuals are $0.3,0.1,-1.1,0.7$. Thus the residual sum of squares is $0.3^2+0.1^2+(-1.1)^2+0.7^2=0.09+0.01+1.21+0.49=1.80$. With $n-p=2$, the dispersion estimate is $\hat{\sigma}^2=1.80/2=0.90$. Hence

\widehat{\operatorname{Var}}(\hat{\beta}) =0.90\cdot\frac{1}{20} \begin{pmatrix} 14&-6\\ -6&4 \end{pmatrix} = \begin{pmatrix} 0.63&-0.27\\ -0.27&0.18 \end{pmatrix}.

The estimated standard error of the intercept is $\sqrt{0.63}\approx0.794$; that of the slope is $\sqrt{0.18}\approx0.424$. The off-diagonal entry $-0.27$ is their estimated covariance, not another standard error.

Worked Example: A Weighted Intercept-Only Model

Suppose three independent responses follow a common-mean model with relative weights $1,2,3$. Then $X$ is a column of ones and $X^\mathsf{T}WX=1+2+3=6$. Therefore $\hat{\beta}$ is the weighted mean, and its variance is

\operatorname{Var}(\hat{\beta})=\frac{\sigma^2}{6}.

If $\sigma^2=9$, the coefficient variance is $9/6=1.5$, so its standard error is $\sqrt{1.5}\approx1.225$. For responses $2,4,5$, the estimate itself is $\hat{\beta}=(1\cdot2+2\cdot4+3\cdot5)/6=25/6\approx4.167$. The coefficient’s standard error describes uncertainty around this weighted-mean estimate; it is not the standard deviation of an individual response, which depends on that response’s weight.

Interpreting and Checking the Result

The inverse matrix $(X^\mathsf{T}WX)^{-1}$ captures how the design and the relative weights combine to determine coefficient precision. In general, more information in the weighted design reduces coefficient variances, but correlations among columns of $X$ can also affect them. The off-diagonal entries matter whenever uncertainty in combinations of coefficients is considered; the diagonal alone is not the whole covariance structure.

A useful check is the units. The entries of $(X^\mathsf{T}WX)^{-1}$ contribute the appropriate design scaling, and multiplying by $\sigma^2$ gives coefficient variances. Standard errors are in the units of their corresponding coefficients. Also check that $X^\mathsf{T}WX$ is invertible and that $n>p$ before using the residual estimate with denominator $n-p$.

Common Mistakes

Inverting the weights. The covariance is $\sigma^2(X^\mathsf{T}WX)^{-1}$, not $\sigma^2(X^\mathsf{T}W^{-1}X)^{-1}$. The model gives response covariance $\sigma^2W^{-1}$, while $W$ itself is the precision-weight matrix in the estimator.
Forgetting the dispersion factor. The inverse weighted cross-product matrix alone is not the coefficient covariance when $\sigma^2$ is not one.
Calling an estimated standard error an exact standard deviation. With unknown dispersion, replacing $\sigma^2$ by $\hat{\sigma}^2$ produces an estimated covariance and estimated standard errors.
Using the wrong residual degrees of freedom. The residual estimate here divides by $n-p$, where $p$ is the number of fitted coefficients, not by $n$.
Confusing a response standard deviation with a coefficient standard error. A response’s conditional variance is $\sigma^2/w_i$; a coefficient’s variance comes from the full matrix $\sigma^2(X^\mathsf{T}WX)^{-1}$.

Key takeaway: Under known relative precision weights, the coefficient covariance is $\sigma^2(X^\mathsf{T}WX)^{-1}$. When dispersion is unknown, estimate it from weighted residuals and multiply that estimate by the same inverse matrix; take square roots of the diagonal entries for coefficient standard errors.

Check Your Understanding

Use the model assumptions and formulas above to answer the following.

In $Y\sim N_n(X\beta,\sigma^2W^{-1})$, what is the variance of observation $i$ when its weight is $w_i$?
Which entries of $\sigma^2(X^\mathsf{T}WX)^{-1}$ give the coefficient variances, and which describe covariances between distinct coefficients?
If the estimated dispersion is $2$ and the relevant diagonal entry of $(X^\mathsf{T}WX)^{-1}$ is $0.25$, what is the estimated standard error of that coefficient?
Why does estimating dispersion require $n>p$ for the residual degrees-of-freedom formula?
Explain why a coefficient standard error is not the same as the standard deviation of an individual response.