Background: We’ve been building up a gamma mixture model with censoring to help understand some interval data. Initially, our fit had serious problems (multiple modes, lots of divergences), but a combination of changes to the model and better priors seems to have fixed most of them. Unfortunately, we’re still getting a divergent transition about every 5000 steps or so. Also, sometimes one chain lags behind the other ones, taking several times longer to run (we think that this is related to the same bad geometry of the model that’s causing divergences, however). We have not finished adding parameters to this model but would like to fix the divergent transition before continuing further.
Attempts to fix issue: In addition to tweaking our priors, we’ve tried raising adapt_delta as high as 0.99, and made sure that both gamma distributions had a log link, as suggested in response to an earlier question. We’ve also tried looking at the divergences in various pairs-plots and in shinystan, but didn’t see any clear funnels, so we’re not really sure which part of the model is at fault.
I’ve included a copy of our model code below, along with a data frame. Do you see anything we could do to eliminate the bad geometry that’s causing the last few divergences (and the occasional super-slow chain)? Thanks in advance for any help you can provide!
Relevant data frame:
brms-data.csv (30.2 KB)
intervals <- read_csv(‘brms-data.csv’)
set.seed(5) # With this seed, we get one divergence
chains = 4
prior1 <- c(
set_prior("gamma(50, 10)", class="shape1"),
set_prior("gamma(100, 50)", class="shape2"),
set_prior("student_t(10, 2, .5)", class = "Intercept", dpar = "mu1"),
set_prior("student_t(10, 5, 0.5)", class = "Intercept", dpar = "mu2"),
set_prior("student_t(10, -1, 1)", class = "Intercept", dpar = "theta2"),
set_prior("student_t(10, 0, .5)", class = "b", dpar = "mu1"),
set_prior("student_t(10, 0, .5)", class = "b", dpar = "mu2"),
set_prior("student_t(10, 0, 1)", class = "b", dpar = "theta2"),
set_prior("lkj(2)", class = "cor"),
set_prior("lognormal(-2, .25)", class = "sd", dpar = "mu1"),
set_prior("lognormal(-2, .25)", class = "sd", dpar = "mu2"),
set_prior("lognormal(-2, .25)", class = "sd", dpar = "theta2")
)
fit1 <- brm(
brmsformula(
intervals | cens(censored) ~ 1 + phase + (1|male|male_id) + (1|trial),
theta2 ~ 1 + phase + (1|male|male_id) + (1|trial)
),
data = x,
family = mixture(
Gamma("log"),
Gamma("log"),
# Enforces that mixture elements (short vs long)
# have the same interpretation across chains
order = TRUE
),
prior = prior1,
control = list(adapt_delta = 0.95),
chains = chains,
cores = chains,
inits = "random",
sample_prior = FALSE,
init_r = 0.25,
sample_file = "chains",
iter = 2000
)