From 45f93c8abec5fa57142f12e9897c42b7ef44f0ad Mon Sep 17 00:00:00 2001 From: Imad Ali Date: Tue, 15 Oct 2019 18:03:26 -0500 Subject: [PATCH 1/2] adding ode forcing function and event examples --- knitr/ode-examples/data/two_comp.data.R | 13 + knitr/ode-examples/data/two_comp_e.data.R | 20 + knitr/ode-examples/data/two_comp_ef.data.R | 20 + .../ode-examples/data/two_comp_ef_fit.data.R | 50 + knitr/ode-examples/data/two_comp_f.data.R | 13 + knitr/ode-examples/models/one_comp.stan | 39 + knitr/ode-examples/models/two_comp.stan | 30 + knitr/ode-examples/models/two_comp_e.stan | 42 + knitr/ode-examples/models/two_comp_ef.stan | 64 ++ .../ode-examples/models/two_comp_ef_fit.stan | 71 ++ knitr/ode-examples/models/two_comp_f.stan | 52 ++ knitr/ode-examples/ode_examples.Rmd | 437 +++++++++ knitr/ode-examples/ode_examples.html | 884 ++++++++++++++++++ 13 files changed, 1735 insertions(+) create mode 100644 knitr/ode-examples/data/two_comp.data.R create mode 100644 knitr/ode-examples/data/two_comp_e.data.R create mode 100644 knitr/ode-examples/data/two_comp_ef.data.R create mode 100644 knitr/ode-examples/data/two_comp_ef_fit.data.R create mode 100644 knitr/ode-examples/data/two_comp_f.data.R create mode 100644 knitr/ode-examples/models/one_comp.stan create mode 100644 knitr/ode-examples/models/two_comp.stan create mode 100644 knitr/ode-examples/models/two_comp_e.stan create mode 100644 knitr/ode-examples/models/two_comp_ef.stan create mode 100644 knitr/ode-examples/models/two_comp_ef_fit.stan create mode 100644 knitr/ode-examples/models/two_comp_f.stan create mode 100644 knitr/ode-examples/ode_examples.Rmd create mode 100644 knitr/ode-examples/ode_examples.html diff --git a/knitr/ode-examples/data/two_comp.data.R b/knitr/ode-examples/data/two_comp.data.R new file mode 100644 index 000000000..f4975ce97 --- /dev/null +++ b/knitr/ode-examples/data/two_comp.data.R @@ -0,0 +1,13 @@ +K <- 2 +S <- 2 +T <- 91 +t0 <- 0.9 +theta <- +c(6, 0.6) +ts <- +c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, +3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, +6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, +8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10) +y0 <- +c(3, 6) diff --git a/knitr/ode-examples/data/two_comp_e.data.R b/knitr/ode-examples/data/two_comp_e.data.R new file mode 100644 index 000000000..7efe7091a --- /dev/null +++ b/knitr/ode-examples/data/two_comp_e.data.R @@ -0,0 +1,20 @@ +E <- 9 +end <- +c(11, 21, 31, 41, 51, 61, 71, 81, 91) +K <- 2 +S <- 2 +start <- +c(1, 12, 22, 32, 42, 52, 62, 72, 82) +T <- 91 +t0 <- +c(0.9, 2, 3, 4, 5, 6, 7, 8, 9) +theta <- +c(6, 0.6) +ts <- +c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, +3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, +6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, +8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10) +y0 <- +structure(c(0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0), +.Dim = c(9, 2)) diff --git a/knitr/ode-examples/data/two_comp_ef.data.R b/knitr/ode-examples/data/two_comp_ef.data.R new file mode 100644 index 000000000..7efe7091a --- /dev/null +++ b/knitr/ode-examples/data/two_comp_ef.data.R @@ -0,0 +1,20 @@ +E <- 9 +end <- +c(11, 21, 31, 41, 51, 61, 71, 81, 91) +K <- 2 +S <- 2 +start <- +c(1, 12, 22, 32, 42, 52, 62, 72, 82) +T <- 91 +t0 <- +c(0.9, 2, 3, 4, 5, 6, 7, 8, 9) +theta <- +c(6, 0.6) +ts <- +c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, +3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, +6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, +8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10) +y0 <- +structure(c(0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0), +.Dim = c(9, 2)) diff --git a/knitr/ode-examples/data/two_comp_ef_fit.data.R b/knitr/ode-examples/data/two_comp_ef_fit.data.R new file mode 100644 index 000000000..0f36b4959 --- /dev/null +++ b/knitr/ode-examples/data/two_comp_ef_fit.data.R @@ -0,0 +1,50 @@ +E <- 9 +end <- +c(11, 21, 31, 41, 51, 61, 71, 81, 91) +K <- 2 +S <- 2 +start <- +c(1, 12, 22, 32, 42, 52, 62, 72, 82) +T <- 91 +t0 <- +c(0.9, 2, 3, 4, 5, 6, 7, 8, 9) +theta <- +c(6, 0.6) +ts <- +c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, +3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, +5.8, 5.9, 6, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, +8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10) +y <- +structure(c(0.0501581870130331, 0.171390081674708, 0.154594268549219, 0.062065259109517, 0.0456829424379411, 0, +0.208754305578589, 0, 0.0789583981577013, 0, 0, 1.13180612537173, 0.488747295719004, 0.175938234736721, +0.267220370371945, 0.0596273820158646, 0.114899295074526, 0.272005811828668, 0.00321200420427842, 0.0389646648821001, 0, +1.20316038685346, 0.904118171078214, 0.98223423254061, 0.565295590488372, 0.522927186207073, 0.450513348921578, +0.620365646294086, 0.4678699568707, 0.37057962028136, 0.555095261235546, 1.74251402715741, 1.14647265763887, +0.908613731853189, 0.784366309313526, 0.679373169916111, 0.417099037784891, 0.653246649036365, 0.675125974398131, +0.716075219103793, 0.278019100183284, 1.35021531291007, 0.759361286445196, 0.465848106411765, 0.303854355878931, 0, +0.044578906469163, 0, 0, 0, 0, 1.08801243790365, 0.559402702889078, 0.365565008029735, 0.161785136396892, +0.168557908358781, 0.0127705580341626, 0.04979082649047, 0.0881241905230884, 0, 0, 1.49806956362223, 0.884173967451606, +0.86524922688302, 0.502676374667655, 0.562566837439883, 0.476620165921046, 0.519539509440848, 0.218662655660154, +0.507152267513044, 0.513108897994652, 1.49996791821863, 1.17971436361413, 0.772735263706271, 0.630423972389203, +0.601900948517316, 0.582514140732464, 0.377620333755696, 0.635375982811837, 0.490251481506544, 0.569144292781708, +1.46665189789528, 0.651212389778649, 0.314880821297867, 0.133348004125639, 0.117135004199271, 0.0739408633489209, 0, +0.0789415707251588, 0.0717642520375456, 0, 0, 0, 0, 0, 0.0138889574457143, 0.0699245986925465, 0, 0.0317266826234932, +0.0563205469874709, 0.051618034181338, 0, 0.841733488377438, 1.25733816758221, 1.29762986984564, 1.43653036429298, +1.49264448556212, 1.47832947165913, 1.39076145170894, 1.44881690940123, 1.29443317998131, 1.1183347031059, +2.05961490452341, 2.565206781713, 2.95672006361366, 3.17636228719576, 3.25413343388872, 3.5490595350775, +3.51759017530972, 3.6994981311587, 3.57849480197177, 3.86912766768025, 4.8210520239222, 5.31095419533127, +5.47286924999716, 5.63455729766036, 5.6922340334613, 5.66667261852888, 5.89035701971403, 5.69128413842573, +5.5653802812262, 5.73401806024937, 6.25614629211007, 6.54052967860694, 6.40081961001448, 6.17014388043333, +6.16265043518996, 5.46368006613957, 5.56019481866796, 5.20963902503991, 4.8578719579305, 4.650638197639, +5.05538464031898, 5.1948775633154, 5.3505287015876, 5.1999262717203, 4.96108694353612, 4.7424579450011, 4.25759754692586, +4.23333121655404, 4.02914711154008, 3.59490771430027, 4.50205709132745, 4.92517194331467, 5.12484869975415, +5.12002491226594, 5.2107712706115, 5.06519828158445, 5.30513662266828, 5.21038342568496, 5.17183165177851, +5.29201340365431, 6.1766443939395, 6.63712849299992, 6.49368187291553, 6.84264585701854, 6.71851118465657, +6.71572385763204, 6.6640121268618, 6.43922293809072, 6.40336702821227, 6.10826396009771, 7.10335309546733, +7.36477065392934, 7.0881024310677, 6.8154777231462, 6.61436329940234, 6.11434070682372, 5.93348503591031, +5.46484119601089, 5.33824209235697, 5.20892459607444), +.Dim = c(91, 2)) +y0 <- +structure(c(0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0), +.Dim = c(9, 2)) diff --git a/knitr/ode-examples/data/two_comp_f.data.R b/knitr/ode-examples/data/two_comp_f.data.R new file mode 100644 index 000000000..f4975ce97 --- /dev/null +++ b/knitr/ode-examples/data/two_comp_f.data.R @@ -0,0 +1,13 @@ +K <- 2 +S <- 2 +T <- 91 +t0 <- 0.9 +theta <- +c(6, 0.6) +ts <- +c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, +3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, +6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, +8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10) +y0 <- +c(3, 6) diff --git a/knitr/ode-examples/models/one_comp.stan b/knitr/ode-examples/models/one_comp.stan new file mode 100644 index 000000000..4d146a776 --- /dev/null +++ b/knitr/ode-examples/models/one_comp.stan @@ -0,0 +1,39 @@ +functions { + real[] one_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) { + real dydt[1]; + dydt[1] = -theta[1] * y[1]; + return dydt; + } +} +data { + int T; + int E; + int S; + real y0[E,S]; + real t0[E]; + real ts[T]; + real theta[1]; + int start[E]; + int end[E]; +} +transformed data { + real x[0]; + int x_int[0]; +} +model { +} +generated quantities { + real y_hat[T,S]; + // 1. stop the integrator at t + // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state) + // 3. combine y0 with the event at t + // 4. restart the integrator + for (i in 1:E) { + real y_state[S]; + if (i == 1) + y_state = y0[1,]; + else + y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],])); + y_hat[start[i]:end[i],] = integrate_ode_bdf(one_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int); + } +} diff --git a/knitr/ode-examples/models/two_comp.stan b/knitr/ode-examples/models/two_comp.stan new file mode 100644 index 000000000..570c26a34 --- /dev/null +++ b/knitr/ode-examples/models/two_comp.stan @@ -0,0 +1,30 @@ +// two compartment (sim) +functions { + real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) { + real dydt[2]; + real a = theta[1]; + real b = theta[2]; + dydt[1] = -a*y[1]; + dydt[2] = a*y[1] - b*y[2]; + return dydt; + } +} +data { + int T; + int S; + int K; + real y0[S]; + real t0; + real ts[T]; + real theta[K]; +} +transformed data { + real x[0]; + int x_int[0]; +} +model { +} +generated quantities { + real y_hat[T,S]; + y_hat = integrate_ode_bdf(two_comp, y0, t0, ts, theta, x, x_int); +} diff --git a/knitr/ode-examples/models/two_comp_e.stan b/knitr/ode-examples/models/two_comp_e.stan new file mode 100644 index 000000000..08d62d955 --- /dev/null +++ b/knitr/ode-examples/models/two_comp_e.stan @@ -0,0 +1,42 @@ +// two compartment event (sim) +functions { + real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) { + real dydt[2]; + dydt[1] = -theta[1] * y[1]; + dydt[2] = theta[1] * y[1] - theta[2] * y[2]; + return dydt; + } +} +data { + int T; + int E; + int S; + int K; + real y0[E,S]; + real t0[E]; + real ts[T]; + real theta[K]; + int start[E]; + int end[E]; +} +transformed data { + real x[0]; + int x_int[0]; +} +model { +} +generated quantities { + real y_hat[T,S]; + // 1. stop the integrator at t + // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state) + // 3. combine y0 with the event at t + // 4. restart the integrator + for (i in 1:E) { + real y_state[S]; + if (i == 1) + y_state = y0[1,]; + else + y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],])); + y_hat[start[i]:end[i],] = integrate_ode_bdf(two_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int); + } +} diff --git a/knitr/ode-examples/models/two_comp_ef.stan b/knitr/ode-examples/models/two_comp_ef.stan new file mode 100644 index 000000000..c223e47b4 --- /dev/null +++ b/knitr/ode-examples/models/two_comp_ef.stan @@ -0,0 +1,64 @@ +// two compartment event and forcing (sim) +functions { + int max_indx (real t, real[] x) { + int N = num_elements(x); + int indx = 1; + while (indx <= N && x[indx] <= t) { + indx = indx + 1; + } + indx = indx - 1; + return indx; + } + + real force_constant (real t, real[] x_step, real[] y_step) { + int indx; + real value; + indx = max_indx(t, x_step); + if (indx == 0) + value = 0; + else + value = y_step[indx]; + return value; + } + + real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) { + real dydt[2]; + real c = force_constant(t, {0.0,3,5,7,9}, {0.0,3,0,3,0}); + dydt[1] = -theta[1] * y[1] + c; + dydt[2] = theta[1] * y[1] - theta[2] * y[2]; + return dydt; + } +} +data { + int T; + int E; + int S; + int K; + real y0[E,S]; + real t0[E]; + real ts[T]; + real theta[K]; + int start[E]; + int end[E]; +} +transformed data { + real x[0]; + int x_int[0]; +} +model { +} +generated quantities { + real y_hat[T,S]; + // 1. stop the integrator at t + // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state) + // 3. combine y0 with the event at t + // 4. restart the integrator + for (i in 1:E) { + real y_state[S]; + if (i == 1) + y_state = y0[1,]; + else + y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],])); + y_hat[start[i]:end[i],] = integrate_ode_bdf(two_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int); + } +} diff --git a/knitr/ode-examples/models/two_comp_ef_fit.stan b/knitr/ode-examples/models/two_comp_ef_fit.stan new file mode 100644 index 000000000..53bdfdc18 --- /dev/null +++ b/knitr/ode-examples/models/two_comp_ef_fit.stan @@ -0,0 +1,71 @@ +// two compartment event and forcing (fit) +functions { + int max_indx (real t, real[] x) { + int N = num_elements(x); + int indx = 1; + while (indx <= N && x[indx] <= t) { + indx = indx + 1; + } + indx = indx - 1; + return indx; + } + + real force_constant (real t, real[] x_step, real[] y_step) { + int indx; + real value; + indx = max_indx(t, x_step); + if (indx == 0) + value = 0; + else + value = y_step[indx]; + return value; + } + + real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) { + real dydt[2]; + real c = force_constant(t, {0.0,3,5,7,9}, {0.0,3,0,3,0}); + dydt[1] = -theta[1] * y[1] + c; + dydt[2] = theta[1] * y[1] - theta[2] * y[2]; + return dydt; + } +} +data { + int T; + int E; + int S; + int K; + real y0[E,S]; + real t0[E]; + real ts[T]; + int start[E]; + int end[E]; + real y[T,S]; +} +transformed data { + real x[0]; + int x_int[0]; +} +parameters { + real theta[K]; + real sigma; +} +model { + real y_hat[T,S]; + // 1. stop the integrator at t + // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state) + // 3. combine y0 with the event at t + // 4. restart the integrator + for (i in 1:E) { + real y_state[S]; + if (i == 1) + y_state = y0[1,]; + else + y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],])); + y_hat[start[i]:end[i],] = integrate_ode_bdf(two_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int); + } + for (t in 1:T) + target+= normal_lpdf(y[t] | y_hat[t], sigma); + target+= normal_lpdf(theta[1] | 5, 1); + target+= normal_lpdf(theta[2] | 0, 1); + target+= normal_lpdf(sigma | 0, 0.1); +} diff --git a/knitr/ode-examples/models/two_comp_f.stan b/knitr/ode-examples/models/two_comp_f.stan new file mode 100644 index 000000000..105772a8f --- /dev/null +++ b/knitr/ode-examples/models/two_comp_f.stan @@ -0,0 +1,52 @@ +// two compartment forcing (sim) +functions { + int max_indx (real t, real[] x) { + int N = num_elements(x); + int indx = 1; + while (indx <= N && x[indx] <= t) { + indx = indx + 1; + } + indx = indx - 1; + return indx; + } + + real force_constant (real t, real[] x_step, real[] y_step) { + int indx; + real value; + indx = max_indx(t, x_step); + if (indx == 0) + value = 0; + else + value = y_step[indx]; + return value; + } + + real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) { + real dydt[2]; + real a = theta[1]; + real b = theta[2]; + real c = force_constant(t, {0.0,3,5,7,9}, {0.0,3,0,3,0}); + dydt[1] = -a*y[1] + c; + dydt[2] = a*y[1] - b*y[2]; + return dydt; + } +} +data { + int T; + int K; + int S; + real y0[S]; + real t0; + real ts[T]; + real theta[K]; +} +transformed data { + real x[0]; + int x_int[0]; +} +model { +} +generated quantities { + real y_hat[T,2]; + y_hat = integrate_ode_bdf(two_comp, y0, t0, ts, theta, x, x_int); +} diff --git a/knitr/ode-examples/ode_examples.Rmd b/knitr/ode-examples/ode_examples.Rmd new file mode 100644 index 000000000..7c6dd9e8b --- /dev/null +++ b/knitr/ode-examples/ode_examples.Rmd @@ -0,0 +1,437 @@ +--- +title: "ODEs with Forcing Functions and Events in Stan" +author: "Imad Ali" +date: "10/7/2019" +output: + html_document: + theme: sandstone + toc: TRUE +--- + +```{r setup, include=FALSE} +library(deSolve) +library(rstan) +library(bayesplot) +rstan_options(auto_write = TRUE) +knitr::opts_chunk$set(echo = TRUE) +``` + +## Introduction + +The purpose of this set of examples is to go over how to implement events and forcing functions in ordinary differential equations (ODEs) using [RStan](https://mc-stan.org/rstan/). Recall that an ODE is an equation that relates a function of an independent variable to the derivative of that function. The solution of an ODE involves finding the function that satifies the equation. + +An ODE can have external variables and events that impact the system and therefore the solution. A **forcing function** is a function that allows an external variable to influence the model as a function of time. Typically, a series of points define a forcing function so the evaluation of the function at each time step defined by the ODE needs to be interpolated (e.g. constant or linear). An **event** occurs when the dependent variable suddenly changes (e.g. a value is added/multiplied to the dependent variable). This causes a discontinuity in the variable's sequence. + +We will consider the following ODE from Soetaert et al (2012). This is a stylized example from pharmacokinetics where the concentration of a drug is modeled in the intestines and blood. The first equation $\frac{dy_1}{dt}$ governs the elimination of a drug from the intestines. The second equation $\frac{dy_2}{dt}$ governs the intake of that drug from the intestines and into the blood and the subsequent elimination from the blood. The ODE is defined as, + +$$ +\begin{align*} +\frac{dy_1}{dt} &= -ay_1 \\ +\frac{dy_2}{dt} &= ay_1 - by_2\\ +\end{align*} +$$ + +where the independent variable is $t$ and the parameters $a$ and $b$ are known. + +In the sections that follow we will find numerical solutions to the ODE above, + +1. without any event or forcing function +2. with a forcing function that defines a continuous external effect at various time intervals +3. with a set of events that define the uptake of the drug +4. with both the forcing function defined in (2) and the events defined in (3). + +In the final section we will also show how to treat the ODE as a density where we know the numerical solution and would like to infer the parameters that generated that solution. + +In each case we will find the solution using the R package [deSolve](https://cran.r-project.org/web/packages/deSolve/index.html) (a popular package used to numerically solve ODEs) and then find the equivalent solution using the integrators provided in the Stan language. We use the same example throughout so that the user can clearly understand the changes made to the Stan program when implementing forcing functions and events, as it is not trivial. (Note that by means do we claim that the use of forcing functions and events in this exposition is appropriate given this or any pharmacokinetic model.) + + + +## No Event or Forcing Function + +Here we model the ODE as is using the initial values $y_1 = 3, y_2 = 6$ and parameters $a = 6, b = 0.6$. Within the context of the pharmacokinitic model one could assume that some concentration of the drug exists in the patient and we would like to simulate its elimination. The solution using deSolve is provided below. + +```{r ode, fig.align='center', fig.width=10, fig.height=5} +# time sequence +times <- seq(from = 1, to = 10, by = 0.1) +times_desolve <- c(times[1] - (times[2]-times[1]), times) +# parameters +pars <- list(a = 6, b = 0.6) +# initial state values +yini <- c(Intestine = 3, Blood = 6) +# ODE +two_comp_fun <- function(t, y, p) { + list2env(p, environment()) + dy1 <- -a*y[1] + dy2 <- a*y[1] - b*y[2] + return(list(c(dy1,dy2))) +} + +sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf") +plot(sim_desolve) +``` + +The comparable solution using Stan is computed below. For details see the Stan program defined in `models/two_comp.stan`. + +```{r} +stan_data <- list(theta = unlist(pars), + ts = times, + T = length(times), + K = length(unlist(pars)), + S = length(yini), + y0 = unname(yini), + t0 = times_desolve[1]) +tc_stan <- stan("models/two_comp.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1) +``` + +Below we overlay the solution using Stan on top of the solution defined by deSolve. + +```{r, fig.align='center', fig.width=10, fig.height=5} +sim_tc_stan <- as.matrix(tc_stan) +sim_tc_stan <- cbind(time = times, + y1 = unname(sim_tc_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_stan))]), + y2 = unname(sim_tc_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_stan))])) + +par(mfrow = c(1,2)) +plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4, + main = "Intestine", + xlab = "time", ylab = "") +lines(sim_tc_stan[,"time"], sim_tc_stan[,"y1"], col = "#CCFF00") +legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4, + main = "Blood", + xlab = "time", ylab = "") +lines(sim_tc_stan[,"time"], sim_tc_stan[,"y2"], col = "#CCFF00") +legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +``` + + + +## Forcing Function + +Now we assume that we have a forcing function that introduces a constant external effect from time steps 3 to 5 and 7 to 9 that increases the concentration in the intestines by 3 (otherwise 0). The solution in deSolve is provided below and is acheived by defining the forcing function using the `approxfun` function provided in R. This function allows us to perform either constant of linear interpolation. + +```{r ode_force, fig.align='center', fig.width=10, fig.height=5} +# forcing function (constant interpolation) +force <- approxfun(c(0,3,5,7,9), c(0,3,0,3,0), method = "constant", rule = 2) +# ODE with forcing function +two_comp_fun <- function(t, y, p) { + list2env(p, environment()) + dy1 <- -a*y[1] + force(t) + dy2 <- a*y[1] - b*y[2] + return(list(c(dy1,dy2))) +} + +sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf") +plot(sim_desolve) +``` + +The solution in Stan is computed below. Note that in order to achieve comparable results we had to implement an identical version of the forcing function in the `functions {}` block in Stan. For details see the Stan program `models/two_comp_f.stan`. + +```{r} +# create stan data +stan_data <- list(theta = unlist(pars), + ts = times, + T = length(times), + K = length(unlist(pars)), + S = length(yini), + y0 = yini, + t0 = times_desolve[1]) +tc_force_stan <- stan("models/two_comp_f.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1) +``` + +Below we overlay the solution using Stan on top of the solution defined by deSolve. + +```{r, fig.align='center', fig.width=10, fig.height=5} +sim_tc_force_stan <- as.matrix(tc_force_stan) +sim_tc_force_stan <- cbind(time = times, + y1 = unname(sim_tc_force_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_force_stan))]), + y2 = unname(sim_tc_force_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_force_stan))])) + +par(mfrow = c(1,2)) +plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4, + main = "Intestine", + xlab = "time", ylab = "") +lines(sim_tc_force_stan[,"time"], sim_tc_force_stan[,"y1"], col = "#CCFF00") +legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4, + main = "Blood", + xlab = "time", ylab = "") +lines(sim_tc_force_stan[,"time"], sim_tc_force_stan[,"y2"], col = "#CCFF00") +legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +``` + + + +## Event + +Now assume that there is no concentration of the drug in the patient at the initial time $t_0$ by specifying the initial state values as $y_1 = 0$ and $y_2 = 0$. Additionally a dose event takes place at each time step in the sequence $2, 3, 4, \ldots, 9$. The dose increases the concentration of the drug in the intestines by a value of 2. The solution in deSolve is provided below and is accomplished by passing an event table to the `ode` function. + +```{r ode_event, fig.align='center', fig.width=10, fig.height=5} +yini <- c(Intestine = 0, Blood = 0) +# event table +dose <- data.frame(var = "Intestine", + time = seq(2,9,by=1), + value = 2, + method = "add") +# ODE +two_comp_fun <- function(t, y, p) { + list2env(p, environment()) + dy1 <- -a*y[1] + dy2 <- a*y[1] - b*y[2] + return(list(c(dy1,dy2))) +} + +sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf", events = list(data = dose)) +plot(sim_desolve) +``` + +The solution in Stan is computed below. The way events are captured in the ODE defined is to stop the integrator at each time step defined by the event, incorporate the event (in this case via the addition operation), and then start the integrator again. This is done in the Stan program `models/two_comp_e.stan`. + +```{r} +# break up time into event segments +seg <- dose$time +seg_indxs <- which(times %in% seg) + +if (tail(seg_indxs,1) != length(times)) { + seg_indxs <- c(seg_indxs, length(times)) +} + +if (seg_indxs[1] != 1) { + seg_indxs <- c(1, seg_indxs) +} + +# for readability create a Stan event table +seg_start <- head(seg_indxs, -1) + 1 +seg_start[1] <- 1 +seg_end <- tail(seg_indxs, -1) +seg_value <- cbind(dose$value, 0) +seg_value <- rbind(yini, seg_value) + +# create stan data +stan_data <- list(theta = unlist(pars), + ts = times, + T = length(times), + E = nrow(seg_value), + S = ncol(seg_value), + K = length(unlist(pars)), + y0 = seg_value, + t0 = times[seg_start] - (times[2]-times[1]), + start = seg_start, + end = seg_end) + +# simulate from ode +tc_event_stan <- stan("models/two_comp_e.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1) + +``` + +Below we overlay the solution using Stan on top of the solution defined by deSolve. + +```{r, fig.align='center', fig.width=10, fig.height=5} +sim_tc_event_stan <- as.matrix(tc_event_stan) +sim_tc_event_stan <- cbind(time = times, + y1 = unname(sim_tc_event_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_event_stan))]), + y2 = unname(sim_tc_event_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_event_stan))])) + +par(mfrow = c(1,2)) + +plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4, + main = "Intestine", + xlab = "time", ylab = "") +lines(sim_tc_event_stan[,"time"], sim_tc_event_stan[,"y1"], col = "#CCFF00") +legend("topleft", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4, + main = "Blood", + xlab = "time", ylab = "") +lines(sim_tc_event_stan[,"time"], sim_tc_event_stan[,"y2"], col = "#CCFF00") +legend("bottomright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +``` + + + +## Event and Forcing Function + +Now we include both the event and forcing function described in the previous sections. The deSolve solution is provided below. + +```{r ode_event_force, fig.align='center', fig.width=10, fig.height=5} +# forcing function +force <- approxfun(c(0,3,5,7,9), c(0,3,0,3,0), method = "constant", rule = 2) +# event table +dose <- data.frame(var = "Intestine", + time = seq(2,9,by=1), + value = 2, + method = "add") +# ODE +two_comp_fun <- function(t, y, p) { + list2env(p, environment()) + dy1 <- -a*y[1] + force(t) + dy2 <- a*y[1] - b*y[2] + return(list(c(dy1,dy2))) +} + +sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf", events = list(data = dose)) +plot(sim_desolve) +``` + +The comparable solution using Stan is provided below. For details see the Stan program `models/two_comp_ef.stan`. + +```{r} +# break up time into event segments +seg <- dose$time +seg_indxs <- which(times %in% seg) + +if (tail(seg_indxs,1) != length(times)) { + seg_indxs <- c(seg_indxs, length(times)) +} + +if (seg_indxs[1] != 1) { + seg_indxs <- c(1, seg_indxs) +} + +# for readability create a Stan event table +seg_start <- head(seg_indxs, -1) + 1 +seg_start[1] <- 1 +seg_end <- tail(seg_indxs, -1) +seg_value <- cbind(dose$value, 0) +seg_value <- rbind(yini, seg_value) + +# create stan data +stan_data <- list(theta = unlist(pars), + ts = times, + T = length(times), + E = nrow(seg_value), + S = ncol(seg_value), + K = length(unlist(pars)), + y0 = seg_value, + t0 = times[seg_start] - (times[2]-times[1]), + start = seg_start, + end = seg_end) + +# simulate from ode +tc_event_force_stan <- stan("models/two_comp_ef.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1) +``` + +Below we overlay the solution using Stan on top of the solution defined by deSolve. + +```{r, fig.align='center', fig.width=10, fig.height=5} +sim_tc_event_force_stan <- as.matrix(tc_event_force_stan) +sim_tc_event_force_stan <- cbind(time = times, + y1 = unname(sim_tc_event_force_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_event_force_stan))]), + y2 = unname(sim_tc_event_force_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_event_force_stan))])) + +par(mfrow = c(1,2)) +plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4, + main = "Intestine", + xlab = "time", ylab = "") +lines(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y1"], col = "#CCFF00") +legend("topleft", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4, + main = "Blood", + xlab = "time", ylab = "") +lines(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y2"], col = "#CCFF00") +legend("bottomright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8) +``` + + + +## Estimation + +All the examples above involve simulation and assume that we know the value of the parameters. Alternatively, we can use Stan to infer the parameters used to generate the numerical solution (once we have observed the solution). This means that we can treat the observed solution as our outcome variable and infer the parameters $a$ and $b$ that generated the solution. + +As an example we use the ODE defined in the previous section along with its numerical solution computed in Stan. To resemble real data with measurment error we add some arbitrary noise to this solution (this could be done in the simulation itself within the `generated quantities {}` block). Below we overlay the true solution along with the solution with noise. + +```{r ode_event_force_fit, fig.width=10, fig.height=5} +# add noise to the solution +y <- sim_tc_event_force_stan +y[,2] <- y[,2] + rnorm(length(times), 0, 0.1) +y[,3] <- y[,3] + rnorm(length(times), 0, 0.1) + +y[,2] <- ifelse(y[,2] < 0, 0, y[,2]) +y[,3] <- ifelse(y[,3] < 0, 0, y[,3]) + +par(mfrow = c(1,2)) +plot(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y1"], type = "l", lwd = 4, + main = "Intestine", + xlab = "time", ylab = "") +lines(y[,"time"], y[,"y1"], col = "#FF6688") +legend("topleft", c("Stan","Stan with noise"), col = c("#000000","#FF6688"), pch = c(15,15), pt.cex = 2, cex = 0.8) + +plot(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y2"], type = "l", lwd = 4, + main = "Blood", + xlab = "time", ylab = "") +lines(y[,"time"], y[,"y2"], col = "#FF6688") +legend("bottomright", c("Stan","Stan with noise"), col = c("#000000","#FF6688"), pch = c(15,15), pt.cex = 2, cex = 0.8) +``` + +To estimate the parameters $a$ and $b$ we move the integrator specification from the `generated quantities {}` block to the `model {}` block. In addition to the integrator we incorporate the following pieces into the `model {}` block of the Stan program: + +**Likelihood** +$$ +y \sim \mathcal{N}(\hat{y}, \sigma) +$$ + +**Priors** +$$ +\begin{align*} +\sigma &\sim \mathcal{N}(0, 0.1) \\ +a &\sim \mathcal{N}(5, 1) \\ +b &\sim \mathcal{N}(0, 1) \\ +\end{align*} +$$ + +Where $y$ is the data we acquired from simulation in the previous section (with noise added) and $\hat{y}$ is the simulation computed by the ODE integrator in the Stan program. Essentially we are modeling $y$ according to the deterministic process that generates $\hat{y}$ along with some measurement error. For details on how this is implemented see the Stan program `models/two_comp_ef_fit.stan`. Below we fit the model to the data. + +```{r ode_event_force_fit_cache, results='hide', warning=FALSE, cache=TRUE} +stan_data$y <- y[,-1] +tc_event_force_fit <- stan("models/two_comp_ef_fit.stan", data = stan_data, chains = 4, iter = 500, cores = 4) +``` + +Looking at the parameter values and diagnostics we find that we appropriately recovered the parameters that generated the data. + +```{r} +print(tc_event_force_fit) +``` + +```{r, fig.align='center', fig.width=15, fig.height=5} +mcmc_trace(as.array(tc_event_force_fit), regex_pars = "theta|sigma") +``` + +Using these results the user can continue with useful Bayesian post-estimation steps, such as generating posterior predictions and constructing credible intervals. + + + +## References + +Gabry, J. & Mahr, T. (2019) _bayesplot: Plotting for Bayesian Models_. R package version 1.7.0, mc-stan.org/bayesplot. + +Mattuck, A., Miller, H., Orloff, J., & Lewis, J. (2011) _Differential Equations_. Massachusetts Institute of Technology: MIT OpenCourseWare. https://ocw.mit.edu. + +Soetaert, K., Cash, J., & Mazzia F. (2012) _Solving Differential Equations in R_. Use R! + +Soetaert, K., Petzoldt, T., & Setzer, R.W. (2010). Solving Differential Equations in R: Package deSolve. _Journal of Statistical Software_. 33(9), 1-25. http://www.jstatsoft.org/v33/i09/. + +Stan Development Team (2019) _Stan Modeling Language User's Guide_. Version 2.20. http://mc-stan.org. + +Stan Development Team (2019) _RStan: the R interface to Stan_. R package version 2.19.1, http://mc-stan.org/. + + + +## Computing Environment + +```{r} +devtools::session_info("rstan") +``` + + + +## Licences + +Code (c) 2019 Imad Ali, BSD 3. + +Text (c) 2019 Imad Ali, licensed under CC-BY-NC 4.0. \ No newline at end of file diff --git a/knitr/ode-examples/ode_examples.html b/knitr/ode-examples/ode_examples.html new file mode 100644 index 000000000..3d3b5f204 --- /dev/null +++ b/knitr/ode-examples/ode_examples.html @@ -0,0 +1,884 @@ + + + + + + + + + + + + + + + +ODEs with Forcing Functions and Events in Stan + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + +
+

Introduction

+

The purpose of this set of examples is to go over how to implement events and forcing functions in ordinary differential equations (ODEs) using RStan. Recall that an ODE is an equation that relates a function of an independent variable to the derivative of that function. The solution of an ODE involves finding the function that satifies the equation.

+

An ODE can have external variables and events that impact the system and therefore the solution. A forcing function is a function that allows an external variable to influence the model as a function of time. Typically, a series of points define a forcing function so the evaluation of the function at each time step defined by the ODE needs to be interpolated (e.g. constant or linear). An event occurs when the dependent variable suddenly changes (e.g. a value is added/multiplied to the dependent variable). This causes a discontinuity in the variable’s sequence.

+

We will consider the following ODE from Soetaert et al (2012). This is a stylized example from pharmacokinetics where the concentration of a drug is modeled in the intestines and blood. The first equation \(\frac{dy_1}{dt}\) governs the elimination of a drug from the intestines. The second equation \(\frac{dy_2}{dt}\) governs the intake of that drug from the intestines and into the blood and the subsequent elimination from the blood. The ODE is defined as,

+

\[ +\begin{align*} +\frac{dy_1}{dt} &= -ay_1 \\ +\frac{dy_2}{dt} &= ay_1 - by_2\\ +\end{align*} +\]

+

where the independent variable is \(t\) and the parameters \(a\) and \(b\) are known.

+

In the sections that follow we will find numerical solutions to the ODE above,

+
    +
  1. without any event or forcing function
    2. +
  2. with a forcing function that defines a continuous external effect at various time intervals
    3. +
  3. with a set of events that define the uptake of the drug
    4. +
  4. with both the forcing function defined in (2) and the events defined in (3).
    5. +
+

In the final section we will also show how to treat the ODE as a density where we know the numerical solution and would like to infer the parameters that generated that solution.

+

In each case we will find the solution using the R package deSolve (a popular package used to numerically solve ODEs) and then find the equivalent solution using the integrators provided in the Stan language. We use the same example throughout so that the user can clearly understand the changes made to the Stan program when implementing forcing functions and events, as it is not trivial. (Note that by means do we claim that the use of forcing functions and events in this exposition is appropriate given this or any pharmacokinetic model.)

+
+
+

No Event or Forcing Function

+

Here we model the ODE as is using the initial values \(y_1 = 3, y_2 = 6\) and parameters \(a = 6, b = 0.6\). Within the context of the pharmacokinitic model one could assume that some concentration of the drug exists in the patient and we would like to simulate its elimination. The solution using deSolve is provided below.

+
# time sequence
+times <- seq(from = 1, to = 10, by = 0.1)
+times_desolve <- c(times[1] - (times[2]-times[1]), times)
+# parameters
+pars <- list(a = 6, b = 0.6)
+# initial state values
+yini <- c(Intestine = 3, Blood = 6)
+# ODE
+two_comp_fun <- function(t, y, p) {
+  list2env(p, environment())
+  dy1 <- -a*y[1]
+  dy2 <- a*y[1] - b*y[2]
+  return(list(c(dy1,dy2)))
+}
+
+sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf")
+plot(sim_desolve)
+

+

The comparable solution using Stan is computed below. For details see the Stan program defined in models/two_comp.stan.

+
stan_data <- list(theta = unlist(pars),
+                  ts = times,
+                  T = length(times),
+                  K = length(unlist(pars)),
+                  S = length(yini),
+                  y0 = unname(yini),
+                  t0 = times_desolve[1])
+tc_stan <- stan("models/two_comp.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
+
## Warning in strptime(xx, f <- "%Y-%m-%d %H:%M:%OS", tz = tz): unknown
+## timezone 'zone/tz/2019b.1.0/zoneinfo/America/Chicago'
+
## 
+## SAMPLING FOR MODEL 'two_comp' NOW (CHAIN 1).
+## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
+## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
+## Chain 1:                0.000701 seconds (Sampling)
+## Chain 1:                0.000701 seconds (Total)
+

Below we overlay the solution using Stan on top of the solution defined by deSolve.

+
sim_tc_stan <- as.matrix(tc_stan)
+sim_tc_stan <- cbind(time = times,
+                           y1 = unname(sim_tc_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_stan))]),
+                           y2 = unname(sim_tc_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_stan))]))
+
+par(mfrow = c(1,2))
+plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
+     main = "Intestine",
+     xlab = "time", ylab = "")
+lines(sim_tc_stan[,"time"], sim_tc_stan[,"y1"], col = "#CCFF00")
+legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+
+plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
+     main = "Blood",
+     xlab = "time", ylab = "")
+lines(sim_tc_stan[,"time"], sim_tc_stan[,"y2"], col = "#CCFF00")
+legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+

