# =============================================
# Auto data
# =============================================
#
# Summary: 
# Data used in the book:
#    "Data anlysis and data mining" by A.Azzalini and B.Scarpa,
#    © Oxford University Press, 2012 (ISBN 978-0-19-976710-6).
#
# Data obtained by manipulation of the data available at
# ftp://ftp.ics.uci.edu/pub/machine-learningdatabases/autos
#
# Data description: 
# The data set consists of the following 18 variables:
# Response variable:
#  city.distance: city distance covered (km/l)
# Covariates:
#  make:                manufacturer (factor, 22 levels)    
#  fuel.type:           type of engine fuel (factor, 2 levels: diesel, gasoline)   
#  aspiration:          type of engine aspiration (factor, 2 levels: standard, turbo)   
#  body.style:          type of body style (factor, 5 levels: hardtop, wagon, sedan, hatchback, convertible)  
#  drive.wheels:        type of drive wheels (factor,  3 levels: 4wd, fwd, rwd)   
#  engine.location:     location of engine (factor, 2 levels: front, rear)  
#  wheel.base:          distance between axes (cm)   
#  length:              length (cm)    
#  width:               width (cm)  
#  height:              height (cm)   
#  curb.weight:         weight (kg)   
#  engine.size:         engine size (l)   
#  compression.rate:    compression rate   
#  hp:                  horse power   
#  peak-rpm:            number of peak revolutions per minute   
#  highway.distance:    highway distance (km/l)  
#  n.cylinders:         number of cylinders
#
#======================================================

## Read the data
auto <- read.table("http://azzalini.stat.unipd.it/Book-DM/auto.dat", header = TRUE)

## Fit a an additive model
library(gam)
fit.gam <- gam(city.distance~s(engine.size)+s(curb.weight)+s(HP)+factor(aspiration),data=auto)
par(mfrow=c(1,3))
plot(fit.gam,se=TRUE,terms=c("s(engine.size)","s(curb.weight)","s(HP)"))
summary(fit.gam)

## Fit a linear model and compare
fit.lm <- lm(city.distance ~ engine.size+curb.weight+HP+factor(aspiration),data=auto)
summary(fit.lm) 
anova(fit.lm,fit.gam,test="F")

## Fit a regression tree
library(tree)
set.seed(34534)
N = nrow(auto)
train = sample(1:N,101,replace=FALSE)

N.tr <- length(train)
fit.t1 <- tree(city.distance~engine.size+curb.weight,data=auto[train,],split="deviance",control=tree.control(mincut=2,minsize=4,nobs=N.tr))
plot(fit.t1)
text(fit.t1, digits=2, pretty=3,cex=0.5)
pred.t1 = predict(fit.t1,auto[-train,])
err.t1 <- mean((auto$city.distance[-train]-pred.t1)^2)

# Prune the tree with 10-fold cross validation
K <- 10
t1.size <- sum(fit.t1$frame$var == "<leaf>")
index <- sample(rep(1:K, ceiling(N.tr/ K))[1:N.tr])# split the data in approximately equisize folds
err.CV <- rep(0,t1.size)
fit.p1 <- prune.tree(fit.t1,newdata=auto[train,],method="deviance")
for(k in 1:K)
{
  fit.t1.k <- tree(city.distance~engine.size+curb.weight,data=auto[train[index != k],],split="deviance",control=tree.control(mincut=2,minsize=4,nobs=sum(index != k)))
  fit.p1.k <- prune.tree(fit.t1.k,,newdata=auto[train[index == k],],method="deviance")
  err.CV[fit.p1.k$size] <- err.CV[fit.p1.k$size]+fit.p1.k$dev
}
plot(1:t1.size,err.CV[1:t1.size],type="s")
# Alternatively, we use the built-in function in R
fit.p1.CV <- cv.tree(fit.t1,,FUN=prune.tree,K=10)
plot(fit.p1.CV)
best.size <- fit.p1.CV$size[which.min(fit.p1.CV$dev)]
best.size
# Fit the best tree
fit.p1.best <- prune.tree(fit.t1, newdata=auto[train,],best=best.size) 
plot(fit.p1.best)
text(fit.p1.best, digits=2, pretty=3,cex=0.5)
pred.p1 = predict(fit.p1.best,auto[-train,])
err.p1 <- mean((auto$city.distance[-train]-pred.p1)^2)

c(err.t1,err.p1)

# =============================================
# E-mail spam data
# =============================================
#
# Summary: 
# Information from 4601 e-mail messages with the purpose of 
# constructing an automatic spam detector
#
# Data description: 
# The data set consists of the following 59 variables:
# Response variable:
#  y: indicator of whether the e-mail is spam True/False
# Covariates:
#  x1-x57: 57 of the most  commonly occuring words and 
#          punctuation marks in the e-mails
# Indicator of training and test sets:
#  train: indicator of whether to use as training or test sample
#
#===============================================

## Read data
spam = read.table("http://www.uio.no/studier/emner/matnat/math/STK2100/data/spam_data.txt",header=TRUE)

# Transform input data due to long-tailed distribution
spam[,1:57] = log(spam[,1:57]+0.1)

## Fit logistic GAM with all covariates
library(gam)
nam = names(spam)[1:57]
k = 57
formula = as.formula(paste("y",paste(paste("s(",nam[1:k],")",sep=""),collapse="+"), sep="~"))
fit.gam = gam(formula,data=spam,subset=spam$train,family=binomial)
pred.gam = predict(fit.gam,spam[!spam$train,],type="response")>0.5
err.gam=mean(spam$y[!spam$train]!=pred.gam)

