One of the toughest problems in predictive model occurs when the classes have a severe imbalance. We spend an entire chapter on this subject itself. One consequence of this is that the performance is generally very biased against the class with the smallest frequencies. For example, if the data have a majority of samples belonging to the first class and very few in the second class, most predictive models will maximize accuracy by predicting everything to be the first class. As a result there's usually great sensitivity but poor specificity.
As a demonstration will use a simulation system described here. By default it has about a 50-50 class frequency but we can change this by altering the function argument called intercept:
(All the code used here was formatted by Pretty R at inside-R.org.)
In the training set the class frequency breakdown looks like this:
> table(training$Class)Class1 Class2 899 101
There is almost a 9:1 imbalance in these data.
Let's use a standard random forest model with these data using the default value of mtry. We'll also use 10-fold cross validation to get a sense of performance:
> set.seed(949) > mod0 <- train(Class ~ ., data = training, + method = "rf", + metric = "ROC", + tuneGrid = data.frame(mtry = 3), + trControl = trainControl(method = "cv", + classProbs = TRUE, + summaryFunction = twoClassSummary)) > getTrainPerf(mod0)TrainROC TrainSens TrainSpec method 1 0.9288911 0.99 0.4736364 rf
The area under the ROC curve is very high, indicating that the model has very good predictive power for these data. Here's a test set ROC curve for this model:
The plot shows the default probability cut off value of 50%. The sensitivity and specificity values associated with this point indicate that performance is not that good when an actual call needs to be made on a sample.
One of the most common ways to deal with this is to determine an alternate probability cut off using the ROC curve. But to do this well, another set of data (not the test set) is needed to set the cut off and the test set is used to validate it. We don't have a lot of data this is difficult since we will be spending some of our data just to get a single cut off value.
Alternatively the model can be tuned, using resampling, to determine any model tuning parameters as well as an appropriate cut off for the probabilities.
The latest update to the caret package allows users to define their own modeling and prediction components. This also gives us a huge amount of flexibility for creating your own models or doing some things that were originally intended by the package. This page shows a lot of the details for creating custom models.
Suppose the model has one tuning parameter and we want to look at four candidate values for tuning. Suppose we also want to tune the probability cut off over 20 different thresholds. Now we have to look at 20×4=80 different models (and that is for each resample). One other feature that has been opened up his ability to use sequential parameters: these are tuning parameters that don't require a completely new model fit to produce predictions. In this case, we can fit one random forest model and get it's predicted class probabilities and evaluate the candidate probability cutoffs using these same hold-out samples. Again, there's a lot of details on this page and, without going into them, our code for these analyses can be found here.
Basically, we define a list of model components (such as the fitting code, the prediction code, etc.) and feed this into the train function instead of using a pre-listed model string (such as method = "rf"). For this model and these data, there was an 8% increase in training time to evaluate 20 additional values of the probability cut off.
How do we optimize this model? Normally we might look at the area under the ROC curve as a metric to choose our final values. In this case the ROC curve is independent of the probability threshold so we have to use something else. A common technique to evaluate a candidate threshold is see how close it is to the perfect model where sensitivity and specificity are one. Our code will use the distance between the current model's performance and the best possible performance and then have train minimize this distance when choosing it's parameters. Here is the code that we use to calculate this:
fourStats <- function (data, lev = levels(data$obs), model = NULL) { ## This code will get use the area under the ROC curve and the ## sensitivity and specificity values using the current candidate ## value of the probability threshold. out <- c(twoClassSummary(data, lev = levels(data$obs), model = NULL)) ## The best possible model has sensitivity of 1 and specifity of 1. ## How far are we from that value? coords <- matrix(c(1, 1, out["Spec"], out["Sens"]), ncol = 2, byrow = TRUE) colnames(coords) <- c("Spec", "Sens") rownames(coords) <- c("Best", "Current") c(out, Dist = dist(coords)[1]) }
Now let's run our random forest model and see what it comes up with for the best possible threshold:
set.seed(949) mod1 <- train(Class ~ ., data = training, ## 'modelInfo' is a list object found in the linked ## source code method = modelInfo, ## Minimize the distance to the perfect model metric = "Dist", maximize = FALSE, tuneLength = 20, trControl = trainControl(method = "cv", classProbs = TRUE, summaryFunction = fourStats))
The resulting model output notes that:
Tuning parameter 'mtry' was held constant at a value of 3 Dist was used to select the optimal model using the smallest value. The final values used for the model were mtry = 3 and threshold = 0.887.
Using ggplot(mod1) will show the performance profile. Instead here is a plot of the sensitivity, specificity, and distance to the perfect model:
You can see that as we increase the probability cut off for the first class it takes more and more evidence for a sample to be predicted as the first class. As a result the sensitivity goes down when the threshold becomes very large. The upside is that we can increase specificity in the same way. The blue curve shows the distance to the perfect model. The value of 0.887 was found to be optimal.
Now we can use the test set ROC curve to validate the cut off we chose by resampling. Here the cut off closest to the perfect model is 0.881. We were able to find a good probability cut off value without setting aside another set of data for tuning the cut off.
One great thing about this code is that it will automatically apply the optimized probability threshold when predicting new samples. Here is an example:
Class1 Class2 Class Note 1 0.874 0.126 Class2 * 2 1.000 0.000 Class1 3 0.930 0.070 Class1 4 0.794 0.206 Class2 * 5 0.836 0.164 Class2 * 6 0.988 0.012 Class1
However we should be careful because the probability values are not consistent with our usual notion of a 50-50 cut off.
caret webinar on Feb 25
I"ll be doing a webinar with the Orange County R User Group on the caret package on Tue, Feb 25, 2014 1:00 PM - 2:00 PM EST.
Here is the url in case you are interested: https://www3.gotomeeting.com/register/673845982
Thanks to Ray DiGiacomo for setting this up.
Calibration Affirmation
In the book, we discuss the notion of a probability model being "well calibrated". There are many different mathematical techniques that classification models use to produce class probabilities. Some of values are "probability-like" in that they are between zero and one and sum to one. This doesn't necessarily mean that the probability estimates are consistent with the true event rate seen with similar samples. Suppose a doctor told you that you have a model predicts that you have a 90% chance of having cancer. Is that number consistent with the actual likelihood of you being sick?
As an example, we can look at the cell image segmentation data use in the text (from this paper). The data are in the caret package. The outcome is whether or not a cell in an image is well segmented (WS) or poorly segments (PS). We try to predict this using measured characteristics of the supposed cell.
We can fit a partial least squares discriminant analysis model to these data. This model can use the basic PLS model for regression (based on binary dummy variables for the classes) and we can post-process the data into class probabilities and class predictions. One method to do this is the softmax function. This simple mathematical formula can used but, in my experience, doesn't lead to very realistic predictions. Let's fit that model:
library(caret) data(segmentationData) ## Retain the original training set segTrain <- subset(segmentationData, Case == "Train") ## Remove the first three columns (identifier columns) segTrainX <- segTrain[, -(1:3)] segTrainClass <- segTrain$Class ## Test Data segTest <- subset(segmentationData, Case != "Train") segTestX <- segTest[, -(1:3)] segTestClass <- segTest$Class
set.seed(1) useSoftmax <- train(segTrainX, segTrainClass, method = "pls", preProc = c("center", "scale"), ## Tune over 15 values of ncomp using 10 fold-CV tuneLength = 15, trControl = trainControl(method = "cv"), probMethod = "softmax") ## Based on cross-validation, how well does the model fit? getTrainPerf(useSoftmax)
## TrainAccuracy TrainKappa method ## 1 0.8108 0.592 pls
## Predict the test set probabilities SoftmaxProbs <- predict(useSoftmax, segTestX, type = "prob")[,"PS",] ## Put these into a data frame: testProbs <- data.frame(Class = segTestClass, Softmax = SoftmaxProbs)
(All the code used here was formatted by Pretty R at inside-R.org.)
The accuracy and Kappa statistics reflect a reasonable model.
We are interested in knowing if cells are poorly segmented. What is the distribution of the test set class probabilities for the two classes?
Looking at the bottom panel, the mode of the distribution is around 40%. Also, very little of the truly well segmented samples are not confidently predicted as such since most probabilities are greater than 30%. The same is true for the poorly segmented cells. Very few values are greater than 80%.
What does the calibration plot look like?
This isn't very close to the 45 degree reference line so we shouldn't expect the probabilities to be very realistic.
Another approach for PLS models is to use Bayes' rule to estimate the class probabilities. This is a little more complicated but at least it is based on probability theory. We can fit the model this way using the option probMethod = "Bayes":
## TrainAccuracy TrainKappa method ## 1 0.8108 0.5984 pls
getTrainPerf(useSoftmax)
## TrainAccuracy TrainKappa method ## 1 0.8108 0.592 pls
BayesProbs <- predict(useBayes, segTestX, type = "prob")[,"PS",] ## Put these into a data frame as before: testProbs$Bayes <- BayesProbs
Performance is about the same. The test set ROC curves both have AUC values of 0.87.
What do these probabilities look like?
Even though the accuracy and ROC curve results show the models to be very similar, the class probabilities are very different. Here, the most confidently predicted samples have probabilities closer to zero or one. Based on the calibration plot, is this model any better?
Much better! Looking at the right-hand side of the plot, the event rate of samples with Bayesian generated probabilities is about 97% for class probabilities between 91% and 100%. The softmax model has an event rate of zero.
This is a good example because, by most measures,the model performance is almost identical. However, one generates much more realistic values than the other.
Applied Predictive Modeling
Applied Predictive Modeling is a book on the practice of modeling when accuracy is the primary goal.