---
title: "Introduction to the MachineShop Package"
subtitle: "Package Version `r packageVersion('MachineShop')`"
author:
- name: "Brian J Smith"
  affiliation: "University of Iowa"
  email: "brian-j-smith@uiowa.edu"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{Introduction to the MachineShop Package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
source("setup.R")
rdoc_url <- function(name) name
```


# Description

`MachineShop` is a meta-package for statistical and machine learning with a unified interface for model fitting, prediction, performance assessment, and presentation of results.  Support is provided for predictive modeling of numerical, categorical, and censored time-to-event outcomes and for resample (bootstrap, cross-validation, and split training-test sets) estimation of model performance.  This vignette introduces the package interface with a survival data analysis example, followed by supported methods of variable specification; applications to other response variable types; available performance metrics, resampling techniques, and graphical and tabular summaries; and modeling strategies.


# Features

```{r overview_modelinfo, echo=FALSE}
library("MachineShop")
info <- modelinfo()
```

* Unified and concise interface for model fitting, prediction, and performance assessment.
* Current support for `r length(info)` established models from `r length(unique(sapply(info, function(x) x$packages[1])))` **R** packages.
* Dynamic model parameters.
* Ensemble modeling with stacked regression and super learners.
* Modeling of response variables types: binary factors, multi-class nominal and ordinal factors, numeric vectors and matrices, and censored time-to-event survival.
* Model specification with traditional formulas, design matrices, and flexible pre-processing [recipes](https://cran.r-project.org/package=recipes).
* Resample estimation of predictive performance, including cross-validation, bootstrap resampling, and split training-test set validation.
* Parallel execution of resampling algorithms.
* Choices of performance metrics: accuracy, areas under ROC and precision recall curves, Brier score, coefficient of determination (R^2^), concordance index, cross entropy, F score, Gini coefficient, unweighted and weighted Cohen's kappa, mean absolute error, mean squared error, mean squared log error, positive and negative predictive values, precision and recall, and sensitivity and specificity.
* Graphical and tabular performance summaries: calibration curves, confusion matrices, partial dependence plots, performance curves, lift curves, and variable importance.
* Model tuning over automatically generated grids of parameter values and randomly sampled grid points.
* Model selection and comparisons for any combination of models and model parameter values.
* User-definable models and performance metrics.


# Getting Started


## Installation

```{r overview_install, eval = FALSE}
# Current release from CRAN
install.packages("MachineShop")

# Development version from GitHub
# install.packages("devtools")
devtools::install_github("brian-j-smith/MachineShop")

# Development version with vignettes
devtools::install_github("brian-j-smith/MachineShop", build_vignettes = TRUE)
```


## Documentation

Once installed, the following **R** commands will load the package and display its help system documentation.  Online documentation and examples are available at the [MachineShop website](https://brian-j-smith.github.io/MachineShop).

```{r overview_docs, eval = FALSE, message = FALSE}
library(MachineShop)

# Package help summary
?MachineShop

# Vignette
RShowDoc("Introduction", package = "MachineShop")
```


# Melanoma Example

The package is illustrated in the following sections with an overall survival analysis example in which the response variable is a time to event outcome.  Since survival outcomes are a combination of numerical (time to event) and categorical (event) variables, package features for both variable types are illustrated in the example.  Outcomes other than survival, including nominal and ordinal factors as well as numeric vectors and matrices, are supported by `MachineShop` and will be discussed.

Survival analysis is performed with the `Melanoma` dataset from the `MASS` package [@andersen:1993:SMB].  The dataset provides survival time, in days, from disease treatment to (1) death from disease, (2) alive at end of study, or (3) death from other causes for 205 Denmark patients with malignant melanomas.  Also provided are potential predictors of the survival outcomes.  The analysis begins by loading required packages `MachineShop`, `survival`, and `MASS` as well as `magrittr` [@bache:2014:MFP] for its pipe (`%>%`) operator to simplify some of the code syntax.  For the analysis, a binary overall survival outcome is created by combining the two death categories (1 and 3) into one. 

```{r using_example_melanoma}
## Analysis libraries
library(MachineShop)
library(survival)
library(MASS)
library(magrittr)

## Malignant melanoma analysis dataset
surv_df <- within(Melanoma, status <- as.numeric(status != 2))
```

Descriptive summaries of the study variables are given below in Table 1, followed by a plot of estimated overall survival probabilities and 95% confidence intervals.

<center>
Table 1. Variable summaries for the Melanoma survival analysis example.
</center>

```{r using_example_summary, echo=FALSE}
surv_stats <- list(
  list("Number of subjects" = ~ length(status)),
  "time" = list("Median (Range)" = ~ median_range(time)),
  "status" = list("1 = Dead" = ~ n_perc(status == 1),
                  "0 = Alive" = ~ n_perc(status == 0)),
  "sex" = list("1 = Male" = ~ n_perc(sex == 1),
               "0 = Female" = ~ n_perc(sex == 0)),
  "age" = list("Median (Range)" = ~ median_range(age)),
  "year" = list("Median (Range)" = ~ median_range(year)),
  "thickness" = list("Median (Range)" = ~ median_range(thickness)),
  "ulcer" = list("1 = Presence" = ~ n_perc(ulcer == 1),
                 "0 = Absence" = ~ n_perc(ulcer == 0))
)

summary_kbl(surv_stats, surv_df)
```

```{r using_example_survfit, echo=FALSE}
library(ggplot2)

col <- "#F8766D"
survfit(Surv(time, status) ~ 1, data = surv_df) %>%
  with(data.frame(time, surv, lower, upper, censor = ifelse(n.censor > 0, time, NA))) %>%
  ggplot(aes(x = time, y = surv)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = col, alpha = 0.2) +
  geom_step(color = col) +
  geom_point(aes(x = censor), shape = 3, color = col) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = "Follow-Up Time (Days)", y = "Overall Survival Probability",
       title = "Kaplan-Meier survival plot")
```

For the analyses, the dataset is split into a training set on which a survival model will be fit and a test set on which predictions will be made and performance evaluated.  A global formula `surv_fo` is defined to relate the predictors on the right hand side to the overall survival outcome on the left and will be used for all subsequent survival models.

```{r using_example_datasets}
## Training and test sets
set.seed(123)
train_indices <- sample(nrow(surv_df), nrow(surv_df) * 2 / 3)
surv_train <- surv_df[train_indices, ]
surv_test <- surv_df[-train_indices, ]

## Global formula for the analysis
surv_fo <- Surv(time, status) ~ sex + age + year + thickness + ulcer
```


# Model Fit and Prediction


## Model Information

Model fitting requires user specification of a `MachineShop` compatible model.  A named list of package-supplied models can be obtained interactively with the `modelinfo` function, and includes the following components for each.

`label`
  : Character descriptor for the model.
  
`packages`
  : Character vector of source packages required to use the model.  These need only be installed with the `install.packages` function or by equivalent means; but need not be loaded with, for example, the `library` function.
  
`response_types`
  : Character vector of response variable types supported by the model.
  
`arguments`
  : Closure with the argument names and corresponding default values of the model function.
  
`grid`
  : Logical indicating whether automatic generation of tuning parameter grids is implemented for the model.
  
`varimp`
  : Logical indicating whether variable importance is defined for the model.

Function `modelinfo` may be called without arguments, with one or more metric functions, an observed response variable, an observed and predicted response variable pair, response variable types, or resampled output; and will return information on all matching metrics.

```{r using_fit_modelinfo}
## All available models
modelinfo() %>% names

