04-results-discussion.Rmd

# Results and Discussion

Below, the results from data preprocessing, model evaluation/-selection and predictions are shown.

## Preprocessing With Package kdd98

The self-written package `kdd98` ensures consistent data provisioning. It handles downloading and preprocessing of the data set for both the learning and validation set. All preprocessing steps are trained on the learning data set. The individual trained transformations are persisted on disk. After training, the transformations can be applied on the validation data. This process is transparent to the user. It is enough to instantiate the data provider for either learning or validation data set and request the data. Examples for usage can be found in the Jupyter notebooks.

The data sets can be obtained at intermediate steps from `kdd98.data_handler.KDD98DataProvider`:

* **raw**, as imported through `pandas.read_csv()`
* **preprocessed**, input errors removed, correct data types for all features, missing at random (MAR) imputations applied
* **numeric**, after feature engineering (encoded categories, date and zip code transformations)
* **imputed**, with missing values imputed
* **all-relevant**, filtered down to a set of relevant features

For some transformers, behavior can be controlled by specifying parameters (which has to be done in the code). The package's architecture makes it easy to implement additional transformation steps.

The source code, along with a short introduction, is available online^[`r make_github_link("kdd98")`].



## Imputation

The evaluation of several imputation strategies led to a straightforward approach: Categorical features had a *missing* level added during encoding (Section \@ref(methods-feature-engineering)). All other features were imputed by their median value to account for the skewed distributions. The notebook *4_Imputation.ipynb*^[`r make_github_link("notebooks", "4_Imputation.ipynb")`] contains details on the other approaches studied.

In concordance with the cup documentation's requirements (Appendix \@ref(data-set-documentation)), constant and *sparse* features were dropped from the data set before imputation. The approach used in the R-package `caret` was implemented in a scikit-learn transformer for this purpose. Work from a blog post [@goddard2017variance] was adapted. The beauty of the method used in `caret` is that it is data type agnostic. It works on the number of unique values per feature and the frequency ratio between the top and second values by count.

The features with most and least missing values are shown in Figure \@ref(fig:most-fewest-missing). It is not surprising to find the *MONTHS_TO_DONATION_\** features among those with most missing because few examples respond to the promotions with a donation.

Among the incomplete features with least missing values, there are several of the US census features. The *RFA_\** features give the status in reference to promotion $i$. All examples who did donate at some point before have an RFA status. Thus, the missingness indicates when new members were added because newly added members do not have an RFA status yet.

(ref:m-f-missing-legend) Features with most (left) and fewest (right) missing values.

```{r most-fewest-missing, fig.cap="(ref:m-f-missing-legend)", out.width="49%", fig.show="hold", echo=F}
include_graphics(c("figures/imputation/most-missing.png", "figures/imputation/fewest-missing.png"))
```

## Feature Selection

After preprocessing and feature engineering, `r ncol(py$numeric$data)` features were present in the data. Using boruta, `r ncol(py$all_relevant$data)` features were identified as important, resulting in a `r round(1 - ncol(py$all_relevant$data) / ncol(py$numeric$data),2)*100` % reduction of the number of features. For details, refer to notebook *5_Feature_Extraction.ipynb*^[`r make_github_link("notebooks", "5_Feature Extraction.ipynb")`].

Three groups of features were selected. By crosschecking with the dataset dictionary (Appendix \@ref(data-set-dictionary)), their meaning could be identified:

Features from the **giving history** broadly correspond to those used in classical RFM models mentioned in literature (Section \@ref(intro)). It is reassuring to find them among the all-relevant features:

* Donation amount for promotion 14
* Summary features: All-time donation amount, all-time number of donations, smallest, average and largest donation, donation amount of most recent donation
* 24 Features on frequency and amount of donations as per the date of past donations
* Time since first donation, Time since largest donation, time since last donation
* Number of donations in response to card promotions
* Number of months between first and second donation
* An indicator for *star* donors

The **promotion history** features can be interpreted as a measure of the importance of the examples to the organization. Those who receive many promotions are deemed valuable:

