florianhartig / DHARMa

Diagnostics for HierArchical Regession Models
http://florianhartig.github.io/DHARMa/
214 stars 22 forks source link

Check compatibility with brms #33

Open florianhartig opened 7 years ago

florianhartig commented 7 years ago

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 https://github.com/paul-buerkner/brms/issues/281

simschul commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 3 years ago

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

florianhartig commented 3 years ago

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

StaffanBetner commented 3 years ago

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.

melina-leite commented 2 weeks ago

DHARMa.helpers, github page: https://pakillo.github.io/DHARMa.helpers/ https://frodriguezsanchez.net/post/using-dharma-to-check-bayesian-models-fitted-with-brms/