# 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