6.4 Weighted least squares

Suppose the errors are uncorrelated but have unequal variances, \[ \operatorname {Var}(\varepsilon _i) = \frac {\sigma ^2}{w_i} , \qquad i = 1, \dots , n , \] with known positive weights \(w_i\). Observations with large \(w_i\) are precise; those with small \(w_i\) are noisy. Ordinary least squares treats them all alike. It seems right to give the precise observations more say, and to minimise the weighted sum of squares \[ S_w(\boldsymbol \beta ) = \sum _{i=1}^n w_i\bigl (y_i - \mathbf {x}_i^\top \boldsymbol \beta \bigr )^2 = (\mathbf {y} - \mathbf {X}\boldsymbol \beta )^\top \mathbf {W}(\mathbf {y} - \mathbf {X}\boldsymbol \beta ) , \] where \(\mathbf {W} = \operatorname {diag}(w_1, \dots , w_n)\). In matrix form, \(\operatorname {Var}(\boldsymbol \varepsilon ) = \sigma ^2\mathbf {W}^{-1}\).

Theorem 6.4.1 (Weighted least squares). Under A1, A3 and \(\operatorname {Var}(\boldsymbol \varepsilon ) = \sigma ^2\mathbf {W}^{-1}\), with \(\mathbf {X}\) of full column rank:

(i): \(S_w\) is minimised by \[ \hat {\boldsymbol \beta }_W = (\mathbf {X}^\top \mathbf {W}\mathbf {X})^{-1}\mathbf {X}^\top \mathbf {W}\mathbf {y} ; \]
(ii): \(\hat {\boldsymbol \beta }_W\) is unbiased, with \(\operatorname {Var}(\hat {\boldsymbol \beta }_W) = \sigma ^2(\mathbf {X}^\top \mathbf {W}\mathbf {X})^{-1}\);
(iii): \(\hat {\boldsymbol \beta }_W\) has the smallest variance, for every coefficient, among all linear unbiased estimators of \(\boldsymbol \beta \);
(iv): \(\hat \sigma ^2_W = S_w(\hat {\boldsymbol \beta }_W)/(n - p)\) is unbiased for \(\sigma ^2\), and under normality all the tests and intervals of Section 3 hold with \(\mathbf {X}^\top \mathbf {W}\mathbf {X}\) in place of \(\mathbf {X}^\top \mathbf {X}\) and \(\hat \sigma ^2_W\) in place of \(MS_{\mathrm {Res}}\).

Proof. Let \(\mathbf {W}^{1/2} = \operatorname {diag}(\sqrt {w_1}, \dots , \sqrt {w_n})\), and transform the model by multiplying through by it: \[ \mathbf {z} = \mathbf {W}^{1/2}\mathbf {y}, \qquad \mathbf {B} = \mathbf {W}^{1/2}\mathbf {X}, \qquad \mathbf {z} = \mathbf {B}\boldsymbol \beta + \boldsymbol \delta , \quad \boldsymbol \delta = \mathbf {W}^{1/2}\boldsymbol \varepsilon . \] That is, \(z_i = \sqrt {w_i}\,y_i\), and each row of \(\mathbf {X}\) is multiplied by \(\sqrt {w_i}\). By Lemma 3.1.3, \(E(\boldsymbol \delta ) = \mathbf {0}\) and \[ \operatorname {Var}(\boldsymbol \delta ) = \mathbf {W}^{1/2}(\sigma ^2\mathbf {W}^{-1})\mathbf {W}^{1/2} = \sigma ^2\mathbf {I} . \] So the transformed model satisfies A1–A3, and \(\mathbf {B}\) has full column rank because \(\mathbf {W}^{1/2}\) is invertible. Everything proved for ordinary least squares applies to it.

(i) \(S_w(\boldsymbol \beta ) = \|\mathbf {W}^{1/2}(\mathbf {y} - \mathbf {X}\boldsymbol \beta )\|^2 = \|\mathbf {z} - \mathbf {B}\boldsymbol \beta \|^2\), the ordinary sum of squares of the transformed model, minimised by Theorem 3.3.1 at \((\mathbf {B}^\top \mathbf {B})^{-1}\mathbf {B}^\top \mathbf {z} = (\mathbf {X}^\top \mathbf {W}\mathbf {X})^{-1}\mathbf {X}^\top \mathbf {W}\mathbf {y}\).

(ii) By Theorem 3.5.1 applied to the transformed model, with \((\mathbf {B}^\top \mathbf {B})^{-1} = (\mathbf {X}^\top \mathbf {W}\mathbf {X})^{-1}\).

(iii) A linear estimator \(\mathbf {Ay}\) of \(\boldsymbol \beta \) is \(\mathbf {AW}^{-1/2}\mathbf {z}\), a linear function of \(\mathbf {z}\), and conversely; unbiasedness means the same in both. So the linear unbiased estimators are the same whether written in terms of \(\mathbf {y}\) or of \(\mathbf {z}\), and Theorem 3.5.2 for the transformed model says that its least-squares estimator, \(\hat {\boldsymbol \beta }_W\), is best among them.

