Can't shake last few divergences in brms multilevel model for modest dataset


Hi there,

I’ve been working on fitting some multilevel models using brms and I can’t seem to shake the last few divergences, despite my best efforts. I’m trying to model the effects of two particular treatments on aphid population growth over time. I’m trying to predict aphid population size at t=2 using the treatment and the population size at t=1, with varying effects across 3 trial blocks (30 total reps per block).

I’ve raised adapt_delta to > .9999 and max_treedepth to 30, tweaked my priors a bunch (I’ve tried Student’s t, cauchy, and normal priors with lots of different widths), and I still get 2-8 diversions each time. The model still runs pretty quickly on my machine, even with the high control values. Since brms uses the non-centered parameterization by default, I don’t think there’s anything I can do on that front.

I’ve examined the pairs plots and I can’t seem to find any really strong patterns for the divergent iterations, but I’ve included the plot here in case it’s helpful.

Here is the model as I’ve written it in R:

model_3.1 <- brm(
  data = data,
  family = negbinomial,
  aphid_total_count_2 ~ 1 + treatment + aphid_total_count_1 +
    (1 + treatment + aphid_total_count_1 |
  prior = c(
    set_prior("normal(3,1)", class = "Intercept"),
    set_prior("normal(0.5,5)", class = "b", coef = "aphid_total_count_1"),
    set_prior("normal(0,.5)", class = "b"),
    set_prior("cauchy(0,.1)", class = "sd"),
    set_prior("lkj(1)", class = "cor")
  iter = 5000,
  warmup = 1000,
  chains = 3,
  cores = future::availableCores(),
  control = list(adapt_delta = 0.9999999999, max_treedepth = 40)

Here are the data and the little bit of code to read and clean the data for use in the model:

APHIS_Fall_2018_Trials_2_3_4.csv (2.5 KB)

data <- read_csv("raw_data/APHIS_Fall_2018_Trials_2_3_4.csv", skip = 1)
data$trial_number <- as.factor(data$trial_number)

I’ve searched high and low on discourse, the google group, the manual, etc. and I’ve found lots of useful information, but nothing quite seems to work. Unfortunately I’m still learning to write straight Stan, so some of the advice was a bit over my head. I’m hoping perhaps I’m missing something obvious, and I’d really like to improve my application of these types of multilevel models in general. Thank you!


Kudos to your efforts of trying to get rid of divergent transitions!

I see muliple things that may explain your problem:

  1. You may want to bring all predictors on the same scale. That is I would standardize aphid_total_count_1 before specifying the model. After doing this, I got rid of all divergent transitions:
data$z_aphid_total_count_1 <- scale(data$aphid_total_count_1)

model_3.2 <- brm(
  data = data,
  family = negbinomial,
  aphid_total_count_2 ~ 1 + treatment + z_aphid_total_count_1 +
    (1 + treatment + z_aphid_total_count_1 | trial_number),
  prior = c(
    set_prior("normal(0,3)", class = "Intercept"),
    set_prior("normal(0.5,5)", class = "b", coef = "z_aphid_total_count_1"),
    set_prior("normal(0,.5)", class = "b"),
    set_prior("normal(0,.5)", class = "sd"),
    set_prior("lkj(2)", class = "cor")
  iter = 2000,
  warmup = 1000,
  chains = 2,
  cores = 2,
  control = list(adapt_delta = 0.99, max_treedepth = 10)

Some additional thoughts:

  1. Fitting varying effects with only 3 levels may not induce much hierarchical shrinkage and may
    lead to convergence problems in some cases.
  2. I would prefer half-normal or half-student-t over cauchy for SDS especially because your SDs are only identified through 3 levels.
  3. Estimating the correlations between varying effects won’t be very helpful for only 3 levels. Imagine a correlation being estimated on a sample size of 3.


Hi Paul,

Thank you so much for the helpful feedback! I was able to scale the aphid_total_count_1 and shake the last of the divergences as well.

I also really appreciate the additional thoughts, I’ll take them into consideration as I move forward with this project!