1
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ Why have you put such uninformative priors on this model? Have you simulated from the prior predictive to see what that implies? Do you think that is a reasonable prior for the data? $\endgroup$ Mar 4, 2021 at 21:54
  • $\begingroup$ I am trying to do a homework and trying to throw in a uniform prior seems like a quick and dirty way to get it done :) $\endgroup$ Mar 4, 2021 at 22:06

1 Answer 1

1
$\begingroup$

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.

enter image description here

$\endgroup$
4
  • $\begingroup$ Can you please tell what was the rationale behind choosing a Cauchy prior? How does that help reduce the divergences? $\endgroup$ Mar 5, 2021 at 0:22
  • 1
    $\begingroup$ @desert_ranger Long and the short of it is that it allows the parameter to be big if it needs to be, but has most of the mass near 0. It isn't a good default for everything, but it appears to be better than the uniform prior in a lot of cases. $\endgroup$ Mar 5, 2021 at 0:30
  • $\begingroup$ Thank you. Could you please elaborate just a little more. For instance, in a uniform prior, I would be assuming all values between a certain range. What exactly does Cauchy do in this case? Feel free to suggest some reading material :) $\endgroup$ Mar 5, 2021 at 0:32
  • $\begingroup$ Never mind! I got it. Thank you once again. $\endgroup$ Mar 5, 2021 at 0:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.