30-sml.Rmd

# Supervised Learning

## Introduction

In **supervised learning** (SML), the learning algorithm is presented
with labelled example inputs, where the labels indicate the desired
output. SML itself is composed of **classification**, where the output
is qualitative, and **regression**, where the output is quantitative.

When two sets of labels, or classes, are available, one speaks of
**binary classification**. A classical example thereof is labelling an
email as *spam* or *not spam*. When more classes are to be learnt, one
speaks of a **multi-class problem**, such as annotation of a new *Iris*
example as being from the *setosa*, *versicolor* or *virginica*
species. In these cases, the output is a single label (of one of the
anticipated classes). If multiple labels may be assigned to each
example, one speaks of **multi-label classification**.

## Preview

To start this chapter, let's use a simple, but useful classification
algorithm, k-nearest neighbours (kNN) to classify the *iris*
flowers. We will use the `knn` function from the `r CRANpkg("class")`
package.

K-nearest neighbours works by directly measuring the (Euclidean)
distance between observations and inferring the class of unlabelled data
from the class of its nearest neighbours. In the figure below, the
unlabelled instances *1* and *2* will be assigned classes *c1* (blue)
and *c2* (red) as their closest neighbours are red and blue,
respectively.

```{r knnex, echo=FALSE, fig.cap="Schematic illustrating the k nearest neighbors algorithm."}
p1 <- c(0, 0)
p2 <- c(0.7, 0.5)
x1 <- rbind(c(0.2, 0.2),
            c(-0.3, -0.8),
            c(-0.2, 1.3))
x2 <- rbind(c(1, 1),
            c(0.5, 0.7))
x3 <- c(1.5, -.9)
x <- rbind(p1, p2, x1, x2, x3)
col <- c("black", "black",
         rep("steelblue", 3),
         rep("red", 2),
         "darkgreen")

plot(x, pch = 19, col = col,
     cex = 5, xlab = "", ylab = "",
     xaxt = "n", yaxt = "n")
grid()
text(p1[1], p1[2], "1", col = "white", cex = 2)
text(p2[1], p2[2], "2", col = "white", cex = 2)
for (i in 1:3)
    segments(p1[1], p1[2],
             x1[i, 1], x1[i, 2],
             lty = "dotted",
             col = "steelblue")
segments(p2[1], p2[2],
         x1[1, 1], x1[1, 2],
         lty = "dotted",
         col = "steelblue")
for (i in 1:2)
    segments(p2[1], p2[2],
             x2[i, 1], x2[i, 2],
             lty = "dotted",
             col = "red")
legend("topright",
       legend = expression(c[1], c[2], c[3]),
       pch = 19,
       col = c("steelblue", "red", "darkgreen"),
       cex = 2,
       bty = "n")
```

Typically in machine learning, there are two clear steps, where one
first **trains** a model and then uses the model to **predict** new
outputs (class labels in this case). In the kNN, these two steps are
combined into a single function call to `knn`.

Lets draw a set of 50 random iris observations to train the model and
predict the species of another set of 50 randomly chosen flowers. The
`knn` function takes the training data, the new data (to be inferred)
and the labels of the training data, and returns (by default) the
predicted class.

```{r knn1}
set.seed(12L)
tr <- sample(150, 50)
nw <- sample(150, 50)
library("class")
knnres <- knn(iris[tr, -5], iris[nw, -5], iris$Species[tr])
head(knnres)
```

We can now compare the observed kNN-predicted class and the expected
known outcome and calculate the overall accuracy of our model.

```{r knn1acc}
table(knnres, iris$Species[nw])
mean(knnres == iris$Species[nw])
```

We have omitted an important argument from `knn`, which is the
parameter *k* of the classifier. This parameter defines how many
nearest neighbours will be considered to assign a class to a new
unlabelled observation. From the arguments of the function,

```{r knnargs}
args(knn)
```

we see that the default value is 1. But is this a good value? Wouldn't
we prefer to look at more neighbours and infer the new class using a
vote based on more labels?

> Challenge
>
> Repeat the kNN classification above by using another value of k, and
> compare the accuracy of this new model to the one above. Make sure
> to use the same `tr` and `nw` training and new data to avoid any
> biases in the comparison.

