Dear all,
I’m new to Stan and I’m currently following the tutorial on “Bayesian workflow for disease transmission modeling in Stan” to develop an SEIR model for the transmission of hand, foot, and mouth disease (HFMD).
To begin, I generated simulated data using the deSolve package in R and then attempted to estimate the model parameters using the simulated data in Stan. However, I have encountered three issues during this process:
- The transmission rate of HFMD exhibits a seasonal pattern, and I have attempted to incorporate this using a cosine function in both R and Stan. I am uncertain whether I have implemented this correctly in Stan and would appreciate guidance on the proper approach.
- Despite obtaining a good fit for the “Infections” variable, I consistently observe negative values for the “Susceptibles” variable. The parameters I utilized to generate the simulated data in R did not result in negative “Susceptibles”. I am unsure of the cause behind this issue and would appreciate assistance in identifying a potential solution.
- Initially, I employed the negative binomial distribution for the case incidence, but I consistently encountered initialization failures. Consequently, I switched to using the Poisson distribution. I am uncertain about the principles governing the selection of priors for the parameter “phi” or “phi_inv” in a negative binomial distribution. Additionally, I am unsure whether it is appropriate to utilize the Poisson distribution in my particular situation.
Could anyone please help me with the above issues? Thanks so much.
The following is the R code that I used to generate the simulated data and the Stan code for estimation.
library(deSolve)
library(tidyr)
library(dplyr)
library(rstan)
library(gridExtra)
rstan_options (auto_write = TRUE)
options (mc.cores = parallel::detectCores ())
seasonal.sir.model <- function (t, x, params) {
with(
as.list(c(x, params)),
{
beta <- beta0*(1+beta1*cos(2*pi*t/180))
dS <- mu*(N-S)-beta*S*I/N
dI <- beta*S*I/N -(mu+gamma)*I
dR <- gamma*I-mu*R
dI_hat <- beta*S*I/N
res <- c(dS,dI,dR,dI_hat)
list(res)
}
) }
times <- seq(0,365*4,by=1)
params <- c(mu=1/700,N=8000000,beta0=0.25,
beta1=0.5, gamma=1/14)
xstart <- c(S=8000000-10,I=10,R=0,I_hat=10)
out <- as.data.frame(ode(xstart,times,seasonal.sir.model,params,rtol=1e-12,hmax=1/120))
cases=as.integer(lead(out$I_hat,n=1L)-out$I_hat)
cases=cases[which(is.na(cases)==F)]
# times
n_days <- length(cases)
t <- seq(0, n_days, by = 1)
t0 = 0
t <- t[-1]
#initial conditions
N=8000000
i0 <- 10
s0 <- N-i0
r0 <- 0
y0 = c(S = s0, I = i0, R = r0,I_hat=i0)
# data for Stan
data_seir <- list(n_days = n_days, y0 = y0, t0 = t0, ts = t, N = N, cases = cases)
stan model:
real[] seir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i) {
real N = x_i[1];
real beta0 = theta[1];
real beta1 = theta[2];
real S = y[1];
real I = y[2];
real R = y[3];
real dS_dt = - beta0 *(1+ beta1 *cos(2*3.14*t/180)) * I * S / N + 0.0014 * (N-S);
real dI_dt = beta0 *(1+ beta1 *cos(2*3.14*t/180)) * I * S / N - 0.07 * I - 0.0014*I;
real dR_dt = 0.07 * I - 0.0014 * R;
real dI_hat= beta0 *(1+ beta1 *cos(2*3.14*t/180)) * I * S / N;
return {dS_dt, dI_dt, dR_dt, dI_hat};
}
}
data {
int<lower=1> n_days;
real y0[4];
real t0;
real ts[n_days];
int N;
int cases[n_days];
}
transformed data {
real x_r[0]; // x_r being of length 0 because we have nothing to put in it
int x_i[1] = { N };
}
parameters {
real <lower=0> beta0;
real <lower=0> beta1;
}
transformed parameters{
real y[n_days, 4];
real theta[2] = {beta0,beta1};
real incidence [n_days-1];
y = integrate_ode_rk45(seir, y0, t0, ts, theta, x_r, x_i);
for (i in 1:(n_days-1)) {
incidence[i]= y[i+1, 4] - y[i, 4] ;
}
}
model {
//priors
beta0 ~ normal(0.25, 0.1);
beta1 ~ normal(0.5, 0.1);
cases[1:(n_days-1)] ~ poisson(incidence);
}
generated quantities {
real pred_cases1[n_days-1];
real ysir[n_days, 4]=y;
pred_cases1 = poisson_rng(incidence);
}