# ROC and AUC in R!!!

[sourcecode language=”R”]

library(pROC) # install with install.packages(“pROC”)

library(randomForest) # install with install.packages(“randomForest”)

#######################################

##

## Generate weight and obesity datasets.

##

#######################################

set.seed(420) # this will make my results match yours

num.samples <- 100

## genereate 100 values from a normal distribution with

## mean 172 and standard deviation 29, then sort them

weight <- sort(rnorm(n=num.samples, mean=172, sd=29))

## Now we will decide if a sample is obese or not.

## NOTE: This method for classifying a sample as obese or not

## was made up just for this example.

## rank(weight) returns 1 for the lightest, 2 for the second lightest, …

## … and it returns 100 for the heaviest.

## So what we do is generate a random number between 0 and 1. Then we see if

## that number is less than rank/100. So, for the lightest sample, rank = 1.

## This sample will be classified "obese" if we get a random number less than

## 1/100. For the second lightest sample, rank = 2, we get another random

## number between 0 and 1 and classify this sample "obese" if that random

## number is < 2/100. We repeat that process for all 100 samples

obese <- ifelse(test=(runif(n=num.samples) < (rank(weight)/num.samples)),

yes=1, no=0)

obese ## print out the contents of "obese" to show us which samples were

## classified "obese" with 1, and which samples were classified

## "not obese" with 0.

## plot the data

plot(x=weight, y=obese)

## fit a logistic regression to the data…

glm.fit=glm(obese ~ weight, family=binomial)

lines(weight, glm.fit$fitted.values)

#######################################

##

## draw ROC and AUC using pROC

##

#######################################

## NOTE: By default, the graphs come out looking terrible

## The problem is that ROC graphs should be square, since the x and y axes

## both go from 0 to 1. However, the window in which I draw them isn't square

## so extra whitespace is added to pad the sides.

roc(obese, glm.fit$fitted.values, plot=TRUE)

## Now let's configure R so that it prints the graph as a square.

##

par(pty = "s") ## pty sets the aspect ratio of the plot region. Two options:

## "s" – creates a square plotting region

## "m" – (the default) creates a maximal plotting region

roc(obese, glm.fit$fitted.values, plot=TRUE)

## NOTE: By default, roc() uses specificity on the x-axis and the values range

## from 1 to 0. This makes the graph look like what we would expect, but the

## x-axis itself might induce a headache. To use 1-specificity (i.e. the

## False Positive Rate) on the x-axis, set "legacy.axes" to TRUE.

roc(obese, glm.fit$fitted.values, plot=TRUE, legacy.axes=TRUE)

## If you want to rename the x and y axes…

roc(obese, glm.fit$fitted.values, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive Percentage", ylab="True Postive Percentage")

## We can also change the color of the ROC line, and make it wider…

roc(obese, glm.fit$fitted.values, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive Percentage", ylab="True Postive Percentage", col="#377eb8", lwd=4)

## If we want to find out the optimal threshold we can store the

## data used to make the ROC graph in a variable…

roc.info <- roc(obese, glm.fit$fitted.values, legacy.axes=TRUE)

str(roc.info)

## and then extract just the information that we want from that variable.

roc.df <- data.frame(

tpp=roc.info$sensitivities*100, ## tpp = true positive percentage

fpp=(1 – roc.info$specificities)*100, ## fpp = false positive precentage

thresholds=roc.info$thresholds)

head(roc.df) ## head() will show us the values for the upper right-hand corner

## of the ROC graph, when the threshold is so low

## (negative infinity) that every single sample is called "obese".

## Thus TPP = 100% and FPP = 100%

tail(roc.df) ## tail() will show us the values for the lower left-hand corner

## of the ROC graph, when the threshold is so high (infinity)

## that every single sample is called "not obese".

## Thus, TPP = 0% and FPP = 0%

## now let's look at the thresholds between TPP 60% and 80%

roc.df[roc.df$tpp > 60 & roc.df$tpp < 80,]

## We can calculate the area under the curve…

roc(obese, glm.fit$fitted.values, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive Percentage", ylab="True Postive Percentage", col="#377eb8", lwd=4, print.auc=TRUE)

## …and the partial area under the curve.

roc(obese, glm.fit$fitted.values, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive Percentage", ylab="True Postive Percentage", col="#377eb8", lwd=4, print.auc=TRUE, print.auc.x=45, partial.auc=c(100, 90), auc.polygon = TRUE, auc.polygon.col = "#377eb822")

#######################################

##

## Now let's fit the data with a random forest…

##

#######################################

rf.model <- randomForest(factor(obese) ~ weight)

## ROC for random forest

roc(obese, rf.model$votes[,1], plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive Percentage", ylab="True Postive Percentage", col="#4daf4a", lwd=4, print.auc=TRUE)

#######################################

##

## Now layer logistic regression and random forest ROC graphs..

##

#######################################

roc(obese, glm.fit$fitted.values, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive Percentage", ylab="True Postive Percentage", col="#377eb8", lwd=4, print.auc=TRUE)

plot.roc(obese, rf.model$votes[,1], percent=TRUE, col="#4daf4a", lwd=4, print.auc=TRUE, add=TRUE, print.auc.y=40)

legend("bottomright", legend=c("Logisitic Regression", "Random Forest"), col=c("#377eb8", "#4daf4a"), lwd=4)

#######################################

##

## Now that we're done with our ROC fun, let's reset the par() variables.

## There are two ways to do it…

##

#######################################

par(pty = "m")

[/sourcecode]