Fitting seir in r using stan

Modeling
,
1

Hi, I am having trouble fitting the SEIR model via Stan. I get the following error: ror in mod$fit_ptr() : Exception: variable does not exist; processing stage=data initialization; variable name=y0; base type=double (in ‘model3fe828ac77ef_SEIR’ at line 30. But the SIR code runs just fine. I just do not understand when we reconstruct y0 in the SEIR in what do we list back in R script. Your help is highly appreciated. I copied the code from Bayesian workflow for disease transmission modeling in Stan

functions {
  real[] sir(real t, real[] y, real[] theta, 
             real[] x_r, int[] x_i) {
             
      real N = x_i[1];
      
      real lambda = theta[1];
      real mu = theta[2];
      real sigma = theta[3];
      real i0 = theta[4]; // initial number of infected people
      real e0 = theta[5]; // initial number of infected people
      
      real init[4] = {N - i0 - e0, e0, i0, 0}; // 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 = - lambda * I * S / N;
      real dE_dt =  lambda * I * S / N - sigma * E;
      real dI_dt =  sigma * 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];
  int N;
  int cases[n_days];
}

transformed data {

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

parameters {
  real<lower=0> mu;
  real<lower=0> lambda;
  real<lower=0> sigma;
  real<lower=0> phi_inv;
  real<lower=0> i0; // number of infected people inititally
  real<lower=0> e0; // number of exposed people inititally
  real<lower=0, upper=1> p_reported;
  }

transformed parameters{
  real incidence[n_days - 1];
  real y[n_days, 4];
  real phi = 1. / phi_inv;
  real theta[5] = {lambda, mu, sigma, i0, e0};
    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))
}


model {
  //priors
  lambda ~ lognormal(log(0.25), 0.8); 
  mu ~ lognormal(log(0.25), 0.5); 
  sigma ~ lognormal(log(0.2), 0.5);
  phi_inv ~ exponential(5);
  i0 ~ 	cauchy(0, 20);
  e0 ~ cauchy(0, 20);
  p_reported ~ beta(1, 2);

  //compute likelihood to the data
  cases[1:(n_days-1)] ~ neg_binomial_2(incidence, phi);
}

My R script is

#----FITTING SEIR MODEL----
#importing the data
Data = read.csv(“C:/Users/Lebo/Desktop/thesis/Data.csv”, header = T)
attach(Data)
Data.csv (152.7 KB)

#initial conditions
N = 2374697;

times

n_days = length(new_cases)
t = seq(0, n_days, by = 1)
t0 = 0
t = t[-1]

number of MCMC steps

niter = 2000

data_seir = list(n_days = n_days, t0 = t0,
ts = t, N = N, cases = new_cases)

model_seir = stan_model(“SEIR.stan”)

fit_seir = sampling(model_seir,
data = data_seir,
iter = niter,
warmup = 1000,
chains = 4,
seed = 0)

traceplot(fit_seir, pars = c(“lambda”, “sigma”, “mu”, “lp__”))

2

@andrjohns Please assist.

3

The error message means that the model is expecting a data variable y0, but you are not submitting this variable.

Note that the Stan model you posted does not have a data variable y0, neither does the data-list you specify in the R code:

However, some of the models in the blog post you link to have a variable y0. So it seems that your file “SEIR.stan” has a different model than the one you quoted above, specifically a model that has a variable y0. So I think a first step is to check that you have the right model in your file “SEIR.stan”.

1 Like
4

Thank you for the response. The model that has y0 defined is the SIR and it runs just fine for me. However, the transformed parameters block for the SEIR model we have y0 replaced with

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

and y0 is not listed in R script according to the paper I am following and without listing I get that error; variable does not exist.

functions {
  real[] sir(real t, real[] y, real[] theta, 
             real[] x_r, int[] x_i) {
             
      real N = x_i[1];
      
      real lambda = theta[1];
      real mu = theta[2];
      real sigma = theta[3];
      real i0 = theta[4]; // initial number of infected people
      real e0 = theta[5]; // initial number of infected people
      
      real init[4] = {N - i0 - e0, e0, i0, 0}; // 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 = - lambda * I * S / N;
      real dE_dt =  lambda * I * S / N - sigma * E;
      real dI_dt =  sigma * 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 y0[4];
  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> mu;
  real<lower=0> lambda;
  real<lower=0> sigma;
  real<lower=0> phi_inv;
  real<lower=0> i0; // number of infected people inititally
  real<lower=0> e0; // number of exposed people inititally
  real<lower=0, upper=1> p_reported;
  }

transformed parameters{
  real incidence[n_days - 1];
  real y[n_days, 4];
  real phi = 1. / phi_inv;
  real theta[5] = {lambda, mu, sigma, i0, e0};
    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))
}

I did include y0 in the data block here and listed it as follows in R script

data_seir = list(n_days = n_days, y0 = y0, t0 = t0, 
                 ts = t, N = N, cases = new_cases)

But i get the error; Error: object ‘y0’ not found
so i am wondering how then is the rep_array(0.0, 4) called in R script

5

sorry for the late reply

This is an error from R, which suggests that you did not define a variable y0 (e.g. y0 = 1) before you made the list data_seir.

Note that I didn’t look at the vignette again, so there might be other issues.

6

Thank you so much and i am so sorry for replying so late. Could you kindly assist on how i define the variable y0