## Model-specific information
modelinfo(GBMModel)
```

Information is displayed above for the `GBMModel` function --- a generalized boosted regression model that can be applied to survival outcomes.


### Type-Specific

When data objects are supplied as arguments, information is returned on all models applicable to response variables of the same data types.  If model functions are additionally supplied as arguments, information on the subset matching the data types is returned. 

```{r using_fit_modelinfo_type}
## All survival response-specific models
modelinfo(Surv(0)) %>% names

## Identify survival response-specific models
modelinfo(Surv(0), CoxModel, GBMModel, SVMModel) %>% names
```


### Response Variable-Specific

As a special case of type-specific arguments, existing response variables to be used in analyses may be given as arguments to identify applicable models.

```{r using_fit_modelinfo_response}
## Models for a responses variable
modelinfo(Surv(surv_df$time, surv_df$status)) %>% names
```


## Fit Function

Package models, like `GBMModel`, can be specified in the `model` argument of the `fit` function to estimate a relationship (`surv_fo`) between predictors and an outcome based on a set of data (`surv_train`).  Argument specifications may be in terms of the model function, function name, or a function call.

```{r using_fit_function, results="hide"}
## Generalized boosted regression fit

## Model function
surv_fit <- fit(surv_fo, data = surv_train, model = GBMModel)

## Model function name
fit(surv_fo, data = surv_train, model = "GBMModel")

## Model function call
fit(surv_fo, data = surv_train, model = GBMModel(n.trees = 100, interaction.depth = 1))
```


## Dynamic Model Parameters

*Dynamic model parameters* are model function arguments defined as expressions to be evaluated at the time of model fitting.  As such, their values can change based on the number of observations in the dataset supplied to the fit function.  Expressions to dynamic parameters are specified within the package-supplied quote operator `.()` and can include the following objects:

`nobs`
  : number of observations in `data`.
   
`nvars`
  : number of predictor variables in `data`.

`y`
  : response variable.
  
In the example below, Bayesian information criterion (BIC) based stepwise variable selection is performed by creating a `CoxStepAICModel` with dynamic parameter `k` to be calculated as the log number of observations in the fitted dataset.

```{r using_fit_dynamic, results="hide"}
## Dynamic model parameter k = log number of observations

## Number of observations: nobs
fit(surv_fo, data = surv_train, model = CoxStepAICModel(k = .(log(nobs))))

## Response variable: y
fit (surv_fo, data = surv_train, model = CoxStepAICModel(k = .(log(length(y)))))
```


## Predict Function

A `predict` function is supplied for application to model fit results to obtain predicted values on a dataset specified with its `newdata` argument or on the original dataset if not specified.  Survival means are predicted for survival outcomes by default.  Estimates of the associated survival distributions are needed to calculate the means.  For models, like `GBMModel`, that perform semi- or non-parametric survival analysis, Weibull approximations to the survival distributions are the default for mean estimation.  Other choices of distributional approximations are exponential, Rayleigh, and empirical.  Empirical distributions are applicable to Cox proportional hazards-based models and can be calculated with the method of Breslow [-@breslow:1972:DPC], Efron [-@efron:1977:ECL, default], or Fleming and Harrington [-@fleming:1984:NES].  Note, however, that empirical survival means are undefined mathematically if an event does not occur at the longest follow-up time.  In such situations, a restricted survival mean is calculated by changing the longest follow-up time to an event, as suggested by Efron [-@efron:1967:PFB], which will be negatively biased.  

```{r using_predict_function}
## Predicted survival means (default: Weibull distribution)
predict(surv_fit, newdata = surv_test) %>% head

## Predicted survival means (empirical distribution)
predict(surv_fit, newdata = surv_test, dist = "empirical") %>% head
```

In addition to survival means, predicted survival probabilities (`type = "prob"`) or 0-1 survival events (default: `type = "response"`) can be obtained with the follow-up `times` argument.  The cutoff probability for classification of survival events (or other binary responses) can be set optionally with the `cutoff` argument (default: `cutoff = 0.5`).  As with mean estimation, distributional approximations to the survival functions may be specified for the predictions, with the default for survival probabilities being the empirical distribution.

```{r using_predict_function_times}
## Predict survival probabilities and events at specified follow-up times
surv_times <- 365 * c(5, 10)

predict(surv_fit, newdata = surv_test, times = surv_times, type = "prob") %>% head

predict(surv_fit, newdata = surv_test, times = surv_times, cutoff = 0.5) %>% head
```


# Variable Specifications

Variable specification defines the relationship between response and predictor variables as well as the data used to estimate the relationship.  Four main types of specifications are supported by the `fit`, `resample`, and `tune` functions: traditional formula, design matrix, model frame, and recipe.


## Traditional Formula

Variables may be specified with a traditional formula and data frame pair, as was done at the start of the survival example.  This specification allows for allows for crossing (`*`), interaction (`:`), and removal (`-`) of predictors in the formula; `.` substitution of variables not already appearing in the formula; in-line functions of response variables; and in-lining of base **R** operators and math functions of predictors.

```{r using_variables_formula, results="hide"}
## Dataset library
library(MASS)

## Formula specification
fit(medv ~ ., data = Boston, model = GBMModel)
```


## Design Matrix

Variables stored separately in a design matrix of predictors and object of responses can be supplied to the fit functions directly.  Fitting with design matrices has less computational overhead than traditional formulas and allows for greater numbers of predictor variables in some models, including `GBMModel`, `GLMNetModel`, and `RandomForestModel`.

```{r using_variables_matrix, results="hide"}
## Example design matrix and response object
x <- model.matrix(medv ~ . - 1, data = Boston)
y <- Boston$medv

## Design matrix specification
fit(x, y, model = GBMModel)
```


## Model Frame

A `ModelFrame` class is provided by the package for specification of predictor and response variables along with other attributes to control model fitting.  Model frames can be specified in calls to the `ModelFrame` constructor function with a syntax similar to the traditional formula or design matrix.

```{r using_variables_modelframe, results="hide"}
## Model frame specification

## Formula
mf <- ModelFrame(medv ~ ., data = Boston)

fit(mf, model = GBMModel)

## Design matrix
mf <- ModelFrame(x, y)

fit(mf, model = GBMModel)
```

The model frame approach has a few advantages over model fitting directly with a traditional formula.  One is that cases with missing values on any of the response or predictor variables are excluded from the model frame by default.  This is often desirable for models that do not handle missing values.  Conversely, missing values can be retained in the model frame by setting its argument `na.rm = FALSE` for models, like `GBMModel`, that do handle them.  A second advantage is that case weights can be included in the model frame to be passed on to the model fitting functions.

```{r using_variables_modelframe_weights, results="hide"}
## Model frame specification with case weights
mf <- ModelFrame(ncases / (ncases + ncontrols) ~ agegp + tobgp + alcgp, data = esoph,
                 weights = with(esoph, ncases + ncontrols))

fit(mf, model = GBMModel)
```

A third, which will be illustrated later, is user-specification of a variable for stratified resampling via the constructor's `strata` argument.


## Preprocessing Recipe

The `recipes` package [@kuhn:2018:RPT] provides a flexible framework for defining predictor and response variables as well as preprocessing steps to be applied to them prior to model fitting.  Using recipes helps ensure that estimation of predictive performance accounts for all modeling step.  They are also a convenient way of consistently applying preprocessing to new data.  A basic recipe is given below in terms of the formula and data frame ingredients needed for the analysis.

```{r using_variables_recipe, results="hide"}
## Recipe specification
library(recipes)

