---
title: "cutpointr benchmarks"
author: "Christian Thiele"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{cutpointr benchmarks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(fig.width = 6, fig.height = 5, fig.align = "center")
options(rmarkdown.html_vignette.check_title = FALSE)
load("vignettedata/vignettedata.Rdata")
```


To offer a comparison to established solutions,
**cutpointr** will be benchmarked against `optimal.cutpoints` 
from the **OptimalCutpoints** package, **ThresholdROC** and custom functions based on
the **ROCR** and **pROC** packages. By generating data of different sizes 
the benchmarks will offer a comparison of the scalability of the different 
solutions. 

Using `prediction` and `performance` from the **ROCR** package and `roc` from the
**pROC** package, we can write functions for computing the cutpoint that maximizes the sum of sensitivity and
specificity. **pROC** has a built-in function to optimize a few metrics:

```{r, eval = FALSE}
# Return cutpoint that maximizes the sum of sensitivity and specificiy
# ROCR package
rocr_sensspec <- function(x, class) {
    pred <- ROCR::prediction(x, class)
    perf <- ROCR::performance(pred, "sens", "spec")
    sens <- slot(perf, "y.values")[[1]]
    spec <- slot(perf, "x.values")[[1]]
    cut <- slot(perf, "alpha.values")[[1]]
    cut[which.max(sens + spec)]
}

# pROC package
proc_sensspec <- function(x, class) {
    r <- pROC::roc(class, x, algorithm = 2, levels = c(0, 1), direction = "<")
    pROC::coords(r, "best", ret="threshold", transpose = FALSE)[1]
}
```

The benchmarking will be carried out using the **microbenchmark** package and randomly
generated data. The values of the `x` predictor variable are drawn from a normal distribution
which leads to a lot more unique values than were encountered before in the 
`suicide` data. Accordingly, the search for an optimal cutpoint is much more 
demanding, if all possible cutpoints are evaluated.

Benchmarks are run for sample sizes of 100, 1000, 1e4, 1e5, 1e6, and 1e7.
For low sample sizes **cutpointr** is slower than the other
solutions. While this should be of low practical importance, **cutpointr** scales
more favorably with increasing sample size. The speed disadvantage in small
samples that leads to the lower limit of around 25ms is mainly due to the nesting
of the original data and the results that makes the compact output of `cutpointr`
possible. This observation is emphasized by the fact that `cutpointr::roc` is 
quite fast also in small samples. For sample sizes > 1e5 **cutpointr**
is a little faster than the function based on **ROCR** and **pROC**. Both of these 
solutions are generally faster than **OptimalCutpoints** and **ThresholdROC** with the exception of
small samples. **OptimalCutpoints** and **ThresholdROC** had to be excluded from benchmarks with 
more than 1e4 observations due to high memory requirements and/or excessive run times, rendering
the use of these packages in larger samples impractical.


```{r, eval = FALSE, echo = FALSE}
library(OptimalCutpoints)
library(ThresholdROC)
library(dplyr)
n <- 100
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
x_pos <- dat$x[dat$y == 1]
x_neg <- dat$x[dat$y == 0]
bench_100 <- microbenchmark::microbenchmark(
    cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
              direction = ">=", metric = youden, break_ties = mean),
    rocr_sensspec(dat$x, dat$y),
    proc_sensspec(dat$x, dat$y),
    optimal.cutpoints(X = "x", status = "y", tag.healthy = 0, methods = "Youden",
                      data = dat),
    thres2(k1 = x_neg, k2 = x_pos, rho = 0.5,
           method = "empirical", ci = FALSE),
    times = 100, unit = "ms"
)

n <- 1000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
x_pos <- dat$x[dat$y == 1]
x_neg <- dat$x[dat$y == 0]
bench_1000 <- microbenchmark::microbenchmark(
    cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
              direction = ">=", metric = youden, break_ties = mean),
    rocr_sensspec(dat$x, dat$y),
    proc_sensspec(dat$x, dat$y),
    optimal.cutpoints(X = "x", status = "y", tag.healthy = 0, methods = "Youden",
                      data = dat),
    thres2(k1 = x_neg, k2 = x_pos, rho = 0.5,
           method = "empirical", ci = FALSE),
    times = 100, unit = "ms"
)

