Generalized Least Squares for Correlated Cases

Why Correlated Cases Change Least Squares

Ordinary least squares treats the case errors as having the same variance and being uncorrelated across cases. When cases are correlated, or have unequal variances, that geometry no longer matches the error structure: directions with more error variation should count less, and correlated directions should not be counted as independent information. Generalized least squares (GLS) incorporates this structure through a weight matrix. Here the row covariance matrix is known and positive definite, so the weighting can be derived directly.

The focus is the coefficient estimator for a multivariate response. The response covariance across outcomes is kept distinct from the covariance across cases. This distinction matters: the GLS weighting in this model acts on the rows of the data matrices, while the outcomes remain columns.

Model and Generalized Least Squares

Definition: Let \(Y\) be an \(n\times p\) response matrix, \(X\) an \(n\times q\) design matrix of full column rank, and \(B\) a \(q\times p\) coefficient matrix. The multivariate linear model is \[ Y=XB+E. \] Assume \(E\) has mean zero and covariance structure \(\operatorname{Cov}(\operatorname{vec}E)=\Sigma\otimes\Omega\), where \(\Omega\) is a known positive-definite \(n\times n\) row covariance matrix and \(\Sigma\) is a positive-definite \(p\times p\) response covariance matrix. The vectorization stacks columns. GLS minimizes \[ \operatorname{tr}\!\left\{\Sigma^{-1}(Y-XB)'\Omega^{-1}(Y-XB)\right\} \] over \(B\). Since \(\Sigma\) is positive definite, this minimizer is the same as the minimizer of \(\operatorname{tr}\{(Y-XB)'\Omega^{-1}(Y-XB)\}\).

The separable covariance assumption says that case dependence is represented by \(\Omega\) and response dependence by \(\Sigma\). The estimator below requires only the specified second moments, not normality. If the errors are matrix normal, the estimator also has a matrix-normal sampling distribution.

Whitening the Model and Deriving the Estimator

Because \(\Omega\) is positive definite, it has a symmetric positive-definite inverse square root \(\Omega^{-1/2}\), satisfying \(\Omega^{-1/2}\Omega\Omega^{-1/2}=I_n\). Premultiply the model by this matrix:

Y_w=\Omega^{-1/2}Y,\qquad X_w=\Omega^{-1/2}X,\qquad Y_w=X_wB+E_w,\quad E_w=\Omega^{-1/2}E.

The transformed errors have identity row covariance. In fact, using column-stacking vectorization, \(\operatorname{Cov}(\operatorname{vec}E_w)=\Sigma\otimes(\Omega^{-1/2}\Omega\Omega^{-1/2})=\Sigma\otimes I_n\). Thus GLS in the original model is ordinary least squares in the whitened model.

Theorem: Under the model conditions above, the unique GLS estimator is \[ \hat B_{\mathrm{GLS}}=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}Y. \] Its expectation is \(B\), and its covariance is \[ \operatorname{Cov}(\operatorname{vec}\hat B_{\mathrm{GLS}}) =\Sigma\otimes(X'\Omega^{-1}X)^{-1}. \] If \(E\) is matrix normal, then \[ \hat B_{\mathrm{GLS}}\sim MN_{q,p}\!\left(B,\,(X'\Omega^{-1}X)^{-1},\,\Sigma\right). \]