rec <- recipe(medv ~ ., data = Boston)

fit(rec, model = GBMModel)
```

Case weights and stratified resampling are also supported for recipes via the designations of `"case_weight"` and `"case_strata"` roles, respectively.  As an example, an initial step is included in the recipe below to update (alter) the role of variable `weights` to be case weights.  That is followed by a step to convert three ordinal factors to integer scores.

```{r using_variables_recipe_weights, results="hide"}
## Recipe specification with case weights
df <- within(esoph, {
  y <- ncases / (ncases + ncontrols)
  weights <- ncases + ncontrols
})

rec <- recipe(y ~ agegp + tobgp + alcgp + weights, data = df) %>%
  update_role(weights, new_role = "case_weight") %>%
  step_ordinalscore(agegp, tobgp, alcgp)

fit(rec, model = GBMModel)
```


## Summary

The variable specification approaches differ with respect to support for preprocessing, in-line functions, case weights, resampling strata, and computational overhead, as summarized in the table below.  Only recipes apply preprocessing steps automatically during model fitting and should be used when it is important to account for such steps in the estimation of model predictive performance.  Preprocessing would have to be done manually and separately otherwise.  Design matrices have the lowest computational overhead and can enable analyses involving larger numbers of predictors than the other approaches.  Both recipes and model frames allow for user-defined case weights (default: equal) and resampling strata (default: none).  The remaining approaches default to equal weights and to strata defined by the response variable.  Syntax ranges from simplest to most complex for design matrices, traditional formulas, model frames, and recipes, respectively.  The relative strengths of each approach should be considered within the context of a given analysis when deciding upon which one to use.

```{r using_variables_summary, echo=FALSE}
df <- data.frame(
  "Specification" = c("Traditional Formula", "Design Matrix",
                      "Traditional Formula", "Design Matrix", "Recipe"),
  "Preprocessing" = factor(c("manual", "manual", "manual", "manual",
                             "automatic"), levels = c("manual", "automatic")),
  "In-line Functions" = factor(c("yes", "no", "yes", "no", "no"),
                               levels = c("no", "yes")),
  "Case Weights" = factor(c("equal", "equal", "user", "user", "user"),
                          levels = c("equal", "user")),
  "Resampling Strata" = factor(c("response", "response", "user", "user",
                                 "user"), levels = c("response", "user")),
  "Computational Overhead" = factor(c("medium", "low", "medium", "low", "high"),
                                    levels = c("high", "medium", "low")),
  check.names = FALSE
)

bg_colors <- c("orange", "blue", "green")
df[-1] <- lapply(df[-1], function(x) {
  bg_colors <- if (nlevels(x) == 2) bg_colors[c(1, 3)] else bg_colors
  cell_spec(x, color = "white", background = bg_colors[x])
})

kable(df, align = c("l", rep("c", ncol(df) - 1)), escape = FALSE,
      caption = "Table 2. Characteristics of available variable specification approaches.") %>%
  kable_styling(c("striped", "condensed"), full_width = FALSE,
                position = "center") %>%
  column_spec(1, bold = TRUE) %>%
  kableExtra::group_rows("Model Frame", 3, 4)
```


# Response Variable Types

The **R** class types of response variables play a central role in their analysis with the package.  They determine, for example, the specific models that can be fit, fitting algorithms employed, predicted values produced, and applicable performance metrics and analyses.  As described below, factors, ordered factors, numeric vectors and matrices, and survival responses are supported by the package.


## Factors

Categorical responses with two or more levels should be coded as `factor` variables for analysis.

```{r using_responses_factor, results="hide"}
## Iris flowers species (3-level factor)
fit(Species ~ ., data = iris, model = GBMModel)
```

```{r using_responses_factor_binary, results="hide"}
## Pima Indians diabetes statuses (binary factor)
library(MASS)

fit(type ~ ., data = Pima.tr, model = GBMModel)
```


## Ordered Factors

Categorical responses can be designated as having ordered levels by storing them as `ordered` factor variables.  For categorical vectors, this can be accomplished with the `factor` function and its argument `ordered = TRUE` or more simply with the `ordered` function.  Numeric vectors can be converted to ordered factors with the `cut` function.

```{r using_responses_ordered, results="hide"}
## Boston housing prices (ordered factor)
library(MASS)

df <- within(Boston, {
  medv <- cut(medv, breaks = 3, ordered_result = TRUE)
})

fit(medv ~ ., data = df, model = GBMModel)
```


## Numeric Vectors

Code univariate numerical responses as `numeric` variables.

```{r using_responses_numeric, results="hide"}
## Boston housing prices
library(MASS)

fit(medv ~ ., data = Boston, model = GBMModel)
```


## Numeric Matrices

Store multivariate numerical responses as numeric `matrix` variables for model fitting with traditional formulas or model frames.

```{r using_responses_matrix, results="hide"}
## Anscombe's multiple regression models dataset

## Numeric matrix response formula
fit(cbind(y1, y2, y3) ~ x1, data = anscombe, model = LMModel)
```

For recipes, include the multiple response variables on the left hand side of the formula specification.

```{r using_responses_matrix_recipe, results="hide"}
## Numeric matrix response recipe
rec <- recipe(y1 + y2 + y3 ~ x1, data = anscombe)

fit(rec, model = LMModel)
```


## Survival Objects

Censored time-to-event survival responses should be stored as `Surv` variables for model fitting with traditional formulas or model frames.

```{r using_responses_surv, results="hide"}
## Survival response formula
library(survival)

fit(Surv(time, status) ~ ., data = surv_train, model = GBMModel)
```

For recipes, survival responses should be specified with the individual survival time and event variables given on the left hand side of the formula and with their roles designated as `"surv_time"` and `"surv_event"`.

```{r using_responses_surv_recipe, results="hide"}
## Survival response recipe
rec <- recipe(time + status ~ ., data = surv_train) %>%
  add_role(time, new_role = "surv_time") %>%
  add_role(status, new_role = "surv_event")

fit(rec, model = GBMModel)
```


# Model Performance Metrics

Performance metrics quantify associations between observed and predicted responses and provide a means of assessing and comparing the predictive performances of models.


## Performance Function

Metrics can be computed with the `performance` function applied to observed responses and responses predicted with the `predict` function.  In the case of observed versus predicted survival probabilities or events, metrics will be calculated at each survival time and returned along with their time-integrated mean.

```{r using_performance_function}
## Survival performance metrics

## Observed responses
obs <- response(surv_fit, newdata = surv_test)

## Predicted survival means
pred_means <- predict(surv_fit, newdata = surv_test)
performance(obs, pred_means)

## Predicted survival probabilities
pred_probs <- predict(surv_fit, newdata = surv_test, times = surv_times, type = "prob")
performance(obs, pred_probs)

## Predicted survival events
pred_events <- predict(surv_fit, newdata = surv_test, times = surv_times)
performance(obs, pred_events)
```

Function `performance` computes a default set of metrics according to the observed and predicted response types, as indicated and in the order given in the table below.


Table 3. Default performance metrics by response types.

Response                 | Default Metrics
:----------------------- | :------------------------------------------------------
Factor                   | Brier Score, Accuracy, Cohen's Kappa
Binary Factor            | Brier Score, Accuracy, Cohen's Kappa, Area Under ROC Curve, Sensitivity, Specificity
Numeric Vector or Matrix | Root Mean Squared Error, R^2^, Mean Absolute Error
Survival Means           | Concordance Index
Survival Probabilities   | Area Under ROC Curve, Brier Score, Accuracy
Survival Events          | Accuracy


These defaults may be changed by specifying one or more package-supplied metric functions to the `metrics` argument of `performance`.  Specification of the `metrics` argument can be in terms of a single metric function, function name, or list of metric functions.  List names, if specified, will be displayed as metric labels in graphical and tabular summaries; otherwise, the function names will be used as labels for unnamed lists.

```{r using_performance_function_metrics, eval=FALSE}
## Single metric function
performance(obs, pred_means, metrics = cindex)

