## On lecture 9 February 2015
# How does pH in Nowegian lakes depend on sulfate, nitrate, calcium, aluminium and organic content (x1, ... x5), area of lake (x6) and location (x7 = 0, Telemark, or x7 = 1, Trøndelag)? Data from Statens forurensningstilsyn (1986). Here 26 random lakes from Telemark and Trøndelag out of 1005 lakes have been drawn 

acidrain <- read.table("http://www.math.ntnu.no/~mettela/TMA4267/Data/acidrain.txt",header=TRUE)

fit <- lm(y~.,data=acidrain) # lm: linear model
# or
attach(acidrain)
fit<-lm(y~x1+x2+x3+x4+x5+x6+x7)
# 1 is added by R as a covariate for both alternatives

summary(fit)

n <- length(y)

x <- cbind(rep(1,n),acidrain[,2:8])
names(x)[1] <- 1
x <- as.matrix(x)
p <- dim(x)[2]

# LS estimates of beta
bhat <- solve(t(x)%*%x)%*%t(x)%*%y
# t: transpose; %*%: matrix multiplication; solve: invert
fit$coefficients

# SSE
t(y-x%*%bhat)%*%(y-x%*%bhat) # or
sum((y-x%*%bhat)^2) # or
sse <- sum(fit$residuals^2)

# ML estimate of sigma^2
sse/n # or
summary(fit)$sigma^2*(n-p)/n # summary(fit)$sigma^2 is an unbiased estimate

# QR decomposition of X
qrdec <- qr(x)
q <- qr.Q(qrdec)
r <- qr.R(qrdec)
# view q, r, q%*%*r, x
solve(r)%*%t(q)%*%y
qr.solve(x,y) # "solves" the over-determined system x%*%betahat = y by QR decomposition of x and least squares - we get the same solution as bhat

## 17 February 2015:

h <- x%*%solve(t(x)%*%x)%*%t(x)
h%*%h-h # zero matrix since h is idempotent

h%*%y # fitted values
fit$fitted.values # computed by R

# From (c) p. 80 in B/F - NOT to be confused by SS and SSR:

i <- diag(n) # identity matrix
j <- matrix(1,nrow=n,ncol=n) # matrix of 1s, from lecture

sstot <- sum((y-mean(y))^2)
t(y)%*%(i-j/n)%*%y # the same in matrix language

ssreg <- sum((fit$fitted.values-mean(y))^2)
t(y)%*%(h-j/n)%*%y # matrix language

t(y)%*%(i-h)%*%y # sse in matrix language

(i-h)%*%(h-j/n) # 0 matrix, as it should be

sstot
sse+ssreg # the same since the model ('fit' above)  has an intercept

sum(diag(h)) # p = 8
sum(diag(i-h)) # n-p = 18

1-sse/sstot
summary(fit)$r.squared # the same again computed by lm.summary

sum(y^2)
sum(fit$fitted.values^2)+sum(fit$residuals^2) # they are equal - Pythagoras

# t statistic for testing whether beta1 = 0 (coeff. of x1, not intercept):
tstat<-fit$coefficients[2]/(summary(fit)$sigma*sqrt(solve(t(x)%*%x)[2,2]))
# Compare with t values in summary(fit):
summary(fit)
# p-value:
2*pt(-abs(tstat),n-p) # again compare with p-value from summary(fit)