---
title: "Estimating predictive models"
author: "David Kepplinger"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimating predictive models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", error = TRUE, warning = TRUE)
load("computing_adapense_fits.RData")
```

In this guide we will estimate predictive models by means of robust adaptive PENSE estimates for high-dimensional linear regression.
These estimates can tolerate up to 50\% of contamination, i.e., the adaptive PENSE estimates are reliable even if up to half the observations in the data set contain anomalous values.
Compute adaptive PENSE estimates is implemented in the function `adapense_cv()` in the pense package.

While the following guide computes adaptive PENSE estimates, everything also applies to other estimates implemented in the pense package: non-adaptive PENSE estimates and regularized M-estimates.
Non-adaptive PENSE estimates (computed by `pense_cv()`) are typically better at identifying all relevant predictors than adaptive PENSE.
However, this comes at the price of often including a large number of irrelevant predictors as well.
Regularized M-estimates (computed by `pensem_cv()`) can be more accurate than either PENSE or adaptive PENSE estimates, but may be unreliable in presence of highly detrimental contamination.

# Computing adaptive PENSE estimates

First we need to load the pense package:
```{r}
library(pense)
```

Computing robust, regularized estimates for high-dimensional linear regression models can take a long time.
To save time, many steps can be done in parallel.
If your computer has more than 1 CPU core, you can harness more computing power by creating a cluster of R processes:

```{r}
library(parallel)
# If you don't know how many CPU cores are available, first run `detectCores(logical = FALSE)`
cluster <- makeCluster(3)
```

This guide uses the following simulated data with 50 observations and 40 available predictors.
The error distribution is a heavy-tailed *t*-distribution and only the first 3 predictors are truly relevant for predicting the response *y*:
```{r}
set.seed(1234)
x <- matrix(rweibull(50 * 40, 2, 3), ncol = 40)
y <- 1 * x[, 1] - 1.1 * x[, 2] + 1.2 * x[, 3] + rt(nrow(x), df = 2)
```

To make the scenario more realistic, let's add some contamination to the response value of the first 3 observations and to some predictors:
```{r}
y[1:3] <- 5 * apply(x[1:3, ], 1, max)
x[3:6, 4:6] <- 1.5 * max(x) + abs(rcauchy(4 * 3))
```

## Step 1: Computing the estimates

The first step is to compute adaptive PENSE estimates for a fixed value for hyper-parameter `alpha` and many different penalization levels (hyper-parameter `lambda`).
The `adapense_cv()` function automatically determines a grid of penalization levels, with parameter `nlambda=` controlling the number of different penalization levels to be used (default 50).
We are going to choose the penalization level which leads to a good balance between prediction accuracy and model size.
The pense package can automatically evaluate prediction accuracy of the adaptive PENSE estimates via cross-validation.

In its simplest form, computing adaptive PENSE estimates and estimating their prediction accuracy is done with the code below.
Here, prediction accuracy is estimated via 5-fold cross-validation, replicated 3 times:
```{r fit_075, eval=FALSE}
set.seed(1234)
fit_075 <- adapense_cv(x, y, alpha = 0.75, cv_k = 5, cv_repl = 3, cl = cluster)
```
```{r, eval=FALSE, include=FALSE}
fit_075_3repl <- fit_075
```

You should always set a seed prior to computing adaptive PENSE or PENSE estimates to ensure reproducibility of the internal cross-validation.

By default, adaptive PENSE estimates are computed with a breakdown point of ~25\%.
This means, the estimates are reliable if up to 25\% of observations contain arbitrary contamination.
If you suspect that a larger proportion of observations may be affected by contamination, you can increase the breakdown point of the estimates with argument `bdp=` to up to 50\%.
Note, however, that a higher breakdown point also leads to less accurate estimates.

## Step 2: Assessing prediction performance

The `plot()` function for the object `fit_075` shows the estimated prediction accuracy for all fitted models.

```{r, include=FALSE}
fit_075 <- fit_075_3repl
```
```{r, fig.width=7, fig.height=5, fig.cap="Estimated prediction accuracy using 3 replications of 5-fold CV."}
plot(fit_075)
```

The plot shows the estimated scale of the prediction error for all 50 models.
The penalization level leading to the best prediction performance is highlighted by a dark blue dot.
If more than one CV replication was performed, the plot also shows a light blue dot, marking the most parsimonious model with prediction performance "indistinguishable" from the best model.
The plot uses the "one-standard-error" rule using the minimum average scale of the prediction error plus 1 standard error of this estimated scale.
You can adjust this rule to the "*m*-standard-error" with `plot(fit_075, lambda = "m-se")`, where `m` is any positive number (e.g., `lambda = "1.5-se"`).

In our case, the estimated prediction performance is a bit unstable and has large variability.
With such unstable estimates of prediction performance it is difficult to reliably select the penalization level.
This is not unusual for robust estimators and can be improved by increasing the number of CV replications.

The more CV replications, the more accurate the estimates of prediction accuracy, but the longer the computing time.
If we repeat step #1, but with 10 CV replications instead of 3, we get a more stable evaluation of prediction performance:

```{r, eval=FALSE}
set.seed(1234)
fit_075 <- adapense_cv(x, y, alpha = 0.75, cv_k = 5, cv_repl = 10, cl = cluster)
plot(fit_075)
```
```{r, eval=FALSE, include=FALSE}
fit_075_10repl <- fit_075
```
```{r, fig.width=7, fig.height=5, echo=FALSE, fig.cap="Estimated prediction accuracy using 10 replications of 5-fold CV."}
fit_075 <- fit_075_10repl
plot(fit_075)
```

With 10 replications, the error bars are still large, a testament to the contamination in the sample, but the general trend of the prediction performance is much clearer.
Of course, you could increase the number of CV replications further to get an even smoother plot.

## Step 3: Extracting coefficients

Once we are happy with the stability of the estimated prediction performance, we can extract summary information from the predictive model with
```{r}
summary(fit_075)
```

This model corresponds to the model with smallest scale of the prediction error (the blue dot in the plot above).
There are a total of `r sum(abs(coef(fit_075)) > 0) - 1` predictors in the model.
If you think a sparser model may be more appropriate for your application, you can also apply the *q*-standard-error rule as in the plots.
The default, the one-standard-error rule leads to the following predictive model:
```{r}
summary(fit_075, lambda = "2-se")
```

In this fit, only `r sum(abs(coef(fit_075, lambda = "2-se")) > 0) - 1` out of the 40 predictors are relevant, including the 3 truly relevant predictors.
But maybe different values for hyper-parameters `alpha` and `exponent` lead to even better prediction?

## Step 4: Exploring different hyper-parameters

The choice for hyper-parameters `alpha` and `exponent` (which was kept at its default value of 1) are rather arbitrary.
The effects of these two hyper-parameters on the estimates are in general less pronounced than of the penalization level.
But you may still want to explore different values for `alpha` and `exponent`.

While `alpha=0.75` is a good value for many applications, `alpha=1` may also be of interest, particularly in applications where correlation between predictors is not an issue.
In applications with high correlation between predictors, lower values of `alpha` (e.g., `alpha=0.5`) may lead to more stability in variable selection.

The hyper-parameter `exponent` generally has an effect on the sparsity of the models.
With higher values for `exponent`, typically only predictors with the largest (in absolute magnitude) standardized coefficients will be non-zero.
While this helps to screen out many or most of the truly irrelevant, it also risks missing some of the truly relevant predictors.

Let us compute adaptive PENSE estimates for different values of hyper-parameters `alpha` and `exponent`.
In the code below, this is done for two values: `alpha=0.75` and `alpha=1` as well as `exponent=1`, `exponent=2`, and `exponent=3`.

```{r fit_all, eval=FALSE}
set.seed(1234)
fit_exp_1 <- adapense_cv(x, y, alpha = c(0.75, 1), exponent = 1, cv_k = 5, cv_repl = 10, cl = cluster)

