Bayesian Mixed-Model – How to Compare DHARMa Residuals Plot vs. Binned Residuals Using stan_glm Object

bayesianlogisticmixed modelresidualsseparation

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:

residuals(simulationOutput)
plot(simulationOutput)

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

Best Answer

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

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.

Related Question