Diagnostic probability plots in logistic regression

Question

There is some discussion on StackExchange about diagnostic plots for logistic regression, but all are focusing on "residuals", for which there is not even a consensus how to define them for logistic regression. And whether they are useful seems to be another can of worms.

I wonder, however, why dignostic plots are not simply based on comparing the predicted probabilities with the observed probabilities dierectly estimated form the data. Two obvious approaches come to my mind ($x_i$ is the linear predictor, i.e. $x_i=\sum_j\beta_j x_{ij}$ and $y_i$ the binary response):

1. Compare the predited probability $P(y=1|x=x_i)$ with its non-parametric estimate from a conditional density plot on basis of the $x_i$, as e.g. computed by the R function cdplot.
2. Compare the predicted cumulative probability $P(y=1|x\leq x_i)$ with its empirical value computed from the data.

As I have not found these dignostics discussed in text books on logistic regression, there must be strong objections to these plots or their usefulness. Does someone know why these diagnostics are not useful?

PS: From its abstract, it seems that this article suggests method 1), but unfortunately I cannot check this because the article is behind a paywall:

Fowlkes, Edward B. “Some Diagnostics for Binary Logistic Regression via Smoothing.” Biometrika 74.3, pp. 503–15, 1987

Answer 1

Eventually, I have found a comprehensive description of the algorithm for creating a calibration plot in

J. Esarey, A. Pierce: "Assessing Fit Quality and Testing for Misspecification in Binary Dependent Variable Models." Political Analysis 20.4, pp. 480-500, 2012

The article compares it with classification based evaluation. Here is a summary of the ideas together with my comments and R code for creating a calibration plot.

When comparing the probability predicted by the model with the "observed" probability, there is the problem that no probabilities are observed but only zeros or ones, i.e. (non-) occurences of the response. These values can be smoothed out to probabilities by a distance weighted average in the "neighborhood" of each value, e.g. with a LOESS local regression. The distance for establishing the "neighborhood" and the weights can be measured in different spaces. Two obvious possible choices are

1. The distance on the link scale, i.e. on $\eta_i=\beta_0 + \langle\vec{\beta},\vec{x}_i\rangle$, where $x_i$ are the predictor variable values for the $i$-th observation, and $\beta$ are the model parameters.
2. The distance on the probability scale, i.e. on $p_i=P(Y=1|\vec{x}_i) = 1 / (1+e^{-\eta_i})$

A LOESS fit through the points $(y_i,\eta_i)$ or $(y_i,p_i)$ will then yield an estimator $\hat{p}_i$ for each $y_i$, which can be compared to the probability $p_i$ predicted by the model: There are two caveats, however:

• For degrees greater than zero, the LOESS fit can yield values outside [0,1]. For this reason, the first value is missing in both of the above plots: its estimated probability $\hat{p}_i$ is negative. This can be easily corrected by cutting off the probabilities at zero and one.
• LOESS only takes a certain percentage (parameter span) of neighbors into account.

The above plots have been created with the default span=0.75. Esarey & Pierce suggest two different optimization methods and link to a reference implementation in a footnote, but that link is meanwhile stalled. I have therefore implemented a very simple optimization criterion: the minimum MSE between $\hat{p}_i$ and $p_i$, i.e. $\sum_i(\hat{p}_i - p_i)^2$. The result on the Challenger Space Shuttle O-Ring dataset can be seen here: I have also included the 95% prediction interval for $p_i$ as predicted by the model. Esarey & Pierce also compute the percentage of values that lie outside an 80% confidence interval by means of a parametric bootstrap, but this might easier be computed directly from the confidence intervals for $p_i$. Here is the code to produce the calibration plot on the link level (right hand side):

# Challenger Space Shuttle O-ring data:  ok vs temp
data <- data.frame(y=c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1,
                       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1),
                   x=c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 
                       70, 70, 70, 70, 72, 73, 75, 75, 76, 76, 
                       78, 79, 81))
fit <- glm(y ~ x, data=data, family=binomial)

#
# calibration plot on link level
#
link.model <- predict(fit, data, type="link", se.fit=TRUE)
sort.key <- order(link.model$fit)
x <- link.model$fit[sort.key]

# prediction interval for probability
plot(link.model$fit, data$y, main="link level")
p.lower <- plogis(link.model$fit - qnorm(1-0.05/2) * 
            link.model$se.fit)[sort.key]