# Compare with the ordinary logistic regression model
fit = glm(y~.-train,data=spam,family=binomial,subset=spam$train)
pred = predict(fit,spam[!spam$train,],type="response")>0.5
err.logist=mean(spam$y[!spam$train]!=pred)

c(err.gam,err.logist)
summary(fit.gam)

# Using variables selected in textbook
ind = c(5:8,16:17,25:27,37,45:46,52:53,56:57)

formula = as.formula(paste("y",paste(paste("s(",nam[ind],")",sep=""),collapse="+"), sep="~"))
fit.gam.sub = gam(formula,data=spam,subset=spam$train,family=binomial)
pred.gam.sub = predict(fit.gam.sub,spam[!spam$train,],type="response")>0.5
err.gam.sub=mean(spam$y[!spam$train]!=pred.gam.sub)

formula = as.formula(paste("y",paste(paste(nam[ind],sep=""),collapse="+"), sep="~"))
fit.glm.sub = glm(formula,data=spam,subset=spam$train,family=binomial)
pred.glm.sub = predict(fit.glm.sub,spam[!spam$train,],type="response")>0.5
err.glm.sub=mean(spam$y[!spam$train]!=pred.glm.sub)

c(err.gam.sub,err.glm.sub)
summary(fit.gam.sub)

# Replace some of the terms with linear effects
ind.smooth <- c(6:7,17,27,37,45:46,52:53,56)
ind.lin <- c(5,8,16,25:26,57)
formula = as.formula(paste("y",paste(c(paste("s(",nam[ind.smooth],")",sep=""),nam[ind.lin]),collapse="+"), sep="~"))
fit.gam.sub.lin = gam(formula,data=spam,subset=spam$train,family=binomial)
pred.gam.sub.lin = predict(fit.gam.sub.lin,spam[!spam$train,],type="response")>0.5
err.gam.sub.lin = mean(spam$y[!spam$train]!=pred.gam.sub.lin)

c(err.gam.sub.lin,err.glm.sub)
summary(fit.gam.sub.lin)

## Fit a classification tree to all covariates
library(tree)
formula = as.formula(paste("factor(y)",paste(nam,collapse="+"), sep="~"))
fit.t1 <- tree(formula,data=spam,subset=spam$train,split="deviance",control=tree.control(mincut=2,minsize=4,nobs=sum(spam$train)))
plot(fit.t1)
text(fit.t1, digits=2, pretty=3,cex=0.5)
pred.t1 = predict(fit.t1,spam[!spam$train,],type="class")
err.t1 = mean(spam$y[!spam$train]!=pred.t1)
err.t1

# Prune the tree with 10-fold cross validation
fit.p1.CV <- cv.tree(fit.t1,,FUN=prune.tree,K=10)
plot(fit.p1.CV) # No gain by pruning
best.size <- fit.p1.CV$size[which.min(fit.p1.CV$dev)]
best.size

# Using variables selected in textbook
ind = c(5:8,16:17,25:27,37,45:46,52:53,56:57)
formula = as.formula(paste("factor(y)",paste(nam[ind],collapse="+"), sep="~"))
fit.t2 <- tree(formula,data=spam,subset=spam$train,split="deviance",control=tree.control(mincut=2,minsize=4,nobs=sum(spam$train)))
plot(fit.t2)
text(fit.t2, digits=2, pretty=3,cex=0.5)
pred.t2 = predict(fit.t2,spam[!spam$train,],type="class")
err.t2 = mean(spam$y[!spam$train]!=pred.t2)
err.t2 #the tree t1 is better

# Prune the tree with 10-fold cross validation
fit.p2.CV <- cv.tree(fit.t2,,FUN=prune.tree,K=10)
plot(fit.p2) # No gain by pruning
best.size <- fit.p2.CV$size[which.min(fit.p2.CV$dev)]
best.size


## Compare prediction performance of logistic regression, logistic GAM and 
## classification tree with ROC-curves and the area under the curve (AUC)
library(ROCR)
pred.glm.sub <- predict(fit.glm.sub,newdata=spam[!spam$train,],type="response")
prediction.glm.sub <- prediction(pred.glm.sub,spam$y[!spam$train]) #prediction object with predicted probabilities and true ys
perf.glm.sub <- performance(prediction.glm.sub,"tpr","fpr")
pred.gam.sub.lin <- predict(fit.gam.sub.lin,newdata=spam[!spam$train,],type="response")
prediction.gam.sub.lin <- prediction(pred.gam.sub.lin,spam$y[!spam$train]) #prediction object with predicted probabilities and true ys
perf.gam.sub.lin <- performance(prediction.gam.sub.lin,"tpr","fpr")
prediction.t1 <- prediction(as.numeric(pred.t1),spam$y[!spam$train])
perf.t1 <- performance(prediction.t1,"tpr","fpr")

# Plot ROC-curves
plot(perf.glm.sub,col="blue")
plot(perf.gam.sub.lin,col="red",add=TRUE)
plot(perf.t1,col="green",add=TRUE)

# Compute the AUCs
auc.glm.sub <- performance(prediction.glm.sub,"auc")
auc.gam.sub.lin <- performance(prediction.gam.sub.lin,"auc")
auc.t1 <- performance(prediction.t1,"auc")
c(auc.glm.sub@y.values[[1]],auc.gam.sub.lin@y.values[[1]],auc.t1@y.values[[1]])