beyond_regression_and_classification.qmd

---
aliases:
  - "/beyond_regression_and_classification.html"
  - "/survival.html"
---

# Beyond Regression and Classification {#sec-special}

{{< include ../../common/_setup.qmd >}}

`r chapter = "Beyond Regression and Classification"`
`r authors(chapter)`

So far, this book has only considered two tasks.
In @sec-basics we introduced deterministic regression as well as deterministic and probabilistic single-label classification (@tbl-alltasks).
But our infrastructure also works well for many other tasks, some of which are available in extension packages (@fig-mlr3verse) and some are available by creating pipelines with `r mlr3pipelines`.
In this chapter, we will take you through just a subset of these new tasks, focusing on the ones that have a stable API.
As we work through this chapter we will refer to the 'building blocks' of `r mlr3`, this refers to the base classes that must be extended to create new tasks, these are `r ref("Prediction")`, `r ref("Learner")`, `r ref("Measure")`, and `r ref("Task")`.
@tbl-alltasks summarizes available extension tasks, including the package(s) they are implemented in and a brief description of the task.

| Task | Package | Description |
| ---- | -- | ---- |
| Deterministic regression | `r mlr3`  | Point prediction of a continuous variable. |
| Deterministic single-label classification | `r mlr3`  | Prediction of a single class for each observation. |
| Probabilistic single-label classification | `r mlr3`  | Prediction of the probability of an observation falling into one or more mutually exclusive categories. |
| Cost-sensitive classification | `r mlr3` and `r mlr3pipelines`  | Classification predictions with unequal costs associated with misclassifications. |
| Survival analysis | `r mlr3proba`  | Time-to-event predictions with possible 'censoring'. |
| Density estimation | `r mlr3proba`  | Unsupervised estimation of probability density functions. |
| Spatiotemporal analysis | `r mlr3spatiotempcv` and `r mlr3spatial`  | Supervised prediction of data with spatial (e.g., coordinates) and/or temporal outcomes. |
| Cluster analysis | `r mlr3cluster`  | Unsupervised estimation of homogeneous clusters of data points. |

