[R] residuals from lm

[Rolf Turner]

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