library(mlbench)

##load Bostong Housing data
data(BostonHousing)
plot(BostonHousing,pch=20,cex=0.5,col=rgb(0.2,0.2,0.2,0.5))




## prepare Y and X matrix
Y <- BostonHousing$medv
X <- BostonHousing[,1:13]
## make 4th variable a numeric
X[,4] <- as.numeric(X[,4])
n <- nrow(X)

## fit a linear model
linmod <- lm(Y~ X)
summary(linmod)

## add intercept term (first column) by hand for the following
X <- cbind(rep(1,n),as.matrix(X))

## Projection matrix P=Proj
Proj <- X %*% solve( t(X) %*% X) %*% t(X)

## compare fitted values from linear model and "direct" projection PY
plot(fitted(linmod), Proj %*% Y)
## compare residuals with QY
plot(residuals(linmod), (diag(n) - Proj) %*% Y)


## compute eigenvalue decomposition of P=Proj
eigProj <- eigen(Proj)
## look at eigenvalues
plot(eigProj$values)