Davis=read.table(
"http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]

summary(lm(height~X1F+X1H,data=Davis))


X=cbind(rep(1,200),Davis$weight)
Y=Davis$height
solve(t(X) %*% X) %*% (t(X) %*% Y)

X=cbind(rep(1,200),Davis$X1F,Davis$X1M)
solve(t(X) %*% X)

X1=X[,c(2,3)]
solve(t(X1) %*% X1)
solve(t(X1) %*% X1) %*% (t(X1) %*% Y)

X1=X[,c(1,3)]
solve(t(X1) %*% X1)
solve(t(X1) %*% X1) %*% (t(X1) %*% Y)

X=cbind(rep(1,200),Davis$weight)
H=X %*% solve(t(X) %*% X) %*% t(X)


beta=solve(t(X) %*% X) %*% (t(X) %*% Y)
Yhat=X%*%beta
plot(X[,2],Y,xlim=c(0,120),ylim=c(130,200))
points(X[,2],Yhat,col="green")
abline(beta[1],beta[2],col="red")
epsilon=Y-Yhat
s2=var(epsilon)
V=as.numeric(s2)*solve(t(X) %*% X)
sqrt(diag(V))
beta[2]
u=seq(0,1,by=.001)
plot(u,dnorm(u,beta[2],sd=sqrt(diag(V))[2]),type="l")


X=X[150:200,]
Y=Y[150:200]
beta=solve(t(X) %*% X) %*% (t(X) %*% Y)
Yhat=X%*%beta
epsilon=Y-Yhat
s2=var(epsilon)
V=as.numeric(s2)*solve(t(X) %*% X)
sqrt(diag(V))
beta[2]
u=seq(0,1,by=.001)
lines(u,dnorm(u,beta[2],sd=sqrt(diag(V))[2]),type="l",col="blue")