Block-Diagonal GLS for Correlated Observation Groups

Why Block-Diagonal GLS?

Many regression data sets have a natural grouping: repeated measurements on people, observations from clinics, or records collected in separate time periods. Errors for observations in the same group may be correlated, while errors from different groups are independent. Ordinary least squares does not use that covariance pattern in its weighting. Generalized least squares (GLS) does: it weights residual combinations according to their covariance, while a block-diagonal structure lets us carry out that weighting group by group.

The key distinction is that independence between groups does not imply independence within groups. A correct GLS calculation preserves the within-group covariance and uses the zero cross-group covariances to simplify the matrix operations.

Model and Covariance Structure

Suppose there are \(G\) groups. Group \(g\) has \(n_g\) observations, response vector \(y_g\), and regressor matrix \(X_g\), with \(X_g\) having \(k\) columns. Assume the regression coefficient vector \(\beta\) is common across groups.

Definition: A grouped linear regression model is \(y_g=X_g\beta+u_g\), for \(g=1,\ldots,G\). Conditional on the regressor matrices, assume \(E(u_g\mid X)=0\), \(Var(u_g\mid X)=\Omega_g\), and errors from distinct groups are independent. Each \(\Omega_g\) is a known, symmetric, positive-definite \(n_g\times n_g\) covariance matrix. Stack the groups in the same order in \(y\), \(X\), and \(u\). Then \(y=X\beta+u\) and \(Var(u\mid X)=\Omega=\operatorname{diag}(\Omega_1,\ldots,\Omega_G)\).

Here \(y\) is the stacked response vector and \(X\) is the stacked regressor matrix. The zero off-diagonal blocks in \(\Omega\) represent zero covariance across groups; the matrices \(\Omega_g\) retain any covariance within a group. The independence assumption is stronger than merely having zero cross-group covariances, although the block-diagonal covariance follows from either independence or zero cross-group covariance.

The positive-definiteness assumption ensures that every nonzero residual vector has positive variance and that \(\Omega^{-1}\) exists. We also require \(X\) to have full column rank, so the regression coefficients are identified by the design matrix. In this tutorial, the covariance matrices are treated as known. If covariance parameters have to be estimated, their estimation adds a separate layer to the procedure.

Worked Example: Assembling and Inverting Group Covariances

Suppose a hypothetical data set has two independent groups, each with two observations. In response-squared units, their within-group covariance matrices are \(\Omega_1=\begin{pmatrix}2&1\\1&2\end{pmatrix}\) and \(\Omega_2=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\). Find the stacked covariance matrix and its inverse.

Because the groups are independent, the covariance matrix is block diagonal:

\Omega= \begin{pmatrix} 2&1&0&0\\ 1&2&0&0\\ 0&0&1&0.5\\ 0&0&0.5&1 \end{pmatrix}.

Invert each block separately. For the first block, the determinant is \(2(2)-1(1)=3\), so \(\Omega_1^{-1}=\frac{1}{3}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\). For the second, the determinant is \(1(1)-0.5(0.5)=0.75\), so \(\Omega_2^{-1}=\frac{1}{0.75}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix} =\begin{pmatrix}1.333&-0.6667\\-0.6667&1.333\end{pmatrix}\). Thus:

\Omega^{-1}=\operatorname{diag}(\Omega_1^{-1},\Omega_2^{-1}).

The inverse also has zero off-diagonal blocks. In particular, inverting the full matrix does not require mixing observations from different groups.

GLS Estimator and Its Derivation

With the covariance matrix specified, GLS chooses coefficients to minimize the residual quadratic form weighted by the inverse covariance. A residual direction with high variance receives less weight; a direction with low variance receives more. The relevant criterion is:

Q(b)=(y-Xb)'\Omega^{-1}(y-Xb) =\sum_{g=1}^G (y_g-X_gb)'\Omega_g^{-1}(y_g-X_gb).

The second equality is exactly where the block structure helps: the full quadratic form is a sum of independent group contributions. Differentiating with respect to \(b\), setting the derivative to zero, and solving gives the GLS estimator.

\frac{\partial Q(b)}{\partial b} =-2X'\Omega^{-1}(y-Xb)=0 \quad\Longrightarrow\quad \hat{\beta}_{GLS}=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y.

The same expression can be computed from group-level contributions, without forming the full stacked matrices:

\hat{\beta}_{GLS} = \left(\sum_{g=1}^G X_g'\Omega_g^{-1}X_g\right)^{-1} \left(\sum_{g=1}^G X_g'\Omega_g^{-1}y_g\right).

To verify the estimator's conditional covariance, write its estimation error as \(\hat{\beta}_{GLS}-\beta=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}u\). Using \(Var(u\mid X)=\Omega\) gives:

Var(\hat{\beta}_{GLS}\mid X) =(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}\Omega\Omega^{-1}X(X'\Omega^{-1}X)^{-1} =(X'\Omega^{-1}X)^{-1}.

This variance result uses the stated covariance model and conditional mean assumption. Under those assumptions, GLS is unbiased conditional on \(X\); normal errors are not needed for this algebraic result. If the errors are also conditionally normal, the estimator has a normal conditional distribution.

Worked Example: GLS Regression Across Three Groups

Consider hypothetical observations from three independent groups. Each group contains two records, with response \(y\) measured in score points and regressor \(x\) measured in study hours. Fit a common intercept and slope, so each row of \(X_g\) is \((1,x)\). Suppose every group has error covariance in score-points-squared units \(\Omega_g=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\). The observations are:

Group	\(x\) (hours)	\(y\) (points)
1	0	1
1	1	3
2	0	2
2	1	2
3	0	4
3	1	5

The assumptions match the setup: the specified covariance matrix is symmetric and positive definite because its determinant is \(1-0.25=0.75>0\); errors are assumed independent across groups and have mean zero conditional on the regressors. Its inverse is \(\Omega_g^{-1}=\frac{1}{0.75}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix} =\frac{4}{3}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix}\). For each group, \(X_g=\begin{pmatrix}1&0\\1&1\end{pmatrix}\), and direct multiplication gives:

X_g'\Omega_g^{-1}X_g = \begin{pmatrix}4/3&2/3\\2/3&4/3\end{pmatrix}.

Adding this contribution across the three groups gives \(\sum_gX_g'\Omega_g^{-1}X_g=\begin{pmatrix}4&2\\2&4\end{pmatrix}\). For the right-hand side, group 1 contributes \((8/3,10/3)'\), group 2 contributes \((8/3,4/3)'\), and group 3 contributes \((6,4)'\). Therefore:

\sum_gX_g'\Omega_g^{-1}y_g = \begin{pmatrix}34/3\\26/3\end{pmatrix}, \qquad \begin{pmatrix}4&2\\2&4\end{pmatrix}^{-1} =\frac{1}{12}\begin{pmatrix}4&-2\\-2&4\end{pmatrix}.

Multiplication yields \(\hat{\beta}_{GLS}=(7/3,1)'\), or an estimated intercept of \(2.333\) points and an estimated slope of \(1.000\) point per additional study hour. The slope describes the fitted difference in mean score associated with one more study hour, under this regression model; it is not, by itself, a causal effect.

Whitening and Efficient Computation

A second way to understand GLS is to transform the model so the errors have identity covariance. For each group, factor its covariance as \(\Omega_g=L_gL_g'\), where \(L_g\) is a nonsingular lower-triangular Cholesky factor. Premultiply the group equation by \(L_g^{-1}\):

L_g^{-1}y_g=L_g^{-1}X_g\beta+L_g^{-1}u_g, \qquad Var(L_g^{-1}u_g\mid X)=I_{n_g}.

Stacking these transformed equations and applying ordinary least squares gives the GLS estimator. This is an algebraic equivalence: the transformed least-squares criterion is the original inverse-covariance-weighted criterion. In practice, numerical software typically uses factorizations and linear solves rather than explicitly forming matrix inverses, which is generally more stable.

Worked Example: Checking a Groupwise Whitening Transformation

For one group in the preceding setup, use \(\Omega_g=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\), \(X_g=\begin{pmatrix}1&0\\1&1\end{pmatrix}\), and \(y_g=(1,3)'\). A Cholesky factor is \(L_g=\begin{pmatrix}1&0\\0.5&\sqrt{0.75}\end{pmatrix}\), since \(L_gL_g'=\Omega_g\). Its inverse is \(L_g^{-1}=\begin{pmatrix}1&0\\-1/\sqrt{3}&2/\sqrt{3}\end{pmatrix}\).

For the candidate coefficient vector \(b=(2,1)'\), the original residual is \(r=y_g-X_gb=(-1,0)'\). The GLS contribution to the objective is:

r'\Omega_g^{-1}r = \begin{pmatrix}-1&0\end{pmatrix} \frac{4}{3}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix} \begin{pmatrix}-1\\0\end{pmatrix} =\frac{4}{3}.

The transformed residual is \(L_g^{-1}r=(-1,1/\sqrt{3})'\), whose squared length is \((-1)^2+(1/\sqrt{3})^2=1+1/3=4/3\). Thus whitening produces the same weighted residual contribution. Applying this transformation to every group gives a computation that can use standard least squares on the transformed data.

What the Block Structure Buys

When the number of groups is large, block structure can make both computation and interpretation clearer. The estimator needs only the groupwise matrices \(X_g'\Omega_g^{-1}X_g\) and vectors \(X_g'\Omega_g^{-1}y_g\), which are added across groups. There is no need to calculate or store cross-group covariance terms under the model. Groups may also have different sizes and different covariance matrices, provided each block has the stated dimensions and is positive definite.

A common pitfall is to reason that independent groups mean all observations are independent, and therefore to replace every \(\Omega_g\) by a diagonal matrix. That drops within-group dependence and changes the GLS weights. The correction is to represent independence across groups with zero off-diagonal blocks, while retaining the specified off-diagonal entries inside each \(\Omega_g\). A related practical check is to use \(\Omega_g^{-1}\), not \(\Omega_g\), in the GLS criterion.

Key takeaway: Block-diagonal GLS is ordinary GLS with the covariance structure organized by group. The estimator is obtained by summing each group's weighted design cross-product and weighted response cross-product, then solving the resulting normal equations.

Check Your Understanding

Use the model and assumptions in this tutorial to answer the following.

What does a zero off-diagonal block in the stacked covariance matrix say about errors in two different groups?
If a group has covariance matrix \(\Omega_g\), which matrix appears in that group's GLS quadratic form?
Why can the stacked GLS estimator be calculated from sums of groupwise matrices and vectors?
What conditions on \(\Omega_g\) ensure that its inverse and a Cholesky factor are available?
In the regression example, what are the units of the slope, and what association does it describe?