## Single metric function name
performance(obs, pred_means, metrics = "cindex")

## List of metric functions
performance(obs, pred_means, metrics = c(cindex, rmse, rmsle))

## Named list of metric functions
performance(obs, pred_means, metrics = c("CIndex" = cindex,
                                         "RMSE" = rmse,
                                         "RMSLE" = rmsle))
```

Metrics based on classification of two-level class probabilities, like sensitivity and specificity, optionally allow for specification of the classification cutoff probability (default: `cutoff = 0.5`).

```{r using_performance_function_cutoff}
## User-specified survival probability metrics
performance(obs, pred_probs, metrics = c(sensitivity, specificity), cutoff = 0.5)
```


## Metric Functions

Whereas multiple package-supplied metrics can be calculated simultaneously with the `performance` function, each exists as a stand-alone function that can be called individually.

```{r using_metrics_functions}
## Metric functions for survival means
cindex(obs, pred_means)

rmse(obs, pred_means)

rmsle(obs, pred_means)

## Metric functions for survival probabilities
sensitivity(obs, pred_probs)

specificity(obs, pred_probs)
```


## Metric Information

A named list of available metrics can be obtained interactively with the `metricinfo` function, and includes the following components for each one.

`label`
  : Character descriptor for the metric.
  
`maximize`
  : Logical indicating whether higher values of the metric correspond to better predictive performance.
  
`arguments`
  : Closure with the argument names and corresponding default values of the metric function.
  
`response_types`
  : Data frame of the observed and predicted response variable types supported by the metric.


Function `metricinfo` may be called without arguments, with one or more metric functions, an observed response variable, an observed and predicted response variable pair, response variable types, or resampled output; and will return information on all matching metrics.

```{r using_metrics_metricinfo}
## All available metrics
metricinfo() %>% names

## Metric-specific information
metricinfo(cindex)
```


### Type-Specific

When data objects are supplied as arguments, information is returned on all metrics applicable to response variables of the same data types.  Observed response variable type is inferred from the first data argument and predicted type from the second, if given.  For survival responses, predicted types may be `numeric` for survival means, `SurvEvents` for 0-1 survival events at specified follow-up times, or `SurvProbs` for follow-up time survival probabilities.  If model functions are additionally supplied as arguments, information on the subset matching the data types is returned. 

```{r using_metrics_metricinfo_type}
## Metrics for observed and predicted response variable types
metricinfo(Surv(0)) %>% names

metricinfo(Surv(0), numeric(0)) %>% names

metricinfo(Surv(0), SurvEvents(0)) %>% names

metricinfo(Surv(0), SurvProbs(0)) %>% names

## Identify survival-specific metrics
metricinfo(Surv(0), auc, cross_entropy, cindex) %>% names
```


### Response Variable-Specific

Existing response variables observed and those obtained from the `predict` function may be given as arguments to identify metrics that are applicable to them.

```{r using_metrics_metricinfo_response}
## Metrics for observed and predicted responses from model fits
metricinfo(obs, pred_means) %>% names

metricinfo(obs, pred_probs) %>% names
```


## Factors

Metrics applicable to multi-level factor response variables are summarized below.

`accuracy`
  : Proportion of correctly classified responses.

`brier`
  : [Brier score](https://en.wikipedia.org/wiki/Brier_score).
  
`cross_entropy`
  : [Cross entropy](https://en.wikipedia.org/wiki/Cross_entropy) loss averaged over the number of cases.
  
`kappa2`
  : [Cohen's kappa](https://en.wikipedia.org/wiki/Cohen%27s_kappa) statistic measuring relative agreement between observed and predicted classifications.

`weighted_kappa2`
  : [Weighted Cohen's kappa](https://en.wikipedia.org/wiki/Cohen%27s_kappa#Weighted_kappa).  This metric is only available for ordered factor responses.
  
Brier score and cross entropy loss are computed directly on predicted class probabilities.  The other metrics are computed on predicted class membership, defined as the factor level with the highest predicted probability.


## Binary Factors

Metrics for binary factors include those given for multi-level factors as well as the following.

`auc`
  : Area under a performance curve.

`cindex`
  : Concordance index computed as rank order agreement between predicted probabilities for paired event and non-event cases.  This metric can be interpreted as the probability that a randomly selected event case will have a higher predicted value than a randomly selected non-event case, and is the same as area under the ROC curve.

`f_score`
  : [F score](https://en.wikipedia.org/wiki/Precision_and_recall#F-measure), $F_\beta = (1 + \beta^2) \frac{\text{precision} \times \text{recall}}{\beta^2 \times \text{precision} + \text{recall}}$.  F1 score $(\beta = 1)$ is the package default.

`fnr`
  : False negative rate, $FNR = \frac{FN}{TP + FN} = 1 - TPR$.

```{r using_metrics_conf, echo=FALSE}
conf <- matrix(c("True Negative (TN)", "False Positive (FP)",
                 "False Negative (FN)", "True Positive (TP)"),
               2, 2,
               dimnames = list("Predicted Response" = c("Negative", "Positive"),
                               "Observed Response" = c("Negative", "Positive")))
kable(conf,
      caption = "Table 4. Confusion matrix of observed and predicted response classifications.",
      align = c("c", "c")) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  add_header_above(c("Predicted Response" = 1, "Observed Response" = 2))