<details>
```{r knn5}
knnres5 <- knn(iris[tr, -5], iris[nw, -5], iris$Species[tr], k = 5)
mean(knnres5 == iris$Species[nw])
table(knnres5, knnres)
```
</details>

> Challenge
>
> Rerun the kNN classifier with a value of *k* > 1, and specify `prob
> = TRUE` to obtain the proportion of the votes for the winning class.

<details>
```{r knn5prob}
knnres5prob <- knn(iris[tr, -5], iris[nw, -5], iris$Species[tr], k = 5, prob = TRUE)
table(attr(knnres5prob, "prob"))
```
</details>


This introductory example leads to two important and related questions
that we need to consider:

- How can we do a good job in training and testing data? In the
  example above, we choose random training and new data.

- How can we estimate our model parameters (*k* in the example above)
  so as to obtain good classification accuracy?

## Model performance

### In-sample and out-of-sample error

In supervised machine learning, we have a desired output and thus know
precisely what is to be computed. It thus becomes possible to directly
evaluate a model using a quantifiable and objective metric. For
regression, we will use the **root mean squared error** (RMSE), which
is what linear regression (`lm` in R) seeks to minimise. For
classification, we will use **model prediction accuracy**.

Typically, we won't want to calculate any of these metrics using
observations that were also used to calculate the model. This
approach, called **in-sample error** leads to optimistic assessment of
our model. Indeed, the model has already *seen* these data upon
construction, and is considered optimised for these observations
in particular; it is said to **over-fit** the data. We prefer to
calculate an **out-of-sample error**, on new data, to gain a better
idea of how to model performs on unseen data, and estimate how well
the model **generalises**.