Derivation of the estimator. Write \(A=\Omega^{-1/2}\), so \(X_w=AX\) and \(Y_w=AY\). Since \(X\) has full column rank and \(A\) is invertible, \(X_w\) has full column rank. The ordinary least-squares normal equation for the whitened model is \(X_w'X_w\hat B=X_w'Y_w\). Its coefficient matrix is \(X_w'X_w=X'\Omega^{-1}X\), which is positive definite and hence invertible. Also \(X_w'Y_w=X'\Omega^{-1}Y\). Solving gives the stated estimator. This is the unique minimizer because the quadratic criterion has positive-definite coefficient matrix \(X'\Omega^{-1}X\).

Proof of unbiasedness and covariance. Define \(C=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}\), so that \(\hat B_{\mathrm{GLS}}=CY\). Since \(CX=I_q\), \[ E[\hat B_{\mathrm{GLS}}]=C(XB+E)=B. \] For the covariance, first calculate \[ C\Omega C' =(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}\Omega\Omega^{-1}X(X'\Omega^{-1}X)^{-1} =(X'\Omega^{-1}X)^{-1}. \] Further, \(\operatorname{vec}(CE)=(I_p\otimes C)\operatorname{vec}(E)\). Therefore \[ \begin{aligned} \operatorname{Cov}(\operatorname{vec}(CE)) &=(I_p\otimes C)(\Sigma\otimes\Omega)(I_p\otimes C')\\ &=\Sigma\otimes(C\Omega C')\\ &=\Sigma\otimes(X'\Omega^{-1}X)^{-1}. \end{aligned} \] Since \(\hat B_{\mathrm{GLS}}-B=CE\), this is the covariance of the estimator. Under matrix-normal errors, a linear transformation of \(E\) remains matrix normal, with the stated mean and row and response covariance matrices. \(\square\)

Worked Calculations

Worked Example: Diagonal Case Weights and the Corrected Cross-Products

Take a scalar response with \(n=3\), design rows \((1,0),(1,1),(1,2)\), response vector \(Y=(1,2,2)'\), and \(\Omega=\operatorname{diag}(1,4,1)\). Thus \(W=\Omega^{-1}=\operatorname{diag}(1,1/4,1)\). Calculate \(\hat B_{\mathrm{GLS}}\).

The design and response are \[ X=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix},\qquad Y=\begin{pmatrix}1\\2\\2\end{pmatrix}. \] The weighted cross-product is \[ X'WX= \begin{pmatrix} 1+1/4+1&1/4+2\\ 1/4+2&1/4+4 \end{pmatrix} = \begin{pmatrix}9/4&9/4\\9/4&17/4\end{pmatrix}. \] Its determinant is \((9/4)(17/4)-(9/4)^2=153/16-81/16=9/2\). Hence \[ (X'WX)^{-1}=\frac{2}{9} \begin{pmatrix}17/4&-9/4\\-9/4&9/4\end{pmatrix} = \begin{pmatrix}17/18&-1/2\\-1/2&1/2\end{pmatrix}. \] The weighted response product is \[ X'WY= \begin{pmatrix}1+2/4+2\\2/4+2(2)\end{pmatrix} = \begin{pmatrix}7/2\\9/2\end{pmatrix}. \] Consequently, \[ \hat B_{\mathrm{GLS}} = \begin{pmatrix}17/18&-1/2\\-1/2&1/2\end{pmatrix} \begin{pmatrix}7/2\\9/2\end{pmatrix} = \begin{pmatrix}19/18\\1/2\end{pmatrix}. \] The fitted intercept is \(19/18\) and the slope is \(1/2\). The calculation illustrates why each entry of \(X'WX\) must be recomputed from the given \(X\) and \(W\); an incorrect cross-product changes the estimator.

Worked Example: A Non-Diagonal Case Covariance

Suppose \(Y=(1,3)'\), \(X=(1,1)'\), and \[ \Omega=\begin{pmatrix}1&1/2\\1/2&1\end{pmatrix}. \] This matrix is positive definite because its eigenvalues are \(3/2\) and \(1/2\). Find the GLS estimate of the common mean.

