Could anyone tell me why I am getting different results from R weighted least squares and manual solution by matrix operation?

Specifically, I am trying to manually solve $\mathbf W \mathbf A\mathbf x=\mathbf W \mathbf b$, where $\mathbf W$ is the diagonal matrix on weights, $\mathbf A$ is the data matrix, $\mathbf b$ is the response vector.

I am trying to compare the results with the R lm function using the weights argument.

enter image description here

  • $\begingroup$ I edited tags: this was definitely not [self-study]. It's also not really about GLS (but about a very special case), so I removed that one as well. $\endgroup$ – amoeba Oct 6 '16 at 0:02

As you can see from the mathematical expressions for your calculations, you are obtaining

$$((WA)^\prime (WA))^{-1} \; ((WA)^\prime (Wb)) = (A^\prime W^2 A)^{-1} (A^\prime W^2 b).$$

Evidently your weights are $W^2$, not $W$. Thus you should be comparing your answer to the output of

> lm(form, mtcars, weights=w^2)
      wt        hp      disp  
14.12980   0.08391  -0.16446 

The agreement is perfect (to within floating point error--internally, R uses a numerically more stable algorithm.)

  • 1
    $\begingroup$ Arguably, we're just talking about software conventions here: where the software expects "weights," does it want you to give it $W$ or $W^2$? I thought this was a valuable question because the problem could affect any statistical package. Regardless of the conventions, the brief analysis in this answer suggests what alternative interpretations of "weights" might be reasonable and worth experimenting with in any circumstance. $\endgroup$ – whuber Jun 9 '16 at 20:25
  • $\begingroup$ Yes, I think it is confusing, I got the expression from Gilbert Strang's linear algebra book Chapter 8.6, where he says weighted least square is just an adjustment from $Ax=b$ to $WAx=Wb$ $\endgroup$ – Haitao Du Jun 9 '16 at 20:28
  • 8
    $\begingroup$ Strang is correct, but he has the pedagogical orientation backwards: he starts with the answer rather than the problem. The problem concerns how to perform the analog of a least-squares procedure when the variances of the residuals have known, but different, values. For various (but simple) theoretical reasons, the data should be weighted by the inverse variances (sometimes called the "precisions"). From that one can work out that $W$ must be the square root of the weights. $\endgroup$ – whuber Jun 9 '16 at 20:35

Your Answer

