# =============================================
# Orange juice data
# =============================================
#
# Summary: 
# These are records of 1070 purchases where the customer either 
# purchased Citrus Hill or Minute Maid Orange Juice. A number of 
# characteristics of the customer and product are recorded.
#
# Data description: 
# The data set consists of the following 18 variables
# (the first colums contains the number of the individual):
# Response variable:
#  Purchase: A factor with levels CH and MM indicating whether 
#            the customer purchased Citrus Hill or Minute 
#            Maid Orange Juice
# Covariates:
#  WeekofPurchase: Week of purchase
#  StoreID: Store ID (1,2,3 or 7)
#  PriceCH: Price charged for CH
#  PriceMM: Price charged for MM
#  DiscCH: Discount offered for CH
#  DiscMM: Discount offered for MM
#  SpecialCH: Indicator of special on CH
#  SpecialMM: Indicator of special on MM
#  LoyalCH: Customer brand loyalty for CH
#  SalePriceMM: Sale price for MM
#  SalePriceCH: Sale price for CH
#  PriceDiff: Sale price of MM less sale price of CH
#  Store7: A factor with levels No and Yes indicating 
#          whether the sale is at Store 7
#  PctDiscMM: Percentage discount for MM
#  PctDiscCH: Percentage discount for CH
#  ListPriceDiff: List price of MM less list price of CH
#  STORE: Which of 5 possible stores the sale occured at
#
#======================================================

library(ISLR2)
data(OJ)

# Divide into a training and a test set
set.seed(1001)
N <- nrow(OJ)
train <- sample(1:N,round(N/2))

# Remove obsolete covariates, define factors for categorical covariates and
# convert response to 0-1
oj <- OJ[,-c(4:5,13:14,17:18)]
oj$StoreID <- as.factor(oj$StoreID)
oj$SpecialCH <- as.factor(oj$SpecialCH)
oj$SpecialMM <- as.factor(oj$SpecialMM)
purchase <- rep(0,N)
purchase[OJ$Purchase == "MM"] <- 1
oj$Purchase <- purchase

# Fit a logistic regression model to all the covariates
fit.glm <- glm(Purchase~.,data=oj,subset=train,family=binomial)
summary(fit.glm)
pred.glm <- predict(fit.glm,newdata=oj[-train,],type="response")
yhat.glm <- pred.glm>0.5
err.glm <- mean(oj$Purchase[-train] != yhat.glm)
err.glm
# Compute AUC
library(ROCR)
prediction.glm <- prediction(pred.glm,oj$Purchase[-train])
auc.glm <- performance(prediction.glm,"auc")@y.values[[1]]
auc.glm
# Try with a subset of the covariates
fit.glm.sub <- glm(Purchase~DiscMM+SpecialCH+LoyalCH+SalePriceMM+SalePriceCH+PctDiscMM,data=oj,subset=train,family=binomial)
summary(fit.glm.sub)
pred.glm.sub <- predict(fit.glm.sub,newdata=oj[-train,],type="response")
yhat.glm.sub <- pred.glm.sub>0.5
err.glm.sub <- mean(oj$Purchase[-train] != yhat.glm.sub)
err.glm.sub
prediction.glm.sub <- prediction(pred.glm.sub,oj$Purchase[-train])
auc.glm.sub <- performance(prediction.glm.sub,"auc")@y.values[[1]]
auc.glm.sub

# Fit an LDA model 
library(MASS)
fit.lda <- lda(Purchase~.,data=oj[train,],method="mle")
pred.lda <- predict(fit.lda,newdata=oj[-train,])$posterior[,2]
yhat.lda <- pred.lda>0.5
err.lda <- mean(oj$Purchase[-train] != yhat.lda)
err.lda
prediction.lda <- prediction(pred.lda,oj$Purchase[-train])
auc.lda <- performance(prediction.lda,"auc")@y.values[[1]]
auc.lda

# Fit a QDA model 
fit.qda <- qda(Purchase~.,data=oj[train,],method="mle")
pred.qda <- predict(fit.qda,newdata=oj[-train,])$posterior[,2]
yhat.qda <- pred.qda>0.5
err.qda <- mean(oj$Purchase[-train] != yhat.qda)
err.qda
prediction.qda <- prediction(pred.qda,oj$Purchase[-train])
auc.qda <- performance(prediction.qda,"auc")@y.values[[1]]
auc.qda