(iv) \(S_w(\hat {\boldsymbol \beta }_W)\) is the residual sum of squares of the transformed model, so Theorems 3.6.1 and 3.7.2 apply to it directly, and the tests and intervals are those of the transformed model. □

The proof is also the practical recipe: multiply each \(y_i\) and each row of \(\mathbf {X}\), including the column of ones, by \(\sqrt {w_i}\), and fit by ordinary least squares without an intercept (the old intercept column has become \(\sqrt {w_i}\), which plays its part). Residual analysis should use the weighted residuals \(\sqrt {w_i}\,(y_i - \hat y_i)\), the residuals of the transformed model, which have constant variance if the weights are right.

Remark (Generalised least squares). The same argument works when the errors are correlated, \(\operatorname {Var}(\boldsymbol \varepsilon ) = \sigma ^2\mathbf {V}\) with \(\mathbf {V}\) any known positive definite matrix. By the spectral theorem \(\mathbf {V} = \mathbf {P}\boldsymbol \Lambda \mathbf {P}^\top \), and \(\mathbf {K} = \mathbf {P}\boldsymbol \Lambda ^{1/2}\mathbf {P}^\top \) is a symmetric matrix with \(\mathbf {K}^2 = \mathbf {V}\). Multiplying through by \(\mathbf {K}^{-1}\) gives a model with covariance \(\sigma ^2\mathbf {I}\), and the best linear unbiased estimator is \[ \hat {\boldsymbol \beta }_{\mathrm {GLS}} = (\mathbf {X}^\top \mathbf {V}^{-1}\mathbf {X})^{-1}\mathbf {X}^\top \mathbf {V}^{-1}\mathbf {y} . \] Weighted least squares is the case \(\mathbf {V} = \mathbf {W}^{-1}\). The case of errors correlated in time is taken up in the section on autocorrelated errors.

In simple linear regression the weighted estimates have the familiar form, with every mean and sum weighted.

Corollary 6.4.2. For \(y = \beta _0 + \beta _1 x + \varepsilon \) with weights \(w_i\), let \(\bar x_w = \sum w_i x_i/\sum w_i\) and \(\bar y_w = \sum w_i y_i/\sum w_i\). Then \[ \hat \beta _1 = \frac {\sum w_i(x_i - \bar x_w)(y_i - \bar y_w)}{\sum w_i(x_i - \bar x_w)^2} , \qquad \hat \beta _0 = \bar y_w - \hat \beta _1\bar x_w , \] and \(\operatorname {Var}(\hat \beta _1) = \sigma ^2\big /\sum w_i(x_i - \bar x_w)^2\).

Proof. Setting the partial derivatives of \(S_w = \sum w_i(y_i - \beta _0 - \beta _1 x_i)^2\) to zero gives \(\sum w_i(y_i - \beta _0 - \beta _1 x_i) = 0\), which is \(\beta _0 = \bar y_w - \beta _1\bar x_w\), and \(\sum w_i x_i(y_i - \beta _0 - \beta _1 x_i) = 0\). Substituting the first in the second, \(\sum w_i x_i\bigl [(y_i - \bar y_w) - \beta _1(x_i - \bar x_w)\bigr ] = 0\). Since \(\sum w_i(y_i - \bar y_w) = \sum w_i(x_i - \bar x_w) = 0\), \(x_i\) may be replaced by \(x_i - \bar x_w\), which gives \(\hat \beta _1\). The variance is the \((2,2)\) entry of \(\sigma ^2(\mathbf {X}^\top \mathbf {W}\mathbf {X})^{-1}\), found exactly as in Example 3.3.2 with every sum weighted. □

What happens if the weights are ignored

Ordinary least squares remains unbiased when the variances are unequal, since unbiasedness used only A1. But it is no longer the best linear unbiased estimator, and, more seriously, its usual standard errors are wrong. By Lemma 3.1.3, with \(\hat {\boldsymbol \beta } = \mathbf {C}\mathbf {X}^\top \mathbf {y}\), \[ \operatorname {Var}(\hat {\boldsymbol \beta }) = \sigma ^2\,\mathbf {C}\mathbf {X}^\top \mathbf {W}^{-1}\mathbf {X}\mathbf {C} , \] not \(\sigma ^2\mathbf {C}\), and \(MS_{\mathrm {Res}}\) no longer estimates a common \(\sigma ^2\). The \(t\)-tests and intervals computed as if the variance were constant can then be badly misleading.

Choosing the weights

Weighted least squares needs the weights, at least up to a common constant. They come from three sources.

