diff --git a/.DS_Store b/.DS_Store index b4c28ec..7099bf3 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/.Rprofile b/.Rprofile new file mode 100644 index 0000000..dea09f7 --- /dev/null +++ b/.Rprofile @@ -0,0 +1 @@ +options(pkg.build_extra_flags = FALSE) diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index 16a5d07..235fbf6 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 1.1.0 -Date: 2025-05-16 08:32:30 UTC -SHA: 5bffcb29721b4c7f8d2a876aadb3bbf99b197511 +Version: 1.2.1 +Date: 2025-06-02 07:51:53 UTC +SHA: 4f2bc2b974fe84b46895036a04980c02584faee9 diff --git a/DESCRIPTION b/DESCRIPTION index 61b05b2..5dc0af2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: denim Type: Package Title: Generate and Simulate Deterministic Discrete-Time Compartmental Models -Version: 1.1.0 -Date: 2025-05-16 +Version: 1.2.1 +Date: 2025-06-02 Authors@R: c( person("Thinh", "Ong", , "thinhop@oucru.org", role = c("aut", "cph"), comment = c(ORCID = "0000-0001-6772-9291")), @@ -20,9 +20,12 @@ License: MIT + file LICENSE URL: https://drthinhong.com/denim/, https://github.com/thinhong/denim BugReports: https://github.com/thinhong/denim/issues +Depends: R (>= 4.1.0) Imports: Rcpp (>= 1.0.6), - viridisLite + colorspace, + rlang, + glue Suggests: covr, knitr, @@ -31,6 +34,7 @@ Suggests: waldo, xml2, deSolve, + tidyverse, DiagrammeR LinkingTo: Rcpp, diff --git a/NAMESPACE b/NAMESPACE index 312d3e3..da94b16 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,8 +6,12 @@ export(d_exponential) export(d_gamma) export(d_lognormal) export(d_weibull) +export(denim_dsl) export(mathexpr) export(nonparametric) export(sim) import(Rcpp) +import(colorspace) +import(glue) +import(rlang) useDynLib(denim, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index 6b25014..2e873ce 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# denim 1.2.1 +* Minor runtime improvement +* Add option to select color palette for plot + +# denim 1.2.0 +* Add denim DSL support +* Add error handlers + # denim 1.1.0 * Fix multinomial vs competing risks transition * Add distribute initial value feature diff --git a/R/distribution.R b/R/distribution.R index 276bb29..dafea8e 100644 --- a/R/distribution.R +++ b/R/distribution.R @@ -1,4 +1,15 @@ # Constructors for distributions +evaluate_par <- function(par){ + + if(is.symbol(par) | is.character(par)) { + par <- as.character(par) + } else{ + par <- as.numeric(eval(par)) + } + + par +} + #' Discrete gamma distribution #' @@ -9,8 +20,19 @@ #' #' @examples #' transitions <- list("S -> I" = d_gamma(rate = 1, shape = 5)) +#' transitions_dsl <- denim_dsl({S -> I = d_gamma(rate = 1, shape = 5)}) +#' # define model parameters as distributional parameters +#' transitions_dsl <- denim_dsl({S -> I = d_gamma(rate = i_rate, shape = i_shape)}) #' @export d_gamma <- function(rate, shape, dist_init = FALSE) { + # capture expression for evaluation + rate <- substitute(rate) + shape <- substitute(shape) + + # check whether it is a fixed value (numeric) or a model parameter (string or expression) + rate <- evaluate_par(rate) + shape <- evaluate_par(shape) + distr <- list( distribution = "gamma", rate = rate, @@ -31,8 +53,17 @@ d_gamma <- function(rate, shape, dist_init = FALSE) { #' #' @examples #' transitions <- list("I -> D" = d_weibull(0.6, 2)) +#' transitions <- denim_dsl({ I -> D = d_weibull(0.6, 2) }) #' @export d_weibull <- function(scale, shape, dist_init = FALSE) { + # capture expression for evaluation + scale <- substitute(scale) + shape <- substitute(shape) + + # check whether it is a fixed value (numeric) or a model parameter (string or expression) + scale <- evaluate_par(scale) + shape <- evaluate_par(shape) + distr <- list( distribution = "weibull", scale = scale, @@ -52,10 +83,14 @@ d_weibull <- function(scale, shape, dist_init = FALSE) { #' #' @examples #' transitions <- list("I -> D" = d_exponential(0.3)) -#' -#' +#' transitions <- denim_dsl({I -> D = d_exponential(0.3)}) #' @export d_exponential <- function(rate, dist_init = FALSE) { + # capture expression for evaluation + rate <- substitute(rate) + # check whether it is a fixed value (numeric) or a model parameter (string or expression) + rate <- evaluate_par(rate) + distr <- list( distribution = "exponential", rate = rate, @@ -75,8 +110,17 @@ d_exponential <- function(rate, dist_init = FALSE) { #' #' @examples #' transitions <- list("I -> D" = d_lognormal(3, 0.6)) +#' transitions <- denim_dsl({I -> D = d_lognormal(3, 0.6)}) #' @export d_lognormal <- function(mu, sigma, dist_init = FALSE) { + # capture expression for evaluation + mu <- substitute(mu) + sigma <- substitute(sigma) + + # check whether it is a fixed value (numeric) or a model parameter (string or expression) + mu <- evaluate_par(mu) + sigma <- evaluate_par(sigma) + distr <- list( distribution = "lognormal", mu = mu, @@ -97,7 +141,8 @@ d_lognormal <- function(mu, sigma, dist_init = FALSE) { #' @return a Distribution object for simulator #' #' @examples -#' transitions <- list("S->I"=mathexpr("beta*S/N")) +#' transitions <- list("S->I"="beta*S/N") +#' transitions <- denim_dsl({S->I=beta*S/N}) #' # definition for parameters in the expression required #' params <- c(N = 1000, beta = 0.3) #' @export @@ -136,17 +181,49 @@ constant <- function(x) { #' #' Convert a vector of frequencies, percentages... into a distribution #' -#' @param ... a vector of values +#' @param x a vector of values #' @return a Distribution object for simulator #' @param dist_init whether to distribute initial value across subcompartments following this distribution. (default to FALSE, meaning init value is always in the first compartment)) #' #' @examples -#' transitions <- list("S->I"=nonparametric(0.1, 0.2, 0.5, 0.2)) +#' transitions <- list("S->I"=nonparametric( c(0.1, 0.2, 0.5, 0.2) )) +#' transitions <- denim_dsl({S->I=nonparametric( c(0.1, 0.2, 0.5, 0.2) )}) +#' # you can also define a model parameter for the distribution +#' transitions <- denim_dsl({S->I=nonparametric( dwelltime_dist )}) #' @export -nonparametric <- function(..., dist_init = FALSE) { +nonparametric <- function(x, dist_init = FALSE) { + # x <- substitute(c(...)) + # + # # remove "c" + # x <- as.character(x)[-1] + # + # # try parsing ... + # failed_parse <- FALSE + # try_eval <- sapply(parse(text = x), \(expr){ + # tryCatch( + # eval(expr), + # error = \(e) { + # failed_parse <<- TRUE + # } + # ) + # }) + # x <- if(!failed_parse) { + # try_eval + # } else{ + # x + # } + # + # # if ... is a parameter and more than 1 parameter is given -> throw error + # if(all(is.character(x)) && length(x) > 1){ + # stop("nonparametric only takes 1 parameter as input") + # } + # + + x <- substitute(x) + x <- evaluate_par(x) distr <- list( distribution = "nonparametric", - waitingTime = c(...), + waitingTime = x, dist_init = as.numeric(dist_init)) class(distr) <- c("Distribution", class(distr)) diff --git a/R/dsl.R b/R/dsl.R new file mode 100644 index 0000000..8ceeac0 --- /dev/null +++ b/R/dsl.R @@ -0,0 +1,75 @@ +#' Define transitions using denim's domain-specific language (DSL) +#' +#' This function parses model transitions defined in denim's DSL syntax +#' +#' @param x - an expression written in denim's DSL syntax. Each line should be a transition written in the +#' format `compartment -> out_compartment = expression` where expression can be either a math expression +#' or one of denim's built-in dwell time distribution function +#' +#' @return denim_transition object +#' @import rlang +#' @export +#' +#' @examples +#' transitions <- denim_dsl({ +#' S -> I = beta * (I/N) * S * timeStep +#' I -> R = d_gamma(rate = 1/4, shape = 3) +#' }) +denim_dsl <- function(x) { + # Capture the expression as a quoted expression + expr <- substitute(x) + # Convert the expression into a list of individual expressions + expr_list <- as.list(expr)[-1] + + transitions <- list() + + # helper function to flip the -> assigment which is non-standard in R + flip_assignment <- function(expr) { + # Split the string at "<-" + comps <- strsplit(expr, "<-")[[1]] + # Trim whitespace from both sides + comps <- trimws(comps) + + # Check that the expression is valid + if (length(comps) != 2) { + stop("Transition must be of the form 'comp_A -> comp_B' or 'weight*comp_A -> comp_B'") + } + + # Rearrange the parts and return + paste(comps[2], "->", comps[1]) + } + + for (expr in expr_list) { + # Each DSL argument should be an assignment (i.e. a call to `=`) + if (!is_call(expr, "=")) { + stop("Each argument must be an assignment using '='.") + } + + expr <- as.character(expr) + + # Extract lhs and rhs of the assignment + lhs_expr <- expr[2] + rhs_expr <- expr[3] + + # Convert the LHS to a string (e.g., "S -> I" or "36 * I -> R") + # lhs_str <- paste(deparse(lhs_expr), collapse = " ") + lhs_str <- flip_assignment(lhs_expr) + + # Check whether rhs is mathematical formula or call to d_* function + # If the rhs is a math formula, deparse it into a string. + # Otherwise (e.g., a call to d_gamma or d_exponential) leave it as an expression. + rhs_val <- if (grepl("^d_", rhs_expr) | grepl("^nonparametric", rhs_expr)) { + # If the call is to a function starting with "d_", leave it as an expression. + rhs_expr <- eval(parse_expr(rhs_expr)) + } else { + # Otherwise, deparse the expression to convert it to a string. + rhs_expr + } + + transitions[[lhs_str]] <- rhs_val + } + + class(transitions) <- "denim_transition" + + transitions +} diff --git a/R/model.R b/R/model.R index a62535e..75e91a0 100644 --- a/R/model.R +++ b/R/model.R @@ -1,4 +1,5 @@ # Constructor +#' @import glue newModel <- function(simulationDuration, errorTolerance, initialValues, parameters=NULL, transitions, timeStep = 1) { @@ -13,14 +14,40 @@ newModel <- function(simulationDuration, errorTolerance, initialValues, transitions[[i]] <- constant(transitions[[i]]) } - # if transition -> get parameters from Distribution object if( (class(transitions[[i]])=="Distribution")[[1]] ){ - if((transitions[[i]]$distribution) == "nonparametric" | - (transitions[[i]]$distribution) == "constant"| + if((transitions[[i]]$distribution) == "constant"| (transitions[[i]]$distribution) == "multinomial"){ - # dont do anything to non parametric model next + }else if((transitions[[i]]$distribution) == "nonparametric"){ + # check whether the distribution + # parameters are numeric or string, if it is a string (i.e. model input) + # make sure that value is provided as one of the parameters + if(is.character(transitions[[i]]$waitingTime)){ + par_val <- transitions[[i]]$waitingTime + + if(!(transitions[[i]]$waitingTime %in% names(parameters))){ + stop(glue::glue("Value for {par_val} of transition {names(transitions)[i]} must provided as one of the model parameters")) + } + transitions[[i]]$waitingTime <- parameters[[par_val]] + # remove waiting dist from parameters list + parameters[[par_val]] <- NULL + } }else if((transitions[[i]]$distribution) != "mathExpression"){ + # check whether the distribution + # parameters are numeric or string, if it is a string (i.e. model input) + # make sure that value is provided as one of the parameters + sapply(names(transitions[[i]][-1]), \(par_name){ + par_val <- transitions[[i]][[par_name]] + + if(is.character(par_val)){ + if(!(par_val %in% names(parameters))){ + stop(glue::glue("Value for {par_val} of transition {names(transitions)[i]} must provided as one of the model parameters")) + } + transitions[[i]][[par_name]] <<- parameters[[par_val]] + } + }) + + # also get parameters from Distribution objects parameters <- c(parameters, transitions[[i]][-1]) }else{ has_math_dist <- TRUE @@ -86,45 +113,3 @@ modelToJson <- function(mod) { } - -# experimental DSL -# TODO: test across various test case -# denim_transitions <- function(...) { -# exprs <- enexprs(...) -# transitions <- list() -# -# for (expr in exprs) { -# # Each DSL argument should be an assignment (i.e. a call to `=`) -# if (!is_call(expr, "=")) { -# stop("Each argument must be an assignment using '='.") -# } -# -# # Extract lhs and rhs of the assignment -# lhs_expr <- expr[[2]] -# rhs_expr <- expr[[3]] -# -# # Convert the LHS to a string (e.g., "S -> I" or "36 * I -> R") -# lhs_str <- paste(deparse(lhs_expr), collapse = " ") -# -# # Check whether rhs is mathematical formula or call to d_* function -# # If the rhs is a math formula, deparse it into a string. -# # Otherwise (e.g., a call to d_gamma or d_exponential) leave it as an expression. -# rhs_val <- if (is.call(rhs_expr)) { -# fn <- as.character(rhs_expr[[1]]) -# if (grepl("^d_", fn) | grepl("^nonparametric", fn)) { -# # If the call is to a function starting with "d_", leave it as an expression. -# rhs_expr -# } else { -# # Otherwise, deparse the expression to convert it to a string. -# paste(deparse(rhs_expr), collapse = " ") -# } -# } else { -# # For non-call objects (names, numbers), deparse them to a string. -# paste(deparse(rhs_expr), collapse = " ") -# } -# -# transitions[[lhs_str]] <- rhs_val -# } -# -# class(transitions) <- "denim_transition" -# } diff --git a/R/simulators.R b/R/simulators.R index 7f72560..27af19b 100644 --- a/R/simulators.R +++ b/R/simulators.R @@ -33,7 +33,7 @@ #' #' Simulation function that call the C++ simulator #' -#' @param transitions a list of transitions follows this format `"transition" = distribution()` +#' @param transitions output of function `denim_dsl()` or a list of transitions follows this format `"transition" = distribution()` #' @param initialValues a vector contains the initial values of all compartments defined #' in the **transitions**, follows this format `compartment_name = initial_value` #' @param parameters a vector contains values of any parameters that are not compartments, @@ -49,6 +49,13 @@ #' @export #' #' @examples +#' # model can be defined using denim DSL +#' transitions <- denim_dsl({ +#' S -> I = beta * S * I / N +#' I -> R = d_gamma(1/3, 2) +#' }) +#' +#' # or as a list #' transitions <- list( #' "S -> I" = "beta * S * I / N", #' "I -> R" = d_gamma(1/3, 2) @@ -94,16 +101,24 @@ sim <- function(transitions, initialValues, parameters=NULL, df } -# Overloaded plot function for denim object +#' Overloaded plot function for denim object +#' @param x - output of denim::sim function +#' +#' @param ... - additional parameter for `plot()` function +#' @param color_palette - a palette name from the colorspace package. You can view available palettes with `colorspace::hcl_palettes("qualitative", plot = TRUE)`. +#' +#' @import colorspace +#' @method plot denim #' @export -plot.denim <- function(x, ...) { +plot.denim <- function(x, ..., color_palette=NULL) { # Set color codes and compartment names - col_codes <- viridisLite::viridis(ncol(x) - 1) - # Test different color scale to be discrete instead of continuous - # col_codes <- hue_pal(l = 65, c = 70)(ncol(x) - 1) - comp_names <- colnames(x)[-1] + # col_codes <- viridisLite::viridis(ncol(x) - 1) + col_codes <- colorspace::qualitative_hcl(ncol(x) - 1, palette = color_palette) + + comp_names <- colnames(x)[-1] + # Plot the first compartment cmd1 <- paste0("with(x, { plot(Time, ", comp_names[1], ", type = \"l\", lwd = 3, col = \"", col_codes[1], diff --git a/README.Rmd b/README.Rmd index 939d41c..dd9aac2 100644 --- a/README.Rmd +++ b/README.Rmd @@ -24,7 +24,13 @@ An R package for building and simulating deterministic discrete-time compartment ## Installation -You can install the development version of denim from [GitHub](https://github.com/) with: +You can install denim from CRAN with: + +``` r +install.packages("denim") +``` + +Or install the development version of denim from [GitHub](https://github.com/) with: ``` r # install.packages("devtools") @@ -38,13 +44,13 @@ This is a basic example to illustrate the specification of a simple SIR model, w ```{r example} library(denim) -transitions <- list( - "S -> I" = "beta * S * I / N", - "I -> R" = d_gamma(rate = 1/3, shape = 2) -) +transitions <- denim_dsl({ + S -> I = beta * S * I / N * timeStep + I -> R = d_gamma(rate = 1/3, shape = 2) +}) parameters <- c( - beta = 0.12, + beta = 1.2, N = 1000 ) @@ -54,7 +60,7 @@ initialValues <- c( R = 0 ) -simulationDuration <- 10 +simulationDuration <- 20 timeStep <- 0.01 mod <- sim(transitions = transitions, initialValues = initialValues, @@ -71,5 +77,5 @@ head(mod) We can plot the output with: ```{r example-plot, dpi=300, fig.width=5} -plot(mod) +plot(mod, ylim = c(1, 1000)) ``` diff --git a/README.md b/README.md index b53850f..53eb8a9 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,13 @@ compartmental models with memory. ## Installation -You can install the development version of denim from +You can install denim from CRAN with: + +``` r +install.packages("denim") +``` + +Or install the development version of denim from [GitHub](https://github.com/) with: ``` r @@ -36,13 +42,13 @@ are gamma distributed in this example: ``` r library(denim) -transitions <- list( - "S -> I" = "beta * S * I / N", - "I -> R" = d_gamma(rate = 1/3, shape = 2) -) +transitions <- denim_dsl({ + S -> I = beta * S * I / N * timeStep + I -> R = d_gamma(rate = 1/3, shape = 2) +}) parameters <- c( - beta = 0.12, + beta = 1.2, N = 1000 ) @@ -52,7 +58,7 @@ initialValues <- c( R = 0 ) -simulationDuration <- 10 +simulationDuration <- 20 timeStep <- 0.01 mod <- sim(transitions = transitions, initialValues = initialValues, @@ -66,17 +72,17 @@ The output is a data frame with 4 columns: `Time`, `S`, `I` and `R` head(mod) #> Time S I R #> 1 0.00 999.0000 1.000000 0.000000e+00 -#> 2 0.01 998.8801 1.119874 5.543225e-06 -#> 3 0.02 998.7459 1.254092 2.278823e-05 -#> 4 0.03 998.5956 1.404364 5.306419e-05 -#> 5 0.04 998.4273 1.572606 9.785981e-05 -#> 6 0.05 998.2389 1.760961 1.588423e-04 +#> 2 0.01 998.9880 1.011982 5.543225e-06 +#> 3 0.02 998.9759 1.024097 2.219016e-05 +#> 4 0.03 998.9636 1.036346 5.000038e-05 +#> 5 0.04 998.9512 1.048730 8.903457e-05 +#> 6 0.05 998.9386 1.061252 1.393545e-04 ``` We can plot the output with: ``` r -plot(mod) +plot(mod, ylim = c(1, 1000)) ``` ![](man/figures/README-example-plot-1.png) diff --git a/_pkgdown.yml b/_pkgdown.yml index bbd238a..03d90ba 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -14,11 +14,12 @@ articles: - title: denim features navbar: ~ contents: + - denim_dsl - multinomial + - nonparametric - title: denim comparison navbar: denim comparison contents: - - denim_vs_deSolve - denim_benchmark - title: migrate to denim navbar: migrate to denim @@ -48,9 +49,13 @@ reference: # desc: Define a set of probabilities of transition from one compartment to multiple compartments # contents: # - multinomial +- title: DSL parser + contents: + - denim_dsl - title: Simulator contents: - sim + - plot.denim - title: internal contents: - denim-package diff --git a/cran-comments.md b/cran-comments.md index 1059dc3..6605c8d 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,8 +1,7 @@ ## Resubmission This is a resubmission. In this version I have: -* Change parameter scale of d_gamma to rate -* Fix bug in multinomial and competing risks modeling -* Add the functionality to distribute the initial values according to the specified distribution +* Improve the run time for sim() function +* Add option to select color palette for plot() ## R CMD check results @@ -26,5 +25,3 @@ We checked 0 reverse dependencies, comparing R CMD check results across CRAN and * We saw 0 new problems * We failed to check 0 packages - -* This is a new release. diff --git a/docs/404.html b/docs/404.html index dd3c82f..5dc9cd9 100644 --- a/docs/404.html +++ b/docs/404.html @@ -31,7 +31,7 @@ denim - 1.0.1 + 1.2.1 + + + + + +
+ + + + +
+
+

denim DSL

+ + + Source: vignettes/denim_dsl.Rmd +
denim_dsl.Rmd
+
+ + + + +
+

Model definition in denim +

+

In denim, model is defined by a set of transitions between +compartments. Each transition is provided in the form of a +key-value pair, where:

+
    +

  • key show the transition direction between 2 +compartments.

    • +

  • value is an expression that describe the transition, it +can be a math expression or a +built-in dwell-time distribution +function.

    • +
+

These key-value pairs can be provided in 2 ways

+
    +

  • Using denim domain-specific language (DSL).

    • +

  • Define as a list in R.

    • +
+
+
+

Denim DSL +

+

In denim, each line of code must be a transition. The syntax for +defining a transition in denim DSL is as followed:

+

compartment_A -> compartment_B = transition

+

Model definition written in denim DSL must be parsed by the function +denim_dsl()

+

Math expression

+

For math expression, some basic supported operators include: ++ for addition, - for minus, * +for multiplication, / for division, ^ for +power. Users can also define additional model parameters in the math +expression.

+

Math expressions in denim are parsed using muparser. For a full list +of operators, visit the muparser website at https://beltoforion.de/en/muparser/features.php.

+

Distribution functions

+

Several built-in functions are provided to describe transitions based +on the distribution of dwell time:

+ +

Each of these functions accepts either fixed numeric values or +model parameters as inputs for their distributional parameters. +Note that current version of denim does not +accept mathematical expression as a distribution’s +parameters.

+
+

Define a classic SIR model +

+

A classic SIR model can be defined in denim as followed

+
+sir_model <- denim_dsl({
+  S -> I = beta * (I/N) * S * timeStep
+  I -> R = d_exponential(rate = gamma)
+})
+

Any variable on the right hand side (RHS) of the transitions +definition (that are not compartments) will be considered model +parameters and their values must be provided later on.

+

In this example, model parameters are: N, +beta, S, gamma.

+

timeStep is a special variable in denim that +automatically uses the time step defined in simulation configuration as +it’s value. Note that this special variable can ONLY be used within math +expression.

+

Users can also choose to provide fixed values as the distributional +parameter as followed

+
+sir_model <- denim_dsl({
+  S -> I = beta*(I/N)*S*timeStep
+  I -> R = d_exponential(rate = 1/4)
+})
+

Similar to R, users can also add comments in denim DSL by starting +the comment with # sign.

+
+sir_model <- denim_dsl({
+  # this is a comment
+  S -> I = beta*(I/N)*S*timeStep
+  I -> R = d_exponential(rate = 1/4) # this is another comment
+})
+

Run model

+

To run the model, users must provide:

+
    +

  • Values for model parameters (in this example, N, +beta, S and gamma).

    • +

  • Initial population for the compartments.

    • +

  • Simulation configurations.

    • +
+

Parameters and initial values can be defined as named vectors +or named lists in R.

+
+# parameters for the model
+parameters <- c(
+  beta = 0.4,
+  N = 1000,
+  gamma = 1/7
+)
+# initial population for each compartment 
+initValues <- c(
+  S = 999, 
+  I = 50,
+  R = 0
+)
+

Simulation configurations are provided as parameters for +sim() function which runs the model, where:

+
    +

  • timeStep is the duration of each time step in the +model.

    • +

  • simulationDuration is the duration to run the +simulation.

    • +
+
+mod <- sim(sir_model, 
+    parameters = parameters, 
+    initialValues = initValues, 
+    timeStep = 0.01,
+    simulationDuration = 40)
+
+plot(mod, ylim = c(1, 1000))
+

+
+
+

Time varying transition +

+

Aside from timeStep, denim also have another special +variable time for time varying transition (e.g. for +modeling seasonality). Note that this variable can ONLY be used within +math expression.

+

Example: time varying transition

+
+time_varying_mod <- denim_dsl({
+  A -> B = 20 * (1+cos(omega * time)) * timeStep
+})
+
+# parameters for the model
+parameters <- c(
+  omega = 2*pi/10
+)
+# initial population for each compartment 
+initValues <- c(A = 1000, B = 0)
+
+mod <- sim(time_varying_mod, 
+    parameters = parameters, 
+    initialValues = initValues, 
+    timeStep = 0.01,
+    simulationDuration = 40)
+
+plot(mod, ylim = c(0, 1000))
+

+
+
+
+

R list +

+

Users can also define the model structure as a list in R. For +example, the SIR model from previous example can be represented as +followed.

+
+sir_model_list <- list(
+  "S -> I" = "beta * (I/N) * S * timeStep",
+  "I -> R" = d_exponential(rate = "gamma")
+)
+
+sir_model_list
+#> $`S -> I`
+#> [1] "beta * (I/N) * S * timeStep"
+#> 
+#> $`I -> R`
+#> Discretized exponential distribution
+#> Rate = gamma
+

Note that the transitions (S -> I, +I -> R), mathematical expression +(beta * (I/N) * S * timeStep), and the model parameter +(gamma) must now be provided as strings.

+

We can then run the model in the same manner as previously +demonstrated.

+
+# parameters for the model
+parameters <- c(
+  beta = 0.4,
+  N = 1000,
+  gamma = 1/7
+)
+# initial population for each compartment 
+initValues <- c(
+  S = 999, 
+  I = 50,
+  R = 0
+)
+# run the simulation
+mod <- sim(sir_model_list, 
+    parameters = parameters, 
+    initialValues = initValues, 
+    timeStep = 0.01,
+    simulationDuration = 40)
+# plot output
+plot(mod, ylim = c(1, 1000))
+

+

When should I define model as a list in +R?

+

While denim DSL offers cleaner and more readable syntax to define +model structure, using R list may be more familiar to R users and better +suited for integration to a more R-centric workflow.

+

For example, consider a use case below, where we explore how model +dynamics change under three different I -> R dwell time +distributions (d_gamma, d_weibull, +d_lognormal) using map2.

+
+library(tidyverse)
+#> Warning: package 'ggplot2' was built under R version 4.3.1
+#> Warning: package 'tidyr' was built under R version 4.3.1
+#> Warning: package 'readr' was built under R version 4.3.1
+#> Warning: package 'purrr' was built under R version 4.3.3
+#> Warning: package 'dplyr' was built under R version 4.3.1
+#> Warning: package 'stringr' was built under R version 4.3.1
+#> Warning: package 'lubridate' was built under R version 4.3.3
+# configurations for 3 different I->R transitions
+model_config <- tibble(
+  IR_dists = c(d_gamma, d_weibull, d_lognormal),
+  IR_pars = list(c(rate = 0.1, shape = 3), c(scale = 5, shape = 0.3), c(mu = 0.3, sigma = 2))
+)
+
+model_config %>% 
+  mutate(
+    plots = map2(IR_dists, IR_pars, \(dist, par){
+      transitions <- list(
+        "S -> I" = "beta * S * (I / N) * timeStep",
+        # This is not applicable when using denim_dsl()
+        "I -> R" = do.call(dist, as.list(par))
+      )
+      
+      # model settings
+      denimInitialValues <- c(S = 980, I = 20, R = 0)
+      parameters <- c(
+        beta = 0.4,
+        N = 1000
+      )
+      
+      # compare output 
+      mod <- sim(transitions = transitions, 
+                 initialValues = denimInitialValues, 
+                 parameters = parameters, 
+                 simulationDuration = 60, 
+                 timeStep = 0.05)
+      
+      plot(mod, ylim = c(0,1000))
+    })
+  ) %>% 
+  pull(plots)
+

+
#> [[1]]
+#> [[1]]$rect
+#> [[1]]$rect$w
+#> [1] 8.5166
+#> 
+#> [[1]]$rect$h
+#> [1] 324.2026
+#> 
+#> [[1]]$rect$left
+#> [1] 53.8834
+#> 
+#> [[1]]$rect$top
+#> [1] 662.1013
+#> 
+#> 
+#> [[1]]$text
+#> [[1]]$text$x
+#> [1] 60.30808 60.30808 60.30808
+#> 
+#> [[1]]$text$y
+#> [1] 581.0507 500.0000 418.9493
+#> 
+#> 
+#> 
+#> [[2]]
+#> [[2]]$rect
+#> [[2]]$rect$w
+#> [1] 8.5166
+#> 
+#> [[2]]$rect$h
+#> [1] 324.2026
+#> 
+#> [[2]]$rect$left
+#> [1] 53.8834
+#> 
+#> [[2]]$rect$top
+#> [1] 662.1013
+#> 
+#> 
+#> [[2]]$text
+#> [[2]]$text$x
+#> [1] 60.30808 60.30808 60.30808
+#> 
+#> [[2]]$text$y
+#> [1] 581.0507 500.0000 418.9493
+#> 
+#> 
+#> 
+#> [[3]]
+#> [[3]]$rect
+#> [[3]]$rect$w
+#> [1] 8.5166
+#> 
+#> [[3]]$rect$h
+#> [1] 324.2026
+#> 
+#> [[3]]$rect$left
+#> [1] 53.8834
+#> 
+#> [[3]]$rect$top
+#> [1] 662.1013
+#> 
+#> 
+#> [[3]]$text
+#> [[3]]$text$x
+#> [1] 60.30808 60.30808 60.30808
+#> 
+#> [[3]]$text$y
+#> [1] 581.0507 500.0000 418.9493
+
+
+
+ + + +
+

Developed by Thinh Ong, Anh Phan, Marc Choisy.

+
+ +
+

Site built with pkgdown 2.0.9.

+
+ +
+
+ + + + + + + diff --git a/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-10-1.png new file mode 100644 index 0000000..f0d0618 Binary files /dev/null and b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-10-1.png differ diff --git a/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-1.png b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-1.png new file mode 100644 index 0000000..8d3499c Binary files /dev/null and b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-1.png differ diff --git a/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-2.png b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-2.png new file mode 100644 index 0000000..01c5645 Binary files /dev/null and b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-2.png differ diff --git a/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-3.png b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-3.png new file mode 100644 index 0000000..2694e87 Binary files /dev/null and b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-11-3.png differ diff --git a/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-7-1.png b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-7-1.png new file mode 100644 index 0000000..f0d0618 Binary files /dev/null and b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-7-1.png differ diff --git a/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-8-1.png b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-8-1.png new file mode 100644 index 0000000..104e344 Binary files /dev/null and b/docs/articles/denim_dsl_files/figure-html/unnamed-chunk-8-1.png differ diff --git a/docs/articles/denim_files/figure-html/unnamed-chunk-11-1.png b/docs/articles/denim_files/figure-html/unnamed-chunk-11-1.png new file mode 100644 index 0000000..7f8f84d Binary files /dev/null and b/docs/articles/denim_files/figure-html/unnamed-chunk-11-1.png differ diff --git a/docs/articles/denim_files/figure-html/unnamed-chunk-16-1.png b/docs/articles/denim_files/figure-html/unnamed-chunk-16-1.png index 4bed4da..9e3871b 100644 Binary files a/docs/articles/denim_files/figure-html/unnamed-chunk-16-1.png and b/docs/articles/denim_files/figure-html/unnamed-chunk-16-1.png differ diff --git a/docs/articles/denim_files/figure-html/unnamed-chunk-18-1.png b/docs/articles/denim_files/figure-html/unnamed-chunk-18-1.png new file mode 100644 index 0000000..2f44f8c Binary files /dev/null and b/docs/articles/denim_files/figure-html/unnamed-chunk-18-1.png differ diff --git a/docs/articles/index.html b/docs/articles/index.html index f0158a4..e6d7d94 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -10,7 +10,7 @@ denim - 1.0.1 + 1.2.1 + + + + + +
+ + + + +
+
+

nonparametric distribution

+ + + Source: vignettes/nonparametric.Rmd +
nonparametric.Rmd
+
+ + + + +
+

Non-parametric vs parametric +

+

In denim, users have 2 options to define dwell-time distribution:

+
    +

  • As a parametric distribution: using d_*() +functions.

    • +

  • As a non-parametric distribution: using +nonparametric() function and user must provide the +histogram of distribution where bin-width matches +timeStep.

    • +
+
+

Example 1 +

+

To demonstrate the difference between 2 approaches, we can try +modeling an SIR model with Weibull distributed infectious period using +nonparametric() and d_weibull().

+

Model definition using +d_weibull()

+
+sir_parametric <- denim_dsl({
+  S -> I = beta * (I/N) * S * timeStep
+  I -> R = d_weibull(scale = r_scale, shape = r_shape)
+})
+

Model parameters that must be defined are: beta, +N, r_scale, r_shape

+

Model definition using +nonparametric()

+
+sir_nonparametric <- denim_dsl({
+  S -> I = beta * (I/N) * S * timeStep
+  I -> R = nonparametric(dwelltime_dist)
+})
+

Model parameters that must be defined are: beta, +N, dwelltime_dist (the discrete dwell time +distribution)

+

Run model

+

We will run both models under the following model settings

+
+# parameters
+mod_params <- list(
+  beta = 0.4,
+  N = 1000,
+  r_scale = 4,
+  r_shape = 3
+)
+# initial population
+init_vals <- c(S = 950, I = 50, R = 0)
+# simulation duration and timestep
+sim_duration <- 30
+timestep <- 0.05
+

Running the model with d_weibull() is straight +forward

+
+parametric_mod <- sim(sir_parametric,
+    initialValues = init_vals,
+    parameters = mod_params,
+    simulationDuration = sim_duration,
+    timeStep = timestep) 
+
+plot(parametric_mod, ylim = c(0, 1000))
+

+

However, to run the model using nonparametric(), we +first need to compute the discrete dwell time distribution +(dwelltime_dist).

+

Since all parametric distributions are asymptotic to 1, we will set +the maximal dwell time as the time point where the cumulative +probability is sufficiently close to 1 (i.e. above the threshold +1 - error_tolerance).

+

A helper function to compute discrete dwell time distribution from a +distribution function in R is provided below.

+
+Helper function +
+# Compute discrete distribution of dwell-tinme
+# dist_func - R distribution function for dwell time (pexp, pgamma, etc.)
+# ... - parameters for dist_func
+compute_dist <- function(dist_func,..., timestep=0.05, error_tolerance=0.0001){
+  maxtime <- timestep
+  prev_prob <- 0
+  prob_dist <- numeric()
+  
+  while(TRUE){
+     # get current cumulative prob and check whether it is sufficiently close to 1
+     temp_prob <-  ifelse(
+       dist_func(maxtime, ...) < (1 - error_tolerance), 
+       dist_func(maxtime, ...), 
+       1);
+
+     # get f(t)
+     curr_prob <- temp_prob - prev_prob
+     prob_dist <- c(prob_dist, curr_prob)
+     
+     prev_prob <- temp_prob
+     maxtime <- maxtime + timestep
+     
+     if(temp_prob == 1){
+       break
+     }
+  }
+  
+  prob_dist
+}
+

We can then run the model as followed

+
+# Compute the discrete distribution
+dwelltime_dist <- compute_dist(pweibull, 
+                               scale = mod_params$r_scale, shape = mod_params$r_shape,
+                               timestep = timestep)
+
+# Compute the discrete distribution
+nonparametric_mod <- sim(sir_nonparametric,
+    initialValues = init_vals,
+    parameters = list(
+      beta = mod_params$beta,
+      N = mod_params$N,
+      dwelltime_dist = dwelltime_dist
+    ),
+    simulationDuration = sim_duration,
+    timeStep = timestep) 
+plot(nonparametric_mod, ylim = c(0, 1000))
+

+
+
+

Example 2 +

+

By using nonparametric(), we can run the model with any +dwell time distribution shape

+

Consider the following multimodal distribution.

+
+timestep <- 0.05
+plot(seq(0, by = 0.05, length.out = length(multimodal_dist)), 
+     multimodal_dist, 
+     type = "l", col = "#374F77", lty = 1, lwd = 3,
+     xlab = "Length of stay (days)", ylab = "", yaxt = 'n')
+

+

We can also run the sir_nonparametric model from last +example with this dwell time distribution

+
+# model parameter
+parameters <- list(beta = 0.4, N = 1000,
+                   dwelltime_dist = multimodal_dist)
+# initial population
+init_vals <- c(S = 950, I = 50, R = 0)
+# simulation duration and timestep
+sim_duration <- 30
+timestep <- 0.05
+
+# Run the model with multimodel distribution
+nonparametric_mod <- sim(
+  sir_nonparametric,
+  initialValues = init_vals,
+  parameters = parameters,
+  simulationDuration = sim_duration,
+  timeStep = timestep) 
+plot(nonparametric_mod, ylim = c(0, 1000))
+

+
+
+
+
+ + + +
+

Developed by Thinh Ong, Anh Phan, Marc Choisy.

+
+ +
+

Site built with pkgdown 2.0.9.

+
+ +
+
+ + + + + + + diff --git a/docs/articles/nonparametric_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-10-1.png new file mode 100644 index 0000000..d8b9d4d Binary files /dev/null and b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-10-1.png differ diff --git a/docs/articles/nonparametric_files/figure-html/unnamed-chunk-5-1.png b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-5-1.png new file mode 100644 index 0000000..ad75367 Binary files /dev/null and b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-5-1.png differ diff --git a/docs/articles/nonparametric_files/figure-html/unnamed-chunk-7-1.png b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-7-1.png new file mode 100644 index 0000000..89afd6d Binary files /dev/null and b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-7-1.png differ diff --git a/docs/articles/nonparametric_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-9-1.png new file mode 100644 index 0000000..900145d Binary files /dev/null and b/docs/articles/nonparametric_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/authors.html b/docs/authors.html index 1fb5619..04d48b2 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -10,7 +10,7 @@ denim - 1.0.1 + 1.2.1 + + + + + +
+
+
+

Define transitions using denim's domain-specific language (DSL)

+ Source: R/dsl.R +
denim_dsl.Rd
+
+ +
+

This function parses model transitions defined in denim's DSL syntax

+
+ +
+

Usage

+ 
denim_dsl(x)
+
+ +
+

Arguments

+
x
+

  • an expression written in denim's DSL syntax. Each line should be a transition written in the +format compartment -> out_compartment = expression where expression can be either a math expression +or one of denim's built-in dwell time distribution function

    • +
+ +
+
+

Value

+ + +

denim_transition object

+
+ +
+

Examples

+ 
transitions <- denim_dsl({
+  S -> I = beta * (I/N) * S * timeStep
+  I -> R = d_gamma(rate = 1/4, shape = 3)
+})
+
+
+
+ + +
+

Developed by Thinh Ong, Anh Phan, Marc Choisy.

+
+ +
+

Site built with pkgdown 2.0.9.

+
+ +
+ + + + + + + diff --git a/docs/reference/figures/README-example-plot-1.png b/docs/reference/figures/README-example-plot-1.png index 1efcd30..cfe1dc1 100644 Binary files a/docs/reference/figures/README-example-plot-1.png and b/docs/reference/figures/README-example-plot-1.png differ diff --git a/docs/reference/index.html b/docs/reference/index.html index a273a0d..bbcca0a 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -10,7 +10,7 @@ denim - 1.0.1 + 1.2.1 + + + + + +
+
+
+

Overloaded plot function for denim object

+ Source: R/simulators.R +
plot.denim.Rd
+
+ +
+

Overloaded plot function for denim object

+
+ +
+

Usage

+ 
# S3 method for denim
+plot(x, ..., color_palette = NULL)
+
+ +
+

Arguments

+
x
+

  • output of denim::sim function

    • +
+ + +
...
+

  • additional parameter for plot() function

    • +
+ + +
color_palette
+

  • a palette name from the colorspace package. You can view available palettes with colorspace::hcl_palettes("qualitative", plot = TRUE).

    • +
+ +
+ +
+ + +
+

Developed by Thinh Ong, Anh Phan, Marc Choisy.

+
+ +
+

Site built with pkgdown 2.0.9.

+
+ +
+ + + + + + + diff --git a/docs/reference/sim.html b/docs/reference/sim.html index 1ee7691..7672256 100644 --- a/docs/reference/sim.html +++ b/docs/reference/sim.html @@ -10,7 +10,7 @@ denim - 1.0.1 + 1.2.1