+
+
+

Forcing Function

+

Now we assume that we have a forcing function that introduces a constant external effect from time steps 3 to 5 and 7 to 9 that increases the concentration in the intestines by 3 (otherwise 0). The solution in deSolve is provided below and is acheived by defining the forcing function using the approxfun function provided in R. This function allows us to perform either constant of linear interpolation.

+
# forcing function (constant interpolation)
+force <- approxfun(c(0,3,5,7,9), c(0,3,0,3,0), method = "constant", rule = 2)
+# ODE with forcing function
+two_comp_fun <- function(t, y, p) {
+  list2env(p, environment())
+  dy1 <- -a*y[1] + force(t)
+  dy2 <- a*y[1] - b*y[2]
+  return(list(c(dy1,dy2)))
+}
+
+sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf")
+plot(sim_desolve)
+

+

The solution in Stan is computed below. Note that in order to achieve comparable results we had to implement an identical version of the forcing function in the functions {} block in Stan. For details see the Stan program models/two_comp_f.stan.

+
# create stan data
+stan_data <- list(theta = unlist(pars),
+                  ts = times,
+                  T = length(times),
+                  K = length(unlist(pars)),
+                  S = length(yini),
+                  y0 = yini,
+                  t0 = times_desolve[1])
+tc_force_stan <- stan("models/two_comp_f.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
+
## 
+## SAMPLING FOR MODEL 'two_comp_f' NOW (CHAIN 1).
+## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
+## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
+## Chain 1:                0.007601 seconds (Sampling)
+## Chain 1:                0.007601 seconds (Total)
+

Below we overlay the solution using Stan on top of the solution defined by deSolve.

+
sim_tc_force_stan <- as.matrix(tc_force_stan)
+sim_tc_force_stan <- cbind(time = times,
+                           y1 = unname(sim_tc_force_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_force_stan))]),
+                           y2 = unname(sim_tc_force_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_force_stan))]))
+
+par(mfrow = c(1,2))
+plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
+     main = "Intestine",
+     xlab = "time", ylab = "")
+lines(sim_tc_force_stan[,"time"], sim_tc_force_stan[,"y1"], col = "#CCFF00")
+legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+
+plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
+     main = "Blood",
+     xlab = "time", ylab = "")
+lines(sim_tc_force_stan[,"time"], sim_tc_force_stan[,"y2"], col = "#CCFF00")
+legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+

