Use ODE in STAN with the simplest example.

Since there seems to be almost no relevant simple example for using ODE ("integrate_ode" function) in STAN, then I write it here. And it would allow one to write a code just like writing in NONMEM with ADVAN6 or ADVAN13. Let's light it on FIRE!

The simplest 1-compartment model

Consider the following 1-compartment model with log normal error. Let\(\hat{C}(t)\) be the concentration at time \(t\) and \(\displaystyle k\left(:=\frac{CL}{V}\right)\) be the elimination constant, then we have: \begin{align} \frac{\partial \hat{C}(t)}{\partial t}=-k\hat{C}(t) \label{eq:ode} \end{align}

Data generation

Analytical solution for Eq.(\ref{eq:ode}) is: \begin{align} \hat{C}(t)=\frac{D}{V}\exp \left( -kt\right), \label{eq:anasol} \end{align} where \(D\) is amount of dose (bioavailability is not considered here.) and \(V\) is volume of distribution. We can generate simulation dataset using Eq.(\ref{eq:anasol}) as follows: \begin{align} C(t)\sim \mbox{lognormal} \left( \log(\hat{C}(t)), \sigma \right), \label{eq:gene} \end{align} where \(\sigma\) is standard deviation of the measurement error. Suppose we measured from \(t=1.0\) to \(24.0\) by \(0.5\) hour and \(D=20\), \(V=30\), \(CL=6\) and \(\sigma=0.05\).

D = 20
theta=c(30,6)
t = seq(1,24,0.5)
sigma = 0.05

FOSD = function(theta,t){
  V  = theta[1]
  k  = theta[2]/theta[1]
  res = D/V*(exp(-k*t))
  return(res)
}

SCFOSD = function(theta,t,sigma){
  yhat = FOSD(theta,t)
  res = yhat*exp(rnorm(length(yhat),0,sigma))
  return(res)
}

res = c()
for (i in 1:length(t)){
  res = c(res,SCFOSD(theta,t[i],sigma))
}

p = ggplot(data.frame(t, C=res))
p + geom_point(aes(x = t, y = C))

The generated dataset is shown as below.

Implementation in STAN

The following STAN code is using the function "integrate_ode_rk45" which allows us to apply ODE expression without having analytical solution for it. The part of "transformed data" is a kind of mantra.

library(rstan)

model_code_1_comp='
functions {
  real[] onecom(real t,
                real[] y,
                real[] theta,
                real[] x_r,
                int[] x_i) {
    real dydt[1];
    real k;
    k = theta[2]/theta[1];
    dydt[1] = -k*y[1];
    return dydt;
  }
}

data {
  int T;
  real t0;
  real ts[T];
  real y[T];
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters{
  real sigma;
  real theta[2]; // V, CL
}

transformed parameters{
  real y0[1];
  real y_hat[T, 1];

  y0[1] = 20/theta[1];
  y_hat = integrate_ode_rk45(onecom, y0, t0, ts, theta, x_r, x_i);
}

model {
  for (i in 1:T)
    y[i] ~ lognormal(log(y_hat[i]),sigma);
}
'

model_1_comp = stan_model(model_code = model_code_1_comp)

Then, we can do sampling by the following script.

fit = sampling(model_1_comp,
                 data = list(T = length(t),
                             ts = t,
                             t0 = 0,
                             y = res),
                 chains = 1, iter = 1000
)

We get the following estimates and all the parameters are almost accurately recovered.

> print(fit)
Inference for Stan model: c11033a128464b6cb6a9e92a6927784c.
1 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=500.

              mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
sigma         0.05    0.00 0.00   0.04   0.04   0.05   0.05   0.06   227 1.00
theta[1]     30.83    0.04 0.46  29.90  30.55  30.83  31.13  31.71   141 1.00
theta[2]      6.11    0.01 0.07   5.99   6.07   6.11   6.16   6.24   135 1.00
y0[1]         0.65    0.00 0.01   0.63   0.64   0.65   0.65   0.67   142 1.00
y_hat[1,1]    0.53    0.00 0.01   0.52   0.53   0.53   0.54   0.55   140 1.00
y_hat[2,1]    0.48    0.00 0.01   0.47   0.48   0.48   0.49   0.50   139 1.00
y_hat[3,1]    0.44    0.00 0.01   0.43   0.43   0.44   0.44   0.45   139 1.00
y_hat[4,1]    0.40    0.00 0.01   0.39   0.39   0.40   0.40   0.41   138 1.00
                                 (snip)
y_hat[42,1]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500 1.00
y_hat[43,1]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500 1.00
y_hat[44,1]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500 1.00
y_hat[45,1]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500 1.00
y_hat[46,1]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500 1.00
y_hat[47,1]   0.01    0.00 0.00   0.01   0.01   0.01   0.01   0.01   500 1.00
lp__        124.84    0.13 1.28 121.57 124.31 125.11 125.74 126.29    98 1.02

Note that even when we changed the value of \(\sigma\) to larger value, almost accurate estimates are obtained. Perhaps it would be better to sample the posterior distribution and draw them as a Posterior Predictive Check by adding "generated quantity" part in STAN code.