Generate ROC Curve Charts for Print and Interactive Use

Michael C Sachs

2023-10-06

Introduction

About ROC Curves

The Receiver Operating Characteristic (ROC) curve is used to assess the accuracy of a continuous measurement for predicting a binary outcome. In medicine, ROC curves have a long history of use for evaluating diagnostic tests in radiology and general diagnostics. ROC curves have also been used for a long time in signal detection theory.

The accuracy of a diagnostic test can be evaluated by considering the two possible types of errors: false positives, and false negatives. For a continuous measurement that we denote as \(M\), convention dictates that a test positive is defined as \(M\) exceeding some fixed threshold \(c\): \(M > c\). In reference to the binary outcome that we denote as \(D\), a good outcome of the test is when the test is positive among an individual who truly has a disease: \(D = 1\). A bad outcome is when the test is positive among an individual who does not have the disease \(D = 0\).

Formally, for a fixed cutoff \(c\), the true positive fraction is the probability of a test positive among the diseased population:

\[ TPF(c) = P\{ M > c | D = 1 \} \]

and the false positive fraction is the probability of a test positive among the healthy population:

\[ FPF(c) = P\{ M > c | D = 0 \} \]

Since the cutoff \(c\) is not usually fixed in advance, we can plot the TPF against the FPF for all possible values of \(c\). This is exactly what the ROC curve is, \(FPF(c)\) on the \(x\) axis and \(TPF(c)\) along the \(y\) axis.

Motivation

In the medical literature, ROC curves are commonly plotted without the cutoff values displayed. Other problems with ROC curve plots are abundant in the medical literature. We aim to solve some of these problems by providing a plotting interface for the ROC curve that comes with sensible defaults. It is easy to create interactive ROC curves for local or web-based use. The next section details the usage of the plotROC package.

Usage

Shiny application

I created a shiny application in order to make the features more accessible to non-R users. A limited subset of the functions of the plotROC package can be performed on an example dataset or on data that users upload to the website. Resulting plots can be saved to the users’ machine as a pdf or as a stand-alone html file. It can be used in any modern web browser with no other dependencies at the website here: https://sachsmc.shinyapps.io/plotROC.

Installation and loading

plotROC can be installed from github or CRAN. It requires a recent version of ggplot2 (>2.0.0).

install.packages("ggplot2")
install.packages("plotROC")
library(plotROC)

Quick start

After installing, the interactive Shiny application can be run locally.

shiny_plotROC()

Command line basic usage

I start by creating an example data set. There are 2 markers, one that is moderately predictive and one that is not as predictive.

## Loading required package: ggplot2
set.seed(2529)
D.ex <- rbinom(200, size = 1, prob = .5)
M1 <- rnorm(200, mean = D.ex, sd = .65)
M2 <- rnorm(200, mean = D.ex, sd = 1.5)

test <- data.frame(D = D.ex, D.str = c("Healthy", "Ill")[D.ex + 1], 
                   M1 = M1, M2 = M2, stringsAsFactors = FALSE)

The Roc Geom

Next I use the ggplot function to define the aesthetics, and the geom_roc function to add an ROC curve layer. The geom_roc function requires the aesthetics d for disease status, and m for marker. The disease status need not be coded as 0/1, but if it is not, stat_roc assumes (with a warning) that the lowest value in sort order signifies disease-free status. stat_roc and geom_roc are linked by default, with the stat doing the underlying computation of the empirical ROC curve, and the geom consisting of the ROC curve layer.

basicplot <- ggplot(test, aes(d = D, m = M1)) + geom_roc()
basicplot

The disease status aesthetic can be specified as a string or factor, but with a warning.

ggplot(test, aes(d = D.str, m = M1)) + geom_roc()
## Warning in verify_d(data$d): D not labeled 0/1, assuming Healthy = 0 and Ill =
## 1!

The geom_roc layer includes the ROC curve line combined with points and labels to display the values of the biomarker at the different cutpoints. It accepts the argument n.cuts to define the number of cutpoints to display along the curve. Labels can be suppressed by using n.cuts = 0 or labels = FALSE. The size of the labels and the number of significant digits can be adjusted with labelsize and labelround, respectively.