+
+
+

Event

+

Now assume that there is no concentration of the drug in the patient at the initial time \(t_0\) by specifying the initial state values as \(y_1 = 0\) and \(y_2 = 0\). Additionally a dose event takes place at each time step in the sequence \(2, 3, 4, \ldots, 9\). The dose increases the concentration of the drug in the intestines by a value of 2. The solution in deSolve is provided below and is accomplished by passing an event table to the ode function.

+
yini <- c(Intestine = 0, Blood = 0)
+# event table
+dose <- data.frame(var = "Intestine",
+                   time = seq(2,9,by=1),
+                   value = 2,
+                   method = "add")
+# ODE
+two_comp_fun <- function(t, y, p) {
+  list2env(p, environment())
+  dy1 <- -a*y[1]
+  dy2 <- a*y[1] - b*y[2]
+  return(list(c(dy1,dy2)))
+}
+
+sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf", events = list(data = dose))
+plot(sim_desolve)
+

+

The solution in Stan is computed below. The way events are captured in the ODE defined is to stop the integrator at each time step defined by the event, incorporate the event (in this case via the addition operation), and then start the integrator again. This is done in the Stan program models/two_comp_e.stan.

+
# break up time into event segments
+seg <- dose$time
+seg_indxs <- which(times %in% seg)
+
+if (tail(seg_indxs,1) != length(times)) {
+  seg_indxs <- c(seg_indxs, length(times))
+}
+
+if (seg_indxs[1] != 1) {
+  seg_indxs <- c(1, seg_indxs)
+}
+
+# for readability create a Stan event table
+seg_start <- head(seg_indxs, -1) + 1
+seg_start[1] <- 1
+seg_end <- tail(seg_indxs, -1) 
+seg_value <- cbind(dose$value, 0)
+seg_value <- rbind(yini, seg_value)
+
+# create stan data
+stan_data <- list(theta = unlist(pars),
+                  ts = times,
+                  T = length(times),
+                  E = nrow(seg_value),
+                  S = ncol(seg_value),
+                  K = length(unlist(pars)),
+                  y0 = seg_value,
+                  t0 = times[seg_start] - (times[2]-times[1]),
+                  start = seg_start,
+                  end = seg_end)
+
+# simulate from ode
+tc_event_stan <- stan("models/two_comp_e.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
+
## 
+## SAMPLING FOR MODEL 'two_comp_e' NOW (CHAIN 1).
+## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
+## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
+## Chain 1:                0.001852 seconds (Sampling)
+## Chain 1:                0.001852 seconds (Total)
+

Below we overlay the solution using Stan on top of the solution defined by deSolve.

+
sim_tc_event_stan <- as.matrix(tc_event_stan)
+sim_tc_event_stan <- cbind(time = times,
+                           y1 = unname(sim_tc_event_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_event_stan))]),
+                           y2 = unname(sim_tc_event_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_event_stan))]))
+
+par(mfrow = c(1,2))
+
+plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
+     main = "Intestine",
+     xlab = "time", ylab = "")
+lines(sim_tc_event_stan[,"time"], sim_tc_event_stan[,"y1"], col = "#CCFF00")
+legend("topleft", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+
+plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
+     main = "Blood",
+     xlab = "time", ylab = "")
+lines(sim_tc_event_stan[,"time"], sim_tc_event_stan[,"y2"], col = "#CCFF00")
+legend("bottomright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+

+
+
+

Event and Forcing Function

+

Now we include both the event and forcing function described in the previous sections. The deSolve solution is provided below.

+
# forcing function
+force <- approxfun(c(0,3,5,7,9), c(0,3,0,3,0), method = "constant", rule = 2)
+# event table
+dose <- data.frame(var = "Intestine",
+                   time = seq(2,9,by=1),
+                   value = 2,
+                   method = "add")
+# ODE
+two_comp_fun <- function(t, y, p) {
+  list2env(p, environment())
+  dy1 <- -a*y[1] + force(t)
+  dy2 <- a*y[1] - b*y[2]
+  return(list(c(dy1,dy2)))
+}
+
+sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf", events = list(data = dose))
+plot(sim_desolve)
+

+

The comparable solution using Stan is provided below. For details see the Stan program models/two_comp_ef.stan.

+
# break up time into event segments
+seg <- dose$time
+seg_indxs <- which(times %in% seg)
+
+if (tail(seg_indxs,1) != length(times)) {
+  seg_indxs <- c(seg_indxs, length(times))
+}
+
+if (seg_indxs[1] != 1) {
+  seg_indxs <- c(1, seg_indxs)
+}
+
+# for readability create a Stan event table
+seg_start <- head(seg_indxs, -1) + 1
+seg_start[1] <- 1
+seg_end <- tail(seg_indxs, -1) 
+seg_value <- cbind(dose$value, 0)
+seg_value <- rbind(yini, seg_value)
+
+# create stan data
+stan_data <- list(theta = unlist(pars),
+                  ts = times,
+                  T = length(times),
+                  E = nrow(seg_value),
+                  S = ncol(seg_value),
+                  K = length(unlist(pars)),
+                  y0 = seg_value,
+                  t0 = times[seg_start] - (times[2]-times[1]),
+                  start = seg_start,
+                  end = seg_end)
+
+# simulate from ode
+tc_event_force_stan <- stan("models/two_comp_ef.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
+
## 
+## SAMPLING FOR MODEL 'two_comp_ef' NOW (CHAIN 1).
+## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
+## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
+## Chain 1:                0.009013 seconds (Sampling)
+## Chain 1:                0.009013 seconds (Total)
+

Below we overlay the solution using Stan on top of the solution defined by deSolve.

+
sim_tc_event_force_stan <- as.matrix(tc_event_force_stan)
+sim_tc_event_force_stan <- cbind(time = times,
+                           y1 = unname(sim_tc_event_force_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_event_force_stan))]),
+                           y2 = unname(sim_tc_event_force_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_event_force_stan))]))
+
+par(mfrow = c(1,2))
+plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
+     main = "Intestine",
+     xlab = "time", ylab = "")
+lines(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y1"], col = "#CCFF00")
+legend("topleft", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+
+plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
+     main = "Blood",
+     xlab = "time", ylab = "")
+lines(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y2"], col = "#CCFF00")
+legend("bottomright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+

+
+
+

Estimation

+

All the examples above involve simulation and assume that we know the value of the parameters. Alternatively, we can use Stan to infer the parameters used to generate the numerical solution (once we have observed the solution). This means that we can treat the observed solution as our outcome variable and infer the parameters \(a\) and \(b\) that generated the solution.

+

As an example we use the ODE defined in the previous section along with its numerical solution computed in Stan. To resemble real data with measurment error we add some arbitrary noise to this solution (this could be done in the simulation itself within the generated quantities {} block). Below we overlay the true solution along with the solution with noise.

+
# add noise to the solution
+y <- sim_tc_event_force_stan
+y[,2] <- y[,2] + rnorm(length(times), 0, 0.1)
+y[,3] <- y[,3] + rnorm(length(times), 0, 0.1)
+
+y[,2] <- ifelse(y[,2] < 0, 0, y[,2])
+y[,3] <- ifelse(y[,3] < 0, 0, y[,3])
+
+par(mfrow = c(1,2))
+plot(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y1"], type = "l", lwd = 4,
+     main = "Intestine",
+     xlab = "time", ylab = "")
+lines(y[,"time"], y[,"y1"], col = "#FF6688")
+legend("topleft", c("Stan","Stan with noise"), col = c("#000000","#FF6688"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+
+plot(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y2"], type = "l", lwd = 4,
+     main = "Blood",
+     xlab = "time", ylab = "")
+lines(y[,"time"], y[,"y2"], col = "#FF6688")
+legend("bottomright", c("Stan","Stan with noise"), col = c("#000000","#FF6688"), pch = c(15,15), pt.cex = 2, cex = 0.8)
+

+

To estimate the parameters \(a\) and \(b\) we move the integrator specification from the generated quantities {} block to the model {} block. In addition to the integrator we incorporate the following pieces into the model {} block of the Stan program:

+

Likelihood \[ +y \sim \mathcal{N}(\hat{y}, \sigma) +\]

+

Priors \[ +\begin{align*} +\sigma &\sim \mathcal{N}(0, 0.1) \\ +a &\sim \mathcal{N}(5, 1) \\ +b &\sim \mathcal{N}(0, 1) \\ +\end{align*} +\]

+

Where \(y\) is the data we acquired from simulation in the previous section (with noise added) and \(\hat{y}\) is the simulation computed by the ODE integrator in the Stan program. Essentially we are modeling \(y\) according to the deterministic process that generates \(\hat{y}\) along with some measurement error. For details on how this is implemented see the Stan program models/two_comp_ef_fit.stan. Below we fit the model to the data.

+
stan_data$y <- y[,-1]
+tc_event_force_fit <- stan("models/two_comp_ef_fit.stan", data = stan_data, chains = 4, iter = 500, cores = 4)
+

Looking at the parameter values and diagnostics we find that we appropriately recovered the parameters that generated the data.

+
print(tc_event_force_fit)
+
## Inference for Stan model: two_comp_ef_fit.
+## 4 chains, each with iter=500; warmup=250; thin=1; 
+## post-warmup draws per chain=250, total post-warmup draws=1000.
+## 
+##            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
+## theta[1]   6.07    0.00 0.08   5.91   6.02   6.07   6.13   6.21   405 1.00
+## theta[2]   0.60    0.00 0.00   0.60   0.60   0.60   0.60   0.60  1044 1.00
+## sigma      0.08    0.00 0.00   0.08   0.08   0.08   0.09   0.09   264 1.01
+## lp__     186.79    0.08 1.30 183.42 186.19 187.14 187.74 188.31   270 1.01
+## 
+## Samples were drawn using NUTS(diag_e) at Tue Oct 15 16:11:05 2019.
+## For each parameter, n_eff is a crude measure of effective sample size,
+## and Rhat is the potential scale reduction factor on split chains (at 
+## convergence, Rhat=1).
+
mcmc_trace(as.array(tc_event_force_fit), regex_pars = "theta|sigma")
+

+

Using these results the user can continue with useful Bayesian post-estimation steps, such as generating posterior predictions and constructing credible intervals.

+
+
+

References

+

Gabry, J. & Mahr, T. (2019) bayesplot: Plotting for Bayesian Models. R package version 1.7.0, mc-stan.org/bayesplot.

+

Mattuck, A., Miller, H., Orloff, J., & Lewis, J. (2011) Differential Equations. Massachusetts Institute of Technology: MIT OpenCourseWare. https://ocw.mit.edu.

+

Soetaert, K., Cash, J., & Mazzia F. (2012) Solving Differential Equations in R. Use R!

+

Soetaert, K., Petzoldt, T., & Setzer, R.W. (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software. 33(9), 1-25. http://www.jstatsoft.org/v33/i09/.

+

Stan Development Team (2019) Stan Modeling Language User’s Guide. Version 2.20. http://mc-stan.org.

+

Stan Development Team (2019) RStan: the R interface to Stan. R package version 2.19.1, http://mc-stan.org/.

+
+
+

Computing Environment

+
devtools::session_info("rstan")
+
## Session info -------------------------------------------------------------
+
##  setting  value                       
+##  version  R version 3.4.0 (2017-04-21)
+##  system   x86_64, darwin15.6.0        
+##  ui       X11                         
+##  language (EN)                        
+##  collate  en_US.UTF-8                 
+##  tz       <NA>                        
+##  date     2019-10-15
+
## Packages -----------------------------------------------------------------
+
##  package      * version    date       source                         
+##  assertthat     0.2.0      2017-04-11 CRAN (R 3.4.0)                 
+##  backports      1.1.2      2017-12-13 cran (@1.1.2)                  
+##  base64enc      0.1-3      2015-07-28 CRAN (R 3.4.0)                 
+##  BH             1.66.0-1   2018-02-13 cran (@1.66.0-)                
+##  callr          3.0.0.9000 2018-09-06 Github (r-pkgs/callr@10c1269)  
+##  cli            1.0.1      2018-09-25 cran (@1.0.1)                  
+##  colorspace     1.3-2      2016-12-14 CRAN (R 3.4.0)                 
+##  compiler       3.4.0      2017-04-21 local                          
+##  crayon         1.3.4      2017-09-16 cran (@1.3.4)                  
+##  desc           1.2.0      2018-05-01 cran (@1.2.0)                  
+##  digest         0.6.18     2018-10-10 cran (@0.6.18)                 
+##  ggplot2      * 3.0.0      2018-07-03 cran (@3.0.0)                  
+##  glue           1.3.0      2018-07-17 cran (@1.3.0)                  
+##  graphics     * 3.4.0      2017-04-21 local                          
+##  grDevices    * 3.4.0      2017-04-21 local                          
+##  grid           3.4.0      2017-04-21 local                          
+##  gridExtra      2.3        2017-09-09 cran (@2.3)                    
+##  gtable         0.2.0      2016-02-26 CRAN (R 3.4.0)                 
+##  inline         0.3.15     2018-05-18 cran (@0.3.15)                 
+##  labeling       0.3        2014-08-23 CRAN (R 3.4.0)                 
+##  lattice        0.20-35    2017-03-25 CRAN (R 3.4.0)                 
+##  lazyeval       0.2.1      2017-10-29 cran (@0.2.1)                  
+##  loo            2.0.0      2018-04-11 CRAN (R 3.4.4)                 
+##  magrittr       1.5        2014-11-22 CRAN (R 3.4.0)                 
+##  MASS           7.3-47     2017-02-26 CRAN (R 3.4.0)                 
+##  Matrix         1.2-14     2018-04-09 CRAN (R 3.4.4)                 
+##  matrixStats    0.52.2     2017-04-14 CRAN (R 3.4.0)                 
+##  methods      * 3.4.0      2017-04-21 local                          
+##  mgcv           1.8-17     2017-02-08 CRAN (R 3.4.0)                 
+##  munsell        0.5.0      2018-06-12 cran (@0.5.0)                  
+##  nlme           3.1-131    2017-02-06 CRAN (R 3.4.0)                 
+##  parallel       3.4.0      2017-04-21 local                          
+##  pillar         1.1.0      2018-01-14 cran (@1.1.0)                  
+##  pkgbuild       1.0.0      2018-09-06 Github (r-lib/pkgbuild@b748041)
+##  plyr           1.8.4      2016-06-08 CRAN (R 3.4.0)                 
+##  processx       3.2.0      2018-08-16 cran (@3.2.0)                  
+##  ps             1.1.0      2018-08-10 cran (@1.1.0)                  
+##  R6             2.3.0      2018-10-04 cran (@2.3.0)                  
+##  RColorBrewer   1.1-2      2014-12-07 CRAN (R 3.4.0)                 
+##  Rcpp           1.0.0      2018-11-07 cran (@1.0.0)                  
+##  RcppEigen      0.3.3.4.0  2018-02-07 cran (@0.3.3.4)                
+##  reshape2       1.4.3      2017-12-11 cran (@1.4.3)                  
+##  rlang          0.3.1      2019-01-08 CRAN (R 3.4.4)                 
+##  rprojroot      1.3-2      2018-01-03 cran (@1.3-2)                  
+##  rstan        * 2.18.1     2018-09-07 local                          
+##  scales         1.0.0      2018-08-09 cran (@1.0.0)                  
+##  StanHeaders  * 2.18.0     2018-09-07 local                          
+##  stats        * 3.4.0      2017-04-21 local                          
+##  stats4         3.4.0      2017-04-21 local                          
+##  stringi        1.2.4      2018-07-20 cran (@1.2.4)                  
+##  stringr        1.3.1      2018-05-10 cran (@1.3.1)                  
+##  tibble         1.4.2      2018-01-22 cran (@1.4.2)                  
+##  tools          3.4.0      2017-04-21 local                          
+##  utf8           1.1.3      2018-01-03 cran (@1.1.3)                  
+##  utils        * 3.4.0      2017-04-21 local                          
+##  viridisLite    0.3.0      2018-02-01 cran (@0.3.0)                  
+##  withr          2.1.2      2018-03-15 CRAN (R 3.4.4)
+
+
+

Licences

+

Code (c) 2019 Imad Ali, BSD 3.

+

Text (c) 2019 Imad Ali, licensed under CC-BY-NC 4.0.

+
+ + + + +
+ + + + + + + + + + + + + + + From 50ba237b7e9aa8c9312414f80dfbf6bebcab894b Mon Sep 17 00:00:00 2001 From: Imad Ali Date: Tue, 15 Oct 2019 23:02:51 -0500 Subject: [PATCH 2/2] removing html render --- knitr/ode-examples/ode_examples.html | 884 --------------------------- 1 file changed, 884 deletions(-) delete mode 100644 knitr/ode-examples/ode_examples.html diff --git a/knitr/ode-examples/ode_examples.html b/knitr/ode-examples/ode_examples.html deleted file mode 100644 index 3d3b5f204..000000000 --- a/knitr/ode-examples/ode_examples.html +++ /dev/null @@ -1,884 +0,0 @@ - - - - - - - - - - - - - - - -ODEs with Forcing Functions and Events in Stan - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
-

Introduction

-

The purpose of this set of examples is to go over how to implement events and forcing functions in ordinary differential equations (ODEs) using RStan. Recall that an ODE is an equation that relates a function of an independent variable to the derivative of that function. The solution of an ODE involves finding the function that satifies the equation.

-

An ODE can have external variables and events that impact the system and therefore the solution. A forcing function is a function that allows an external variable to influence the model as a function of time. Typically, a series of points define a forcing function so the evaluation of the function at each time step defined by the ODE needs to be interpolated (e.g. constant or linear). An event occurs when the dependent variable suddenly changes (e.g. a value is added/multiplied to the dependent variable). This causes a discontinuity in the variable’s sequence.

-

We will consider the following ODE from Soetaert et al (2012). This is a stylized example from pharmacokinetics where the concentration of a drug is modeled in the intestines and blood. The first equation \(\frac{dy_1}{dt}\) governs the elimination of a drug from the intestines. The second equation \(\frac{dy_2}{dt}\) governs the intake of that drug from the intestines and into the blood and the subsequent elimination from the blood. The ODE is defined as,

-

\[ -\begin{align*} -\frac{dy_1}{dt} &= -ay_1 \\ -\frac{dy_2}{dt} &= ay_1 - by_2\\ -\end{align*} -\]

-

where the independent variable is \(t\) and the parameters \(a\) and \(b\) are known.

-

In the sections that follow we will find numerical solutions to the ODE above,

-
    -
  1. without any event or forcing function
    2. -
  2. with a forcing function that defines a continuous external effect at various time intervals
    3. -
  3. with a set of events that define the uptake of the drug
    4. -
  4. with both the forcing function defined in (2) and the events defined in (3).
    5. -
-

In the final section we will also show how to treat the ODE as a density where we know the numerical solution and would like to infer the parameters that generated that solution.

-

In each case we will find the solution using the R package deSolve (a popular package used to numerically solve ODEs) and then find the equivalent solution using the integrators provided in the Stan language. We use the same example throughout so that the user can clearly understand the changes made to the Stan program when implementing forcing functions and events, as it is not trivial. (Note that by means do we claim that the use of forcing functions and events in this exposition is appropriate given this or any pharmacokinetic model.)

-
-
-

No Event or Forcing Function

-

Here we model the ODE as is using the initial values \(y_1 = 3, y_2 = 6\) and parameters \(a = 6, b = 0.6\). Within the context of the pharmacokinitic model one could assume that some concentration of the drug exists in the patient and we would like to simulate its elimination. The solution using deSolve is provided below.

-
# time sequence
-times <- seq(from = 1, to = 10, by = 0.1)
-times_desolve <- c(times[1] - (times[2]-times[1]), times)
-# parameters
-pars <- list(a = 6, b = 0.6)
-# initial state values
-yini <- c(Intestine = 3, Blood = 6)
-# ODE
-two_comp_fun <- function(t, y, p) {
-  list2env(p, environment())
-  dy1 <- -a*y[1]
-  dy2 <- a*y[1] - b*y[2]
-  return(list(c(dy1,dy2)))
-}
-
-sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf")
-plot(sim_desolve)
-

-

The comparable solution using Stan is computed below. For details see the Stan program defined in models/two_comp.stan.

-
stan_data <- list(theta = unlist(pars),
-                  ts = times,
-                  T = length(times),
-                  K = length(unlist(pars)),
-                  S = length(yini),
-                  y0 = unname(yini),
-                  t0 = times_desolve[1])
-tc_stan <- stan("models/two_comp.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
-
## Warning in strptime(xx, f <- "%Y-%m-%d %H:%M:%OS", tz = tz): unknown
-## timezone 'zone/tz/2019b.1.0/zoneinfo/America/Chicago'
-
## 
-## SAMPLING FOR MODEL 'two_comp' NOW (CHAIN 1).
-## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
-## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
-## Chain 1:                0.000701 seconds (Sampling)
-## Chain 1:                0.000701 seconds (Total)
-

Below we overlay the solution using Stan on top of the solution defined by deSolve.

-
sim_tc_stan <- as.matrix(tc_stan)
-sim_tc_stan <- cbind(time = times,
-                           y1 = unname(sim_tc_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_stan))]),
-                           y2 = unname(sim_tc_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_stan))]))
-
-par(mfrow = c(1,2))
-plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
-     main = "Intestine",
-     xlab = "time", ylab = "")
-lines(sim_tc_stan[,"time"], sim_tc_stan[,"y1"], col = "#CCFF00")
-legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-
-plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
-     main = "Blood",
-     xlab = "time", ylab = "")
-lines(sim_tc_stan[,"time"], sim_tc_stan[,"y2"], col = "#CCFF00")
-legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-

