Dear Stan-community,
I am modelling a simple SEIR model in stan.
\frac{dS}{dt} = - \beta S I \\ \frac{dE}{dt} = \beta S I - \tau E \\ \frac{dI}{dt} = \tau E - \gamma I \\ \frac{dR}{dt} = \gamma I
Due to limitations in the accuracy of the solver ( I am using ode_rk45_tol, but had similar issues with other solvers ) I sometimes get slightly negative values for the number of individuals in the I comparmtent. This leads to errors in my program, because I use these values as the rate in a poisson process. I though this could be solved by log transforming the model. Then, the number of infected individuals would be the exponential of the output of the model and always be positive. I was hoping this would make my model more efficient.
Log transformation of the SEIR model:
Step 1, divide by variables:
\frac{1}{S} \frac{dS}{dt} = - \beta S I \frac{1}{S} \\
\frac{1}{E} \frac{dE}{dt} = (\beta S I - \tau E) \frac{1}{E} \\
\frac{1}{I} \frac{dI}{dt} = (\tau E - \gamma I) \frac{1}{I} \\
\frac{1}{R} \frac{dR}{dt} = \gamma I \frac{1}{R} \\
Step 2, simplify:
\frac{d logS}{dt} = - \beta I \\
\frac{d logE}{dt} = \frac{\beta S I}{E} - \tau \\
\frac{d logI}{dt} = \frac{\tau E}{I} - \gamma \\
\frac{d logR}{dt} = \frac{\gamma I}{R}
When I analyse this model and the original SEIR model in r, I find that both are equivalent. I also included a model in the code that is scaled by the population size and models the whole population istead of the fraction of infected individuals. The initial conditions cannot be zero, because I cannot logtransform zero. Thus, I choose very small values, this does not really influence the output of the models.
# Load library
library(deSolve)
# Set parameters and time sequence
p <- c(beta = 0.08, gamma = 2/7, tau = 5/7, contact = 77, popsize = 10000)
times <- seq(0, 45, by = 0.01)
# Set initial coditions
s_norm <- c(S = (1-1/10000), E=0, I = 1/10000, R = 0)
s_log <- c(logS = log(1-1/10000), logE=-50, logI = log(1/10000), logR = -50)
s_pop <- c(logS = log(10000-1), logE=log(0.00001), logI = log(1), logR = log(0.00001))
# if initial values for E and R are too high, the result is not the same
s_pop_stan <- c(logS = 9.21024,logE=-1,logI=0,logR=-1)
# Traditional way to code SEIR model
gfull <- function(time, y, params) {
g <- with(as.list(c(y,params)),
c( -beta*contact*S*I,
beta*contact*S*I - tau*E,
tau*E - gamma*I,
gamma *I )
)
return(list(g))
}
# Log transformed SEIR model
gfull_log <- function(time, y, params) {
g <- with(as.list(c(y,params)),
c( - beta*contact*(exp(logI)),
beta*contact*(exp(logI)*(exp(logS)))/(exp(logE)) - tau,
tau*(exp(logE))/(exp(logI)) - gamma,
gamma*(exp(logI))/(exp(logR)) )
)
return(list(g))
}
# Log transformed and scaled with population size
gpop <- function(time, y, params) {
g <- with(as.list(c(y,params)),
c( - beta*contact*(exp(logI))/popsize,
(beta*contact*(exp(logI)*(exp(logS)))/(exp(logE)))/popsize - tau,
tau*(exp(logE))/(exp(logI)) - gamma,
gamma*(exp(logI))/(exp(logR)) )
)
return(list(g))
}
# Run all versions of the model
result_log <- ode(y = s_log, times = times, func = gfull_log, parms = p)
result_norm <- ode(y = s_norm, times = times, func = gfull, parms = p)
result_pop <- ode(y = s_pop, times = times, func = gpop, parms = p)
# compare results for I compartment
par(mfrow=c(1,1))
plot(times, exp(result_log[,3]), col="blue", xlab="time in weeks", ylab="Number of infected individuals", type="l", main="R deSolve model")
lines(times, result_norm[,3], col="red",type="l")
lines(times, exp(result_pop[,3])/10000,col="violet")
legend("topright", legend=c("Traditional model", "Logtransformed model", "Logtransformed population scaled"), col=c("red", "blue","violet"),lty=1)
However, when I run the exact same models in Stan, I get nan values for the ODE output if the initial values are smaller than -11 for the logtransformed model and smaller than log(0.01) for the version that is scaled by the population size.
functions {
// ODE for basic SEIR model
vector gfull(real t, vector y, data real tau, data real gamma, real beta,
data real contact, data int num_comp) {
// define compartments
vector[num_comp] dydt;
// define ODEs
dydt[1] = - beta * contact * y[1] * y[3];
dydt[2] = beta * contact * y[1] * y[3] - tau * y[2];
dydt[3] = tau * y[2] - gamma * y[3];
dydt[4] = gamma * y[3];
return dydt;
}
// ODE for logtransformed SEIR
vector gfull_logtrans(real t, vector y, real I0, data real tau, data real gamma, real beta,
data real contact, data int num_comp) {
// define compartments
vector[num_comp] dydt;
// define ODEs
dydt[1] = - beta*contact*(exp(y[3]));
dydt[2] = (beta*contact*exp(y[3])*exp(y[1]))/exp(y[2]) - tau;
dydt[3] = tau*exp(y[2])/exp(y[3]) - gamma;
dydt[4] = gamma*exp(y[3])/exp(y[4]);
return dydt;
}
// ODE for basic SEIR model logtransformed and scaled with populationsize
vector gfull_pop(real t, vector y, real I0, data real tau, data real gamma, real beta,
data real contact, data int num_comp, data int popsize) {
// define compartments
vector[num_comp] dydt;
dydt[1] = -beta*contact*(exp(y[3]))/popsize;
dydt[2] = (beta*contact*exp(y[3])*exp(y[1]))/(exp(y[2])*popsize) - tau;
dydt[3] = tau*(exp(y[2]))/(exp(y[3])) - gamma;
dydt[4] = gamma*(exp(y[3]))/(exp(y[4]));
return dydt;
}
}
// load data objects
data {
// Load parameters from data
// data objects that are generic independent of the model specifications
int num_t;
int popsize;
int contact;
array[num_t] real ts;
// true values for simulated data
real beta_sim;
real I0_raw_sim;
// fixed quantities
real incubation_period;
real infectious_period;
// control parameters
real rtol;
real atol;
int max_num_steps;
}
transformed data {
// generic transformed data block
int num_comp = 4; // total compartments
int num_prev = 3; // index of the I compartment
real t0 = 0.;
// fixed parameters
real tau = 1./incubation_period;
real gamma = 1./infectious_period;
// scale initial seed
real I0_sim = I0_raw_sim/popsize;
print("initial values traditional model:");
print(to_vector(append_array( {(1-I0_sim), 0, I0_sim}, rep_array(0, num_comp-3))));
// get reference simulation
array[num_t] vector[num_comp] y_sim_traditional = ode_rk45_tol(
gfull,
to_vector(append_array( {(1-I0_sim), 0, I0_sim}, rep_array(0, num_comp-3))), // initial values work well
t0, // initial time = 0
ts, // evaluation times
rtol, atol, max_num_steps, // tolerances
tau, gamma, beta_sim, contact, // data
num_comp // metadata
);
print("initial values log transformed:");
print(to_vector(append_array( {log(1-I0_sim), -11, log(I0_sim)}, rep_array(-11, num_comp-3))));
array[num_t] vector[num_comp] y_sim_log = ode_rk45_tol(
gfull_logtrans,
to_vector(append_array( {log(1-I0_sim), -11, log(I0_sim)}, rep_array(-11, num_comp-3))), // values more negative then -10 lead to Nan result for y
t0, // initial time = 0
ts, // evaluation times
rtol, atol, max_num_steps, // tolerances
I0_sim, tau, gamma, beta_sim, contact, // data
num_comp // metadata
);
print("initial values log transformed population scaled:");
print(to_vector(append_array( {log(popsize-I0_sim*popsize), log(0.1), log(I0_sim*popsize)}, rep_array(log(0.1), num_comp-3))));
array[num_t] vector[num_comp] y_sim_poplog = ode_rk45_tol(
gfull_pop,
to_vector(append_array( {log(popsize-I0_sim*popsize), log(0.1), log(I0_sim*popsize)}, rep_array(log(0.1), num_comp-3))), // values smaller than log(1) lead to Nan result for y
t0, // initial time = 0
ts, // evaluation times
rtol, atol, max_num_steps, // tolerances
I0_sim, tau, gamma, beta_sim, contact, // data
num_comp, popsize // metadata
);
}
parameters {
}
transformed parameters {
}
generated quantities {
// posterior predictive check
array[num_t] vector[num_comp] y_traditional = y_sim_traditional;
array[num_t] vector[num_comp] y_log = y_sim_log;
array[num_t] vector[num_comp] y_poplog = y_sim_poplog;
}
With these initial conditions, the output is not exactly the same for the three models:
Could anyone explain me why there is an issue with the logtransformed model in Stan but not in R? Is there something about how rounding is done in Stan that causes this issue? When I print the ODE output for every step, I see that the E and R compartment reach very high values in only a few time steps and result in inf values and nan in a next step. Is this due to the division by exp(E) and exp(R)?
I would be very keen on getting a better understanding of Stan and how the ODE solver in Stan is different from the one in R. Any pointers would be greatly appreciated.
Best regards,
Judith
This is the code I used to run and analyse the stan models:
# Load stan file that contains these functions with same parameters as used above.
# define parameters for the spline
fixed_pars = list(incubation_period = (5/7),
infectious_period=(2/7),
contact=7*11)
simul_pars = list(I0_raw_sim=1,
beta_sim=0.08,
ts=1:45,
num_t=45,
popsize=10000)
data_list = c(fixed_pars,simul_pars)
data_list = within(data_list,{
rtol = 1e-5;
atol = 1e-5;
max_num_steps = 1000})
# compile
mod = cmdstan_model(stan_file = "/various_tests/MRE_logtransform.stan")
samples_prior = mod$sample(
data = data_list,
seed = 123,
chains = 1, parallel_chains = 1,
iter_warmup = 1, iter_sampling = 1,
fixed_param=TRUE)
y_traditional <- samples_prior$summary(c("y_traditional"))
y_log <- samples_prior$summary(c("y_log"))
y_poplog <- samples_prior$summary(c("y_poplog"))
y_anthony <- samples_prior$summary(c("y_anthony"))
y_traditional_forPlot <- y_traditional[91:135,]
y_traditional_forPlot$t <- 1:45
y_log_forPlot <- y_log[91:135,]
y_log_forPlot$t <- 1:45
y_poplog_forPlot <- y_poplog[91:135,]
y_poplog_forPlot$t <- 1:45
y_anthony_forPlot <- y_anthony[91:135,]
y_anthony_forPlot$t <- 1:45
# compare results of the three models
ggplot() + geom_line(aes(x=t, y=mean, colour="Traditional model"), data = y_traditional_forPlot) +
geom_line(aes(x=t, y= exp(mean), colour="Logtransformed model" ), data = y_log_forPlot) +
geom_line(aes(x=t, y= exp(mean)/10000, colour="Logtransformed population scaled"), data = y_poplog_forPlot) +
geom_line(aes(x=t, y= exp(mean), colour="Anthony special"), data = y_anthony_forPlot) +
xlab("time (weeks)")+ylab("Number of infected individuals") + labs(title="Results stan model") +
scale_colour_manual(values=c("Traditional model"="red", "Logtransformed model"="blue",
"Logtransformed population scaled"="violet", "Anthony special"="green")) +
theme(legend.position = "bottom")
