### R code from vignette source 'nnet.Rnw'

###################################################
### code chunk number 1: nnet.Rnw:12-13
###################################################
options(show.signif.stars=FALSE,digits=4,width=80)


###################################################
### code chunk number 2: f1
###################################################
data(ozone, package="faraway")
head(ozone)
# transform ibt to get correct units, deg. C
ozone$ibt <- 32 + 0.18 * ozone$ibt
# use same variables as faraway:
# O3, ozone concentration measured at Sandburg AFB E of Los Angeles
# temp, Daily (maximum?) temerature
# ibh, Inversion base height, feet, a measure of atmospheric conditions
# ibt, inversion base temperature, degrees F
library(car)
spm(~O3 + temp + ibh + ibt, ozone)


###################################################
### code chunk number 3: nnet.Rnw:82-85
###################################################
(summary(p1 <- powerTransform(cbind(ibh, ibt, temp) ~ 1, ozone)))
# ibh is apparently truncated at 5000 feet; lets take logs to lessen the
#    effect of truncation


###################################################
### code chunk number 4: f2
###################################################
spm(~ O3 + temp + log(ibh) + ibt, ozone)


###################################################
### code chunk number 5: nnet.Rnw:95-104
###################################################
# Let's do cross-validation by default:
set.seed(50)
construct <- rep(FALSE, 330)
construct[sample(1:330, 250)] <- TRUE
check.fit <- function(response=ozone$O3, fit, group=construct) {
 error <- (response - fit)^2
 data.frame(Construction = sqrt(sum(error[group])/sum(group)),
            Validation= sqrt(sum(error[!group])/sum(!group)))
 }


###################################################
### code chunk number 6: nnet.Rnw:109-112
###################################################
library(nnet)
nmd1 <- nnet(O3 ~ temp + log(ibh) + ibt, ozone, subset=construct,
    size=2, linout=TRUE )


###################################################
### code chunk number 7: nnet.Rnw:115-119
###################################################
# maximum RSS by group
check.fit(ozone$O3, rep(mean(ozone$O3[construct]), 330))
# RSS of fitted model
check.fit(ozone$O3, predict(nmd1, new=ozone))


###################################################
### code chunk number 8: nnet.Rnw:128-131
###################################################
ozone$logibh <- log(ozone$ibh)
sx <- scale(ozone) # centered and scaled
sx1 <- as.data.frame(sx)


###################################################
### code chunk number 9: f3
###################################################
maxrss <- (199) * var(sx1$O3)
bestrss <- dim(sx1)[1] -1
ntries <- 100
rss <- rep(0, ntries)
for (i in 1:ntries){
  nnmd1 <- nnet(O3 ~ temp + logibh + ibt, sx1, size=2, 
  subset=construct, linout=TRUE, trace=FALSE)
  rss[i] <- nnmd1$value
  if (nnmd1$value < bestrss){
     bestrss <- nnmd1$value
     bestnn <- nnmd1}
  }
hist(rss, main=paste("RSS min =", round(bestrss,3), ", R^2 =",
   round(1-bestrss/maxrss, 2)))


###################################################
### code chunk number 10: nnet.Rnw:152-155
###################################################
summary(bestnn)
check.fit(sx1$O3, rep(mean(sx1$O3[construct]), 330) )
check.fit(sx1$O3, predict(bestnn, new=sx1) )


###################################################
### code chunk number 11: nnet.Rnw:161-163
###################################################
m1 <- lm(O3 ~ temp + logibh + ibt, sx1, subset=construct)
check.fit(sx1$O3, predict(m1, newdata=sx1))


###################################################
### code chunk number 12: f4
###################################################
library(mgcv)
summary(g1 <- gam( O3 ~ s(temp) + s(logibh) + s(ibt), data= sx1, subset=construct))
check.fit(sx1$O3, predict(g1, newdata=sx1))
plot(g1, pages=1)


###################################################
### code chunk number 13: f5
###################################################
library(randomForest)
r1 <- randomForest(O3 ~ temp + logibh + ibt, sx1, subset=construct)
check.fit(sx1$O3, predict(r1, newdata=sx1))
op <- par(mfrow=c(1, 3))
partialPlot(r1, sx1, "temp")
partialPlot(r1, sx1, "logibh")
partialPlot(r1, sx1, "ibt")
par(op)


###################################################
### code chunk number 14: f6
###################################################
par(mfrow=c(1,3))
xx <- expand.grid(temp=seq(-3, 3, .1), logibh=0, ibt=0)
plot(xx$temp, predict(bestnn, new=xx), xlab="temp", ylab="predicted O3", lwd=3,type="l")
xx <- expand.grid(temp=seq(-3, 3, .1), logibh=-1, ibt=0)
lines(xx$temp, predict(bestnn, new=xx), lwd=1, type="l")
xx <- expand.grid(temp=seq(-3, 3, .1), logibh=1, ibt=0)
lines(xx$temp, predict(bestnn, new=xx), lwd=2, type="l")
xx <- expand.grid(temp=0, logibh=seq(-3, 3, .1), ibt=0)
plot(xx$logibh, predict(bestnn, new=xx), xlab="logibh", ylab="predicted O3")
xx <- expand.grid(temp=0, logibh=0, ibt=seq(-3, 3, .1))
plot(xx$ibt, predict(bestnn, new=xx), xlab="ibt", ylab="predicted O3")


###################################################
### code chunk number 15: f7
###################################################
bestrss1 <- dim(sx1)[1] -1
rss <- rep(0, ntries)
for (i in 1:ntries){
  nnmd1 <- nnet(O3 ~ temp + logibh + ibt, sx1[construct, ] , size=2, linout=TRUE,
           trace=FALSE, decay=.03)
  rss[i] <- nnmd1$value
  if (nnmd1$value < bestrss1){
     bestrss1 <- nnmd1$value
     bestnn1 <- nnmd1}
  }
check.fit(sx1$O3, predict(bestnn1, new=sx1))
par(mfrow=c(1, 3))
xx <- expand.grid(temp=seq(-3, 3, .1), logibh=0, ibt=0)
plot(xx$temp, predict(bestnn1, new=xx), xlab="temp", ylab="predicted O3", lwd=3,type="l")
xx <- expand.grid(temp=seq(-3, 3, .1), logibh=-1, ibt=0)
lines(xx$temp, predict(bestnn1, new=xx), lwd=1, type="l")
xx <- expand.grid(temp=seq(-3, 3, .1), logibh=1, ibt=0)
lines(xx$temp, predict(bestnn1, new=xx), lwd=2, type="l")
xx <- expand.grid(temp=0, logibh=seq(-3, 3, .1), ibt=0)
plot(xx$logibh, predict(bestnn1, new=xx), xlab="logibh", ylab="predicted O3")
xx <- expand.grid(temp=0, logibh=0, ibt=seq(-3, 3, .1))
plot(xx$ibt, predict(bestnn1, new=xx), xlab="ibt", ylab="predicted O3")


###################################################
### code chunk number 16: nnet.Rnw:232-236
###################################################
mean.predictions <-
  (as.vector(predict(bestnn, new=sx1)) + predict(m1, new=sx1) + predict(g1, newdata=sx1) +
  predict(r1, new=sx1) + as.vector(predict(bestnn1, new=sx1)))/5
check.fit(sx1$O3, mean.predictions)