ggplot(test, aes(d = D, m = M1)) + geom_roc(n.cuts = 0)

ggplot(test, aes(d = D, m = M1)) + geom_roc(n.cuts = 5, labelsize = 5, labelround = 2)

ggplot(test, aes(d = D, m = M1)) + geom_roc(n.cuts = 50, labels = FALSE)

We provide a function style_roc that can be added to a ggplot that contains an ROC curve layer. This adds a diagonal guideline, sets the axis labels, and adjusts the major and minor grid lines. The direct_label function operates on a ggplot object, adding a direct label to the plot. It attempts to intelligently select an appropriate location for the label, but the location can be adjusted with nudge_x, nudge_y and label.angle. If the labels argument is NULL, it will take the name from the mapped aesthetic.

styledplot <- basicplot + style_roc()
styledplot

basicplot + style_roc(theme = theme_grey, xlab = "1 - Specificity")

direct_label(basicplot) + style_roc()

direct_label(basicplot, labels = "Biomarker", nudge_y = -.1) + style_roc()

Confidence regions and the Rocci Geom

It is common to compute confidence regions for points on the ROC curve using the Clopper and Pearson (1934) exact method. Briefly, exact confidence intervals are calculated for the \(FPF\) and \(TPF\) separately, each at level \(1 - \sqrt{1 - \alpha}\). Based on result 2.4 from Pepe (2003), the cross-product of these intervals yields a \(100 * (1 - \alpha)\) percent rectangular confidence region for the pair.

This is implemented in the stat_rocci and displayed as a geom_rocci layer. These both require the same aesthetics as the ROC geom, d for disease status and m for marker. By default, a set of 3 evenly spaced points along the curve are chosen to display confidence regions. You can select points by passing a vector of values in the range of m to the ci.at argument. By default, the significance level \(\alpha\) is set to 0.05, this can be changed using the sig.level option.

styledplot + geom_rocci()

styledplot + geom_rocci(sig.level = .01)

ggplot(test, aes(d = D, m = M1)) + geom_roc(n.cuts = 0) +
  geom_rocci(ci.at = quantile(M1, c(.1, .4, .5, .6, .9))) + style_roc()

Interactive Plots

Ggplot objects that contain a GeomRoc layer can be used to create an interactive plot and display it in the Rstudio viewer or default web browser by passing it to the plot_interactive_roc, or export_interactive_roc function. The style_roc function is applied by default. Give the function an optional path to an html file as an argument called file to save the interactive plot as a complete web page. By default, any existing Rocci layers are removed and replaced with a dense layer of confidence regions so that the user can click anywhere for a confidence region. This can be suppressed by add.cis = FALSE. Furthermore, the points layer of the Roc geom can be hidden by using the hide.points option.

Hovering over the display shows the cutoff value at the point nearest to the cursor. Clicking makes the cutoff label stick until the next click, and if confidence regions are available, clicks will also display those as grey rectangles. The confidence regions are automatically detected. When the user clicks on the ROC curve, the confidence region for the TPF and FPF is overlaid using a grey rectangle. The label and region stick until the next click.

plot_interactive_roc(basicplot)

An interactive ROC plot can be exported by using the export_interactive_roc function, which returns a character string containing the necessary HTML and JavaScript. The character string can be copy-pasted into an html document, or better yet, incorporated directly into a dynamic document using knitr (knitr homepage).

In a knitr document, it is necessary to use the cat function on the results and use the chunk options results = 'asis' and fig.keep='none' so that the interactive plot is displayed correctly. For documents that contain multiple interactive plots, it is necessary to assign each plot a unique name using the prefix argument of export_interactive_roc. This is necessary to ensure that the JavaScript code manipulates the correct svg elements. The next code block shows an example knitr chunk that can be used in an .Rmd document to display an interactive plot.

```{r int-no, fig.keep='none', results = 'asis'}
cat(
  export_interactive_roc(basicplot, 
                        prefix = "a")
  )
```

The result is shown below:

cat(
  export_interactive_roc(basicplot, 
                        prefix = "a")
  )
