On 03/07/12 03:58, chuck.01 wrote:
Hi,
I was playing around with something else and I noticed this matrix code for
residuals in a linear model doesn't say what lm() says. Please tell me if I
am completely misguided here.
data(mtcars)
Y <- as.matrix(mtcars[,1])
X <- as.matrix(mtcars[,c(2:11)])
# shouldnt this:
H <- X %*% solve(t(X) %*% X) %*% t(X)
(diag(dim(H)[1]) - H) %*% Y
# be equal to this:
residuals(lm(Y~X))
# ???
# thanks
You are forgetting about the constant term.
Try
X <- cbind(1,X)
H <- X %*% solve(t(X) %*% X) %*% t(X)
R1 <-(diag(dim(H)[1]) - H) %*% Y
R2 <-residuals(lm(Y~X))
range(R1-R2)
I get:
[1] -3.886808e-12 1.262768e-11
OMMMMMMMMMMMMMMMM!!! :-)
cheers,
Rolf Turner
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