# Fit a logistic GAM
library(gam)
fit.gam.1 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+s(DiscMM)+SpecialCH+SpecialMM+s(LoyalCH)+s(SalePriceMM)+s(SalePriceCH)+s(PctDiscMM)+s(PctDiscCH),data=oj,subset=train,family=binomial)
summary(fit.gam.1)
anova(fit.glm,fit.gam.1,test="Chisq")
# See if some covariates can be dropped
# SpecialMM
fit.gam.2 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+s(DiscMM)+SpecialCH+s(LoyalCH)+s(SalePriceMM)+s(SalePriceCH)+s(PctDiscMM)+s(PctDiscCH),data=oj,subset=train,family=binomial)
summary(fit.gam.2)
anova(fit.gam.2,fit.gam.1,test="Chisq")
# PctDiscCH
fit.gam.3 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+s(DiscMM)+SpecialCH+s(LoyalCH)+s(SalePriceMM)+s(SalePriceCH)+s(PctDiscMM),data=oj,subset=train,family=binomial)
summary(fit.gam.3)
anova(fit.gam.3,fit.gam.2,test="Chisq")
# SpecialCH
fit.gam.4 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+s(DiscMM)+s(LoyalCH)+s(SalePriceMM)+s(SalePriceCH)+s(PctDiscMM),data=oj,subset=train,family=binomial)
summary(fit.gam.4)
anova(fit.gam.4,fit.gam.3,test="Chisq")
# Try a linear effect for DiscMM
fit.gam.5 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+DiscMM+s(LoyalCH)+s(SalePriceMM)+s(SalePriceCH)+s(PctDiscMM),data=oj,subset=train,family=binomial)
summary(fit.gam.5)
anova(fit.gam.5,fit.gam.4,test="Chisq")
# Try a linear effect for PctDiscMM
fit.gam.6 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+DiscMM+s(LoyalCH)+s(SalePriceMM)+s(SalePriceCH)+PctDiscMM,data=oj,subset=train,family=binomial)
summary(fit.gam.6)
anova(fit.gam.6,fit.gam.5,test="Chisq")
# Try a linear effect for LoyalCH
fit.gam.7 <- gam(Purchase~s(WeekofPurchase)+StoreID+s(DiscCH)+DiscMM+LoyalCH+s(SalePriceMM)+s(SalePriceCH)+PctDiscMM,data=oj,subset=train,family=binomial)
summary(fit.gam.7)
anova(fit.gam.7,fit.gam.6,test="Chisq")
# Try a linear effect for WeekofPurchase
fit.gam.8 <- gam(Purchase~WeekofPurchase+StoreID+s(DiscCH)+DiscMM+LoyalCH+s(SalePriceMM)+s(SalePriceCH)+PctDiscMM,data=oj,subset=train,family=binomial)
summary(fit.gam.8)
anova(fit.gam.8,fit.gam.7,test="Chisq")
# Try a linear effect for DiscCH
fit.gam.9 <- gam(Purchase~WeekofPurchase+StoreID+DiscCH+DiscMM+LoyalCH+s(SalePriceMM)+s(SalePriceCH)+PctDiscMM,data=oj,subset=train,family=binomial)
summary(fit.gam.9)
anova(fit.gam.9,fit.gam.8,test="Chisq")
# Try a linear effect for SalePriceCH
fit.gam.10 <- gam(Purchase~WeekofPurchase+StoreID+DiscCH+DiscMM+LoyalCH+s(SalePriceMM)+SalePriceCH+PctDiscMM,data=oj,subset=train,family=binomial)
summary(fit.gam.10)
anova(fit.gam.10,fit.gam.9,test="Chisq")
pred.gam <- predict(fit.gam.10,newdata=oj[-train,],type="response")
yhat.gam <- pred.gam>0.5
err.gam <- mean(oj$Purchase[-train] != yhat.gam)
err.gam
prediction.gam <- prediction(pred.gam,oj$Purchase[-train])
auc.gam <- performance(prediction.gam,"auc")@y.values[[1]]
auc.gam

# Fit a classification tree
library(tree)
cov.names <- names(oj)[-1] 
formula <- as.formula(paste(paste("factor(",names(oj)[1],")",sep=""),paste(cov.names,collapse="+"), sep="~"))
fit.t1 <- tree(formula,data=oj,subset=train,split="deviance",control=tree.control(mincut=2,minsize=4,nobs=length(train)))
plot(fit.t1)
text(fit.t1, digits=2, pretty=3,cex=0.5)
pred.t1 <- predict(fit.t1,oj[-train,],type="vector")[,2]
yhat.t1 <- pred.t1>0.5
err.t1 <- mean(oj$Purchase[-train] != yhat.t1)
err.t1
prediction.t1 <- prediction(pred.t1,oj$Purchase[-train])
auc.t1 <- performance(prediction.t1,"auc")@y.values[[1]]
auc.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) 
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=oj[train,],best=best.size) 
plot(fit.p1.best)
text(fit.p1.best, digits=2, pretty=3,cex=0.5)
pred.p1.best <- predict(fit.p1.best,oj[-train,],type="vector")[,2]
yhat.p1.best <- pred.p1.best>0.5
err.p1.best <- mean(oj$Purchase[-train] != yhat.p1.best)
err.p1.best
prediction.p1.best <- prediction(pred.p1.best,oj$Purchase[-train])
auc.p1.best <- performance(prediction.p1.best,"auc")@y.values[[1]]
auc.p1.best