Inf 1.7 1.2 0.9 0.7 0.3 0.2 -0.1 -0.5 -1.8 -1.1 -1 -0.9 -0.8 -0.7 -0.7 -0.7 -0.6 -0.6 -0.6 -0.5 -0.5 -0.5 -0.5 -0.5 -0.4 -0.4 -0.4 -0.3 -0.3 -0.2 -0.2 -0.2 -0.1 -0.1 -0.1 -0.1 0 0 0 0 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.5 0.5 0.5 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.9 0.9 0.9 0.9 0.9 1 1 1 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.3 1.3 1.3 1.3 1.3 1.3 1.5 1.5 1.5 1.6 1.7 1.7 1.7 1.8 1.8 1.8 1.9 1.9 1.9 1.9 2.1 2.4 2.6 2.8 0.00 0.10 0.25 0.50 0.75 0.90 1.00 0.00 0.10 0.25 0.50 0.75 0.90 1.00 False positive fraction True positive fraction

Click for confidence regions.

Multiple ROC curves

If you have grouping factors in your dataset, or you have multiple markers measured on the same subjects, you may wish to plot multiple ROC curves on the same plot. plotROC fully supports faceting and grouping done by ggplot2. In out example dataset, we have 2 markers measured in a paired manner:

head(test)
##   D   D.str         M1          M2
## 1 1     Ill 1.48117155 -2.50636605
## 2 1     Ill 0.61994478  1.46861033
## 3 0 Healthy 0.57613345  0.07532573
## 4 1     Ill 0.85433197  2.41997703
## 5 0 Healthy 0.05258342  0.01863718
## 6 1     Ill 0.66703989  0.24732453

These data are in wide format, with the 2 markers going across 2 columns. ggplot requires long format, with the marker result in a single column, and a third variable identifying the marker. We provide the function melt_roc to perform this transformation. The arguments are the data frame, a name or index identifying the disease status column, and a vector of names or indices identifying the the markers. Optionally, the names argument gives a vector of names to assign to the marker, replacing their column names. The result is a data frame in long format.

longtest <- melt_roc(test, "D", c("M1", "M2"))
head(longtest)
##     D          M name
## M11 1 1.48117155   M1
## M12 1 0.61994478   M1
## M13 0 0.57613345   M1
## M14 1 0.85433197   M1
## M15 0 0.05258342   M1
## M16 1 0.66703989   M1

Then, the dataset can be passed to the ggplot function, with the marker name given as a grouping or faceting variable.

ggplot(longtest, aes(d = D, m = M, color = name)) + geom_roc() + style_roc()

ggplot(longtest, aes(d = D, m = M, color = name)) + geom_roc(n.cuts = 0) + style_roc()

ggplot(longtest, aes(d = D, m = M)) + geom_roc() + facet_wrap(~ name) + style_roc()

ggplot(longtest, aes(d = D, m = M, linetype = name)) + geom_roc() + geom_rocci()

ggplot(longtest, aes(d = D, m = M, color = name)) + geom_roc() + style_roc()

pairplot <- ggplot(longtest, aes(d = D, m = M, color = name)) + 
  geom_roc(show.legend = FALSE) + style_roc()
direct_label(pairplot)

pairplot + geom_rocci()

pairplot + geom_rocci(linetype = 1)

pairplot + geom_rocci()

Interactive versions of the plots are fully supported.

