---
output: 
    github_document:
        toc: true
        toc_depth: 2
---

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

```{r, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "bench/fig-"
)
```

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

# Benchmarks

These benchmarks compare pROC with competing ROC analysis packages in R. They can 
serve as a way to detect performance bottleneck that must be fixed, and 
possible regressions in performance.

The benchmarking are carried out with the **microbenchmark** package and randomly
generated data. The values of the `x` predictor variable are drawn from a normal distribution, resulting in every value being 
essentially unique. Predictor values for positive examples are increased to have a mean of 1, resulting in ROC curves with an AUC of 0.76.

<!-- TODO: bench with rounded or fewer unique values -->

The benchmark code is adapted from the [cutpointr vignette by Christian Thiele](https://github.com/Thie1e/cutpointr/blob/master/vignettes/cutpointr.Rmd),
released under a GPL-3 license.

## Building the ROC curve

This first benchmark looks at the time needed to building the ROC curve only, and getting sensitivities, specificities and thresholds. Only packages allowing turn off the
calculation of the AUC, or not computing it by default, were tested.

```{r, echo = FALSE}
# Simply compute sensitivity, specificity and thresholds

rocr_roc <- function(predictor, response) {
  pred <- ROCR::prediction(predictor, response)
  perf <- ROCR::performance(pred, "sens", "spec")
  se <- slot(perf, "y.values")[[1]]
  sp <- slot(perf, "x.values")[[1]]
  thr <- slot(perf, "alpha.values")[[1]]
}

proc_roc <- function(response, predictor) {
  r <- pROC::roc(response, predictor, algorithm = 2, levels = c(0, 1), direction = "<",
                 auc = FALSE)
  se <- r$sensitivities
  sp <- r$specificities
  thr <- r$thresholds
}
```


```{r, echo = FALSE}
n <- 1000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
roc_bench_1000 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x)
)

n <- 10000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
roc_bench_10000 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    times=50
)

n <- 1e5
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
roc_bench_1e5 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    times = 20
)

n <- 1e6
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
roc_bench_1e6 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    times = 15
)

n <- 1e7
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
roc_bench_1e7 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    times = 10
)

roc_results <- rbind(
    data.frame(time = summary(roc_bench_1000)$median,
               solution = summary(roc_bench_1000)$expr, 
               n = 1000),
    data.frame(time = summary(roc_bench_10000)$median,
               solution = summary(roc_bench_10000)$expr, 
               n = 10000),
    data.frame(time = summary(roc_bench_1e5)$median,
               solution = summary(roc_bench_1e5)$expr, 
               n = 1e5),
    data.frame(time = summary(roc_bench_1e6)$median,
               solution = summary(roc_bench_1e6)$expr, 
               n = 1e6),
    data.frame(time = summary(roc_bench_1e7)$median,
               solution = summary(roc_bench_1e7)$expr, 
               n = 1e7)
)
roc_results$solution <- as.character(roc_results$solution)
roc_results$solution[grep(pattern = "rocr", x = roc_results$solution)] <- "ROCR"
roc_results$solution[grep(pattern = "proc", x = roc_results$solution)] <- "pROC"
```



```{r, echo = FALSE}
ggplot(roc_results, aes(x = n, y = time, col = solution, shape = solution)) +
    geom_point(size = 3) + geom_line() +
    scale_y_log10(breaks = c(3, 5, 10, 25, 100, 250, 1000, 5000, 1e4, 15000)) +
    scale_x_log10(breaks = c(1000, 1e4, 1e5, 1e6, 1e7)) +
    ggtitle("ROC building benchmark results", "n = 1000, 10000, 1e5, 1e6, 1e7") +
    ylab("Median time (milliseconds, log scale)") + xlab("n (log scale)")
```
```{r, echo = FALSE}
res_table <- tidyr::spread(roc_results, solution, time)
knitr::kable(res_table)
```
## AUC

This benchmark tests how long it takes to calculate the ROC curve and the area under the ROC curve (AUC).

```{r, echo = FALSE}
# Calculate the AUC

rocr_auc <- function(predictor, response) {
    pred <- ROCR::prediction(predictor, response)
    perf <- ROCR::performance(pred, measure = "auc")
    perf@y.values[[1]]
}

proc_auc <- function(response, predictor) {
      r <- pROC::roc(response, predictor, algorithm = 2, levels = c(0, 1), direction = "<")
      r$auc
}

prroc_auc <- function(positives, negatives) {
  r <- PRROC::roc.curve(positives, negatives)
  r$auc
}

epi_auc <- function(predictor, response) {
  e <- Epi::ROC(predictor, response, plot=FALSE)
  e$AUC
}
```

```{r, echo = FALSE}
n <- 1000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
negatives <- dat$x[dat$y == 0]
positives <- dat$x[dat$y == 1]
auc_bench_1000 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    prroc_auc(positives, negatives),
    epi_auc(dat$x, dat$y)
)

n <- 10000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
negatives <- dat$x[dat$y == 0]
positives <- dat$x[dat$y == 1]
auc_bench_10000 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    prroc_auc(positives, negatives),
    epi_auc(dat$x, dat$y),
    times=50
)

n <- 1e5
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
negatives <- dat$x[dat$y == 0]
positives <- dat$x[dat$y == 1]
auc_bench_1e5 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    prroc_auc(positives, negatives),
    epi_auc(dat$x, dat$y),
    times = 20
)

n <- 1e6
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
negatives <- dat$x[dat$y == 0]
positives <- dat$x[dat$y == 1]
auc_bench_1e6 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    prroc_auc(positives, negatives),
    epi_auc(dat$x, dat$y),
    times = 15
)

n <- 1e7
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
negatives <- dat$x[dat$y == 0]
positives <- dat$x[dat$y == 1]
auc_bench_1e7 <- microbenchmark::microbenchmark(unit = "ms",
    rocr_roc(dat$x, dat$y),
    proc_roc(dat$y, dat$x),
    prroc_auc(positives, negatives),
    times = 10
)

auc_results <- rbind(
    data.frame(time = summary(auc_bench_1000)$median,
               solution = summary(auc_bench_1000)$expr, 
               n = 1000),
    data.frame(time = summary(auc_bench_10000)$median,
               solution = summary(auc_bench_10000)$expr, 
               n = 10000),
    data.frame(time = summary(auc_bench_1e5)$median,
               solution = summary(auc_bench_1e5)$expr, 
               n = 1e5),
    data.frame(time = summary(auc_bench_1e6)$median,
               solution = summary(auc_bench_1e6)$expr,
               n = 1e6),
    data.frame(time = summary(auc_bench_1e7)$median,
               solution = summary(auc_bench_1e7)$expr,
               n = 1e7)
)
auc_results$solution <- as.character(auc_results$solution)
auc_results$solution[grep(pattern = "epi", x = auc_results$solution)] <- "Epi"
auc_results$solution[grep(pattern = "prroc", x = auc_results$solution)] <- "PRROC"
auc_results$solution[grep(pattern = "rocr", x = auc_results$solution)] <- "ROCR"
auc_results$solution[grep(pattern = "proc", x = auc_results$solution)] <- "pROC"
```


```{r, echo = FALSE}
ggplot(auc_results, aes(x = n, y = time, col = solution, shape = solution)) +
    geom_point(size = 3) + geom_line() +
    scale_y_log10(breaks = c(3, 5, 10, 25, 100, 250, 1000, 5000, 1e4, 15000)) +
    scale_x_log10(breaks = c(1000, 1e4, 1e5, 1e6, 1e7)) +
    ggtitle("ROC building benchmark results", "n = 1000, 10000, 1e5, 1e6, 1e7") +
    ylab("Median time (milliseconds, log scale)") + xlab("n (log scale)")
```
```{r, echo = FALSE}
res_table <- tidyr::spread(auc_results, solution, time)
knitr::kable(res_table)
```

## Best threshold

Benchmarks packages that extract the "best" threshold. At the moment they all use the Youden index. This includes building the ROC curve first.

```{r, echo = FALSE}
# Get the best threshold as a numeric value 

proc_best <- function(response, predictor) {
  r <- pROC::roc(response, predictor, algorithm = 2, levels = c(0, 1), direction = "<")
  pROC::coords(r, "best", ret="threshold", drop=TRUE)
}

cutpointr_best <- function(data, predictor_name, response_name) {
  cu <- cutpointr::cutpointr_(data, predictor_name, response_name, 
                       pos_class = 1, neg_class = 0,
                       direction = ">=", metric = cutpointr::youden, 
                       break_ties = mean)
  cu[,"optimal_cutpoint", drop=TRUE]
}

optimalcutpoints_best <- function(data, predictor_name, response_name) {
  o <- OptimalCutpoints::optimal.cutpoints(predictor_name, response_name, data=data,
                                           tag.healthy = 0, methods = "Youden")
  o$Youden$Global$optimal.cutoff$cutoff
}

thresholdroc_best <- function(negatives, positives) {
  tr <- ThresholdROC::thres2(negatives, positives, rho = 0.5, 
           method = "empirical", ci = FALSE)
  tr$T$thres
}
```




```{r, echo = FALSE}
n <- 100
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
positives <- dat$x[dat$y == 1]
negatives <- dat$x[dat$y == 0]
best_bench_100 <- microbenchmark::microbenchmark(
  proc_best(dat$y, dat$x),
  cutpointr_best(dat, "x", "y"),
  optimalcutpoints_best(dat, "x", "y"),
  thresholdroc_best(negatives, positives),
  unit = "ms"
)

n <- 1000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
positives <- dat$x[dat$y == 1]
negatives <- dat$x[dat$y == 0]
best_bench_1000 <- microbenchmark::microbenchmark(
  proc_best(dat$y, dat$x),
  cutpointr_best(dat, "x", "y"),
  optimalcutpoints_best(dat, "x", "y"),
  thresholdroc_best(negatives, positives),
  unit = "ms"
)

n <- 10000
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
positives <- dat$x[dat$y == 1]
negatives <- dat$x[dat$y == 0]
best_bench_10000 <- microbenchmark::microbenchmark(
  proc_best(dat$y, dat$x),
  cutpointr_best(dat, "x", "y"),
  optimalcutpoints_best(dat, "x", "y"),
  thresholdroc_best(negatives, positives),
  times = 20, unit = "ms"
)

n <- 1e5
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
positives <- dat$x[dat$y == 1]
negatives <- dat$x[dat$y == 0]
best_bench_1e5 <- microbenchmark::microbenchmark(
  proc_best(dat$y, dat$x),
  cutpointr_best(dat, "x", "y"),
  times = 20, unit = "ms"
)

n <- 1e6
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
best_bench_1e6 <- microbenchmark::microbenchmark(
  proc_best(dat$y, dat$x),
  cutpointr_best(dat, "x", "y"),
  times = 10, unit = "ms"
)

n <- 1e7
set.seed(123)
dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
dat$x <- dat$x + dat$y
best_bench_1e7 <- microbenchmark::microbenchmark(
  proc_best(dat$y, dat$x),
  cutpointr_best(dat, "x", "y"),
    times = 10, unit = "ms"
)

best_results <- rbind(
    data.frame(time = summary(best_bench_100)$median,
               solution = summary(best_bench_100)$expr, 
               n = 100),
    data.frame(time = summary(best_bench_1000)$median,
               solution = summary(best_bench_1000)$expr, 
               n = 1000),
    data.frame(time = summary(best_bench_10000)$median,
               solution = summary(best_bench_10000)$expr, 
               n = 10000),
    data.frame(time = summary(best_bench_1e5)$median,
               solution = summary(best_bench_1e5)$expr, 
               n = 1e5),
    data.frame(time = summary(best_bench_1e6)$median,
               solution = summary(best_bench_1e6)$expr, 
               n = 1e6),
    data.frame(time = summary(best_bench_1e7)$median,
               solution = summary(best_bench_1e7)$expr, 
               n = 1e7)
)
best_results$solution <- as.character(best_results$solution)
best_results$solution[grep(pattern = "cutpointr", x = best_results$solution)] <- "cutpointr"
best_results$solution[grep(pattern = "optimalcutpoints", x = best_results$solution)] <- "OptimalCutpoints"
best_results$solution[grep(pattern = "proc", x = best_results$solution)] <- "pROC"
best_results$solution[grep(pattern = "thresholdroc", x = best_results$solution)] <- "ThresholdROC"
```


```{r, echo = FALSE}
ggplot(best_results, aes(x = n, y = time, col = solution, shape = solution)) +
    geom_point(size = 3) + geom_line() +
    scale_y_log10(breaks = c(3, 5, 10, 25, 100, 250, 1000, 5000, 1e4, 15000)) +
    scale_x_log10(breaks = c(100, 1000, 1e4, 1e5, 1e6, 1e7)) +
    ggtitle("Benchmark results", "n = 1000, 10000, 1e5, 1e6, 1e7") +
    ylab("Median time (milliseconds, log scale)") + xlab("n (log scale)")
```

```{r, echo = FALSE}
res_table <- tidyr::spread(best_results, solution, time)
knitr::kable(res_table)
```