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

## Code for Workshop 4: Predictive Modeling on Data with Severe 

## Class Imbalance: Applications on Electronic Health Records.  

## The course was conducted for the  International Conference on 

## Health Policy Statistics (ICHPS) on Wed, Oct 7, from 

## 10:15 AM - 12:15 PM.

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

## Example Data

load("emr.RData")

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

## Training/Test Split

library(caret)

set.seed(1732)

in_train <- createDataPartition(emr$Class, p = .75, list = FALSE)

training <- emr[ in_train,]

testing  <- emr[-in_train,]

mean(training$Class == "event")

mean(testing$Class == "event")

table(training$Class)

table(testing$Class)

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

## Overfitting to the Majority Class

library(partykit)

library(rpart)

rpart_small <- rpart(Class ~ ., data = training,

                    control = rpart.control(cp = 0.0062))

plot(as.party(rpart_small))

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

## Subsampling for class imbalances

## Define the resampling method and how we calculate performance

ctrl <- trainControl(method = "repeatedcv",

                     repeats = 5,

                     classProbs = TRUE,

                     savePredictions = TRUE,

                     summaryFunction = twoClassSummary)

## Tune random forest models over this grid

mtry_grid <- data.frame(mtry = c(1:15, (4:9)*5))

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

## The basic random forest model with no adaptations

set.seed(1537)

rf_mod <- train(Class ~ ., 

                data = training,

                method = "rf",

                metric = "ROC",

                tuneGrid = mtry_grid,

                ntree = 1000,

                trControl = ctrl)

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

## This function is used to take the out of sample predictions and

## create an approximate ROC curve from them

roc_train <- function(object, best_only = TRUE, ...) {

  caret:::requireNamespaceQuietStop("pROC")

  caret:::requireNamespaceQuietStop("plyr")

  if(object$modelType != "Classification")

    stop("ROC curves are only availible for classification models")

  if(!any(names(object$modelInfo) == "levels"))

    stop(paste("The model's code is required to have a 'levels' module.",

               "See http://topepo.github.io/caret/custom_models.html#Components"))

  lvs <- object$modelInfo$levels(object$finalModel)

  if(length(lvs) != 2) 

    stop("ROC curves are only implemented here for two class problems")

  ## check for predictions

  if(is.null(object$pred)) 

    stop(paste("The out of sample predictions are required.",

               "See the `savePredictions` argument of `trainControl`"))

  if(best_only) {

    object$pred <- merge(object$pred, object$bestTune)

  }

  ## find tuning parameter names

  p_names <- as.character(object$modelInfo$parameters$parameter)

  p_combos <- object$pred[, p_names, drop = FALSE]

  ## average probabilities across resamples

  object$pred <- plyr::ddply(.data = object$pred, 

                             .variables = c("obs", "rowIndex", p_names),

                             .fun = function(dat, lvls = lvs) {

                               out <- mean(dat[, lvls[1]])

                               names(out) <- lvls[1]

                               out

                             })

  make_roc <- function(x, lvls = lvs, nms = NULL, ...) {

    out <- pROC::roc(response = x$obs,

                     predictor = x[, lvls[1]],

                     levels = rev(lvls))

    out$model_param <- x[1,nms,drop = FALSE]

    out

  }

  out <- plyr::dlply(.data = object$pred, 

                     .variables = p_names,

                     .fun = make_roc,

                     lvls = lvs,

                     nms = p_names)

  if(length(out) == 1)  out <- out[[1]]

  out

}

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

## Some plots of the data

ggplot(rf_mod)

plot(roc_train(rf_mod), 

     legacy.axes = TRUE,

     print.thres = .5,

     print.thres.pattern="   <- default %.1f threshold")

plot(roc_train(rf_mod), 

     legacy.axes = TRUE,

     print.thres.pattern = "Cutoff: %.2f (Sp = %.2f, Sn = %.2f)",

     print.thres = "best")

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

## Internal down-sampling

set.seed(1537)

rf_down_int <- train(Class ~ ., 

                     data = training,

                     method = "rf",

                     metric = "ROC",

                     strata = training$Class,

                     sampsize = rep(sum(training$Class == "event"), 2),

                     ntree = 1000,

                     tuneGrid = mtry_grid,

                     trControl = ctrl)

ggplot(rf_mod$results, aes(x = mtry, y = ROC)) + geom_point() + geom_line() + 

  geom_point(data = rf_down_int$results, aes(x = mtry, y = ROC), col = mod_cols[2]) + 

  geom_line(data = rf_down_int$results, aes(x = mtry, y = ROC), col = mod_cols[2]) + 

  theme_bw() + 

  xlab("#Randomly Selected Predictors") + 

  ylab("ROC (Repeated Cross-Validation)")

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

## External down-sampling

ctrl$sampling <- "down"

set.seed(1537)

rf_down <- train(Class ~ ., 

                 data = training,

                 method = "rf",

                 metric = "ROC",

                 tuneGrid = mtry_grid,

                 ntree = 1000,

                 trControl = ctrl)