```

`fpr`
  : False positive rate, $FPR = \frac{FP}{TN + FP} = 1 - TNR$.

`npv`
  : Negative predictive value, $NPV = \frac{TN}{TN + FN}$.

`ppv`, `precision`
  : Positive predictive value, $PPV = \frac{TP}{TP + FP}$.

`pr_auc`, `auc`
  : Area under a precision recall curve.
  
`roc_auc`, `auc`
  : [Area under an ROC curve](https://en.wikipedia.org/wiki/Receiver_operating_characteristic#Area_under_the_curve).

`roc_index`
  : A tradeoff function of sensitivity and specificity as defined by the `f` argument in this function (default: sensitivity + specificity).  The function allows for specification of tradeoffs [@perkins:2006:IOC] other than the default of Youden's J statistic [@youden:1950:IRD].
  
`rpp`
  : Rate of positive prediction, $RPP = \frac{TP + FP}{TP + FP + TN + FN}$.
  
`sensitivity`, `recall`, `tpr`
  : True positive rate, $TPR =\frac{TP}{TP + FN} = 1 - FNR$.
  
`specificity`, `tnr`
  : True negative rate, $TNR = \frac{TN}{TN + FP} = 1 - FPR$.

Area under the ROC and precision-recall curves as well as the concordance index are computed directly on predicted class probabilities.  The other metrics are computed on predicted class membership.  Memberships are defined to be in the second factor level if predicted probabilities are greater than the function default or user-specified cutoff value.


## Numerics

Performance metrics are defined below for numeric vector responses.  If applied to a numeric matrix response, the metrics are computed separately for each column and then averaged to produce a single value.

`gini`
  : [Gini coefficient](https://en.wikipedia.org/wiki/Gini_coefficient).
  
`mae`
  : Mean absolute error, $MAE = \frac{1}{n}\sum_{i=1}^n|y_i - \hat{y}_i|$, where $y_i$ and $\hat{y}_i$ are the $n$ observed and predicted responses.

`mse`
  : Mean squared error, $MSE = \frac{1}{n}\sum_{i=1}^n(y_i - \hat{y}_i)^2$.

`msle`
  : Mean squared log error, $MSLE = \frac{1}{n}\sum_{i=1}^n(log(1 + y_i) - log(1 + \hat{y}_i))^2$.

`r2`
  : One minus residual divided by total sums of squares, $R^2 = 1 - \sum_{i=1}^n(y_i - \hat{y}_i)^2 / \sum_{i=1}^n(y_i - \bar{y})^2$.

`rmse`
  : Square root of mean squared error.

`rmsle`
  : Square root of mean squared log error.


## Survival Objects

All previously described metrics for binary factor responses---plus accuracy, Brier score and Cohen's kappa---are applicable to survival probabilities predicted at specified follow-up times.  Metrics are evaluated separately at each follow-up time and reported along with a time-integrated mean.  The survival concordance index is computed with the method of Harrell [-@harrell:1982:EYM] and Brier score according to Graf et al. [-@graf:1999:ACP]; whereas, the others are computed according to the confusion matrix probabilities below, in which term $\hat{S}(t)$ is the predicted survival probability at follow-up time $t$ and $T$ is the survival time [@heagerty:2004:TDR].

```{r using_metrics_conf_surv, echo=FALSE}
conf <- matrix(c("$TN = \\Pr(\\hat{S}(t) \\gt \\text{cutoff} \\cap T \\ge t)$",
                 "$FP = \\Pr(\\hat{S}(t) \\le \\text{cutoff} \\cap T \\ge t)$",
                 "$FN = \\Pr(\\hat{S}(t) \\gt \\text{cutoff} \\cap T \\lt t)$",
                 "$TP = \\Pr(\\hat{S}(t) \\le \\text{cutoff} \\cap T \\lt t)$"),
               2, 2,
               dimnames = list("Predicted Response" = c("Non-Event", "Event"),
                               "Observed Response" = c("Non-Event", "Event")))
kable(conf,
      caption = "Table 5. Confusion matrix of observed and predicted survival response classifications.",
      align = c("c", "c"),
      escape = FALSE) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  add_header_above(c("Predicted Response" = 1, "Observed Response" = 2))
```

In addition, all of the metrics described for numeric vector responses are applicable to predicted survival means and are computed using only those cases with observed (non-censored) events.


# Resample Performance Estimation


## Algorithms

Model performance can be estimated with resampling methods that simulate repeated training and test set fits and predictions.  With these methods, performance metrics are computed on each resample to produce an empirical distribution for inference.  Resampling is controlled in the `MachineShop` with the functions:

BootControl
  : Simple bootstrap resampling [@efron:1993:IB].  Models are repeatedly fit with bootstrap resampled training sets and used to predict the full data set.

BootOptimismControl
  : Optimism-corrected bootstrap resampling [@efron:1983:LLB; @harrell:1996:MPM].

CVControl
  : Repeated K-fold cross-validation [@kohavi:1995:SCB].  The full data set is repeatedly partitioned into K-folds.  For a given partitioning, prediction is performed on each of the K folds with models fit on all remaining folds.

OOBControl
  : Out-of-bootstrap resampling.  Models are fit with bootstrap resampled training sets and used to predict the unsampled cases.
  
SplitControl
  : Split training and test sets [@hastie:2009:ESL7].  The data are randomly partitioned into a training and test set.
  
TrainControl
  : Training resubstitution.  A model is fit on and used to predict the full training set in order to estimate training, or apparent, error [@efron:1986:HBA].
  
For the survival example, repeated cross-validation control structures are defined to estimate model performance in predicting survival means and 5 and 10-year survival probabilities.  In addition to arguments controlling the resampling algorithms, a `seed` can be set to ensure reproducibility of resampling results obtained with the structures.

```{r using_resample_control}
## Control parameters for K-fold cross-validation

## Prediction of survival means
surv_means_control <- CVControl(folds = 5, repeats = 3, seed = 123)

## Prediction of survival probabilities
surv_probs_control <- CVControl(folds = 5, repeats = 3, times = surv_times, seed = 123)
```


## Parallel Processing

Resampling is implemented with the `foreach` package [@microsoft:2017:FPF] and will run in parallel if a compatible backend is loaded, such as that provided by the `doParallel` package [@microsoft:2017:DFP].

```{r using_resample_parallel}
## Register multiple cores for parallel computations
library(doParallel)
registerDoParallel(cores = 2)
```


## Resample Function

Resampling is performed by calling the `resample` function with a variable specification, model, and control structure.  Like the `fit` function, variables may be specified in terms of a traditional formula, design matrix, model frame, or recipe.

```{r using_resample_function}
## Resample estimation for survival means and probabilities
(res_means <- resample(surv_fo, data = surv_train, model = GBMModel, control = surv_means_control))

(res_probs <- resample(surv_fo, data = surv_train, model = GBMModel, control = surv_probs_control))
```


## Summary Statistics

The `summary` function when applied directly to output from `resample` computes summary statistics for the default performance metrics described in the *Performance Function* section.

```{r using_resample_summary}
## Summary of survival means metric
summary(res_means)

## Summary of survival probability metrics
summary(res_probs)
```

Other relevant metrics can be identified and summarized with the `metricinfo` and `performance` functions.

```{r using_resample_summary_performance}
## Resample-specific metrics
metricinfo(res_means) %>% names

## User-specified survival means metrics
summary(performance(res_means, metrics = c(cindex, rmse)))
```

Futhermore, summaries can be customized with a user-defined statistics function or list of statistics functions passed to the `stats` argument of `summary`.

```{r using_resample_summary_stats}
## User-defined statistics function
percentiles <- function(x) quantile(x, probs = c(0.25, 0.50, 0.75))
summary(res_means, stats = percentiles)

## User-defined list of statistics functions
summary(res_means, stats = c(Mean = mean, Percentile = percentiles))
```


## Plots

Summary plots of resample output can be obtained with the `plot` function.  Boxplots are the default plot type; but density, errorbar, and violin plots are also available.  Plots are generated with the `ggplot2` package [@wickham:2016:GEG] and returned as `ggplot` objects.  As such, annotation and formatting defined for ggplots can be applied to the returned plots.

```{r using_resample_plots}
## Libraries for plot annotation and fomatting
library(ggplot2)
library(gridExtra)

## Individual ggplots
p1 <- plot(res_means)
p2 <- plot(res_means, type = "density")
p3 <- plot(res_means, type = "errorbar")
p4 <- plot(res_means, type = "violin")

## Grid of plots
grid.arrange(p1, p2, p3, p4, nrow = 2)
```


## Stratified Resampling

Stratification of cases for the construction of resampled training and test sets can be employed to help achieve balance across the sets.  Stratified resampling is automatically performed if variable specification is in terms of a traditional formula and will be done according to the response variable if a numeric vector or factor, the event variable if survival, and the first variable if a numeric matrix.  For model frames and recipes, stratification variables must be defined explicitly with the `strata` argument to the `ModelFrame` constructor or with the `"case_strata"` role designation in a recipe step. 

```{r using_resample_strata, results="hide"}
## Model frame with case status stratification
mf <- ModelFrame(surv_fo, data = surv_train, strata = surv_train$status)

resample(mf, model = GBMModel)

