# Visualizing the calibration of predicted probability of a model

Suppose I have a predictive model that produces, for each instance, a probability for each class. Now I recognize that there are many ways to evaluate such a model if I want to use those probabilities for classification (precision, recall, etc.). I also recognize that an ROC curve and the area under it can be used to determine how well the model differentiates between classes. Those are not what I'm asking about.

I'm interested in assessing the calibration of the model. I know that a scoring rule like the Brier score can be useful for this task. That's OK, and I'll likely incorporate something along those lines, but I'm not sure how intuitive such metrics will be for the lay person. I'm looking for something more visual. I want the person interpreting results to be able to see whether or not when the model predicts something is 70% likely to happen that it actually happens ~70% of the time, etc.

I heard of (but never used) Q-Q plots, and at first I thought this was what I was looking for. However, it seems that is really meant for comparing two probability distributions. That's not directly what I have. I have, for a bunch of instances, my predicted probability and then whether the event actually occurred:

Index    P(Heads)    Actual Result
2          .3            Tails
4         .65            Tails
...         ...              ...


So is a Q-Q plot really what I want, or am I looking for something else? If a Q-Q plot is what I should be using, what is the correct way to transform my data into probability distributions?

I imagine I could sort both columns by predicted probability and then create some bins. Is that the type of thing I should be doing, or am I off in my thinking somewhere? I'm familiar with various discretization techniques, but is there a specific way of discretizing into bins that is standard for this sort of thing?

John Tukey recommended binning by halves: split the data into upper and lower halves, then split those halves, then split the extreme halves recursively. Compared to equal-width binning, this allows visual inspection of tail behavior without devoting too many graphical elements to the bulk of the data (in the middle).

Here is an example (using R) of Tukey's approach. (It's not exactly the same: he implemented mletter a little differently.)

First, let's create some predictions and some outcomes that conform to those predictions:

set.seed(17)
prediction <- rbeta(500, 3/2, 5/2)
actual <- rbinom(length(prediction), 1, prediction)
plot(prediction, actual, col="Gray", cex=0.8)


The plot is not very informative, because all the actual values are, of course, either $0$ (did not occur) or $1$ (did occur). (It appears as the background of gray open circles in the first figure below.) This plot needs smoothing. To do so, we bin the data. Function mletter does the splitting-by-halves. Its first argument r is an array of ranks between 1 and n (the second argument). It returns unique (numeric) identifiers for each bin:

mletter <- function(r,n) {
lower <-  2 + floor(log(r/(n+1))/log(2))
upper <- -1 - floor(log((n+1-r)/(n+1))/log(2))
i <- 2*r > n
lower[i] <- upper[i]
lower
}


Using this, we bin both the predictions and the outcomes and average each within each bin. Along the way, we compute bin populations:

classes <- mletter(rank(prediction), length(prediction))
pgroups <- split(prediction, classes)
agroups <- split(actual, classes)
bincounts <- unlist(lapply(pgroups, length)) # Bin populations
x <- unlist(lapply(pgroups, mean))           # Mean predicted values by bin
y <- unlist(lapply(agroups, mean))           # Mean outcome by bin


To symbolize the plot effectively we should make the symbol areas proportional to bin counts. It can be helpful to vary the symbol colors a little, too, whence:

binprop <- bincounts / max(bincounts)
colors <- -log(binprop)/log(2)
colors <- colors - min(colors)
colors <- hsv(colors / (max(colors)+1))


With these in hand, we now enhance the preceding plot:

abline(0,1, lty=1, col="Gray")                           # Reference curve
points(x,y, pch=19, cex = 3 * sqrt(binprop), col=colors) # Solid colored circles
points(x,y, pch=1, cex = 3 * sqrt(binprop))              # Circle outlines


As an example of a poor prediction, let's change the data:

set.seed(17)
prediction <- rbeta(500, 5/2, 1)
actual <- rbinom(length(prediction), 1, 1/2 + 4*(prediction-1/2)^3)


Repeating the analysis produces this plot in which the deviations are clear:

This model tends to be overoptimistic (average outcome for predictions in the 50% to 90% range are too low). In the few cases where the prediction is low (less than 30%), the model is too pessimistic.

• (+1) Very nice, thanks. I think the colors might distract a bit from the purpose, but the rest was a nice idea and very nice explanation. – Michael McGowan Mar 29 '12 at 20:52
• Michael, I found that some color was necessary to help see the very small circles that appear at either end. A constant color would accomplish this, of course. Just replace col=colors by the color you want, such as col="Red". – whuber Mar 29 '12 at 21:12
• +1, this is very nice. However, I don't quite intuit why the reference line a simple, straight 45 degree line, instead of the proper logistic regression line, or a loess? I should think those would be more appropriate references against which to judge the quality of the predictions. – gung - Reinstate Monica Mar 31 '12 at 3:33
• @gung Think about it: a prediction is accurate when the outcomes conform with the prediction. Therefore, among a collection of predictions all near $p$, it had better be the case that the average outcome is also close to $p$. Any other relationship would reflect inaccuracies. In particular, just how would you propose drawing a "proper logistic regression line"--which sounds like it might extend to $\pm\infty$ along one or both axes--on a plot that must be contained within the unit square $[0,1]\times [0,1]$? – whuber Mar 31 '12 at 15:04
• @gung (Take 2) What I think you might have in mind is improving the visualization to account for expected variation in the residuals. That variation should be proportional to $\sqrt{p(1-p)/n}$ for a point associated with predictions near $p$ consisting of $n$ observations. It is an interesting challenge to think of a useful, preferably gestalt way, to represent that on a statistical graphic like this. (I suppose you could erect prediction limits around each point...) – whuber Mar 31 '12 at 15:06

Another option is isotonic regression. It is similar to whuber's answer except the bins are generated dynamically instead of by splitting in halves, with a requirement that outputs are strictly increasing.

This primary usage of isotonic regression is to recalibrate your probabilities if they are shown to be poorly calibrated, but it can also be used for visualization. Basically, if the isotonic regression line roughly follows the Y=X line, then your probabilities are properly calibrated.

This is Isotonic Regression applied to the problem shown by Whuber.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.isotonic import IsotonicRegression

prediction = np.random.beta(3.0/2.0, 5.0/2.0, size=500)
actual = np.random.binomial(1,prediction, len(prediction))
plt.scatter(prediction, actual,  facecolors='none', edgecolors=[0.3,0.3,0.3], label='Data')

ir = IsotonicRegression()
isotonic = ir.fit_transform(prediction, actual)
plt.plot(prediction, isotonic,'ok', label='Isotonic Fit')

plt.xlabel('Prediction')
plt.ylabel('Actual')
plt.plot([0,1],[0,1], '--k', label='y=x line')
plt.legend(loc = 'center left')


http://fa.bianp.net/blog/2013/isotonic-regression/

http://stat.wikia.com/wiki/Isotonic_regression

You might also want to look at package "verification":

http://cran.r-project.org/web/packages/verification/index.html

There are plots in the vignette that might be useful:

http://cran.r-project.org/web/packages/verification/vignettes/verification.pdf

• The vignette link no longer works, and the first CRAN link doesn't seem to have a vignette anymore. – BLT Mar 28 '17 at 23:00