Fitting the model based on C'(t) data

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);
}
}
``````

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);
}

``````

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