#### R code for section 8.5-2
#### Example: Table 8.5-6 Pavement Durability
data = read.table("table856.dat", head=T)
y = data$Y       # change in rut depth
x1 = data$X1     # viscosity of asphalt
x2 = data$X2     # percentage of asphalt in surface course
x3 = data$X3     # percentage of asphalt in base course
x4 = data$X4     # percentage of fines in surface course
x5 = data$X5     # percentage of voids in surface course
x6 = data$X6     # indicate variable: two sets of runs

## pre-analysis
pairs(data)
round(cor(data),3)

# Box-Cox transformation
library(MASS)
boxcox(y~x1)
boxcox(y~x1+x2+x3+x4+x5+x6)

logy = log(y)
pairs(cbind(logy,x1,x2,x3,x4,x5,x6))
boxcox(x1~logy)
logx1 = log(x1)

## plot y and x1
par(mfrow=c(1,2))
plot(x1, y, main="Scatter Plot")
plot(log(x1), log(y), main="Scatter Plot in Log Scale")

## model 1
fit1 = lm(logy~logx1)
par(mfrow=c(2,2))
summary(fit1)
plot(fit1,1)
plot(fit1,2)
plot(fit1,3)
plot(fit1,4)

## model 2
fit2 = lm(logy~logx1+x6)
summary(fit2)
par(mfrow=c(2,2))
plot(fit2)

## model 3
fit3 = lm(logy~logx1+x2+x3+x4+x5+x6)
summary(fit3)
par(mfrow=c(2,2))
plot(fit3)


## model 4
fit4 = lm(logy~logx1+x2+x5+x6)
summary(fit4)
par(mfrow=c(2,2))
plot(fit4)


## model 5 (interaction)
x16 = logx1*x6
fit5 = lm(logy~logx1+x2+x5+x6+x16)
summary(fit5)
par(mfrow=c(2,2))
plot(fit5)

## Choose a model by AIC (optional)
## AIC: Akaike information criterion
## AIC = 2(p+1) - 2*log(L) = 2(p+1) + n*log(RSS/n), compared with BIC=(p+1)*log(n) - 2*log(L)
step(fit3)