Diagnostic probability plots in logistic regression 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):

*

*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.

*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


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))

