[[http://www.gnu.org/licenses/gpl-3.0.html][http://img.shields.io/:license-gpl3-blue.svg]]
News
The microsimulation package now allows for in-line use and provides a new SummaryReport class. See the in-line example below for a demonstration of both of these facilities.
The C++ application programming interface is now [[https://htmlpreview.github.io/?https://github.com/mclements/microsimulation/blob/master/inst/doc/html/index.html][described]] using Doxygen.
Background [[https://en.wikipedia.org/wiki/Microsimulation][Microsimulation]] is a form of individual-based stochastic simulation. In continuous time, microsimulation is closely related to [[https://en.wikipedia.org/wiki/Discrete_event_simulation][discrete event simulation]], and in discrete time it is closely related to [[https://en.wikipedia.org/wiki/Agent-based_model][agent-based models]]. In econometrics and health care, microsimulation is often used to model policy changes. Our implementation is in continuous time and uses event-based discrete event simulation for the model specification.
The package provides several approaches for microsimulation and event-based, discrete event simulation. The package includes an R implementation of discrete event simulation, building on several R5 classes. This implementation is useful from a pedagogical perspective, but it is slow for larger microsimulations. For speed, we also provide C++ classes for discrete event simulation, building on several well established libraries, including the [[http://www.inf.usi.ch/carzaniga/ssim/index.html][SSIM]] C++ discrete event simulation library, the [[http://www.iro.umontreal.ca/~lecuyer/myftp/streams00/][RngStream]] C library for common random numbers, the [[http://www.boost.org/][Boost]] libraries for making many C++11 tools available to C++98, and [[http://www.rcpp.org/][Rcpp]] for providing glue between R, R's C API and C++.
We specifically developed this package for modelling the cost-effectiveness of cancer screening, where there are many (e.g. 10^7) individuals who are followed from birth to death. Notably, we provide a complete prostate cancer screening model, including tools for cost-effectiveness analysis.
install.packages("microsimulation")
For compiling from source, use steps 2--3.
2 Dependencies: :: The microsimulation package requires [[http://www.rcpp.org/][Rcpp]]. A convenient, but not required, way of installing github-packages in R is to use [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]]. Since both of the dependencies and [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]] are available on [[https://cran.r-project.org/][CRAN]] just run the following in R.
install.packages("Rcpp") install.packages("devtools")
3a Installation with devtools: :: To install the microsimulation using [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]] just run the following in R:
require(devtools) install_github("mclements/microsimulation")
3b Alternative installation from shell: ::
If you prefer the shell over [[https://cran.r-project.org/web/packages/devtools/README.html][devtools]], just run the following to download the microsimulation R-package:
git clone https://github.com/mclements/microsimulation.git
To install the microsimulation R-package run this in your shell:
R CMD INSTALL path_to_microsimulation
** Simple examples
(concat "#> "(mapconcat 'identity (split-string result "\n") "\n#> "))
require(microsimulation) sim1 <- callSimplePerson2(n = 1e5) summary(sim1$events)
: #> state event age number : #> Healthy:338 toOtherDeath :142 Min. : 1.00 Min. : 1.0 : #> Cancer : 0 toCancer :100 1st Qu.: 37.00 1st Qu.: 74.0 : #> Death : 0 toCancerDeath: 96 Median : 58.00 Median : 325.0 : #> Mean : 56.98 Mean : 445.1 : #> 3rd Qu.: 79.00 3rd Qu.: 721.8 : #> Max. :100.00 Max. :1722.0 : #>
sim2 <- callIllnessDeath(n = 1e5, cure = 0.2) summary(sim2$prev)
: #> state age number : #> Healthy:200 Min. : 0.00 Min. : 1 : #> Cancer : 0 1st Qu.: 25.75 1st Qu.: 2030 : #> Median : 50.50 Median : 8096 : #> Mean : 50.49 Mean : 34914 : #> 3rd Qu.: 75.25 3rd Qu.: 79029 : #> Max. :100.00 Max. :100000 : #>
** In-line example: Sick-Sicker example from Krijkamp et al (2018)
We have provided a continuous-time re-implementation of the discrete-time Markov model developed by Krijkamp et al (2018).
The example demonstrates:
For the specific example, the transformation from discrete to continuous time showed that the discrete time formulation could have been improved. In particular, the discrete time formulation assumes no transitions between Healthy and Sicker over one year, while the approximate probability of that event is the one-year probability of moving from Healthy to Sick times the probability of moving from Sick to Sicker. We have included that probability in the transition matrix and using matrix logarithms to calculate the transition probabilities.
library(expm) # logm library(Rcpp) # sourceCpp library(microsimulation) # Rcpp::depends and include files library(ascii); options(asciiType="org")
param <- within(list(), {
p.HD = 0.005 # probability to die when healthy
p.HS1 = 0.15 # probability to become sick when healthy
p.S1H = 0.5 # probability to become healthy when sick
p.S1S2 = 0.105 # probability to become sicker when sick
rr.S1 = 3 # rate ratio of death when sick vs healthy
rr.S2 = 10 # rate ratio of death when sicker vs healthy
r.HD = -log(1-p.HD) # rate of death when healthy
r.S1D = rr.S1 * r.HD # rate of death when sick
r.S2D = rr.S2 * r.HD # rate of death when sicker
p.S1D = 1-exp(-r.S1D) # probability to die when sick
p.S2D = 1-exp(-r.S2D) # probability to die when sicker
## Cost and utility inputs
c_H = 2000 # cost of remaining one cycle healthy
c_S1 = 4000 # cost of remaining one cycle sick
c_S2 = 15000 # cost of remaining one cycle sicker
c_Trt = 12000 # (additional) cost of treatment (per cycle)
u_H = 1 # utility when healthy
u_S1 = 0.75 # utility when sick
u_S2 = 0.5 # utility when sicker
u_Trt = 0.95 # utility when sick (as per the code) and being treated
## new parameters
discountRate = 0.03 # discount rate
partitionBy = 1.0 # partition used in the report
Trt = FALSE # Treatment?
debug = FALSE})
param = within(param, { p.HS2 = p.HS1p.S1S2 }) Pmat = with(param, matrix(c(1-p.HD-p.HS1-p.HS2,p.HS1,p.HS2,p.HD, p.S1H,1-p.S1H-p.S1S2-p.S1D,p.S1S2,p.S1D, 0,0,1-p.S2D,p.S2D, 0,0,0,1), 4, byrow=TRUE)) stopifnot(all(abs(rowSums(Pmat)-1)<10.Machine$double.eps)) Qmat = expm::logm(Pmat) # matrix logarithm stopifnot(all(abs(rowSums(Qmat))<10*.Machine$double.eps))
param = within(param, { r_HS1 = Qmat[1,2]; r_HD = Qmat[1,4] r_S1H = Qmat[2,1]; r_S1S2 = Qmat[2,3]; r_S1D = Qmat[2,4] r_S2D = Qmat[3,4] })
sourceCpp(code=" // [[Rcpp::depends(BH)]] // [[Rcpp::depends(RcppArmadillo)]] // [[Rcpp::depends(microsimulation)]]
enum state_t {Healthy, Sick, Sicker, Dead};
enum event_t {toS1, toS2, toH, toD, toEOF};
typedef ssim::SummaryReport<short,short> Report;
/
Utility: Random exponential using rate parameterisation
*/
template
ascii.proc_time = function(x, include.rownames=FALSE, include.colnames=TRUE, ...) ascii(summary(x), include.rownames, include.colnames, ...)
set.seed(12345) ascii(system.time(sim1 <- simulations(1e3, param=modifyList(param, list(Trt = FALSE)))),header=TRUE) set.seed(12345) ascii(system.time(sim2 <- simulations(1e3, param=modifyList(param, list(Trt = TRUE)))),FALSE,FALSE) cat("\n") ascii(ICER(sim1,sim2), rownames=c("No treatment","Treatment"), caption="Continuous-time Sick-Sicker model for n=10,000 individuals")
:RESULTS: | user | system | elapsed | |------+--------+---------| | 0.10 | 0.00 | 0.10 | | 0.09 | 0.00 | 0.09 |
| | Total | | | | Incremental | | | | | | | Costs | (se) | QALYs | (se) | Costs | (se) | QALYs | (se) | ICER | |--------------+----------+--------+--------+-------+---------------+--------+-------+-------+----------| | No treatment | 117835.7 | 2263.0 | 13.101 | 0.160 | | | | | | | Treatment | 223871.9 | 4097.8 | 13.658 | 0.167 | 106036.1 | 1881.9 | 0.558 | 0.013 | 190110.5 | :END:
library(parallel) set.seed(12345) mcsimulations <- function(n, simulations, ..., mc.cores = getOption("mc.cores", 2L)) { n.seg <- diff(c((0:(mc.cores-1))*floor(n/mc.cores),n)) do.call(rbind, mclapply(n.seg, simulations, ..., mc.cores=mc.cores)) } ascii(system.time(sim1 <- mcsimulations(1e5, simulations, param=modifyList(param, list(Trt = FALSE)))), header=TRUE) set.seed(12345) ascii(system.time(sim2 <- mcsimulations(1e5, simulations, param=modifyList(param, list(Trt = TRUE)))), FALSE,FALSE) cat("\n") ascii(ICER(sim1,sim2),caption="Continuous-time Sick-Sicker model for n=100,000 individuals")
:RESULTS: | user | system | elapsed | |------+--------+---------| | 8.25 | 0.03 | 4.14 | | 4.18 | 0.02 | 4.20 |
| | Total | | | | Incremental | | | | | | | Costs | (se) | QALYs | (se) | Costs | (se) | QALYs | (se) | ICER | |-----------+----------+-------+--------+-------+---------------+-------+-------+-------+----------| | Reference | 118425.3 | 234.4 | 13.133 | 0.016 | | | | | | | Treatment | 225548.4 | 424.3 | 13.720 | 0.017 | 107123.1 | 485.7 | 0.587 | 0.023 | 182387.6 | :END:
There was an appreciable difference in the estimates from the discrete-time case and the continuous-time case. This issue warrants further investigation, particularly given that the discrete time case ignores any transitions between Healthy and Sicker within a one-year period.
| | Total | | | | Incremental | | | | |
| | Costs | (se) | QALYs | (se) | Costs | (se) | QALYs | (se) | ICER |
|-----------+---------+------+-------+-------+---------------+------+-------+-------+--------|
| | | | | |
set.seed(12345) invisible(callSim(5, modifyList(param, list(debug=TRUE)))) set.seed(12345) invisible(callSim(5, modifyList(param, list(Trt=TRUE, debug=TRUE))))
id: 0, state: 0, kind: 0, previous: 0.000000, now: 1.554306 id: 0, state: 1, kind: 1, previous: 1.554306, now: 1.660357 id: 0, state: 2, kind: 3, previous: 1.660357, now: 2.137922 id: 1, state: 0, kind: 0, previous: 0.000000, now: 22.523279 id: 1, state: 1, kind: 2, previous: 22.523279, now: 23.538790 id: 1, state: 0, kind: 0, previous: 23.538790, now: 25.110528 id: 1, state: 1, kind: 2, previous: 25.110528, now: 25.529621 id: 1, state: 0, kind: 4, previous: 25.529621, now: 31.000000 id: 2, state: 0, kind: 0, previous: 0.000000, now: 1.279403 id: 2, state: 1, kind: 2, previous: 1.279403, now: 2.482209 id: 2, state: 0, kind: 0, previous: 2.482209, now: 4.467754 id: 2, state: 1, kind: 2, previous: 4.467754, now: 5.240185 id: 2, state: 0, kind: 0, previous: 5.240185, now: 6.123519 id: 2, state: 1, kind: 2, previous: 6.123519, now: 7.951863 id: 2, state: 0, kind: 0, previous: 7.951863, now: 10.967305 id: 2, state: 1, kind: 2, previous: 10.967305, now: 11.127479 id: 2, state: 0, kind: 0, previous: 11.127479, now: 15.295999 id: 2, state: 1, kind: 1, previous: 15.295999, now: 16.663526 id: 2, state: 2, kind: 4, previous: 16.663526, now: 31.000000 id: 3, state: 0, kind: 0, previous: 0.000000, now: 0.021088 id: 3, state: 1, kind: 2, previous: 0.021088, now: 1.648953 id: 3, state: 0, kind: 0, previous: 1.648953, now: 1.665451 id: 3, state: 1, kind: 1, previous: 1.665451, now: 2.825812 id: 3, state: 2, kind: 3, previous: 2.825812, now: 3.378253 id: 4, state: 0, kind: 0, previous: 0.000000, now: 1.548766 id: 4, state: 1, kind: 2, previous: 1.548766, now: 4.303079 id: 4, state: 0, kind: 0, previous: 4.303079, now: 5.128089 id: 4, state: 1, kind: 2, previous: 5.128089, now: 5.418233 id: 4, state: 0, kind: 0, previous: 5.418233, now: 9.734630 id: 4, state: 1, kind: 2, previous: 9.734630, now: 11.013894 id: 4, state: 0, kind: 0, previous: 11.013894, now: 26.388802 id: 4, state: 1, kind: 2, previous: 26.388802, now: 26.514306 id: 4, state: 0, kind: 0, previous: 26.514306, now: 27.803678 id: 4, state: 1, kind: 2, previous: 27.803678, now: 28.437418 id: 4, state: 0, kind: 4, previous: 28.437418, now: 31.000000 id: 0, state: 0, kind: 0, previous: 0.000000, now: 1.554306 id: 0, state: 1, kind: 1, previous: 1.554306, now: 1.660357 id: 0, state: 2, kind: 3, previous: 1.660357, now: 2.137922 id: 1, state: 0, kind: 0, previous: 0.000000, now: 22.523279 id: 1, state: 1, kind: 2, previous: 22.523279, now: 23.538790 id: 1, state: 0, kind: 0, previous: 23.538790, now: 25.110528 id: 1, state: 1, kind: 2, previous: 25.110528, now: 25.529621 id: 1, state: 0, kind: 4, previous: 25.529621, now: 31.000000 id: 2, state: 0, kind: 0, previous: 0.000000, now: 1.279403 id: 2, state: 1, kind: 2, previous: 1.279403, now: 2.482209 id: 2, state: 0, kind: 0, previous: 2.482209, now: 4.467754 id: 2, state: 1, kind: 2, previous: 4.467754, now: 5.240185 id: 2, state: 0, kind: 0, previous: 5.240185, now: 6.123519 id: 2, state: 1, kind: 2, previous: 6.123519, now: 7.951863 id: 2, state: 0, kind: 0, previous: 7.951863, now: 10.967305 id: 2, state: 1, kind: 2, previous: 10.967305, now: 11.127479 id: 2, state: 0, kind: 0, previous: 11.127479, now: 15.295999 id: 2, state: 1, kind: 1, previous: 15.295999, now: 16.663526 id: 2, state: 2, kind: 4, previous: 16.663526, now: 31.000000 id: 3, state: 0, kind: 0, previous: 0.000000, now: 0.021088 id: 3, state: 1, kind: 2, previous: 0.021088, now: 1.648953 id: 3, state: 0, kind: 0, previous: 1.648953, now: 1.665451 id: 3, state: 1, kind: 1, previous: 1.665451, now: 2.825812 id: 3, state: 2, kind: 3, previous: 2.825812, now: 3.378253 id: 4, state: 0, kind: 0, previous: 0.000000, now: 1.548766 id: 4, state: 1, kind: 2, previous: 1.548766, now: 4.303079 id: 4, state: 0, kind: 0, previous: 4.303079, now: 5.128089 id: 4, state: 1, kind: 2, previous: 5.128089, now: 5.418233 id: 4, state: 0, kind: 0, previous: 5.418233, now: 9.734630 id: 4, state: 1, kind: 2, previous: 9.734630, now: 11.013894 id: 4, state: 0, kind: 0, previous: 11.013894, now: 26.388802 id: 4, state: 1, kind: 2, previous: 26.388802, now: 26.514306 id: 4, state: 0, kind: 0, previous: 26.514306, now: 27.803678 id: 4, state: 1, kind: 2, previous: 27.803678, now: 28.437418 id: 4, state: 0, kind: 4, previous: 28.437418, now: 31.000000
*** Limitations of the in-line approach
One limitation for the in-line code is that common random numbers, which are manipulated in C++ and use R's random number functions, are not available. Common random numbers can be used in a package, which is used by the [[https://github.com/mclements/prostata][prostata]] package.
** Extensive use case For more advance use of the microsimulation framework, please have a look at our prostate cancer natural history model: