---
output:
github_document:
toc: true
toc_depth: 2
---
```{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.
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)
```