# Fit a neural network with one hidden layer
library(nnet)
# Standardise the numerical covariates
num.ind <- c(2,4:5,9:12)
oj.stand <- oj
oj.stand[,num.ind] <- scale(oj[,num.ind])
m <- 2:15
err.nnet.size <- rep(NA,length(m))
for(i in 1:length(m))
{
  print(i)
  fit.nnet.size <- nnet(as.factor(Purchase)~.,data=oj.stand[train,],size=m[i],decay=0,MaxNWts=10000,maxit=500)
  pred <- predict(fit.nnet.size,newdata=oj.stand[-train,],type="class")
  err.nnet.size[i] <- mean(oj.stand$Purchase[-train] != pred)
} 
err.nnet.size
decay <- seq(0.5,5,length=10)
m <- 8
err.nnet.decay <- rep(NA,length(decay))
for(i in 1:length(decay))
{
  print(i)
  fit.nnet.decay <- nnet(as.factor(Purchase)~.,data=oj.stand[train,],size=m,decay=decay[i],MaxNWts=10000,maxit=500)
  pred <- predict(fit.nnet.decay,newdata=oj.stand[-train,],type="class")
  err.nnet.decay[i] <- mean(oj.stand$Purchase[-train] != pred)
}  
err.nnet.decay
# Refit best model
fit.nnet <- nnet(as.factor(Purchase)~.,data=oj.stand[train,],size=m,decay=decay[2],MaxNWts=10000,maxit=500)
pred.nnet <- predict(fit.nnet,newdata=oj.stand[-train,],type="raw")
yhat.nnet <- pred.nnet>0.5
err.nnet <- mean(oj$Purchase[-train] != yhat.nnet)
err.nnet
prediction.nnet <- prediction(pred.nnet,oj$Purchase[-train])
auc.nnet <- performance(prediction.nnet,"auc")@y.values[[1]]
auc.nnet

# Fit gradient boosted trees
library(xgboost)
x.tr <- matrix(unlist(oj[train,-1]),nrow=length(train),ncol=11)
y.tr <- unlist(oj$Purchase[train])
x.te <- matrix(unlist(oj[-train,-1]),nrow=N-length(train),ncol=11)
y.te <- unlist(oj$Purchase[-train])
# Use tree depth 1 and 10-fold cross-validation to select M and the learning rate is set to 0.1
fit.boost.1.cv <- xgb.cv(data=xgb.DMatrix(data=x.tr,label=y.tr),nrounds = 1000,eta=0.1,nfold = 10,metrics=list("error"),max_depth=1,objective="binary:logistic")
best.boost.1 <- which.min(fit.boost.1.cv$evaluation_log$test_error_mean)
fit.boost.1 <- xgboost(data=x.tr,label=y.tr,max.depth=1,nrounds=best.boost.1,eta=0.1,objective="binary:logistic")
pred.boost.1 <- predict(fit.boost.1,newdata=x.te)
yhat.boost.1 <- pred.boost.1>0.5
err.boost.1 <- mean(y.te != yhat.boost.1)
err.boost.1
prediction.boost.1 <- prediction(pred.boost.1,oj$Purchase[-train])
auc.boost.1 <- performance(prediction.boost.1,"auc")@y.values[[1]]
auc.boost.1
# Use tree depth 2 
fit.boost.2.cv <- xgb.cv(data=xgb.DMatrix(data=x.tr,label=y.tr),nrounds = 1000,eta=0.1,nfold = 10,metrics=list("error"),max_depth=2,objective="binary:logistic")
best.boost.2 <- which.min(fit.boost.2.cv$evaluation_log$test_error_mean)
fit.boost.2 <- xgboost(data=x.tr,label=y.tr,max.depth=2,nrounds=best.boost.2,eta=0.1,objective="binary:logistic")
pred.boost.2 <- predict(fit.boost.2,newdata=x.te)
yhat.boost.2 <- pred.boost.2>0.5
err.boost.2 <- mean(y.te != yhat.boost.2)
err.boost.2
prediction.boost.2 <- prediction(pred.boost.2,oj$Purchase[-train])
auc.boost.2 <- performance(prediction.boost.2,"auc")@y.values[[1]]
auc.boost.2

