The purpose of this package is to allow users to simulate and estimate the parameters from a system of ordinary differential equations (ODE) in Stan while not having to write Stan code. Events are supported by providing an event data frame.
Unfortunately, due to the user-defined nature ODE system it is not possible for the Stan model to be precompiled. Thus some compile time will be involve when simulating/fitting the model specified by the user.
NOTE: estimating parameters from the ODE (i.e. sampling = TRUE
) is not yet supported.
Load the stanode package and the deSolve package (for comparison).
library(rstanode) library(deSolve) rstan::rstan_options(auto_write = TRUE)
The mathematical representation of a simple harmonic oscillator ODE is, $$ \begin{align} dy_1 &= y_2 \ dy_2 &= -y_1 -\theta \cdot y_2 \end{align} $$
In deSolve we can fit a simple harmonic oscillator as follows,
# specify the ODE sho <- function(t, y, p) { with(as.list(c(y,p)), { dy1 = y2 dy2 = - y1 - theta * y2 return(list(c(dy1 = dy1, dy2 = dy2))) }) } # declare the parameters, initial conditions, and time steps pars <- c("theta" = 0.15) yini <- c("y1" = 1, "y2" = 0) time_steps <- seq(1, 100, by = 0.001) sho_deSolve <- ode(func = sho, y = yini, times = time_steps, parms = pars)
In stanode we can use the same function arguments to simulate from the same model,
sho_stanode <- stan_ode(sho, state = yini, pars = pars, times = time_steps, integrator = "rk45", sampling = FALSE) sims <- sho_stanode$simulations
plot(sho_deSolve[,2], sho_deSolve[,3], col = "#808080", lwd = 3, type = "l", main = "Simple Harmonic Oscillator", xlab = "y1", ylab = "y2") lines(sims[,2], sims[,3], col = "#FF6688") legend("bottomleft", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))
The stanode package also supports events. Specifically, the user provides an event table to which describes the time the event takes place, the value that each state should take on at the specfied state, and the method to use (currently "add"
and "multiply"
are supported).
Below we fit the two compartment PKPD model in deSolve
two_cpt <- function(t, y, parms) { with(as.list(c(y, parms)), { dy_gut = -ka * y_gut dy_cent = ka * y_gut - (CL/V_cent + Q/V_cent) * y_cent + (Q/V_peri) * y_peri dy_peri = (Q/V_cent) * y_cent - (Q/V_peri) * y_peri return(list(c(dy_gut=dy_gut, dy_cent=dy_cent, dy_peri=dy_peri))) }) } dosing_table <- data.frame(var = "y_gut", time = seq(10, 70, by=10), value = 5, method = "add") pars <- c("CL" = 10, Q = 13, "V_cent" = 20, "V_peri" = 73, ka = 3) yini <- c("y_gut" = 0, "y_cent" = 0, "y_peri" = 0) time_steps <- seq(0, 150, by = 0.005) two_cpt_deSolve <- deSolve::ode(y = yini, times = time_steps, func = two_cpt, parms = pars, events = list(data = dosing_table))
To fit the model in stanode we need to specify the event table,
dosing_table_stan <- data.frame(time = dosing_table$time, "y_gut" = 5, "y_cent" = 0, "y_peri" = 0, method = "add") print(dosing_table_stan)
Now we can simulate the model,
two_cpt_stanode <- stan_ode(two_cpt, state = yini, pars = pars, times = time_steps, events = dosing_table_stan, integrator = "rk45", sampling = FALSE) sims <- two_cpt_stanode$simulations
Below we plot the simulations for each compartment from both packages.
par(mfrow = c(1,3)) plot(two_cpt_deSolve[,1], two_cpt_deSolve[,2], col = "#808080", lwd = 3, type = "l", main = "Gut Compartment", xlab = "Time", ylab = "Concentration") lines(sims[,1], sims[,2], col = "#FF6688") legend("topright", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1)) plot(two_cpt_deSolve[,1], two_cpt_deSolve[,3], col = "#808080", lwd = 3, type = "l", main = "Central Compartment", xlab = "Time", ylab = "Concentration") lines(sims[,1], sims[,3], col = "#FF6688") legend("topright", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1)) plot(two_cpt_deSolve[,1], two_cpt_deSolve[,4], col = "#808080", lwd = 3, type = "l", main = "Peripheral Compartment", xlab = "Time", ylab = "Concentration") lines(sims[,1], sims[,4], col = "#FF6688") legend("topright", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))
Arenstorf <- function(t, y, p) { with(as.list(c(y,p)), { D1 <- ((y[1] + mu1)^2 + y[2]^2)^(3/2) D2 <- ((y[1] - mu2)^2 + y[2]^2)^(3/2) dy1 <- y[3] dy2 <- y[4] dy3 <- y[1] + 2*y[4] - mu2*(y[1]+mu1)/D1 - mu1*(y[1]-mu2)/D2 dy4 <- y[2] - 2*y[3] - mu2*y[2]/D1 - mu1*y[2]/D2 return(list( c(dy1, dy2, dy3, dy4) )) }) }
Simulate in deSolve.
mu1 <- 0.012277471 pars <- c(mu1 = mu1, mu2 = 1 - mu1) yini <- c(y1 = 0.994, y2 = 0, y3 = 0, y4 = -2.00158510637908252240537862224) time_steps <- seq(from = 0, to = 18, by = 0.01) orbit_deSolve <- ode(func = Arenstorf, y = yini, parms = pars, times = time_steps)
Simulate in rstanode.
orbit_rstanode <- stan_ode(Arenstorf, yini, pars, time_steps, integrator = "rk45") sims <- orbit_rstanode$simulations
Overlay both simulations.
plot(orbit_deSolve[,2], orbit_deSolve[,3], type = "l", col = "#808080", lwd = 3) lines(sims[,2], sims[,3], col = "#FF6688") legend("bottomleft", c("deSolve","stanode"), col = c("#808080", "#FF6688"), lwd = c(3, 1))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.