brms_dharma.R

#pacman::p_load(pacman, tidyverse, rio, brms, DHARMa)

# inspired by https://frodriguezsanchez.net/post/using-dharma-to-check-bayesian-models-fitted-with-brms/


# example code commented out, assumes that you have DHARMa and brms installed.
# posterior_epred is more likely a better choice than posterior_linpred, 
# but the former does not work with ordinal models
# due to being a multivariate model in disguise (which shares linear predictor)

# example data
#import("https://stats.idre.ucla.edu/stat/data/ologit.dta", setclass="tbl") %>% 
#  mutate(apply = apply %>% factor(labels=c("unlikely", "somewhat likely", "very likely"), ordered = TRUE),
#         across(c(pared, public), ~.x %>% factor(labels=c("no","yes")))) ->
#  ordinal_dat

# function
model.check <- function(model) {
  require(DHARMa)
  response_variable <- all.vars(model$formula$formula)[1] # extract name of response variable
  observedResponse <- as.numeric(model$data[[response_variable]]) # convert response into numeric
  createDHARMa(
    simulatedResponse = t(posterior_predict(model)),
    observedResponse = observedResponse,
    fittedPredictedResponse = apply(t(posterior_linpred(model)), 1, mean), 
    integerResponse = !any(as.numeric(model$data[[response_variable]])%%1!=0) # checks if response consists only of integers
  )
}

# it is a good idea to memoise model.check
model.check <- memoise::memoise(model.check)

# example model for ordinal data
#brm(
#  formula = apply ~ pared*public,
#  data = ordinal_dat,
#  family = cumulative(),
#  cores = 4,
#  backend = "cmdstanr"
#  ) ->
#model


#model %>%
#  model.check() %>%
#  plot()