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.

enter image description here

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.

Binned residuals plot produced with performance::binned_residuals

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:


And the plots came out perfect.

Residuals plot using DHARMa

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

enter image description here

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.

Boxplot returned by DHARMa

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

Usual-looking scatter plot returned by DHARMa

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


2 Answers 2


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:

enter image description here

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

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

#### 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:

enter image description here

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.

  • $\begingroup$ Hi Shawn, thank you for your response and the extra miles you went for labelling the steps for me. I followed the codes and produced a boxplot (if using median) and a scatterplot (if using mean), please see the update in the main post. Can I use either one? And do they look ok? (please excuse me for these rather newb questions). $\endgroup$ Commented Mar 4, 2023 at 17:20

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.


# 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

enter image description here

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:


p <- plot(r)

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

enter image description here

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.

  • $\begingroup$ Instead of an equal binning, as in the default histogram, the dots plot shows a selected amount of dots depicting the quantiles of predicted values. In my experience, dots plots are good at showing tails, outliers, and possible discontinuities. $\endgroup$ Commented Mar 6, 2023 at 12:13
  • $\begingroup$ Thank you for the suggestion! I have tried it and things come out perfect, too! I have a reputation too low to give you an upvote; but know that I thank you really! $\endgroup$ Commented Mar 6, 2023 at 17:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.