Rejecting Initial values, ESS too low, etc

Modeling
1

Good day house, I am trying to follow the code lines in “Bayesian workflow for disease transmission modeling in Stan” to estimate the parameter values of SIR model on covid for a particular country but I have been getting loads of errors. I am new to RStan. Help Pls. Below are Stan and R codes, alongside the error messages.
Stan code

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

      real S = y[1];
      real I = y[2];
      real R = y[3];
      int N = x_i[1];
      
      real beta = theta[1];
      real gamma = theta[2];
      
      real dS_dt = -beta * I * S / N;
      real dI_dt = beta * I * S / N - gamma * I;
      real dR_dt =  gamma * I;
      
      return {dS_dt, dI_dt, dR_dt};
  }
}
data {
  int<lower=1> n_days;
  real t0;
  real y0[3];
  real ts[n_days];
  int N;
  int cases[n_days];
}
transformed data {
  real x_r[0];
  int x_i[1] = { N };
}
parameters {
  real<lower=0> gamma;
  real<lower=0> beta;
  real<lower=0> phi_inv;
  real<lower=0, upper=1> p_reported; // proportion of infected (symptomatic) people reported
}
transformed parameters{
  real y[n_days, 3];
  real incidence[n_days - 1];
  real phi = 1. / phi_inv;
  // initial compartement values
  real theta[2] = {beta, gamma};
  y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
  for (i in 1:n_days-1){
    incidence[i] =  (y[i, 1] - y[i+1, 1]) * p_reported;
  }
}
model {
  //priors
  beta ~ normal(2, 1);
  gamma ~ normal(0.4, 0.5);
  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-1)] ~ neg_binomial_2(incidence, phi);
}
generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  real pred_cases[n_days-1];
  pred_cases = neg_binomial_2_rng(incidence, phi);
}

R codes

setwd("~/R")

nga <- read.csv("covid_data.csv", header=TRUE, sep=",")

#rename(date = 'Date_reported') 
nga %>%
ggplot() + 
  geom_bar(mapping = aes(x = Date_reported, y = New_cases), stat = "identity") +
  labs(y="Number of reported cases")

# Nigeria population
N <- 211400708;

# time series of cases
cases <- nga$New_cases  # Number of students in bed

# times
n_days <- length(cases) 
t <- seq(0, n_days, by = 1)
t0 = 0 
t <- t[-1]

#initial conditions
i0 <- 5
s0 <- N - i0
r0 <- 0
y0 = c(S = s0, I = i0, R = r0)

# data for Stan
data_sir <- list(n_days = n_days, y0 = y0, t0 = t0, ts = t, N = N, cases = cases)

# number of MCMC steps
niter <- 2000

model <- stan_model("seir.stan")


#pass the data to stan and run the model
fit_seir <- sampling(model,
                           data = data_sir,
                           iter = niter,
                           chains = 4, 
                           seed = 0)
pars=c('beta', 'gamma', "R0", "recovery_time")

print(fit_sir_negbin, pars = pars)

stan_dens(fit_sir_negbin, pars = pars, separate_chains = TRUE)

Error messages
SAMPLING FOR MODEL ‘seir’ NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: neg_binomial_2_lpmf: Location parameter[184] is -9.82331e-09, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: neg_binomial_2_lpmf: Location parameter[17] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: neg_binomial_2_lpmf: Location parameter[27] is -3.54775e-08, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: neg_binomial_2_lpmf: Location parameter[202] is -2.4914e-10, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 1:
Chain 1: Gradient evaluation took 0.003 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 30 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 625.273 seconds (Warm-up)
Chain 1: 901.8 seconds (Sampling)
Chain 1: 1527.07 seconds (Total)
Chain 1:

SAMPLING FOR MODEL ‘seir’ NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.002 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 20 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 22.021 seconds (Warm-up)
Chain 2: 19.652 seconds (Sampling)
Chain 2: 41.673 seconds (Total)
Chain 2:

SAMPLING FOR MODEL ‘seir’ NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: neg_binomial_2_lpmf: Location parameter[121] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: neg_binomial_2_lpmf: Location parameter[16] is -1.46082e-07, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: neg_binomial_2_lpmf: Location parameter[4] is -2.43315e-08, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: neg_binomial_2_lpmf: Location parameter[49] is -2.27157e-08, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: neg_binomial_2_lpmf: Location parameter[6] is -2.59156e-08, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 3: Rejecting initial value:
Chain 3: Error evaluating the log probability at the initial value.
Chain 3: Exception: neg_binomial_2_lpmf: Location parameter[22] is -2.83628e-07, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 3:
Chain 3: Gradient evaluation took 0.003 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 30 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 22.072 seconds (Warm-up)
Chain 3: 18.496 seconds (Sampling)
Chain 3: 40.568 seconds (Total)
Chain 3:

SAMPLING FOR MODEL ‘seir’ NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[10] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[196] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[10] is -1.51708e-08, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[79] is -3.07714e-12, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[146] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[9] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[15] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[4] is -2.19999e-08, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[109] is -2.32483e-09, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[24] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: neg_binomial_2_lpmf: Location parameter[17] is 0, but must be > 0! (in ‘model5e9058fa2d4b_seir’ at line 57)

Chain 4:
Chain 4: Gradient evaluation took 0.001 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 415.724 seconds (Warm-up)
Chain 4: 561.561 seconds (Sampling)
Chain 4: 977.285 seconds (Total)
Chain 4:
Warning messages:
1: There were 1519 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
2: Examine the pairs() plot to diagnose sampling problems
3: The largest R-hat is 2.42, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess

2

Please, where should the command

stan(..., control = list(max_treedepth = 15))

for warming message 1 be implemented. Thanks

3

Hi, you have a severely misspecified model. My guess is, it won’t matter if you increase adapt_delta or max_treedepth.

1 Like
4

Thank you Torkar for your feedback, I know this might sound daft, please can you make clear to me what you mean by " misspecified model" and how?. I am a maths student trying to familiarize Bayesian and this programme. Thanks

5

Aaah, sorry about that.

You have many transitions that exceeded max treedepth (a question of efficiency) and you have very large Rhat values. Then ESS indicates that very few samples were recovered by the model. In short, this all points to some problems with your model (it could be several things showing this behavior). I recommend you to read this which I found helpful myself:
https://mc-stan.org/misc/warnings.html

6

Alright then. Let me keep digging. Thank you.

1 Like