# Fit random forests
library(randomForest)
x.tr <- matrix(unlist(oj[train,-1]),nrow=length(train),ncol=11)
y.tr <- as.factor(unlist(oj$Purchase[train]))
x.te <- matrix(unlist(oj[-train,-1]),nrow=N-length(train),ncol=11)
y.te <- as.factor(unlist(oj$Purchase[-train]))
# Fit a random forest with m=4 and B=2500
m <- 4
B <- 2500
fit.rf.1 <- randomForest(x=x.tr,y=y.tr,ntree=B,mtry=m)
pred.rf.1 <- predict(fit.rf.1,newdata=x.te,type="prob")[,2] 
yhat.rf.1 <- as.numeric(pred.rf.1>0.5)
err.rf.1 <- mean(y.te != yhat.rf.1)
err.rf.1
prediction.rf.1 <- prediction(pred.rf.1,oj$Purchase[-train])
auc.rf.1 <- performance(prediction.rf.1,"auc")@y.values[[1]]
auc.rf.1
plot(fit.rf.1$err.rate[,1],type="l",ylim=c(0.15,0.23),lwd=2,xlab="Number of trees",ylab="OOB error estimate")
# Fit a random forest with m=6 and B=2500
m <- 6
B <- 2500
fit.rf.2 <- randomForest(x=x.tr,y=y.tr,ntree=B,mtry=m)
pred.rf.2 <- predict(fit.rf.2,newdata=x.te,type="prob")[,2] 
yhat.rf.2 <- as.numeric(pred.rf.2>0.5)
err.rf.2 <- mean(y.te != yhat.rf.2)
err.rf.2
prediction.rf.2 <- prediction(pred.rf.2,oj$Purchase[-train])
auc.rf.2 <- performance(prediction.rf.2,"auc")@y.values[[1]]
auc.rf.2
lines(fit.rf.2$err.rate[,1],lwd=2,col=2)
# Fit a random forest with m=8 and B=2500
m <- 8
B <- 2500
fit.rf.3 <- randomForest(x=x.tr,y=y.tr,ntree=B,mtry=m)
pred.rf.3 <- predict(fit.rf.3,newdata=x.te,type="prob")[,2] 
yhat.rf.3 <- as.numeric(pred.rf.3>0.5)
err.rf.3 <- mean(y.te != yhat.rf.3)
err.rf.3
prediction.rf.3 <- prediction(pred.rf.3,oj$Purchase[-train])
auc.rf.2 <- performance(prediction.rf.3,"auc")@y.values[[1]]
auc.rf.3
lines(fit.rf.3$err.rate[,1],lwd=2,col=3) 

## Clustering (forgetting about the known label)

# K-Means clustering with K=2, dropping the categorical covariates and the label
k.means.fit <- kmeans(oj.stand[,num.ind], 2) 
k.means.fit$centers
k.means.fit$size
# Compare the K-means clusters to the true groups
table(oj.stand[,1],k.means.fit$cluster)
# Compute and plot the within group sum of squares for K=1,...,nc
wssplot <- function(data, nc=15, seed=1234){
  wss <- (nrow(data)-1)*sum(apply(data,2,var))
  for (i in 2:nc){
    set.seed(seed)
    wss[i] <- sum(kmeans(data, centers=i)$withinss)}
  plot(1:nc, wss, type="b", xlab="Number of Clusters",
       ylab="Within group sum of squares")}
wssplot(oj.stand[,num.ind])
clusplot(oj.stand[,num.ind], k.means.fit$cluster, main='2D representation of the Cluster solution',
         color=TRUE, shade=TRUE, labels=2, lines=0)

# K-Means clustering with K=4
k.means.fit.2 <- kmeans(oj.stand[,num.ind], 4) 
# Centroids:
k.means.fit.2$centers
# Cluster size:
k.means.fit.2$size
table(oj.stand[,1],k.means.fit.2$cluster)
library(cluster)
clusplot(oj.stand[,num.ind], k.means.fit.2$cluster, main='2D representation of the Cluster solution',
         color=TRUE, shade=TRUE, labels=2, lines=0)

# K-medoids clustering with K=2
k.medoids.fit <- pam(oj.stand[,-1],2)
table(oj.stand[,1],k.medoids.fit$clustering)
clusplot(oj.stand[,num.ind], k.medoids.fit$cluster, main='2D representation of the Cluster solution',
         color=TRUE, shade=TRUE, labels=2, lines=0)

# K-medoids clustering with K=4
k.medoids.fit.2 <- pam(oj.stand[,-1],4)
table(oj.stand[,1],k.medoids.fit.2$clustering)
clusplot(oj.stand[,num.ind], k.medoids.fit.2$cluster, main='2D representation of the Cluster solution',
         color=TRUE, shade=TRUE, labels=2, lines=0)