p.upper <- plogis(link.model$fit + qnorm(1-0.05/2) * 
            link.model$se.fit)[sort.key]
polygon(c(x,rev(x)), c(p.lower, rev(p.upper)), col="#dddddd", 
          border=NA)
points(link.model$fit, data$y) # replot overplotted points
lines(x, plogis(x), col="red")

# LOESS fit
optim.span <- optimize(resub.mse, c(0.1,1.0),
                       y=data$y, x=link.model$fit, 
                       p.model=p.model)
span <- optim.span$minimum
p.fit <- loess(y ~ x, data=data.frame(y=data$y, 
          x=link.model$fit), family="gaussian", degree=1, 
                span=span)
p.cutfit <- predict(p.fit, data.frame(x=x))
p.cutfit[p.cutfit < 0] <- 0
p.cutfit[p.cutfit > 1] <- 1
lines(x, p.cutfit)

legend("topleft", c("model", sprintf("LOESS (span=%4.2f)", 
       span)),
       col=c("red","black"), lty=1)

# the optimization function for estimating span
resub.mse <- function(span, y, x, p.model) {
  fit <- loess(y ~ x, family="gaussian", degree=1, span=span)
  return(sum((fit$fitted - p.model)^2))
}

Answer 2

Another approach, apparently not discussed in the literature, is the conditional density plot as provided out-of-the box by the R function cdplot.

The conditional density plot directly estimates $P(Y=\omega_i|x)$ for an arbitrary number of levels $\omega_i$ non-parametrically without assuming a statistical model. In the case of logistic regression, there are only two levels (0 and 1) and the regression fits a parametric model for $P(Y=1|x)$. The two estimators can thus be directly compared to see whether the logistic model matches the data.

cdplot estimates $P(Y=1|x)$ by means of Bayes' Theorem $$P(Y=1|x) = \frac{f(x|Y=1)\cdot P(Y=1)}{f(x)}$$ where $f$ denotes the probability densities, which are estimated by a kernel density estimator from the data. The only tricky part in this estimation is that both the estimator for $f(x)$ and for $f(x|Y=1)$ must use the same kernel bandwidth. Compared to the LOESS approach, this has two conceptual advantages:

• The result is guaranteed to yield a probability, and it can not happen (like for LOESS) that the value lies outside the range $[0,1]$
• It does not require a numerical interpretation of the levels as 0 and 1 in order to make numerically sense or to be applicable at all.

Like for the LOESS approach, the predictor must be a scalar value, for which in complete analogy the link value $\eta$ can be used. The kernel desity estimator requires to choose a bandwidth, for which the "plugin method" (bw="SJ" in the R function cdplot) is generally recommended in the literature (and in the documentation of density, too, although it uses a different default).

For comparison, I have implemented an additional bandwidth selection method that chooses that bandwidth which makes the cdplot most close to the logistic prediction. This can serve as a baseline what could be the best to be said about the logistic model ;-)

And here the code with the plot of the 95% confidence band from the logistic model omitted for better legibility:

# Space Shuttle Challenger temp vs oring-ok
data <- data.frame(y=c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1,
                       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1),
                   x=c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
                       70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81))
fit <- glm(y ~ x, data=data, family=binomial)

# helper function for finding the bandwidth
# that is closest to the logistic model
resub.mse <- function(bw, y, x, p.model) {
  cdfit <- cdplot(x, y, bw=bw, plot=FALSE)
  return(sum((cdfit[[levels(y.factor)[1]]](x) - p.model)^2))
}

#
# logistic prediction vs. link
#
link.model <- predict(fit, data, type="link", se.fit=TRUE)
p.model <- plogis(link.model$fit)
sort.key <- order(link.model$fit)
x <- link.model$fit[sort.key]
plot(link.model$fit, data$y, main="link level")
lines(x, plogis(x), col="red")

# cdplot vs. link
# note that we must code the level of interest
# as FIRST level (for cdplot)
y.factor <- factor(data$y, levels=c(1,0))
optim.bw <- optimize(resub.mse, c(bw.nrd0(x)/10, (max(x)-min(x))/2),
                     y=y.factor, x=link.model$fit, p.model=p.model)
