Issues with fitting this poisson-regression

Modeling
,
1

Hi all,

I am trying to implement and copy the model in the example ‘1.5 Ache Hunting with Age Trends’ from the book Bayesian Ideas and Data Analysis (Christensen et al., 2010). I have attached the relevant example pages from the book and the dataset used.

1.5 Ache Hunting with Age Trends.pdf (2.3 MB)

Exercise 1.5.csv (528 Bytes)

Some of the warning messages I am getting is below:

  • There were 2963 divergent transitions after warmup.
  • There were 2381 transitions after warmup that exceeded the maximum treedepth.
  • There were 1 chains where the estimated Bayesian Fraction of Missing Information was low.

The rstan code to fit the model is:

stan(file=‘Exercise 1.5.stan’, data=data, chains=4, iter=6000, warmup=2000, thin=1, cores=4,
control=list(adapt_delta=0.99, max_treedepth=15))

I suspect that the issue may be relating to my Stan codes. Can someone please explain to me what I might be doing wrong? And what I can do to fix this?

The Stan codes are:

data {
  int<lower=1> N;
  vector<lower=0>[N] age;
  int<lower=0> kills[N];
}

transformed data {
  real mu_age = mean(age);
  vector[N] ones = rep_vector(1, N);
  vector[N] age_centred = age - mu_age;
  vector[N] age_centred_squared = square(age_centred);
  matrix[N, 3] X;
  X = append_col(ones, append_col(age_centred, age_centred_squared));
}

parameters {
  vector[3] beta;
  real error;
  real<lower=0> sigma_squared;
}

transformed parameters {
  real<lower=0> sigma = sqrt(sigma_squared);
}

model {
  // priors
  beta ~ normal(0, 1000);
  sigma_squared ~ inv_gamma(0.001, 0.001);
  error ~ normal(0, sigma);
  
  // likelihood
  kills ~ poisson_log(X*beta + error);
}

Many thanks.

2

Isn’t the first element in “beta” serving the same role as the “error” parameter?

3

Hi Mike,

I am trying to replicate the model as per that example. The beta1 and error/delta parameter are from different normal distributions.

4

Ah, having looked at the paper now, “error” needs to be a vector of length N

2 Likes
5

Thanks Mike! The model has converged properly now :)

1 Like
6

Actually, I’m still getting similar convergence issues I had before when I re-run the updated code where “error” is now a vector. Sometimes there are no issues, and other times there are. What could be the reason for this? Is it to do with initial values?

Also, when the model converges, for the parameter means I am getting \beta_1=0.23, \beta_2=0.11, \beta_3=0.00. This happens to be very different to means in the attached paper, where the \beta_1=-0.7147, \beta_2=0.01368, \beta_3=0.002683. Is anyone able to replicate the results from the paper?

My updated code is:

data {
  int<lower=1> N;
  vector<lower=0>[N] age;
  int<lower=0> kills[N];
}

transformed data {
  real mu_age = mean(age);
  vector[N] ones = rep_vector(1, N);
  vector[N] age_centred = age - mu_age;
  vector[N] age_centred_squared = square(age_centred);
  matrix[N, 3] X;
  X = append_col(ones, append_col(age_centred, age_centred_squared));
}

parameters {
  vector[3] beta;
  vector[N] error;
  real<lower=0> sigma_squared;
}

transformed parameters {
  real<lower=0> sigma = sqrt(sigma_squared);
}

model {
  // priors
  beta ~ normal(0, 1000);
  sigma_squared ~ inv_gamma(0.001, 0.001);
  error ~ normal(0, sigma);
  
  // likelihood
  kills ~ poisson_log(X*beta + error);
}