## Recipe with case status stratification
rec <- recipe(time + status ~ ., data = surv_train) %>%
  add_role(time, new_role = "surv_time") %>%
  add_role(status, new_role = "surv_event") %>%
  add_role(status, new_role = "case_strata")

resample(rec, model = GBMModel)
```


## Dynamic Model Parameters

As discussed previously in the *Model Fit and Prediction* section, dynamic model parameters are evaluated at the time of model fitting and can depend on the number of observations in the fitted dataset.  In the context of resampling, dynamic parameters are repeatedly evaluated at each fit of the resampled datasets.  As such, their values can change based on the observations selected for training at each iteration of the resampling algorithm.

```{r using_resample_dynamic, results="hide"}
## Dynamic model parameter k = log number of training set observations
resample(surv_fo, data = surv_train, model = CoxStepAICModel(k = .(log(nobs))))
```


## Model Comparisons

Resampled metrics from different models can be combined for comparison with the `Resamples` function.  Optional names given on the left hand side of equal operators within calls to `Resamples` will be used as labels in output from the `summary` and `plot` functions.  For comparisons of resampled output, the same control structure must be used in all associated calls to `resample` to ensure that resulting model metrics are computed on the same resampled training and test sets.  The combined resample output can be summarized and plotted as usual.

```{r using_resample_comparisons}
## Resample estimation
res1 <- resample(surv_fo, data = surv_train, model = GBMModel(n.trees = 25),
                 control = surv_means_control)
res2 <- resample(surv_fo, data = surv_train, model = GBMModel(n.trees = 50),
                 control = surv_means_control)
res3 <- resample(surv_fo, data = surv_train, model = GBMModel(n.trees = 100),
                 control = surv_means_control)

## Combine resample output for comparison
(res <- Resamples(GBM1 = res1, GBM2 = res2, GBM3 = res3))

summary(res)

plot(res)
```

Pairwise model differences for each metric can be calculated with the `diff` function applied to results from a call to `Resamples`.  Resulting differences can be summarized descriptively with the `summary` and `plot` functions and assessed for statistical significance with the `t.test` function.

```{r using_resample_diff}
## Pairwise model comparisons
(perfdiff <- diff(res))

summary(perfdiff)

plot(perfdiff)
```

```{r using_resample_diff_test}
t.test(perfdiff)
```


# Performance Analyses

Calculation of performance metrics on test sets or by resampling, as discussed previously, is one method of assessing model performance.  Others available include measures of predictor variable importance, calibration curves comparing observed and predicted response values, partial dependence plots, and receiver operating characteristic analysis.


## Variable Importance

The importance of variables in a model fit is estimated with the `varimp` function and plotted with `plot`.  Variable importance is a relative measure of the contributions of model predictors and has a default range of 0 to 100, where 0 denotes the least important variables and 100 the most.  Classes of models can differ with respect to how variable importance is defined.  In the case of a `GBMModel`, importance of each predictor is based on the sum of squared empirical improvements over all internal tree nodes created by splitting on that variable [@greenwell:2019:GBM].

```{r using_analyses_vi}
## Predictor variable importance
(vi <- varimp(surv_fit))

plot(vi)
```

Alternatively, importance is based on one minus predictor variable p-values for statistical models, like `CoxModel`, that produce them.  For other classes of models, variable importance is generally defined and calculated by their underlying source packages.


## Calibration Curves

Agreement between model-predicted and observed values can be visualized with calibration curves.  In the construction of these curves, cases are partitioned into equally spaced bins according to their (resampled) predicted responses.  Mean observed responses are then calculated within each of the bins and plotted on the vertical axis against the bin midpoints on the horizontal axis. 

```{r using_analyses_cal, results="hide"}
## Binned calibration curves
cal <- calibration(res_probs, breaks = 10)
plot(cal, se = TRUE)
```

As an alternative to discrete bins, curves can be smoothed over the individual predicted values by setting `breaks = NULL`.

```{r using_analyses_cal_smoothed, results="hide"}
## Smoothed calibration curves
cal <- calibration(res_probs, breaks = NULL)
plot(cal)
```

Calibration curves that are close to the 45$^\circ$ line indicate close agreement between observed and predicted responses and a model that is said to be well calibrated.


## Confusion Matrices

Confusion matrices of cross-classified observed and predicted categorical responses are available with the `confusion` function.  They can be constructed with predicted class membership or with predicted class probabilities.  In the latter case, predicted class membership is derived from predicted probabilities according to a probability cutoff value for binary factors (default: `cutoff = 0.5`) and according to the class with highest probability for factors with more than two levels.


```{r using_analyses_conf}
## Confusion matrices
(conf <- confusion(res_probs, cutoff = 0.5))
```

```{r using_analyses_conf_plot, results="hide"}
plot(conf)
```

Confusion matrices are the data structure upon which many of the performance metrics described earlier for factor predictor variables are based.  Metrics commonly reported for confusion matrices are generated by the `summary` function.

```{r using_analyses_conf_summary}
## Summary performance metrics
summary(conf)
```

Summaries can also be obtained with `performance` function but for select metrics specified by the user.

```{r using_analyses_conf_performance}
## Confusion matrix-specific metrics
metricinfo(conf) %>% names

## User-specified metrics
performance(conf, metrics = c("Accuracy" = accuracy,
                              "Sensitivity" = sensitivity,
                              "Specificity" = specificity))
```


## Partial Dependence Plots

Partial dependence plots display the marginal effects of predictors on a response variable.  Dependence for a select set of one or more predictor variables $X_S$ is computed as
$$
\bar{f}_S(X_S) = \frac{1}{N}\sum_{i=1}^N f(X_S, x_{iS'}),
$$
where $f$ is a fitted prediction function and $x_{iS'}$ are values of the remaining predictors in a dataset of $N$ cases.  The response scale displayed in dependence plots will depend on the response variable type: probability for predicted factors and survival probabilities, original scale for numerics, and survival time for predicted survival means.  By default, dependence is computed for each select predictor individually over a grid of 10 approximately evenly spaced values and averaged over the dataset on which the prediction function was fit.

```{r using_analyses_pd, results = "hide"}
## Partial dependence plots
pd <- dependence(surv_fit, select = c(thickness, age))
plot(pd)
```

Averaging may be performed over different datasets to estimate marginal effects in other populations of cases, over different numbers of predictor values, and over quantile spacing of the values.

```{r using_analyses_pd_data, results = "hide"}
pd <- dependence(surv_fit, data = surv_test, select = thickness, n = 20,
                 intervals = "quantile")