n <- 10000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
x_pos <- dat$x[dat$y == 1]
x_neg <- dat$x[dat$y == 0]
bench_10000 <- microbenchmark::microbenchmark(
    cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
              direction = ">=", metric = youden, break_ties = mean, silent = TRUE),
    rocr_sensspec(dat$x, dat$y),
    optimal.cutpoints(X = "x", status = "y", tag.healthy = 0, methods = "Youden",
                      data = dat),
    proc_sensspec(dat$x, dat$y),
    thres2(k1 = x_neg, k2 = x_pos, rho = 0.5,
           method = "empirical", ci = FALSE),
    times = 100
)

n <- 1e5
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1e5 <- microbenchmark::microbenchmark(
    cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
              direction = ">=", metric = youden, break_ties = mean),
    rocr_sensspec(dat$x, dat$y),
    proc_sensspec(dat$x, dat$y),
    times = 100, unit = "ms"
)

n <- 1e6
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1e6 <- microbenchmark::microbenchmark(
    cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
              direction = ">=", metric = youden, break_ties = mean),
    rocr_sensspec(dat$x, dat$y),
    proc_sensspec(dat$x, dat$y),
    times = 30, unit = "ms"
)

n <- 1e7
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1e7 <- microbenchmark::microbenchmark(
    cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
              direction = ">=", metric = youden, break_ties = mean),
    rocr_sensspec(dat$x, dat$y),
    proc_sensspec(dat$x, dat$y),
    times = 30, unit = "ms"
)

results <- rbind(
    data.frame(time = summary(bench_100)$median,
               Solution = summary(bench_100)$expr,
               n = 100),
    data.frame(time = summary(bench_1000)$median,
               Solution = summary(bench_1000)$expr,
               n = 1000),
    data.frame(time = summary(bench_10000)$median,
               Solution = summary(bench_10000)$expr,
               n = 10000),
    data.frame(time = summary(bench_1e5)$median,
               Solution = summary(bench_1e5)$expr,
               n = 1e5),
    data.frame(time = summary(bench_1e6)$median,
               Solution = summary(bench_1e6)$expr,
               n = 1e6),
    data.frame(time = summary(bench_1e7)$median,
               Solution = summary(bench_1e7)$expr,
               n = 1e7)
)
results$Solution <- as.character(results$Solution)
results$Solution[grep(pattern = "cutpointr", x = results$Solution)] <- "cutpointr"
results$Solution[grep(pattern = "rocr", x = results$Solution)] <- "ROCR"
results$Solution[grep(pattern = "optimal", x = results$Solution)] <- "OptimalCutpoints"
results$Solution[grep(pattern = "proc", x = results$Solution)] <- "pROC"
results$Solution[grep(pattern = "thres", x = results$Solution)] <- "ThresholdROC"

results$task <- "Cutpoint Estimation"
```

```{r, echo = FALSE}
# These are the original results on our system
# dput(results)
results <- structure(list(time = c(4.5018015, 1.812802, 0.662101, 2.2887015, 
1.194301, 4.839401, 2.1764015, 0.981001, 45.0568005, 36.2398515, 
8.5662515, 5.667101, 2538.612001, 4.031701, 2503.8012505, 45.384501, 
43.118751, 37.150151, 465.003201, 607.023851, 583.0950005, 5467.332801, 
7850.2587, 7339.356101), Solution = c("cutpointr", "ROCR", "pROC", 
"OptimalCutpoints", "ThresholdROC", "cutpointr", "ROCR", "pROC", 
"OptimalCutpoints", "ThresholdROC", "cutpointr", "ROCR", "OptimalCutpoints", 
"pROC", "ThresholdROC", "cutpointr", "ROCR", "pROC", "cutpointr", 
"ROCR", "pROC", "cutpointr", "ROCR", "pROC"), n = c(100, 100, 
100, 100, 100, 1000, 1000, 1000, 1000, 1000, 10000, 10000, 10000, 
10000, 10000, 1e+05, 1e+05, 1e+05, 1e+06, 1e+06, 1e+06, 1e+07, 
1e+07, 1e+07), task = c("Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation", 
"Cutpoint Estimation")), row.names = c(NA, -24L), class = "data.frame")
```


```{r, eval = FALSE}
# ROCR package
rocr_roc <- function(x, class) {
    pred <- ROCR::prediction(x, class)
    perf <- ROCR::performance(pred, "sens", "spec")
    return(NULL)
}

