rocplot.single <- function(grp, pred, title = "ROC Plot", p.value = FALSE){
  require(ggplot2)
  plotdata <- rocdata(grp, pred)
 
  if (p.value == TRUE){
    annotation <- with(plotdata$stats, paste("AUC=",signif(auc, 2), " (P=", signif(p.value, 2), ")", sep=""))
  } else {
    annotation <- with(plotdata$stats, paste("AUC=",signif(auc, 2), " (95%CI ", signif(ci.upper, 2), " - ", signif(ci.lower, 2), ")", sep=""))
  }
 
  p <- ggplot(plotdata$roc, aes(x = x, y = y)) +
      geom_line(aes(colour = "")) +
      geom_abline (intercept = 0, slope = 1) +
      theme_bw() +
      scale_x_continuous("False Positive Rate (1-Specificity)") +
      scale_y_continuous("True Positive Rate (Sensitivity)") +
      scale_colour_manual(labels = annotation, values = "#000000") +
      opts(title = title,
           plot.title = theme_text(face="bold", size=14), 
           axis.title.x = theme_text(face="bold", size=12),
           axis.title.y = theme_text(face="bold", size=12, angle=90),
           panel.grid.major = theme_blank(),
           panel.grid.minor = theme_blank(),
           legend.justification=c(1,0), 
           legend.position=c(1,0),
           legend.title=theme_blank(),
           legend.key = theme_blank()
           )
  return(p)
}

Let’s make an example dataset and apply this function to it. The following function generates some random data in two groups.

exampledata <- function(n = 50, sd = 1) {
  # Control group
  ctrl <- data.frame(grp = rep("Control", n),
                    res = -1 + rnorm(n, mean = 0, sd = sd)
                    )
  # Disease group
  dis <- data.frame(grp = rep("Disease", n),
                    res = 1 + rnorm(n, mean = 0, sd = sd)
                    )
  df <- rbind(ctrl,dis)
  return(df)
}

We can then take the data frame that this produces and pass it to the rocplot.simple function. It processes the data through our rocdata function from part 1 and then passes the x and y co-ordinates to ggplot2. I’ve added a line of non-discrimination (geom_abline) and used the legend to display the statistics (AUC and confidence intervals). In order to get the legend to display for a single plot line I’ve added an aes colour mapping to geom_line and then used the labels argument of scale_colour_manual to add the annotation to the plot. The annotation can then be customised to display whatever I need – in this case either the AUC and CI or the AUC and P value.

Test1 <- exampledata()
p <- rocplot.single(Test1$grp, Test1$res, title = "")
print(p)

Let’s try a slightly more complicated situation. How about a hypothetical situation where you apply the same test in a set of different situations. How does the test perform? Well, let’s make up 3 test scenarios and build up 3 sets of test data. We can then put those together into a list and apply our function. Because we’re dealing with a list of data objects our rocplot function is slightly different:

rocplot.multiple <- function(test.data.list, groupName = "grp", predName = "res", title = "ROC Plot", p.value = TRUE){
  require(plyr)
  require(ggplot2)
  plotdata <- llply(test.data.list, function(x) with(x, rocdata(grp = eval(parse(text = groupName)), pred = eval(parse(text = predName)))))
  plotdata <- list(roc = ldply(plotdata, function(x) x$roc),
                   stats = ldply(plotdata, function(x) x$stats)
                   )
 
  if (p.value == TRUE){
    annotation <- with(plotdata$stats, paste("AUC=",signif(auc, 2), " (P=", signif(p.value, 2), ")", sep=""))
  } else {
    annotation <- with(plotdata$stats, paste("AUC=",signif(auc, 2), " (95%CI ", signif(ci.upper, 2), " - ", signif(ci.lower, 2), ")", sep=""))
  }
 
  p <- ggplot(plotdata$roc, aes(x = x, y = y)) +
       geom_line(aes(colour = .id)) +
       geom_abline (intercept = 0, slope = 1) +
       theme_bw() +
       scale_x_continuous("False Positive Rate (1-Specificity)") +
       scale_y_continuous("True Positive Rate (Sensitivity)") +
       scale_colour_brewer(palette="Set1", breaks = names(test.data.list), labels = paste(names(test.data.list), ": ", annotation, sep = "")) +
       opts(title = title,
            plot.title = theme_text(face="bold", size=14), 
            axis.title.x = theme_text(face="bold", size=12),
            axis.title.y = theme_text(face="bold", size=12, angle=90),
            panel.grid.major = theme_blank(),
            panel.grid.minor = theme_blank(),
            legend.justification=c(1,0), 
            legend.position=c(1,0),
            legend.title=theme_blank(),
            legend.key = theme_blank()
            )
  return(p)
}

We can generate some test data like this

