Initialization error in multi-season SIR model

Modeling
1

I am trying to fit an SIR model with prior distributions on the model parameters as well as on the initial value for the proportion of the population already infected (the I compartment).

functions {
  
  real[] SIR(real t,
             real[] y,
             real[] theta,
             data real[] x_r,
             data int [] x_i) {
             
             real S = y[1];
             real I = y[2];
             real R = y[3];
             real N = x_i[1];
             
             real beta = theta[1];
             real gamma = theta[2];
             
             real dSdt = -beta * S * I / N;
             real dIdt = beta * S * I / N - gamma * I;
             real dRdt = gamma * I;
             
             return {dSdt, dIdt, dRdt};
             
         }
}

data {
int<lower=0> N;                                 // number of points
int<lower=0> n_seasons;                         // number of seasons
int<lower=0> seasons[N];                        // seasons
real<lower=0> ili[N];                           // all ili data
int<lower=0> x[N];                              // week
real<lower=0> n_weeks[43];       


int<lower=0> ps;                                // season to predict
real<lower=0,upper=1> S0;                       // initial susceptible proportion
real t0;                                        // week 0
}

transformed data {
  real x_r[0];
  int x_i[1] = { 1 };
}

parameters {


vector[n_seasons] beta_s;
vector[n_seasons] rho_s;
real<lower=0,upper=.1> I0_s[n_seasons];
real<lower=0> epsilon;


}



transformed parameters {


vector[n_seasons] gamma_s;

for (i in 1:n_seasons) {
                gamma_s[i] = beta_s[i] * rho_s[i];
                R0 = 1 - S0 - I0_s[i];
}



}

model {

for (i in 1:n_seasons) {
  beta_s[i] ~ gamma(2,.02);
  rho_s[i] ~ gamma(2,.02);
  I0_s[i] ~ normal(.05,.02);
}

epsilon ~ exponential(100);



for (i in 1:N) {
                 
                 ili[i] ~ normal(col(to_matrix(integrate_ode_rk45(SIR, 
                 {S0, I0_s[seasons[i]], 1 - S0 - I0_s[seasons[i]]}, 
                              t0, n_weeks, 
                              {beta_s[seasons[i]], gamma_s[seasons[i]]},
                              x_r, x_i)),2),
                              epsilon);
}
                          
                          
}

I run the model in R using rstan and I get the following error

Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: loop iteration: 1

Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: Max number of iterations exceeded (1000000).  (in 'model27678c0a492_af4cfd756eb609d1cf3b6b499e742e03' at line 96)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                              
[2] "  Exception: Max number of iterations exceeded (1000000).  (in 'model27678c0a492_af4cfd756eb609d1cf3b6b499e742e03' at line 96)"
error occurred during calling the sampler; sampling not done

I’m very lost trying to figure this one out. Any ideas what I might be doing wrong? Also included is the code for the data list I’m using in R, as well as a .csv
file of my data.
ILINet_samp.csv (19.1 KB)

dat = list(N = nrow(ILINet_samp),
           n_seasons = length(unique(ILINet_samp$season)),
           seasons = as.numeric(factor(ILINet_samp$season)),
           ili = ILINet_samp$unweighted_ili/100,
           x = ILINet_samp$season_week,
          # n_weeks = ILINet_samp$n_weeks,
           n_weeks = 1:43,
           ps = max(as.numeric(factor(unique(ILINet_samp$season)))),
           S0 = .9,
           t0 = 0)
2

I’m guessing you want to fit a different ODE solution to each of what you disease seasons (e.g. each year), burt you probably shouldn’t be calling the integrator for each time point, and it seems that’s what you are doing in the for (i in 1:N) loop over there.

You may want to call the integrator for each season in this loop (if you have the same prior for all seasons you can skip loop and just list the priors as beta_s ~ gamma(2,.02) for instance).

It seems like the i indices are getting mixed up between season number and data point indices. It may also be helpful to compute the ODE solution in higher up in the model block (or in the transformed parameters block if you want it saved to the output) and write it more concisely in the final lines, it may make it easier to spot any errors.

That said, even after correcting any problems with the code, that error is pretty common with ODEs, so you may also want to choose initial values for the parameters with the init keyword (maybe it’s inits in rstan, I’m not sure).