cat(
  export_interactive_roc(direct_label(pairplot), 
                        prefix = "b")
  )
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
Inf 1.7 1.2 0.9 0.7 0.3 0.2 -0.1 -0.5 -1.8 Inf 2.4 1.9 1.3 0.7 0.3 -0.2 -0.6 -1.4 -3.8 M1 M2 -2.5 -2.4 -2.2 -1.8 -1.5 -1.5 -1.3 -1.2 -1.1 -1.1 -1 -1 -0.9 -0.8 -0.7 -0.5 -0.5 -0.5 -0.4 -0.4 -0.3 -0.3 -0.2 -0.2 -0.2 -0.2 -0.2 -0.1 -0.1 0 0 0 0 0 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.5 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.9 0.9 1 1 1 1 1.1 1.1 1.1 1.1 1.2 1.2 1.3 1.3 1.3 1.3 1.4 1.4 1.5 1.5 1.5 1.5 1.6 1.7 1.7 1.7 1.7 1.8 1.9 2 2.1 2.2 2.3 2.3 2.4 2.4 2.4 2.8 3.7 5.4 -2.5 -2.4 -2.2 -1.8 -1.5 -1.5 -1.3 -1.2 -1.1 -1.1 -1 -1 -0.9 -0.8 -0.7 -0.5 -0.5 -0.5 -0.4 -0.4 -0.3 -0.3 -0.2 -0.2 -0.2 -0.2 -0.2 -0.1 -0.1 0 0 0 0 0 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.5 0.6 0.6 0.7 0.7 0.7 0.7 0.7 0.7 0.9 0.9 1 1 1 1 1.1 1.1 1.1 1.1 1.2 1.2 1.3 1.3 1.3 1.3 1.4 1.4 1.5 1.5 1.5 1.5 1.6 1.7 1.7 1.7 1.7 1.8 1.9 2 2.1 2.2 2.3 2.3 2.4 2.4 2.4 2.8 3.7 5.4 0.00 0.10 0.25 0.50 0.75 0.90 1.00 0.00 0.10 0.25 0.50 0.75 0.90 1.00 False positive fraction True positive fraction

Click for confidence regions.

cat(
  export_interactive_roc(
    ggplot(longtest, aes(d = D, m = M)) + geom_roc() + facet_wrap(~ name), 
                        prefix = "c", width = 10, height = 5)
  )
Inf 1.7 1.2 0.9 0.7 0.3 0.2 -0.1 -0.5 -1.8 -3.6 -2.5 -2.4 -1.9 -1.5 -1.4 -1.3 -1.1 -1 -1 -0.8 -0.8 -0.8 -0.8 -0.7 -0.6 -0.6 -0.6 -0.5 -0.5 -0.5 -0.3 -0.2 -0.2 -0.2 -0.1 -0.1 -0.1 -0.1 0 0 0 0.1 0.1 0.1 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.6 0.6 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.9 0.9 0.9 0.9 1 1 1.1 1.2 1.2 1.3 1.3 1.3 1.4 1.4 1.4 1.4 1.5 1.5 1.5 1.5 1.6 1.6 1.6 1.7 1.7 1.7 1.7 1.9 1.9 1.9 2.1 2.1 2.3 2.3 2.3 2.4 2.4 Inf 2.4 1.9 1.3 0.7 0.3 -0.2 -0.6 -1.4 -3.8 -3.6 -2.5 -2.4 -1.9 -1.5 -1.4 -1.3 -1.1 -1 -1 -0.8 -0.8 -0.8 -0.8 -0.7 -0.6 -0.6 -0.6 -0.5 -0.5 -0.5 -0.3 -0.2 -0.2 -0.2 -0.1 -0.1 -0.1 -0.1 0 0 0 0.1 0.1 0.1 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5 0.6 0.6 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.9 0.9 0.9 0.9 1 1 1.1 1.2 1.2 1.3 1.3 1.3 1.4 1.4 1.4 1.4 1.5 1.5 1.5 1.5 1.6 1.6 1.6 1.7 1.7 1.7 1.7 1.9 1.9 1.9 2.1 2.1 2.3 2.3 2.3 2.4 2.4 M1 M2 0.00 0.10 0.25 0.50 0.75 0.90 1.00 0.00 0.10 0.25 0.50 0.75 0.90 1.00 0.00 0.10 0.25 0.50 0.75 0.90 1.00 False positive fraction True positive fraction

Click for confidence regions.

Showing multiple curves is also useful when there is a factor that affects the classification accuracy of the test. Let’s create another example dataset.

D.cov <- rbinom(400, 1, .5)
gender <- c("Male", "Female")[rbinom(400, 1, .49) + 1]
M.diff <- rnorm(400, mean = D.cov, sd = ifelse(gender == "Male", .5, 1.5))

test.cov <- data.frame(D = D.cov, gender = gender, M = M.diff)
bygend <- ggplot(test.cov, aes(d = D, m = M, color = gender)) + geom_roc(show.legend = FALSE)
direct_label(bygend) + style_roc()