TestData1 <- list(TrialA = exampledata(n=50, sd = 1), 
                  TrialB = exampledata(n=50, sd = 1.5),
                  TrialC = exampledata(n=50, sd = 2)
                  )

Each data frame containing our data is then put together into a list object which we pass to our rocplot.multiple function.

p <- rocplot.multiple(TestData1, title = "", p.value = FALSE)
print(p)

In addition to ggplot2 this function makes use of the excelent tools available from the plyr library also written by Hadley Wickham. The first step is to apply our rocdata function to each of the 3 data frames in our list using the llply function (equivalent to lapply in R base). We then put together the roc elements and the stats elements from each rocdata output by using the ldply function on each of them. This gives us 2 dataframes with an .id column corresponding to the name of the original data frame it came from. We can put those together into a list so just like our first example we supply ggplot with a list with all the ROC co-ordinates in one element and all the stats in the second. We can then plot one line for each set of data defined by the .id column. Like the first example I’ve used the the legend to display our summary statistics.

Code tested using R version 2.14.2 with the following packages:
plyr – version 1.7.1
ggplot2 – version 0.9.0

If you’re not sure what a ROC plot is a good summary is given here. So now that we know what we will be building, just how can we accomplish this in R? Let’s start by building a function that computes the co-ordinates for the ROC curve and gives us some summary statistics.

rocdata <- function(grp, pred){
  # Produces x and y co-ordinates for ROC curve plot
  # Arguments: grp - labels classifying subject status
  #            pred - values of each observation
  # Output: List with 2 components:
  #         roc = data.frame with x and y co-ordinates of plot
  #         stats = data.frame containing: area under ROC curve, p value, upper and lower 95% confidence interval
 
  grp <- as.factor(grp)
  if (length(pred) != length(grp)) {
    stop("The number of classifiers must match the number of data points")
  } 
 
  if (length(levels(grp)) != 2) {
    stop("There must only be 2 values for the classifier")
  }
 
  cut <- unique(pred)
  tp <- sapply(cut, function(x) length(which(pred > x & grp == levels(grp)[2])))
  fn <- sapply(cut, function(x) length(which(pred < x & grp == levels(grp)[2])))
  fp <- sapply(cut, function(x) length(which(pred > x & grp == levels(grp)[1])))
  tn <- sapply(cut, function(x) length(which(pred < x & grp == levels(grp)[1])))
  tpr <- tp / (tp + fn)
  fpr <- fp / (fp + tn)
  roc = data.frame(x = fpr, y = tpr)
  roc <- roc[order(roc$x, roc$y),]
 
  i <- 2:nrow(roc)
  auc <- (roc$x[i] - roc$x[i - 1]) %*% (roc$y[i] + roc$y[i - 1])/2
 
  pos <- pred[grp == levels(grp)[2]]
  neg <- pred[grp == levels(grp)[1]]
  q1 <- auc/(2-auc)
  q2 <- (2*auc^2)/(1+auc)
  se.auc <- sqrt(((auc * (1 - auc)) + ((length(pos) -1)*(q1 - auc^2)) + ((length(neg) -1)*(q2 - auc^2)))/(length(pos)*length(neg)))
  ci.upper <- auc + (se.auc * 0.96)
  ci.lower <- auc - (se.auc * 0.96)
 
  se.auc.null <- sqrt((1 + length(pos) + length(neg))/(12*length(pos)*length(neg)))
  z <- (auc - 0.5)/se.auc.null
  p <- 2*pnorm(-abs(z))
 
  stats <- data.frame (auc = auc,
                       p.value = p,
                       ci.upper = ci.upper,
                       ci.lower = ci.lower
                       )
 
  return (list(roc = roc, stats = stats))
}

So let’s go through how this works. The function requires two arguments: Firstly a numeric vector giving the results of the test being applied, and secondly a factor vector (or a vector that can be coerced to a factor with two levels) containing the status of the instance being tested eg. disease or healthy, case or control etc. We firstly ensure this grouping vector is a factor

grp <- as.factor(grp)

and then run some checks to ensure there is one classifier for each test result and that there are only two levels of the classifier.

if (length(pred) != length(grp)) {
    stop("The number of classifiers must match the number of data points")
} 
 
if (length(levels(grp)) != 2) {
   stop("There must only be 2 values for the classifier")
}

The ROC plot should show the proportion of true positives against the proportion of false positives for any threshold value. However using the dataset we have we can define a set of thresholds at all the unique values of the test result.

cut <- unique(pred)

Then at each cutoff point we want to calculate the number of true positives, false positives, true negatives and false negatives. We can do this using the sapply function in R. Looking at the true positives as an example, for every cutoff value we count the number of entries where the predictor vector is greater than the cutoff AND the group vector is the higher factor level (hence the reason for ensuring it is an ordered factor).

