####################################################################################################
# Functions for generating and plotting ROC curves
####################################################################################################
library(pROC)
#' Generate roc_obj object that is handy for summarizing discrimination statistics and for passing
#' to other functions for presentation
#'
#' @param x Numeric or ordinal vector of predictors
#' @param y Factor, numeric, or character vector of responses
#' @param units Character specifying units for formatting Youden's threshold
#' @param ci Logical indicating whether bootstrapped confidence intervals should be calculated
#' @param rnd_thresh Rounding precision (in decimal places) for Youden's threshold (default is 0)
#' @param rnd_discr Rounding precision for discrimination statistics (sensitivity and specificity; default is 2)
#' @return A roc_obj object
roc_obj <- function(x, y, xlab="", units="", ci=FALSE, rnd_thresh=0, rnd_discr=2) {
if(xlab!="") {
xlab <- paste0(xlab, "=")
}
robj <- roc(predictor=x, response=y)
roc_tbl <- tibble(sens = robj$sens,
spec = robj$spec,
thresh = robj$thresholds,
sum = robj$specificities + (robj$sensitivities)) %>%
arrange(sens)
youden <- roc_tbl %>% arrange(desc(sum)) %>% slice(1)
if(ci) {
sens_ci = ci.se(roc=robj, specificities=youden$spec)[c(1,3)]
spec_ci = ci.sp(roc=robj, sensitivities=youden$sens)[c(1,3)]
youden %<>% mutate(se_lb = sens_ci[1],
se_ub = sens_ci[2],
sp_lb = spec_ci[1],
sp_ub = spec_ci[2])
sens_txt <- with(youden, paste0(rnd(sens, rnd_discr), " (95\\% CI: ", rnd(se_lb, rnd_discr), ", ", rnd(se_ub, rnd_discr), ")"))
spec_txt <- with(youden, paste0(rnd(spec, rnd_discr), " (95\\% CI: ", rnd(sp_lb, rnd_discr), ", ", rnd(sp_ub, rnd_discr), ")"))
} else {
sens_txt <- spec_txt <- NA
}
youden_txt <- c(paste0(xlab, rnd(youden$thresh, rnd_thresh), units),
paste0("(Spec=", rnd(youden$spec, rnd_discr), ", Sens=", rnd(youden$sens, rnd_discr), ")"))
return(list(obj = robj,
tbl = roc_tbl,
youden = list(tbl = youden,
txt = youden_txt,
sens_txt = sens_txt,
spec_txt = spec_txt)))
}
#' Plot a ROC curve based on a roc_obj object
#'
#' @param roc_obj Object generated by roc_obj()
#' @param thresh_pos Character vector specifying whether the plot should place Youden's criterion above or below its point on the curve. Valid options include "top" or "bottom". Defaults to "top."
#' @param thresh_size Numeric specifying the font size for Youden's criterion. Defaults to 7.
#' @return A ggplot object
plot_roc <- function(roc_obj, thresh_pos="top", thresh_size=7) {
offset <- switch(thresh_pos,
top = c(0.08, 0.03),
bottom = c(-0.04, -0.10))
ggplot(roc_obj$tbl, aes(x=(1-spec), y=sens)) +
geom_line(size=1.5) +
geom_point(x = (1-roc_obj$youden$tbl$spec),
y = roc_obj$youden$tbl$sens,
size=4, color="darkred") +
labs(x="1 - specificity", y="Sensitivity") +
annotate("text",
x = (1-roc_obj$youden$tbl$spec) + 0.01,
y = roc_obj$youden$tbl$sens + offset[1],
label = roc_obj$youden$txt[1],
hjust="left",
size=thresh_size) +
annotate("text",
x = (1-roc_obj$youden$tbl$spec) + 0.01,
y = roc_obj$youden$tbl$sens + offset[2],
label = roc_obj$youden$txt[2],
hjust="left",
size=thresh_size*0.7) +
theme(panel.background = element_rect(fill = "ivory"),
panel.grid.major = element_line(color = "ivory2"),
axis.title = element_text(size=15),
axis.text = element_text(size=10))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.