Weighted Least Squares by lm-function vs. manual calculation

Question

I am trying to manually calculate beta-coefficients using Weighted Least Squares, which are given by:

X should comprise only one variable and the coefficients should include an intercept. I tried it as follows:

set.seed(1)
x = as.matrix(mtcars$wt)
y = as.matrix(mtcars$mpg)
w = runif(length(x))

(t(cbind(1,x)) %*%  diag(w) %*% cbind(1,x))^(-1) %*% 
     t(cbind(1,x)) %*%  diag(w) %*%  y

This leads to:

         [,1]
[1,] 38.92461
[2,] 11.52764

However, the lm-Funktion leads to different results:

lm(y~x, weights = w)

Call:
lm(formula = y ~ x, weights = w)

Coefficients:
(Intercept)            x  
     37.896       -5.437  

I am quite sure it must have something to do with the cbind(1,x)-part, as (t(cbind(x)) %*% diag(w) %*% cbind(x))^(-1) %*% t(cbind(x)) %*% diag(w) %*% y leads to the same results as lm(y~x -1, weights = w). Does anyone see what I did wrong while trying to calculate the coefficients manually?

Answer 1

Both of them return the same results.

Set $A=X^T WX$. When computing the fitted values $\hat{\beta}_{WLS}$, the first half of the expression is the inverse of $A$. However, as in your code above, literally putting an exponent of $(-1)$ will simply perform $\frac{1}{A_{ij}}$ to every element $A_{ij}$ of the matrix.

Instead, to find the inverse of the matrix in R, simply type solve(A). With these changes the results of the Weighted Regression manually and in R will match.

Hope this helps.