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

## Slide 3
fit <- lm(Ozone~Solar+Wind+Temp+Day)
fit$coefficients
fit$fitted.values
fit$residuals
fit$rank
fit$df.residual
n <- nrow(ozone)
p <- fit$rank

## Slide 5
sig2 <- sum(fit$residuals^2)/fit$df.residual
sig <- sqrt(sig2)

## Slide 6
X <- as.matrix(cbind(1,ozone[,-1]))
VarB <- sig2*solve(crossprod(X))
summ <- summary(fit)
summ$sigma
summ$sigma^2*summ$cov.unscaled

## Slide 7
sqrt(diag(VarB))

## Slide 9
lambda <- c(0,0,-5,10,0)
crossprod(lambda,fit$coefficients)
sqrt(t(lambda) %*% VarB %*% lambda)
c(-5,10) %*% VarB[3:4,3:4] %*% c(-5,10)

## Slide 10
25*VarB[3,3] + 100*VarB[4,4]
lambda <- c(0,0,5,10,0)
crossprod(lambda,fit$coefficients)
sqrt(t(lambda) %*% VarB %*% lambda)

## Slide 11
a <- c(1,180,15,70,274)
crossprod(a,fit$coefficients)
sqrt(t(a) %*% VarB %*% a + sig2)
b <- c(1,180,5,90,274)
crossprod(b,fit$coefficients)
sqrt(t(b) %*% VarB %*% b)

## Slide 13
sqrt(t(a) %*% VarB %*% a + sig2)
crossprod(a,fit$coefficients) + c(-1,1)*qt(.975,106)*sqrt(t(a) %*% VarB %*% a + sig2)

## Slide 14
var(Ozone)
var(fit$residuals) + var(fit$fitted.values)
TSS <- crossprod(Ozone-mean(Ozone))
RSS <- crossprod(fit$residuals)
MSS <- crossprod(fit$fitted.values-mean(fit$fitted.values))
MSS/TSS
cor(fit$fitted.values,Ozone)^2

## Slide 17
Wind2 <- Wind + rnorm(n,mean=0,sd=20)
cor(Wind,Wind2)
summ <- summary(lm(Ozone~Solar+Wind+Temp+Day))
summ2 <- summary(lm(Ozone~Solar+Wind+Temp+Day+Wind2))
diag(summ$sigma^2*summ$cov.unscaled)
diag(summ2$sigma^2*summ2$cov.unscaled)

## Slide 18 (simpler code, uglier plot)
SD <- seq(20,.1,len=200)
Correlation <- numeric(200)
VarHat <- matrix(NA,nrow=200,ncol=6)
for (i in 1:200)
  {
    Wind2 <- Wind + rnorm(n,mean=0,sd=SD[i])
    Correlation[i] <- cor(Wind,Wind2)
    summ <- summary(lm(Ozone~Solar+Wind+Temp+Day+Wind2))
    VarHat[i,] <- diag(summ$sigma^2*summ$cov.unscaled)
  }
matplot(Correlation,VarHat[,-1],pch=19,cex=.5,ylim=c(0,10))

## Slide 18 (nicer plot using the lattice package)
SD <- seq(20,.1,len=200)
Results <- NULL
for (i in 1:200)
  {
    Wind2 <- Wind + rnorm(n,mean=0,sd=SD[i])
    summ <- summary(lm(Ozone~Solar+Wind+Temp+Day+Wind2))
    Results <- rbind(Results,data.frame(Correlation=cor(Wind,Wind2),VarHat=diag(summ$sigma^2*summ$cov.unscaled)[-1],Term=attr(summ$terms,"term.labels")))
  }
trellis.par.set(superpose.symbol=list(pch=19,col=c("black","red","green","blue","steelblue")))
xyplot(VarHat~Correlation,Results,group=Term,ylim=c(-0.5,10.5),auto.key=list(columns=5))

## Slide 19
WindSq <- Wind^2
summary(lm(Ozone~Solar+Wind+WindSq+Temp+Day))

## Slide 21
fit <- lm(log(Ozone)~Solar+Wind+Temp+Day)
summary(fit)