Hi there,
I have been trying to fit a hierarchical binomial model on the lottery/surebet choices by rats given the lottery offer. I started with a simple nonhierachical version and synthetic data, it works just fine. I then progressed to a specieslevel model. Interestingly, this model works well when N_SUBJ < 8
, meaning all the 4 chains with default parameters run and give no divergences. However, when I set N_SUBJ >= 8
in the synthetic dataset, only 1 out of 4 chains would run and the others gave this typical error:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Initialization between (2, 2) failed after 100 attempts.
And the remaining one chain that did initialise would often end with 100% divergences.
When I ramp up N_SUBJ
to be a large value, say 20, then none of the chains would initialise…
Note that previously I have completed a bernoulli_logit
version of the same model and it works perfectly fine on all levels, but I want to speed up computation by using binoimal
if possible due to large data size.
Here’s the model, any help is appreciated!
// hierarchical model fitting one species as betarho risky agents with binomial
functions {
real choiceprob(real rho, real beta, real lott_value, real lott_prob, real sb_value, real rew_multi) {
real y;
real u1; // Lottery utility
real u2; // Surebet utility
u1 = lott_prob * (rew_multi * lott_value) ^ rho;
u2 = (rew_multi * sb_value) ^ rho;
y = 1 / (1 + exp(beta * (u1  u2)));
return y;
}
}
data {
int<lower=0> T; // Number of trials types
int<lower=0> K; // number of subjects in each species
int individual[T]; // vector of subjid indexes
vector[T] lott_value; // Lottery value for that choice
vector[T] lott_prob; // Lottery probabilities for that choice
vector[T] sb_value; // Surebet values for that choice
vector[T] total_rew_multi; // total reward multiplier = base_reward * rew_multi
int n_chose_lott[T]; // number chose lottery for this trial type
int n_trials[T]; // total number of trials for this trial type
}
parameters {
real log_rho_s; // specieslevel rho
real log_beta_s; // specieslevel beta
real<lower=0> log_sigma_rho; // standard deviation for species rho
real<lower=0> log_sigma_beta; // standard deviation for species beta
vector[K] log_rho_i; // individuallevel rho
vector[K] log_beta_i; // individuallevel beta
}
transformed parameters{
real rho_s;
real beta_s;
vector[K] rho_i;
vector[K] beta_i;
rho_s = exp(log_rho_s);
beta_s = exp(log_beta_s);
rho_i = exp(log_rho_i);
beta_i = exp(log_beta_i);
}
model {
// set weak priors
log_rho_s ~ std_normal();
log_beta_s ~ std_normal();
// draw individual parameters from the species distribution
log_rho_i ~ normal(log_rho_s, log_sigma_rho);
log_beta_i ~ normal(log_beta_s, log_sigma_beta);
// fit the actual model to each trial
for(t in 1:T){
n_chose_lott[t] ~ binomial(n_trials[t], choiceprob(rho_i[individual[t]], beta_i[individual[t]], lott_value[t], lott_prob[t], sb_value[t], total_rew_multi[t]));
}
}
And some example rows of synthetic data file, note that each individual / subject only has 6 rows of data here (6 different lottery offers).
individual lott_value lott_prob sb_value total_rew_multi n_chose_lott n_trials
1 1 0 0.5 3 1 34 327
2 1 2 0.5 3 1 65 327
3 1 4 0.5 3 1 135 332
4 1 8 0.5 3 1 192 302
5 1 16 0.5 3 1 344 352
6 1 32 0.5 3 1 359 359
7 2 0 0.5 3 1 70 342
8 2 2 0.5 3 1 83 307
9 2 4 0.5 3 1 129 345
10 2 8 0.5 3 1 192 360
11 2 16 0.5 3 1 262 336
12 2 32 0.5 3 1 293 309
13 3 0 0.5 3 1 83 341
Some things I have tried but did not work:
 By using the
shinystan
interface I found out that when the chains have high divergences, thestepsize
is often a ridiculously small number compared to whatstepsize
is when the model runs fine (~0.35). So I addedstepsize = 0.35
in setting, did not help.  Increasing
adapt_delta
even to 0.99 did not help at all.  Increasing both
adapt_delta
and setstepsize
did not help.
Something even more interesting:

N_SUBJ = 8
:only 1 chain runs 
N_SUBJ = 9
: 2 chains run 
N_SUBJ = 10
: 0 chains run 
N_SUBJ = 11
: 1 chain run 
N_SUBJ > 11
: 0 chains run