A Note on Weighted Least Squares

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 115

It is well known that the APL domino function includes a routine to obtain least squares fits of a straight line to a set of observations. Specifically, suppose that the model to be fitted to a set of multivariate data containing \(n\) observations is

where the residuals \(e_i\) (\(i = 1,2,\ldots n\)) are assumed to be identically and independently distributed with an error variance which has to be estimated along with the b’s. Suppose further that the data is presented in the form of a matrix with \(n\) rows and \((p+2)\) columns, the leftmost of which gives the y values, and the remainder the x-values, starting with a column of 1’s to correspond to the fitting of the constant term \(b_0\). The vector of regression coefficients is then given by :

BS : R[;1]⌹0 1↓R

and the residuals by

ES : R+.ׯ1,BS R

The error variance is given by

EVAR : (+/(ES R)*2)÷1+-/⍴R

and the variances of the b’s by

BVAR : 1 1⍉(EVAR R)×IMS R

where IMS stands for “Inverse Matrix Squared”, a function which as \((X'X)^{-1}\) occurs repeatedly in analysing residuals. It is defined as

IMS : ⌹(⍉X)+.×X←0 1↓R

In practice there are situations in which the assumptions that the residuals are identically distributed is not tenable, e.g. if the population is in fact several subpopulations whose residual variances are proportional to one of the x values, or if the residual variances corresponding to particular x values are known as the result of repeated experimentation, as might happen say in a physics experiment. The appropriate analysis is then weighted least squares. Assume that the variance corresponding to the ith. observation is \(w_i\) times the error variance, where \((w_1, w_2,\ldots\, w_n)\) are the items of a known vector \(W\). The appropriate regression coefficients, error variance, and variances of the b’s are the BS, EVAR, and BVAR of the data matrix

(DIAG ÷W*.5)+.×R

where DIAG is a function which constructs a diagonal matrix with the weights \(w_1 \ldots w_n\) as the leading diagonal, e.g.

DIAG : ((2⍴⍴R)⍴R)×T∘.=T←⍳⍴R

To give a specific example, consider the data

as a candidate for a least squares fit

so that R is

 0.8   1    3
 1.2   1    5
 2.5   1    8
 4.9   1   12
 8.3   1   15
12.3   1   17
17.6   1   20

The coefficients \(b_0\) and \(b_1\) are given by

      BS R
¯4.018 0.9466

If, however, there is evidence that the spread of the residuals increases with \(x\) so that a plot of \(e\) vs \(x\) has the appearance

then weighted least squares is appropriate. Assuming weights proportional to \(x\) gives a regression line \(y = -2.562 + 0.8192x\). As the graph below shows the effect is to downgrade the importance of the last two points which are in the relatively more variable x region.

The more the weights are “strengthened”, e.g. by taking them to be the squares of the x values which gives \(y = -1.721 + 0.709\,x\), the more closely does the weighted squares line track the smaller value points.

If the model

is chosen, i.e. the data matrix is changed to

R,R[;3]*2

the fitted quadratics are

and

The effect of weighted least squares is to cause the regression parabola to approach closer to the lowest x-value point as the graph below illustrates.