I am trying to extend the lwr() function of the package McSptial, which fits weigthed regressions as non-parametric estimation. In the core of the lwr() function, it inverts a matrix using solve() instead of a QR decomposition, resulting in numerical instability. I would like to change it but can't figure out how to get the hat matrix (or other derivatives) from the QR decomposition afterward.
With data :
set.seed(0); xmat <- matrix(rnorm(500), nrow=50) ## model matrix
y <- rowSums(rep(2:11,each=50)*xmat) ## arbitrary values to let `lm.wfit` work
w <- runif(50, 1, 2) ## weights
The lwr() function goes as follows :
xmat2 <- w * xmat
xx <- solve(crossprod(xmat, xmat2))
xmat1 <- tcrossprod(xx, xmat2)
vmat <- tcrossprod(xmat1)
I need the value of, for instance :
sum((xmat[1,] %*% xmat1)^2)
sqrt(diag(vmat))
For the moment I use reg <- lm.wfit(x=xmat, y=y, w=w) but cannot manage to get back what seems to me to be the hat matrix (xmat1) out of it.