* Number of promotions received
* Number promotions in last 12 months before the current promotion
* Number of card promotions received
* Number of card promotions in last 12 months before the current promotion

Features from the **US census** data are concerned with the social status and wealth of the neighborhood of donors and by intuition make sense to be deemed relevant.

* Median and average home value
* Percentage of home values above some threshold (5 features)
* Percentage of renters paying more than 500 $
* DMA (designated market area, a geographical grouping)
* Median / average family / house income
* Per capita income
* Percentage of households with interest, rental or dividend income
* Percentage of adults with a Bachelor's degree
* Percentage of people born in state of residence

## Model Evaluation and Selection {#results-models}

Model evaluation and -selection may be studied in detail in notebook \  *6_Model_Evaluation_Selection.ipynb*^[`r make_github_link("notebooks", "6_Model_Evaluation_Selection.ipynb")`]. As will be explained below, all classifiers performed rather weak and were highly influenced by the imbalance in the data. The best results were achieved when using SMOTE resampling. Random over-/undersampling and specifying class weights had an inferior effect on model performance.

Among the classifiers evaluated, GLMnet  showed the best performance and was thus chosen.

For the regression models, RF outperformed the other models, although the differences were less pronounced compared to the classifier results.

### Classifiers

During grid search, models were trained individually for best F1 and recall. The models trained for high recall had only slightly worse precision than those trained for high F1, but at the same time better recall scores (see Jupyter notebook `6_Model_Evaluation_Selection.ipynb`). Therefore, the recall-trained models were considered for selection of the classifier.

Evaluation was based on the recall scores, confusion matrices, receiver operating characteristic (ROC) indicating model performance through an area under the curve score (ROC-AUC) and  precision-recall (PR) curves.

If only the recall score were considered, the decision would be obvious, as shown in Figure \@ref(fig:recall-scores). SVM has ~74 % recall, with the next-best scores at 54 % and 53 %. However, it is important to also consider the false positives as those cost money and decrease net profit.

(ref:recall-scores-legend) Comparison of recall scores for all classifiers evaluated.

```{r recall-scores, fig.cap="(ref:recall-scores-legend)", echo=F}

knitr::include_graphics("figures/learning/recall-scores.png")
```

The ROC-AUC curve (Figure \@ref(fig:roc-auc-pr), left panel) was constructed by evaluating the false positive rate (FPR) against the true positive rate (TPR) at various thresholds for the predicted class probabilities of examples in the training data. The closer the curve is to the top-left corner, the better a model performs (large TPR and at the same time a low FPR for a wide range of thresholds). All classifiers seemingly performed rather weak. In the case of imbalanced data, the majority class dominates this metric. The false positive rate is $FPR = \frac{FP}{FP+TN}$. This means that as the false positives (FP) decrease due to an increasing threshold, FPR does not change a lot.

The PR curve with precision $P = \frac{TP}{TP+FP}$ plotted against recall $R = \frac{TP}{{TP+FN}}$ at different threshold values is sensitive to false positives and, since TN is not  involved, better suited for the imbalanced data at hand. To construct the curve, recall is plotted against precision for various threshold values of the predicted class probabilities. For good models, the curve is close to the top-right corner. The right panel of Figure \@ref(fig:roc-auc-pr) shows the models in direct comparison. All of them suffer from low precision except for the highest threshold values. Again, this is caused by the high imbalance in the data.

(ref:roc-auc-pr-legend) Comparison of ROC-AUC for the evaluated classifiers (left) and PR curves (right) for the classifiers.

```{r roc-auc-pr, fig.cap="(ref:roc-auc-pr-legend)", fig.show="hold", out.width="49%", echo=F}

knitr::include_graphics(c("figures/learning/roc_auc_compared_refit_recall.png", "figures/learning/prec_rec_compared_refit_recall.png"))
```

More insight can be gained from the confusion matrices, shown in Figure \@ref(fig:conf-matrices). The best classifier should have a high recall, predicting many donors correctly, and at the same time a low False Positive Rate (FPR). 

