Profiling the {SpeedyMarkov} implementation

Background

SpeedyMarkov uses a functional approach with the Markov model first being specified as a series of functions. These functions are then sampled using package functionality. Once samples from the model have been taken for each intervention SpeedyMarkov can then be used to simulate the Markov model and internally estimate total costs and QALYs, discounted forward in time. Finally, a separate function can be applied to produce cost effectiveness summary measures. These package steps are wrapped into a single function markov_ce_pipeline - exploring the documentation of this function and the code and documentation of the functions it itself wraps may be a useful approach to understanding the underlying mechanisms of the package.

In principle this functional approach may lead to slightly less efficient code as steps are separated from each other. On the other hand, it may aid generalisability, readability, and make optimising areas of the code easier (as individual functions can be profiled and optimised). A major aim of SpeedyMarkov is to provide an easy to use end-to-end tool chain. However, all code should be highly reusable if users wish to use it in a more modular fashion or as part of more complex workflows.

Profiling

R implementation

The R implementation of the SpeedyMarkov approach is roughly 25% faster than the reference implementation for this example with a substantially reduced memory footprint. However, some of this memory saving may be due to not saving individual simulations which may be desirable behaviour. The inner Markov for loop dominates the run-time - taking roughly 75% of the overall run-time. However, a substantial amount of time is also spent sampling data, summarising the findings and allocating samples to be simulated meaning that some optimisations may be possible here as well. Profiling is implemented using the profvis package (Note: Profiling output will be more detailed in an interactive session post package install).

Strategies for optimisation suggested for the reference application also hold here with the most pressing being the optimisation of the inner Markov loop.

library(profvis)
# Load SpeedyMarkov from local source.
devtools::load_all("..")
#> Loading SpeedyMarkov

profvis({
  markov_ce_pipeline(example_two_state_markov(), duration = 100, samples = 100000,
                     sim_type = "base")
})

Flame Graph

Data

Options ▾

<expr>		Memory				Time
	library(profvis)
	# Load SpeedyMarkov from local source.
	devtools::load_all("..")
	profvis({
	markov_ce_pipeline(example_two_state_markov(), duration = 100, samples = 100000,
	sim_type = "base")
	})

SpeedyMarkov/R/example_two_state_markov.R		Memory				Time
	#' Reference Two State Markov Model
	#' @description This is a two state Markov model - modelling smoking cessation - it was adapted from `reference_two_state_markov`
	#' to use the `SpeedyMarkov` framework. It essentially contains a list of functions that are then used to sample a markov model
	#' that can then be simulated and analysed. Unlike `reference_two_state_markov` this is not a standalone analysis pipeline
	#' but instead represents a model definition.
	#'
	#' @return A named list of functions that all require a samples argument and pass additional arguments (using ...).
	#' The list contains:
	#' * transitions_list: a list of transition functions, with the first taking the number of samples as an argument
	#' and the following being dependent on the a previous transition.
	#' * qalys: a function that samples the qaly cost for each intervention.
	#' * intervention_costs: a function that returns the costs for each intervention.
	#' * state_costs: a function that returns the state costs for each intervention.
	#' * cohorts: a function that returns the initial state for each intervention.
	#'
	#' Please see the code for more details on each required list item.
	#' @export
	#' @importFrom VGAM rdiric
	#' @importFrom stats rnorm
	#' @importFrom purrr map map2 transpose
	#' @author Sam Abbott
	#' @examples
	#' ## Example model run
	#' example_two_state_markov()
	#'
	example_two_state_markov <- function() {
	# Transitions -------------------------------------------------------------
	# 1. Specify transition matrices for each intervention
	# Baseline - Soc
	# Pass additional arguments internally
	soc_transition <- function(samples = NULL, ...) {
	# Sample transitions
	tmp <- list(VGAM::rdiric(samples, c(88, 12)),
	VGAM::rdiric(samples, c(8, 92)))
	# Arrange as matrices
	tmp <- SpeedyMarkov::matrix_arrange(tmp, ...)
	return(tmp)
	}
	# Intervention - Soc with website
	# Depends on Soc
	soc_with_website_transition <- function(baseline = NULL, ...) {
	#Sample transitions for each baseline matrix
	samples <- length(baseline)
	new_row_sample <- VGAM::rdiric(samples,c(85,15))
	# Update baseline
	updated <- purrr::map2(baseline, 1:samples, function(transition, sample) {
	transition[1, ] <- new_row_sample[sample, ]
	return(transition)
	})
	return(updated)
	}
	## Test
	#soc_trans_sample <- soc_transition()
	# soc_trans_sample
	#soc_with_website_trans_sample <- soc_with_website_transition(soc_trans_sample)
	# soc_with_website_trans_sample
	#Set up transition list
	transitions_list <- list(soc_transition,
	soc_with_website_transition)
	names(transitions_list) <- c("SoC", "Soc with Website")
	# Qualies -----------------------------------------------------------------
	# 2. Specify qaly costs per intervention (random sampling)
	qalys <- function(samples = NULL, ...) {
	qaly <- function(samples = 1, ...) {
	## Sample
	tmp <- list(stats::rnorm(samples, mean = 0.95,sd = 0.01) / 2,
	rep(1 / 2, samples))
	out <- SpeedyMarkov::vector_arrange(tmp)
	return(out)
	}
	soc <- qaly(samples = samples)
	soc_with_website <- soc
	out <- list(soc, soc_with_website)
	names(out) <- list("SoC", "Soc with Website")
	out <- purrr::transpose(out)
	return(out)
	}
	# qalys()
	# Costs -------------------------------------------------------------------
	# 3. Specify costs per intervention (random sampling)
	intervention_costs <- function(samples = NULL, ...) {
	## Sample
	tmp <- list(rep(0, samples),
	rep(50, samples))
	out <- SpeedyMarkov::vector_arrange(tmp)
	return(out)
	}
	# intervention_costs()
	state_costs <- function(samples = NULL, ...) {
	state_cost <- function(samples = 1) {
	tmp <- list(rep(0, samples),
	rep(0, samples))
	out <- SpeedyMarkov::vector_arrange(tmp)
	return(out)
	}
	soc <- state_cost(samples = samples)
	soc_with_website <- soc
	out <- list(soc, soc_with_website)
	names(out) <- list("SoC", "Soc with Website")
	out <- purrr::transpose(out)
	return(out)
	}
	# state_costs()
	# Cohort ------------------------------------------------------------------
	#4. Define cohort
	cohorts <- function(samples = NULL, ...) {
	cohort <- function(samples = 1) {
	tmp <- list(rep(1, samples),
	rep(0, samples))
	out <- SpeedyMarkov::vector_arrange(tmp)
	return(out)
	}
	soc <- cohort(samples = samples)
	soc_with_website <- soc
	out <- list(soc, soc_with_website)
	names(out) <- list("SoC", "Soc with Website")
	out <- purrr::transpose(out)
	return(out)
	}
	#cohorts()
	model <- list(
	transitions_list = transitions_list,
	qalys = qalys,
	intervention_costs = intervention_costs,
	state_costs = state_costs,
	cohorts = cohorts
	)
	class(model) <- c("SpeedyMarkov", class(model))
	return(model)
	}

