Hi everyone,
I’m quite new to modeling and Stan, so this might be a stupid question, but I’m still on a learning curve!
I’m trying to model transmission of an infectious disease with yearly rates defining the flow between states (compartments). However, I’d like the ODE integrator (rk45) to solve the equations at a monthly time step to have smoother curves.
A particularity of my model is also that I’m fitting to 1 prevalence sample (1 data point).
I can’t find any examples of this, and when I try doing it on my own, I run into quite a few errors.
Here are some relevant parts of my Stan code (i’ve removed the paramters and equations etc to make it easier to read), with the integrator solving yearly for 40 yrs if I’m understanding and coding correctly:
functions {
real[] dtb1(real t, real[] y, real[] theta,
real[] x_r, int[] x_i) {
//define interger data//
int n_cases = x_i[1];
int n_samples = x_i[2];
int N = x_i[3];
//define real data//
//define paramaters needed to solve model//
//define compartments & init conditions //
//define states & ODE//
return {dS_dt, dE_dt, dL_dt, dI_dt, dR_dt};
}
}
data {
int<lower = 1> n_ts; // number of times to solve for
int<lower = 1> N; // population size
int<lower = 1> n_cases; // observed cases
int<lower = 1> n_samples; // number of data points to fit to
int<lower = 1> n_comp; // number of compartments
real t0; //initial time point 0
real ts[n_ts];
int y[n_samples, 1]; // data, reported prev in ea sample
}
transformed data {
real x_r[2] = {f, r};
int x_i[3] = {n_cases, n_samples, N};
}
parameters {
}
transformed parameters {
real y_hat[n_ts, n_comp]; // Matrix to store the integrated ODE solutions
real theta[9] = {beta, mu, gamma, epsilon, delta, omega, chi, pi, rho};
real init[n_comp] = {1231, 524, 1057, 150, 2036};
simplex [n_comp] output[n_ts];
// Solve the ODE using integrate_ode_rk45
y_hat = integrate_ode_rk45(dtb1, init, t0, ts, theta, x_r, x_i);
for(i in 1:n_ts) output[i] = ( to_vector(y_hat[i,1:5]) ) ./ sum( to_vector(y_hat[i,1:5]));
// get the proportion in each category
}
model {
// sampling distribution
for(i in 1:n_ts) n_cases ~ binomial(N, output[i,4]);
}
generated quantities {
array[n_ts] real pred_cases;
for (i in 1:n_ts)
pred_cases[i] = binomial_rng(N, output[i,4]);
}
',
file = "dtb1_binom.stan", sep="", fill=T)
#data
n_cases<-c(47)
N<-1080
t<-1
ts<-seq(1, 40, by=1)
n_ts<-40
y <- array(c(I = n_cases), dim = c(1, 1))
stan_d = list(n_samples =1,
N=N,
n_comp = 5,
n_ts= n_ts,
n_cases=n_cases,
ts=ts,
y=y,
f=f,
r=r,
t0 = 0)
#compile stan model
model <- stan_model("dtb1_binom.stan")
fit <- sampling(model,
data = stan_d,
iter = 1000,
chains = 4,
seed = 0)
Thank you in advance for any help!
Morgana