Eliminating divergent transitions in Stan

Question

I have the following dataset -

d = data.frame(year=1976:1985,
               fatal_accidents = c(24,25,31,31,22,21,26,20,16,22),
               passenger_deaths = c(734,516,754,877,814,362,764,809,223,1066),
               death_rate = c(0.19,0.12,0.15,0.16,0.14,0.06,0.13,0.13,0.03,0.15))
d$miles_flown = d$passenger_deaths/d$death_rate

I'd like to hierarchically model the passenger deaths($Y_{i}$) using a Poisson distribution with the parameter consisting of miles_flown ($x_{i}$), i.e - $Y_{i} = Poisson(x_{i}\lambda_{i})$

Given the above information, my mid-level model follows a Gamma distribution ($Gamma(\alpha, \beta)$), since it is conjugate to the Poisson.

In order to model, my $\alpha$ and $\beta$, I decided to use a uniform distribution, i.e - $\alpha \sim Uniform(0, a_{0})$ $\beta \sim Uniform(0, b_{0})$,

wherein I decided to manually pick $a_{0}, b_{0}$ values. Here is my Stan code -

library("rstan")
set.seed(8889)

model_default_prior = "
data {
int<lower=0> N;
vector[N] miles_flown;
int<lower=0> fatal_accidents[N];
}
parameters {
real<lower=0> lambda[N];
real<lower=0> alpha;
real<lower=0> beta;
}

model {
// Uninformative prior
alpha ~ uniform(0, 100);
beta ~ uniform(0, 100);

// implicit joint distributions
lambda ~ gamma(alpha,beta);
fatal_accidents~poisson(lambda);
}
"

d.dat = list(fatal_accidents = d$fatal_accidents, miles_flown = d$miles_flown, N = nrow(d))
m = stan_model(model_code=model_default_prior)
r.d = sampling(m, d.dat, c("alpha","beta","lambda"), iter=10000, control = list(adapt_delta = 0.9))
r.d

However, after running the above code, I am running into 14372 divergent transitions. This seems like I have made some serious error with my code or model.

Answer 1

The model has several issues.

The first is that your priors for alpha and beta are far far too wide. The alpha and beta priors you have specified lead to a lambda which is concentrated at 0 and often goes as large as 15,000. Try generating draws from your prior and see for yourself.

Second, the model you have written is not the model you are estimating. In your model, each $y$ has its own $\lambda$ which is independent of the number of miles flown. I've written the model up in Stan myself. I've placed hald cauchy priors on alpha and beta. They are still non-informative, but are more preferable to uniform priors (additionally, there is good geometric reason not to put uniform priors on unbounded quantites like alpha and beta, which I just don't have the time to get into at the moment).

model_code = '
data{
  int N;
  vector[N] x;
  int y[N];
}
parameters{
  real<lower=0> lam;
  real<lower=0> a;
  real<lower=0> b;
}
model{
  y ~ poisson(lam * x);
  lam ~ gamma(a, b);
  a ~ cauchy(0,1);
  b ~ cauchy(0,1);
}
generated quantities{
  real yppc[N] = poisson_rng(lam*x);
}
'

When I fit this model I get seven divergences, which is a sign it is a bad model and in need of additional prior information. For the time being, these divergences can be resolved by increasing the adapt_delta argument to 0.99 (again, it isn't preferable to make the sampler take such small steps, but it gets rid of the divergences).

With these changes, you can now do model checking. A posterior predictive check confirms this is indeed a poor model to fit to the data. The posterior predictive distribution does not look much like the original data.