# Estimation of time-varying infection parameter for a simple SEIR model

Dear all,

I have been stuck on this topic for some time now. I am trying to estimate time-varying infection (with other time-independent epidemiological) parameters for a complex multi-patch SEIRS model. But I decided to first implement the idea for a simple single-patch SEIRS model. I am not quite sure how to achieve this, but I imagine that I need to use a for loop at the transformed parameters block, to evaluate the ODE integrator at each time point. I tried implementing this, but on compiling, I got the following error:

86:
-------------------------------------------------
123:  	//params[4] = nu;
124:
125:  	y[1,] = ode_rk45(SEIRS, y0, -0.000001, dt[1],beta, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff
^
126:
127:
-------------------------------------------------

Ill-typed arguments supplied to function 'ode_rk45':
(<F1>, vector, real, real, vector, real, real, real, real, real, real)
where F1 = (real, vector, real, real, real, real, real, real, real) => vector
Available signatures:
(<F1>, vector, real, array[] real, real, real, real, real, real, real, real) => vector
The fourth argument must be array[] real but got real

I am using the new ODE interface, which requires that the times be given as array real. How to provide a single real time at each iteration is not very clear for me. I will appreciate it if you could help me debug the following Stan code, which I use cmstanpy to run.

functions {
array[] real diff(array[] real x){
int N = num_elements(x);
array[N-1] real result;
for (i in 1:(N-1)){
result[i] = x[i+1] - x[i];
}
return result;
}
vector SEIRS( real t,        // time
vector y,      // State S, E, I, R, N
real beta, // infection parameter beta
real kappa,
real gama,
real Lambda,
real mu,
real tau,
real phi
//vector n_obs
//vector n_fake
//real[] x_r,    // real ODE argument that depend only on data
//int[] x_i
){  // integer ODE argument that depend only on data

vector[6] dydt;

dydt[1] = Lambda*y[6] - beta*y[1]*y[3]/y[6] - mu*y[1] + tau*y[4];
dydt[2] = beta*y[1]*y[3]/y[6] - (kappa + mu) * y[2];
dydt[3] = kappa * y[2] - (gama + phi + mu) * y[3];
dydt[4] = gama * y[3] - (tau + mu) * y[4];
dydt[5] = kappa * y[2];
dydt[6] = Lambda*y[6] - (mu * y[6] + phi * y[3]);

return dydt;

}
}

data{
int<lower = 1> n_obs;              // Number of days sampled
int<lower = 1> n_fake;             // number of fake days sampled for prediction
array[n_obs] int<lower = 0> cases; // The observed data  (daily incidence)
vector<lower = 0>[6] y0;       //initial conditions
real<lower = 0> t0;                // Initial time point (zero)
array[n_obs+1] real<lower = t0> dt;    // Time points that were sampled
array[n_fake] real<lower = t0> fake_ts; // Fake prediction time points that were sampled
array[n_obs]  real<lower = 0> mob_data; //mobility data
real<lower = 0> N0;          // Initial pop size
real<lower = 0> R0;          //
real<lower = 0> Lambda;      // We considere known
real<lower = 0> mu;          // We considere known
real<lower = 0> tau;         // We considere known
real<lower = 0> phi;         // We considere known
}

//transformed data {
//  real x_r[4] = {Lambda, mu, tau, phi};
//  int x_i[0];
//}

parameters{
//Support of parameters
//real<lower =  0.5, upper = 2.0> beta_t;
vector<lower = 0.8, upper = 1.60>[n_obs] beta;
//vector<lower = 0.5, upper=1.6>[n_obs] beta0;
//real<lower = 0.5, upper=1.6> beta0;
//real<lower = 0.5, upper=1.6> beta1;
real<lower =  0.5, upper = 1.5> kappa;
real<lower =  0.5, upper = 1.5> gama;
//real<lower =  0.1> nu_inv;
real<lower =  0.1> nu_inv;
//real<lower =  0, upper = 15> E0;
//real<lower =  0, upper = 30> I0;
//real<lower = 5e-7> nu;
}

transformed parameters{
array[n_obs+1] vector<lower = 0.>[6] y;    // Output from the ODE solver
array[n_fake] vector<lower = 0.>[6] y_pred_Inc; // Output from the ODE solver for prediction
real nu = 1./nu_inv;

//vector[n_obs+1] beta =
//{
//  vector[6] y0;                 // initial conditions
//  y0[1]   = N0 - ( E0 + I0 + R0);
//  y0[2]   = E0;
// y0[3]   = I0;
//  y0[4]   = R0;
//y0[5]   = E0 + I0 + R0;
//y0[6]   = N0;

{
//real params[3];
//params[1] = beta;
//params[2] = kappa;
//params[3] = gama;
//params[4] = nu;

y[1,] = ode_rk45(SEIRS, y0, -0.000001, dt[1],beta, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff

y_pred_Inc[1,] = ode_rk45(SEIRS, y0, -0.000001, fake_ts[1], beta, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff (for predictions)

for (t in 2:n_fake){
//real beta_t = exp(beta0[t] + beta1[t]*mob_data[t]);
real beta_t = beta[t]; //+ beta1*mob_data[t];
//real dt[1] = t - 1;
//real fake_ts[1] = t - 1;
//}

y[t,] = ode_rk45(SEIRS, y[t-1,], -0.000001, dt[t],beta_t, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff

// y0 = y[t];

y_pred_Inc[t,] = ode_rk45(SEIRS, y_pred_Inc[t-1,], -0.000001, fake_ts[t], beta_t, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff (for predictions)

// y0 = y_pred_Inc[t];

}

// }
}
}

model{
real epsilon = 1e-4;
// Prior distributions
for (t in 1:n_obs){
beta[t]    ~ normal(-0.1115975, 0.2966405);
}
//beta0    ~ lognormal(-0.1115975, 0.2966405);
//beta1    ~ lognormal(-0.1115975, 0.2966405);
kappa   ~ gamma(4.25173, 6.09586);
gama    ~ gamma(4.25173, 6.09586);
//nu_inv ~ gamma(9.0045, 2.2508);
nu ~ gamma(0.64906, 0.28971);
//E0      ~ uniform(3.1699, 11.8301);
//I0      ~ uniform(6.3397, 23.6603);

// Likelihood

cases ~ neg_binomial_2(diff(y[ ,5]), nu+epsilon);
}
generated quantities {
real epsilon = 1e-4;
array[n_fake-1] real pred_Inc_cases;

pred_Inc_cases = neg_binomial_2_rng(diff(y_pred_Inc[ ,5]), nu+epsilon);
}