tp <- sapply(cut, function(x) length(which(pred > x & grp == levels(grp)[2])))
fn <- sapply(cut, function(x) length(which(pred < x & grp == levels(grp)[2])))
fp <- sapply(cut, function(x) length(which(pred > x & grp == levels(grp)[1])))
tn <- sapply(cut, function(x) length(which(pred < x & grp == levels(grp)[1])))

Each of the above produces a vector with the same length as cut. We can then calculate the true positive rate as the vector of true positives divided by the sum of the true positives and false negatives vector

tpr <- tp / (tp + fn)

and the false positive rate as the false positive vector divided by the sum of the false positive vector and the true negative vector.

fpr <- fp / (fp + tn)

these two vectors give our y and x co-ordinates for our ROC plot. It’s then a simple matter of writing them into a data frame and ordering them for plotting to produce our basic ROC plot data.

roc = data.frame(x = fpr, y = tpr)
roc <- roc[order(roc$x, roc$y),]

Summary Statistics

Probably the commonest way of summarising the ROC plot is to compute the area under the curve (AUC). This reduces the ROC plot to a single scalar value representing the performance of the test. A perfect test has an AUC of 1 and a test that is equivalent to random guessing has an AUC of 0.5. There are a number of ways of estimating the AUC of a plot; the trapezoid rule is employed here. The function used here is the same as the trapz function from the caTools package.

i <- 2:nrow(roc)
auc <- (roc$x[i] - roc$x[i - 1]) %*% (roc$y[i] + roc$y[i - 1])/2

An interesting statistical property of the AUC is that it represents the probability that the test value of a positive case will rank higher than a negative case. This is essentially equivalent to the Wilcoxon-Mann-Whitney statistic. Hence an alternative way of computing the AUC would be

pos <- pred[grp == levels(grp)[2]]
neg <- pred[grp == levels(grp)[1]]
auc.wilcox <- wilcox.test(pos,neg,exact=0)$statistic /(length(pos)*length(neg))

The standard error of the ROC curve is given by: (Hanley and McNeil, 1982)

(1)  

where is the AUC, is the number of positive cases, is the number of negative cases, and . Using the standard error we can calculate the upper and lower 95% confidence intervals for the AUC

pos <- pred[grp == levels(grp)[2]]
neg <- pred[grp == levels(grp)[1]]
q1 <- auc/(2-auc)
q2 <- (2*auc^2)/(1+auc)
se.auc <- sqrt(((auc * (1 - auc)) + ((length(pos) -1)*(q1 - auc^2)) + ((length(neg) -1)*(q2 - auc^2)))/(length(pos)*length(neg)))
ci.upper <- auc + (se.auc * 0.96)
ci.lower <- auc - (se.auc * 0.96)

We can also use equation (1) to calculate a p value for the AUC – ie. the probabality that the AUC = 0.5 (no discrimination). Substituting 0.5 into (1) and simplifying gives us the standard error for our null hypothesis.

(2)  

We can then calculate a Z score and from that the p value using the normal distribution

pos <- pred[grp == levels(grp)[2]]
neg <- pred[grp == levels(grp)[1]]
se.auc.null <- sqrt((1 + length(pos) + length(neg))/(12*length(pos)*length(neg)))
z <- (auc - 0.5)/se.auc.null
p <- 2*pnorm(-abs(z))

To finish off our rocdata function we just need to output our roc curve co-ordinates, auc and p value. I’ve chosen to put the stats measures into one data frame and the roc curve data into another and then put the two dataframes into a list.

stats <- data.frame (auc = auc,
                     p.value = p,
                     ci.upper = ci.upper,
                     ci.lower = ci.lower
                     )
 
return (list(roc = roc, stats = stats))

References
Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology 1982;143:29-36.

Setting up shared folders on Dropbox is pretty straightforward but if you’ve never done it before there’s a guide here. The next problem we found was that we had to make sure the script could correctly reference the path to data files and other scripts. Fortunately if your Dropbox folder is located in the same place as the R home directory (is is for a default install on Windows and Mac machines – I haven’t tested this on a Linux box yet) then it should be straightforward process.

You can check your R home directory using (example for Windows XP)

Sys.getenv("R_USER")

[1] "C:/Documents and Settings/Kate/My Documents"

Your Dropbox folder should be in My Documents too. Then in your R script set your working directory or file locations as follows:

setwd(path.expand("~/Dropbox/SharedRScript"))

That way you can define the location of any other called scripts or data relative to your shared folder and they should work irrespective of which machine you are on.  So far Dropbox is working pretty well for our needs. It’s free (at least for 2Gb storage), synchronises quickly and is easy to use. You might run into problems if several people are working on a file at the same time – it will save multiple versions rather than try to merge the files.  However for us it’s a good and easy solution to collaborative coding.