#########################
#                       #
# STAT 599 Spring 2013  #
#                       #
# Example R code        #
#                       #
# Chapter 12            #
#                       #
#########################

### Installing the add-on packages needed for this course:

# If you haven't already done so, 
# type at the command line (without the preceding # ):

# install.packages(c('faraway','mgcv','acepack','mda','brlr','sm','wavethresh','SuppDists','lme4'))

# and choose a U.S. mirror site when the window pops up.
# This will install the needed packages for you.




######################################
##
#  Additive Models
##
######################################


# load the 'faraway' package:

library(faraway)

# loading the 'ozone' dataset:

data(ozone)

# attaching the 'ozone' dataset:

# attach(ozone)



# Ordinary linear model:

olm <- lm(O3 ~ temp + ibh + ibt, data=ozone)
summary(olm)

## Additive model, using loess to estimate the f_j functions:

# load the 'gam' package:
library(gam)

ozone.gam.loess <- gam(O3 ~ lo(temp) + lo(ibh) + lo(ibt), data=ozone)
summary(ozone.gam.loess)

# A deviance F-test to check significance of 'ibt' term:
ozone.gam.loess.reduced <- gam(O3 ~ lo(temp) + lo(ibh), data=ozone)
anova(ozone.gam.loess.reduced, ozone.gam.loess, test="F")

# Visualizing the f_j(x_j) functions:
plot(ozone.gam.loess, residuals=TRUE, se=TRUE, pch=".", ask=T)

### We could easily do something similar using smoothing splines:

ozone.gam.ss <- gam(O3 ~ s(temp) + s(ibh) + s(ibt), data=ozone)
summary(ozone.gam.ss)

# A deviance F-test to check significance of 'ibt' term:
ozone.gam.ss.reduced <- gam(O3 ~ s(temp) + s(ibh), data=ozone)
anova(ozone.gam.ss.reduced, ozone.gam.ss, test="F")


### Additive model, using the 'mgcv' package:

# load the 'mgcv' package:
library(mgcv)

ozone.mgcv.ss <- gam(O3 ~ s(temp) + s(ibh) + s(ibt), data=ozone)

plot(ozone.mgcv.ss, residuals=TRUE, se=TRUE, pch=".")

## We may safely drop 'ibt' from the model...


# A deviance-based F-test to see whether temperature actually enters the model in a linear way:

ozone.mgcv.1 <- gam(O3 ~ s(temp) + s(ibh), data=ozone)
ozone.mgcv.lin <- gam(O3 ~ temp + s(ibh), data=ozone)
anova(ozone.mgcv.lin, ozone.mgcv.1, test="F")

#### A test for temp-by-ibh interaction:

# Fitting an additive model that allows for interaction:
ozone.mgcv.int <- gam(O3 ~ s(temp,ibh) + s(ibt), data=ozone)

# Fitting an additive model that does not allow for interaction:
ozone.mgcv.ss <- gam(O3 ~ s(temp) + s(ibh) + s(ibt), data=ozone)

anova(ozone.mgcv.ss, ozone.mgcv.int, test="F")

## A graphical examination of possible interaction:

plot(ozone.mgcv.int)

vis.gam(ozone.mgcv.int, theta=-45, color="gray")

## Using the additive-model results to inform a (piecewise linear) parametric model:

# Defining "hockey-stick" basis functions:

rhs <- function(x,c) ifelse(x>c,x-c,0)
lhs <- function(x,c) ifelse(x<c,c-x,0)

# From the f_j(x_j) function in the 'gam' fit, it appears that the 'temp' trend
# has a changepoint at 60, and the 'ibh' trend has a changepoint at 1000:

ozone.piecewise <- lm(O3 ~ rhs(temp,60) + lhs(temp,60) + rhs(ibh,1000) + lhs(ibh,1000), data=ozone)
summary(ozone.piecewise)

#### Predictions:

# Predict the O3 for an observation with temp=70, ibh=3000, ibt=200:

predict(ozone.mgcv.ss, data.frame(temp=70, ibh=3000, ibt=200), se=T)

## Diagnostics based on residuals:

plot(predict(ozone.mgcv.ss), residuals(ozone.mgcv.ss), xlab='Predicted', ylab='Residuals')

qqnorm(residuals(ozone.mgcv.ss), main='Normal Q-Q plot of Residuals')

############

## Using all the predictors (not just 3):

ozone.all <- gam(O3 ~ s(vh) + s(wind) + s(humidity) + s(temp) + s(ibh) + s(ibt) + s(dpg) + s(vis) + s(doy), data=ozone)
summary(ozone.all)

