Bayesian Modelling of SEIR

Modeling
1

Hello, I have the following stan code for SEIR modelling. I need help with setting the domain for the model priors in the parameters block since my priors assume negative values i.e lambda ~ lognormal(log(0.4), 0.5); and i have in the parameters block, real<lower=0> mu; which i have to change to include negative values. The reference code is from Bayesian workflow for disease transmission modeling in Stan

functions {
  real[] seir(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; //population size
  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> reporting_D;
  }

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(seir, 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]) * reporting_D  ; //-(E(t+1) - E(t) + S(t+1) - S(t))
}

model {
  //priors
  lambda ~ lognormal(log(0.4), 0.5); 
  mu ~ lognormal(log(0.125), 0.2); 
  sigma ~ lognormal(log(0.2), 0.5);
  phi_inv ~ exponential(5);
  i0 ~ 	cauchy(0, 20);
  e0 ~ cauchy(0, 20);
  reporting_D ~ lognormal(log(8), 0.2);
  
  //compute likelihood to the data
  cases[1:(n_days-1)] ~ neg_binomial_2(incidence, phi);
}


generated quantities {
  real R0 = lambda / mu;
  real recovery_time = 1 / mu;
  real incubation_time = 1 / sigma;
  real pred_cases[n_days - 1];
  pred_cases = neg_binomial_2_rng(incidence, phi);
}
1 Like
2

In your model mu (gamma in the case study you mention, and more commonly in the literature) is the absolute value of the recovery rate , and it carries the negative sign in the structure of the model, so I believe it really should be constrained to being positive only.

Most things that are rates should be strictly positive for biological reasons, i.e. there are no negative infections, returning from the infected to the exposed compartment, or “unrecovering”, etc.

If you are having issues with the ranges of parameters the causes are probably other than the ranges themselves.

2 Likes
3

Thank you so much for your help.

1 Like