# pROC package
proc_roc <- function(x, class) {
    r <- pROC::roc(class, x, algorithm = 2, levels = c(0, 1), direction = "<")
    return(NULL)
}
```



```{r, eval = FALSE, echo = FALSE}
n <- 100
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_100 <- microbenchmark::microbenchmark(
    cutpointr::roc(dat, "x", "y", pos_class = 1,
                   neg_class = 0, direction = ">="),
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$x, dat$y),
    times = 100, unit = "ms"
)
n <- 1000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1000 <- microbenchmark::microbenchmark(
    cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
                   direction = ">="),
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$x, dat$y),
    times = 100, unit = "ms"
)
n <- 10000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_10000 <- microbenchmark::microbenchmark(
    cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
                   direction = ">="),
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$x, dat$y),
    times = 100, unit = "ms"
)
n <- 1e5
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1e5 <- microbenchmark::microbenchmark(
    cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
                   direction = ">="),
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$x, dat$y),
    times = 100, unit = "ms"
)
n <- 1e6
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1e6 <- microbenchmark::microbenchmark(
    cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
                   direction = ">="),
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$x, dat$y),
    times = 30, unit = "ms"
)
n <- 1e7
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
bench_1e7 <- microbenchmark::microbenchmark(
    cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
                   direction = ">="),
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$x, dat$y),
    times = 30, unit = "ms"
)

results_roc <- rbind(
    data.frame(time = summary(bench_100)$median,
               Solution = summary(bench_100)$expr,
               n = 100),
    data.frame(time = summary(bench_1000)$median,
               Solution = summary(bench_1000)$expr,
               n = 1000),
    data.frame(time = summary(bench_10000)$median,
               Solution = summary(bench_10000)$expr,
               n = 10000),
    data.frame(time = summary(bench_1e5)$median,
               Solution = summary(bench_1e5)$expr,
               n = 1e5),
    data.frame(time = summary(bench_1e6)$median,
               Solution = summary(bench_1e6)$expr,
               n = 1e6),
    data.frame(time = summary(bench_1e7)$median,
               Solution = summary(bench_1e7)$expr,
               n = 1e7)
)
results_roc$Solution <- as.character(results_roc$Solution)
results_roc$Solution[grep(pattern = "cutpointr", x = results_roc$Solution)] <- "cutpointr"
results_roc$Solution[grep(pattern = "rocr", x = results_roc$Solution)] <- "ROCR"
results_roc$Solution[grep(pattern = "proc", x = results_roc$Solution)] <- "pROC"
results_roc$task <- "ROC curve calculation"
```

```{r, echo = FALSE}
# Our results
results_roc <- structure(list(time = c(0.7973505, 1.732651, 0.447701, 0.859301, 
2.0358515, 0.694802, 1.878151, 5.662151, 3.6580505, 11.099251, 
42.8208515, 35.3293005, 159.8100505, 612.471901, 610.4337005, 
2032.693551, 7806.3854515, 7081.897251), Solution = c("cutpointr", 
"ROCR", "pROC", "cutpointr", "ROCR", "pROC", "cutpointr", "ROCR", 
"pROC", "cutpointr", "ROCR", "pROC", "cutpointr", "ROCR", "pROC", 
"cutpointr", "ROCR", "pROC"), n = c(100, 100, 100, 1000, 1000, 
1000, 10000, 10000, 10000, 1e+05, 1e+05, 1e+05, 1e+06, 1e+06, 
1e+06, 1e+07, 1e+07, 1e+07), task = c("ROC curve calculation", 
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation", 
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation", 
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation", 
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation", 
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation", 
"ROC curve calculation", "ROC curve calculation")), row.names = c(NA, 
-18L), class = "data.frame")
```

```{r, echo = FALSE}
library(ggplot2)
results_all <- dplyr::bind_rows(results, results_roc)

ggplot(results_all, aes(x = n, y = time, col = Solution, shape = Solution)) +
  geom_point(size = 3) + geom_line() +
  scale_y_log10(breaks = c(0.5, 1, 2, 3, 5, 10, 25, 100, 250, 1000, 5000, 1e4, 15000)) +
  scale_x_log10(breaks = c(100, 1000, 1e4, 1e5, 1e6, 1e7)) +
  ylab("Median Time (milliseconds, log scale)") + xlab("Sample Size (log scale)") +
  theme_bw() +
  theme(legend.position = "bottom", 
        legend.key.width = unit(0.8, "cm"), 
        panel.spacing = unit(1, "lines")) +
  facet_grid(~task)
```


```{r, echo = FALSE}
library(tidyr)
res_table <- tidyr::spread(results_all, Solution, time)
res_table <- dplyr::arrange(res_table, task)
knitr::kable(res_table)
```