-
-
-

Forcing Function

-

Now we assume that we have a forcing function that introduces a constant external effect from time steps 3 to 5 and 7 to 9 that increases the concentration in the intestines by 3 (otherwise 0). The solution in deSolve is provided below and is acheived by defining the forcing function using the approxfun function provided in R. This function allows us to perform either constant of linear interpolation.

-
# forcing function (constant interpolation)
-force <- approxfun(c(0,3,5,7,9), c(0,3,0,3,0), method = "constant", rule = 2)
-# ODE with forcing function
-two_comp_fun <- function(t, y, p) {
-  list2env(p, environment())
-  dy1 <- -a*y[1] + force(t)
-  dy2 <- a*y[1] - b*y[2]
-  return(list(c(dy1,dy2)))
-}
-
-sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf")
-plot(sim_desolve)
-

-

The solution in Stan is computed below. Note that in order to achieve comparable results we had to implement an identical version of the forcing function in the functions {} block in Stan. For details see the Stan program models/two_comp_f.stan.

-
# create stan data
-stan_data <- list(theta = unlist(pars),
-                  ts = times,
-                  T = length(times),
-                  K = length(unlist(pars)),
-                  S = length(yini),
-                  y0 = yini,
-                  t0 = times_desolve[1])
-tc_force_stan <- stan("models/two_comp_f.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
-
## 
-## SAMPLING FOR MODEL 'two_comp_f' NOW (CHAIN 1).
-## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
-## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
-## Chain 1:                0.007601 seconds (Sampling)
-## Chain 1:                0.007601 seconds (Total)
-

Below we overlay the solution using Stan on top of the solution defined by deSolve.

-
sim_tc_force_stan <- as.matrix(tc_force_stan)
-sim_tc_force_stan <- cbind(time = times,
-                           y1 = unname(sim_tc_force_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_force_stan))]),
-                           y2 = unname(sim_tc_force_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_force_stan))]))
-
-par(mfrow = c(1,2))
-plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
-     main = "Intestine",
-     xlab = "time", ylab = "")
-lines(sim_tc_force_stan[,"time"], sim_tc_force_stan[,"y1"], col = "#CCFF00")
-legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-
-plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
-     main = "Blood",
-     xlab = "time", ylab = "")
-lines(sim_tc_force_stan[,"time"], sim_tc_force_stan[,"y2"], col = "#CCFF00")
-legend("topright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-

-
-
-

Event

-

Now assume that there is no concentration of the drug in the patient at the initial time \(t_0\) by specifying the initial state values as \(y_1 = 0\) and \(y_2 = 0\). Additionally a dose event takes place at each time step in the sequence \(2, 3, 4, \ldots, 9\). The dose increases the concentration of the drug in the intestines by a value of 2. The solution in deSolve is provided below and is accomplished by passing an event table to the ode function.

-
yini <- c(Intestine = 0, Blood = 0)
-# event table
-dose <- data.frame(var = "Intestine",
-                   time = seq(2,9,by=1),
-                   value = 2,
-                   method = "add")
-# ODE
-two_comp_fun <- function(t, y, p) {
-  list2env(p, environment())
-  dy1 <- -a*y[1]
-  dy2 <- a*y[1] - b*y[2]
-  return(list(c(dy1,dy2)))
-}
-
-sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf", events = list(data = dose))
-plot(sim_desolve)
-

-

The solution in Stan is computed below. The way events are captured in the ODE defined is to stop the integrator at each time step defined by the event, incorporate the event (in this case via the addition operation), and then start the integrator again. This is done in the Stan program models/two_comp_e.stan.

-
# break up time into event segments
-seg <- dose$time
-seg_indxs <- which(times %in% seg)
-
-if (tail(seg_indxs,1) != length(times)) {
-  seg_indxs <- c(seg_indxs, length(times))
-}
-
-if (seg_indxs[1] != 1) {
-  seg_indxs <- c(1, seg_indxs)
-}
-
-# for readability create a Stan event table
-seg_start <- head(seg_indxs, -1) + 1
-seg_start[1] <- 1
-seg_end <- tail(seg_indxs, -1) 
-seg_value <- cbind(dose$value, 0)
-seg_value <- rbind(yini, seg_value)
-
-# create stan data
-stan_data <- list(theta = unlist(pars),
-                  ts = times,
-                  T = length(times),
-                  E = nrow(seg_value),
-                  S = ncol(seg_value),
-                  K = length(unlist(pars)),
-                  y0 = seg_value,
-                  t0 = times[seg_start] - (times[2]-times[1]),
-                  start = seg_start,
-                  end = seg_end)
-
-# simulate from ode
-tc_event_stan <- stan("models/two_comp_e.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
-
## 
-## SAMPLING FOR MODEL 'two_comp_e' NOW (CHAIN 1).
-## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
-## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
-## Chain 1:                0.001852 seconds (Sampling)
-## Chain 1:                0.001852 seconds (Total)
-

Below we overlay the solution using Stan on top of the solution defined by deSolve.

-
sim_tc_event_stan <- as.matrix(tc_event_stan)
-sim_tc_event_stan <- cbind(time = times,
-                           y1 = unname(sim_tc_event_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_event_stan))]),
-                           y2 = unname(sim_tc_event_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_event_stan))]))
-
-par(mfrow = c(1,2))
-
-plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
-     main = "Intestine",
-     xlab = "time", ylab = "")
-lines(sim_tc_event_stan[,"time"], sim_tc_event_stan[,"y1"], col = "#CCFF00")
-legend("topleft", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-
-plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
-     main = "Blood",
-     xlab = "time", ylab = "")
-lines(sim_tc_event_stan[,"time"], sim_tc_event_stan[,"y2"], col = "#CCFF00")
-legend("bottomright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-

-
-
-

Event and Forcing Function

-

Now we include both the event and forcing function described in the previous sections. The deSolve solution is provided below.

-
# forcing function
-force <- approxfun(c(0,3,5,7,9), c(0,3,0,3,0), method = "constant", rule = 2)
-# event table
-dose <- data.frame(var = "Intestine",
-                   time = seq(2,9,by=1),
-                   value = 2,
-                   method = "add")
-# ODE
-two_comp_fun <- function(t, y, p) {
-  list2env(p, environment())
-  dy1 <- -a*y[1] + force(t)
-  dy2 <- a*y[1] - b*y[2]
-  return(list(c(dy1,dy2)))
-}
-
-sim_desolve <- ode(func = two_comp_fun, times = times_desolve, y = yini, parms = pars, method = "bdf", events = list(data = dose))
-plot(sim_desolve)
-

-

The comparable solution using Stan is provided below. For details see the Stan program models/two_comp_ef.stan.

-
# break up time into event segments
-seg <- dose$time
-seg_indxs <- which(times %in% seg)
-
-if (tail(seg_indxs,1) != length(times)) {
-  seg_indxs <- c(seg_indxs, length(times))
-}
-
-if (seg_indxs[1] != 1) {
-  seg_indxs <- c(1, seg_indxs)
-}
-
-# for readability create a Stan event table
-seg_start <- head(seg_indxs, -1) + 1
-seg_start[1] <- 1
-seg_end <- tail(seg_indxs, -1) 
-seg_value <- cbind(dose$value, 0)
-seg_value <- rbind(yini, seg_value)
-
-# create stan data
-stan_data <- list(theta = unlist(pars),
-                  ts = times,
-                  T = length(times),
-                  E = nrow(seg_value),
-                  S = ncol(seg_value),
-                  K = length(unlist(pars)),
-                  y0 = seg_value,
-                  t0 = times[seg_start] - (times[2]-times[1]),
-                  start = seg_start,
-                  end = seg_end)
-
-# simulate from ode
-tc_event_force_stan <- stan("models/two_comp_ef.stan", data = stan_data, algorithm = "Fixed_param", iter = 1, chains = 1)
-
## 
-## SAMPLING FOR MODEL 'two_comp_ef' NOW (CHAIN 1).
-## Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
-## Chain 1:  Elapsed Time: 0 seconds (Warm-up)
-## Chain 1:                0.009013 seconds (Sampling)
-## Chain 1:                0.009013 seconds (Total)
-

Below we overlay the solution using Stan on top of the solution defined by deSolve.

-
sim_tc_event_force_stan <- as.matrix(tc_event_force_stan)
-sim_tc_event_force_stan <- cbind(time = times,
-                           y1 = unname(sim_tc_event_force_stan[,grep("^y_hat\\[[[:digit:]]+,1\\]$", colnames(sim_tc_event_force_stan))]),
-                           y2 = unname(sim_tc_event_force_stan[,grep("^y_hat\\[[[:digit:]]+,2\\]$", colnames(sim_tc_event_force_stan))]))
-
-par(mfrow = c(1,2))
-plot(sim_desolve[,"time"], sim_desolve[,"Intestine"], type = "l", lwd = 4,
-     main = "Intestine",
-     xlab = "time", ylab = "")
-lines(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y1"], col = "#CCFF00")
-legend("topleft", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-
-plot(sim_desolve[,"time"], sim_desolve[,"Blood"], type = "l", lwd = 4,
-     main = "Blood",
-     xlab = "time", ylab = "")
-lines(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y2"], col = "#CCFF00")
-legend("bottomright", c("deSolve","Stan"), col = c("#000000","#CCFF00"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-

-
-
-

Estimation

-

All the examples above involve simulation and assume that we know the value of the parameters. Alternatively, we can use Stan to infer the parameters used to generate the numerical solution (once we have observed the solution). This means that we can treat the observed solution as our outcome variable and infer the parameters \(a\) and \(b\) that generated the solution.

-

As an example we use the ODE defined in the previous section along with its numerical solution computed in Stan. To resemble real data with measurment error we add some arbitrary noise to this solution (this could be done in the simulation itself within the generated quantities {} block). Below we overlay the true solution along with the solution with noise.

-
# add noise to the solution
-y <- sim_tc_event_force_stan
-y[,2] <- y[,2] + rnorm(length(times), 0, 0.1)
-y[,3] <- y[,3] + rnorm(length(times), 0, 0.1)
-
-y[,2] <- ifelse(y[,2] < 0, 0, y[,2])
-y[,3] <- ifelse(y[,3] < 0, 0, y[,3])
-
-par(mfrow = c(1,2))
-plot(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y1"], type = "l", lwd = 4,
-     main = "Intestine",
-     xlab = "time", ylab = "")
-lines(y[,"time"], y[,"y1"], col = "#FF6688")
-legend("topleft", c("Stan","Stan with noise"), col = c("#000000","#FF6688"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-
-plot(sim_tc_event_force_stan[,"time"], sim_tc_event_force_stan[,"y2"], type = "l", lwd = 4,
-     main = "Blood",
-     xlab = "time", ylab = "")
-lines(y[,"time"], y[,"y2"], col = "#FF6688")
-legend("bottomright", c("Stan","Stan with noise"), col = c("#000000","#FF6688"), pch = c(15,15), pt.cex = 2, cex = 0.8)
-

-

To estimate the parameters \(a\) and \(b\) we move the integrator specification from the generated quantities {} block to the model {} block. In addition to the integrator we incorporate the following pieces into the model {} block of the Stan program:

-

Likelihood \[ -y \sim \mathcal{N}(\hat{y}, \sigma) -\]

-

Priors \[ -\begin{align*} -\sigma &\sim \mathcal{N}(0, 0.1) \\ -a &\sim \mathcal{N}(5, 1) \\ -b &\sim \mathcal{N}(0, 1) \\ -\end{align*} -\]

-

Where \(y\) is the data we acquired from simulation in the previous section (with noise added) and \(\hat{y}\) is the simulation computed by the ODE integrator in the Stan program. Essentially we are modeling \(y\) according to the deterministic process that generates \(\hat{y}\) along with some measurement error. For details on how this is implemented see the Stan program models/two_comp_ef_fit.stan. Below we fit the model to the data.

-
stan_data$y <- y[,-1]
-tc_event_force_fit <- stan("models/two_comp_ef_fit.stan", data = stan_data, chains = 4, iter = 500, cores = 4)
-

Looking at the parameter values and diagnostics we find that we appropriately recovered the parameters that generated the data.

-
print(tc_event_force_fit)
-
## Inference for Stan model: two_comp_ef_fit.
-## 4 chains, each with iter=500; warmup=250; thin=1; 
-## post-warmup draws per chain=250, total post-warmup draws=1000.
-## 
-##            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
-## theta[1]   6.07    0.00 0.08   5.91   6.02   6.07   6.13   6.21   405 1.00
-## theta[2]   0.60    0.00 0.00   0.60   0.60   0.60   0.60   0.60  1044 1.00
-## sigma      0.08    0.00 0.00   0.08   0.08   0.08   0.09   0.09   264 1.01
-## lp__     186.79    0.08 1.30 183.42 186.19 187.14 187.74 188.31   270 1.01
-## 
-## Samples were drawn using NUTS(diag_e) at Tue Oct 15 16:11:05 2019.
-## For each parameter, n_eff is a crude measure of effective sample size,
-## and Rhat is the potential scale reduction factor on split chains (at 
-## convergence, Rhat=1).
-
mcmc_trace(as.array(tc_event_force_fit), regex_pars = "theta|sigma")
-

-

Using these results the user can continue with useful Bayesian post-estimation steps, such as generating posterior predictions and constructing credible intervals.

-
-
-

References

-

Gabry, J. & Mahr, T. (2019) bayesplot: Plotting for Bayesian Models. R package version 1.7.0, mc-stan.org/bayesplot.

-

Mattuck, A., Miller, H., Orloff, J., & Lewis, J. (2011) Differential Equations. Massachusetts Institute of Technology: MIT OpenCourseWare. https://ocw.mit.edu.

-

Soetaert, K., Cash, J., & Mazzia F. (2012) Solving Differential Equations in R. Use R!

-

Soetaert, K., Petzoldt, T., & Setzer, R.W. (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software. 33(9), 1-25. http://www.jstatsoft.org/v33/i09/.

-

Stan Development Team (2019) Stan Modeling Language User’s Guide. Version 2.20. http://mc-stan.org.

-

Stan Development Team (2019) RStan: the R interface to Stan. R package version 2.19.1, http://mc-stan.org/.

-
-
-

Computing Environment

-
devtools::session_info("rstan")
-
## Session info -------------------------------------------------------------
-
##  setting  value                       
-##  version  R version 3.4.0 (2017-04-21)
-##  system   x86_64, darwin15.6.0        
-##  ui       X11                         
-##  language (EN)                        
-##  collate  en_US.UTF-8                 
-##  tz       <NA>                        
-##  date     2019-10-15
-
## Packages -----------------------------------------------------------------
-
##  package      * version    date       source                         
-##  assertthat     0.2.0      2017-04-11 CRAN (R 3.4.0)                 
-##  backports      1.1.2      2017-12-13 cran (@1.1.2)                  
-##  base64enc      0.1-3      2015-07-28 CRAN (R 3.4.0)                 
-##  BH             1.66.0-1   2018-02-13 cran (@1.66.0-)                
-##  callr          3.0.0.9000 2018-09-06 Github (r-pkgs/callr@10c1269)  
-##  cli            1.0.1      2018-09-25 cran (@1.0.1)                  
-##  colorspace     1.3-2      2016-12-14 CRAN (R 3.4.0)                 
-##  compiler       3.4.0      2017-04-21 local                          
-##  crayon         1.3.4      2017-09-16 cran (@1.3.4)                  
-##  desc           1.2.0      2018-05-01 cran (@1.2.0)                  
-##  digest         0.6.18     2018-10-10 cran (@0.6.18)                 
-##  ggplot2      * 3.0.0      2018-07-03 cran (@3.0.0)                  
-##  glue           1.3.0      2018-07-17 cran (@1.3.0)                  
-##  graphics     * 3.4.0      2017-04-21 local                          
-##  grDevices    * 3.4.0      2017-04-21 local                          
-##  grid           3.4.0      2017-04-21 local                          
-##  gridExtra      2.3        2017-09-09 cran (@2.3)                    
-##  gtable         0.2.0      2016-02-26 CRAN (R 3.4.0)                 
-##  inline         0.3.15     2018-05-18 cran (@0.3.15)                 
-##  labeling       0.3        2014-08-23 CRAN (R 3.4.0)                 
-##  lattice        0.20-35    2017-03-25 CRAN (R 3.4.0)                 
-##  lazyeval       0.2.1      2017-10-29 cran (@0.2.1)                  
-##  loo            2.0.0      2018-04-11 CRAN (R 3.4.4)                 
-##  magrittr       1.5        2014-11-22 CRAN (R 3.4.0)                 
-##  MASS           7.3-47     2017-02-26 CRAN (R 3.4.0)                 
-##  Matrix         1.2-14     2018-04-09 CRAN (R 3.4.4)                 
-##  matrixStats    0.52.2     2017-04-14 CRAN (R 3.4.0)                 
-##  methods      * 3.4.0      2017-04-21 local                          
-##  mgcv           1.8-17     2017-02-08 CRAN (R 3.4.0)                 
-##  munsell        0.5.0      2018-06-12 cran (@0.5.0)                  
-##  nlme           3.1-131    2017-02-06 CRAN (R 3.4.0)                 
-##  parallel       3.4.0      2017-04-21 local                          
-##  pillar         1.1.0      2018-01-14 cran (@1.1.0)                  
-##  pkgbuild       1.0.0      2018-09-06 Github (r-lib/pkgbuild@b748041)
-##  plyr           1.8.4      2016-06-08 CRAN (R 3.4.0)                 
-##  processx       3.2.0      2018-08-16 cran (@3.2.0)                  
-##  ps             1.1.0      2018-08-10 cran (@1.1.0)                  
-##  R6             2.3.0      2018-10-04 cran (@2.3.0)                  
-##  RColorBrewer   1.1-2      2014-12-07 CRAN (R 3.4.0)                 
-##  Rcpp           1.0.0      2018-11-07 cran (@1.0.0)                  
-##  RcppEigen      0.3.3.4.0  2018-02-07 cran (@0.3.3.4)                
-##  reshape2       1.4.3      2017-12-11 cran (@1.4.3)                  
-##  rlang          0.3.1      2019-01-08 CRAN (R 3.4.4)                 
-##  rprojroot      1.3-2      2018-01-03 cran (@1.3-2)                  
-##  rstan        * 2.18.1     2018-09-07 local                          
-##  scales         1.0.0      2018-08-09 cran (@1.0.0)                  
-##  StanHeaders  * 2.18.0     2018-09-07 local                          
-##  stats        * 3.4.0      2017-04-21 local                          
-##  stats4         3.4.0      2017-04-21 local                          
-##  stringi        1.2.4      2018-07-20 cran (@1.2.4)                  
-##  stringr        1.3.1      2018-05-10 cran (@1.3.1)                  
-##  tibble         1.4.2      2018-01-22 cran (@1.4.2)                  
-##  tools          3.4.0      2017-04-21 local                          
-##  utf8           1.1.3      2018-01-03 cran (@1.1.3)                  
-##  utils        * 3.4.0      2017-04-21 local                          
-##  viridisLite    0.3.0      2018-02-01 cran (@0.3.0)                  
-##  withr          2.1.2      2018-03-15 CRAN (R 3.4.4)
-
-
-

Licences

-

Code (c) 2019 Imad Ali, BSD 3.

-

Text (c) 2019 Imad Ali, licensed under CC-BY-NC 4.0.

-
- - - - -
- - - - - - - - - - - - - - -