---
title: "Extreme Class Imbalance"
subtitle: "Caravan Contest Dataset"
author: "Ian McLeod"
date: "April 7, 2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(fig.pos = 'H')
knitr::opts_chunk$set(message=FALSE)
knitr::opts_chunk$set(warning=FALSE)
library(caret)
library("AppliedPredictiveModeling") #install.packages("AppliedPredictiveModeling")
library(mda) #install.packages("mda")
library(bestglm) #mcdonald data
library(glmnet) #required by caret
library(Matrix) #required by caret when using glmnet
library(knitr)
library(kableExtra) #install.packages("kableExtra")
library(tidyverse)
library("ISLR") #install.packages("ISLR")
library("pROC") #install.packages("pROC")
library("RRF")  #install.packages("RRF")
library("doParallel")
#might need this at end to close down
#stopImplicitCluster()
#get Train/Test Data
set.seed(757513) #reproducible
data(Caravan) #Remark: not exported
Xy <- Caravan
Xy$Purchase <- factor(Xy$Purchase, levels=c("Yes", "No"))
ind <- createDataPartition(Xy$Purchase, p=0.7)[[1]]
XyTr <- Xy[ind,]
XyTe <- Xy[-ind,]
yTr <- XyTr[, ncol(XyTr)]
yTe <- XyTe[, ncol(XyTe)]
XTr <- XyTr[, -ncol(XyTr)]
XTe <- XyTe[, -ncol(XyTe)]
ROCACC <- function(...) c(twoClassSummary(...), defaultSummary(...))
```


#Introduction

Many datasets in medicine and insurance pose difficulties in class prediction 
due to extreme imbalance in the classes.
Chapter 16 of *Applied Predictive Modeling* discusses strategies 
for dealing with this problem.
In the case of binary classification, one strategy is to use ROC Curve 
analysis to select a more useful cutoff than 0.5.

##Poorly Calibrated Probabilities

We noted that probabilities generated using the softmax or 
the logistic function by classifiers such as
Support Vector Machines (SVM) and Multilayer Perceptrons (MLP) 
have a Bayesian interpretation.
But to be useful in practice, there probabilities must be **well-calibrated**.
Fawcett (2006) reports an interesting example 
involving a Naive Bayes (NB) classifier.
As we know the NB classifier makes very strong assumptions 
so it may be not too surprising that
the probabilities it generates may not be realistic.
In Fawcett's example, the NB classifier produced the posterior 
probabilities shown in the Table below:

```{r}
p <- c(0.99999, 0.99999, 0.99993, 0.99986, 0.99964, 0.99955, 0.68139, 0.50961, 0.48880, 0.44951)
Y <- factor(c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0), levels=c(1, 0), labels=c("Yes", "No"))
tb <- list(Y=Y, prediction=p)
kable(tb, caption="Fawcett - Not well-Calibrated Probabilities") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
```
The ROC curve is shown below. We see that AUC=1 so this is a perfect
classifier!

```{r rocFaw}
ROCFawcett <- roc(Y, p)
plot(ROCFawcett)
```
But the default method of predicting produces the confusion matrix:

```{r DEFAULT}
YH <- factor(ifelse(p>0.5, "Yes", "No"), levels=c("Yes", "No"))
table(Y, YH, dnn=list("Truth", "Predicted"))
```
So the mis-classification rate is 2/10 which is not perfect!
Also specificity=50% and sensitivity=100%.  
However simply adjusting the cutoff results in perfect classification.

```{r ADJ}
YH <- ifelse(p>0.8, "Yes", "No")
table(Y, YH, dnn=list("Truth", "Predicted"))
```

Many ML classifiers such as the Naive Bayes, in the above example, 
produce "posterior probabilities" that are useful for ranking 
and provide a measure of relative plausibility 
but they are not true probabilities!

A formal method of using the ROC curve to choose an optimal
classifier will be discussed later in this note.

##Training and contest test data

The Caravan Dataset was used in data mining contest.
The contest details, results and many submissions are described here:
(http://www.liacs.nl/~putten/library/cc2000/)


The dataset itself is described here:
(http://www.liacs.nl/~putten/library/cc2000/data.html)

**TICDATA2000.txt**: 
Dataset to train and validate prediction models and 
build a description (5822 customer records). 
Each record consists of 86 attributes, 
containing sociodemographic data (attribute 1-43) 
and product ownership (attributes 44-86).
The sociodemographic data is derived from zip codes. 
All customers living in areas with the same zip code have 
the same sociodemographic attributes.  
Attribute 86, "CARAVAN:Number of mobile home policies", is the target variable. 

**TICEVAL2000.txt**: 
Dataset for predictions (4000 customer records). 
It has the same format as TICDATA2000.txt, only the target is missing. 
Participants supposed to predict which customers have purchased insurance.
Accurate prediction is valuable because in actual practice 
customers characteristic may be used to 
predict which customers are likely to buy insurance.

**TICTGTS2000.txt**
Targets for the evaluation set.

The dataset `Caravan` available in the package **ISLR** is 
TICDATA2000 that used for the contest.
There are $n=5822$ observations and $p=85$ inputs.
The output variable `Purchase` is `Yes` or `No` 
according as whether or not customer purchased insurance or not. 
We randomly partition this data into two parts 
one for training and one for testing using 70% for training.

The challenges with the dataset include:

* with 85 predictor, it is relatively high-dimensionsal
* the dependence between inputs and output is weak
* extreme class imbalance since only about 6% of the customers are `Yes`

If only 6% are `Yes` this means that simply predicting `No` 
for all new customers will have an expected accuracy of about 94%.

##Summary of Techniques for Dealing with Extreme Imbalance

* recalibrate the probabilties
  * use left-most rule in pROC
  * select top 10 or 20 etc.
  * adjust priors
  * adjust weights
  * adjust tuning parameters
* stratified sampling with caret::upSample() or downSample()
* instead of Accuracy or AUC, cost-sensitive training with rpart() or C50() or other suitable algorithm
* feature selection may be helpful
* LASSO regression does feature selection and regularization
* simple parsimonious methods such as logistic or LDA may better than more complex methods such as RF

##Logistic Regression and ROC Curve Analysis

We use 7 replications of 10-fold CV, it takes less 
than 30 seconds using 16 cores.

```{r FITLor, cache=TRUE}
set.seed(1441414113)
ncores <- parallel::detectCores()
cl <- makeCluster(ncores)
registerDoParallel(cl)
startTime <- proc.time()[3] #elapsed time
#
ctrl <- trainControl(method="repeatedcv", number=10, repeats=7, classProbs = TRUE, summaryFunction = ROCACC)
ansLor <- train(Purchase ~ ., data = XyTr,
               method = "glm", family=binomial,
               trControl = ctrl,
               metric = "ROC")