set.seed(1234)
fit_exp_2 <- adapense_cv(x, y, alpha = c(0.75, 1), exponent = 2, cv_k = 5, cv_repl = 10, cl = cluster)

set.seed(1234)
fit_exp_3 <- adapense_cv(x, y, alpha = c(0.75, 1), exponent = 3, cv_k = 5, cv_repl = 10, cl = cluster)
```

Note that we set the seed before each call to `adapense_cv()`.
This is again to ensure reproducibility of the CV, but also to make the estimated prediction performance more comparable across different values of the `exponent` hyper-parameter.

After checking that the cross-validated prediction performance of the fitted models is smooth enough to reliably select the penalization level, we can compare all the estimates.
For this, the package includes the function `prediction_performance()`, which extracts and prints the prediction performance of all given objects:

```{r}
prediction_performance(fit_exp_1, fit_exp_2, fit_exp_3)
```

Here we see the combination of hyper-parameters `alpha=0.75` and `exponent=3` (in object `fit_exp_3`) leads to the best prediction accuracy, but all models are based on 12 or more predictors.
Instead, we can also compare sparser models with almost the same prediction performance:

```{r}
prediction_performance(fit_exp_1, fit_exp_2, fit_exp_3, lambda = '1-se')
```

These models are still quite large.
If we are interested in even sparser models, we can increase the tolerance level to, for instance, 3 standard errors:

```{r}
prediction_performance(fit_exp_1, fit_exp_2, fit_exp_3, lambda = '3-se')
```

This relaxation still leads to the combination of hyper-parameters `alpha=0.75` and `exponent=3`, but now with only 7 relevant predictors; including all truly relevant predictors.
We can see that the estimated coefficients and estimated relevant predictors are close to the truth:

```{r}
summary(fit_exp_3, alpha = 0.75, lambda = '3-se')
```

## Using different measures of prediction performance

By default, `adapense_cv()` uses the τ-scale of the prediction errors to assess prediction accuracy.
This can be changed by specifying a different metric via `adapense_cv(cv_metric=)`.
The package supports also the median absolute prediction error (`cv_metric = "mape"`) or the classical root mean squared prediction error (`cv_metric = "rmspe"`).
You should, however, not use the RMSPE to evaluate prediction performance in the potential presence of contamination.
Robust methods are not designed to predict contaminated observations well and the RMSPE may be artificially inflated by poor prediction of a few contaminated response values.
You can also specify your own function which takes as input the vector of prediction errors and returns a single number, measuring the prediction performance.
For example, to use the *mean* absolute prediction error, you would write

```{r, eval=FALSE}
mae <- function (prediction_errors) {
  mean(abs(prediction_errors))
}

set.seed(1234)
fit_075_mae <- adapense_cv(x, y, alpha = 0.75, cv_k = 5, cv_repl = 5, cl = cluster, cv_metric = mae)
```

A matrix with estimates of the prediction performance are accessible as slot `$cv_replications` in the object returned by `adapense_cv()`.
The rows correspond to the different penalization levels, and columns correspond to the individual CV replications.

```{r, eval=FALSE, include=FALSE}
save(fit_075_10repl, fit_075_3repl, fit_075_mae, fit_exp_1, fit_exp_2, fit_exp_3,
     file = 'computing_adapense_fits.RData')
```