Dear All,
I want to estimate parameters (beta and alpha) of an SIRD Model whose equations are given by
dS_dt = -beta * I * S / N,
dI_dt = beta * I * S / N - gammaI - alpha * I,
dR_dt = gamma I,
dD_dt = alpha * I.
I already know the value of gamma from literature, it is gamma=0.07. How do I capture this value in my Rstan code? I have attached my code and the dataset that I am using. I am following Bayesian workflow for disease transmission modeling in Stan. Thank you in advance.
Use transformed data, so that’d just be:
transformed data {
real gamma = 0.07;
}
Then you can use it just like any variable, such as in the definition of the diff eq system. Or you can just roll it into your code. Where you have
gamma = 0.07
you can just use
real gamma = 0.07;
and everything should work.
P.S. Pasting your stan code into the chat is more helpful than links which need to be downloaded. I’d also recommend keeping it in its own file so that you can evaluate line numbers in error messages and you don’t have to escape either " or ' in the code. Anyway, here it is for the next person.
functions {
real[] sir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i){
real S = y[1];
real I = y[2];
real R = y[3];
real D = y[4];
real N = x_i[1];
real beta = theta[1];
real alpha = theta[2];
gamma = 0.07; # How best can I incorperate this value of gamma????
real dS_dt = -beta * I * S / N;
real dI_dt = beta * I * S / N - gamma*I - alpha * I;
real dR_dt = gamma* I;
real dD_dt = alpha * I;
return {dS_dt, dI_dt, dR_dt, dD_dt};
}
}
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];
int x_i[1] = { N };
}
parameters {
real<lower=0, upper = 2> beta;
real<lower=0, upper = 1> alpha;
real<lower=0, upper = 10000> phi_inv;
}
transformed parameters{
real y[n_days, 4];
real phi = 1. / phi_inv;
{
real theta[3];
theta[1] = beta;
theta[2] = alpha;
y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
}
}
model {
//priors
beta ~ uniform(0, 2);
alpha ~ uniform(0, 1);
phi_inv ~ exponential(5);
//sampling distribution
//col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
}
generated quantities {
real R0 = beta / (alpha+gamma);
real recovery_time = 1 / gamma;
real pred_cases[n_days];
pred_cases = neg_binomial_2_rng(col(to_matrix(y), 2), phi);
}
@Bob_Carpenter Thank you very much for your response. It did work but R_0 is not picking the value of gamma. There is an error coming from the Stan code which says, gamma is not defined.
Is there anything else I need to do? I have pasted below my Stan and R code.
functions {
real[] sir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i){
real S = y[1];
real I = y[2];
real R = y[3];
real D = y[4];
real N = x_i[1];
real beta = theta[1];
real alpha = theta[2];
real gamma = 0.07;
real dS_dt = -beta * I * S / N;
real dI_dt = beta * I * S / N - gamma*I - alpha * I;
real dR_dt = gamma* I;
real dD_dt = alpha * I;
return {dS_dt, dI_dt, dR_dt, dD_dt};
}
}
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];
int x_i[1] = { N };
}
parameters {
real<lower=0, upper = 2> beta;
real<lower=0, upper = 1> alpha;
real<lower=0, upper = 10000> phi_inv;
}
transformed parameters{
real y[n_days, 4];
real phi = 1. / phi_inv;
{
real theta[2];
theta[1] = beta;
theta[2] = alpha;
y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
}
}
model {
//priors
beta ~ uniform(0, 2);
alpha ~ uniform(0, 1);
phi_inv ~ exponential(5);
//sampling distribution
//col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
}
generated quantities {
real R0 = beta / (alpha + 0.07);
real recovery_time = 1 / 0.07;
real pred_cases[n_days];
pred_cases = neg_binomial_2_rng(col(to_matrix(y), 2), phi);
}
library(ggplot2)
library(rstan)
library(gridExtra)
#library(tidyverse)
setwd("~/Desktop/Rstan")
datalist <- read.csv(file = 'my_dataset.csv')
head(datalist)
ggplot(data = datalist) +
geom_point(mapping = aes(x = day, y = cases0)) +
labs(y = "Number of daily infected people")
rstan_options (auto_write = TRUE)
options (mc.cores = parallel::detectCores ())
# total count
N <- 99999;
# times
cases <- datalist$cases0
n_days <- length(cases)
t <- seq(0, n_days, by = 1)
t0 = 0
t <- t[-1]
#initial conditions
i0 <- 1
s0 <- N - i0
r0 <- 0
d0 <- 0
y0 = c(S = s0, I = i0, R = r0, D=d0)
# data for Stan
data_sir <- list(n_days = n_days, y0 = y0, t0 = t0, ts = t, N = N, cases = cases)
# number of MCMC steps
niter <- 2000
model <- stan_model("sird.stan")
fit_sir_negbin <- sampling(model,data = data_sir,iter = niter,chains = 4,seed = 0)
pars=c('beta', 'alpha', "R0", "recovery_time")
print(fit_sir_negbin, pars = pars)
stan_dens(fit_sir_negbin, pars = pars, separate_chains = TRUE)
smr_pred <- cbind(as.data.frame(summary(
fit_sir_negbin, pars = "pred_cases", probs = c(0.05, 0.5, 0.95))$summary), t, cases)
colnames(smr_pred) <- make.names(colnames(smr_pred)) # to remove % in the col names
# We fit the model to data below:
ggplot(smr_pred, mapping = aes(x = t)) +
geom_ribbon(aes(ymin = X5., ymax = X95.), fill = "red", alpha = 0.6) +
geom_line(mapping = aes(x = t, y = X50.)) +
geom_point(mapping = aes(y = cases)) +
labs(x = "Day", y = "Number of students in bed")
params <- lapply(t, function(i){sprintf("y[%s,2]", i)}) #number of infected for each day
smr_y <- as.data.frame(summary(fit_sir_negbin,
pars = params, probs = c(0.05, 0.5, 0.95))$summary)
colnames(smr_y) <- make.names(colnames(smr_y)) # to remove % in the col names
ggplot(smr_y, mapping = aes(x = t)) +
geom_ribbon(aes(ymin = X5., ymax = X95.), fill = "orange", alpha = 0.6) +
geom_line(mapping = aes(x = t, y = X50.), color = "orange") +
labs(x = "Day", y = "Number of infected students")
traceplot(fit_sir_negbin, pars = c('beta', 'alpha', "R0", "recovery_time"))
check_hmc_diagnostics(fit_sir_negbin)
pairs(fit_sir_negbin, pars = c("beta", "alpha", "R0", "recovery_time"))