v0.2.2

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: odemodeling
 Title: Building, fitting, and validating Bayesian ODE models in 'Stan'
-Version: 0.2.1
+Version: 0.2.2
 Authors@R: 
     person(given = "Juho",
            family = "Timonen",
@@ -15,7 +15,7 @@ License: BSD_3_clause + file LICENCE
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.3
 SystemRequirements: CmdStan (https://mc-stan.org/users/interfaces/cmdstan)
 Imports: 
     methods,

NAMESPACE

@@ -9,6 +9,7 @@ export(bdf_list)
 export(ckrk)
 export(compute_reliability_metrics)
 export(euler)
+export(example_ode_model)
 export(log_ratios)
 export(max_abs_loglik_diff)
 export(max_abs_odesol_diff)

NEWS.md

@@ -1,5 +1,9 @@
 # Version 0.2
 
+## odemodeling 0.2.2
+  * Improves documentation.
+  * Exports `example_ode_model()`.
+
 ## odemodeling 0.2.0
   * Adds the `$diagnose()` method for the `OdeModel` class.
   * Stores the results of `diagnose` and `stansummary`in the

R/classes-fit.R

@@ -85,7 +85,6 @@ OdeModelMCMC <- R6::R6Class("OdeModelMCMC",
                    solver = NULL,
                    fitted_params = NULL,
                    ...) {
-
       # Handle input
       t0 <- replace_if_null(t0, self$t0)
       t <- replace_if_null(t, self$t)

R/classes-model.R

@@ -126,7 +126,6 @@ OdeModel <- R6::R6Class("OdeModel", list(
                  solver = rk45(),
                  params = NULL,
                  ...) {
-
     # Full Stan data
     model <- self
     sd <- create_standata(model, t0, t, solver)
@@ -166,7 +165,6 @@ OdeModel <- R6::R6Class("OdeModel", list(
                     data = list(),
                     solver = rk45(),
                     ...) {
-
     # Check and handle input
     sd <- create_standata(self, t0, t, solver)
     full_data <- c(sd, data)
@@ -212,7 +210,6 @@ OdeModel <- R6::R6Class("OdeModel", list(
                       error = Inf,
                       epsilon = 1e-6,
                       ...) {
-
     # Check and handle input
     sd <- create_standata(self, t0, t, solver)
     full_data <- c(sd, data)
@@ -273,7 +270,7 @@ OdeModel <- R6::R6Class("OdeModel", list(
       fit <- model$sample(
         t0 = t0,
         t = t,
-        data  = data,
+        data = data,
         solver = solver,
         chains = chains,
         ...

R/classes-standecl.R

@@ -1,4 +1,3 @@
-
 #' A Stan variable declaration (an abstract base class)
 #'
 #' @field name name of the variable

R/main-examples.R

@@ -9,6 +9,7 @@
 #'   \item `"lv"` - Lotka-Volterra model
 #' }
 #' @return An object of class `OdeModel`.
+#' @export
 #' @family model constructor functions
 example_ode_model <- function(name, prior_only = FALSE, ...) {
   choices <- c("gsir", "tmdd", "lv")
@@ -29,7 +30,6 @@ example_ode_model <- function(name, prior_only = FALSE, ...) {
 
 # Group-stratified SIR example model
 example_ode_model_gsir <- function(prior_only, ...) {
-
   # Time points
   N <- stan_dim("N", lower = 0) # number of timepoints
 
@@ -134,7 +134,6 @@ example_ode_model_gsir <- function(prior_only, ...) {
 
 # TMDD example model
 example_ode_model_tmdd <- function(prior_only, ...) {
-
   # Dimensions and other data
   N <- stan_dim("N", lower = 1) # number of time points
   D <- stan_dim("D", lower = 1) # ODE system dimension
@@ -205,7 +204,6 @@ example_ode_model_tmdd <- function(prior_only, ...) {
 
 # Lotka-Volterra example model
 example_ode_model_lv <- function(prior_only = FALSE, ...) {
-
   # Dimensions and other data
   N <- stan_dim("N", lower = 1) # number of time points
   D <- stan_dim("D", lower = 1) # ODE system dimension

R/main-model.R

@@ -52,7 +52,6 @@ ode_model <- function(N,
                       verbose = FALSE,
                       compile = TRUE,
                       sig_figs = 18) {
-
   # Argument checks
   choices_vars <- c("StanParameter", "StanTransformation", "StanDeclaration")
   checkmate::assert_class(N, "StanDimension")

R/odemodeling-package.R

@@ -45,6 +45,33 @@
 #'   \item See \code{\link{compare_odefits}} for functions to compare
 #'   different ODE model simulations and fits.
 #' }
+#' @section Importance sampling for reliable and efficient inference in Bayesian ODE models:
+#' Our proposed workflow is to
+#' \enumerate{
+#' \item Select an initial ODE solver M.
+#' \item Sample the parameter posterior using MCMC with M as the ODE solver.
+#' \item  Compute certain metrics using a more accurate solver M∗.
+#' \item Increase the accuracy of M∗ and repeat Step 3 until the metrics converge.
+#'  If the Pareto-k metric converges to a value larger than 0.7, increase the accuracy
+#' of M and go back to Step 2.
+#' \item Compute any posterior estimates using final importance weights. See
+#' Timonen et al. (2022) below.
+#' }
+#' The algorithm can be used to validate the reliability, and correct the errors
+#' of a given method M, which can be for example a software default. On the other hand,
+#' a smart initial selection of M can provide speed gains
+#' compared to often rather conservatively set software defaults, while still maintaining
+#' reliability of the inferences.
+#'
+#' We generally recommend selecting M initially so that sampling is as fast as possible.
+#' For example, for non-adaptive ODE solvers, one can first try using the smallest
+#' sensible number of steps that does not result in immediate failure.
+#' Selecting a good M is more difficult in the case of adaptive solvers.
+#' We have
+#' observed that tolerances on the order of 1e-4 to 1e-3 generally work well.
+#'
+#' See the \code{reliability()} method of the \code{\link{OdeModelMCMC}} class.
+#'
 #' @section Tutorial:
 #' See the tutorial vignette.
 #' @references

R/utils-create_input.R

@@ -1,4 +1,3 @@
-
 # Create full Stan data for sampling or gq
 create_standata <- function(model, t0, t, solver) {
   checkmate::assert_class(model, "OdeModel")