totTimeLor <- proc.time()[3] -  startTime #about 23 sec w 16 cores
stopCluster(cl)
ansLor
```
From the summary we see the Specificity and Accuracy are very high
but the Sensitivity, the most important criterion in this
application is extremely low. 
Our main focus is to be able to predict which type of
customers are likely to buy insurance.
The low sensitivity informs us that there is apparently little
predictive power in predictions.

The ROC curve does indicate that the model is capable
of more than random predictions.
Indeed is works well in informing us that very few customers
purchase insurance!

```{r}
yHTr <- predict(ansLor, XyTr, type = "prob")[,1]
yHTe <- predict(ansLor, XyTe, type = "prob")[,1]
LorROC <- roc(yTe, yHTe)
LorROC
plot(LorROC, main="Performance of Logistic Regression on Test Data")
```

```{r}
CM <- confusionMatrix(predict(ansLor, XyTe), yTe)
CM
```
The 'closest.topleft' rule is used to select an optimal cutoff
based on the idea of having balanced good performance for both
sensitivity and specificity. We see that we get about 71%
correct in predicting customers with insurance.

```{r}
rfThresh <- coords(LorROC, x = "best", ret="threshold", best.method="closest.topleft")
plot(LorROC, print.thres = rfThresh, type = "S",
     print.thres.pattern = "%.3f (Spec = %.2f, Sens = %.2f)",
     print.thres.cex = .8, legacy.axes = TRUE)
```

The confusion matrix correspond to the 'closest.topleft' rule
is shown below.

```{r}
yHTeBest <- factor(ifelse(yHTe>rfThresh, "Yes", "No"), levels=c("Yes", "No"))
CMBest <- confusionMatrix(yHTeBest, yTe)
CMBest
```


##Random Forest

We investigate the use of Random Forest for prediction
with the Caravan data.
We use 10-fold cross-validation with 7 replications to
select the tuning parameter **mtry**.
Best Accuracy was achieved with **mtry=15** but better sensivity
was with **mtry=85**.
We could investigate if larger number values of the tuning parameter
continue to improve sensitivity since that is most important
in this application.

```{r FITRF, cache=TRUE}
set.seed(1441414113)
ncores <- parallel::detectCores()
cl <- makeCluster(ncores)
registerDoParallel(cl)
startTime <- proc.time()[3] #elapsed time
#
ctrl <- trainControl(method="repeatedcv", number=10, repeats=7, classProbs = TRUE, summaryFunction = ROCACC)
ansRF <- train(Purchase ~ ., data = XyTr,
               method = "rf",
               trControl = ctrl,
               ntree = 1000,
               tuneLength = 7,
               metric = "ROC")