From the confusion matrices, it becomes obvious that for SVM, the trade off for high recall is also a high FPR. The false positives can be directly translated to cost: the 12'985 false positives in this case would amount to 8'830 $ at a unit cost of 0.68 $, while the 715 true positives generate an expected profit of only 8'808 $ (with a mean net profit of 12.32 $).

Evidently, GLMnet an NNet have a relatively good balance of recall and FPR. GBM performs well for FPR, which means less money lost due to unit costs, but has a very low recall.

GLMnet and NNet were also combined into a voting classifier. This classifier creates an ensemble that predicts through a majority vote, therefore compensating for the individual classifier's weaknesses. It exhibits a slightly lower recall for a slight decrease in FPR. Since *recall* was seen as being more important, it was not investigated further.

(ref:conf-mat-legend) Confusion matrices for the 6 classifiers studied.

```{r conf-matrices, fig.cap="(ref:conf-mat-legend)", fig.show="hold", out.width="30%", echo=F}
knitr::include_graphics(c(
	"figures/learning/confusion_matrix_model_RF_refit_recall.png",
	"figures/learning/confusion_matrix_model_GBM_refit_recall.png",
	"figures/learning/confusion_matrix_model_GLMnet_refit_recall.png",
	"figures/learning/confusion_matrix_model_NNet_refit_recall.png",
	"figures/learning/confusion_matrix_model_SVM_refit_recall.png",
	"figures/learning/confusion_matrix_voting_classifier.png"
	))
```

The non-null coefficients of the GLMnet model are shown in \@ref(fig:glmnet-coefficients). 22 coefficients are non-null. The absolute values of the coefficients indicate importance of the respective feature. The first five features are from the promotion- and giving history and are concerned with the amount and frequency of donations.

```{r glmnet-coefficients, fig.cap="Coefficient values for GLMnet.", echo=F, out.width="100%"}

knitr::include_graphics("figures/learning/glmnet_coefficients.png")
```

Although RF performed poorly, it can be used to have a second source of information on the important features. The measure implemented in scikit-learn is the *Gini importance* as described in @breiman1984classification. It represents the total decrease in impurity (see Section \@ref(methods-rf)) due to nodes split on feature $f$, averaged over all trees in the forest. The most important features again are all from the giving history, followed by the promotion history. US census features are not important (see Figure \@ref(fig:importances)).

```{r importances, fig.cap="Feature importances determined with the RF classifier. Impurity is measured by Gini importance / mean decrease impurity. Error bars give bootstrap error on 50 repetitions.", echo=F, out.width="100%"}

knitr::include_graphics("figures/learning/feature-importance-rf-classification.png")
```

### Regressors

```{python regressor-stuff, include=F, cache=F}
import pickle
with open("models/target_d_transformer.pkl", "rb") as f:
    target_d_transformer = pickle.load(f)
box_cox_lambda = target_d_transformer.lambdas_[0]
```
As explained in Section \@ref(methods-prediction), the regressors were learned on a subset of the training data comprised of all donors: $\{\{x_i, y_i\}|y_{b,i} = 1\}$. Before learning, *TARGET_D* was Box-Cox transformed for normalization with  $\lambda=`r py$box_cox_lambda`$. The goal of the transformation was to improve regression models' performance. The transformed data somewhat resembles a normal distribution, although there are several modes to be made out (see Figure \@ref(fig:reg-targ-transform)).

```{r reg-targ-transform, fig.cap="Target before transformation (left) and after a Box-Cox transformation (right).", echo=F}

knitr::include_graphics(c("figures/learning/target_d-distributions-before-after-transformation.png"))
```

The resulting distribution of predicted donation amounts on the training data is shown in Figure \@ref(fig:reg-distrib). Except for SVR, the models produce very similar results. Again, the multi-modal distribution is found (refer to Figure \@ref(fig:reg-targ-transform)). RF and SVR are relatively symmetric, while BR and ElasticNet produce right-skewed distributions that predict very large donation amounts for some examples.

