Good R-hat, good n_eff but large se_mean and sd!

Modeling
, ,
1

Hi everyone,

I am trying to fit an SEIR model to influenza data. In the data, I have the incidence of newly reported cases per week. In my Stan model, I simulated the SEIR model per day and then aggregated the data per week for the incidence and fitted it to the data.

My priors are the transmissibility \beta, the initial number of exposed people to the influenza e_{0}, the initial number of infected people i_{0}, and the initial number of recovered people from the previous influenza season r_{0}. I have considered a normal prior for transmissibility \beta and flat priors for the initial conditions, which is a uniform distribution between 0 and an upper limit, which is the population of the country that I am using for the influenza data.

However, when I run the model, R-hat is almost 1 (1.01) and n_{eff} is more than 400, but with 4000 iterations, I get 103 divergent transitions and also low Bulk and Tail Effective Sample Sizes (ESS). Additionally, the se_mean and standard deviation (sd) are very large.

Why does this happen? Why are R-hat and n_{eff} good, but se_mean and sd are very large? Below, I have included my Stan model. I really appreciate it if anyone helps me.

functions {
real sir(real t, real y, real theta,
real x_r, int x_i) {

  real mu = x_r[1];
  real epsilon = x_r[2];
  real N = x_i[1];


  
  real beta = theta[1];
  real e0 = theta[2];
  real i0 = theta[3];
  real r0 = theta[4];
  
  real init[4] = {N - i0 - e0 - r0, e0, i0, r0}; // initial values
  real S = y[1] + init[1];
  real E = y[2] + init[2];
  real I = y[3] + init[3];
  real R = y[4] + init[4];
  
  real dS_dt = -beta * I * S / N;
  real dE_dt =  beta * I * S / N - epsilon * E;
  real dI_dt = epsilon * E - mu * I;
  real dR_dt =  mu * I;
  
  return {dS_dt, dE_dt, dI_dt, dR_dt};

}
}
data {
int<lower=1> n_days;
real t0;
real ts[n_days];
real mu;
real epsilon;
int N;
int cases[n_days/7];
}
transformed data {
real x_r[2]={mu, epsilon};
int x_i[1] = { N };
}
parameters {
real<lower=0> beta;
real<lower=0> e0;
real<lower=0> i0;
real<lower=0> r0;
real<lower=0> phi_inv;
real<lower=0, upper=1> p_reported; // proportion of infected (symptomatic) people reported
}
transformed parameters{
real y[n_days, 4];
real incidence[n_days - 1];
real aggregated_incidence[n_days / 7]; // New aggregated incidence array
real phi = 1. / phi_inv;

real theta[4] = {beta, e0, i0, r0};

y = integrate_ode_rk45(sir, rep_array(0.0, 4), t0, ts, theta, x_r, x_i);

for (i in 1:n_days-1){
incidence[i] = -(y[i+1, 2] - y[i, 2] + y[i+1, 1] - y[i, 1]) * p_reported; //-(E(t+1) - E(t) + S(t+1) - S(t))
}
// Aggregate incidence for every seven days
for (j in 1:(n_days / 7)) {
int start_index = (j - 1) * 7 + 1;
int end_index = min(j * 7, n_days-1);
aggregated_incidence[j] = sum(incidence[start_index:end_index]);
}
}
model {
//priors
beta ~ normal(2, 1);
real upper_limit = N; // Upper limit as 1% of the population
e0 ~ uniform(0, upper_limit); // Flat prior for e0
i0 ~ uniform(0, upper_limit); // Flat prior for i0
r0 ~ uniform(0, upper_limit); // Flat prior for r0
phi_inv ~ exponential(5);
p_reported ~ beta(1, 2);

//sampling distribution
//col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
cases[1:(n_days/7)] ~ neg_binomial_2(aggregated_incidence, phi);
}
generated quantities {
real R0 = beta / mu; // Calculate R0 using fixed mu value
real pred_cases[n_days/7];
pred_cases = neg_binomial_2_rng(aggregated_incidence, phi);
}