I am modelling a dataset from a neuroscience experiment. When I’m trying to incorporate hierarchy in the model, I’m getting a few divergent transitions after warmup, for a variety of distribution families and model definitions. I would appreciate your advice on this; bellow I provide more information about my data and what I tried.
Experimental mice (N=4; total of ~11K trials) from a specific genetic strain were allowed to interact with mice from different genetic strains (
CueID: a factor variable with six levels), and the dependent measure is the time they spent interacting (
DwellTime: a strictly positive and heavily skewed variable). This first figure demonstrates that each mouse has a noticeable number of extreme observations:
This second figure summarizes how
DwellTimelooks in different conditions for each of the four mice.
My first step was to fit a regression model that predicts
DwellTime using the
CueID factor variable without any animal-specific terms:
post_DT_lognormal <- brm(formula = DwellTime ~ CueID, data = big_dframe, family = lognormal(), warmup = 2000, iter = 3000, chains = 4, control = list(adapt_delta = 0.98))
The diagnostics were fine and it seemed like the lognormal distribution did a good job; I tried also
exgaussian() families, but
lognormal() mixed significantly fastest and had better loo log-likelihood.
When I introduced the next level of complexity, i.e., animal-specific intercepts:
post_DT_lognormal_mixeff_int_only <- brm(formula = DwellTime ~ CueID + (1 | AnimalID), data = big_dframe, family = lognormal(), warmup = 2000, iter = 3000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 15))
I encountered a few divergent transitions after warmup (<10 across different attempts). Note that the model with the animal-specific intercepts had a better loo log-likelihood as compared to a population-effects only model.
LOOIC SE post_DT_lognormal_mixeff_int_only 69595.29 312.75 post_DT_lognormal 70063.17 304.59 post_DT_lognormal_mixeff_int_only - post_DT_lognormal -467.88 45.17
When looking at
pairs(post_DT_lognormal_mixeff_int_only) () we see that the
sd_AnimalID_Intercept term introduces a complex geometry with the other parameters (second column from right). If I understand correctly, brms implements a non-centered parametrization so this is probably not the problem.
And this isn’t solved when I try the other families that I mentioned above; It persists even when the model is in the simplified hierarchical form:
formula = DwellTime ~ 1 + (1 | AnimalID), and moreover, when I trimmed outlier observations as well (and we wouldn’t like to throw these observations in any case).
Note also that when I log the data, it doesn’t look like a classic normal distribution, so not sure whether the log-normal distribution is the best here.
Similarly and as expected, when I introduce even more complexity and add animal-specific
post_DT_lognormal_mixeff <- brm(formula = DwellTime ~ CueID + (CueID || AnimalID), data = big_dframe, family = lognormal(), warmup = 2000, iter = 3000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 15))
I get 3 divergent transitions after warmup and similar geometric issues. Again, the model does a better job than the varying-intercepts model and ideally I would prefer it.
LOOIC SE post_DT_lognormal_mixeff 69460.41 313.12 post_DT_lognormal_mixeff_int_only 69595.29 312.75 post_DT_lognormal_mixeff - post_DT_lognormal_mixeff_int_only -134.88 22.25
Do you have any suggestions?
- Operating System: macOS Mojave 10.14.2,
- brms Version: 2.01