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

**URL:** <https://discourse.mc-stan.org/t/cant-shake-last-few-divergences-in-brms-multilevel-model-for-modest-dataset/6781>\
**Category:** brms\
**Created:** [November 29, 2018, 12:50am UTC](https://discourse.mc-stan.org/t/cant-shake-last-few-divergences-in-brms-multilevel-model-for-modest-dataset/6781 "2018-11-29T00:50:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Michael\_C-M](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/michael_c-m/32/5900_2.png) [@Michael\_C-M](https://discourse.mc-stan.org/u/Michael_C-M)\
**Post date:** [November 29, 2018, 12:50am UTC](https://discourse.mc-stan.org/t/cant-shake-last-few-divergences-in-brms-multilevel-model-for-modest-dataset/6781/1 "2018-11-29T00:50:34Z")

</div>

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.

 ![mcmc_pairs_plot_negbin_multi](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/02f5da777d933f944d7fd9c7d55f4a2f22c0f54b.jpeg)

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 |
       trial_number),
  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](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/089ef6c78d6e1c6230b197b0170e05f07984f63c.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!

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [November 29, 2018, 9:06am UTC](https://discourse.mc-stan.org/t/cant-shake-last-few-divergences-in-brms-multilevel-model-for-modest-dataset/6781/2 "2018-11-29T09:06:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Michael\_C-M](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/michael_c-m/32/5900_2.png) [@Michael\_C-M](https://discourse.mc-stan.org/u/Michael_C-M)\
**Post date:** [December 7, 2018, 12:26am UTC](https://discourse.mc-stan.org/t/cant-shake-last-few-divergences-in-brms-multilevel-model-for-modest-dataset/6781/3 "2018-12-07T00:26:20Z")

</div>

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!