Merge pull request #174 from stan-dev/case-study/liza-sir-model

corrections

Author: bob-carpenter

knitr/disease_transmission/biblio.bib

@@ -0,0 +1,199 @@
+@misc{chatzilena2019contemporary,
+    title={Contemporary statistical inference for infectious disease models using Stan},
+    author={Anastasia Chatzilena and Edwin van Leeuwen and Oliver Ratmann and Marc Baguelin and Nikolaos Demiris},
+    year={2019},
+    eprint={1903.00423},
+    archivePrefix={arXiv},
+    primaryClass={stat.AP}
+}
+
+@article{kermack_contributions_1991,
+	title = {Contributions to the mathematical theory of epidemics{I}},
+	volume = {53},
+	issn = {1522-9602},
+	url = {https://doi.org/10.1007/BF02464423},
+	doi = {10.1007/BF02464423},
+	number = {1},
+	journal = {Bulletin of Mathematical Biology},
+	author = {Kermack, W. O. and McKendrick, A. G.},
+	month = mar,
+	year = {1991},
+	pages = {33--55}
+}
+
+@misc{carpenter2015stan,
+    title={The Stan Math Library: Reverse-Mode Automatic Differentiation in C++},
+    author={Bob Carpenter and Matthew D. Hoffman and Marcus Brubaker and Daniel Lee and Peter Li and Michael Betancourt},
+    year={2015},
+    eprint={1509.07164},
+    archivePrefix={arXiv},
+    primaryClass={cs.MS}
+}
+
+@article {PMID:12995455,
+	Title = {The analysis of equilibrium in malaria},
+	Author = {MACDONALD, G},
+	Number = {9},
+	Volume = {49},
+	Month = {September},
+	Year = {1952},
+	Journal = {Tropical diseases bulletin},
+	ISSN = {0041-3240},
+	Pages = {813—829},
+	URL = {http://europepmc.org/abstract/MED/12995455}
+}
+
+@conference{margossian_differential_2017,
+	Author = {Margossian, Charles C and Gillespie, William R},
+	Booktitle = {StanCon 2017},
+	Journal = {Journal of Pharmacokinetics and Pharmacodynamics},
+	Month = {January},
+	Title = {Differential Equation Based Models in Stan},
+	Year = {2017},
+	doi = {10.5281/zenodo.1284264}}
+
+@misc{carpenter_predator-prey_2018,
+	title = {Predator-{Prey} {Population} {Dynamics}: the {Lotka}-{Volterra} model in {Stan}},
+	url = {https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html},
+	author = {Carpenter, Bob},
+	year = {2018}
+}
+@misc{prior_choice,
+	title = {Prior Choice Recommendations},
+	url = {https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations},
+	author = {Stan Team},
+	year = {2020}
+}
+
+@article{hall_infectious_1992,
+	title = {Infectious {Diseases} of {Humans}. {R}. {M}. {Anderson} \& {R}. {M}. {May}. {Oxford} etc.: {Oxford} {University} {Press}, 1991. viii+757 pp. {Price} £50. {ISBN} 0-19-854599-1},
+	volume = {86},
+	issn = {0035-9203},
+	url = {https://doi.org/10.1016/0035-9203(92)90276-I},
+	doi = {10.1016/0035-9203(92)90276-I},
+	number = {4},
+	journal = {Transactions of The Royal Society of Tropical Medicine and Hygiene},
+	author = {Hall, A.J.},
+	year = {1992},
+	pages = {461--461}
+}
+
+@article{Riou_covid_2020,
+              Author = {Riou, Julien and Hauser, Anthony and Counotte, Michel J and Margossian, Charles C
+                              and Konstantinoudis, Garyfallos and Low, Nicola and Althaus, Christian L},
+              Title = {Estimation of SARS-CoV-2 mortality during the early stages of an epidemic: 
+                          a modeling study in Hubei, China and northern Italy},
+               journal = {medRxiv:2020.03.04.20031104},
+               year = {2020}}
+
+@article{Betancourt_hmc_2018,
+	Author = {Betancourt, Michael},
+	Journal = {arXiv:1701.02434v1},
+	Month = {January},
+	Title = {A Conceptual Introduction to Hamiltonian Monte Carlo},
+	Year = {2018}}
+	
+@article{Hoffman_nuts_2014,
+	Author = {Hoffman, Matthew D. and Gelman, Andrew},
+	Date-Added = {2016-10-07 17:06:45 +0000},
+	Date-Modified = {2017-08-17 09:27:44 +0000},
+	Journal = {Journal of Machine Learning Research},
+	Month = {April},
+	Pages = {1593-1623},
+	Title = {The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo},
+	Year = {2014}}
+
+@techreport{Hindmarsh_ode_2020,
+	Author = {Hindmarsh, Alan and Serban, Radu},
+	Date-Modified = {2020-01-16 22:47:24 -0500},
+	Institution = {Lawrence Livermore National Laboratory},
+	Month = {January},
+	Title = {User Documentation for CVODES v5.1.0},
+	Year = {2020}}
+
+@book{Griewank_ad_2008,
+	Author = {Griewank, Andreas and Walther, Andrea},
+	Edition = {Second},
+	Publisher = {Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA},
+	Title = {Evaluating derivatives},
+	Year = {2008}}
+
+@article{Margossian_ad_2019,
+	Author = {Charles C. Margossian},
+	Month = {3},
+	Title = {A Review of automatic differentiation and its efficient implementation},
+	journal = {Wiley interdisciplinary reviews: data mining and knowledge discovery},
+	volume = {9},
+        issue = {4},
+        doi = {10.1002/WIDM.1305},
+	Year = {2019}}
+
+@article{margossian_mixed_solver_2017,
+	Author = {Margossian, Charles C and Gillespie, William R},
+	Booktitle = {American },
+	Journal = {Journal of Pharmacokinetics and Pharmacodynamics},
+	Month = {October},
+	Title = {Gaining Efficiency by Combining Analytical and Numerical Methods to Solve ODEs: Implementation in Stan and Application to Bayesian PK/PD },
+	Volume = {44},
+	Year = {2017}}
+
+@article{Flaxman_covid_2020,
+        Author = {Seth Flaxman and Swapnil Mishra and Axel Gandy and  H Juliette T Unwin and Helen Coupland and ... and Samir Bhatt},
+        title = {Report 13 - Estimating the number of infections and the impact of non-pharmaceutical interventions on COVID-19 in 11 European countries},
+        journal = {arXiv:2004.11342},
+        year = {2020},
+        month = {April}}
+
+@article{Carpenter_stan_2017,
+         Title = {Stan: A Probabilistic Programming Language},
+         Author = {Carpenter, Bob and Gelman, Andrew and Hoffman, Matt and Lee, Daniel and Goodrich, Ben and Betancourt, Michael and Brubaker, Marcus A. and Guo, Jiqiang and Li, Peter and Riddel, Allen},
+         Journal = {Journal of Statistical Software},
+         volume = {76},
+         issue = {1},
+         pages = {1 --32},
+         Year = {2017},
+         doi = {10.18637/jss.v076.i01}}
+
+@article{Chatzilena_tutorial_2019,
+  author = {Chatzilena, Anastasia and van Leeuwen, Edwin and Ratmann, Olivier
+    and Baguelin, Olivier and Demiris, Nikolaos},
+  title =  {Contemporary statistical inference for infectious disease models using Stan},
+  journal = {Epidemics},
+  volume = {29},
+  year = {2019},
+  doi = {https://doi.org/10.1016/j.epidem.2019.100367}}
+  
+@article{Mihaljevic_tutorial_2016,
+  title = {Estimating transmission by fitting mechanistic models in Stan},
+  author = {Mihaljevic, Joseph},
+  year = {2016},
+  url = {https://jrmihalj.github.io/estimating-transmission-by-fitting-mechanistic-models-in-Stan/}
+  }
+  
+@article{Weber_ode_2018,
+  title = {Solving ODEs in the wild: Scalable pharmacometrics with Stan},
+  author = {Weber, Sebastian},
+  journal = {StanCon Helsinki 2018},
+  year = {2018}}
+
+@article{Carpenter_predatory-prey_2018,
+  title = {Predator-Prey Population Dynamics: the Lotka-Volterra model in Stan},
+  author = {Carpenter, Bob},
+  url = {https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html},
+  year = {2018}}
+
+@article{Talts_sbc_2018,
+  author = {Talts, Sean and Betancourt, Michael and Simpson, Daniel and Vehtari, Aki and Gelman, Andrew},
+  title = {Validating Bayesian inference algorithms with simulation-based calibration},
+  journal = {arXiv:1804.06788v1},
+  year = {2018}}
+
+@article{Roberts_mcmc_2004,
+  author = {Roberts, Gareth O and Rosenthal, Jeffrey S},
+  title = {General state space Markov chains and MCMC algorithms},
+  journal = {Probability survey},
+  volume = {1},
+  pages = {20 - 71},
+  doi = {10.1214/154957804100000024},
+  year = {2004}}
+    

knitr/disease_transmission/boarding_school_case_study.Rmd

@@ -305,21 +305,25 @@ We motivate this signature in our discussion on scaling ODEs (section 3).
 
 In our example, the ODEs for the SIR model is defined as follows:
 ```{stan, eval=F, output.var="md"}
-  //theta[1] = beta
-  //theta[2] = gamma
-  //x_i[1] = N, total population
-  // y = S, I, R
-  real[] sir(real t,
-             real[] y,
-             real[] theta,
-             real[] x_r,
-             int[] x_i) {
-    real dydt[3];
-    dydt[1] = - theta[1] * y[1] * y[2] / x_i[1];
-    dydt[2] = theta[1] * y[1] * y[2] / x_i[1] - theta[2] * y[2];
-    dydt[3] = theta[2] * y[2];
-    return dydt;
+functions {
+  real[] sir(real t, real[] y, real[] theta, 
+             real[] x_r, int[] x_i) {
+
+      real S = y[1];
+      real I = y[2];
+      real R = y[3];
+      real N = x_i[1];
+      
+      real beta = theta[1];
+      real gamma = theta[2];
+      
+      real dS_dt = -beta * I * S / N;
+      real dI_dt =  beta * I * S / N - gamma * I;
+      real dR_dt =  gamma * I;
+      
+      return {dS_dt, dI_dt, dR_dt};
   }
+}
 ```
 
 We evaluate the solution numerically by using one of Stan's numerical integrators.
@@ -343,6 +347,13 @@ with
 
 We now have all the ingredients to solve our ODE.
 
+Note that in the given example, when we assume that the total population remains constant, the three derivatives $\frac{dS}{dt}$, $\frac{dI}{dt}$, $\frac{dR}{dt}$ sum up to $0$:  We can use this fact to improve computational efficiency of the `sir` function by deriving the value of $\frac{dI}{dt}$ from $\frac{dS}{dt}$ and $\frac{dR}{dt}$:
+```{stan, eval=F, output.var="md"}
+      real dS_dt = -beta * I * S / N;
+      real dR_dt =  gamma * I;
+      real dI_dt =  -(dS_dt + dR_dt);
+```
+
 ### Building the model in Stan
 
 We next code the model in Stan, working through the various coding blocks.
@@ -553,24 +564,24 @@ ggplot(smr_y, mapping = aes(x = t)) +
 ## Complete Stan program
 ```{stan, eval=F, output.var="md"}
 functions {
-  //theta[1] = beta
-  //theta[2] = gamma
-  //x_i[1] = N, total population
-  // y = S, I, R
-  real[] sir(real t,
-             real[] y,
-             real[] theta,
-             real[] x_r,
-             int[] x_i) {
-    real dydt[3];
-    dydt[1] = - theta[1] * y[1] * y[2] / x_i[1];
-    dydt[2] = theta[1] * y[1] * y[2] / x_i[1] - theta[2] * y[2];
-    dydt[3] = theta[2] * y[2];
-    return dydt;
-
+  real[] sir(real t, real[] y, real[] theta, 
+             real[] x_r, int[] x_i) {
+
+      real S = y[1];
+      real I = y[2];
+      real R = y[3];
+      real N = x_i[1];
+      
+      real beta = theta[1];
+      real gamma = theta[2];
+      
+      real dS_dt = -beta * I * S / N;
+      real dI_dt =  beta * I * S / N - gamma * I;
+      real dR_dt =  gamma * I;
+      
+      return {dS_dt, dI_dt, dR_dt};
   }
 }
-
 data {
   int<lower=1> n_days;
   real y0[3];
@@ -579,18 +590,15 @@ data {
   int N;
   int cases[n_days];
 }
-
 transformed data {
   real x_r[0];
   int x_i[1] = { N };
 }
-
 parameters {
   real<lower=0> gamma;
   real<lower=0> beta;
   real<lower=0> phi_inv;
 }
-
 transformed parameters{
   real y[n_days, 3];
   real phi = 1. / phi_inv;
@@ -602,8 +610,6 @@ transformed parameters{
     y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
   }
 }
-
-
 model {
   //priors
   beta ~ normal(2, 1);
@@ -614,14 +620,12 @@ model {
   //col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people 
   cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
 }
-
 generated quantities {
   real R0 = beta / gamma;
   real recovery_time = 1 / gamma;
   real pred_cases[n_days];
   pred_cases = neg_binomial_2_rng(col(to_matrix(y), 2), phi);
 }
-
 ```
 
 # 2 Using simulated data to understand our model
@@ -716,7 +720,7 @@ colnames(smr_pred) <- make.names(colnames(smr_pred)) # to remove % in the col na
 ggplot(smr_pred, mapping=aes(x=t)) +
   geom_ribbon(aes(ymin = X5., ymax = X95.), fill = "orange", alpha = 0.6) +
   geom_line(mapping=aes(x=t, y=X50.)) + 
-  geom_abline(intercept=577, color="red", ) +
+  geom_abline(intercept=577, color="red" ) +
   geom_text(x=1.8, y=560, label="Population size", color="red") +
   labs(x = "Day", y="Number of students in bed")
 ```