Check compatibility with brms

[florianhartig]

https://cran.r-project.org/web/packages/brms/index.html

See also https://github.com/florianhartig/DHARMa/wiki/Adding-new-R-packages-to-DHARMA

See also paul-buerkner/brms#281

[simschul]

Are there any updates on that?

I am rather new to both multilevel modeling and Bayesian stats, thus my approach might be too naive or just wrong:

library(brms)
library(DHARMa)

# example from brms
bprior1 <- prior(student_t(5,0,10), class = b) +
  prior(cauchy(0,2), class = sd)
fit1 <- brm(count ~ zAge + zBase * Trt + (1|patient),
            data = epilepsy, family = poisson(), prior = bprior1)

# sample from the Posterior Predictive Distribution
preds <- posterior_predict(fit1, nsamples = 250, summary = FALSE)
preds <- t(preds)

res <- createDHARMa(simulatedResponse = preds, 
                    fittedPredictedResponse = apply(preds, 1, median), 
                    observedResponse = epilepsy$count, 
                    integerResponse = "poisson")

plot(res, quantreg = FALSE)

Happy to hear if that makes sense to you!

[florianhartig]

Hi Simon,

without looking at the details, this makes sense to me. Only interResonse should be T/F - the only reason that this doesn't throw an error is that I recently changed DHARMa to a new residual definition, which doesn't use the integerResponse info any more.

In general, my intention for this ticket was more to support brm like I do for the other packages, so that you could do simulateResiduals(fit1). In your example, this looks definitely doable, but I don't know how generalizeable this is, i.e. if users can define models where this would break.

A tricky thing is definitely how to extract the response from the model. You use the variable directly, but this wouldn't work for me. I suppose that given posterior_predict knows that the response is, it should be possible to extract this from fit1 automatically as well, but I'm not sure how fast this would run into problems. I think I will probably have to ask the brm developers, because I'm not using brm enough myself to have a good overview.

[simschul]

Thanks for your response!

The function get_response from the insight-package extracts the response from (brm) models. However, brms also support models with more than one response variable which makes things more complicated...

[florianhartig]

Hi Simon,

hmm ... yes, tricky. At the moment, I'm tending to think that brms is probably too flexible to create a reliable out-of-the-box support, but I suppose I should have another look, and the least I could do is to add some examples along the lines of what you do in the vignette.

Thanks for reviving this issue, I'll put this on the list for the next release!

[florianhartig]

https://twitter.com/StaffanBetner/status/1404071597837795335

[florianhartig]

https://gist.github.com/StaffanBetner/6241c5b816a8a78b20bd300743b3a66d

[StaffanBetner]

Regarding brms, don't forget this blog post. https://frodriguezsanchez.net/post/using-dharma-to-check-bayesian-models-fitted-with-brms/

The latest link you added is for the mgcv::ocat family.