geom_point(data = rf_down$results, aes(x = mtry, y = ROC), col = mod_cols[1]) + 

  geom_line(data = rf_down$results, aes(x = mtry, y = ROC), col = mod_cols[1]) + 

  theme_bw() + 

  xlab("#Randomly Selected Predictors") + 

  ylab("ROC (Repeated Cross-Validation)")

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

## Up-sampling

ctrl$sampling <- "up"

set.seed(1537)

rf_up <- train(Class ~ ., 

               data = training,

               method = "rf",

               tuneGrid = mtry_grid,

               ntree = 1000,

               metric = "ROC",

               trControl = ctrl)

ggplot(rf_mod$results, aes(x = mtry, y = ROC)) + geom_point() + geom_line() + 

  geom_point(data = rf_up$results, aes(x = mtry, y = ROC), col = mod_cols[3]) + 

  geom_line(data = rf_up$results, aes(x = mtry, y = ROC), col = mod_cols[3]) + 

  theme_bw() + 

  xlab("#Randomly Selected Predictors") + 

  ylab("ROC (Repeated Cross-Validation)")

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

## Up-sampling done **wrong**

ctrl2 <- trainControl(method = "repeatedcv",

                      repeats = 5,

                      classProbs = TRUE,

                      savePredictions = TRUE,

                      summaryFunction = twoClassSummary)

upped <- upSample(x = training[, -1], y = training$Class)

set.seed(1537)

rf_wrong <- train(Class ~ ., 

                  data = upped,

                  method = "rf",

                  tuneGrid = mtry_grid,

                  ntree = 1000,

                  metric = "ROC",

                  trControl = ctrl2)

ggplot(rf_mod$results, aes(x = mtry, y = ROC)) + geom_point() + geom_line() + 

  geom_point(data = rf_wrong$results, aes(x = mtry, y = ROC), col = mod_cols[3]) + 

  geom_line(data = rf_wrong$results, aes(x = mtry, y = ROC), col = mod_cols[3]) + 

  theme_bw() + 

  xlab("#Randomly Selected Predictors") + 

  ylab("ROC (Repeated Cross-Validation)")

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

## SMOTE 

ctrl$sampling <- "smote"

set.seed(1537)

rf_smote <- train(Class ~ ., 

                  data = training,

                  method = "rf",

                  tuneGrid = mtry_grid,

                  ntree = 1000,

                  metric = "ROC",

                  trControl = ctrl)

ggplot(rf_mod$results, aes(x = mtry, y = ROC)) + geom_point() + geom_line() + 

  geom_point(data = rf_smote$results, aes(x = mtry, y = ROC), col = mod_cols[4]) + 

  geom_line(data = rf_smote$results, aes(x = mtry, y = ROC), col = mod_cols[4]) + 

  theme_bw() + 

  xlab("#Randomly Selected Predictors") + 

  ylab("ROC (Repeated Cross-Validation)")

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

## Make code to measure performance for cost-sensitive learning

fourStats <- function (data, lev = levels(data$obs), model = NULL) {

  accKapp <- postResample(data[, "pred"], data[, "obs"])

  out <- c(accKapp,

           sensitivity(data[, "pred"], data[, "obs"], lev[1]),

           specificity(data[, "pred"], data[, "obs"], lev[2]))

  names(out)[3:4] <- c("Sens", "Spec")

  out

}

ctrl_cost <- trainControl(method = "repeatedcv",

                          repeats = 5,

                          classProbs = FALSE,

                          savePredictions = TRUE,

                          summaryFunction = fourStats)

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

## Setup a custom tuning grid by first fitting a rpart model and

## getting the unique Cp values

rpart_init <- rpart(Class ~ ., data = training, cp = 0)$cptable

cost_grid <- expand.grid(cp = rpart_init[, "CP"],

                         Cost = 1:5)

set.seed(1537)

rpart_costs <- train(Class ~ ., data = training, 

                     method = "rpartCost",

                     tuneGrid = cost_grid,

                     metric = "Kappa",

                     trControl = ctrl_cost)

ggplot(rpart_costs) + 

  scale_x_log10(breaks = 10^pretty(log10(rpart_costs$results$cp), n = 5)) + 

  theme(legend.position = "top")

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

## C5.0 with costs

cost_grid <- expand.grid(trials = 1:3,

                         winnow = FALSE,

                         model = "tree",

                         cost = seq(1, 10, by = .25))

set.seed(1537)

c5_costs <- train(Class ~ ., data = training, 

                  method = "C5.0Cost",

                  tuneGrid = cost_grid,

                  metric = "Kappa",

                  trControl = ctrl_cost)

c5_costs_res <- subset(c5_costs$results, trials <= 3)

c5_costs_res$trials <- factor(c5_costs_res$trials)

ggplot(c5_costs_res, aes(x = cost, y = Kappa, group = trials)) +

  geom_point(aes(color = trials)) + 

  geom_line(aes(color = trials)) + 

  ylab("Kappa (Repeated Cross-Validation)")+ 

  theme(legend.position = "top")