# Looks like the 'ibh' and 'ibt' terms can be dropped (confirm with F-test):

ozone.most <- gam(O3 ~ s(vh) + s(wind) + s(humidity) + s(temp) + s(dpg) + s(vis) + s(doy), data=ozone)

anova(ozone.most, ozone.all, test="F")

summary(ozone.most)

## The analogous linear model fit does worse (lower adjusted R^2), but could be improved 
## through well-chosen transformations of predictors:

ozone.most.lm <- gam(O3 ~ vh + wind + humidity + temp + dpg + vis + doy, data=ozone)
summary(ozone.most.lm)


#################################################
##
## Section 12.3: Generalized Additive Models
##
#################################################

## Using a Poisson GAM for the 'O3' response:


# load the 'mgcv' package:
library(mgcv)

ozone.pois.gam <- gam(O3 ~ s(temp) + s(ibh) + s(ibt), family=poisson, scale=-1, data=ozone)
summary(ozone.pois.gam)

plot(ozone.pois.gam, residuals=TRUE, se=TRUE, pch=".")



#################################################
##
## Section 12.6: Generalized Additive Mixed Models
##
#################################################


# load the 'faraway' package:

library(faraway)

# loading the 'epilepsy' dataset:

data(epilepsy)

# attaching the 'epilepsy' dataset:

# attach(epilepsy)


# load the 'mgcv' package:
library(mgcv)


# Fitting a GAMM that assumes a Poisson distribution for the response, 
# and assumes correlation within "id" (within subject)

# We include the factors 'treat' and 'expind' (and their interaction)
# in the model, and we model age with a smoothing-spline-based f_j

ep.gamm <- gamm(seizures ~ treat*expind + s(age), family=poisson, random=list(id=~1), data=epilepsy, subset=(id!=49))

# observation 49 is removed because of unusually high seizure counts.

summary(ep.gamm$gam)

# The age predictor appears non-significant.

# Note the essentially flat linear trend in the spline fit:
 plot(ep.gamm$gam,pages=1)


#################################################
##
## Section 12.7: Mulitvariate Adaptive Regression Splines (MARS)
##
#################################################

# load the 'faraway' package:

library(faraway)

# loading the 'ozone' dataset:

data(ozone)

# attaching the 'ozone' dataset:

# attach(ozone)


# load the 'mda' package:
library(mda)

# A MARS model with 7 terms (and no interactions):

a <- mars(ozone[,-1], ozone[,1], nk=7)
summary(lm(ozone[,1] ~ a$x-1))

# A MARS model with 10 terms (and no interactions):

a <- mars(ozone[,-1], ozone[,1], nk=10)
summary(lm(ozone[,1] ~ a$x-1))

# A MARS model with 10 terms (allowing interactions):

a <- mars(ozone[,-1], ozone[,1], nk=10, degree=2)
summary(lm(ozone[,1] ~ a$x-1))

# Looking at the indicator matrix from the last fit:

a$factor[a$selected.terms,]

# A "1" indicates a right-hockey-stick basis function is included in the model for that variable.
# A "-1" indicates a left-hockey-stick basis function is included in the model for that variable.

# If two variables have a nonzero entry on a row, it indicates an interaction between those two variables.

# Plotting the effect of 'ibh' (base height) and
# plotting the effect of 'doy' (day of year):

plot(ozone[,6], a$x[,4]*a$coef[4], xlab="ibh", ylab="Contribution of ibh")
x11()
plot(ozone[,10], a$x[,7]*a$coef[7] + a$x[,6]*a$coef[6], xlab="Day", ylab="Contribution of Day of Year")

## Showing graphically the interaction effect of Temperature and Humidity:

humidity <- seq(10,100,length=20)  # reasonable grid of 'humidity' values
temp <- seq(20,100,length=20)  # reasonable grid of 'temp' values
# setting the other variables to be at their median value:
medians <- apply(ozone,2, median) 


pdf1<-matrix(medians,nrow=400,ncol=10,byrow=T)
pdf1[,4] <- rep(humidity,20)
pdf1[,5] <- rep(temp,rep(20,20))
pdf1 <- as.data.frame(pdf1)
names(pdf1) <- names(medians)

z<- predict(a,pdf1[,-1])
zm <- matrix(z,ncol=20,nrow=20)


### Some diagnostic plots:

qqnorm(a$res, main='Normal Q-Q plot')
x11()
plot(a$fit, a$res, xlab='Fitted', ylab='Residuals')