README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# modelrecon

The goal of modelrecon is to apply thresholds to predicted probabilities and calculate net benefit in the presence of resource constraints.

## Installation

You can install the released version of modelrecon from GitHub with:

``` r
remotes::install_github('ML4LHS/modelrecon')
```

## Get started

Let's load the package and generate and example dataset containing the probability of an adverse outcome and whether or not that outcome was experienced (`TRUE`) or not experienced (`FALSE`).

```{r example}
library(modelrecon)
library(dplyr)
library(tidyr)
library(ggplot2)

example_data = data.frame(probability = c(0.8, 0.7, 0.6, 0.5, 0.3, 0.2, 0.1, 0.1, 0.05, 0.01),
                          outcome = c(T, T, F, T, T, F, F, F, T, F))
```


## What does our example dataset look like?

```{r}
example_data
```


## Let's apply a threshold of 0.2

This means we will call all predictions with a probability >= 0.2 as `TRUE`.

```{r}
example_data %>% 
  apply_threshold(0.2)
```


## Let's calculate a net benefit with a threshold of 0.2


```{r}
example_data %>%
  apply_threshold(0.2) %>%
  calculate_net_benefit()
```

## What is going on behind the scenes?

Behind the scenes, the `calculate_net_benefit()` function is calculating the number of true and false positives, and then using that along with the previously applied threshold to calculate the net benefit.

### How did `calculate_net_benefit()` know about the threshold?

This information is captured in the `thresholds` attribute.

```{r}
example_data %>%
  apply_threshold(0.2) %>%
  attributes() %>% 
  .$thresholds
```
### Want more information about the number of true and false positives?

Set the `verbose` argument of `calculate_net_benefit()` to `TRUE`. This will print, *not* return, a data frame with the information it used to calculate the net benefit. The value returned is still the net benefit.

```{r}
example_data %>%
  apply_threshold(0.2) %>%
  calculate_net_benefit(verbose = TRUE)
```

## What happens when you apply an absolute constraint?

Two of the five predicted `TRUE` values are converted to `FALSE` because only the first 3 `TRUE` values (those with the highest predicted probability) are able to be acted upon.

```{r}
example_data %>%
  apply_threshold(0.2) %>%
  apply_constraint(3)
```


## Calculate a realized net benefit with a threshold of 0.2 and an capacity of 3

This is an example of an *absolute* constraint.

```{r}
example_data %>%
  apply_threshold(0.2) %>%
  apply_constraint(3) %>%
  calculate_net_benefit(verbose = TRUE)

```


## Calculate a realized net benefit with an absolute threshold of 0.2 and capacity of 3, and *then* a relative constraint of 0.5:

```{r}
example_data %>%
  apply_threshold(0.2) %>%
  apply_constraint(3) %>%
  apply_threshold(0.5) %>%
  calculate_net_benefit(verbose = TRUE)
```

## The default assumption when we set a threshold without a subsequent constraint is that the capacity is infinite.

You can also explicitly note the infinite capacity, which will be applied only to the immediate prior threshold.

```{r}
example_data %>%
  apply_threshold(0.2) %>%
  apply_constraint(3) %>%
  apply_threshold(0.5) %>%
  apply_constraint(Inf) %>% 
  calculate_net_benefit()
```


Using this mechanism, you can construct multiple layers of absolute and relative constraints as the piped functions retain metadata about prior constraints and thus know that each constraint applies to only the prior threshold.

You *cannot* apply a threshold that is *lower* than a prior threshold because it would make no sense to apply a permissive criterion *before* a more restrictive one.

## Setting a new threshold that is lower than the prior one will generate an error.

```{r error = TRUE}
example_data %>%
  apply_threshold(0.5) %>%
  apply_constraint(3) %>%
  apply_threshold(0.2) %>%
  calculate_net_benefit()
```


## Setting a new threshold that is the same as a prior one will generate a warning.

In a future version, this may be upgraded to an error.

```{r error = TRUE}
example_data %>%
  apply_threshold(0.2) %>%
  apply_constraint(3) %>%
  apply_threshold(0.2) %>%
  calculate_net_benefit()
```

# Let's plot a decision curve for an absolute constraint, and an absolute + relative constraint

```{r, warning = FALSE}
plot_data =
  expand_grid(constraint = c(0, 1, 3, 5, 7, Inf),
              threshold = seq(from = 0, to = 1, by = 0.05)) %>%
  group_by(constraint, threshold) %>%
  mutate(net_benefit = example_data %>%
           apply_threshold(threshold) %>%
           apply_constraint(constraint) %>%
           calculate_net_benefit()) %>%
  ungroup()

# Vary absolute constraint and add relative constraint (up to threshold of 0.5)

plot_data_2 =
  expand_grid(constraint = c(0, 1, 3, 5, 7, Inf),
              threshold = seq(from = 0, to = 0.5, by = 0.05)) %>%
  group_by(constraint, threshold) %>%
  mutate(net_benefit = example_data %>%
           apply_threshold(threshold) %>%
           apply_constraint(constraint) %>%
           apply_threshold(pmax(0.5, threshold)) %>%
           calculate_net_benefit()) %>%
  ungroup()

bind_rows(
  plot_data %>% mutate(constraint_type = 'Absolute constraint'),
  plot_data_2 %>% mutate(constraint_type = paste0('Absolute constraint\n',
                                                  'relaxed by relative\n',
                                                  'constraint at threshold\n',
                                                  'of 0.5'))
) %>%
  mutate(constraint = if_else(constraint == Inf, 'Infinity', as.character(constraint))) %>%
  mutate(constraint = as.factor(paste('Capacity =',constraint))) %>%
  filter((constraint == 'Capacity = 3' & threshold == 0.2) |
           (constraint == 'Capacity = Infinity' & threshold == 0.2)) %>% 
  slice(1:3) %>% 
  mutate(text = c('Case study 2', 'Case study 1', 'Case study 3')) %>% 
  mutate(x = c(0.3, 0.4, 0.3), y = c(0.05, 0.4, 0.33)) -> 
  point_data

bind_rows(
  plot_data %>% mutate(constraint_type = 'Absolute constraint'),
  plot_data_2 %>% mutate(constraint_type = paste0('Absolute constraint\n',
                                                  'relaxed by relative\n',
                                                  'constraint at threshold\n',
                                                  'of 0.5'))
  ) %>%
  mutate(constraint = if_else(constraint == Inf, 'Infinity', as.character(constraint))) %>%
  mutate(constraint = as.factor(paste('Capacity =',constraint))) %>%
  ggplot(aes(x = threshold, y = net_benefit,
             linetype = constraint_type)) +
  geom_line() +
  geom_point(data = point_data) +
  geom_text(data = point_data,
            aes(label = text, x = x, y = y), size = 3) +
  facet_wrap(~constraint) +
  coord_cartesian(ylim = c(0, 0.5)) +
  theme_bw() +
  theme(axis.text = element_text(size = 6)) +
  labs(x = 'Threshold probability',
       y = 'Realized net benefit',
       linetype = 'Constraint')

# ggsave('Figure 2.pdf',
#        width = 6.5, height = 4, units = 'in')

```