DHARMa residuals plot vs. binned residuals using stan_glm object

Question

I am quite a newbie in R and even more so in Bayesian regression. I have fit a stan_glm binomial model with 1689 observations, 12 variables and two interaction terms. All predictors are categorical. One of the main predictors suffers from quasi-complete separation.

I also followed the procedures detailed in Gelman, Hill and Vehtari (2020) to fit the model. One of the diagnostic tools introduced in the book is the binned residuals plot.

To produce a binned residuals plot, I tried performance::binned_residuals and arm::binnedplot. The two plots are identical. Unfortunately, it turned out that only 32% of the residuals fall within the error bound, with the rest on the two extremes, outside the bounds.

I don't know if it is due to the problem of separation or other issues in the model. I stumbled across this post and someone recommended the package DHARMa. Interpreting a binned residual plot in logistic regression

And so, I tried it. When I ran

simulationOutput <- simulateResiduals(fittedModel = fit.final, plot=F, integerResponse = NULL)

This warning message pops up:

Warning message:
In checkModel(fittedModel) :
  DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work!

I continued running the following code:

residuals(simulationOutput)
plot(simulationOutput)

And the plots came out perfect.

My questions are thus:

(i) What is the problem with the binned residuals plot? Anyway to fix it? (ii) Is the use of DHARMa residuals warranted and appropriate despite the warning?

============ UPDATE ============

Thank you, Shawn, for the extra efforts in labelling the steps.

Following your advice by plugging my model into createDHARMa directly:

I transformed the response into integers from factors

DHARMaRes = createDHARMa( simulatedResponse = t(posterior_predict(fit.final)), observedResponse = x$adj_pla1, fittedPredictedResponse = apply(t(posterior_predict(fit.final)), 1, median), integerResponse=T ) plot(DHARMaRes)

I got this as plot. Using median (recommended for Bayesian models) in fittedPredictedResponse = returns a boxplot.

If I use mean, I get the usual scatter plots.

In any case, can I use either of the plots? (Do excuse me for this rather dumb question. Please bear with me.)

Answer 1

Bayesian DHARMA Residuals

I recommend reading through this specific section in the DHARMA package vignette to understand why it is saying this, along with the supplementary vignette they mention:

The package creator has some code below that section for modeling with Bayesian data, which I've slightly modified to label what's going on and to include a random seed for reproducibility:

#### Load DHARMa Library and Set Random Seed ####
library(DHARMa)
set.seed(123)

#### Create Test Data and Fit to GLM ####
testData <- createData(sampleSize = 200, 
                      overdispersion = 0.5,
                      family = poisson())

fittedModel <- glm(observedResponse 
                   ~ Environment1,
                   family = "poisson",
                   data = testData)

#### Create Simulation Function ####
simulatePoissonGLM <- function(fittedModel, n){
  pred = predict(fittedModel, type = "response")
  nObs = length(pred)
  sim = matrix(nrow = nObs, ncol = n)
  for(i in 1:n) sim[,i] = rpois(nObs, pred)
  return(sim)
}

#### Use FUnction for Fitted Model ####
sim <- simulatePoissonGLM(fittedModel, 100)

#### Create DHARMa Residuals ####
DHARMaRes <- createDHARMa(simulatedResponse = sim, observedResponse = testData$observedResponse, 
                         fittedPredictedResponse = predict(fittedModel), integerResponse = T)

#### Plot Them ####
plot(DHARMaRes, quantreg = F)

You can see the residuals plotted below:

For your specific case, I think you just have to include your model into the createDHARMa function and then this will use your residuals in the way you prescribe, rather than using the simulateResiduals function typically used.

Binned Residual Plot

As for the binned residual plot, notice this section in the same vignette:

One reason why GL(M)Ms residuals are harder to interpret is that the expected distribution of the data (aka predictive distribution) changes with the fitted values. Reweighting with the expected dispersion, as done in Pearson residuals, or using deviance residuals, helps to some extent, but it does not lead to visually homogenous residuals, even if the model is correctly specified. As a result, standard residual plots, when interpreted in the same way as for linear models, seem to show all kind of problems, such as non-normality, heteroscedasticity, even if the model is correctly specified. Questions on the R mailing lists and forums show that practitioners are regularly confused about whether such patterns in GL(M)M residuals are a problem or not.

Basically your standard residual plots can be severely inaccurate, and I know from personal experience this is the case.

Answer 2

Another good option for posterior predictive checks for binomial predictions is a reliability diagram by Dimitriadis et al. (2021). Here, the mean predictive probabilities (x-axis) are compared to conditional event probabilities (CEP) on the y-axis. CEP are computed by fitting a monotonic step function to the observed events, when they are sorted by predictive probabilities. That is, the diagram is designed to assess the posterior calibration and instead of a even binning, uses a monotonicity assumption.

library(reliabilitydiag)
library(rstanarm)
library(ggplot2)
set.seed(123)

# We fit an example model.
fittedModel <- stan_glm(switch ~ arsenic + dist + educ, data = wells, 
                        chains = 4, iter = 1000, family = "binomial")

# Compute the reliabilitydiagram
r <- reliabilitydiag(x = fittedModel$fitted.values, y = wells$switch)

# Plot the diagram
plot(r)

The resulting plot is a ggplot-object, so we can easily edit the diagram to for example show ggdist:stat_dots instead of the histogram of predictive probabilities:

library(ggdist)

p <- plot(r)

p$layers[[1]] <- stat_dots(aes(x = fittedModel$fitted.values),
                           quantiles = 100, alpha = .5, scale = .5)
p

Now, if the red line of CEP values would fall under confidence interval, the predicted event probabilities are too high, and respectively, if the CEP values are above the intervals, the model is predicting too low probabilities.