New Features

AUC tables

test.cov2 <- data.frame(D = D.cov, gender = gender, M = M.diff, 
                        group = c("A", "B")[rbinom(400, 1, .4) + 1])
bygend2 <- ggplot(test.cov2, aes(d = D, m = M, color = gender)) + 
           geom_roc() + facet_wrap(~ group)
calc_auc(bygend2)
##   PANEL group group.1 gender       AUC
## 1     1     1       1 Female 0.6374434
## 3     1     2       2   Male 0.9531326
## 2     2     1       1 Female 0.6328904
## 4     2     2       2   Male 0.8732526

Custom labels and locations

ggplot(test, aes(d = D, m = M1)) + geom_roc(cutoffs.at = c(2, 1, .5, 0, -.5, -1))

ggplot(test, aes(d = D, m = M1)) + geom_roc(cutoffs.at = c(2, 1, .5, 0, -.5, -1), cutoff.labels = letters[1:6])

Increasing parameter

By default, the ROC calculation assumes that larger M is associated with a larger Pr(D = 1), you can now change this by setting increasing = FALSE.

ggplot(test, aes(d = D, m = M1)) + geom_roc(increasing = FALSE)

ggplot(test, aes(d = D, m = -M1)) + geom_roc(increasing = FALSE)

Advanced options

Themes and annotations

plotROC uses the ggplot2 package to create the objects to be plotted. Therefore, themes and annotations can be added in the usual ggplot2 way, or with external libraries such as ggthemes. The area under the ROC curve can be calculated by passing the ggplot object to the calc_auc function.

basicplot + 
  style_roc(theme = theme_grey) +
  theme(axis.text = element_text(colour = "blue")) +
  ggtitle("Themes and annotations") + 
  annotate("text", x = .75, y = .25, 
           label = paste("AUC =", round(calc_auc(basicplot)$AUC, 2))) +
  scale_x_continuous("1 - Specificity", breaks = seq(0, 1, by = .1))
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

Other estimation methods

By default calculate_roc computes the empirical ROC curve. There are other estimation methods out there, as I have summarized in the introduction. Any estimation method can be used, as long as the cutoff, the TPF and the FPF are returned. Then you can simply pass those values in a data frame to the ggroc function, and overriding the default stat to identity.

D.ex <- test$D
M.ex <- test$M1
mu1 <- mean(M.ex[D.ex == 1])
mu0 <- mean(M.ex[D.ex == 0])
s1 <- sd(M.ex[D.ex == 1])
s0 <- sd(M.ex[D.ex == 0])
c.ex <- seq(min(M.ex), max(M.ex), length.out = 300)

binorm.roc <- data.frame(c = c.ex, 
                             FPF = pnorm((mu0 - c.ex)/s0), 
                             TPF = pnorm((mu1 - c.ex)/s1)
                             )

binorm.plot <- ggplot(binorm.roc, aes(x = FPF, y = TPF, label = c)) + 
  geom_roc(stat = "identity") + style_roc(theme = theme_grey)
binorm.plot

Interactive plots with stat = "identity" are not currently supported.

Another potential use of this approach is for plotting time-dependent ROC curves for time-to-event outcomes estimated as described in Heagerty, Lumley, and Pepe (2000).

library(survivalROC)
survT <- rexp(350, 1/5)
cens <- rbinom(350, 1, .1)
M <- -8 * sqrt(survT) + rnorm(350, sd = survT)

sroc <- lapply(c(2, 5, 10), function(t){
  stroc <- survivalROC(Stime = survT, status = cens, marker = M, 
                       predict.time = t, method = "NNE", 
                       span = .25 * 350^(-.2))
  data.frame(TPF = stroc[["TP"]], FPF = stroc[["FP"]], 
             c = stroc[["cut.values"]], 
             time = rep(stroc[["predict.time"]], length(stroc[["FP"]])))
})

sroclong <- do.call(rbind, sroc)

ggplot(sroclong, aes(x = FPF, y = TPF, label = c, color = time)) + 
  geom_roc(labels = FALSE, stat = "identity") + style_roc(theme = theme_gray)

Acknowledgements

This package would not be possible without the following: