# Good R-hat, good n\_eff but large se\_mean and sd!

**URL:** <https://discourse.mc-stan.org/t/good-r-hat-good-n-eff-but-large-se-mean-and-sd/35274>\
**Category:** Modeling\
**Tags:** stan, fitting-issues, rstan\
**Created:** [May 27, 2024, 2:32pm UTC](https://discourse.mc-stan.org/t/good-r-hat-good-n-eff-but-large-se-mean-and-sd/35274 "2024-05-27T14:32:35Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![rmrmasoomi](https://avatars.discourse-cdn.com/v4/letter/r/47e85d/32.png) [@rmrmasoomi](https://discourse.mc-stan.org/u/rmrmasoomi)\
**Post date:** [May 27, 2024, 2:32pm UTC](https://discourse.mc-stan.org/t/good-r-hat-good-n-eff-but-large-se-mean-and-sd/35274/1 "2024-05-27T14:32:35Z")

</div>

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