plot(pd)
```

In addition, dependence may be computed for combinations of multiple predictors to examine interaction effects and for summary statistics other than the mean.


## Performance Curves

Tradeoffs between correct and incorrect classifications of binary outcomes, across the range of possible cutoff probabilities, can be studied with performance curves.


### ROC

Receiver operating characteristic (ROC) curves are one example in which true positive rates (sensitivity) are plotted against false positive rates (1 - specificity) [@fawcett:2006:IRA].  Higher ROC curves are indicative of better predictive performance.

```{r using_analyses_roc}
## ROC curves
roc <- performance_curve(res_probs)
plot(roc, diagonal = TRUE)
```

ROC curves show the relation between the two rates being plotted but not their relationships with specific cutoff values.  The latter may be helpful for the selection of a cutoff value to apply in practice.  Accordingly, separate plots of each rate versus the range of possible cutoff values are available with the `type = "cutoffs"` option.

```{r using_analyses_roc_cutoffs}
plot(roc, type = "cutoffs")
```

Area under resulting ROC curves can be computed as an overall measure of model predictive performance and interpreted as the probability that a randomly selected event case will have a higher predicted value than a randomly selected non-event case.

```{r using_analyses_roc_auc}
auc(roc)
```


### Precision Recall

In general, any two binary response metrics may be specified for the construction of a performance curve.  Precision recall curves are another example [@davis:2006:RPR].

```{r using_analyses_pr}
## Precision recall curves
pr <- performance_curve(res_probs, metrics = c(precision, recall))
plot(pr)
```

```{r using_analyses_pr_auc}
auc(pr)
```


### Lift

Lift curves depict the rate at which observed binary responses are identifiable from (resampled) predicted response probabilities.  In particular, they plot the true positive findings (sensitivity) against the positive test rates for all possible classification probability cutoffs.  Accordingly, a lift curve can be interpreted as the rate at which positive responses are found as a function of the positive test rate among cases.

```{r using_analyses_lift}
## Lift curves
lf <- lift(res_probs)
plot(lf, find = 0.75)
```


# Modeling Strategies


## Model Tuning

Many of the modeling functions have arguments, or tuning parameters, that control aspects of their model fitting algorithms.  For example, `GBMModel` parameters `n.trees` and `interaction.depth` control the number of decision trees to fit and the maximum depth of variable interactions.  The `tune` function performs model fitting over a grid of parameter values and returns the model with the most optimal values.  Optimality is determined based on the first performance metric of the `metrics` argument of `tune` if supplied or the first default metric of the `performance` function otherwise.  Argument `grid` additionally controls the construction of grid values and can be a single numeric value giving the grid length in each parameter dimension.  As shown in the output below, tuning results include the tuning parameter grid values, the names of models fit to each, all calculated metrics, the final model selected, the metric upon which its selection was based, and its tuning parameters.

```{r using_strategies_tune}
## Tune over automatic grid of model parameters
(tuned_model <- tune(surv_fo, data = surv_train, model = GBMModel,
                     grid = 3,
                     control = surv_means_control,
                     metrics = c("CIndex" = cindex, "RMSE" = rmse)))
```

Grid values may also be a call to `Grid` with the grid `length` and number of grid points to sample at `random` or a user-specified data frame of grid points.

```{r using_strategies_tune_grid, eval=FALSE}
## Tune over randomly sampled grid points
tune(surv_fo, data = surv_train, model = GBMModel,
     grid = Grid(length = 100, random = 10),
     control = surv_means_control)

## Tune over user-specified grid points
tune(surv_fo, data = surv_train, model = GBMModel,
     grid = expand.grid(n.trees = c(25, 50, 100),
                        interaction.depth = 1:3),
     control = surv_means_control)
```

Statistics summarizing the resampled performance metrics across all tuning parameter combinations can be obtained with the `summary` function.
  
```{r using_strategies_tune_summary}
summary(tuned_model)
```

Line plots of tuning results display the resampled metric means, or another statistic specified with the `stat` argument, versus the first tuning parameter values and with lines grouped according to the remaining parameters, if any.

```{r using_strategies_tune_plot}
plot(tuned_model, type = "line")
```

The return value of `tune` is a model object with the optimal tuning parameters and not a model fit object.  The returned model can be fit subsequently to a set of data with the `fit` function.  Having the final model allows for fitting to the same dataset on which tuning was performed or to datasets from other case populations.

```{r using_strategies_tune_fit}
## Fit the tuned model
surv_fit <- fit(surv_fo, data = surv_train, model = tuned_model)
```

As illustrated in the examples above, the `tune` function allows for comparisons of all models fit to the parameter grid points and returns the optimal model for subsequent fitting.  For applications in which fitting, performance assessment, and prediction with the optimal model are of primary interest, the `TunedModel` function is provided as a model interface to the `tune` function.  The model interface has the advantage of enabling predictive performance estimation, with the `resample` function, that accounts for the full tuning grid process.

```{r using_strategies_tune_model}
## Model interface for grid tuning
tuned_model <- TunedModel(GBMModel, grid = 3, control = surv_means_control,
                          metrics = c("CIndex" = cindex, "RMSE" = rmse))
surv_fit <- fit(surv_fo, data = surv_train, model = tuned_model)
(vi <- varimp(surv_fit))

## Predictive performance of the tuning process
res_tuned <- resample(surv_fo, data = surv_train, model = tuned_model,
                      control = surv_means_control)
summary(performance(res_tuned))
```


## Model Selection

Model selection can be conducted with the `tune` or `SelectedModel` function to automatically choose from any combination of models and model parameters.  Selection has as a special case the just-discussed tuning of a single model over a grid of parameter values.  Combinations of model functions, function names, or function calls can be supplied as lists to the `models` argument of `tune` or as arguments to `SelectedModel` in order to define sets of candidate models from which to select.  An `expand.model` helper function is additionally available to expand a model over a grid of tuning parameters for inclusion in the candidate set if so desired.  In this general form of model selection, the `grid` argument in the tuning functions is not used.

```{r using_strategies_selection, results="hide"}
## Model interface for model selection
selected_model <- SelectedModel(
  expand.model(GBMModel, n.trees = c(50, 100), interaction.depth = 1:2),
  GLMNetModel(lambda = 0.01),
  CoxModel,
  SurvRegModel
)

## Fit the selected model
fit(surv_fo, data = surv_train, model = selected_model)
```

Selection may also be performed over candidate sets that include tuned models.  For instance, the `SelectedModel` function is applicable to sets containing different classes of models each individually tuned over a grid of parameters.

```{r using_strategies_selection_tuned, results="hide"}
## Model interface for selection among tuned models
selected_tuned_model <- SelectedModel(
  TunedModel(GBMModel, control = surv_means_control),
  TunedModel(GLMNetModel, control = surv_means_control),
  TunedModel(CoxModel, control = surv_means_control)
)

## Fit the selected tuned model
fit(surv_fo, data = surv_train, model = selected_tuned_model)
```


## Ensemble Models

Ensemble modelling methods combine $m = 1, \ldots, M$ base learning models as a strategy to improve predictive performance.  Two methods implemented in `Machineshop` are *stacked regression* [@breiman:1996:SR] and *super learners* [@vanderLann:2007:SL].  Stacked regression fits a linear combination of predictions from specified base learners to produce a prediction function of the form
$$
\hat{f}(x) = \sum_{m=1}^M \hat{w}_m \hat{f}_m(x).
$$
Stacking weights $w$ are estimated by (constrained) least squares regression of case responses $y_i$ on predictions $\hat{f}^{-\kappa(i)}(x_i)$ from learners fit to data subsamples $-\kappa(i)$ not containing the corresponding cases. In particular, they are obtained as the solution
$$
\hat{w} = \underset{w}{\operatorname{argmin}} \sum_{i=1}^{N}\left(y_i - \sum_{m=1}^{M} w_m \hat{f}^{-\kappa(i)}(x_i) \right)^2
$$
subject to the constraints that all $w_m \ge 0$ and $\sum_m w_m = 1$.  K-fold cross-validation is the default subsampling method employed in the estimation, with the other resampling methods provided by the package available as other options.  Survival outcomes are handled with a modified version of the stacked regression algorithm in which

* minimization of least squares is replaced by maximization of Harrell's concordance index [-@harrell:1982:EYM] to accommodate censoring, and
* prediction can only be performed on the same response type used for the model fit; i.e., either survival means or survival probabilities at set follow-up times.

Super learners are a generalization of stacked regression that fit a specified model, such as `GBMModel`, to case responses $y_i$, base learner predictions $\hat{f}^{-\kappa(i)}(x_i)$, and optionally also to the original predictor variables $x_i$.  Given below are examples of a stacked regression and super learner each fit with gradient boosted, random forest, and Cox regression base learners.  A separate gradient boosted model is used as the super learner in the latter.

```{r using_strategies_ensembles}
## Stacked regression
stackedmodel <- StackedModel(GLMBoostModel, CForestModel, CoxModel)
res_stacked <- resample(surv_fo, data = surv_train, model = stackedmodel)
summary(res_stacked)

