In xrobin/pROC: Display and Analyze ROC Curves
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "bench/fig-"
)
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,
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.
# 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
}
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"
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)")
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).
# 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
}
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"
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)")
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.
# 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
}
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"
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)")
res_table <- tidyr::spread(best_results, solution, time)
knitr::kable(res_table)
xrobin/pROC documentation built on Nov. 7, 2023, 2:34 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "bench/fig-" )
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, 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.
# 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 }
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"
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)")
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).
# 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 }
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"
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)")
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.
# 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 }
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"
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)")
res_table <- tidyr::spread(best_results, solution, time) knitr::kable(res_table)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.