03-methods.Rmd

# Experimental Setup and Methods

## Tools Used

The problem itself was solved using the python language, making use of established packages. From the scipy ecosystem: `numpy`, a library for multidimensional arrays and vectorized operations on them [@oliphant2006guide]; `scipy`, a library for scientific computing [@mckinney-proc-scipy-2010]; `pandas`, a library for data analysis [@mckinney2011pandas]; `matplotlib`, for high-quality plots [@hunter2007matplotlib]. Furthermore: the machine learning library`scikit-learn` [@scikit-learn] and the high-level data visualization library `seaborn` [ @waskom2014seaborn].

Calculations were run in the cloud on a linux virtual machine^[For details, refer to: [https://docs.microsoft.com/en-us/azure/machine-learning/data-science-virtual-machine/overview](https://docs.microsoft.com/en-us/azure/machine-learning/data-science-virtual-machine/overview), accessed on 5.6.2019].

Except for the development of a helper package, most programming was performed in interactive Jupyter notebooks [@Kluyver:2016aa].

The report was written in `rmarkdown` [@allaire2016rmarkdown] using `knitr` [@xie2015] and `bookdown` [@xie2016bookdown] to render the document into several output formats.

All work was tracked in version control.


## Data Handling

The learning data set was used to establish the complete analysis pipeline, i.e. preprocessing, imputation and feature selection, model evaluation and -selection and prediction method. The test data was not touched until making the final prediction, subjecting it to the complete pipeline determined on the learning data.

In accordance with recommendations in @friedman2001elements, the learning data set was further split 80/20 into training and validation sets (see Figure \@ref(fig:data-splitting)). The training set was used to train different models while the validation set served to tune hyperparameters. The split was performed using a stratified sampling algorithm to preserve target class frequencies based on the response indicator.

```{r data-splitting, fig.cap="Data set use for training and predictions.", echo=F, results='asis', out.width="50%"}
data_graph <- grViz(
  "
  digraph d_split {
  
    graph [overlap=false, fonsize=18, fontname=Lato, layout=dot, nodesep=1.4]
    
    node [fillcolor=white, shape=circle]
    C [label=<<b>Complete Data Set</b><br />191'779 examples>, width=3.0]
    L [label=<<b>Learning</b><br />95'412 examples>, width=2.0]
    T [label=<<b>Test</b><br />n96'367 examples>, width=2.0]
    LT [label=<<b>Training</b>>, width=1.5]
    LV [label=<<b>Validation</b>>, width=0.5]
    
    edge [color = black, penwidth=0.95]
    C -> L [label='used for\npreprocessing,\n training models']
    C -> T [label='kept back\n for final predicitons']
    L -> LT [label='  80 % of data']
    L -> LV [label='  20 % of data']
  }
  ")

rsvg_pdf(charToRaw(DiagrammeRsvg::export_svg(data_graph)),
         file="figures/methods/data-splitting.pdf")
rsvg_png(charToRaw(DiagrammeRsvg::export_svg(data_graph)),
         file="figures/methods/data-splitting.png")
knitr::include_graphics("figures/methods/data-splitting.png")
```


## Data Preprocessing

The necessary preprocessing was guided by practical necessity (input errors, inconsistent categories), the requirements of the algorithms that were examined (Section \@ref(eval-and-select)) and the requirements set out in the cup documentation:

* Only numeric features (required by some algorithms)
* Imputation of missing values (required by cup documentation)
* Removal of constant and *sparse* features (required by cup documentation)

The transformations were established interactively in Jupyter notebooks. Once finalized, transformations were implemented in the python package `kdd98`^[Available from `r make_github_link("kdd98")`].



### Cleaning

The transformations applied can be studied in the Jupyter notebook *1_Preprocessing.ipynb*^[`r make_github_link("notebooks", "1_Preprocessing.ipynb")`].

The cleaning stage of preprocessing encompassed the following transformations:

* Removing *noise*: Input errors, inconsistent encoding of binary / categorical features
* Dropping constant and sparse (i.e. those where only few examples have a value set) features
* Imputation of values missing at random (MAR)

MAR values in the sense of @rubin1976inference are missing conditionally on other features in the data. For example, there are three related features from the promotion and giving history: *ADATE*, the date of mailing a promotion, *RDATE*, the date of receiving a donation in response to the promotion and *RAMOUNT*, the amount received. For missing *RAMOUNT* values, it can be checked if *RDATE* is non-missing. If *RDATE* is missing, then the example most likely has not donated and *RAMOUNT* can be set to zero. If, on the other hand, both date features have a value, *RAMOUNT* is truly missing.



### Feature Engineering {#methods-feature-engineering}

The transformations applied during feature engineering are described in detail in the Jupyter notebook *2_Feature_Engineering.ipynb*^[`r make_github_link("notebooks", "2_Feature Engineering.ipynb")`]. The result of this step was an all-numeric data set usable for downstream learning.

All non-numeric (i.e. categorical) features were encoded as numeric representations. For ordinal features, manual mappings from alphanumeric levels to integer numbers were specified. For nominal features, two encoding techniques provided in package `categorical-encoding`^[Available at: [https://github.com/scikit-learn-contrib/categorical-encoding/](https://github.com/scikit-learn-contrib/categorical-encoding/), accessed on 5.6.2019] were employed, depending on the number of levels:

* One-hot encoding for $\leq 10$ levels: For each level of a categorical feature, a new feature is created. An additional feature may be added to indicate missing values. Exactly one of these new features is set to $1$, indicating the original level.

* Binary encoding, $> 10$ levels: The levels of the categorical are first transformed to ordinal (i.e. to a sequence of integer numbers). Then, these numbers are taken to the base of 2. (A $5$, for example, becomes `101`). According to the number of levels, new features for the binary digits are created. As an example: To represent 60 levels, 6 features are required ($2^6=64$).

Also, several features were transformed to better usable formats (i.e. converting dates to time deltas, converting zip codes to coordinates). Care was taken to keep the dimensionality of the data set as low as possible. 


### Imputation

Three different approaches were evaluated. The details are shown in Jupyter notebook *4_Imputation.ipynb*^[`r make_github_link("notebooks", "4_Imputation.ipynb")`].

#### K-Nearest Neighbors

kNN imputation uses those examples that are "near" an example with a missing value in feature $f$ to impute the value. In short, the kNN algorithm by @troyanskaya2001missing implemented in `missingpy`^[Available at: [https://github.com/epsilon-machine/missingpy](https://github.com/epsilon-machine/missingpy), accessed on 11.6.2019] works as follows:

1. Construct the distance matrix $D$ with distances between examples.
2. Order all features descending by number of missing values.
3. Starting with the feature with most missing, for each example with a missing value, use the $k$ nearest neighbors' mean or median to impute.

The algorithm runs until all values are imputed.

While very attractive because of the intuitive approach and because it preserves data characteristics (continuous / discrete, categoricals) in the features, the distance matrix is very memory-intensive for large data sets. Also, it is required to remove features with > 80 % missing values first.

#### Iterative imputation

Iterative imputation, implemented in package `fancyimpute`^[Available at: [https://pypi.org/project/fancyimpute/](https://pypi.org/project/fancyimpute/), accessed on 30.06.2019], works similar to the R-package `mice` (see @buuren2011mice). Before imputation, all features have to be transformed to numerical data types and normalized because the underlying model is a linear regresssion.

In short, the algorithm works as follows:

1. Features are ordered by the fraction of missing values
2. Starting with the feature with most missing values, use the other features to build a linear model, using the current feature as the dependent variable and predict missing values.
3. Repeat step 2 until all features are complete
4. Repeat steps 2 -- 3 $n$ times, $n=5$ was chosen

#### Simple Imputation and Categorical Indicator

This approach is straightforward: Numeric features are imputed by their median value to make the imputation robust to skewed distributions.

As the algorithm used (`sklearn.impute.SimpleImputer`) only supports numerical data types, categorical features were treated separately during feature engineering (Section \@ref(methods-feature-engineering)): The one-hot or binary encoded categoricals had one more feature added, indicating missing values.



### Feature Selection {#methods-feature-selection}

One of the biggest caveats in dealing with high-dimensional data is the infamous "Curse of Dimensionality" coined by @bellman1966dynamic. The curse comes from the fact that with an increasing number of dimensions of the feature space, the number of possible combinations grows exponentially. In order to cover all possible combinations with several examples, a huge amount of data would be required as a result. In the area of machine learning, high dimensionality frequently manifests in the form of overfitting, which leads to an unacceptably large generalization error @goodfellow2016deep. @hughes1968mean showed that for a fixed number of examples, model performance first increases with increasing number of dimensions but then decreases again. 

Boruta, introduced by @kursa2010boruta in the form of an R package, performs feature selection by solving the so-called all-relevant feature problem. The algorithm was found to perform very well regarding selection of relevant features in @kursa2011boruta.

In short, the algorithm works as follows:

1. The input matrix $\mathbf{X}$ of dimension $n \times p$ is extended with $p$ so-called *shadow features*. The shadow features are permuted copies of the features in $\mathbf{X}$. They are therefore de-correlated with the target.
2. On the resulting matrix $\mathbf{X^*}$, a random forest classifier is trained and the Z-scores $\frac{\bar{\text{loss}}}{sd}$ for each of the $2p$ features calculated.
3. The highest Z-score among the shadow features $MZSA$ is determined.
4. All original features are compared against $MZSA$ and those features with a higher score selected as important.
5. With the remaining features, a two-sided test for equality of the Z-scores with $MZSA$ is performed and all features with significantly lower score are deemed unimportant.
6. All shadow copies are removed, go to step 1.

The algorithm terminates when all attributes are marked as either important or not important or when the maximum number of iterations is reached.

For this thesis, a python implementation^[Available from [https://github.com/scikit-learn-contrib/boruta_py](https://github.com/scikit-learn-contrib/boruta_py), accessed on 5.6.2019] was used.

The models were all learned on the same data set resulting from the boruta selection.


## Prediction {#methods-prediction}

The desired quantity to estimate is net profit $\hat{\Pi}$. In order to predict this quantity, a two-step procedure was implemented, utilizing the response indicator $\text{TARGET}_B$ and the continuous target *$\text{TARGET}_D$*, respectively. For each step, a model was trained. One is a classifier, predicting $\hat{y}_b$, the probability of donating. The other is a regressor, predicting the donation amount $\hat{y}_d$. The classifier was trained on the complete learning data set, while the regressor was trained on $\mathbf{X}_d = \{x_i|y_{b,i} = 1, i=1 \ldots n\}$ and $y_d = \{y_{d,i}|y_{b,i} = 1, i=1 \ldots n\}$. The non-random sample leads to a conditionally predicted donation amount, thereby introducing bias. A correction for the bias is necessary as a consequence. This approach resembles @heckman1976common's two-stage procedure which is widely used in econometrics. Heckman's procedure was presented for the question of wage offerings for women. Data on wages was only available for working women, thereby introducing bias if wage offerings were only predicted on this sample. The question of whether an example works or not is seen as an unobserved feature. In @heckman1976common, a probit model is used for this first *selection* stage. The inverse Mills ratio of the probit, $\frac{-\phi(\hat{y_i}-\mathbf{X}\beta)}{\Phi(\hat{y_i}-\mathbf{X}\beta}$, is calculated and included in the data for predictions in the second *observation* stage with OLS regression.

### Setup of the Two-Stage Prediction

For the first stage some classifier is used, predicting the probability for example $x_i$ to donate. This is not a probit ($P(Y=1|X) = \Phi(X^T\beta)$), of course. Instead, the resulting distribution depends on the classifier.

\begin{equation}
\hat{y}_b = f(\mathbf{X})
(\#eq:y-b)
\end{equation}

where $\hat{y}_b$ is the vector of predicted probabilities of donating and $f$ is the classifier.

The second stage is performed on $\mathbf{X}_d$ and consists in predicting the donation amount using a regression model:

 \begin{equation}
\hat{y}_{dt} = g(\mathbf{X_d})
(\#eq:y-d)
\end{equation}

where $\mathbf{\hat{y}_{dt}}$ is the vector of predicted donation amounts and $g$ is the conditionally learned regression model. $g$ is learned with Box-Cox transformed target $\mathbf{y_d}$, so $\hat{y}_{dt}$ is also Box-Cox transformed with parameter $\lambda$.

The decision of whether to include example $i$ in the promotion is governed by the following indicator function. In it, $\alpha^*$ accounts for the introduced bias. Every example that has a predicted donation amount of more than the unit cost is included.

\begin{equation}
 \mathbb{1}_{\hat{y}_{i,b} * \exp(\hat{y}_{i,dt}) * \alpha^* > \exp(u_t)}(\hat{y}_{i,dt})
(\#eq:indicator)
\end{equation}

where $\alpha^* \in [0,1]$ is a factor to correct for bias introduced due to the non-randomness of $\mathbf{X}_d$, $\hat{y}_{i,dt}$ is the predicted donation amount, transformed so as to normalize the distribution learned beforehand and $u_t$ is the unit cost, Box-Cox transformed with parameter $\lambda$. The exponential is used to deal with negative values resulting from the Box-Cox transformation.

Finally, the quantity estimated is net profit $\hat{\Pi}$. It is calculated by summing over the product of the indicator function \@ref(eq:indicator) and the net profit for examples $1 \ldots n$. For unseen data, \@eq:pi-alpha is altered by using the estimated net profit $\hat{y}_{d,i}$ in the product.

\begin{equation}
\hat{\Pi}_\alpha = \sum_{i=1}^n  \mathbb{1}_{\hat{y}_{i,b} * \exp(\hat{y}_{i,dt}) * \alpha^* > \exp(u_t)}(\hat{y}_{i,dt})*(y_{i,d} - u)
(\#eq:pi-alpha)
\end{equation}



### Optimization of $\alpha^*$

With equation \@ref(eq:pi-alpha), the estimated profit $\Pi$ was calculated on the learning data for a grid of $\alpha$ values, $\alpha \in [0,1]$. The optimal value was then $\alpha^{*} = \underset{\alpha}{\operatorname{argmax}} f(\alpha)$ where $f$ is a function that was fit to $\Pi$.

For $f$, a cubic spline $s$ was used. $\alpha^*$ was determined as follows:

1. Fit $s(\Pi)$, the cubic spline on the estimated profits for the grid of $\alpha$ values
2. Derive $ds = \frac{\delta}{\delta \alpha} s$
3. Find the finite roots of $ds$, $\alpha_{\text{candidates}}$, representing candidates for $\alpha^*$
4. Determine $\alpha^* = \underset{\alpha}{\operatorname{argmax}} s(\alpha_{\text{candidates}})$


## Model Evaluation and -Selection {#eval-and-select}

Several algorithms (estimators) were trained in parallel. Once good hyperparameters were found using randomized grid search, performance of the estimators was compared using a common metric in order to select the best (see Figure \@ref(fig:evaluation-selection) for a schematic process overview). This was done independently for classifiers (predicting the binary target) and regressors (predicting the continuous target). The process is documented in notebook *6_Model_Evaluation_Selection.ipynb*^[`r make_github_link("notebooks", "6_Model_Evaluation_Selection.ipynb")`].

The pipeline functionality of `scikit-learn` was used to combine preliminary data transformations, i.e. scaling of features (where necessary) and resampling with the estimator-algorithm. This allows to jointly tune hyperparameters for transformations and the estimator.

```{r evaluation-selection, fig.cap="Learning process schematic.", echo=F, out.width="80%"}

workflow_graph <- grViz(
  " 
    digraph workflow {

      labeljust='l';

      graph [overlap=false,
             fontname=Lato,
             fonsize=30,
             layout=dot,
             nodesep=1.1,
             rankdir=LR]

      node[fontname=Lato, fillcolor=white, shape=circle, width=1.3]
      
      scale[group='a', label='Scaler', style='dashed'];
      samp[group='a', label='Resampler', style='dashed'];
      estim[group='a', label='Estimator'];

      pip[group='a', label='Pipeline']
      param[group='a', label='Set Hyperparameter\nDistributions'];
      gs[group='a', label='Randomized\nGrid Search'];
      eval[group='a', label='Evaluation\nof results by metric,\nhyperparameters chosen'];
      select[group='a', label='Select best model'];

      edge[style=solid, color=black, fontname=Lato];
      {scale;samp}->pip [style='dashed'];
      estim -> pip;
      pip -> param;
      param -> gs;
      gs -> eval;
      eval -> param [label='using best\n cv iteration \n for metric'];
      eval -> select [label='hyperparameters\nconverged'];
    }
  ")

rsvg_pdf(charToRaw(DiagrammeRsvg::export_svg(workflow_graph)),
         file="figures/methods/evaluation-process.pdf")
rsvg_png(charToRaw(DiagrammeRsvg::export_svg(workflow_graph)),
         file="figures/methods/evaluation-process.png")
knitr::include_graphics("figures/methods/evaluation-process.png")
```

A common random seed was used for all algorithms and random number generators.

### Evaluation

Randomized grid search was run with 10-fold cross-validation (CV). The best-performing pipeline for each algorithm was stored in a python dictionary during evaluation. The dictionary was persisted to disk and only updated when an algorithm's metric improved. This ensured that the best hyperparameter settings were always retained during the extensive model evaluation phase.

#### Randomized Grid Search
In randomized grid search, [@bergstra2012random], probability distributions for hyperparameter values are specified. The algorithm then runs a defined number (10 were used) of random combinations by sampling from the distributions. Compared to the usual grid search, this can greatly speed up the learning process because good hyperparameter settings are generally identified with less iterations.

After one round of learning, the hyperparameter distributions were adjusted before the next iteration as follows: When the best value was found near the limits of the domain, the distribution was shifted in this direction. For values falling inside the domain, the distribution was narrowed down towards the found value. This procedure was repeated until the hyperparameters converged.

The implementation of randomized grid search in `scikit-learn` returns a summary table with the CV results for each hyperparameter combination. By comparing the mean test score for the chosen metric, the best hyperparameters can be determined.

#### Cross-Validation
CV splits the training data into several *folds* of equal size. The algorithm is trained as many times as there are folds, holding back one of the folds at each training step for validation using some specified performance metric and training with the rest of the data. This procedure enables quantification of the generalization error and the calculation of statistics that indicate the variance of the model. Following guidance in @kohavi1995study, 10-fold CV was used, which trades higher bias for lower variance compared to fewer folds.


### Selection

#### Classifiers

Both recall and F1 were initially selected and calculated during model evaluation. In the end, recall was used to select the best classifier as it produced better estimators for the task. Reasoning for the choice of these metrics is given below.

For classification problems, the confusion matrix (see Figure \@ref(fig:conf-mat-plot)) can be used to construct various performance metrics. A true positive (TP) indicates a correctly predicted 1, a false negative (FN) is a falsely predicted 0, a false positive is a falsely predicted 1 and finally a true negative (TN) is a correctly predicted 0.

For the data analyzed here, *1* means an example has donated, *0* means the example has not donated.

```{python conf-mat, include=F}

cm = np.array([[0.0, 0.0],[0.0, 0.0]])
labs = np.array([["TP", "FN"],["FP", "TN"]])
classes = ["1", "0"]

fig, ax = plt.subplots()
im = ax.imshow(cm, interpolation='nearest', cmap=sns.light_palette(sns.husl_palette(8)[4], n_colors=256, as_cmap=True))
# We want to show all ticks...
ax.set(xticks=np.arange(cm.shape[1]),
       yticks=np.arange(cm.shape[0]),
       # ... and label them with the respective list entries
       xticklabels=classes, yticklabels=classes,
       ylabel='True',
       xlabel='Predicted')

# Loop over data dimensions and create text annotations.
fmt = "s"
for i in range(labs.shape[0]):
    for j in range(labs.shape[1]):
        ax.text(j, i, format(labs[i, j], fmt),
                ha="center", va="center",
                color="black")
save_fig("methods/conf-mat")
plt.close('all')
```

(ref:conf-mat-cap) Definition of the confusion matrix for a two-class problem. 

```{r conf-mat-plot, fig.cap="(ref:conf-mat-cap)", echo=F, out.width="30%"}
knitr::include_graphics("figures/methods/conf-mat.png")
```

The definitions of some often-used metrics are given below^[taken from [https://en.wikipedia.org/wiki/Confusion_matrix](https://en.wikipedia.org/wiki/Confusion_matrix), accessed on 28.05.2019]. The choice of metric depends on the goal of the prediction and the data at hand.

\begin{align*}
\text{Recall / Sensitivity / True Positive Rate TPR} &= \frac{TP}{TP + FN} \\
\text{Specificity / True Negative Rate TNR} &= \frac{TN}{TN+FP} \\ 
\text{Precision / Positive Predictive Value PPV} &= \frac{TP}{TP + FP}\\
\text{Negative Predictive Value NPV} &= \frac{TN}{TN+FN}\\
\text{False Negative Rate FNR} &= \frac{FN}{FN+TP}\\
\text{False Positive Rate FPR} &= \frac{FP}{FP+TN}\\
\text{Accuracy} &= \frac{TP + TN}{TP + FP + FN + TN} \\
\text{F1 score} &= \frac{2TP}{2TP+FP+FN}
\end{align*}

The goal, as mentioned earlier, is to maximize net profit. To achieve this, a balance between predicting as many TP's as possible while keeping the number of FP's low has to be found. One FP costs `r 0.68` $ (sending a letter to a non-donor). Keeping in mind the distribution of $\text{TARGET}_D$ (Section \@ref(targets)), one FN means loosing at least `r 0.32` $ of possible profit (not sending a promotion to a donor, smallest donation amount is 1 $). The expected loss in profit for one FN is even `r 15` $ (corresponding to the mean donation amount), which means that with each TP, the cost of  `r floor(15/0.68)` FP's can be covered on average.

The default metric for classification is accuracy, but in the case of imbalanced targets, it is not a desirable metric (TN is present both in the nominator and denominator, dominating the score). The metrics that could be used beneficially because they involve TP an not TN are F1, recall and precision. Since in precision, the FP of the majority class have the highest influence from among the three, it was discarded from the candidate list.

#### Regressors

For regression, $R^2 = 1- \frac{\sum_{i=1}^n(y_i-\hat{y}_i)^2}{\sum_{i=1}^n(y_i-\bar{y})^2}$ was used, mainly because it is the default metric for regression algorithms in `scikit-learn`. $R^2$ has the drawback of depending on the variance of the data used to fit the model and therefore is different for other data. It was however assumed that because learning and test data have the same generating function, $R^2$ can be used to select a regression model.

### Dealing With Imbalanced Data {#imblearn}

Several different approaches were explored. The following over-/undersampling techniques available in package `imblearn` by @lemaitre2017imblearn were studied:

* Random oversampling of the minority class
* Random undersampling of the majority class
* SMOTE (synthetic minority oversampling technique), variant borderline-1

The random sampling algorithms either draw from the minority class with repetition or draw random samples from the majority class without repetition until the labels are balanced. SMOTE does not sample from the data but instead generates synthetic samples from the minority class, thereby countering the danger of overfitting by learning on a small number of repeated observations (random oversampling) or a small subset of the training data available (random undersampling). The SMOTE variant *borderline-1* chosen generates samples that are close to the optimal decision boundary, where misclassification is likely.

Additional experiments were run with class- and sample-weights set on the data without resampling. For class weights, the ratio of non-donors vs. donors was used. For sample weights, the donation amount was employed, rescaled to the interval $[0,1]$.

### Algorithms

A short introduction of each algorithm is given below. For each algorithm, the hyperparameters that were tuned during learning are given. The choice of algorithms was made so as to cover a wide range of underlying concepts.

#### Random Forest {#methods-rf}

Random forest (RF) belongs to the family of so-called ensemble learners and was introduced by @breiman2001rf. Predictions are made by majority vote of an ensemble of decision trees (CART, @breiman1984classification). The RF can be employed both for regression and classification tasks. RF's are insensitive towards scale differences in the individual features, and, depending on the implementation, can deal with missing values. Another important feature of RF is the assessment of *variable importance* by summing the loss improvement for each split in every tree per feature [@friedman2001elements].

During learning, a random sample of the available features is drawn with replacement for each tree (bagging, or bootstrap aggregating, see @breiman1996bagging), thereby reducing the variance of the ensemble estimator. Furthermore, splits within a tree are determined again on a random subset of the features. These sources of randomness tend to increase bias of the forest, yet the decrease in variance due to averaging through majority vote outweighs the bias increase @friedman2001elements. @breiman2001rf shows that as the forest grows, the generalization error converges almost surely. This means that random forests are insensitive to overfitting and perform better the more trees are grown.

As explained in @friedman2001elements, trees are grown as follows:

For data with n examples and p features $D = \{\{x_i,y_i\}, i=1 \ldots n, x_i=\{x_{i,1}, x_{i,2}, \ldots x_{i,p}\}\}$, the CART algorithm decides on the structure of the tree, the splitting features and the split points. As a result, the data is partitioned into $M$ regions $R_1, R_2, \ldots, R_M$.

**Regression**
The response $\hat{y}$ is modeled as a constant $c_m$ for each region:
\begin{equation}
f(x) = \sum_{m=1}^M c_m\mathbb{1}(x \in R_m)
(\#eq:cart-const)
\end{equation}

Using as loss criterion the sum of squares $\sum_{i=1}^n (y_i - f(x_i))^2$, the best $\hat{c}_m$ is the average of $y_i$ in the region:
\begin{equation}
\hat{c}_m=ave(y_i|x_i \in R_m).
(\#eq:cart-hatc)
\end{equation}

The algorithm greedily decides on the best partition. Starting with all data, a splitting feature $j$ and a split point $s$ is considered, creating two regions $R_1(j,s) = \{X|X_j \leq s\}, R_2(j,s) = \{X|X_j > s\}$. Feature $j$ and split point $s$ are chosen by solving
\begin{equation}
\min_{j,s}\left[\min_{c1} \sum_{x_i \in R_1(j,s)}(y_i-c_1)^2 + min_{c2} \sum_{x_i \in R_2(j,s)}(y_i-c_2)^2\right]
(\#eq:cart-opt)
\end{equation}
with the inner minimization solved using \@ref(eq:cart-hatc).


**Classification**
For a binary classification problem with outcomes $\{0,1\}$, predictions are made through the proportion of the positive class in a region $R_m$ with $N_m$ examples $x_i$ inside, which is given by:
\begin{equation}
\hat{p}_m = \frac{1}{N_m} \sum_{x_i \in R_m} \mathbb{1}(y_i = 1)
(\#eq:cart-class)
\end{equation}

When $\hat{p}_m>0.5$, the positive class is chosen.

The loss is described by impurity. When making a split, the feature $j$ resulting in the highest impurity decrease is selected. Impurity is measured by the *Gini index*. For binary classification: A node $m$, representing region $R_m$ with $N_m$ observations has as proportion of the positive class $\hat{p}_m = \frac{1}{N_m} \sum_{x_i \in R_m} \mathbb{1}(y_i = 1)$. The Gini index is then defined as $2p(1-p)$. So the decision function is:
\begin{equation}
\min_{j,s}(\min_{Gini_l} \hat{p}_l + \min_{Gini_r} \hat{p}_r)
(\#eq:cart-class-dec)
\end{equation}
for regions $R_l, R_r$ below the node.

The `RandomForestClassifier` and `RandomForestRegressor` included in `scikit-learn` were used for learning.

**Hyperparameters**

* `max_depth`, $\{1,2,3, \ldots\}$: depth of the trees, $2^n$ leafs maximum. Controls tree size.
* `min_samples_split`, $\{2,3,4,\ldots\}$: Minimum number of samples required to split a node, controls tree size.
* `max_features`, $\{1, 2, \ldots, m\}$: Maximum number of the $m$ features to consider when searching for a split. @friedman2001elements recommend values in  $m = \{1, 2, \ldots \sqrt{m}\}$, but for high dimensional data with few relevant features, larger $m$ can lead to better results because the probability of including relevant features increases.
* `n_estimators`, $\{1, 2, \ldots\}$: Number of trees to grow. In combination with early stopping, this can be set to a high value since learning will stop when the loss converges.
* `class_weight` $\{\text{balanced}, 1,2,\ldots\}$: Weights on target classes: "balanced" calculates weights according to class frequencies, integer values specify weight on majority class relative to minority


#### Gradient Boosting Machine {#methods-gbm}

The main idea behind boosting is to sequentially train an ensemble of weak learners which on their own are only slightly better than a random decision. The predictions of the individual weak learners are then combined into a majority vote [@kearns1988thoughts].

Gradient boosting machine (GBM) extends on this idea. Like a random forest, GBM learns an ensemble of decision trees. However, trees are learned in an additive manner. At each iteration, the tree that improves the model most (i.e. in the direction of the gradient of the loss function) is added. For this thesis, the package `XGBoost` by @chen2016xgboost was used. They describe the algorithm as follows:

Assume a data set with $n$ examples and $p$ features: $D = \{\{x_i,y_i\}, i=1\ldots n, x_i=\{x_{i,1}, x_{i,2}, \ldots x_{i,p}\}\}$. The implementation uses a tree ensemble with $K$ regression trees to predict the outcome for an example in the data by summing up the weights predicted by each tree:

\begin{equation}
\hat{y}_i = \phi(x_i) = \sum_{k=1}^K f_k(x_i), f_k \in F
(\#eq:gbm-ensemble)
\end{equation}

where $F = \{f(x) = w_{q(x)}\} (q: \mathbb{R}^m \rightarrow T, w \in \mathbb{R}^T)$ is the space of regression trees. $T$ is the number of leaves in a tree, $q$ is the structure of each tree, mapping an example to the corresponding leaf index. Each tree $f_k$ has an independent structure $q$ and weights $w$ at the terminal leafs.

For growing the trees in $F$, the following loss function is minimized:

\begin{equation}
L(\phi) = \sum_{i=1}^n l(y_i, \hat{y_i}) + \sum_k^K \Omega(f_k)
(\#eq:gbm-loss)
\end{equation}

Here, $l$ is a differentiable, convex loss function that measures the difference between predictions and true values. Since $l$ is convex, a global minimum is guaranteed. $\Omega(f) = \gamma T + \frac{1}{2}\lambda||w||^2$ is a penalty on the complexity of the trees to counter over-fitting. The algorithm thus features integrated regularization. 

Now, at each iteration $t$, the tree $f_t$ that improves the model most is added. For this, $f_t(\mathbf{x_i})$ is added to the predictions at $t-1$ and compute the loss:

\begin{equation}
L^{(t)} = \sum_{i=1}^n l(y_i, \hat{y_i}^{t-1} + f_t(\mathbf{x_i})) +\Omega(f_t)
(\#eq:gbm-iterate)
\end{equation}

To find the best $f_t$ to add, the gradients finally come into play. With $g_i$ and $h_i$ the first- and second-order gradient statistics of $l$, the loss function becomes:

\begin{equation}
\tilde{L}^{(t)} = \sum_{i=1}^n l(y_i, \hat{y_i}^{t-1} + g_if_t(x_i) + h_i f_t^2(x_i)) +\Omega(f_t)
(\#eq:gbm-grad)
\end{equation}

**Hyperparameters**

* `learning_rate`, $[0, 1]$: *Shrinkage*, decreases step size for the gradient descent when $\eta < 1.0$, helping convergence. The number of estimators $f_k$ has to be increased for small learning rates in order for the algorithm to converge.
* `min_child_weight`, $[0, \inf)$: Minimum sum of weights of the hessian in a node. When close to zero, the node is pure. Controls regularization.
* `subsample`, $[0, 1]$: A random sample from the $n$ examples of size $s, s < n$ is drawn for each iteration, countering overfitting and speeding up learning.
* `colsample_by_tree`, $[0, 1]$: A random sample of the $m$ features is drawn for growing each tree.
* `n_iter_no_change`, $\{0,1,2,3, \ldots\}$: Early stopping. Based on an evaluation set, learning stops when no improvement on the performance metric (misclassification error was chosen) is made for a fixed number of steps.



#### GLMnet

The GLMnet is an implementation of a generalized linear model (GLM) with penalized maximum likelihood by @hastie2014glmnet. Regularization is achieved through $L^2$ (ridge)  and $L^1$ (lasso) penalties or their combination known as elastic net [@zou2005regularization].

The loss functions are described in @hastie2014glmnet. For the binary classification task at hand, a logistic regression was performed. The logistic regression model for a two-class response $G = \{0,1\}$ with target $y_i = \mathbb{1}(g_i=1)$ is:

$P(G=1|X=x) = \frac{e^{\beta_0+\beta^Tx}}{1+e^{\beta_0+\beta^Tx}}$ or, in the log-odds transformation: $log\frac{P(G=1|X=x)}{P(G=0|X=x)}=\beta_0+\beta^Tx$.

The GLMnet loss function for logistic regression:

\begin{equation}
\min_{\beta_0,\beta} \frac{1}{N} \sum_{i=1}^{N} w_i l(y_i,\beta_0+\beta^T x_i) + \lambda\left[(1-\alpha)||\beta||_2^2/2 + \alpha ||\beta||_1\right]
(\#eq:glmnet-logit)
\end{equation}

where $w_i$ are individual sample weights, $l(\cdot)$ the (negative) log-likelihood of the parameters $\mathbf{\beta}$ given the data, $\lambda$ the amount of penalization and $\alpha \in [0,1]$ the elastic net parameter, for $\alpha=0$ pure ridge and for $\alpha=1$ pure lasso.

For the regression task, a gaussian family model was used, having loss function:

\begin{equation}
\min_{(\beta_0, \beta) \in \mathbb{R}^{p+1}}\frac{1}{2N} \sum_{i=1}^N (y_i -\beta_0-x_i^T \beta)^2+\lambda \left[ (1-\alpha)||\beta||_2^2/2 + \alpha||\beta||_1\right]
(\#eq:glmnet-gaussian)
\end{equation}

where $l()$ the (negative) log-likelihood of the parameters $\mathbf{\beta}$ given the data, $\lambda$ the amount of penalization and $\alpha \in [0,1]$ the elastic net parameter, for $\alpha=0$ pure ridge and for $\alpha=1$ pure lasso.

For learning, GLMnet evaluates many different $\lambda$ values for a given $\alpha$ through cross validation. Because GLMnet is sensitive to scale differences in the features, input data (features and target) should be transformed to mean zero and unit variance.

**Hyperparameters**

* `n_splits`, $\{3,4,5, \ldots\}$: Number of CV-splits. Typical values are 3, 5 and 10.
* $\alpha$, $[0, 1]$: parametrizes the elastic net. For $\alpha = 0$ pure ridge, for $\alpha = 1$ pure lasso.
* `scoring`: Scoring method for cross-validation (log-loss, classification error, accuracy, precision, recall, average precision, roc-auc)



#### Multilayer Perceptron

The multilayer perceptron (MLP) is a so-called feed-forward neural network. The network consists of at least three layers: an input layer, an arbitrary number of hidden layers and an output layer. Each layer is made up of units. The term feed-forward means that information flows from the input layer through intermediary steps and then to the output. The goal is to approximate the function $f^*$. For a classifier, $y = f^*(x)$ maps an example $x$ to a category $y$. A feed-forward network defines a mapping $y=f(x, w)$ and learns the weights $w$ by approximating the function $f$ [@goodfellow2016deep].

For a binary classification problem on a dataset with $n$ examples and $p$ features $D = \{\{x_i, y_i\}\}, x_i = \{x_{i,1},x_{i,2},\ldots,x_{i,p}\}, y_i \in \{0,1\}, i = 1 \ldots n\}$, the input layer has $p$ units, the output layer has $1$ unit. The hidden layers each have an arbitrary number of *hidden* units.

Each unit, except for the input layer, consists of a perceptron, which is in effect a linear model with some non-linear activation function applied:

\begin{equation}
y = \phi(w^Tx+b)
(\#eq:perceptron)
\end{equation}

where $\phi$ is a non-linear activation function, $w$ is the vector of weights, $x$ is the vector of inputs and b is the bias. For $\phi$, typical functions are the hyperbolic tangens $tanh(\cdot)$, the logistic sigmoid $\sigma(x) = \frac{1}{1+e^{-x}}$, or, more recently, the rectified linear unit $relu(x) = max(0,x)$ [@hahnloser2000digital; @goodfellow2016deep].

During learning, the training examples are fed to the network sequentially. For each example, the prediction error is calculated using a loss function, which is typically the negative log-likelihood.

Then, the partial derivatives of the loss function with respect to the weights are computed for each unit and the parameters updated using stochastic gradient descent. This process is called back-propagation.

The complete network is a chain of functions, for the network trained here, it is:
\begin{equation}
\mathbf{y} = \phi^{(3)}(\mathbf{W}^{(3)T} \mathbf{\phi}^{(2)}(\mathbf{W}^{(2)T} \mathbf{\phi}^{(1)}(\mathbf{W}^{(1)T} \mathbf{x} + \mathbf{b}^{(1)}) + \mathbf{b}^{(2)}) + \mathbf{b}^{(3)})
\end{equation}
where  $x$ is the vector of input features, $y$ is the vector of outputs, $\mathbf{W}^{(1)}, \mathbf{W}^{(2)}, \mathbf{W}^{(3)}$ are the weight matrices for each layer and $b^{(1)}, b^{(2)}, b^{(3)}$ are the bias vectors for each layer and $\phi^{(1)}, \phi^{(2)}, \phi^{(3)}$ are the sets of perceptrons in the corresponding layer.

(ref:mlp-g-legend) Neural network topology used. Two hidden layers $\mathbf{h^{(1)}, h^{(2)}}$ are contained. $\mathbf{b^{(2)}, b^{(3)}}$ and $\mathbf{b^{(4)}}$ are the bias vectors for the respective layers.

```{r mlp-graph, fig.cap="(ref:mlp-g-legend)", echo=F, results="asis"}
data_graph <- grViz(
  "
  digraph G {

    graph [overlap=false,
           fontname=Lato,
           fontsize=28,
           nodesep=0.5,
           rankdir=LR,
           ranksep= 1.1;]

    node [fontname=Lato, fonsize=24, fillcolor=white, shape=circle, height=0.8]

    splines=false;
    edge [color=red, style=invis]
    rankdir = LR;
    {
      rank=same;
        labels [group='l', label='.', style = invis];
        network [group='n', label='.', style = invis];
        labels -> network
    }

    subgraph layer_labels{
      l0 [shape=plaintext, label='layer 1\n(input)', group='l'];
      l1 [shape=plaintext, label='layer 2\n(hidden)', group='l'];
      l2 [shape=plaintext, label='layer 3\n(hidden)', group='l'];
      l3 [shape=plaintext, label='layer 4\n(output)', group='l'];
      l0->l1->l2->l3;
    }

    labels->l0;

    subgraph network_nodes{
      node [shape=circle];
      b1 [label=<b<sup>(2)</sup>>, group='n'];
      b2 [label=<b<sup>(3)</sup>>, group='n'];
      b3 [label=<b<sup>(4)</sup>>, group='n'];

      x1 [label=<x<sub>1</sub>>, group='n'];
      x2 [label=<x<sub>2</sub>>, group='n'];
      x3 [shape=plain, label='...', fontsize=36, group='n']
      x4 [label=<x<sub>n</sub>>, group='n'];
      h12 [label=<h<sub>1</sub><sup>(2)</sup>>, group='n'];
      h22 [label=<h<sub>2</sub><sup>(2)</sup>>, group='n'];
      h32 [label=<h<sub>3</sub><sup>(2)</sup>>, group='n'];
      h42 [shape=plain, label='...', fontsize=36, group='n'];
      h52 [label=<h<sub>n</sub><sup>(2)</sup>>, group='n'];
      h13 [label=<h<sub>1</sub><sup>(3)</sup>>, group='n'];
      h23 [label=<h<sub>2</sub><sup>(3)</sup>>, group='n'];
      h33 [label=<h<sub>3</sub><sup>(3)</sup>>, group='n'];
      h43 [shape=plain, label='...', fontsize=36, group='n']
      h53 [label=<h<sub>n</sub><sup>(3)</sup>>, group='n'];
      O1 [label=y, group='n'];

      h12->h13;
      {
          rank=same;
          l0->b1->x1->x2->x3->x4;
      }
      {
          rank=same;
          l1->b2->h12->h22->h32->h42->h52;
      }
      {
          rank=same;
          l2->b3->h13->h23->h33->h43->h53;
      }
      {
        rank=same;
        l3->O1
      }

      edge[color=black, style=solid, tailport=e, headport=w, arrowhead=vee];
      
      {x1;x2;x3;x4} -> {h12;h22;h32;h42;h52};
      {h12;h22;h32;h42;h52} -> {h13;h23;h33;h43;h53};
      {h13;h23;h33;h43;h53} -> {O1};

      edge[style=dashed]
      b1 -> {h12;h22;h32;h42;h52};
      b2 -> {h13;h23;h33;h43;h53};
      b3 -> O1;
    }

    network -> x1
  }
")
rsvg_pdf(charToRaw(DiagrammeRsvg::export_svg(data_graph)),
         file="figures/methods/mlp-structure.pdf")
rsvg_png(charToRaw(DiagrammeRsvg::export_svg(data_graph)),
         file="figures/methods/mlp-structure.png")
knitr::include_graphics("figures/methods/mlp-structure.png")
```

For this thesis, the scikit-learn implementation `sklearn.neural_network.MLPClassifier` was used with default settings (rectified linear unit as activation function, *adam* (stochastic gradient descent) solver). The network topology was determined by treating it as a hyperparameter. Best results were achieved with a network featuring two hidden layers with 28 hidden units each. The network is shown symbolically in Figure \@ref(fig:mlp-graph).

**Hyperparameters**

* `hidden_layer_sizes`: Tuples, i.e. (100,50,) for a network with 2 hidden layers of 100 and 50 units and one output unit.
* $\alpha$, [0,1]: $L^2$ regularization amount
* `learning_rate_init` $>0$: Controls step size for the solver during back-propagation (regularization).


#### Support Vector Machine

In the case of binary classification, support vector classifiers (SVC) find a separating hyperplane $\{x: f(x)=x^T\beta+\beta_0 = 0\}$ with classification rule $G(x)=sign(x^T\beta+\beta_0)$. The best hyperplane is such that a margin $M$ defined by parallel hyperplanes on either side is maximized. The larger the margin, the lower the generalization error. The margin planes contain the examples of the classes that are nearest to each other. In the case of linearly not separable classes (when the classes overlap, see Figure \@ref(fig:svm-schematic-plot)), a soft-margin SVC may be employed. For examples on the correct side of the hyperplane, the loss is zero. For wrongly classified examples, the loss is proportional to the distance from the hyperplane. A global budget for loss is defined as a constraint and the hyperplane found subject to the constraint [@friedman2001elements].

```{python svm-schematic, include=F}
import matplotlib.pyplot as plt
from kdd98.config import Config
from scipy.stats import uniform
import numpy as np
palette = Config.get("qual_palette")

def f1(x):
    return -0.2*x + 0.4



def f2(x):
    return -0.2*x + 1.2

def hp(x):
    return -0.2*x + 0.8

np.random.seed(seed=42)
x2 = [0.6,1.2,2.1]
x1 = [0.2,1.5,1.8]
plt.plot([0, 2.1], [f1(x) for x in [0, 2.1]], color=palette[1])
plt.plot([0, 2.1], [f2(x) for x in [0, 2.1]], color=palette[1])
plt.plot(x1,[f1(x) for x in x1], 'o', color=palette[0])
plt.plot(x2,[f2(x) for x in x2], 'o', color=palette[5])
plt.plot([0.0,2.1], [hp(x) for x in [0, 2.1]], '-', color="black")
plt.plot(uniform(0.0,2.0).rvs(7), uniform(0.4,1.5).rvs(7), 'o', color=palette[5])
plt.plot(uniform(0.2,1.8).rvs(13), uniform(-0.5,1.3).rvs(13), 'o', color=palette[0])
plt.tick_params(axis='both', which='both', bottom=False, top=False, labelbottom=False, right=False, left=False, labelleft=False)
save_fig("methods/svm-schematic")
plt.close('all')
```

(ref:svm-schematic-legend) Schematic display of an SVM hyperplane (in black), separating two overlapping classes. The margins are shown around the hyperplane, with support vectors falling on the margins. Misclassifications (examples on the wrong side of the hyperplane) have a total budget for distance from the separating plane. The margins are determined by respecting the budget. Adapted from @friedman2001elements.

```{r svm-schematic-plot,  fig.cap="(ref:svm-schematic-legend)", echo=F}

knitr::include_graphics("figures/methods/svm-schematic.png")

```

Support Vector Machines (SVM), [@boser1992svc; @cortes1995support] are another approach to the overlapping class case. By mapping the original input space to a  high- or infinite-dimensional feature space using nonlinear transformations, a linear hyperplane separating the classes can be found. Calculation of the SVM involves a dot product between the $x_i$ [@friedman2001elements]:

\begin{equation}
f(x) = \sum_{i=1}^N \alpha_i*y_i\langle h(x), h(x_i)\rangle + \beta_0
(\#eq:svm)
\end{equation}

where $h(x)$ are the transformations mapping from the input space to the high-dimensional space.

Since the transformations can be prohibitively expensive, the so-called kernel trick [@aizerman1964theoretical] is used [@boser1992svc]. The trick lies in the fact that knowledge of the kernel function $K(x,x') = \langle h(x), h(x') \rangle$ suffices, without need to compute the dot products.

Several kernels are possible. In `scikit-learn`, linear, polynomial, radial basis function and sigmoid are available.

For regression, the Support Vector Regression Machine (SVR) was introduced by @drucker1997support.

`scikit-learn` implementations were used for both classification and regression.

**Hyperparameters**

* `C`, $(0,\inf)$: Penalty on the margin size
* `kernel`: One of linear, poly, rbf or sigmoid
* `degree` $\{1,2, \ldots \}$: Degree of polynomial kernel
* `class_weight`: Automatic balancing or a dictionary with weights per class

#### Bayesian Ridge Regression

Bayesian Ridge Regression (BR) can be seen as a Bayesian approach to a linear model with $l_2$ regularization. The `sklearn.linear_model.BayesianRidge` algorithm was used, which is implemented as described in @tipping2001sparse. The linear model is:

\begin{equation}
y_i = f(x_n, w) + \epsilon_i
(\#eq:bayesion-linear-model-sk)
\end{equation}

where $\epsilon_i$ are iid error terms $\sim \mathcal{N}(0, \sigma^2)$.

The probabilistic model is then:

\begin{equation}
p(y|x) =\mathcal{N}(y|f(x,w), \alpha)
(\#eq:bayesian-prob-model)
\end{equation}

In order to regularize the model in a ridge manner, a prior probability on the weights in $f(x,w)$ is defined as:

\begin{equation}
p(w|\lambda) = \prod_{i=0}^p \mathcal{N}(w_i|0, \lambda_i^-1)
(\#eq:prior-alpha)
\end{equation}

The algorithm was trained with default settings.