#### Figure 8.3-1
## (a)
set.seed(123)
n = 23
x = runif(n, 0, 1.5)
y = 0.5 + 2*x - x^2 + rnorm(n)*0.05
fit1 = lm(y~x)
par(mfrow=c(2,3))
plot(x, y, main="(a)")
abline(fit1$coef[1],fit1$coef[2])

## (c)
plot(fit1$fitted, fit1$residual, xlab=expression(hat(y)), 
     ylab="e", main="(c)")
abline(0,0)

## (e)
plot(x, fit1$residual, xlab="x", ylab="e", main="(e)")
abline(0,0)

## (b)
set.seed(456)
n = 41
x = abs(rnorm(n))
y = 2 - 0.5*x + rnorm(n)*0.5*(max(x)-x)^1.8
fit2 = lm(y~x)
plot(x, y, main="(b)")
abline(fit2$coef[1],fit2$coef[2])

## (d)
plot(fit2$fitted, fit2$residual, xlab=expression(hat(y)), 
     ylab="e", main="(c)")
abline(0,0)

## (f)
plot(x, fit2$residual, xlab="x", ylab="e", main="(e)")
abline(0,0)


#### Example 8.3-1 on page 459
temp = read.table("exercise836.dat", head=T)
x = temp$Weight
y = temp$MPG
fit3 = lm(y~x)

## Figure 8.3-2 on page 460
par(mfrow=c(1,3))
plot(x, y, main="(a)", xlab="Weight", ylab="MPG")
abline(fit3$coef[1],fit3$coef[2])

plot(fit3$fitted, fit3$residual, xlab=expression(hat(y)), 
     ylab="Residual", main="(b)")
abline(0,0)

plot(x, fit3$residual, xlab="Weight", ylab="Residual", main="(c)")
abline(0,0)

## add x^2 into consideration *** multiple linear regression
x2 = x^2
fit4 = lm(y~x+x2)
par(mfrow=c(1,3))
plot(x, y, main="(a)", xlab="Weight", ylab="MPG")
points(x[order(x)], fit4$fitted[order(x)], type="l")
plot(fit4$fitted, fit4$residual, xlab=expression(hat(y)), 
     ylab="Residual", main="(b)")
abline(0,0)
plot(x, fit4$residual, xlab="Weight", ylab="Residual", main="(c)")
abline(0,0)

## do transformation on y
gpm = 100/y
fit5 = lm(gpm~x)
## Figure 8.3-3 on page 461
par(mfrow=c(1,3))
plot(x, gpm, main="(a)", xlab="Weight", ylab="GPM")
abline(fit5$coef[1],fit5$coef[2])
plot(fit5$fitted, fit5$residual, xlab="Fitted values", 
     ylab="Residuals", main="(b)")
abline(0,0)
plot(x, fit5$residual, xlab="Weight", ylab="Residual", main="(c)")
abline(0,0)


### check output of linear regression
summary(fit5)
summary(fit5$residual)
n = length(x)
y2bar = mean(gpm^2)
ybar = mean(gpm)
x2bar = mean(x^2)
xbar = mean(x)
xybar = mean(x*gpm)
SSTO = n*(y2bar - ybar^2)
SSR = n*(xybar - xbar*ybar)^2/(x2bar - xbar^2)
SSE = SSTO - SSR
R2 = SSR/SSTO
# adjusted R^2
1 - (1-R2)*(n-1)/(n-1-1)


### use Box-Cox transformation, updated on 02/15/2010
library("MASS")
boxcox(y~x)