(ref:reg-distrib-legend) Distribution of (Box-Cox transformed) donation amounts for the four regressors evaluated.

 ```{r reg-distrib, fig.cap="(ref:reg-distrib-legend)", echo=F}

knitr::include_graphics("figures/learning/regressor-predictions-comparison.png")
```
The regressor for final prediction was selected by $R^2$ score on a test set (20% of the learning data was used). RF was the best performing model with $R^2 = 0.72$ (see Figure \@ref(fig:reg-eval)).

(ref:reg-eval-legend) Evaluation metric $R^2$ for all regression models evaluated. The domain for $R^2$ is $(-\inf, 1]$.)

```{r reg-eval, fig.cap="(ref:reg-eval-legend)", echo=F}

knitr::include_graphics("figures/learning/regressor-score-comparison.png")
```

Again, RF enables to interpret the importance of features. The results (Figure \@ref(fig:reg-importance)) show that mainly the amount of the last donation, and to a lesser amount the average donation amount and the amounts of the last three donations were important for predicting donation amounts. This would indicate that donors tend to always donate the same sums.

(ref:reg-importance-legend) Feature importances for the RF regressor.

```{r reg-importance, fig.cap="(ref:reg-importance-legend)", echo=F}

knitr::include_graphics("figures/learning/feature-importance-rf-regression.png")
```

## Prediction

The intermediate steps for arriving at the final prediction will be examined in this section. For details on the steps, see *7_Predictions.ipynb*^[`r make_github_link("notebooks", "7_Predictions.ipynb")`]. Up until the final prediction, the results reported are on the learning data set. The test data was only used for the final prediction for comparison with the results of the participants of the cup.

### Prediction of Donation Probability

Using the GLMnet regressor on the learning data set, $\hat{y}_b$ is distributed close to normal with a slight right skew (see Figure \@ref(fig:y-b-predict)). 33'559 examples are selected, which amounts to `r 33559/95412*100` % of the training data set.

(ref:y-b-predict-legend) Distribution of predicted donation probabilities $\hat{y}_b$ on the learning data set.

```{r y-b-predict, fig.cap="(ref:y-b-predict-legend)", echo=F}
knitr::include_graphics("figures/predictions/y_b_predicted.png")
```

### Conditional Prediction of the Donation Amount

Compared to the true distribution of $y_d$ (see Figure \@ref(fig:reg-targ-transform)), a similar distribution for $\hat{y}_d$ results from predictions with the RF regressor on the complete learning data, shown in Figure \@ref(fig:y-d-predict). The distribution is again multi-modal, with approximately the same median value. The predicted donation amounts are however strictly positive with a minimum of 3.2 \$. The reason is that the predictions are biased due to the non-random sample of donors used to learn the regressor. Unfortunately, the high-dollar donations are missing. The highest prediction is 82.1 \$, the median is 30.7 \$, which is almost double the median of the complete learning data set.

(ref:y-d-predict-legend) Conditionally predicted donation amounts, Box-Cox transformed (left) and on the original scale (right).

```{r y-d-predict, fig.cap="(ref:y-d-predict-legend)", echo=F}

knitr::include_graphics("figures/predictions/hat_y_d-distributions-before-after-transformation.png")
```

### Profit Optimization

The correction factor $\alpha$, used to account for bias introduced by learning the regression models on a non-random sample, was found as described in section \@ref(methods-prediction). First, estimated net profit $\hat{\Pi}$ was calculated using equation \@ref(eq:pi-alpha) for a grid of $\alpha$ values in $[0,1]$. Then, a cubic spline was fitted to the data and $\alpha$ optimized subsequently.

As can be observed in Figure \@ref(fig:alpha-grid), the region of $\alpha$ for high profit is narrow. This is not surprising given the distribution of $\hat{y}_d$ (Figure \@ref(fig:y-d-predict)), which is narrowly concentrated around ~15 \$. Furthermore, the curve is constant over much of the domain, meaning that all examples were selected from $\alpha \approx 0.15$. The cubic spline fits the data very well and finding the maximum at $\alpha^*=0.02$. For reference, a polynomial of degree 12 is also shown in \@ref(fig:alpha-grid), highlighting the difference in the fit for the two approaches.