bw <- optim.bw$minimum
p.kernel <- cdplot(link.model$fit, y.factor, bw="SJ", plot=FALSE)
lines(x, p.kernel$'1'(x))
p.kernel <- cdplot(link.model$fit, y.factor, bw=bw, plot=FALSE)
lines(x, p.kernel$'1'(x), col="blue", lty=2)

legend("topleft",
       c("model", "cdplot (bw='SJ')", sprintf("closest cdplot (bw=%4.2f)", bw)),
       col=c("red", "black", "blue"), lty=c(1,1,2))

Answer 3

Calibration

For completeness sake, here are two other ways to produce calibration plots: The first is using calibration belts, introduced by Nattino et al. (2014)$^{1}$, Nattino et al. (2016)$^{2}$ and Nattino et al. (2017)$^{3}$. Briefly, they fit an $m$-th-order polynomial logistic function (with $m\geq 2$) to the observed outcomes using the predicted probabilities of the model to be assessed. The parameter $m$ is selected using a standard forward selection procedure controlled by a likelihood-ratio statistic that accounts for the forward process used to select $m$. The calibration belt can be used for internal and external calibration. The procedure is implemented in Stata (calibrationbelt) and R (package givitiR). Here is the example using the Challenger data:

# Challenger Shuttle Challenger temp vs oring-ok
dat <- data.frame(y=c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1,
                       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1),
                   x=c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
                       70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81))

fit <- glm(y ~ x, data=dat, family=binomial)

preds_p <- predict(fit, type = "response")

cb <- givitiCalibrationBelt(o = dat$y, e = preds_p , devel = "internal", confLevels = c(0.95, 0.8))

plot(cb, main = "", las = 1, ylab = "Observed probabilities", xlab = "Model probabilities", table = FALSE)

The identity line is displayed in red. The light gray area is the 80% confidence calibration interval whereas the dark gray is a 95% confidence interval. Ideally, the red line is inside the belt over the whole range of probabilities. In the example, the confidence intervals are huge: The calibration belt shows a large uncertainty with respect to calibration. As the red line lies within the interval, we cannot reject the hypothesis of a well calibrated model.

To better illustrate the calibration belt, let's look at a calibration belt of a well calibrated model:

Here, the confidence intervals are much narrower. Because the red identity line lies within the belt over the whole range, it offers little evidence for miscalibration.

As some of the other answers, the second method implemented in the R package rms relies on a nonparametric smoother fitted to the predicted and observed probabilities. It also plots bias-corrected estimates based on the bootstrap. Details can be found in Harrel (2015)$^{4}$.

library(rms)

mod <- lrm(y~x, dat = dat, x = TRUE, y = TRUE)

res <- calibrate(mod, B = 10000)

plot(res)

The model seems to underestimate probabilities lower than $0.75$ and overestimate probabilities in over $0.75$. But again, due to the small sample size, the uncertainty is large.

Residuals

There are many types of residuals for generalized linear models but their interpretation is often difficult. One possibility is to look at simulation-based quantile residuals as implemented in the DHARMa package for R. Here is the example using the same data as above:

fit <- glm(y ~ x, data=dat, family=binomial)

simres <- simulateResiduals(fit, n = 1e4, seed = 142857)

plot(simres)

The nice thing about these residuals is that they can be interpreted as the "usual" residuals from linear regression models. On the left, a Q-Q-plot of the residuals is shown. On the right, the residuals are plotted against the predicted values. In both cases, there seems to be little evidence for a problem.

$[1]:$ Nattino, G., Finazzi, S., & Bertolini, G. (2014). A new calibration test and a reappraisal of the calibration belt for the assessment of prediction models based on dichotomous outcomes. Statistics in medicine, 33(14), 2390-2407.

$[2]:$ Nattino, G., Finazzi, S., & Bertolini, G. (2016). A new test and graphical tool to assess the goodness of fit of logistic regression models. Statistics in medicine, 35(5), 709-720.

$[3]:$ Nattino, G., Lemeshow, S., Phillips, G., Finazzi, S., & Bertolini, G. (2017). Assessing the calibration of dichotomous outcome models with the calibration belt. The Stata Journal, 17(4), 1003-1014.

$[4]:$ Harrell, F. E. (2015). Regression modeling strategies: with applications to linear models, logistic and ordinal regression, and survival analysis (Vol. 3). New York: springer.