## Super learner
supermodel <- SuperModel(GLMBoostModel, CForestModel, CoxModel,
                         model = GBMModel)
res_super <- resample(surv_fo, data = surv_train, model = supermodel)
summary(res_super)
```


# Package Extensions

Custom models and metrics can be defined with `MLModel` and `MLMetric` for use with the model fitting, prediction, and performance assessment tools provided by the package.


## Custom Models

The `MLModel` function creates a model object that can be used with the previously described fitting functions.  It take the following arguments.

`name`
  : Character name of the object to which the model is assigned.
  
`label`
  : Optional character descriptor for the model (default: `name`).
  
`packages`
  : Character vector of source packages required to use the model.
  
`response_types`
  : Character vector of response variable types to which the model can be fit.  Supported types are `"binary"`, `"factor"`, `"matrix"`, `"numeric"`, `"ordered"`, and `"Surv"`.
  
`fit`
  : Model fitting function whose arguments are a `formula`, a `ModelFrame` named `data`, case `weights`, and an ellipsis.  Argument `data` may be converted to a data frame with the `as.data.frame` function as commonly needed.  The fit function should return the object resulting from the model fit.
  
`predict`
  : Prediction function whose arguments are the `object` returned by `fit`, a `ModelFrame` named
`newdata` of predictor variables, optional vector of `times` at which to predict survival, and an ellipsis.  Argument `data` may be converted to a data frame with the `as.data.frame` function as needed.  Values returned by the function should be formatted according to the response variable types below.

* `factor` : Vector or column matrix of probabilities for the second level of binary factors or a matrix whose columns contain the probabilities for factors with more than two levels.
* `matrix` : Matrix of predicted responses.
* `numeric` : Vector or column matrix of predicted responses.
* `Surv` : Matrix whose columns contain survival probabilities at `times` if supplied or a vector of predicted survival means otherwise.
  
`varimp`
  : Variable importance function whose arguments are the `object` returned by `fit`, optional arguments passed from calls to `varimp`, and an ellipsis.  The function should return a vector of importance values named after the predictor variables or a matrix or data frame whose rows are named after the predictors.

```{r using_extensions_mlmodel}
## Logistic regression model extension
LogisticModel <- MLModel(
  name = "LogisticModel",
  label = "Logistic Model",
  response_types = "binary",
  fit = function(formula, data, weights, ...) {
    glm(formula, data = as.data.frame(data), weights = weights,
        family = binomial, ...)
  },
  predict = function(object, newdata, ...) {
    predict(object, newdata = as.data.frame(newdata), type = "response")
  },
  varimp = function(object, ...) {
    pchisq(coef(object)^2 / diag(vcov(object)), 1)
  }
)
```


## Custom Metrics

The `MLMetric` function creates a metric object that can be used as previously described for the model performance metrics.  Its first argument is a function to compute the metric, defined to accept `observed` and `predicted` as the first two arguments and with an ellipsis to accommodate others.  Its remaining arguments are as follows.

`name`
  :  Character name of the object to which the metric is assigned.
  
`label`
  : Optional character descriptor for the metric (default: `name`).
  
`maximize`
  : Logical indicating whether higher values of the metric correspond to better predictive performance.

```{r using_extensions_mlmetric}
## F2 score metric extension
f2_score <- MLMetric(
  function(observed, predicted, ...) {
    f_score(observed, predicted, beta = 2, ...)
  },
  name = "f2_score",
  label = "F2 Score",
  maximize = TRUE
)
```


## Usage

Once created, model and metric extensions can be used with the package-supplied fitting and performance functions.

```{r using_extensions_usage}
## Logistic regression analysis
library(MASS)
res <- resample(type ~ ., data = Pima.tr, model = LogisticModel)
summary(performance(res, metric = f2_score))
```


# Model Constructor Functions

```{r reference_models, echo = FALSE}
library(MachineShop)

info <- modelinfo()
types <- c("binary" = "b", "factor" = "f", "matrix" = "m", "numeric" = "n",
           "ordered" = "o", "Surv" = "S")
x <- lapply(names(info), function(modelname) {
  c(modelname, ifelse(names(types) %in% info[[modelname]]$response_types, types, NA))
})
df <- as.data.frame(do.call(rbind, x), stringsAsFactors = FALSE)
names(df) <- c("Function", names(types))

toString2 <- function(x) toString(na.omit(x))
df_classes <- data.frame(
  Function = rdoc_url(df$Function),
  Label = sapply(info, getElement, name = "label"),
  Categorical = apply(df[c("binary", "factor", "ordered")], 1, toString2),
  Continuous = apply(df[c("matrix", "numeric")], 1, toString2),
  Survival = apply(df["Surv"], 1, toString2)
)
names(df_classes)[3:5] <- paste0(names(df_classes)[3:5], footnote_marker_number(1:3))

kable(df_classes,
      caption = "Package-supplied model constructor functions and supported response variable types.",
      align = c("l", "l", "c", "c", "c"), row.names = FALSE,
      escape = FALSE) %>%
  kable_styling(c("striped", "condensed"), full_width = FALSE, position = "center") %>%
  add_header_above(c(" " = 1, " " = 1, "Response Variable Types" = 3)) %>%
  footnote(number = c("b = binary factor, f = factor, o = ordered factor",
                      "m = matrix, n = numeric",
                      "S = Surv"))
```


# Metric Functions

```{r reference_metrics, table_metrics, echo=FALSE}
library(MachineShop)

f <- function(x) {
  types <- x$response_types
  
  is_type <- function(observed, predicted) {
    any(types$observed == observed & types$predicted == predicted)
  }
  
  categorical <- if (is_type("factor", "matrix")) "f" else
    if (is_type("factor", "numeric")) "b" else NULL
  if (is_type("ordered", "ordered")) categorical <- c(categorical, "o")

  continuous <- NULL
  if (is_type("matrix", "matrix")) continuous <- "m"
  if (is_type("numeric", "numeric")) continuous <- c(continuous, "n")
  
  survival <- NULL
  if (any(mapply(is_type, "Surv", c("numeric", "SurvEvents", "SurvProbs")))) {
    survival <- "S"
  }

  data.frame(
    Label = x$label,
    Categorical = toString(categorical),
    Continuous = toString(continuous),
    Survival = toString(survival)
  )
}

info <- metricinfo()
df <- cbind("Function" = rdoc_url(names(info)),
            do.call(rbind, lapply(info, f)))
names(df)[3:5] <- paste0(names(df)[3:5], footnote_marker_number(1:3))

kable(df, caption = "Package-supplied performance metric functions and supported response variable types.",
      align = c("l", "l", "c", "c", "c"), row.names = FALSE,
      escape = FALSE) %>%
  kable_styling(c("striped", "condensed"), full_width = FALSE, position = "center") %>%
  add_header_above(c(" " = 1, " " = 1, "Response Variable Types" = 3)) %>%
  footnote(number = c("b = binary factor, f = factor, o = ordered factor",
                      "m = matrix, n = numeric",
                      "S = Surv"))
```


# References