In this course, we will focus on the `r CRANpkg("caret")` package for
Classification And REgression Training (see also
https://topepo.github.io/caret/index.html). It provides a common and
consistent interface to many, often repetitive, tasks in supervised
learning.


```{r caretlib}
library("caret")
```

The code chunk below uses the `lm` function to model the price of
round cut diamonds and then predicts the price of these very same
diamonds with the `predict` function.


```{r}
data(diamonds)
model <- lm(price ~ ., diamonds)
p <- predict(model, diamonds)
```

> Challenge
>
> Calculate the root mean squared error for the prediction above

<details>
```{r}
## Error on prediction
error <- p - diamonds$price
rmse_in <- sqrt(mean(error^2)) ## in-sample RMSE
rmse_in
```
</details>

Let's now repeat the exercise above, but by calculating the
out-of-sample RMSE. We prepare a 80/20 split of the data and use
80% to fit our model, and predict the target variable (this is called the
**training data**), the price, on the 20% of unseen data (the **testing
data**).

> Challenge
>
> 1. Let's create a **random** 80/20 split to define the test and
>    train subsets.
> 2. Train a regression model on the training data.
> 3. Test the model on the testing data.
> 4. Calculating the out-of-sample RMSE.

<details>
```{r}
set.seed(42)
ntest <- nrow(diamonds) * 0.80
test <- sample(nrow(diamonds), ntest)
model <- lm(price ~ ., data = diamonds[test, ])
p <- predict(model, diamonds[-test, ])
error <- p - diamonds$price[-test]
rmse_out <- sqrt(mean(error^2)) ## out-of-sample RMSE
rmse_out
```
</details>

The values for the out-of-sample RMSE will vary depending on what
exact split was used. The diamonds is a rather extensive data set, and
thus even when building our model using a subset of the available data
(80% above), we manage to generate a model with a low RMSE, and
possibly lower than the in-sample error.

When dealing with datasets of smaller sizes, however, the presence of
a single outlier in the train and test data split can substantially
influence the model and the RMSE. We can't rely on such an approach and
need a more robust one where we can generate and use multiple,
different train/test sets to sample a set of RMSEs, leading to a
better estimate of the out-of-sample RMSE.

### Cross-validation

Instead of doing a single training/testing split, we can systematise
this process, produce multiple, different out-of-sample train/test
splits, that will lead to a better estimate of the out-of-sample RMSE.

The figure below illustrates the cross validation procedure, creating
3 folds. One would typically do a 10-fold cross validation (if the
size of the data permits it). We split the data into 3 *random* and
complementary folds, so that each data point appears exactly once in
each fold. This leads to a total test set size that is identical to
the size of the full dataset but is composed of out-of-sample
predictions.

![Schematic of 3-fold cross validation producing three training (blue) and testing (white) splits.](./figure/xval.png)

After cross-validation, all models used within each fold are
discarded, and a new model is built using the whole dataset, with the
best model parameter(s), i.e those that generalised over all folds.

This makes cross-validation quite time consuming, as it takes *x+1*
(where *x* in the number of cross-validation folds) times as long as
fitting a single model, but is essential.

Note that it is important to maintain the class proportions within the
different folds, i.e. respect the proportion of the different classes
in the original data. This is also taken care when using the
`r CRANpkg("caret")` package.


The procedure of creating folds and training the models is handled by
the `train` function in `r CRANpkg("caret")`. Below, we apply it to
the diamond price example that we used when introducing the model
performance.

- We start by setting a random seed to be able to reproduce the example.
- We specify the method (the learning algorithm) we want to use. Here,
  we use `"lm"`, but, as we will see later, there are many others to
  choose from[^1].
- We then set the out-of-sample training procedure to 10-fold cross
  validation (`method = "cv"` and `number = 10`). To simplify the
  output in the material for better readability, we set the verbosity
  flag to `FALSE`, but it is useful to set it to `TRUE` in interactive
  mode.

[^1]: There are exactly `r length(names(getModelInfo()))` available
    methods. See
    http://topepo.github.io/caret/train-models-by-tag.html for
    details.

```{r trxval}
set.seed(42)
model <- train(price ~ ., diamonds,
               method = "lm",
               trControl = trainControl(method = "cv",
                                        number = 10,
                                        verboseIter = FALSE))
model
```

Once we have trained our model, we can directly use this `train`
object as input to the `predict` method:

```{r texval}
p <- predict(model, diamonds)
error <- p - diamonds$price
rmse_xval <- sqrt(mean(error^2)) ## xval RMSE
rmse_xval
```

> Challenge
>
> Train a linear model using 10-fold cross-validation and then use it
> to predict the median value of owner-occupied homes in Boston from
> the `Boston` dataset as described above. Then calculate the RMSE.


<details>
```{r xvalsol}
library("MASS")
data(Boston)
model <- train(medv ~ .,
               Boston,
               method = "lm",
               trControl = trainControl(method = "cv",
                                        number = 10))
model
p <- predict(model, Boston)
sqrt(mean((p - Boston$medv)^2))
```
</details>

## Classification performance

Above, we have used the RMSE to assess the performance of our
regression model. When using a classification algorithm, we want to
assess its accuracy to do so.

### Confusion matrix

Instead of calculating an error between predicted value and known
value, in classification we will directly compare the predicted
class matches with the known label. To do so, rather than calculating the
mean accuracy as we did above, in the introductory kNN example, we can
calculate a **confusion matrix**.

A confusion matrix contrasts predictions to actual results. Correct
results are *true positives* (TP) and *true negatives* (TN) are found
along the diagonal. All other cells indicate false results, i.e *false
negatives* (FN) and *false positives* (FP).

```{r cmat, echo=FALSE}
cmat <- data.frame(c("TP", "FN"),
                   c("FP", "TN"))
rownames(cmat) <- c("Predicted Yes", "Predicted No")
colnames(cmat) <- c("Reference Yes", "Reference No")
knitr::kable(cmat)
```

The values that populate this table will depend on the cutoff that
we set to define whether the classifier should predict *Yes* or
*No*. Intuitively, we might want to use 0.5 as a threshold, and assign
every result with a probability > 0.5 to *Yes* and *No* otherwise.

Let's experiment with this using the `Sonar` dataset, and see if we
can differentiate mines from rocks using a logistic classification
model use the `glm` function from the `r CRANpkg("stats")` package.

```{r sonarex0}
library("mlbench")
data(Sonar)
## 60/40 split
tr <- sample(nrow(Sonar), round(nrow(Sonar) * 0.6))
train <- Sonar[tr, ]
test <- Sonar[-tr, ]
```

```{r sonarex1, warning = FALSE}
model <- glm(Class ~ ., data = train, family = "binomial")
p <- predict(model, test, type = "response")
summary(p)
```

```{r sonarex2}
cl <- ifelse(p > 0.5, "M", "R")
table(cl, test$Class)
```

The caret package offers its own, more informative function to
calculate a confusion matrix:

```{r soncmat}
confusionMatrix(factor(cl), test$Class)
```

We get, among others

- the accuracy: $\frac{TP + TN}{TP + TN + FP + FN}$
- the sensitivity (recall, TP rate): $\frac{TP}{TP + FN}$
- the specificity: $\frac{TN}{TN + FP}$
- positive predictive value (precision): $\frac{TP}{TP + FP}$
- negative predictive value: $\frac{TN}{TN + FN}$
- FP rate (fall-out): $\frac{FP}{FP + TN}$

> Challenge
>
> Compare the model accuracy (or any other metric) using thresholds of
> 0.1 and 0.9.

<details>
```{r confmatsol}
confusionMatrix(factor(ifelse(p > 0.9, "M", "R")), test$Class)
confusionMatrix(factor(ifelse(p > 0.1, "M", "R")), test$Class)
```
</details>

### Receiver operating characteristic (ROC) curve

There is no reason to use 0.5 as a threshold. One could use a low
threshold to catch more mines with less certainty or or higher
threshold to catch fewer mines with more certainty.

This illustrates the need to adequately balance TP and FP rates. We
need to have a way to do a cost-benefit analysis, and the solution
will often depend on the question/problem.

One solution would be to try with different classification
thresholds. Instead of inspecting numerous confusion matrices, it is
possible to automate the calculation of the TP and FP rates at each
threshold and visualise all results along a ROC curve.

This can be done with the `colAUC` function from the
`r CRANpkg("caTools")` package:


```{r}
caTools::colAUC(p, test[["Class"]], plotROC = TRUE)
```

- x: FP rate (1 - specificity)
- y: TP rate (sensitivity)
- each point along the curve represents a confusion matrix for a given
  threshold

In addition, the `colAUC` function returns the area under the curve
(AUC) model accuracy metric. This is single number metric, summarising
the model performance along all possible thresholds:

- an AUC of 0.5 corresponds to a random model
- values > 0.5 do better than a random guess
- a value of 1 represents a perfect model
- a value 0 represents a model that is always wrong

### AUC in `caret`

When using `r CRANpkg("caret")`'s `trainControl` function to train a
model, we can set it so that it computes the ROC and AUC properties
for us.

```{r trctrlroc, warning = FALSE}
## Create trainControl object: myControl
myControl <- trainControl(
    method = "cv", ## cross validation
    number = 10,   ## 10-fold
    summaryFunction = twoClassSummary, ## NEW
    classProbs = TRUE, # IMPORTANT
    verboseIter = FALSE
)
## Train glm with custom trainControl: model
model <- train(Class ~ ., Sonar,
               method = "glm", ## to use glm's logistic regression
               trControl = myControl)

## Print model to console
print(model)
```

> Challenge
>
> Define a `train` object that uses the AUC and 10-fold cross
> validation to classify the Sonar data using a logistic regression,
> as demonstrated above.

## Random forest

Random forest models are accurate and non-linear models and robust to
over-fitting and hence quite popular. They however require
hyperparameters to be tuned manually, like the value *k* in the
example above.

Building a random forest starts by generating a high number of
individual decision trees. A single decision tree isn't very accurate,
but many different trees built using different inputs (with
bootstrapped inputs, features and observations) enable us to explore a
broad search space and, once combined, produce accurate models, a
technique called *bootstrap aggregation* or *bagging*.

### Decision trees

A great advantage of decision trees is that they make a complex
decision simpler by breaking it down into smaller, simpler decisions
using a divide-and-conquer strategy. They basically identify a set of
if-else conditions that split the data according to the value of the
features.


```{r rpart, fig.cap="Descision tree with its if-else conditions"}
library("rpart") ## recursive partitioning
m <- rpart(Class ~ ., data = Sonar,
           method = "class")
library("rpart.plot")
rpart.plot(m)
p <- predict(m, Sonar, type = "class")
table(p, Sonar$Class)
```

Decision trees choose splits based on most homogeneous partitions, and
lead to smaller and more homogeneous partitions over their iterations.

An issue with single decision trees is that they can grow, and become
large and complex with many branches, which corresponds to
over-fitting. Over-fitting models noise, rather than general patterns in
the data, focusing on subtle patterns (outliers) that won't
generalise.

To avoid over-fitting, individual decision trees are pruned. Pruning
can happen as a pre-condition when growing the tree, or afterwards, by
pruning a large tree.

- *Pre-pruning*: stop growing process, i.e stops divide-and-conquer
  after a certain number of iterations (grows tree to a certain
  predefined level), or requires a minimum number of observations in
  each mode to allow splitting.

- *Post-pruning*: grow a large and complex tree, and reduce its size;
  nodes and branches that have a negligible effect on the
  classification accuracy are removed.

### Training a random forest

Let's return to random forests and train a model using the `train`
function from `r CRANpkg("caret")`:

```{r loadrange, echo=FALSE, message=FALSE}
suppressPackageStartupMessages(library("ranger"))
```

```{r rftrain, fig.cap="", cache=TRUE}
set.seed(12)
model <- train(Class ~ .,
               data = Sonar,
               method = "ranger")
print(model)
```

```{r rfplotmodel, fig.cap=""}
plot(model)
```

The main hyperparameter is *mtry*, i.e. the number of randomly selected
variables used at each split. Two variables produce random models, while
hundreds of variables tend to be less random, but risk
over-fitting. The `caret` package can automate the tuning of the hyperparameter using a
**grid search**, which can be parametrised by setting `tuneLength`
(that sets the number of hyperparameter values to test) or directly
defining the `tuneGrid` (the hyperparameter values), which requires
knowledge of the model.

```{r tuneLength, eval = FALSE}
model <- train(Class ~ .,
               data = Sonar,
               method = "ranger",
               tuneLength = 5)
```

```{r tuneGrid, fig.cap="", cache=TRUE}
set.seed(42)
myGrid <- expand.grid(mtry = c(5, 10, 20, 40, 60),
                      splitrule = c("gini", "extratrees"),
                      min.node.size = 1) ## Minimal node size; default 1 for classification
model <- train(Class ~ .,
               data = Sonar,
               method = "ranger",
               tuneGrid = myGrid,
               trControl = trainControl(method = "cv",
                                       number = 5,
                                       verboseIter = FALSE))
print(model)
plot(model)
```

> Challenge
>
> Experiment with training a random forest model as described above,
> by using 5-fold cross validation, and setting a `tuneLength` of 5.

<details>
```{r rftrainsol, cache=TRUE, fig.cap=""}
set.seed(42)
model <- train(Class ~ .,
               data = Sonar,
               method = "ranger",
               tuneLength = 5,
               trControl = trainControl(method = "cv",
                                        number = 5,
                                        verboseIter = FALSE))
plot(model)
```
</details>

## Data pre-processing

### Missing values

Real datasets often come with missing values. In R, these should be
encoded using `NA`. There are basically two approaches to deal with
such cases.

- Drop the observations with missing values, or, if one feature
  contains a very high proportion of NAs, drop the feature
  altogether. These approaches are only applicable when the proportion
  of missing values is relatively small. Otherwise, it could lead to
  losing too much data.

- Impute (replace) missing values.

Data imputation can however have critical consequences depending on the
proportion of missing values and their nature. From a statistical
point of view, missing values are classified as *missing completely at
random* (MCAR), *missing at random* (MAR) or *missing not at random*
(MNAR), and the type of the missing values will influence the
efficiency of the imputation method.

The figure below shows how different imputation methods perform
depending on the proportion and nature of missing values
(from [Lazar *et al.*](https://www.ncbi.nlm.nih.gov/pubmed/26906401),
on quantitative proteomics data).

![Normalised RMSE (RMSE-observation standard deviation ration) describing the effect of different imputation methods depending on the nature and proportion of the missing values: kNN (a), SVDimpute (b), MLE (c), MinDet (d), and MinProb (e).](./figure/imp.png)


Let's start by simulating a dataset containing missing values using
the `mtcars` dataset. Below, we will want to predict the `mpg`
variable using `cyl`, `disp`, and `hp`, with the latter containing 10
missing values.

```{r makena}
data(mtcars)
mtcars[sample(nrow(mtcars), 10), "hp"] <- NA
Y <- mtcars$mpg    ## target variable
X <- mtcars[, 2:4] ## predictors
```

If we now wanted to train a model (using the non-formula interface):

```{r rflib, echo=FALSE}
suppressPackageStartupMessages(library("randomForest"))

```

```{r trainna, warning=FALSE}
try(train(X, Y))
```

(Note that the occurrence of the error will depend on the model
chosen.)

We could perform imputation manually, but `r CRANpkg("caret")`
provides a whole range of pre-processing methods, including imputation
methods, that can directly be passed when training the model.

### Median imputation

Imputation using median of features. This method works well if the
data are missing at random.

```{r, eval=TRUE}
train(X, Y, preProcess = "medianImpute")
```

Imputing using caret also allows us to optimise the imputation based on
the cross validation splits, as `train` will do median imputation
inside each fold.

### kNN imputation

If there is a systematic bias in the missing values, then median
imputation is known to produce incorrect results. kNN imputation will
impute missing values using other, similar non-missing rows. The
default value is 5.

```{r, eval=TRUE}
train(X, Y, preProcess = "knnImpute")
```

## Scaling and centering

We have seen in the *Unsupervised learning* chapter how data at
different scales can substantially disrupt a learning
algorithm. Scaling (division by the standard deviation) and centering
(subtraction of the mean) can also be applied directly during model
training by setting. Note that they are set to be applied by default
prior to training.

```{r, eval=FALSE}
train(X, Y, preProcess = "scale")
train(X, Y, preProcess = "center")
```

As we have discussed in the section on Principal Component
Analysis, PCA can be used as pre-processing method, generating a set
of high-variance and perpendicular predictors, preventing
collinearity.

```{r, eval=FALSE}
train(X, Y, preProcess = "pca")
```

### Multiple pre-processing methods

It is possible to chain multiple processing methods: imputation,
center, scale, pca.

```{r, eval=TRUE, warning=FALSE}
train(X, Y, preProcess = c("knnImpute", "center", "scale", "pca"))
```

The pre-processing methods above represent a classical order of
operations, starting with data imputation to remove missing values,
then centering and scaling, prior to PCA.

<!-- ## Low information predictors -->

<!-- To remove, for example constant or random variables, or variables with -->
<!-- low variance. -->

<!-- ```{r, eval=FALSE} -->
<!-- train(X, Y, preProcess = c("zv", ...))  ## remove constant (zero-variance) columns -->
<!-- train(X, Y, preProcess = c("nzv", ...)) ## nearly constant columns -->
<!-- ``` -->

For further details, see `?preProcess`.

## Model selection

In this final section, we are going to compare different predictive
models and choose the best one using the tools presented in the
previous sections.

To to so, we are going to first create a set of common training
controller object with the same train/test folds and model evaluation
metrics that we will re-use. This is important to guarantee fair
comparison between the different models.

For this section, we are going to use the `churn` data. Below, we see
that about 15% of the customers churn. It is important to maintain
this proportion in all of the folds.

```{r churndata}
library("C50")
data(churn)
table(churnTrain$churn)/nrow(churnTrain)
```

Previously, when creating a train control object, we specified the
method as `"cv"` and the number of folds. Now, as we want the same
folds to be re-used over multiple model training rounds, we are going
to pass the train/test splits directly. These splits are created with
the `createFolds` function, which creates a list (here of length 5)
containing the element indices for each fold.

```{r createFolds}
myFolds <- createFolds(churnTrain$churn, k = 5)
str(myFolds)
```

> Challenge
>
> Verify that the folds maintain the proportion of yes/no results.

<details>
```{r foldprop}
sapply(myFolds, function(i) {
    table(churnTrain$churn[i])/length(i)
})
```
</details>

We can now a train control object to be reused consistently for
different model trainings.

```{r trctrol}
myControl <- trainControl(
  summaryFunction = twoClassSummary,
  classProb = TRUE,
  verboseIter = FALSE,
  savePredictions = TRUE,
  index = myFolds
)
```

### `glmnet` model

The `glmnet` is a linear model with built-in variable selection and
coefficient regularisation.

```{r glmnetmodel, fig.cap=""}
glm_model <- train(churn ~ .,
                   churnTrain,
                   metric = "ROC",
                   method = "glmnet",
                   tuneGrid = expand.grid(
                       alpha = 0:1,
                       lambda = 0:10/10),
                   trControl = myControl)
print(glm_model)
plot(glm_model)
```


Below, we are going to repeat this same modelling with a variety of
different classifiers, some of which we haven't looked at. This
illustrates another advantage of of using **meta-packages** such as
`r CRANpkg("caret")`, that provide a consistant interface to different
backends (in this case for machine learning). Once we have mastered
the interface, it becomes easy to apply it to a new backend.

Note that some of the model training below will take some time to run,
depending on the tuning parameter settings.

### random forest model

> Challenge
>
> Apply a random forest model, making sure you reuse the same train
> control object.

<details>
```{r rfmodel, cache=TRUE, fig.cap=""}
rf_model <- train(churn ~ .,
                  churnTrain,
                  metric = "ROC",
                  method = "ranger",
                  tuneGrid = expand.grid(
                      mtry = c(2, 5, 10, 19),
                      splitrule = c("gini", "extratrees"),
                      min.node.size = 1),
                  trControl = myControl)
print(rf_model)
plot(rf_model)
```
</details>

### kNN model

> Challenge
>
> Apply a kNN model, making sure you reuse the same train
> control object.

<details>
```{r knnmodel, cache=TRUE, fig.cap=""}
knn_model <- train(churn ~ .,
                   churnTrain,
                   metric = "ROC",
                   method = "knn",
                   tuneLength = 20,
                   trControl = myControl)
print(knn_model)
plot(knn_model)
```
</details>

### Support vector machine model

> Challenge
>
> Apply a svm model, making sure you reuse the same train control
> object. Hint: Look at `names(getModelInfo())` for all possible model
> names.

<details>
```{r svmmodel, cache=TRUE, fig.cap=""}
svm_model <- train(churn ~ .,
                   churnTrain,
                   metric = "ROC",
                   method = "svmRadial",
                   tuneLength = 10,
                   trControl = myControl)
print(svm_model)
plot(svm_model)
```
</details>

### Naive Bayes

> Challenge
>
> Apply a naive Bayes model, making sure you reuse the same train
> control object.

<details>
```{r nbmodel, fig.cap=""}
nb_model <- train(churn ~ .,
                  churnTrain,
                  metric = "ROC",
                  method = "naive_bayes",
                  trControl = myControl)

print(nb_model)
plot(nb_model)
```
</details>

### Comparing models

We can now use the `caret::resamples` function that will compare the
models and pick the one with the highest AUC and lowest AUC standard
deviation.

```{r resamples}
model_list <- list(glmmet = glm_model,
                   rf = rf_model,
                   knn = knn_model,
                   svm = svm_model,
                   nb = nb_model)
resamp <- resamples(model_list)
resamp
summary(resamp)
```

```{r plotresam, fig.cap = "Comparing distributions of AUC values for various models."}
lattice::bwplot(resamp, metric = "ROC")
```

### Pre-processing

The random forest appears to be the best one. This might be related to
its ability to cope well with different types of input and require
little pre-processing.

> Challenge
>
> If you haven't done so, consider pre-processing the data prior to
> training for a model that didn't perform well and assess whether
> pre-processing affected the modelling.

<details>
```{r svmmodel2, cache=TRUE, fig.cap=""}
svm_model1 <- train(churn ~ .,
                    churnTrain,
                    metric = "ROC",
                    method = "svmRadial",
                    tuneLength = 10,
                    trControl = myControl)

svm_model2 <- train(churn ~ .,
                    churnTrain[, c(2, 6:20)],
                    metric = "ROC",
                    method = "svmRadial",
                    preProcess = c("scale", "center", "pca"),
                    tuneLength = 10,
                    trControl = myControl)

model_list <- list(svm1 = svm_model1,
                   svm2 = svm_model2)
resamp <- resamples(model_list)
summary(resamp)
bwplot(resamp, metric = "ROC")
```
</details>

### Predict using the best model

> Challenge
>
> Choose the best model using the `resamples` function and comparing
> the results and apply it to predict the `churnTest` labels.

<details>
```{r}
p <- predict(rf_model, churnTest)
confusionMatrix(p, churnTest$churn)
```
</details>