## Read in data
ozone <- read.delim("../../data/ozone.txt")
attach(ozone)

## Scatterplot matrix
pairs(ozone)
panel.histogram <- function(x,...)
  {
    par(new = TRUE)
    fit <- hist(x, plot=FALSE)
    hist(x, main = "", axes = FALSE, col = "gray",ylim=c(0,max(fit$counts)*1.5))
  }
pairs(ozone,pch=19,cex=.8,diag.panel=panel.histogram)

## Slide 10
fit <- lm(Ozone~Solar+Wind+Temp+Day)

## Slide 11
y <- Ozone
X <- cbind(1,as.matrix(ozone[,-1]))
dim(X)
beta <- solve(t(X) %*% X) %*% t(X) %*% y

## Slide 12
univariate.beta <- c(coef(lm(Ozone~Solar))[2],
                     coef(lm(Ozone~Wind))[2],
                     coef(lm(Ozone~Temp))[2],
                     coef(lm(Ozone~Day))[2])
cbind(coef(fit)[-1],univariate.beta)

## Slide 15-16
scaled.beta <- fit$coef[-1] * apply(ozone[,-1],2,sd) / sd(Ozone)
scaled.univ.beta <- c(coef(lm(scale(Ozone)~0+scale(Solar))),
                      coef(lm(scale(Ozone)~0+scale(Wind))),
                      coef(lm(scale(Ozone)~0+scale(Temp))),
                      coef(lm(scale(Ozone)~0+scale(Day)))) ## Same as correlation
lm(Ozone~0+Solar+Wind+Temp+Day,
   data=as.data.frame(scale(ozone)))
cbind(coef(fit)[-1],univariate.beta,scaled.beta,scaled.univ.beta)

## Slide 18
x <- c(1,180,15,70,274)
crossprod(x,fit$coef)
x <- c(1,180,5,90,274)
crossprod(x,fit$coef)

## Slide 21
Extra <- 3*Solar - 2*Wind + 0.5*Temp
lm(Ozone~Solar+Wind+Temp+Day+Extra)
XX <- cbind(1,Solar,Wind,Temp,Day,Extra)
solve(crossprod(XX))