The inverse is \[ \Omega^{-1}=\frac{4}{3}\begin{pmatrix}1&-1/2\\-1/2&1\end{pmatrix} =\begin{pmatrix}4/3&-2/3\\-2/3&4/3\end{pmatrix}. \] Then \(\Omega^{-1}X=(2/3,2/3)'\), so \(X'\Omega^{-1}X=4/3\) and \(X'\Omega^{-1}Y=8/3\). Therefore \[ \hat B_{\mathrm{GLS}}=(4/3)^{-1}(8/3)=2. \] Here the GLS estimate happens to equal the ordinary average. The equality follows because \(\Omega^{-1}\) has equal row sums, so \(\Omega^{-1}\mathbf 1\) is proportional to \(\mathbf 1\). If the scalar error covariance is \(\sigma^2\Omega\), the estimator variance is \(\sigma^2(X'\Omega^{-1}X)^{-1}=3\sigma^2/4\).

Worked Example: Two Responses and Their Coefficient Covariance

Let \(X=(1,1)'\), \(\Omega=\operatorname{diag}(1,4)\), and \[ Y=\begin{pmatrix}1&2\\5&4\end{pmatrix},\qquad \Sigma=\begin{pmatrix}4&1\\1&9\end{pmatrix}. \] Find the GLS coefficient row and its covariance.

The weight matrix is \(\Omega^{-1}=\operatorname{diag}(1,1/4)\), giving \(X'\Omega^{-1}X=1+1/4=5/4\). The weighted response product is \[ X'\Omega^{-1}Y=(1,1/4) \begin{pmatrix}1&2\\5&4\end{pmatrix} =(9/4,3). \] Thus \(\hat B_{\mathrm{GLS}}=(4/5)(9/4,3)=(9/5,12/5)\). The row covariance factor is \((X'\Omega^{-1}X)^{-1}=4/5\), so the covariance matrix of the two estimated coefficients is \[ \frac45\Sigma= \begin{pmatrix}16/5&4/5\\4/5&36/5\end{pmatrix}. \] The off-diagonal entries retain the response correlation represented by \(\Sigma\); the case weights determine the scalar row factor \(4/5\).

Interpretation and a Common Pitfall

The matrix \(\Omega^{-1}\) determines how case directions contribute to the fit. If \(\Omega\) is diagonal, its diagonal entries encode unequal case variances and GLS uses their reciprocals as weights. With non-diagonal \(\Omega\), the weighting also accounts for correlations between cases; in general, it cannot be reproduced by applying separate scalar weights to each observation. Whitening makes this precise: it transforms the error covariance across cases to the identity before fitting.

Pitfall: treating \(\Omega\) as the response covariance. In this model \(\Omega\) is \(n\times n\), acts on rows, and describes dependence across cases. \(\Sigma\) is \(p\times p\), acts on columns, and describes dependence across responses.
Pitfall: using \(\Omega\) instead of \(\Omega^{-1}\) in the criterion. GLS minimizes the residual quadratic form weighted by \(\Omega^{-1}\). A direction with larger error variance is downweighted, not upweighted.
Pitfall: assuming GLS must change every coefficient. Special designs or covariance structures can make some GLS estimates equal to their ordinary least-squares counterparts. Equality in one example does not make the two methods generally identical.

The full-column-rank assumption ensures that the displayed inverse exists and the coefficient matrix is uniquely determined. The next issue is what can be estimated when that rank condition fails; it is not needed for the derivation here.

Key takeaway: For known positive-definite case covariance \(\Omega\), whitening by \(\Omega^{-1/2}\) converts the model to ordinary least squares. The resulting coefficient estimate is \((X'\Omega^{-1}X)^{-1}X'\Omega^{-1}Y\), with covariance \(\Sigma\otimes(X'\Omega^{-1}X)^{-1}\).

Check Your Understanding

Use the model dimensions and covariance roles carefully when answering.

Why does positive definiteness of \(\Omega\) allow the model to be whitened?
In \(Y=XB+E\), which covariance matrix describes dependence across cases, and which describes dependence across responses?
For scalar-response errors with covariance \(\sigma^2\Omega\), write the variance of \(\hat B_{\mathrm{GLS}}\).
Explain why the GLS objective uses \(\Omega^{-1}\), rather than \(\Omega\).
In the two-response example, identify what determines the off-diagonal covariance between the two estimated coefficients.