: Table of extension tasks that can be used with `mlr3` infrastructure. As we have a growing community of contributors, this list is far from exhaustive and many 'experimental' task implementations exist; this list just represents the tasks that have a functioning interface. {#tbl-alltasks}

## Cost-Sensitive Classification {#sec-cost-sens}

We begin by discussing a task that does not require any additional packages or infrastructure, only the tools we have already learned about from earlier chapters.
In 'regular' classification, the aim is to optimize a metric (often the misclassification rate) while assuming all misclassification errors are deemed equally severe.
A more general approach is `r index("cost-sensitive classification", "cost-sensitive", parent = "classification", aside = TRUE)`, in which costs caused by different kinds of errors may not be equal.
The objective of cost-sensitive classification is to minimize the expected costs.
We will use `tsk("german_credit")` as a running example.

Imagine you are trying to calculate if giving someone a loan of \$5K will result in a profit after one year, assuming they are expected to pay back \$6K.
To make this calculation, you will need to predict if the person will have good credit.
This is a deterministic classification problem where we are predicting whether someone will be in class 'Good' or 'Bad'.
Now let us consider some potential costs associated with each prediction and the eventual truth.
As cost-sensitive classification is a minimization problem, we assume lower costs correspond to higher profits/positive outcomes, hence we write profits as negative values and losses as positive values:

```{r special-002}
costs = matrix(c(-1, 0, 5, 0), nrow = 2, dimnames =
  list("Predicted Credit" = c("good", "bad"),
    Truth = c("good", "bad")))
costs
```

In this example, if the model predicts that the individual has bad credit (bottom row) then there is no profit or loss, the loan is not provided.
If the model predicts that the individual has good credit and indeed the customer repays the loan with interest (top left), then you will make a \$1K profit.
On the other hand, if they default (top right), you will lose \$5K.

### Cost-Sensitive Measure

We will now see how to implement a more nuanced approach to classification errors with `msr("classif.costs")`.
This measure takes one argument, which is a matrix with row and column names corresponding to the class labels in the task of interest.
Let us put our insurance example into practice, notice that we have already named the cost matrix as required for the measure:

```{r special-003}
library(mlr3verse)

tsk_german = tsk("german_credit")

msr_costs = msr("classif.costs", costs = costs)
msr_costs

learners = lrns(c("classif.log_reg", "classif.featureless",
  "classif.ranger"))
bmr = benchmark(benchmark_grid(tsk_german, learners,
  rsmp("cv", folds = 3)))
bmr$aggregate(msr_costs)[, c(4, 7)]
```

In this experiment, we find that the logistic regression learner happens to perform best as it minimizes the expected costs (and maximizes expected profits) and the featureless learner performs the worst.
All losses result in positive costs, which means each model results in us losing money.
To improve our models, we will now turn to thresholding.

### Thresholding

As we have discussed in @sec-basics, `r index("thresholding")` is a method to fine-tune the probability at which an observation will be predicted as one class label or another.
Currently in our running example, the models above will predict a customer has good credit (in the class 'Good') if the probability of good credit is greater than 0.5.
Here, this might not be a sensible approach as we would likely act more conservatively and reject more credit applications with a higher threshold due to the non-uniform costs.
This is highlighted in the `"threshold"` `autoplot` (@fig-costsens-threshold), which plots `msr("classif.costs")` over all possible thresholds.

```{r special-004}
#| output: false
#| cache: false
prediction = lrn("classif.log_reg",
  predict_type = "prob")$train(tsk_german)$predict(tsk_german)
autoplot(prediction, type = "threshold", measure = msr_costs)
```

```{r}
#| echo: false
#| warning: false
#| message: false
#| label: fig-costsens-threshold
#| fig-cap: Changing values of cost-sensitive measure as the prediction threshold is changed.
#| fig-alt: "Line graph with x-axis labeled 'Probability Threshold' ranging between 0-1, and y-axis labeled 'classif.costs' ranging between -0.3 and 0.7'. Line starts at (0,0.7) decreases linearly to around (0.8,-0.3) then increases linearly to (1, 0)."
plt = ggplot2::last_plot()
plt$layers[[1]]$aes_params$colour = "grey30"
print(plt)
```

As expected, the optimal threshold is greater than 0.5 which means the optimal model should predict 'bad' credit more often than not.

The optimal threshold can be automated by making use of `r mlr3tuning` (@sec-optimization) and `r mlr3pipelines` (@sec-pipelines) to tune `po("tunethreshold")`.
Continuing the same example:

```{r special-005}
po_cv = po("learner_cv", lrn("classif.log_reg", predict_type = "prob"))
graph =  po_cv %>>% po("tunethreshold", measure = msr_costs)

learners = list(as_learner(graph), lrn("classif.log_reg"))
bmr = benchmark(benchmark_grid(tsk_german, learners,
  rsmp("cv", folds = 3)))
bmr$aggregate(msr_costs)[, c(4, 7)]
```

By using `po("learner_cv")` for internal resampling and `po("tunethreshold")` to find the optimal threshold we have improved our model performance considerably and can now even expect a profit.

## Survival Analysis {#sec-survival}

`r index("Survival analysis")` is a field of statistics concerned with trying to predict/estimate the time until an event takes place.
This predictive problem is unique as survival models are trained and tested on data that may include 'censoring', which occurs when the event of interest does *not* take place.
Survival analysis can be hard to explain in the abstract, so as a working example consider a marathon runner in a race.
Here the 'survival problem' is trying to predict the time when the marathon runner finishes the race.
However, if the event of interest does not take place (e.g., the marathon runner gives up and does not finish the race), they are said to be censored.
Instead of throwing away information about censored events, survival analysis datasets include a status variable that provides information about the 'status' of an observation.
So in our example, we might write the runner's outcome as $(4, 1)$ if they finish the race at four hours, otherwise, if they give up at two hours we would write $(2, 0)$.

The key to modeling in survival analysis is that we assume there exists a hypothetical time the marathon runner would have finished if they had not been censored, it is then the job of a survival learner to estimate what the true survival time would have been for a similar runner, assuming they are *not* censored (see @fig-censoring).
Mathematically, this is represented by the hypothetical event time, $Y$, the hypothetical censoring time, $C$, the observed outcome time, $T = \min(Y, C)$, the event indicator $\Delta := (T = Y)$, and as usual some features, $X$.
Learners are trained on $(T, \Delta)$ but, critically, make predictions of $Y$ from previously unseen features.
This means that unlike classification and regression, learners are trained on two variables, $(T, \Delta)$, which, in R, is often captured in a `r ref("survival::Surv")` object.
Relating to our example above, the runner's outcome would then be $(T = 4, \Delta = 1)$ or $(T = 2, \Delta = 0)$.
Another example is in the code below, where we randomly generate six survival times and six event indicators, an outcome with a `+` indicates the outcome is censored, otherwise, the event of interest occurred.

```{r special-009}
library(survival)
Surv(runif(6), rbinom(6, 1, 0.5))
```

Readers familiar with survival analysis will recognize that the description above applies specifically to 'right censoring'.
Currently, this is the only form of censoring available in the `mlr3` universe, hence restricting our discussion to that setting.
For a good introduction to survival analysis see @Collett2014 or for machine learning in survival analysis specifically see @MLSA.

For the remainder of this section, we will look at how `r mlr3proba` [@mlr3proba] extends the building blocks of `mlr3` for survival analysis. We will begin by looking at objects used to construct machine learning tasks for survival analysis, then we will turn to the learners we have implemented to solve these tasks, before looking at measures for evaluating survival analysis predictions, and then finally we will consider how to transform prediction types.

```{r special-010, echo=FALSE}
#| fig-alt: Figure shows give horizontal lines at 1,2,3,4,5 on the y-axis and a vertical line at 8 on the x-axis. Top line (subject 5) has a circle at x=8 and a diamond at x=9, second line (subject 4) has a circle at x=1 and a diamond at x=9, subject 3 has a circle at x=4 and a diamond at x=6, subject 2 has a diamond at x=8, and subject 1 has a diamond at x=7.
#| fig-cap: Plot illustrating different censoring types. Dead and censored subjects (y-axis) over time (x-axis). Black diamonds indicate true death times and white circles indicate censoring times. Vertical line is the study end time. Subjects 1 and 2 die in the study time. Subject 3 is censored in the study and (unknown) dies within the study time. Subject 4 is censored in the study and (unknown) dies after the study. Subject 5 dies after the end of the study. Figure and caption from @Sonabend2021b.
#| label: fig-censoring
library(ggplot2)
ggplot(data.frame(x = c(1, 2, 5, 7, 10), y = 1:5), aes(x = x, y = y)) +
  geom_blank() +
  geom_vline(xintercept = 8) +
  geom_segment(y = 1, yend = 1, x = 0, xend = 7) +
  geom_segment(y = 2, yend = 2, x = 0, xend = 8) +
  geom_segment(y = 3, yend = 3, x = 0, xend = 4) +
  geom_segment(y = 3, yend = 3, x = 4, xend = 6, linetype = 2) +
  geom_segment(y = 4, yend = 4, x = 0, xend = 1) +
  geom_segment(y = 4, yend = 4, x = 1, xend = 9, linetype = 2) +
  geom_segment(y = 5, yend = 5, x = 0, xend = 8) +
  geom_segment(y = 5, yend = 5, x = 8, xend = 9, linetype = 2) +
  geom_point(x = 4, y = 3, shape = 21, color = "black", fill = "white", size = 4) +
  geom_point(x = 6, y = 3, shape = 18, size = 4) +
  geom_point(x = 7, y = 1, shape = 18, size = 4) +
  geom_point(x = 8, y = 2, shape = 18, size = 4) +
  geom_point(x = 9, y = 5, shape = 18, size = 4) +
  geom_point(x = 8, y = 5, shape = 21, size = 4, fill = "white") +
  geom_point(x = 1, y = 4, shape = 21, size = 4, fill = "white") +
  geom_point(x = 9, y = 4, shape = 18, size = 4) +
  labs(y = "Subject", x = "Time") +
  scale_x_continuous(labels = as.character(1:10), breaks = 1:10) +
  theme_minimal() + theme(panel.grid.minor = element_blank())
```

### TaskSurv

As we saw in the introduction to this section, survival algorithms require two targets for training, this means the new `r ref("TaskSurv")` object expects two targets.
The simplest way to create a survival task is to use `r ref("as_task_surv()")`, as in the following code chunk.
Note this has more arguments than `r ref("as_task_regr()")` to reflect multiple target and censoring types, `time` and `event` arguments expect strings representing column names where the 'time' and 'event' variables are stored, `type` refers to the censoring type (currently only right censoring supported so this is the default).
`as_task_surv()` coerces the target columns into a `r ref("survival::Surv")` object.
In this section we will use the `rats` dataset as a running example, this dataset looks at predicting if a drug treatment was successful in preventing 150 rats from developing tumors.
The dataset, by its own admission, is not perfect and should generally be treated as 'dummy' data, which is good for examples but not real-world analysis.

```{r special-011, warning=FALSE, message=FALSE}
library(mlr3verse)
library(mlr3proba)
library(survival)

tsk_rats = as_task_surv(survival::rats, time = "time",
  event = "status", type = "right", id = "rats")

tsk_rats$head()
```

Plotting the task with `autoplot` results in a `r index('Kaplan-Meier', lower = FALSE)` plot (@fig-autokm) which is a non-parametric estimator of the probability of survival for the average observation in the training set.

```{r, warning=FALSE, message=FALSE, output = FALSE, cache = FALSE}
autoplot(tsk_rats)
```
```{r, warning=FALSE, message=FALSE, echo = FALSE}
#| fig-cap: Kaplan-Meier plot of `tsk("rats")`. x-axis is time variable and y-axis is survival function, S(T), defined by $1 -$ F(T) where F is the cumulative distribution function. Crosses indicate points where censoring takes place.
#| fig-alt: Figure shows a line plot with "Time" on the x-axis from 0 to 100 and "Survival" on the y-axis from 0.80 to 1.00. The line plot is a black line from (0, 1) to (25, 1) then starts to drop slowly and then quickly down to (100, 0.80).
#| label: fig-autokm
p = autoplot(tsk_rats)
p$layers[[2]]$aes_params$colour = "grey30"
p
```

As well as creating your own tasks, you can load any of the tasks shipped with `mlr3proba`:

```{r special-013}
as.data.table(mlr_tasks)[task_type == "surv"]
```

### LearnerSurv, PredictionSurv and Predict Types

The interface for `r ref("LearnerSurv")` and `r ref("PredictionSurv")` objects is identical to the regression and classification settings discussed in @sec-basics.
Similarly to these settings, survival learners are constructed with `r ref("lrn()")`.

`mlr3proba` has a different predict interface to `mlr3` as all possible types of prediction ('predict types') are returned when possible for all survival models -- i.e., if a model *can* compute a particular predict type then *it will be* returned in `PredictionSurv`.
The reason for this design decision is that all these predict types can be transformed to one another and it is therefore computationally simpler to return all at once instead of rerunning models to change predict type.
In survival analysis, the following predictions can be made:

* `response` -- Predicted survival time.
* `distr` -- Predicted survival distribution, either discrete or continuous.
* `lp` -- Linear predictor calculated as the fitted coefficients multiplied by the test data.
* `crank` -- Continuous risk ranking.

We will go through each of these prediction types in more detail and with examples to make them less abstract.
We will use `lrn("surv.coxph")`\index{Cox Proportional Hazards} trained on `tsk("rats")` as a running example, since for this model, all predict types except `response` can be computed.

```{r special-015}
tsk_rats = tsk("rats")
split = partition(tsk_rats)
prediction_cph = lrn("surv.coxph")$train(tsk_rats, split$train)$
  predict(tsk_rats, split$test)
prediction_cph
```

#### predict_type = "response" {.unnumbered .unlisted}

Counterintuitively for many, the `response` prediction of predicted survival times is the least common predict type in survival analysis.
The likely reason for this is due to the presence of censoring.
We rarely observe the true survival time for many observations and therefore it is unlikely any survival model can confidently make predictions for survival times.
This is illustrated in the code below.

In the example below we train and predict from a survival SVM\index{support vector machine!survival} (`lrn("surv.svm")`), note we use `type = "regression"` to select the algorithm that optimizes survival time predictions and `gamma.mu = 1e-3` is selected arbitrarily as this is a required parameter (this parameter should usually be tuned).
We then compare the predictions from the model to the true data.

```{r special-016}
library(mlr3extralearners)
prediction_svm = lrn("surv.svm", type = "regression", gamma = 1e-3)$
  train(tsk_rats, split$train)$predict(tsk_rats, split$test)
data.frame(pred = prediction_svm$response[1:3],
  truth = prediction_svm$truth[1:3])
```

As can be seen from the output, our predictions are all less than the true observed time, which means we know our model underestimated the truth.
However, because each of the true values are censored times, we have absolutely no way of knowing if these predictions are slightly bad or absolutely terrible, (i.e., the true survival times could be $105, 99, 92$ or they could be $300, 1000, 200$).
Hence, with no realistic way to evaluate these models, survival time predictions are rarely useful.

#### predict_type = "distr" {.unnumbered .unlisted}

Unlike regression in which deterministic/point predictions are most common, in survival analysis distribution predictions are much more common.
You will therefore find that the majority of survival models in `mlr3proba` will make distribution predictions by default.
These predictions are implemented using the `r ref_pkg("alan-turing-institute/distr6")` package, which allows visualization and evaluation of survival curves (defined as $1 -$ cumulative distribution function).
Below we extract the first three `$distr` predictions from our running example and calculate the probability of survival at $t = 77$.

```{r special-017}
prediction_cph$distr[1:3]$survival(77)
```

The output indicates that there is a `r paste0(round(prediction_cph$distr[1:3]$survival(77)*100, 1), "%")`, chance of the first three predicted rats being alive at time 77 respectively.

#### predict_type = "lp" {.unnumbered .unlisted}

`lp`, often written as $\eta$ in academic writing, is computationally the simplest prediction and has a natural analog in regression modeling.
Readers familiar with linear regression will know that when fitting a simple linear regression model, $Y = X\beta$, we are estimating the values for $\beta$, and the estimated `r index("linear predictor")` (lp) is then $X\hat{\beta}$, where $\hat{\beta}$ are our estimated coefficients.
In simple survival models, the linear predictor is the same quantity (but estimated in a slightly more complicated way).
The learner implementations in `mlr3proba` are primarily machine-learning focused and few of these models have a simple linear form, which means that `lp` cannot be computed for most of these.
In practice, when used for prediction, `lp` is a proxy for a relative risk/continuous ranking prediction, which is discussed next.

#### predict_type = "crank" {.unnumbered .unlisted}

The final prediction type, `crank`, is the most common in survival analysis and perhaps also the most confusing.
Academic texts will often refer to 'risk' predictions in survival analysis (hence why survival models are often known as 'risk prediction models'), without defining what 'risk' means.
Often, risk is defined as $\exp(\eta)$ as this is a common quantity found in simple linear survival models.
However, sometimes risk is defined as $\exp(-\eta)$, and sometimes it can be an arbitrary quantity that does not have a meaningful interpretation.
To prevent this confusion in `mlr3proba`, we define the predict type `crank`, which stands for **c**ontinuous **rank**ing.
This is best explained by example; continuing from the previous we output the first three `crank` predictions.

```{r special-018}
prediction_cph$crank[1:3]
```

The output tells us that the first rat is at the lowest risk of death (smaller values represent lower risk) and the third rat is at the highest risk.
The distance between predictions also tells us that the difference in risk between the second and third rats is smaller than the difference between the first and second.
The actual values themselves are meaningless and therefore comparing `crank` values between samples (or papers or experiments) is not meaningful.

The `crank` prediction type is informative and common in practice because it allows identifying observations at lower/higher risk to each other, which is useful for resource allocation, e.g., which patient should be given an expensive treatment, and clinical trials, e.g., are people in a treatment arm at lower risk of disease X than people in the control arm.

:::{.callout-warning}
## Interpreting Survival Risk
The interpretation of 'risk' for survival predictions differs across R packages and sometimes even between models in the same package.
In `mlr3proba` there is one consistent interpretation of `crank`: lower values represent a lower risk of the event taking place and higher values represent higher risk.
:::

### MeasureSurv

Survival models in `mlr3proba` are evaluated with `r ref("MeasureSurv")` objects, which are constructed in the usual way with `msr()`.

In general survival measures can be grouped into the following:

1. Discrimination measures -- Quantify if a model correctly identifies if one observation is at higher risk than another. Evaluate `crank` and/or `lp` predictions.
2. Calibration measures -- Quantify if the average prediction is close to the truth (all definitions of calibration are unfortunately vague in a survival context). Evaluate `crank` and/or `lp` predictions.
3. Scoring rules -- Quantify if probabilistic predictions are close to true values. Evaluate `distr` predictions.

```{r special-020}
as.data.table(mlr_measures)[
  task_type == "surv", c("key", "predict_type")][1:5]
```

There is not a consensus in the literature around the 'best' survival measures to use to evaluate models.
We recommend ISBS (Integrated Survival Brier Score) (`msr("surv.graf")`) to evaluate the quality of `distr` predictions, concordance index (`msr("surv.cindex")`) to evaluate a model's discrimination, and D-Calibration (`msr("surv.dcalib")`) to evaluate a model's calibration.

Using these measures, we can now evaluate our predictions from the previous example.

```{r}
prediction_cph$score(msrs(c("surv.graf", "surv.cindex", "surv.dcalib")))
```

The model's performance seems okay as the ISBS and DCalib are relatively low and the C-index is greater than 0.5 however it is very hard to determine the performance of any survival model without comparing it to some baseline (usually the Kaplan-Meier).

### Composition {#sec-surv-comp}

Throughout `mlr3proba` documentation we refer to "native" and "composed" predictions.
We define a 'native' prediction as the prediction made by a model without any post-processing, whereas a 'composed' prediction is returned after post-processing.

#### Internal Composition

`mlr3proba` makes use of composition internally to return a `"crank"` prediction for every learner.
This is to ensure that we can meaningfully benchmark all models according to at least one criterion (discrimination performance).
The package uses the following rules to create `"crank"` predictions:

1. If a model returns a 'risk' prediction then `crank = risk` (we may multiply this by $-1$ to ensure the 'low-value low-risk' interpretation).
1. Else if a model returns a `response` prediction then we set `crank = -response`.
1. Else if a model returns a `lp` prediction then we set `crank = lp` (or `crank = -lp` if needed).
1. Else if a model returns a `distr` prediction then we set `crank` as the sum of the cumulative hazard function (see @Sonabend2022 for full discussion as to why we picked this method).

#### Explicit Composition and Pipelines

At the start of this section, we mentioned that it is possible to transform prediction types between each other.
In `mlr3proba` this is possible with 'compositor' pipelines (@sec-pipelines).
There are several pipelines implemented in the package but three in particular focus on predict type transformation:

1. `r ref("pipeline_crankcompositor()")` -- Transforms a `"distr"` prediction to `"crank"`
2. `r ref("pipeline_distrcompositor()")` -- Transforms a `"lp"` prediction to `"distr"`
3. `r ref("pipeline_responsecompositor()")` -- Transforms a `"distr"` prediction to `"response"` (survival time)

We internally use a version of the first pipeline whenever we return predictions from survival models so that every model has a `"crank"` prediction type - so only use the first pipeline to overwrite these ranking predictions.
In practice, the second pipeline is more common as Cox or Accelerated Failure Time (AFT) type models always return a linear predictor (`"lp"`), but sometimes the internal `predict()` functions don't provide a transformation to a survival distribution prediction (`"distr"`).
The third pipeline summarizes the predicted survival curves to a single number (expected survival time), and as previously mentioned, are rarely useful for evaluating the performance of survival machine learning models.

In the example below we load the `rats` dataset, remove factor columns, and then partition the data into training and testing.
We construct the `distrcompositor` pipeline around a survival XGBoost Accelerated Failure Time (AFT) learner (`lrn("surv.xgboost.aft")`) which by default makes predictions for `"lp"`, `"crank"` and `"response"`.
In the pipeline, we specify that we will estimate the baseline distribution with a `r index("Kaplan-Meier", lower = FALSE)` estimator (`estimator = "kaplan"`) and that we want to assume an AFT form for our estimated distribution (`form = "aft"`).
We then train and predict in the usual way and in our output we can now see a `distr` prediction.

```{r special-021, warning=FALSE}
library(mlr3verse)
library(mlr3extralearners)

tsk_rats = tsk("rats")$select(c("litter", "rx"))
split = partition(tsk_rats)

learner = lrn("surv.xgboost.aft", nrounds = 10)

# no distr output
learner$train(tsk_rats, split$train)$predict(tsk_rats, split$test)

graph_learner = ppl(
  "distrcompositor",
  learner = learner,
  estimator = "kaplan",
  form = "aft",
  graph_learner = TRUE
)

# now with distr
graph_learner$train(tsk_rats, split$train)$predict(tsk_rats, split$test)
```

Mathematically, we have done the following:

1. Assume our estimated distribution will have the form $S(t) = S_0(\frac{t}{\exp(\eta)})$ where $S$ is the survival function, $S_0$ is the baseline survival function and $\eta$ is the linear predictor.
1. Estimate $\hat{\eta}$ prediction using XGBoost
1. Estimate $\hat{S}_0(t)$ with the Kaplan-Meier estimator
1. Put this all together as $S(t) = \hat{S}_0(\frac{t}{\exp(\hat{\eta})})$

For more detail about prediction types and composition we recommend @Kalbfleisch2011.


### Putting It All Together {#sec-survival-all}

Finally, we will put all the above into practice in a small benchmark experiment.
We first load `tsk("grace")` (which only has numeric features) and sample 500 rows randomly.
We then select the ISBS, D-Calibration, and C-index to evaluate predictions, set up the same pipeline we used in the previous experiment, and load a Cox PH and Kaplan-Meier estimator.
We run our experiment with three-fold CV and aggregate the results.

```{r special-022, warning=FALSE}
set.seed(42)
library(mlr3extralearners)

tsk_grace = tsk("grace")
tsk_grace$filter(sample(tsk_grace$nrow, 500))
msr_txt = c("surv.graf", "surv.cindex", "surv.dcalib")
measures = msrs(msr_txt)

graph_learner = ppl(
  "distrcompositor",
  learner = lrn("surv.xgboost.aft", nrounds = 10),
  estimator = "kaplan",
  form = "aft",
  graph_learner = TRUE,
  scale_lp = TRUE
)
graph_learner$id = "XGBoost-AFT"
learners = c(lrns(c("surv.coxph", "surv.kaplan")), graph_learner)

bmr = benchmark(benchmark_grid(tsk_grace, learners,
  rsmp("cv", folds = 3)))
bmr$aggregate(measures)[, c("learner_id", ..msr_txt)]
```

In this small experiment, XGBoost-AFT and Cox PH have the best discrimination, the Kaplan-Meier baseline has the best calibration, and Cox PH has the best overall predictive accuracy (with the lowest ISBS).

## Density Estimation {#sec-density}

`r index('Density estimation')` is a learning task to estimate the unknown distribution from which a univariate dataset is generated or put more simply to estimate the probability density (or mass) function for a single variable.
As with survival analysis, density estimation is implemented in `mlr3proba`, as both can make probability distribution predictions (hence the name "**mlr3proba**bilistic").
Unconditional density estimation (i.e. estimation of a target without any covariates) is viewed as an unsupervised task, which means the 'truth' is never known.
For a good overview of density estimation see @Silverman1986.

The package `mlr3proba` extends `mlr3` with the following objects for density estimation:

* `r ref("TaskDens")` to define density tasks.
* `r ref("LearnerDens")` as the base class for density estimators.
* `r ref("PredictionDens")` for density predictions.
* `r ref("MeasureDens")` as a specialized class for density performance measures.

We will consider each in turn.

### TaskDens

As density estimation is an unsupervised task, there is no target for prediction.
In the code below we construct a density task using `r ref("as_task_dens()")` which takes one argument, a `data.frame` type object with exactly one column (which we will use to estimate the underlying distribution).

```{r special-023}
tsk_dens = as_task_dens(data.table(x = rnorm(1000)))
tsk_dens
```

As with other tasks, we have included a couple of tasks that come shipped with `mlr3proba`:

```{r}
as.data.table(mlr_tasks)[task_type == "dens", c(1:2, 4:5)]
```

### LearnerDens and PredictionDens

Density learners may return the following prediction types:

1. `distr` -- probability distribution
1. `pdf` -- probability density function
1. `cdf` -- cumulative distribution function

All learners will return a `distr` and `pdf` prediction but only some can make `cdf` predictions.
Again, the `distr` predict type is implemented using `distr6`.
In the code below we train and 'predict' with a histogram learner and then plot the estimated probability density function (@fig-dens-hist), which closely matches the underlying Normally-distributed data.

```{r special-024}
#| label: fig-dens-hist
#| fig-cap: Predicted density from the histogram learner, which closely resembles the underlying N(0, 1) data.
#| fig-alt: Image shows a line plot with x-axis between (-2,2) and y-axis between (0,0.4). The plot closely resembles a Normal(0, 1) distribution with a peak at 0.4.
lrn_hist = lrn("dens.hist")
prediction = lrn_hist$train(tsk_dens, 1:900)$predict(tsk_dens, 901:1000)
x = seq.int(-2, 2, 0.01)
df = data.frame(x = x, y = prediction$distr$pdf(x))
ggplot(df, aes(x = x, y = y)) + geom_line() + theme_minimal()
```

The `pdf` and `cdf` predict types are simply wrappers around `distr$pdf` and `distr$cdf` respectively:

```{r special-025}
prediction = lrn_hist$train(tsk_dens, 1:10)$predict(tsk_dens, 11:13)
# pdf and cdf columns in output
prediction
# comparing cdf from prediction to $cdf method from distr
cbind(prediction$distr$cdf(tsk_dens$data()$x[11:13]),
  prediction$cdf[1:3])
```

### MeasureDens and Putting It All Together

At the time of publication, the only measure implemented in `mlr3proba` for density estimation is logloss, which is defined in the same way as in classification, $L(y) = -\log(\hat{f}_Y(y))$, where $\hat{f}_Y$ is our estimated probability density function.
Putting this together with the above we are now ready to train a density learner, estimate a distribution, and evaluate our estimation:

```{r special-026}
msr_logloss = msr("dens.logloss")
msr_logloss

prediction$score(msr_logloss)
```

This output is most easily interpreted when compared to other learners in a benchmark experiment, so let us put everything together to conduct a small benchmark study on `tsk("faithful")` task using some of the integrated density learners:

```{r special-028, message=FALSE, warning=FALSE, results='hide'}
library(mlr3extralearners)
tsk_faithful = tsk("faithful")
learners = lrns(c("dens.hist", "dens.pen", "dens.kde"))
measure = msr("dens.logloss")
bmr = benchmark(benchmark_grid(tsk_faithful, learners,
  rsmp("cv", folds = 3)))
bmr$aggregate(measure)
```

```{r}
#| output: false
#| cache: false
autoplot(bmr, measure = measure)
```

```{r}
#| echo: false
#| warning: false
#| message: false
#| label: fig-beyond-density
#| fig-cap: Three boxplots comparing performance of dens.hist, dens.pen, and dens.kde on `tsk("faithful")`.
#| fig-alt: Three boxplots labeled 'dens.hist', 'dens.pen', 'dens.kde'. y-axis is 'dens.logloss' between 0.9 and 1.2. 'dens.hist' is a narrow boxplot between 1.0 and 1.1. 'dens.pen' is very wide between 0.95 and 2.1. 'dens.kde' is narrow between 0.95 and 1.0.
plt = ggplot2::last_plot()
plt$layers[[1]]$mapping$fill = NULL
print(plt)
```

The results (@fig-beyond-density) of this experiment indicate that the sophisticated Penalized Density Estimator does not outperform the baseline histogram, but the Kernel Density Estimator has at least consistently better (i.e. lower) logloss results.

## Cluster Analysis {#sec-cluster}

`r index("Cluster analysis")` is another unsupervised task implemented in `mlr3`.
The objective of cluster analysis is to group data into clusters, where each cluster contains similar observations.
The similarity is based on specified metrics that are task and application-dependent.
Unlike classification where we try to predict a class for each observation, in cluster analysis there is no 'true' label or class to predict.

The package `r mlr3cluster` extends `mlr3` with the following objects for cluster analysis:

* `r ref("TaskClust")` to define clustering tasks
* `r ref("LearnerClust")` as the base class for clustering learners
* `r ref("PredictionClust")` as the specialized class for `r ref("Prediction")` objects
* `r ref("MeasureClust")` as the specialized class for performance measures

We will consider each in turn.

### TaskClust

Similarly to density estimation (@sec-density), there is no target for prediction and so no `truth` field in `r ref("TaskClust")`.
By example, we will look at the `r ref("cluster::ruspini")` dataset, which has 75 rows and two columns and was first introduced in @Ruspini1970 to illustrate different clustering techniques.
The observations in the dataset form four natural clusters (@fig-beyond-clust-ruspini).
In the code below we construct a cluster task using `r ref("as_task_clust()")` which only takes one argument, a `data.frame` type object.
```{r special-030, warning=FALSE, message=FALSE}
#| label: fig-beyond-clust-ruspini
#| fig-cap: Distribution of the `ruspini` dataset.
#| fig-alt: "Four grids. Top-left shows a curve increasing sharply between (0,0.003) and (30,0.012) then decreasing to (120, 0.003). Top-right just says 'Corr: 0.065'. Bottom-left shows four distinct clusters of points. Bottom-right increases from (0, 50) to (140, 150) then decreases to (155, 100)."
library(mlr3verse)
library(cluster)
tsk_ruspini = as_task_clust(ruspini)
tsk_ruspini
tsk_ruspini$data(1:3) # print first 3 rows

autoplot(tsk_ruspini)
```

Technically, we did not need to create a new task for the `ruspini` dataset since it is already included in the package, along with one other task:

```{r special-031}
as.data.table(mlr_tasks)[task_type == "clust", c(1:2, 4:5)]
```

### LearnerClust and PredictionClust

As with density estimation, we refer to `training` and `predicting` for clustering to be consistent with the `mlr3` interface, but strictly speaking, this should be `clustering` and `assigning` (the latter we will return to shortly).
Two `predict_types` are available for clustering learners:

1. `"partition"` -- Estimate of which cluster an observation falls into
1. `"prob"` -- Probability of an observation belonging to each cluster

Similarly to classification, prediction types of clustering learners are either deterministic (`"partition"`) or probabilistic (`"prob"`).

Below we construct a C-Means clustering learner with `"prob"` prediction type and three clusters (`centers = 3`), train it on the `ruspini` dataset and then return the cluster assignments (`$assignments`) for six random observations.

```{r special-032}
lrn_cmeans = lrn("clust.cmeans", predict_type = "prob", centers = 3)
lrn_cmeans

lrn_cmeans$train(tsk_ruspini)
lrn_cmeans$assignments[sample(tsk_ruspini$nrow, 6)]
```

As clustering is unsupervised, it often does not make sense to use `predict` for new data however this is still possible using the `mlr3` interface.

```{r special-033, warning = FALSE, message = FALSE}
#| label: fig-beyond-clust-ruspini-estimated
#| fig-cap: Distribution of the estimated clusters.
#| fig-alt: "Four grids. Top-left shows three overlapping curves in purple (cluster 1), blue (cluster 2) and green (cluster 3). The purple and blue curves are zero in most places but then peak at (30, 120) and (60, 120) respectively. The green curve starts at (0,0) then increases slowly to (40, 120) then decreases bumpily to (120, 60). Top-right says '-0.78' in green (cluster 3), '0' in blue (cluster 2), and '-0.05' in purple (cluster 1). Bottom-left shows four distinct clusters of points, two clusters are green, one (bottom) is blue, one (bottom left) is purple. Bottom-right: line graphs that show a similar but inverted shape as top-left."

# using different data for estimation (rare use case)
lrn_cmeans$train(tsk_ruspini, 1:30)$predict(tsk_ruspini, 31:32)

# using same data as for estimation (common use case)
prediction = lrn_cmeans$train(tsk_ruspini)$predict(tsk_ruspini)
autoplot(prediction, tsk_ruspini)
```

While two prediction types are possible, there are some learners where 'prediction' can never make sense, for example in `r index("hierarchical clustering", aside = TRUE)`.
In hierarchical clustering, the goal is to build a hierarchy of nested clusters by either splitting large clusters into smaller ones or merging smaller clusters into bigger ones.
The final result is a tree or `r index('dendrogram')` which can change if a new data point is added.
For consistency, `mlr3cluster` offers a `predict` method for hierarchical clusters but with a warning:

```{r special-034}
#| label: fig-beyond-clust-dend
#| fig-cap: Dendrogram representing hierarchical clustering of the `ruspini` dataset. y-axis is similarity of points such that the lower observations (x-axis) are connected, the greater their similarity. The top split represents the separation of the two clusters.
#| fig-alt: "Plot shows a horizontal line that connects two vertical lines. Each vertical line connects to another horizontal line that splits into two more vertical lines, which continues for up to nine breaks."
lrn_hclust = lrn("clust.hclust", k = 2)
lrn_hclust$train(tsk_ruspini)$predict(tsk_ruspini)
autoplot(lrn_hclust) + theme(axis.text = element_text(size = 5.5))
```

In this case, the `predict` method simply cuts the dendrogram into the number of clusters specified by `k` parameter of the learner.


### MeasureClust

As previously discussed, unsupervised tasks do not have ground truth data to compare to in model evaluation.
However, we can still measure the quality of cluster assignments by quantifying how closely objects within the same cluster are related (`r index("cluster cohesion")`) as well as how distinct different clusters are from each other (`r index("cluster separation")`).
There are a few built-in evaluation metrics available to assess the quality of clustering, which can be found by searching the `r ref("mlr_measures")` dictionary.

Two common measures are the within sum of squares (WSS) measure (`msr("clust.wss")`) and the silhouette coefficient (`msr("clust.silhouette")`).
WSS calculates the sum of squared differences between observations and centroids, which is a quantification of cluster cohesion (smaller values indicate the clusters are more compact).
The silhouette coefficient quantifies how well each point belongs to its assigned cluster versus neighboring clusters, where scores closer to `1` indicate well clustered and scores closer to `-1` indicate poorly clustered.
Note that the silhouette measure in `mlr3cluster` returns the mean silhouette score across all observations and when there is only a single cluster, the measure simply outputs 0.

Putting this together with the above we can now score our cluster estimation (note we must pass the `task` to `$score`):

```{r special-036}
measures = msrs(c("clust.wss", "clust.silhouette"))

prediction$score(measures, task = tsk_ruspini)
```

The very high WSS and middling mean silhouette coefficient indicate that our clusters could do with a bit more work.

Often reducing an unsupervised task to a quantitative measure may not be useful (given no ground truth) and instead visualization (discussed next) may be a more effective tool for assessing the quality of the clusters.

### Visualization

As clustering is an unsupervised task, visualization can be essential not just for 'evaluating' models but also for determining if our learners are performing as expected for our task. This section will look at visualizations for supporting clustering choices and following that we will consider plots for evaluating model performance.

#### Visualizing Clusters

It is easy to rely on clustering measures to assess the quality of clustering however this should be done with care as choosing between models may come down to other decisions such as how clusters are formed.
By example, consider data generated by `r ref("mlbench::mlbench.spirals")`, which results in two individual lines that spiral around each other (@fig-beyond-clust-spirals).

```{r special-037}
#| label: fig-beyond-clust-spirals
#| fig-cap: Distribution of `spirals` data.
#| fig-alt: "Grid of four plots. Top-left: line increasing from (-1,0.1) to (0,0.5) then decreasing to (1,0.1). Top-right: text that says 'Corr: -0.145'. Bottom-left: two lines of dots that are in tight, non-overlapping spirals around each other. Bottom-right: same shape as top-left."
spirals = mlbench::mlbench.spirals(n = 300, sd = 0.01)
tsk_spirals = as_task_clust(as.data.frame(spirals$x))
autoplot(tsk_spirals)
```

Now let us see what happens when fit two clustering learners on this data:

```{r special-038}
learners = list(
  lrn("clust.kmeans"),
  lrn("clust.dbscan", eps = 0.1)
)

bmr = benchmark(benchmark_grid(tsk_spirals, learners, rsmp("insample")))
bmr$aggregate(msr("clust.silhouette"))[, c(4, 7)]
```

We can see that K-means clustering gives us a higher average silhouette score and so we might conclude that a K-means learner with two centroids is a better choice than the DBSCAN method.
However, now take a look at the cluster assignment plots in @fig-beyond-clust-spirals-pred (`autoplot.PredictionClust` is available but we do not use it here so we can highlight two particular plots).

```{r special-039, message=FALSE, warning=FALSE}
#| label: fig-beyond-clust-spirals-pred
#| fig-cap: Comparing estimated clusters from `lrn("clust.kmeans")` and `lrn("clust.dbscan")`. Both create two distinct clusters that are separated in different ways.
#| fig-alt: "Two plots of the same spirals as in the previous plot. Left (K-means): points above the line x=y are purple (cluster 1) and other points are green (cluster 2). Right (DBSCAN): One of the spirals is purple and the other is green."
library(patchwork)
# get K-Means and DBSCAN partitions
pred_kmeans = as.factor(bmr$resample_result(1)$prediction()$partition)
pred_dbscan = as.factor(bmr$resample_result(2)$prediction()$partition)
# plot
df_kmeans = cbind(tsk_spirals$data(), clust = pred_kmeans)
df_dbscan = cbind(tsk_spirals$data(), clust = pred_dbscan)
map = aes(x = V1, y = V2, color = clust)
p_kmeans = ggplot(df_kmeans, map) + ggtitle("K-means")
p_dbscan = ggplot(df_dbscan, map) + ggtitle("DBSCAN")

p_kmeans + p_dbscan + plot_layout(guides = "collect") & geom_point() &
  theme_minimal() & ggplot2::scale_colour_viridis_d(end = 0.8)
```

The two learners arrived at two different results to cleanly separate clusters -- the K-means algorithm assigned points that are part of the same line into two different clusters whereas DBSCAN assigned each line to its own cluster.
Which one of these approaches is correct?
The answer is it depends on your specific task and the goal of cluster analysis.
If we had only relied on the silhouette score, then the details of how the clustering was performed would have been masked and we would have been unable to decide which method was appropriate for the task.

#### PCA and Silhouette Plots

The two most important plots implemented in `r mlr3viz` to support the evaluation of cluster learners are PCA and silhouette plots.

`r index('Principal components analysis')` (PCA) is a commonly used dimension reduction method in ML to reduce the number of variables in a dataset or to visualize the most important 'components', which are linear transformations of the dataset features.
Components are considered more important if they have higher variance (and therefore more predictive power).
In the context of clustering, by plotting observations against the first two components, and then coloring them by cluster, we could visualize our high-dimensional dataset and we would expect to see observations in distinct groups.

Since our running example only has two features, PCA does not make sense to visualize the data.
So we will use a task based on the `USArrests` dataset instead.
By plotting the result of PCA (@fig-beyond-clust-usarrests), we see that our model has done a good job of separating observations into two clusters along the first two principal components.

```{r special-040, message=FALSE, warning=FALSE}
#| label: fig-beyond-clust-usarrests
#| fig-cap: First two principal components using PCA on `tsk("usarrests")`.
#| fig-alt: "Scatter plot of green (cluster 2) and purple (cluster 1) points. x-axis: PC1 (96.55%) between -0.3 and 0.2. y-axis: PC2 (2.78%) between -0.3 and 0.2. Points are cleanly separated into two clusters by color."
tsk_usarrests = tsk("usarrests")
prediction = lrn("clust.kmeans")$train(tsk_usarrests)$
  predict(tsk_usarrests)
autoplot(prediction, tsk_usarrests, type = "pca")
```

Silhouette plots visually assess the quality of the estimated clusters by visualizing if observations in a cluster are well-placed both individually and as a group.
The plots include a dotted line which visualizes the average silhouette coefficient across all data points and each data point's silhouette value is represented by a bar colored by their assigned cluster.
In our particular case, the average silhouette index is `r round(msr("clust.silhouette")$score(prediction, tsk_usarrests), 2)`.
If the average silhouette value for a given cluster is below the average silhouette coefficient line then this implies that the cluster is not well defined.

Continuing with our new example, we find (@fig-beyond-clust-sil) that a lot of observations are actually below the average line and close to zero, and therefore the quality of our cluster assignments is not very good, meaning that many observations are likely assigned to the wrong cluster.

```{r special-041, message=FALSE, warning=FALSE}
#| label: fig-beyond-clust-sil
#| fig-cap: 'Silhouette plot from predictions made by `lrn("clust.kmeans")` on `tsk("usarrests")`.'
#| fig-alt: "Horizontal barplot with 'Silhouette Values' on x-axis between 0 and 1; y-axis is 'Observations' between 0 and 50. Observations between 0-20 are all colored purple (cluster 1) and observations between 21-50 are colored green (cluster 2). A dashed vertical line passes through x=0.59. The majority of bars finish before this line."
autoplot(prediction, tsk_usarrests, type = "sil")
```

### Putting It All Together {#sec-cluster-all}

Finally, we conduct a small benchmark study using `tsk("usarrests")` and a few integrated cluster learners:

```{r special-042, message=FALSE, warning=FALSE}
tsk_usarrests = tsk("usarrests")
learners = list(
  lrn("clust.featureless"),
  lrn("clust.kmeans", centers = 4L),
  lrn("clust.cmeans", centers = 3L)
)
measures = list(msr("clust.wss"), msr("clust.silhouette"))
bmr = benchmark(benchmark_grid(tsk_usarrests, learners,
  rsmp("insample")))
bmr$aggregate(measures)[, c(4, 7, 8)]
```

The C-means and K-means algorithms are both considerably better than the featureless baseline but further analysis (and visualizations) would be required to decide which of those two is suitable for our needs.

## Spatial Analysis {#sec-spatiotemporal}

The final task we will discuss in this book is `r index("spatial analysis")`.
Spatial analysis can be a subset of any other machine learning task (e.g., regression or classification) and is defined by the presence of spatial information in a dataset, usually stored as coordinates that are often named "x" and "y" or "lat" and "lon" (for 'latitude' and 'longitude' respectively.)

Spatial analysis is its own task as spatial data must be handled carefully due to the complexity of 'autocorrelation'.
Where `r index("correlation")` is defined as a statistical association *between two* variables, `r index("autocorrelation", aside = TRUE)` is a statistical association *within one* variable.
In ML terms, in a dataset with features and observations, correlation occurs when two or more features are statistically associated in some way, whereas autocorrelation occurs when two or more observations are statistically associated across one feature.
Autocorrelation, therefore, violates one of the fundamental assumptions of ML that all observations in a dataset are independent, which results in lower confidence about the quality of a trained machine learning model and the resulting performance estimates [@hastie2001].

Autocorrelation is present in spatial data as there is implicit information encoded in coordinates, such as whether two observations (e.g., cities, countries, continents) are close together or far apart.
By example, let us imagine we are predicting the number of cases of a disease two months after an outbreak in Germany (@fig-autocorrelation).
Outbreaks radiate outwards from an epicenter and therefore countries closer to Germany will have higher numbers of cases and countries further away will have lower numbers (@fig-autocorrelation, bottom).
Thus, looking at the data spatially shows clear signs of autocorrelation across nearby observations.
Note in this example the autocorrelation is radial but in practice, this will not always be the case.

```{r special-043, warning=FALSE, echo=FALSE, message=FALSE, out.width = "60%", out.height = "60%", fig.height = 5}
#| label: fig-autocorrelation
#| fig-cap: Heatmaps where darker countries indicate higher number of cases and lighter countries indicate lower number of cases of imaginary Disease X with epicenter in Germany. The top map imagines a world in which there is no spatial autocorrelation and the number of cases of a disease is randomly distributed. The bottom map shows a more accurate world in which the number of cases radiate outwards from the epicenter (Germany).
#| fig-alt: Image shows two separate maps of Europe. Top map has a random distribution of colors from white to dark gray. Bottom map shows darkest color (dark gray) at Germany with increasing lightness as the countries are increasingly further away.

library(dplyr)
library(patchwork)

m = map_data("world")
m = m[m$region == "Germany", ]
x1 = mean(m$lon)
x2 = mean(m$lat)
m = map_data("world")
m = m %>% dplyr::filter(long > -15, long < 40, lat > 36, lat <  60)
m$dist = (-sqrt((m$lon - x1)^2 + (m$lat - x2)^2)) * 10e5
if (min(m$dist) < 0) {
  m$dist = m$dist - min(m$dist)
}
m$rand = runif(nrow(m), min(m$dist), max(m$dist))
cou = m %>% group_by(region) %>% summarize(long = mean(long), lat = mean(lat)) %>% filter(region == "Germany")

p = ggplot(data = m) +
  scale_fill_gradient2(low = "white", mid = "gray", high = "black", midpoint = mean(m$dist)) +
  theme_minimal() +
  xlim(-10, 25) + ylim(37, 58)

t = geom_text(aes(label = region, x = long, y = lat), data = cou, size = 2.5)

p1 = p + geom_map(map = m, aes(long, lat, map_id = region, fill = rand), color = "black", lwd = 0.6) + t +
  theme(legend.position = "n")

p2 = p + geom_map(map = m, aes(long, lat, map_id = region, fill = dist), color = "black", lwd = 0.6) + t +
  labs(fill = "Cases")

p1 / p2 & patchwork::plot_layout(guides = "collect")
```

Unlike other tasks we have looked at in this chapter, there is no underlying difference between the implemented learners or measures.
Instead, we provide additional resampling methods in `r mlr3spatiotempcv` to account for the similarity in the train and test sets during resampling that originates from spatiotemporal autocorrelation.

Throughout this section we will use the `r ref("mlr3spatiotempcv::ecuador")` dataset and task as a working example.

### TaskClassifST and TaskRegrST

To make use of spatial resampling methods, we have implemented two extensions of `r ref("TaskClassif")` and `r ref("TaskRegr")` to accommodate spatial data, `r ref("TaskClassifST")` and `r ref("TaskRegrST")` respectively.
Below we only show classification examples but regression follows trivially.

```{r special-044, message=FALSE, warning=FALSE}
library(mlr3spatial)
library(mlr3spatiotempcv)

# create task from `data.frame`
tsk_ecuador = as_task_classif_st(ecuador, id = "ecuador_task",
  target = "slides", positive = "TRUE",
  coordinate_names = c("x", "y"), crs = "32717")

# or create task from 'sf' object
data_sf = sf::st_as_sf(ecuador, coords = c("x", "y"), crs = "32717")
tsk_ecuador = as_task_classif_st(data_sf, target = "slides",
  positive = "TRUE")
tsk_ecuador
```

Once a task is created, you can train and predict as normal.

```{r special-045}
lrn("classif.rpart")$train(tsk_ecuador)$predict(tsk_ecuador)
```

However as discussed above, it is best to use the specialized resampling methods to achieve bias-reduced estimates of model performance.

### Spatiotemporal Cross-Validation {#spatiotemp-cv}

Before we look at the spatial resampling methods implemented in `mlr3spatiotempcv` we will first show what can go wrong if non-spatial resampling methods are used for spatial data.
Below we benchmark a decision tree on `tsk("ecuador")` using two different repeated cross-validation resampling methods, the first ("NSpCV" (non-spatial cross-validation)) is a non-spatial resampling method from `mlr3`, the second ("SpCV" (spatial cross-validation)) is from `mlr3spatiotempcv` and is optimized for spatial data.
The example highlights how "NSpCV" makes it appear as if the decision tree is performing better than it is with considerably higher estimated performance, however, this is an overconfident prediction due to the autocorrelation in the data.

```{r special-046, warning=FALSE,message=FALSE}
lrn_rpart = lrn("classif.rpart", predict_type = "prob")
rsmp_nsp = rsmp("repeated_cv", folds = 3, repeats = 2, id = "NSpCV")
rsmp_sp = rsmp("repeated_spcv_coords", folds = 3, repeats = 2,
  id = "SpCV")

design = benchmark_grid(tsk_ecuador, lrn_rpart, c(rsmp_nsp, rsmp_sp))
bmr = benchmark(design)
bmr$aggregate(msr("classif.acc"))[, c(5, 7)]
```

In the above example, applying non-spatial resampling results in train and test sets that are very similar due to the underlying spatial autocorrelation.
Hence there is little difference from testing a model on the same data it was trained on, which should be avoided for an honest performance result (see @sec-basics).
In contrast, the spatial method has accommodated autocorrelation and the test data is less correlated (though some association will remain) with the training data.
Visually this can be seen using `r ref("mlr3spatiotempcv::autoplot()")` methods.
In @fig-sprsmp we visualize how the task is partitioned according to the spatial resampling method (@fig-sprsmp, left) and non-spatial resampling method (@fig-sprsmp, right).
There is a clear separation in space for the respective partitions when using the spatial resampling whereas the train and test splits overlap a lot (and are therefore more correlated) using the non-spatial method.

```{r special-047}
#| fig-cap: Scatterplots show separation of train (blue) and test (orange) data for the first fold of the first repetition of the cross-validation. Left is spatial resampling where train and test data are clearly separated. Right is non-spatial resampling where there is overlap in train and test data.
#| fig-alt: Two scatter plots with points in blue (training data) and orange (test data). Left plot (Spatial Resampling) shows a clean separation between orange and blue points. Right plot (Non-spatial Resampling) shows blue and orange dots randomly scattered among each other.
#| label: fig-sprsmp

library(patchwork)

(autoplot(rsmp_sp, tsk_ecuador, fold_id = 1, size = 0.7) +
  ggtitle("Spatial Resampling") +
  autoplot(rsmp_nsp, tsk_ecuador, fold_id = 1, size = 0.7) +
  ggtitle("Non-spatial Resampling")) +
  plot_layout(guides = "collect") &
  theme_minimal() &
  theme(axis.text = element_text(size = 4), legend.position = "bottom")
```

Now we have seen why spatial resampling matters we can take a look at what methods are available in `mlr3spatiotempcv`.
The resampling methods we have added can be categorized into:

* Blocking -- Create rectangular blocks in 2D or 3D space
* Buffering -- Create buffering zones to remove observations between train and test sets
* Spatiotemporal clustering -- Clusters based on coordinates (and/or time points)
* Feature space clustering -- Clusters based on feature space and not necessarily spatiotemporal
* Custom (partitioning) -- Grouped by factor variables

The choice of method may depend on specific characteristics of the dataset and there is no easy rule to pick one method over another, full details of different methods can be found in @Schratz2021 -- the paper deliberately avoids recommending one method over another as the 'optimal' choice is highly dependent on the predictive task, autocorrelation in the data, and the spatial structure of the sampling design.
The documentation for each of the implemented methods includes details of each method as well as references to original publications.

:::{.callout-tip}
## Spatio*temporal* Resampling

We have focused on spatial analysis but referred to "spatiotemporal" and "spatiotemp".
The spatial-only resampling methods discussed in this section can be extended to temporal analysis (or spatial and temporal analysis combined) by setting the `"time"` `col_role` in the task (@sec-row-col-roles) -- this is an advanced topic that may be added in future editions of this book.
See the `mlr3spatiotempcv` visualization vignette at `r link("https://mlr3spatiotempcv.mlr-org.com/articles/spatiotemp-viz.html")` for specific details about 3D spatiotemporal visualization.
:::

### Spatial Prediction {#sec-spatial-prediction}

Until now we have looked at resampling to accommodate spatiotemporal *features*, but what if you want to make spatiotemporal *predictions*?
In this case, the goal is to make classification or regression predictions at the pixel level, i.e., for an area, defined by the geometric resolution, of a raster image.

To enable these predictions we have created a new function, `r ref("mlr3spatial::predict_spatial()")`, to allow spatial predictions using any of the following spatial data classes:

* `stars` (from package `r ref_pkg("stars")`)
* `SpatRaster` (from package `r ref_pkg("terra")`)
* `RasterLayer` (from package `r ref_pkg("raster")`)
* `RasterStack` (from package `r ref_pkg("raster")`)

In the example below we load the `leipzig_points` dataset for training and coerce this to a spatiotemporal task with `r ref("as_task_classif_st")`, and we load the `leipzig_raster` raster.
Both files are included as example data in `r mlr3spatial`.

```{r special-048, warning=FALSE, message=FALSE}
#| cache: false
library(mlr3spatial)
library(sf)
library(terra, exclude = "resample")

# load sample points
leipzig_vector = sf::read_sf(system.file("extdata",
  "leipzig_points.gpkg", package = "mlr3spatial"),
  stringsAsFactors = TRUE)
# create training data
tsk_leipzig = as_task_classif_st(leipzig_vector, target = "land_cover")

# load raster image
leipzig_raster = terra::rast(system.file("extdata", "leipzig_raster.tif",
  package = "mlr3spatial"))
```

Now we can continue as normal to train and predict with a classification learner, in this case, a random forest.

```{r special-049, cache.lazy=FALSE, cache=FALSE}
lrn_ranger = lrn("classif.ranger")$train(tsk_leipzig)
prediction = predict_spatial(leipzig_raster, lrn_ranger,
  format = "terra")
prediction
```

In this example, we specified the creation of a `terra` object, which can be visualized with in-built plotting methods.

```{r special-050, message = FALSE, cache.lazy=FALSE, cache=FALSE, out.width = "50%", out.height = "50%"}
#| fig-cap: Spatial predictions for forest (purple), pasture (blue), urban (green), and water (yellow) categories.
#| fig-alt: Very zoomed-in map with x-axis from 732000 to 733000 and 5692500 to 5693500 on y-axis. Different clusters are colored in green, blue, purple and yellow.
#| label: fig-beyond-raster
plot(prediction, col = c("#440154FF", "#443A83FF", "#31688EFF",
  "#21908CFF", "#35B779FF", "#8FD744FF", "#FDE725FF"))
```

## Conclusion

In this chapter, we explored going beyond regression and classification to see how classes in `mlr3` can be used to implement other ML tasks.
Cost-sensitive classification extends the 'normal' classification setting by assuming that costs associated with `r index('false negatives', 'false negative')` and `r index('false positives', 'false positive')` are unequal.
Running cost-sensitive classification experiments is possible using only features in `mlr3`.
Survival analysis, available in `r mlr3proba`, can be thought of as a regression problem when the outcome may be censored, which means it may never be observed within a given time frame.
The final task in `mlr3proba` is density estimation, the unsupervised task concerned with estimating univariate probability distributions.
Using `r mlr3cluster`, you can perform cluster analysis on observations, which involves grouping observations according to similarities in their variables.
Finally, with `r mlr3spatial` and `r mlr3spatiotempcv`, it is possible to perform spatial analysis to make predictions using coordinates as features and to make spatial predictions.
The `mlr3` interface is highly extensible, which means future ML tasks can (and will) be added to our universe and we will add these to this chapter of the book in future editions.

| Class | Constructor/Function | Fields/Methods |
| ----- | --- | -- |
| `r ref("mlr3::MeasureClassifCosts")` | `msr("classif.costs")` | - |
| `r ref("mlr3pipelines::PipeOpTuneThreshold")` | `po("tunethreshold")` | - |
| `r ref("mlr3proba::TaskSurv")` | `r ref("as_task_surv()")` | `$data()` |
| `r ref("mlr3proba::PipeOpDistrCompositor")` | `po("distrcompose")` | - |
| `r ref("mlr3proba::TaskDens")` | `r ref("as_task_dens()")` | `$data()` |
| `r ref("mlr3cluster::TaskClust")` | `r ref("as_task_clust()")` | `$data()` |
| `r ref("mlr3spatiotempcv::TaskClassifST")` | `r ref("as_task_classif_st()")` | `$data()` |
| - | `r ref("mlr3spatiotempcv::predict_spatial()")` |  |

: Important classes and functions covered in this chapter with underlying class (if applicable), class constructor or function, and important class fields and methods (if applicable). {#tbl-beyond-api}

## Exercises

1. Run a benchmark experiment on `tsk("german_credit")` with `lrn("classif.featureless")`, `lrn("classif.log_reg")`, and `lrn("classif.ranger")`.
  Tune the prediction thresholds of all learners by encapsulating them in a `po("learner_cv")` (with two-fold CV), followed by a `po("tunethreshold")`.
  Use `msr("classif.costs", costs = costs)`, where the `costs` matrix is as follows: true positive is `-10`, true negative is `-1`, false positive is `2`, and false negative is `3`.
  Use this measure in `po("tunethreshold")` and when evaluating your benchmark experiment.
2. Train and test a survival forest using `lrn("surv.rfsrc")` (from `mlr3extralearners`).
  Run this experiment using `tsk("rats")` and `partition()`.
  Evaluate your model with the RCLL measure.
3. Estimate the density of the "precip" task from the `mlr3proba` package using `lrn("dens.hist")`, evaluate your estimation with the logloss measure.
  As a stretch goal, look into the documentation of `distr6` to learn how to analyse your estimated distribution further.
4. Run a benchmark clustering experiment on the "wine" dataset without a label column.
  Compare the performance of k-means learner with `k` equal to `2`, `3` and `4` using the silhouette measure and the insample resampling technique.
  What value of `k` would you choose based on the silhouette scores?

::: {.content-visible when-format="html"}
`r citeas(chapter)`
:::