# RPD -- Exercise 12.7 -- Watershed Flow Data (Box-Cox)

watershed <- read.table("watershed.dat", header=TRUE)

  # use log-transformed independent variables, as suggested in exercise:

watershed.model <- lm(peakrate ~ log(area) + log(areaimperv) + log(slope) +
                                 log(flow) + log(absorbency) + log(capacity) +
                                 log(infiltration) + log(rainfall) + 
                                 log(raintime), data=watershed)

plot(fitted(watershed.model), residuals(watershed.model))


Y <- watershed$peakrate

X <- with(watershed, cbind(1, log(area), log(areaimperv), log(slope),
                           log(flow), log(absorbency), log(capacity),
                           log(infiltration), log(rainfall), log(raintime)))

XtXinv <- solve(t(X)%*%X)

Ydot <- exp(sum(log(Y))/length(Y))

lambdas <- seq(-1.0,1.0,by=0.05)
SSRes <- rep(0,length(lambdas))
for(i in 1:length(lambdas)){

  lambda <- lambdas[i]

  if(round(lambda,digits=2) == 0)
    Ytrans <- Ydot*log(Y)
  else 
    Ytrans <- (Y^lambda - 1)/(lambda*Ydot^(lambda-1))

  Resid <- Ytrans - X%*%(XtXinv%*%(t(X)%*%Ytrans))

  SSRes[i] <- sum(Resid^2)

}

print(data.frame(lambdas, SSRes))

nu <- length(Y) - ncol(X)

CI.95lim <- min(SSRes) * (1 + qt(1-0.05/2,nu)^2/nu)

print(CI.95lim)

dev.new()
plot(lambdas, SSRes, ylim=c(0,1000000), type="b")
abline(h=CI.95lim)

# Log transformation (lambda = 0) is almost within the confidence limit,
# but not quite.


trans.watershed.model <- lm((peakrate)^(-0.15) ~ log(area) + log(areaimperv) +
                                                 log(slope) + log(flow) +
                                                 log(absorbency) +
                                                 log(capacity) +
                                                 log(infiltration) +
                                                 log(rainfall) + log(raintime),
                            data=watershed)

dev.new()
plot(fitted(trans.watershed.model), residuals(trans.watershed.model))


# Note that the transformation is actually a decreasing function.  Therefore,
# a negative coefficient for a variable implies that an increase in that
# variable is associated with increase in peakrate (all else constant).


# There is also a function "boxcox" for fitting Box-Cox transformations
# directly.  It requires the MASS package.