totTimeRF <- proc.time()[3] -  startTime #about 16 min w 16 cores
stopCluster(cl)
ansRF
```

The performance of RF on the ROC curve is very similar to that
of logistic regression.

```{r}
yHTr <- predict(ansRF, XyTr, type = "prob")[,1]
yHTe <- predict(ansRF, XyTe, type = "prob")[,1]
rfROC <- roc(yTe, yHTe)
rfROC
plot(rfROC, main="Performance of RF on Test Dataset")
```
The sensitivity is slightly better but still very unsatisfactory.

```{r}
CM <- confusionMatrix(predict(ansRF, XyTe), yTe)
CM
```
Using the 'closest.topleft` we achieve much better sensitivity
and even slightly better than logistic regression.

```{r}
rfThresh <- coords(rfROC, x = "best", ret="threshold",
                   best.method="closest.topleft")
plot(rfROC, print.thres = rfThresh, type = "S",
     print.thres.pattern = "%.3f (Spec = %.2f, Sens = %.2f)",
     print.thres.cex = .8, legacy.axes = TRUE)
```


```{r}
yHTeBest <- factor(ifelse(yHTe>rfThresh, "Yes", "No"), levels=c("Yes", "No"))
CMBest <- confusionMatrix(yHTeBest, yTe)
CMBest
```

##Regularized Random Forest

We briefly explored a global regularization method 
but it doesn't seem to help.

method = 'RRFglobal'

Tuning parameters:

* mtry (#Randomly Selected Predictors)
* coefReg (Regularization Value)

Required packages: RRF

We briefly explore this methods using 25 bootstrap iterations with vanilla
default settings for the tuning parameters.

```{r FITRRF, cache=TRUE}
set.seed(1441414113)
ncores <- parallel::detectCores()
cl <- makeCluster(ncores)
registerDoParallel(cl)
startTime <- proc.time()[3] #elapsed time
#
ctrl <- trainControl(method="repeatedcv", number=10, repeats=7, classProbs = TRUE, summaryFunction = ROCACC)
ctrl <- trainControl(classProbs = TRUE, summaryFunction = ROCACC) #vanilla version
ansRRF <- train(Purchase ~ ., data = XyTr,
               method = "RRFglobal",
               trControl = ctrl,
               ntree = 1000,
               #tuneGrid =  expand.grid(mtry=seq(3, 123, 40), coefReg=c(0.01, 0.1, 0.3, 1)),
               metric = "ROC")
totTimeRRF <- proc.time()[3] -  startTime #about 15 min w 16 cores
stopCluster(cl)
ansRRF
```


```{r}
yHTr <- predict(ansRRF, XyTr, type = "prob")[,1]
yHTe <- predict(ansRRF, XyTe, type = "prob")[,1]
rfROC <- roc(yTe, yHTe)
rfROC
plot(rfROC, main="Performance of `RRF` on Test Dataset")
```

```{r}
CM <- confusionMatrix(predict(ansRRF, XyTe), yTe)
CM
```

```{r}
rfThresh <- coords(rfROC, x = "best", ret="threshold",
                   best.method="closest.topleft")
plot(rfROC, print.thres = rfThresh, type = "S",
     print.thres.pattern = "%.3f (Spec = %.2f, Sens = %.2f)",
     print.thres.cex = .8, legacy.axes = TRUE)
```


```{r}
yHTeBest <- factor(ifelse(yHTe>rfThresh, "Yes", "No"), levels=c("Yes", "No"))
CMBest <- confusionMatrix(yHTeBest, yTe)
CMBest
```

##LASSO Logistic Regression

Feature selection may be helpful in understanding this data.
Especially for the second part of the contest.
We show how to use LASSO regression to select the most important
predictors.

```{r LASSO, cache=TRUE}
startTime <- proc.time()[3] #elapsed time
xtr <- model.matrix(Purchase ~ ., data=XyTr)
ansL1 <- glmnet(xtr, yTr, alpha=1, family="binomial")
set.seed(77711153)
ans_cv <- cv.glmnet(xtr, y=yTr, family="binomial", alpha=1)
plot(ans_cv)
lambdaHat <- ans_cv$lambda.1se
abline(v=log(lambdaHat), lwd=2.5, col=alpha("magenta", 0.4))
xte <- model.matrix(Purchase ~ ., data=XyTe)
yHat <- factor(predict(ansL1, newx=xte, type="class", s=lambdaHat)[,1], levels=c("Yes", "No"))
totTimeL1 <- proc.time()[3] -  startTime #about 41 sec w 16 cores
confusionMatrix(yTe, yHat)
```

```{r LASSOSelect}
b<-as.matrix(coef(ans_cv))
tibble(predictor=rownames(b), bHat=c(b)) %>%
  filter(bHat !=0)
```

```{r, eval=FALSE, echo=FALSE}
yHTe <- predict(ansL1, xte, type = "response")[,1]
ROCL1 <- roc(yTe, yHTe)
ROCL1
plot(ROCL1, main="Performance of LASSO on Test Dataset")
```
