Fitting the model based on C'(t) data

Modeling
1

I would like to fit SEIR model in Stan, but there is a problem. In fact, I have a problem in model section where my data is C’(t) or dC/dt at each day. How can I fit? The code is attached here.

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

      real S = y[1];
      real E = y[2];
      real I = y[3];
      real R = y[4];
      real C = y[5];
      real N = x_i[1];
      
      real beta = theta[1];
      real gamma = theta[2];
      real kappa = theta[3];     
      real dS_dt = -beta * I * S / N;
      real dE_dt = beta * I * S / N - kappa * E;
      real dI_dt =  kappa * E - gamma * I;
      real dR_dt =  gamma * I;
      real dC_dt = kappa * E;
      
      return {dS_dt, dE_dt, dI_dt, dR_dt, dC_dt};
  }
}

 
data {
  int<lower=1> n_days;
  real y0[5];
  real t0;
  real ts[n_days];
  int N;
}

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

// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
  real<lower=0> gamma;
  real<lower=0> beta;
  real<lower=0> kappa;
  real<lower=0> phi_inv;
}


transformed parameters{
  real y[n_days, 5];
  real phi = 1. / phi_inv;
  {
    real theta[3];
    theta[1] = beta;
    theta[2] = gamma;
    theta[3] = kappa;

    y = integrate_ode_rk45(seir, y0, t0, ts, theta, x_r, x_i);
  }
}

// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
  //priors
  beta ~ normal(2, 1); //truncated at 0
  gamma ~ normal(0.4, 0.5); //truncated at 0
  kappa ~ normal(0.3, 0.5);
  phi_inv ~ exponential(5);
  
  // compute dc/dt
  for (i in 2:n_days) {
    real dc_dt = (y[i, 5] - y[i - 1, 5]) / (ts[i] - ts[i - 1]);
    int dc_dt_int = floor(dc_dt); // cast to integer
    target += neg_binomial_2_lpmf(dc_dt_int | phi, N);
  }
}






generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  real pred_cases[n_days];
  // predict cases based on c(t)
  pred_cases[1] = neg_binomial_2_rng(y[1, 5], phi);
  for (i in 2:n_days) {
    real dc_dt_pred = (y[i, 5] - y[i - 1, 5]) / (ts[i] - ts[i - 1]);
    pred_cases[i] = neg_binomial_2_rng(dc_dt_pred, phi);
  }
}
2

I think the code appended below is great but it does not really work for simulated data.

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

      real S = y[1];
      real E = y[2];
      real I = y[3];
      real R = y[4];
      real C = y[5];
      real N = x_i[1];
      
      real beta = theta[1];
      real gamma = theta[2];
      real kappa = theta[3];
      
      real dS_dt = -beta * I * S / N;
      real dE_dt = beta * I * S / N - kappa * E;
      real dI_dt =  kappa * E - gamma * I;
      real dR_dt =  gamma * I;
      real dC_dt = kappa * E;
      
      return {dS_dt, dE_dt, dI_dt, dR_dt, dC_dt};
  }
}

 
data {
  int<lower=1> n_days;
  real y0[5];
  real t0;
  real ts[n_days];
  int N;
  int cases[n_days];
}

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

// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
  real<lower=0> gamma;
  real<lower=0> beta;
  real<lower=0> kappa;
  real<lower=0> phi_inv;
}



transformed parameters{
  real y[n_days, 5];
  real phi = 1. / phi_inv;
  {
    real theta[3];
    theta[1] = beta;
    theta[2] = gamma;
    theta[3] = kappa;

    y = integrate_ode_rk45(seir, y0, t0, ts, theta, x_r, x_i);
    for (i in 2:n_days)
    {
      y[i, 5] = (y[i, 5] - y[i - 1, 5]);
    }
  }
}

// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
  //priors
  beta ~ normal(2, 1); //truncated at 0
  gamma ~ normal(0.4, 0.5); //truncated at 0
  kappa ~ normal(0.3, 0.5);
  phi_inv ~ exponential(5);
  
  cases ~ neg_binomial_2(col(to_matrix(y), 5), phi);
  }
generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  real pred_cases[n_days];
  // predict cases based on c(t)
  pred_cases = neg_binomial_2_rng(col(to_matrix(y), 5) + 1e-5, phi);
}
3

Hi, Hamed. I don’t know the SEIR model well enough to debug your example, but you might want to start with the models coded in this tutorial:

https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html

1 Like