Hi,

This is my first try with Stan and R. I’m trying to calibrate the contact rate in an SEIR model from incidence data. I assume the incidence data is binomially distributed, where p indicates the probability of a case being reported.

My problem is that the incidence variable is real vector, and the binomial function requires an integer vector.

Is there a way to make this work or should I change my approach?

```
functions {
real[] SEIR(real t,
real[] y, // stocks
real[] params,
real[] x_r,
int[] x_i) {
real dydt[4];
real forceOfInfection;
real population;
real prevaleceOfInfection;
real newInfections;
real completingLatentPeriod;
real recovery;
population = y[1] + y[2] + y[3] + y[4];
prevaleceOfInfection = y[3] / population;
forceOfInfection = params[1] * prevaleceOfInfection * 0.005; //beta
newInfections = y[1] * forceOfInfection;
completingLatentPeriod = y[2] / 0.42857; //Tao
recovery = y[3] / 1; //Mu
dydt[1] = - newInfections;
dydt[2] = newInfections - completingLatentPeriod;
dydt[3] = completingLatentPeriod - recovery;
dydt[4] = recovery;
return dydt;
}
}
data {
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_params; // Number of model parameters
int<lower = 1> n_difeq; // Number of differential equations in the system
int y[n_obs]; // The binomially distributed data
real t0; // Initial time point (zero)
real ts[n_obs + 1]; // Time points that were sampled
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 0> params[n_params]; // Model parameters
real<lower = 0, upper = 1> p;
}
transformed parameters{
real y_hat[n_obs + 1, n_difeq]; // Output from the ODE solver
real y0[n_difeq]; // Initial conditions
real incidence[n_obs];
y0[1] = 999;
y0[2] = 1;
y0[3] = 0;
y0[4] = 0;
y_hat = integrate_ode_rk45(SEIR, y0, t0, ts, params, x_r, x_i);
for (i in 0:n_obs) {
incidence[i] = y_hat[i + 1, 3] - y_hat[i, 4] + y_hat[i + 1, 4];
}
}
model {
p ~ uniform(0, 1);
params[1] ~ uniform(0, 500);
y ~ binomial(incidence, p);
}
```

R code:

```
library(rstan)
stan_d = list(n_obs = 32,
n_params = 1,
n_difeq = 4,
y = c(1,1,2,2,4,7,7,12,22,28,41,49,59,57,60,51,45,27,19,19,8,5,6,2,3,1,1,0,0,0,0,0),
t0 = 0,
ts = 0:31)
# Test / debug the model:
test = stan("SEIR.stan",
data = stan_d,
chains = 1, iter = 10)
```