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);

}