(ref:alpha-grid-legend) Expected profit for a range of $\alpha$ values in $[0,1]$ with overlayed cubic spline and polynomial function of order 12.

```{r alpha-grid, fig.cap="(ref:alpha-grid-legend)", echo=F}
knitr::include_graphics("figures/predictions/comparison-alpha-profit-models.png")
```

### Final Prediction

For the final prediction, an estimator combining the two-step prediction process was implemented in package `kdd98`. The estimator `Kdd98ProfitEstimator` was initialized with the GLMnet classifier and RF regressor determined as the best during model evaluation and selection.

Then, the estimator was fitted on the complete learning data set, thereby fitting both the classifier and regressor and optimizing $\alpha*$, which was computed as $\alpha^* = 0.023$.

The fitted estimator could then be used for prediction on the test data set. The results are shown in Table \@ref(tab:prediction-results) together with the winners of the cup. The model learned in this thesis selects less donors compared to the top-ranked participants and has a higher mean donation amount. Nevertheless, net profit is lower, resulting in the $4^{th}$ rank.

The winner, Urban Science Applications with their proprietary software Gainsmarts, chose a two-step approach as well. In the first step, a logistic regression was used. The second step consisted of linear regression.
Their software automated feature-engineering (by trying several different transformations for each feature) and feature selection through an expert system^[see [https://www.kdnuggets.com/meetings-past/kdd98/gain-kddcup98-release.html](https://www.kdnuggets.com/meetings-past/kdd98/gain-kddcup98-release.html)]. They selected a distinct subset of features for each step.

No information is available on SAS's approach. Decisionhouse's Quadstone on $3^{rd}$ rank also invested heavily in feature engineering. The proprietary software was specially designed for customer behavior modeling. They used decision trees for conditioned profit estimation using 6 features and an additive scorecard model for selecting examples based on 10 features in combination. For the estimation of donation amount, they also found the amount of the last donation and average donation amounts as important.

Surely, if domain knowledge had influenced feature engineering and feature selection for the model developed here, performance would have been better. Given the radically data-driven approach, the result can be considered good.

(ref:prediction-results-legend) Prediction results for the test data set (in color) and the results of the cup-winners. $N^*$ denotes number of examples selected. The theoretical maximum was calculated with: \$ ($\sum_{i=1}^n \mathbb{1}_{\{\text{TARGET}_{D,i} > 0.0\}}*(\text{TARGET}_{D,i} - u)$ with $u$ the unit cost of 0.68 $ per mailing).


```{r prediction-results, echo=F}

predictions <- t(data.frame(GainSmarts = c(56330, -0.68, 0.26, 5.57, 499.32, 14712, 14712/72776*100),
                            SAS = c(55838, -0.68, 0.26, 5.64, 499.32, 14662, 14662/72776*100),
                            Quadstone = c(57836, -0.68, 0.24, 5.66, 499.32, 13954, 13954/72776*100),
                            Own = c(40984, -0.68, 0.338, 6.28, 499.32, 13877, 13877/72776*100),
                            CARRL = c(55650, -0.68, 0.25, 5.61, 499.32, 13825, 13825/72776*100)))
colnames(predictions) <- c('N*','Min','Mean','Std', 'Max','Net Profit' ,'Percent of Maximum')

knitr::kable(predictions,
             booktabs = T, digits=2,
             caption = "(ref:prediction-results-legend)") %>%
    kableExtra::kable_styling(latex_options=c("hold_position", position="Center")) %>%
    kableExtra::add_header_above(c(" "=1," "=1,"Amount, $"=5, " "=1)) %>%
    kableExtra::row_spec(4, bold=TRUE, color="#39a8d1")
```

However, these predictions were made using the true *TARGET_D*. When treating the data as unseen, predicted net profit was unrealistically high, at `r 861562` \$, which translates to an error of `r 6208` %. In this case, net profit was calculated by substituting $(y_d-u)$ with $(\hat{y}_d-u)$ in Equation \@ref(eq:pi-alpha). However bad this result is, would the model developed be used to select examples for a promotion, it would still lead to relatively good net profit.