Merge pull request #185 from stan-dev/convert-odes

Added model code and data to convert-odes case study (design-doc #19)

Author: bbbales2

knitr/convert-odes/convert_odes.Rmd

@@ -1,7 +1,7 @@
 ---
 title: "Upgrading to the new ODE interface"
 author: "Ben Bales & Sebastian Weber"
-date: "26 July 2020"
+date: "27 July 2020"
 output: html_document
 ---
 
@@ -294,10 +294,10 @@ transformed parameters {
 ### Full Model
 
 The full model with the new interface is
-[here](https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/sir/sir.stan) and the data is [here](https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/sir/sir.data.R).
+[here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.stan) and the data is [here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.data.R).
 
 The full model with the old interface is
-[here](https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/sir/sir.stan), and the data is [here](https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/sir/sir.data.R).
+[here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.old.stan) and the data is [here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.old.data.R).
 
 ## PKPD Model
 
@@ -401,8 +401,8 @@ transformed parameters {
 ### Full Model
 
 The full model with the new interface is
-[here](https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/one_comp_mm_elim_abs/one_comp_mm_elim_abs.stan) and the data is [here](https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/one_comp_mm_elim_abs/one_comp_mm_elim_abs.data.R).
+[here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.stan) and the data is [here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.data.R).
 
 The full model with the old interface is
-[here](https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/pkpd/one_comp_mm_elim_abs.stan) and the data is [here](https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/pkpd/one_comp_mm_elim_abs.data.R).
+[here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.old.stan) and the data is [here](https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.old.data.R).
 

knitr/convert-odes/convert_odes.html

@@ -11,7 +11,7 @@
 
 <meta name="author" content="Ben Bales &amp; Sebastian Weber" />
 
-<meta name="date" content="2020-07-26" />
+<meta name="date" content="2020-07-27" />
 
 <title>Upgrading to the new ODE interface</title>
 
@@ -365,7 +365,7 @@
 
 <h1 class="title toc-ignore">Upgrading to the new ODE interface</h1>
 <h4 class="author">Ben Bales &amp; Sebastian Weber</h4>
-<h4 class="date">26 July 2020</h4>
+<h4 class="date">27 July 2020</h4>
 
 </div>
 
@@ -518,8 +518,8 @@ <h3>Calling the ODE Solver</h3>
 </div>
 <div id="full-model" class="section level3">
 <h3>Full Model</h3>
-<p>The full model with the new interface is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/sir/sir.stan">here</a> and the data is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/sir/sir.data.R">here</a>.</p>
-<p>The full model with the old interface is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/sir/sir.stan">here</a>, and the data is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/sir/sir.data.R">here</a>.</p>
+<p>The full model with the new interface is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.stan">here</a> and the data is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.data.R">here</a>.</p>
+<p>The full model with the old interface is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.old.stan">here</a> and the data is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/sir.old.data.R">here</a>.</p>
 </div>
 </div>
 <div id="pkpd-model" class="section level2">
@@ -602,8 +602,8 @@ <h3>Calling the ODE Solver</h3>
 </div>
 <div id="full-model-1" class="section level3">
 <h3>Full Model</h3>
-<p>The full model with the new interface is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/one_comp_mm_elim_abs/one_comp_mm_elim_abs.stan">here</a> and the data is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/master/benchmarks/one_comp_mm_elim_abs/one_comp_mm_elim_abs.data.R">here</a>.</p>
-<p>The full model with the old interface is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/pkpd/one_comp_mm_elim_abs.stan">here</a> and the data is <a href="https://github.com/stan-dev/stat_comp_benchmarks/blob/64b6597cd17de77e91dcfb1a9977758407017d20/benchmarks/pkpd/one_comp_mm_elim_abs.data.R">here</a>.</p>
+<p>The full model with the new interface is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.stan">here</a> and the data is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.data.R">here</a>.</p>
+<p>The full model with the old interface is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.old.stan">here</a> and the data is <a href="https://github.com/stan-dev/example-models/blob/master/knitr/convert-odes/one_comp_mm_elim_abs.old.data.R">here</a>.</p>
 </div>
 </div>
 </div>

knitr/convert-odes/one_comp_mm_elim_abs.data.R

@@ -0,0 +1,18 @@
+t0 <- 0
+C0 <- 
+structure(c(0),
+.Dim = c(1))
+D <- 30
+V <- 2
+times <- 
+c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6,
+6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10)
+N_t <- 20
+C <- 
+c(5.70812264865215, 7.10126075072086,
+8.38520678426651, 9.79008249883381, 13.4390409239245,
+11.4478987597702, 11.2124282696837, 11.4269217682577,
+12.2432859438401, 13.8201804108938, 13.8408670746042,
+11.422291744251, 10.5031943081843, 11.9121452242965,
+14.0849980781312, 10.5505145523917, 10.1905539351877,
+12.2232272590821, 11.7290653047821, 12.2719396535996)

knitr/convert-odes/one_comp_mm_elim_abs.old.data.R

@@ -0,0 +1,18 @@
+t0 <- 0
+C0 <- 
+structure(c(0),
+.Dim = c(1))
+D <- 30
+V <- 2
+times <- 
+c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6,
+6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10)
+N_t <- 20
+C_hat <- 
+c(5.70812264865215, 7.10126075072086,
+8.38520678426651, 9.79008249883381, 13.4390409239245,
+11.4478987597702, 11.2124282696837, 11.4269217682577,
+12.2432859438401, 13.8201804108938, 13.8408670746042,
+11.422291744251, 10.5031943081843, 11.9121452242965,
+14.0849980781312, 10.5505145523917, 10.1905539351877,
+12.2232272590821, 11.7290653047821, 12.2719396535996)

knitr/convert-odes/one_comp_mm_elim_abs.old.stan

@@ -0,0 +1,75 @@
+functions {
+  real[] one_comp_mm_elim_abs(real t,
+                              real[] y,
+                              real[] theta,
+                              real[] x_r,
+                              int[] x_i) {
+    real dydt[1];
+    real k_a = theta[1]; // Dosing rate in 1/day
+    real K_m = theta[2]; // Michaelis-Menten constant in mg/L
+    real V_m = theta[3]; // Maximum elimination rate in 1/day
+    real D = x_r[1];
+    real V = x_r[2];
+    real dose = 0;
+    real elim = (V_m / V) * y[1] / (K_m + y[1]);
+
+    if (t > 0)
+      dose = exp(- k_a * t) * D * k_a / V;
+
+    dydt[1] = dose - elim;
+
+    return dydt;
+  }
+}
+
+data {
+  real t0;    // Initial time in days;
+  real C0[1]; // Initial concentration at t0 in mg/L
+
+  real D;   // Total dosage in mg
+  real V;   // Compartment volume in L
+
+  int<lower=1> N_t;
+  real times[N_t];   // Measurement times in days
+
+  // Measured concentrations in effect compartment in mg/L
+  real C_hat[N_t];
+}
+
+transformed data {
+  real x_r[2] = {D, V};
+  int x_i[0];
+}
+
+parameters {
+  real<lower=0> k_a; // Dosing rate in 1/day
+  real<lower=0> K_m; // Michaelis-Menten constant in mg/L
+  real<lower=0> V_m; // Maximum elimination rate in 1/day
+  real<lower=0> sigma;
+}
+
+transformed parameters {
+  real C[N_t, 1];
+  {
+    real theta[3] = {k_a, K_m, V_m};
+    C = integrate_ode_bdf(one_comp_mm_elim_abs, C0, t0, times, theta, x_r, x_i);
+  }
+}
+
+model {
+  // Priors
+  k_a ~ cauchy(0, 1);
+  K_m ~ cauchy(0, 1);
+  V_m ~ cauchy(0, 1);
+  sigma ~ cauchy(0, 1);
+
+  // Likelihood
+  for (n in 1:N_t)
+    C_hat[n] ~ lognormal(log(C[n, 1]), sigma);
+}
+
+generated quantities {
+  real C_ppc[N_t];
+  for (n in 1:N_t)
+    C_ppc[n] = lognormal_rng(log(C[n, 1]), sigma);
+}

knitr/convert-odes/one_comp_mm_elim_abs.stan

@@ -0,0 +1,65 @@
+functions {
+  vector one_comp_mm_elim_abs(real t,
+                              vector y,
+                              real k_a, // Dosing rate in 1/day
+                              real K_m, // Michaelis-Menten constant in mg/L
+                              real V_m, // Maximum elimination rate in 1/day
+                              real D,
+                              real V) {
+    vector[1] dydt;
+
+    real dose = 0;
+    real elim = (V_m / V) * y[1] / (K_m + y[1]);
+
+    if (t > 0)
+      dose = exp(- k_a * t) * D * k_a / V;
+
+    dydt[1] = dose - elim;
+
+    return dydt;
+  }
+}
+
+data {
+  real t0;      // Initial time in days;
+  vector[1] C0; // Initial concentration at t0 in mg/L
+
+  real D;   // Total dosage in mg
+  real V;   // Compartment volume in L
+
+  int<lower=1> N_t;
+  real times[N_t];   // Measurement times in days
+
+  // Measured concentrations in effect compartment in mg/L
+  real C[N_t];
+}
+
+parameters {
+  real<lower=0> k_a; // Dosing rate in 1/day
+  real<lower=0> K_m; // Michaelis-Menten constant in mg/L
+  real<lower=0> V_m; // Maximum elimination rate in 1/day
+  real<lower=0> sigma;
+}
+
+transformed parameters {
+  vector[1] mu_C[N_t] = ode_bdf_tol(one_comp_mm_elim_abs, C0, t0, times,
+				 1e-8, 1e-8, 1000,
+				 k_a, K_m, V_m, D, V);
+}
+
+model {
+  // Priors
+  k_a ~ cauchy(0, 1);
+  K_m ~ cauchy(0, 1);
+  V_m ~ cauchy(0, 1);
+  sigma ~ cauchy(0, 1);
+
+  // Likelihood
+  C ~ lognormal(log(mu_C[, 1]), sigma);
+}
+
+generated quantities {
+  real C_ppc[N_t];
+  for (n in 1:N_t)
+    C_ppc[n] = lognormal_rng(log(mu_C[n, 1]), sigma);
+}

knitr/convert-odes/sir.data.R

@@ -0,0 +1,9 @@
+N_t <- 20
+
+t <- c(7,14,21,28,35,42,49,56,63,70,77,84,91,98,105,112,119,126,133,140)
+
+y0 <- c(10000,0,0,10000)
+
+stoi <- c(441,419,553,666,753,692,644,509,371,286,200,147,101,68,49,35,19,12,7,10)
+
+B <- c(6141.73,5801.72,9439.14,13742,11975,16387,15144,14705.6,8736.39,7553.29,6007.09,5979.96,3048.17,1566.77,1275.21,824.85,604.605,516.25,258.046,188.591)