(i): Known from the way the data arose. If \(y_i\) is the mean of \(n_i\) independent observations, each with variance \(\sigma ^2\), then \(\operatorname {Var}(y_i) = \sigma ^2/n_i\) and \(w_i = n_i\).
(ii): A known form of dependence on a regressor. If the standard deviation is proportional to \(x_i\), then \(w_i = 1/x_i^2\) (Problem 6.5.4). The form may be suggested by the funnel in a residual plot against \(x\).
(iii): Estimated. Fit by ordinary least squares, regress the absolute residuals (or the squared residuals) on \(x\) or on \(\hat y\) to estimate how the spread changes, take \(w_i\) as the reciprocal of the fitted variance, and refit. With replicates at each \(x\), the sample variance at each \(x\) can be used. Estimated weights are themselves uncertain, so small samples call for caution.

Theorem 6.4.3. Suppose that at each of \(m\) settings \(\mathbf {x}_1, \dots , \mathbf {x}_m\) there are \(n_i\) observations \(y_{i1}, \dots , y_{in_i}\), with means \(\bar y_i\). The ordinary least-squares fit to all \(n = \sum n_i\) observations is the same as the weighted least-squares fit to the \(m\) means with weights \(w_i = n_i\).

Proof. As in the proof of Theorem 4.6.1, for any \(\boldsymbol \beta \), \[ \sum _{i=1}^m\sum _{j=1}^{n_i}\bigl (y_{ij} - \mathbf {x}_i^\top \boldsymbol \beta \bigr )^2 = \sum _{i=1}^m\sum _{j=1}^{n_i}(y_{ij} - \bar y_i)^2 + \sum _{i=1}^m n_i\bigl (\bar y_i - \mathbf {x}_i^\top \boldsymbol \beta \bigr )^2 , \] the cross term vanishing because deviations from \(\bar y_i\) add to zero. The first sum does not involve \(\boldsymbol \beta \), so the left side and the weighted sum of squares on the right are minimised by the same \(\boldsymbol \beta \). □

So when only the group means survive, the individual-level fit can still be recovered exactly, provided the means are weighted by the group sizes.

Example 6.4.4. A process was run at five temperature settings \(x\), with \(1\), \(4\), \(6\), \(4\) and \(1\) runs, and only the mean response \(\bar y_i\) at each setting was kept:

Setting \(x\)	1	2	3	4	5

Runs \(n_i\)	1	4	6	4	1
Mean \(\bar y_i\)	2	7	8	9	14

Fit the line by weighted least squares, and compare with ordinary least squares on the five means.

Solution. Weighted. \(w_i = n_i\), so \(\sum w_i = 16\), \[ \bar x_w = \frac {1 + 8 + 18 + 16 + 5}{16} = 3, \qquad \bar y_w = \frac {2 + 28 + 48 + 36 + 14}{16} = 8 , \] \begin {align*} \textstyle \sum w_i(x_i - 3)^2 &= 1(4) + 4(1) + 0 + 4(1) + 1(4) = 16, \\ \textstyle \sum w_i(x_i - 3)\bar y_i &= 1(-2)(2) + 4(-1)(7) + 0 + 4(1)(9) + 1(2)(14) = 32 . \end {align*}

So \(\hat \beta _1 = 32/16 = 2\), \(\hat \beta _0 = 8 - 2 \times 3 = 2\), and \[ \hat y_W = 2 + 2x . \] The fitted values are \(4\), \(6\), \(8\), \(10\), \(12\) and the weighted residual sum of squares is \[ S_w = 1(-2)^2 + 4(1)^2 + 6(0)^2 + 4(-1)^2 + 1(2)^2 = 16 , \] so \(\hat \sigma ^2_W = 16/3 = 5.33\), and \(\operatorname {se}(\hat \beta _1) = \sqrt {5.33/16} = 0.577\). By Theorem 6.4.3, this is the line that would have been fitted to all sixteen runs.

Unweighted. \(\bar x = 3\), \(\bar y = 8\), \(S_{xx} = 10\) and \(S_{xy} = -2(2) - 1(7) + 0 + 1(9) + 2(14) = 26\), so \(\hat \beta _1 = 2.6\), \(\hat \beta _0 = 8 - 7.8 = 0.2\): \[ \hat y = 0.2 + 2.6x . \] Comparison. The two end settings, each a single run and so the least reliable, lie furthest from the trend of the well-replicated middle three. Ordinary least squares, giving them full weight, tilts the line towards them. Weighted least squares gives them a sixth of the weight of the middle setting, and follows the precise means. How much precision is at stake? For the slope, Corollary 6.4.2 gives \(\operatorname {Var}(\hat \beta _{1W}) = \sigma ^2/16 = 0.0625\sigma ^2\). The ordinary estimator is \(\hat \beta _1 = \sum (x_i - 3)\bar y_i/10\), with \(\operatorname {Var}(\bar y_i) = \sigma ^2/n_i\), so \[ \operatorname {Var}(\hat \beta _1) = \frac {\sigma ^2}{100}\left (\frac {4}{1} + \frac {1}{4} + 0 + \frac {1}{4} + \frac {4}{1}\right ) = 0.085\sigma ^2 . \] Ignoring the weights costs a quarter of the precision: the unweighted estimate behaves as if it had been based on \(0.0625/0.085 = 74\%\) of the data.

[next] [prev] [prev-tail] [front] [up]

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.