9

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

amoeba
  • 93,463
  • 28
  • 275
  • 317
Haitao Du
  • 32,885
  • 17
  • 118
  • 213
  • 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. – amoeba Oct 06 '16 at 00:02

1 Answers1

13

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)
Coefficients:
      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.)

whuber
  • 281,159
  • 54
  • 637
  • 1,101
  • 1
    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. – whuber Jun 09 '16 at 20:25
  • 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$ – Haitao Du Jun 09 '16 at 20:28
  • 8
    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. – whuber Jun 09 '16 at 20:35