## Default plots
alcohol <- read.delim("alcohol.txt")
fit <- lm(Metabol~Sex*Gastric,data=alcohol)
plot(fit)

## Influence statistics
influence.measures(fit)

## To make an ANOVA table
fit0 <- lm(Metabol~1,data=alcohol)
fit1 <- lm(Metabol~Sex,data=alcohol)
fit2 <- lm(Metabol~Sex+Gastric,data=alcohol)
fit3 <- lm(Metabol~Sex*Gastric,data=alcohol)
anova(fit0,fit1,fit2,fit3)

##
## Extra code used to create plots in lecture
##
## Read if you wish

x <- runif(100,-10,10)
y <- rnorm(100,mean=x^2,sd=10)
fit <- lm(y~x)
plot(fit$fitted,fit$resid,xlab="Fitted values",ylab="Residuals",pch=19)

x <- runif(100,0,20)
y <- rnorm(100,mean=x,sd=x)
fit <- lm(y~x)
plot(fit$fitted,fit$resid,xlab="Fitted values",ylab="Residuals",pch=19)

fit <- lm(Metabol~Sex*Gastric,data=alcohol)
plot(fit$fitted,fit$resid,xlab="Fitted values",ylab="Residuals",pch=19)

hist(fit$resid,main="",xlab="Residual",col="gray")
boxplot(fit$resid,ylab="Residual")

r <- rdex(300,mean=0,var=1)
hist(r,main="",xlab="Residual",col="gray")
boxplot(r,ylab="Residual")

r <- rdex(300,mean=0,var=1)
qqnorm(r,pch=19,main="Thick tails")
qqline(r,col="gray")
r <- runif(300,-3,3)
qqnorm(r,pch=19,ylim=c(-5,5),main="Thin tails")
qqline(r,col="gray")
r <- rexp(300)
qqnorm(r,pch=19,ylim=c(-2,6),main="Skewed right")
qqline(r,col="gray")
r <- -1*rexp(300)
qqnorm(r,pch=19,ylim=c(-6,2),main="Skewed left")
qqline(r,col="gray")

qqnorm(scale(fit$resid),pch=19)
qqline(scale(fit$resid),col="gray")

calibratedQQ <- function(fit,B=20)
  {
    r <- scale(fit$resid)
    n <- length(r)
    X <- matrix(rnorm(n*B),nrow=n)
    X <- apply(X,2,sort)
    theory <- qnorm(ppoints(n))
    matplot(theory,X,lty=1,type="l",col="gray",xlab="Theoretical quantiles",ylab="Sample quantiles",ylim=c(min(X,r),max(X,r)),main="Calibrated Q-Q plot")
    abline(a=0,b=1,col="black")
    lines(theory,sort(r),col="red",lwd=2)    
  }
calibratedQQ(fit)

plot(hatvalues(fit),type="h",lwd=3,xlab="Subject",ylab="Leverage")
abline(h=4/32,col="gray")

plot(alcohol$Gastric,hatvalues(fit),pch=19,xlab="Gastric",ylab="Leverage")

h <- hatvalues(fit)
plot(fit$fitted,fit$resid,pch=19,cex=h/mean(h),xlab="Fitted values",ylab="Residuals")
plot(fit$fitted,studres(fit),pch=19,cex=h/mean(h),xlab="Fitted values",ylab="Studentized deleted residuals")

ylim <- c(-2.5,2.5)
infl <- as.data.frame(influence.measures(fit)$infmat)
plot(infl$dfb.Gstr,type="h",xlab="Subject",ylab=expression(Delta*hat(beta)),lwd=2,ylim=ylim,main="Slope (Females)")
Female <- 1*(alcohol$Sex=="Female")
fit2 <- lm(Metabol~Female*Gastric,data=alcohol)
infl2 <- as.data.frame(influence.measures(fit2)$infmat)
plot(infl2$dfb.Gstr,type="h",xlab="Subject",ylab=expression(Delta*hat(beta)),lwd=2,ylim=ylim,main="Slope (Males)")
plot(infl[,"dfb.SM:G"],type="h",xlab="Subject",ylab=expression(Delta*hat(beta)),lwd=2,ylim=ylim,main="Interaction")

plot(infl$cook.d,type="h",xlab="Subject",ylab="Cook's Distance",lwd=3)

x1 <- runif(100,-10,10)
x2 <- runif(100,-5,5)
y <- rnorm(100,mean=5*x1+x2^2,sd=10)
fit <- lm(y~x1+x2)
par(mfrow=c(1,2))
plot(fit$fitted,fit$resid,xlab="Fitted values",ylab="Residuals",pch=19)
plot(x2,fit$resid,xlab=expression(x[2]),ylab="Residuals",pch=19)

x1 <- runif(100,-10,10)
x2 <- runif(100,-5,10)
y <- rnorm(100,mean=5*x1+x2^2,sd=10)
fit <- lm(y~x1+x2)
plot(x2,y,pch=19)
plot(x2,fit$resid,xlab=expression(x[2]),ylab="Residuals",pch=19)
pr <- fit$residuals+x2*fit$coef[3]
plot(x2,pr,ylab="Partial residuals",pch=19)

ozone <- read.delim("../../data/ozone.txt")
attach(ozone)
fit <- lm(Ozone~Temp)
plot(Temp,Ozone,pch=19,col="gray",main="Marginal",cex=0.7)
abline(a=fit$coef[1],b=fit$coef[2],col="blue",lwd=2)
fit <- lm(Ozone~Solar+Wind+Day)
pr <- fit$residuals
plot(Temp,pr,ylab="Partial residuals",pch=19,col="gray",main="Adjusted for Wind, etc.",cex=0.7)
fit <- lm(pr~Temp)
abline(a=fit$coef[1],b=fit$coef[2],col="blue",lwd=2)

fit <- lm(Ozone~Solar+Wind+Temp+Day)
X <- model.matrix(fit)
beta <- fit$coefficients
n <- nrow(X)
p <- ncol(X)
xmeans <- apply(X,2,mean)
res <- 101
j <- 4
XX <- matrix(xmeans,nrow=res,ncol=p,byrow=TRUE)
colnames(XX) <- colnames(X)
x <- X[,j]
xx <- seq(min(x),max(x),len=res)
XX[,j] <- xx
est <- XX%*%beta
XXX <- matrix(xmeans,nrow=n,ncol=p,byrow=TRUE)
colnames(XXX) <- colnames(X)
XXX[,j] <- x
rr <- fit$residual+XXX%*%beta
se <- sqrt(apply(XX%*%vcov(fit)*XX,1,sum))
plot(x,rr,type="n",xlab=names(beta)[j],ylab="Regression function")
polygon(c(xx,rev(xx)),c(est-1.96*se,rev(est+1.96*se)),col="gray85",border=F)
lines(xx,est,lwd=2,col="blue")
points(x,rr,pch=19,cex=0.7)