Using DHARMa to check Bayesian models fitted with brms

The R package DHARMa is incredibly useful to check many different kinds of statistical models. It can be used with Bayesian models too, although it requires a few more lines of code.

Here I develop an example using DHARMa to check a Bayesian hierarchical generalised linear model fitted with the also fantastic brms package.

Example dataset

First, let’s simulate some count data with a hierarchical structure (random effect) and overdispersion:

library(DHARMa)
testdata <- createData(sampleSize = 1000, 
                       intercept = 0,
                       fixedEffects = 1,
                       numGroups = 20, 
                       family = poisson(), 
                       overdispersion = 0.5)

head(testdata)
##   ID observedResponse Environment1 group time          x         y
## 1  1                9    0.8948431     1    1 0.03217434 0.2598578
## 2  2                4    0.4559604     1    2 0.53322096 0.2235822
## 3  3               18    0.8721283     1    3 0.21682292 0.2619711
## 4  4                0   -0.9403410     1    4 0.07484640 0.9699219
## 5  5                2   -0.3102207     1    5 0.60195011 0.9271296
## 6  6                3    0.8067498     1    6 0.14500979 0.2501301

We have simulated counts (observedResponse) in 20 sites as a function of an environmental variable (Environment1):

library(ggplot2)

ggplot(testdata) +
  geom_point(aes(Environment1, observedResponse)) +
  facet_wrap(~group) 

Fitting the model with brms

First, we try a non-hierarchical Poisson GLM of observed counts as a function of environment. For simplicity, using default priors here.

library(brms)
model <- brm(observedResponse ~ Environment1, 
             data = testdata, 
             family = poisson(),
             cores = 4)
summary(model)
##  Family: poisson 
##   Links: mu = log 
## Formula: observedResponse ~ Environment1 
##    Data: testdata (Number of observations: 1000) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        0.80      0.02     0.75     0.84 1.00     1687     1899
## Environment1     0.98      0.04     0.91     1.06 1.00     1701     1953
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).


Model checking with DHARMa

To use DHARMa with brms (and other Bayesian packages), we must use the createDHARMa function, which takes the posterior simulations and the observed response to generate a DHARMa object, as simulateResiduals do with other supported model types.

model.check <- createDHARMa(
  simulatedResponse = t(posterior_predict(model)),
  observedResponse = testdata$observedResponse,
  fittedPredictedResponse = apply(t(posterior_epred(model)), 1, mean),
  integerResponse = TRUE)

plot(model.check)

We can see there are problems with this model, which we will try to fix below. The posterior predictive check from bayesplot confirms the model is not adequate:

pp_check(model, nsamples = 100) + xlim(0, 20)


Automating DHARMa checks for brms models

To simplify the code, let’s write a function to automate the creation and examination of DHARMa objects from fitted brms models:

check_brms <- function(model,             # brms model
                       integer = FALSE,   # integer response? (TRUE/FALSE)
                       plot = TRUE,       # make plot?
                       ...                # further arguments for DHARMa::plotResiduals 
) {
  
  mdata <- brms::standata(model)
  if (!"Y" %in% names(mdata))
    stop("Cannot extract the required information from this brms model")
  
  dharma.obj <- DHARMa::createDHARMa(
    simulatedResponse = t(brms::posterior_predict(model, nsamples = 1000)),
    observedResponse = mdata$Y, 
    fittedPredictedResponse = apply(
      t(brms::posterior_epred(model, nsamples = 1000, re.form = NA)),
      1,
      mean),
    integerResponse = integer)
  
  if (isTRUE(plot)) {
    plot(dharma.obj, ...)
  }
  
  invisible(dharma.obj)
  
}

Note this function is used here for convenience and will not work with all brms models (e.g. models with more than one response).


Improving the model

Now, let’s try to fix the model. The DHARMa plot above shows that our model is not accounting for potential differences among groups (sites), and there is likely overdispersion (see DHARMa vignette for detailed explanations of how to interpret these plots).

Let’s start adding a varying intercept (random effect) to account for sites/groups:

model2 <- brm(observedResponse ~ Environment1 + (1 | group), 
             data = testdata, 
             family = poisson(), 
             cores = 4,
             iter = 4000)
summary(model2)
##  Family: poisson 
##   Links: mu = log 
## Formula: observedResponse ~ Environment1 + (1 | group) 
##    Data: testdata (Number of observations: 1000) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~group (Number of levels: 20) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.18      0.23     0.84     1.71 1.00      818     1302
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        0.30      0.27    -0.26     0.81 1.00      734     1060
## Environment1     0.95      0.04     0.87     1.02 1.00     3518     4361
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Check with DHARMa using the function we just created:

model2.check <- check_brms(model2, integer = TRUE)

There are still problems, but the residuals show a better pattern now. They are now quite homogeneous among sites:

plot(model2.check, form = testdata$group)

But there are still problems. One of them seems to be overdispersion:

testDispersion(model2.check)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.2343, p-value = 0.008
## alternative hypothesis: two.sided


Improving the model further

Let’s add an observation-level random effect to try to account for overdispersion:

testdata$obs.id <- 1:nrow(testdata)

model3 <- brm(observedResponse ~ Environment1 + (1 | group) + (1 | obs.id), 
             data = testdata, 
             family = poisson(),
             cores = 4, 
             iter = 4000)
summary(model3)
##  Family: poisson 
##   Links: mu = log 
## Formula: observedResponse ~ Environment1 + (1 | group) + (1 | obs.id) 
##    Data: testdata (Number of observations: 1000) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~group (Number of levels: 20) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.18      0.22     0.84     1.68 1.00     1356     2925
## 
## ~obs.id (Number of levels: 1000) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.44      0.03     0.38     0.51 1.00     3190     5113
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        0.20      0.27    -0.34     0.71 1.00      711     1462
## Environment1     0.97      0.05     0.86     1.07 1.00     5896     6340
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Check:

model3.check <- check_brms(model3, integer = TRUE)

There are still some warnings, although the observation-level random effect seems to have corrected for overdispersion:

testDispersion(model3.check)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97885, p-value = 0.808
## alternative hypothesis: two.sided


Improve the model even further

Let’s try using a negative binomial distribution rather than a Poisson:

model4 <- brm(observedResponse ~ Environment1 + (1 | group), 
             data = testdata, 
             family = negbinomial(), 
             cores = 4,
             iter = 4000)

summary(model4)
##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: observedResponse ~ Environment1 + (1 | group) 
##    Data: testdata (Number of observations: 1000) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~group (Number of levels: 20) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.19      0.23     0.83     1.72 1.00      767     1679
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        0.29      0.27    -0.25     0.82 1.01      641     1157
## Environment1     0.96      0.05     0.86     1.06 1.00     4822     4625
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     5.18      0.77     3.90     6.88 1.00     4820     5076
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Check:

model4.check <- check_brms(model4, integer = TRUE)

plot(model4.check, form = testdata$Environment1)

plot(model4.check, form = testdata$group)

The negative binomial model looks better in this case.

Let’s check it with bayesplot too:

pp_check(model4, nsamples = 100) + xlim(0, 20)

library(bayesplot)

ppc_rootogram(y = testdata$observedResponse,
              yrep = posterior_predict(model4, nsamples = 1000)) +
  xlim(0, 20)

Conclusion

DHARMa can greatly facilitate thorough checking and improvement of Bayesian models, as those fitted with brms.

Prof. Florian Hartig kindly revised a draft of this post.

Francisco Rodríguez-Sánchez
Francisco Rodríguez-Sánchez
Researcher

Computational Ecologist & Data Scientist.

Related