SpeedyMarkov/R/sample_markov.R		Memory				Time
	#' Sample a Markov Model
	#'
	#'
	#'@description This model agnostic function samples a markov model specification. It wraps multiple approaches
	#'that may offer various advantages and disadvantages.
	#'
	#' @param markov_model A list of functions that define a markov model across multiple interventions. See `example_two_state_markov`
	#' for the correct format.
	#' @param debug Logical, defaults to \code{FALSE}. Turns on all debug checks - this may impact runtimes.
	#' @export
	#' @inheritParams sample_markov_base
	#' @inherit sample_markov_base return
	#' @examples
	#'
	#' sample <- sample_markov(example_two_state_markov())
	#'
	#' sample
	sample_markov <- function(markov_model = NULL,
	type = "rcpp",
	debug = FALSE,
	samples = 1) {
	if (debug) {
	if (!is.list(markov_model)) {
	stop("The markov model must be supplied as a list.
	See SpeedyMarkov::example_two_state_markov for details of the required data format.")
	}
	}
	out <- sample_markov_base(
	transitions = markov_model[["transitions_list"]],
	cohorts = markov_model[["cohorts"]],
	state_costs = markov_model[["state_costs"]],
	intervention_costs = markov_model[["intervention_costs"]],
	qalys = markov_model[["qalys"]],
	samples = samples,
	type = type
	)
	return(out)
	}
	#' Sample a Markov Model Sample using Base R
	#'
	#'
	#' @description This model agnostic function samples a markov model specification using a base R implementation.
	#' See `example_two_state_markov` for an example of the required input. Alternatively use `sample_markov(type = "base")` passing the
	#' model specification function.
	#'
	#' @param transitions A function that generates a list of transition matrices,
	#' see `example_two_state_markov` for an example of setting this up.
	#' @param cohorts A function that generates a list containing the initial state for each intervention,
	#' see `example_two_state_markov` for an example of setting this up.
	#' @param state_costs A function that generates a list of state costs for each intervention,
	#' see `example_two_state_markov` for an example of setting this up.
	#' @param intervention_costs A function that generates a vector of intervention costs, see `example_two_state_markov` for an example of setting this up.
	#' @param qalys A function that generates a list of QALYs for each intervention, see `example_two_state_markov` for an example of setting this up.
	#' @param samples Numeric, defaults to 1. The number of samples to take from the Markov model
	#' @param type A character string specifying the approach to use to sample the model. Currently implemented
	#' approaches are "base" and "rcpp" with "rcpp" as the default.
	#' @return A data.frame of samples of a model encoded in the `SpeedyMarkov` format (see `example_two_state_markov` for details).
	#' @export
	#' @importFrom purrr map transpose
	#' @importFrom tibble tibble
	#' @examples
	#'
	#'
	#' markov_model <- example_two_state_markov()
	#'
	#' sample_markov_base(
	#' transitions = markov_model$transitions_list,
	#' cohorts = markov_model$cohorts,
	#' state_costs = markov_model$state_costs,
	#' intervention_costs = markov_model$intervention_costs,
	#' qalys = markov_model$qalys
	#' )
	#'
	sample_markov_base <- function(transitions = NULL, state_costs = NULL,
	intervention_costs = NULL, cohorts = NULL,
	qalys = NULL, samples = 1, type = "rcpp") {
	#sample baseline transition matrix
	baseline <- transitions[[1]](samples = samples, type = type)
	#sample all interventions depending on baseline
	interventions <- purrr::map(2:length(transitions), ~ transitions[[.]](baseline, type = type))
	#update transitions as a single sample
	transitions[[1]] <- baseline
	transitions[-1] <- interventions
	## Organise as dataframe
	## Unlist interventions
	samples_df <- tibble::tibble(
	sample = unlist(purrr::map(1:samples,
	~ rep(., length(transitions)))),
	intervention = rep(names(transitions), samples),
	transition = purrr::flatten(purrr::transpose(transitions)),
	state_cost = purrr::flatten(state_costs(samples = samples, type = type)),
	intervention_cost = unlist(intervention_costs(samples = samples, type = type)),
	cohort = purrr::flatten(cohorts(samples = samples, type = type)),
	qalys = purrr::flatten(qalys(samples = samples, type = type)))
	return(samples_df)
	}

SpeedyMarkov/R/markov_simulation_pipeline.R		Memory				Time
	#' Markov Sampling and Simulation Pipeline
	#'
	#' @description This functions wraps multiple modular functions and allows an end-to-end cost effectiveness to
	#' be run, excluding the final analysis of the findings. It may also be used in batch mode to run analyses in
	#' parallel.
	#' @param samples Numeric, defaults to 1. The number of markov model samples to use.
	#' @param sample_type A character string specifying the approach to use to sample the model.
	#' Options and defaults inherited from `sample_markov`.
	#' @param sim_type A character string specifying the approach to use to simulate the model.
	#' Options and defaults inherited from `simulate_markov`.
	#' @param discount Numeric, the discount that should be applied to the costs and QALYs. Defaults to 1.035.
	#' @param batches Numeric, defaults to 1. The number of batches to run simulation/sampling in. When set to
	#' values greater than 1 a `batch_fn` must also be supplied. It is likely that the most efficient option will
	#' be to use a batch number that corresponds to the number of cores being utilised.
	#' @param batch_fn Function, defaults to `NULL`. This is the function to be used to parallise across batches. Potential options
	#' include `parallel::mclapply` (not Windows) or `furrr::future_map` (requires the use of `future::plan` outside the function). When
	#' not given the function will default to using no batching.
	#' @param ... Additional options to pass to `batch_fn`. For example this may be the `mc.cores` argument of `parallel::mclapply`.
	#' @return A list containing the model samples and simulations.
	#' @export
	#' @importFrom data.table rbindlist
	#' @importFrom tibble as_tibble
	#' @importFrom dplyr bind_cols
	#' @importFrom purrr transpose map
	#' @inheritParams simulate_markov
	#' @inheritParams sample_markov
	#' @seealso sample_markov simulate_markov
	#' @examples
	#'
	#'
	#' markov_simulation_pipeline(example_two_state_markov(), duration = 10, samples = 2)
	markov_simulation_pipeline <- function(markov_model = NULL, duration = NULL,
	discount = 1.035, samples = 1,
	sample_type = "rcpp",
	sim_type = "armadillo_all",
	batches = 1,
	batch_fn = NULL,
	debug = FALSE,
	...) {
	# Define pipeline ---------------------------------------------------------
	markov_simulation_pipeline_inner <- function(samples = NULL) {
	# Generate samples --------------------------------------------------------
	model_samples <- SpeedyMarkov::sample_markov(markov_model,
	type = sample_type,
	debug = debug,
	samples = samples)
	# Allocate simulation matrix ----------------------------------------------
	## This optional step gives a small speed boost by making preallocation occur once
	## rather than for every model simulation
	## Pull out a template transition matrix
	template_transition <- model_samples$transition[[1]]
	## Set up simulation preallocation
	sim_storage <- matrix(NA, nrow = duration, ncol = nrow(template_transition))
	colnames(sim_storage) <- colnames(template_transition )
	# Calculate discounting ---------------------------------------------------
	discounting <- SpeedyMarkov::calc_discounting(discount, duration)
	# Simulate model from samples ---------------------------------------------
	## Map samples to a list for efficiency
	samples_list <- purrr::transpose(model_samples)
	## Simulate over samples and interventions
	results <- purrr::map(samples_list, function(sample) {
	SpeedyMarkov::simulate_markov(
	markov_sample = sample,
	duration = duration,
	discounting = discounting,
	type = sim_type,
	sim = sim_storage,
	input_is_list = TRUE,
	debug = debug)
	})
	## Parallel data frame binding for results from data.table
	results <- data.table::rbindlist(results)
	## Combine samples and simulation results
	combined <- dplyr::bind_cols(model_samples, results)
	return(combined)
	}
	# Check batch instructions ------------------------------------------------
	if (is.null(batch_fn) & batches > 1) {
	message("No batch function has been supplied so falling back to not using batching (batches = 1)")
	batches <- 1
	}
	if (batches > samples) {
	stop("The number of batches should not be greater than the number of samples")
	}
	# Simulate/Sample single/batch ----------------------------------------
	if (batches == 1) {
	##Single batch (i.e no parallisation)
	combined <- markov_simulation_pipeline_inner(samples = samples)
	}else{
	## Find samples that won't divide evenly into batches
	div_diff_samples <- samples %% batches
	## Divide divisible samples into batch sizes
	diff_samples <- (samples - div_diff_samples) / batches
	## Make batch sample vector
	batch_samples <- rep(diff_samples, batches)
	## Add non-divisible batches to final batch
	batch_samples[batches] <- batch_samples[batches] + div_diff_samples
	## Run batch simulations
	combined <- batch_fn(batch_samples, function(batch_sample){
	markov_simulation_pipeline_inner(samples = batch_sample)
	}, ...)
	## Combine batch simulations
	combined <- data.table::rbindlist(combined)
	}
	## Convert to tibble for ease of interaction once returned
	combined <- tibble::as_tibble(combined)
	return(combined)
	}

SpeedyMarkov/R/markov_ce_pipeline.R		Memory				Time
	#' Markov Sampling, Simulation, and Cost Effectiveness Analysis Pipeline
	#'
	#' @description This functions wraps multiple modular functions and allows an end-to-end cost effectiveness to
	#' be run, including final analysis of the findings.
	#' @return A list containing the model samples and simulations and cost effectiveness summary measures.
	#' @export
	#' @inheritParams markov_simulation_pipeline
	#' @inheritParams analyse_ce
	#' @seealso markov_simulation_pipeline analyse_ce
	#' @examples
	#'
	#' markov_ce_pipeline(example_two_state_markov(), duration = 10, samples = 5)
	#'
	markov_ce_pipeline <- function(markov_model = NULL, duration = NULL,
	discount = 1.035, samples = 1, baseline = 1,
	willingness_to_pay_threshold = 20000,
	sample_type = "rcpp", sim_type = "armadillo_all", debug = FALSE,
	batches = 1, batch_fn = NULL, ...) {
	# Sample and simulation markov --------------------------------------------
	simulations <- SpeedyMarkov::markov_simulation_pipeline(markov_model = markov_model,
	duration = duration,
	discount = discount,
	samples = samples,
	sample_type = sample_type,
	sim_type = sim_type,
	batches = batches,
	batch_fn = batch_fn,
	debug = debug, ...)
	# Analyse model -----------------------------------------------------------
	sum <- SpeedyMarkov::analyse_ce(simulations, baseline = baseline,
	willingness_to_pay_threshold = willingness_to_pay_threshold)
	return(sum)
	}

SpeedyMarkov/R/RcppExports.R		Memory				Time
	# Generated by using Rcpp::compileAttributes() -> do not edit by hand
	# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
	#' @title An inner Markov loop implemented using RcppArmadillo
	#' @inherit markov_loop
	#' @export
	#' @useDynLib SpeedyMarkov, .registration=TRUE
	#' @examples
	#'
	#' transition <- matrix(rnorm(4), 2, 2)
	#' sim <- matrix(NA, 100, 2)
	#' cohort <- c(1, 0)
	#' duration <- 100
	#'
	#' # Reference R implementation
	#' sim_r <- markov_loop(sim, cohort, transition, duration)
	#'
	#' # RcppArmadillo implementation
	#' sim_rcppArma <- ArmaMarkovLoop(sim, cohort, transition, duration)
	#'
	#' # Check results are within tolerances
	#' all.equal(sim_r, sim_rcppArma)
	#'
	#' # Benchmark
	#' library(microbenchmark)
	#' microbenchmark(markov_loop(sim, cohort, transition, duration),
	#' ArmaMarkovLoop(sim, cohort, transition, duration),
	#' times = 1000)
	ArmaMarkovLoop <- function(sim, cohort, transition, duration) {
	.Call(`_SpeedyMarkov_ArmaMarkovLoop`, sim, cohort, transition, duration)
	}
	#' @title Simulate a Markov Model Sample using RcppArmadillo
	#'
	#' @description This model agnostic function runs a single markov model for the specified duration using a Armadillo implementation.
	#' See `example_two_state_markov` for an example of the required input. Alternatively use `sample_markov(type = "base")`
	#' and the output from `sample_markov`.
	#' @inherit simulate_markov_base
	#' @export
	#' @useDynLib SpeedyMarkov, .registration=TRUE
	#' @examples
	#'
	#' ## Sample model
	#' markov_sample <- sample_markov(example_two_state_markov())
	#'
	#' ## Simulate using R
	#' sim_r <- simulate_markov_base(
	#' ## Specify the storage simulation matrix to maintain consistency
	#' ##here (but not needed for the base implementation).
	#' sim = matrix(NA, nrow = 10, ncol = nrow(markov_sample$transition[[1]])),
	#' transition = markov_sample$transition[[1]],
	#' cohort = markov_sample$cohort[[1]],
	#' state_cost = markov_sample$state_cost[[1]],
	#' intervention_cost = markov_sample$intervention_cost[[1]],
	#' qalys = markov_sample$qalys[[1]],
	#' duration = 10,
	#' discounting = SpeedyMarkov::calc_discounting(1.035, 10),
	#' markov_loop_fn = SpeedyMarkov::markov_loop
	#')
	#'
	#' # RcppArmadillo implementation
	#' sim_rcppArma <- ArmaSimulateMarkov(
	#' sim = matrix(NA, nrow = 10, ncol = nrow(markov_sample$transition[[1]])),
	#' transition = markov_sample$transition[[1]],
	#' cohort = markov_sample$cohort[[1]],
	#' state_cost = markov_sample$state_cost[[1]],
	#' intervention_cost = markov_sample$intervention_cost[[1]],
	#' qalys = markov_sample$qalys[[1]],
	#' duration = 10,
	#' discounting = SpeedyMarkov::calc_discounting(1.035, 10)
	#')
	#'
	#' # Check results are within tolerances
	#' all.equal(sim_r, sim_rcppArma)
	#'
	#' # Benchmark
	#' library(microbenchmark)
	#' microbenchmark(simulate_markov_base(
	#' sim = matrix(NA, nrow = 100, ncol = nrow(markov_sample$transition[[1]])),
	#' transition = markov_sample$transition[[1]],
	#' cohort = markov_sample$cohort[[1]],
	#' state_cost = markov_sample$state_cost[[1]],
	#' intervention_cost = markov_sample$intervention_cost[[1]],
	#' qalys = markov_sample$qalys[[1]],
	#' duration = 100,
	#' discounting = SpeedyMarkov::calc_discounting(1.035, 100),
	#' markov_loop_fn = SpeedyMarkov::markov_loop),
	#' ArmaSimulateMarkov(
	#' sim = matrix(NA, nrow = 100, ncol = nrow(markov_sample$transition[[1]])),
	#' transition = markov_sample$transition[[1]],
	#' cohort = markov_sample$cohort[[1]],
	#' state_cost = markov_sample$state_cost[[1]],
	#' intervention_cost = markov_sample$intervention_cost[[1]],
	#' qalys = markov_sample$qalys[[1]],
	#' duration = 100,
	#' discounting = SpeedyMarkov::calc_discounting(1.035, 100)),
	#' times = 1000)
	ArmaSimulateMarkov <- function(sim, cohort, transition, duration, state_cost, discounting, qalys, intervention_cost) {
	.Call(`_SpeedyMarkov_ArmaSimulateMarkov`, sim, cohort, transition, duration, state_cost, discounting, qalys, intervention_cost)
	}
	#' @title Arrange Vectorised Matrix Samples using Rcpp
	#'
	#' @description A convenience function used to arrange vectorised matrix samples into the correct
	#' matrix format for several functions used to specify Markov model. See `example_two_state_markov`
	#' for an example use case. Implemented using Rcpp.
	#' @inherit matrix_arrange_inner
	#' @export
	#' @useDynLib SpeedyMarkov, .registration=TRUE
	#' @importFrom Rcpp evalCpp
	#' @examples
	#'
	#'
	#' matrix_samples <- list(VGAM::rdiric(1:5, c(88, 12)),
	#' VGAM::rdiric(1:5, c(8, 92)))
	#'
	#' # R implementation
	#' samples_r <- matrix_arrange_inner(matrix_samples)
	#'
	#'
	#' # Rcpp implementation
	#' samples_rcpp <- MatrixArrange(matrix_samples)
	#'
	#' all.equal(samples_r, samples_rcpp)
	#'
	#' # Benchmark
	#' library(microbenchmark)
	#' microbenchmark(matrix_arrange_inner(matrix_samples),
	#' MatrixArrange(matrix_samples),
	#' times = 1000)
	MatrixArrange <- function(samples) {
	.Call(`_SpeedyMarkov_MatrixArrange`, samples)
	}

SpeedyMarkov/R/matrix_arrange.R		Memory				Time
	#' Arrange Vectorised Matrix Samples
	#'
	#' @description A convenience function used to arrange vectorised matrix samples into the correct
	#' matrix format for several functions used to specify Markov model. See `example_two_state_markov`
	#' for an example use case. Both a R and C++ implementation are available.
	#' @param type A character string specifying the approach to use. Currently implemented
	#' approaches are "base", and "rcpp" with "rcpp" as the default.
	#' @export
	#' @inherit matrix_arrange_inner
	#' @examples
	#'
	#' matrix_samples <- list(VGAM::rdiric(1:5, c(88, 12)),
	#' VGAM::rdiric(1:5, c(8, 92)))
	#'
	#' # Default usage
	#' matrix_arrange(matrix_samples)
	#' # R implementation
	#' matrix_arrange(matrix_samples, type = "base")
	#' # Rcpp implementation
	#' matrix_arrange(matrix_samples, type = "rcpp")
	#'
	matrix_arrange <- function(samples, type = "rcpp") {
	#Send findings to inner functions
	if (type == "base") {
	tmp <- SpeedyMarkov::matrix_arrange_inner(samples)
	}else if (type == "rcpp"){
	tmp <- SpeedyMarkov::MatrixArrange(samples)
	}
	return(tmp)
	}
	#' Arrange Vectorised Matrix Samples using R
	#'
	#' @description A convenience function used to arrange vectorised matrix samples into the correct
	#' matrix format for several functions used to specify Markov model. See `example_two_state_markov`
	#' for an example use case. Implemented using R.
	#' @param samples A list of vectorised matrix samples
	#'
	#' @return A list of matrices with each matrix representing a single sample.
	#' @importFrom purrr map
	#' @export
	#'
	#' @examples
	#'
	#' matrix_samples <- list(VGAM::rdiric(1:5, c(88, 12)),
	#' VGAM::rdiric(1:5, c(8, 92)))
	#'
	#' matrix_arrange_inner(matrix_samples)
	matrix_arrange_inner <- function(samples) {
	# Map transitions into matrices
	dim <- length(samples)
	out <- matrix(NA, dim, dim)
	tmp <- purrr::map(1:nrow(samples[[1]]), ~ {
	for (i in 1:dim) {
	out[i, ] <- samples[[i]][., ]
	}
	return(out)
	})
	return(tmp)
	}

SpeedyMarkov/R/vector_arrange.R		Memory				Time
	#' Arrange Vectorised Vector Samples
	#'
	#' @description A convenience function used to arrange vectorised samples into the correct
	#' vector format for several functions used to specify Markov model. See `example_two_state_markov`
	#' for an example use case.
	#' @param samples A list of vectorised samples
	#'
	#' @return A list of vectors with each vector representing a single sample.
	#' @importFrom purrr map
	#' @export
	#'
	#' @examples
	#'
	#' vector_samples <- list(rnorm(10, 1, 2), rnorm(10, 10, 2))
	#'
	#' vector_arrange(vector_samples)
	vector_arrange <- function(samples) {
	## Sample
	tmp <- do.call(cbind, samples)
	## Split into vectors and name
	out <- split(tmp, 1:nrow(tmp))
	names(out) <- NULL
	return(out)
	}

SpeedyMarkov/R/simulate_markov.R		Memory				Time
	#' Simulate a Markov Model Sample
	#'
	#'
	#'@description This model agnostic function runs a single markov model for the specified duration. It wraps multiple approaches
	#'that may offer various advantages and disadvantages.
	#'
	#' @param markov_sample A single row dataframe or a list with no default. See `sample_markov`
	#' for the correct data format.
	#' @param type A character string specifying the approach to use to simulate the model. Currently implemented
	#' approaches are "base", "armadillo_inner" and "armadillo_all with "armadillo_inner" as the default.
	#' "armadillo_all" is likely to be generally faster but has a slightly reduced feature set.
	#' @param debug Logical, defaults to \code{FALSE}. Turns on all debug checks - this may impact runtimes.
	#' @param input_is_list Logical defaults to NULL. What type of input is `markov_sample`? A list or a dataframe.
	#' If not given then the input type will be checked and converted to a list as required.
	#' @importFrom purrr transpose
	#' @export
	#' @inherit simulate_markov_base
	#' @examples
	#'
	#' sample <- sample_markov(example_two_state_markov())
	#'
	#' simulate_markov(sample[1, ], duration = 10, discounting = SpeedyMarkov::calc_discounting(1.035, 10))
	#'
	simulate_markov <- function(markov_sample = NULL,
	duration = NULL,
	discounting = NULL,
	type = "armadillo_all",
	debug = FALSE,
	input_is_list = NULL,
	sim = NULL) {
	if(debug) {
	if (!is.list(markov_sample)) {
	stop("The markov sample must be supplied as a list or dataframe (1 row).
	See SpeedyMarkov::sample_markov for details of the required data format.")
	}
	if (tibble::is.tibble(markov_sample) && nrow(markov_sample) > 1) {
	stop("Only 1 sample may be simulated at a time.")
	}
	if (type == "armadillo_all" & is.null(sim)) {
	stop("This implementation requires the sim matrix to be prespecified.
	See the docs for details.")
	}
	}
	## If the input type is not given check for a dataframe and convert as required
	if (is.null(input_is_list) | isFALSE(input_is_list)) {
	if (is.data.frame(markov_sample)) {
	# Flip input format into a nested list
	markov_sample <- purrr::transpose(markov_sample)
	# Extract the first object from the list
	markov_sample <- markov_sample[[1]]
	}
	}
	## Assign transition
	transition <- markov_sample[["transition"]]
	## Preallocate
	if (is.null(sim)) {
	sim <- matrix(NA, nrow = duration, ncol = nrow(transition))
	}
	if (type == "base") {
	out <- SpeedyMarkov::simulate_markov_base(
	transition = transition,
	cohort = markov_sample[["cohort"]],
	state_cost = markov_sample[["state_cost"]],
	intervention_cost = markov_sample[["intervention_cost"]] ,
	qalys = markov_sample[["qalys"]],
	duration = duration,
	discounting = discounting,
	sim = sim,
	markov_loop_fn = SpeedyMarkov::markov_loop
	)
	}else if (type == "armadillo_inner") {
	out <- SpeedyMarkov::simulate_markov_base(
	transition = transition,
	cohort = markov_sample[["cohort"]],
	state_cost = markov_sample[["state_cost"]],
	intervention_cost = markov_sample[["intervention_cost"]] ,
	qalys = markov_sample[["qalys"]],
	duration = duration,
	discounting = discounting,
	sim = sim,
	markov_loop_fn = SpeedyMarkov::ArmaMarkovLoop
	)
	}else if (type == "armadillo_all") {
	out <- SpeedyMarkov::ArmaSimulateMarkov(
	transition = transition,
	cohort = markov_sample[["cohort"]],
	state_cost = markov_sample[["state_cost"]],
	intervention_cost = markov_sample[["intervention_cost"]] ,
	qalys = markov_sample[["qalys"]],
	duration = duration,
	discounting = discounting,
	sim = sim
	)
	}
	return(out)
	}
	#' Simulate a Markov Model Sample using Base R
	#'
	#'
	#' @description This model agnostic function runs a single markov model for the specified duration using a base R implementation.
	#' See `example_two_state_markov` for an example of the required input. Alternatively use `sample_markov(type = "base")` and
	#' the output from `sample_markov`.
	#' @param state_cost A list of state costs for each intervention,
	#' see `example_two_state_markov` for an example of setting this up.
	#' @param intervention_cost A vector of intervention costs, see `example_two_state_markov` for an example of setting this up.
	#' @param qalys A list of QALYs for each intervention, see `example_two_state_markov` for an example of setting this up.
	#' @param discounting Numeric vector, the discount that should be applied to the costs and QALYs for each time period.
	#' This must be the same length as `duration`.
	#' @param sim Matrix with the same number of rows as the duration of the model and the same number of columns as the number of
	#' states in the model. Used to store model simulations.
	#' @param markov_loop_fn A function, defaults to ] \code{NULL}. The function to use to solve the inner markov loops. Built in examples
	#' are `markov_loop` (using `R`) and `ArmaMarkovLoop` (using `RcppArmadillo`)
	#' @return A list containing total costs and total QALYs as matrices across states
	#' and the duration of the model
	#' @inheritParams markov_loop
	#' @export
	#'
	#' @examples
	#'
	#' markov_sample <- sample_markov(example_two_state_markov())
	#'
	#' simulate_markov_base(
	#' sim = matrix(NA, nrow = 10, ncol = nrow(markov_sample$transition[[1]])),
	#' transition = markov_sample$transition[[1]],
	#' cohort = markov_sample$cohort[[1]],
	#' state_cost = markov_sample$state_cost[[1]],
	#' intervention_cost = markov_sample$intervention_cost[[1]],
	#' qalys = markov_sample$qalys[[1]],
	#' duration = 10,
	#' discounting = SpeedyMarkov::calc_discounting(1.035, 10),
	#' markov_loop_fn = SpeedyMarkov::markov_loop
	#' )
	#'
	#'
	simulate_markov_base <- function(transition = NULL, cohort = NULL, state_cost = NULL,
	intervention_cost = NULL, qalys = NULL, duration = NULL,
	discounting = NULL, sim = NULL,
	markov_loop_fn = NULL) {
	# Simulate model over the loop
	sim <- markov_loop_fn(sim, cohort, transition, duration)
	##Total costs per cycle
	total_costs_cycle <- (sim %*% state_cost) * discounting
	##Total QALYs per cycle
	discounted_qalys <- (sim %*% qalys) * discounting
	total_qalys <- sum(discounted_qalys)
	## Overall costs
	total_costs <- sum(total_costs_cycle) + intervention_cost
	out <- list(total_costs = total_costs, total_qalys = total_qalys)
	return(out)
	}

SpeedyMarkov/R/markov_loop.R		Memory				Time
	#' Inner Markov Loop in base R
	#'
	#' @param sim Matrix with the same number of rows as the duration of the model and the same number of columns as the number of
	#' states in the model.
	#' @param cohort A vector containing the initial state.
	#' @param transition A transition matrix, see `example_two_state_markov` for an example of setting this up.
	#' @param duration Numeric, how many long to run the model for.
	#'
	#' @return A matrix of the same size as the inputted sim matrix.
	#' @export
	#'
	#' @examples
	#'
	#'
	#'transition <- matrix(rnorm(4), 2, 2)
	#'sim <- matrix(NA, 10, 2)
	#'cohort <- c(1, 0)
	#'duration <- 10
	#'
	#'markov_loop(sim, cohort, transition, duration)
	markov_loop <- function(sim = NULL, cohort, transition = NULL, duration = NULL) {
	## Assign initial pop
	sim[1, ] <- cohort
	##Loop over the rest of the model
	for (i in 2:duration) {
	sim[i, ] <- sim[i - 1, ] %*% transition
	}
	return(sim)
	}

SpeedyMarkov/R/analyse_ce.R		Memory				Time
	#' Analyse the Cost Effectiveness of Interventions
	#'
	#'
	#'@description This function produces cost effectiveness summary measures using the output of `markov_simulation_pipeline`
	#' or similar data structures. At least two interventions must be present.
	#' @param markov_simulations A dataframe of markov samples and simulations as produced by `markov_simulation_pipeline`. At least two
	#' interventions must be present.
	#' @param baseline Numeric, the intervention to consider as the baseline for pairwise comparisons.
	#' @param willingness_to_pay_threshold Numeric, defaulting to 20,000. This is the threshold at which an intervention
	#' may be considered cost effective in the UK.
	#' @param type A character string specifying the approach to use to simulate the model. Currently implemented
	#' approaches are "base" with "base" as the default.
	#' @return A list of dataframes including: Cost effectiveness measures for each sample, and summarised cost effectiveness
	#' measures across samples.
	#' @export
	#' @importFrom purrr map
	#' @importFrom data.table as.data.table .N ":="
	#' @importFrom stats sd
	#' @examples
	#'
	#' sims <- markov_simulation_pipeline(example_two_state_markov(),
	#' duration = 10, samples = 10)
	#'
	#' analyse_ce(sims)
	#'
	analyse_ce <- function(markov_simulations = NULL,
	baseline = 1,
	willingness_to_pay_threshold = 20000,
	type = "base") {
	## NULL out variables to deal with package notes
	total_costs <- NULL; total_qalys <- NULL; incremental_qalys <- NULL; incremental_costs <- NULL;
	intervention <- NULL; incremental_net_benefit <- NULL; mean_costs <- NULL; mean_qalys <- NULL;
	total_costs <- NULL; total_costs <- NULL; icer <- NULL; . <- NULL;
	## Convert to data.table
	markov_simulations <- data.table::as.data.table(markov_simulations)
	## Calculate incremental costs and qalys
	incremental_sims <- markov_simulations[, c("incremental_costs", "incremental_qalys") :=
	list(total_costs - total_costs[baseline],
	total_qalys - total_qalys[baseline]),
	by = "sample"][,
	incremental_net_benefit := willingness_to_pay_threshold * incremental_qalys - incremental_costs
	,]
	## Summarise costs
	summarised_sims <- incremental_sims[,
	.( mean_costs = mean(total_costs),
	sd_costs = stats::sd(total_costs),
	mean_qalys = mean(total_qalys),
	sd_qlays = stats::sd(total_qalys),
	mean_incremental_qalys = mean(incremental_qalys),
	sd_incremental_qlays = stats::sd(incremental_qalys),
	mean_incremental_costs = mean(incremental_costs),
	sd_incremental_costs = stats::sd(incremental_costs),
	mean_incremental_net_benefit = mean(incremental_net_benefit),
	sd_incremental_net_benefit = stats::sd(incremental_net_benefit),
	probability_cost_effective = sum(incremental_net_benefit > 0) / .N
	),
	by = "intervention"
	][,
	icer := mean_costs / mean_qalys
	,]
	output <- list(incremental_sims, summarised_sims)
	## Convert output to tibble for nice presentation
	output <- purrr::map(output, as_tibble)
	names(output) <- c("simulations_with_ce", "summarised_ce")
	return(output)
	}

Sample Interval: 10ms

30740ms

☒ Split horizontally

☒ Hide internal function calls

☐ Hide lines of code with zero time

☐ Hide memory results