#**************************************************************************
#*
#* Original work Copyright (C) 2016 Antoine Pierucci
#* Modified work Copyright (C) 2017 Matt Wiener
#* Modified work Copyright (C) 2016 Jordan Amdahl
#*
#* This program is free software: you can redistribute it and/or modify
#* it under the terms of the GNU General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or
#* (at your option) any later version.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#**************************************************************************
#' Run Markov Model
#'
#' Runs one or more strategy. When more than one strategy is
#' provided, all strategies should have the same states and
#' state value names.
#'
#' In order to compute comparisons strategies must be
#' similar (same states and state value names). Thus
#' strategies can only differ through transition matrix cell
#' values and values attached to states (but not state value
#' names).
#'
#' The initial number of individuals in each state and the
#' number of cycle will be the same for all strategies
#'
#' `state_time_limit` can be specified in 3 different ways:
#' 1. As a single value: the limit is applied to all states
#' in all strategies. 2. As a named vector (where names are
#' state names): the limits are applied to the given state
#' names, for all strategies. 3. As a named list of named
#' vectors: the limits are applied to the given state names
#' for the given strategies.
#'
#' @param ... One or more `uneval_model` object.
#' @param parameters Optional. An object generated by
#' [define_parameters()].
#' @param init numeric vector or result of [define_init()],
#' same length as number of states. Number of individuals
#' in each state at the beginning.
#' @param cycles positive integer. Number of Markov Cycles
#' to compute.
#' @param cost Names or expression to compute cost on the
#' cost-effectiveness plane.
#' @param effect Names or expression to compute effect on
#' the cost-effectiveness plane.
#' @param method Counting method.
#' @param uneval_strategy_list List of models, only used by
#' [run_model_()] to avoid using `...`.
#' @param state_time_limit Optional expansion limit for
#' `state_time`, see details.
#' @param central_strategy character. The name of the
#' strategy at the center of the cost-effectiveness plane,
#' for readability.
#' @param inflow numeric vector or result of
#' [define_inflow()], similar to `init`. Number of new
#' individuals in each state per cycle.
#' @param aux_params Optional. an object generated by [define_parameters()]
#' representing object parameters.
#' @param parallel Optional. Whether to use parallel evaluation of strategies.
#'
#' @return A list of evaluated models with computed values.
#' @export
#'
#' @example inst/examples/example_run_model.R
#'
run_model <- function(...,
parameters = define_parameters(),
init = c(1000L, rep(0L, get_state_number(get_states(list(...)[[1]])) - 1)),
cycles = 1,
method = "life-table",
cost = NULL, effect = NULL,
state_time_limit = NULL,
central_strategy = NULL,
inflow = rep(0L, get_state_number(get_states(list(...)[[1]]))),
aux_params = NULL,
parallel = F,
cores = 1,
disc_method = 'start',
create_progress_reporter = create_null_prog_reporter,
progress_reporter = create_progress_reporter(),
state_groups = NULL,
individual_level = F) {
uneval_strategy_list <- list(...)
init <- check_init(init, get_state_names(uneval_strategy_list[[1]]))
inflow <- check_inflow(inflow, get_state_names(uneval_strategy_list[[1]]))
run_model_(
uneval_strategy_list = uneval_strategy_list,
parameters = parameters,
init = init,
cycles = cycles,
method = method,
cost = lazy_(substitute(cost), env = parent.frame()),
effect = lazy_(substitute(effect), env = parent.frame()),
state_time_limit = state_time_limit,
central_strategy = central_strategy,
inflow = inflow,
aux_params = aux_params,
parallel = parallel,
cores = cores,
disc_method = disc_method,
create_progress_reporter = create_progress_reporter,
progress_reporter = progress_reporter,
state_groups = state_groups,
individual_level = individual_level
)
}
#' @export
#' @rdname run_model
run_model_ <- function(uneval_strategy_list,
parameters,
init,
cycles,
method,
cost, effect,
state_time_limit,
central_strategy,
inflow,
aux_params = NULL,
parallel = F,
cores = 1,
disc_method = 'start',
create_progress_reporter = create_null_prog_reporter,
progress_reporter = create_progress_reporter(),
state_groups = NULL,
individual_level = F) {
if (length(uneval_strategy_list) == 0) {
stop("At least 1 strategy is needed.")
}
if (! is.wholenumber(cycles)) {
stop("'cycles' must be a whole number.")
}
if (! all(unlist(lapply(
uneval_strategy_list,
function(x) "uneval_model" %in% class(x))))) {
.x <- names(uneval_strategy_list[! unlist(lapply(
uneval_strategy_list,
function(x) "uneval_model" %in% class(x)))])
stop(sprintf(
"Incorrect model object%s: %s.",
plur(length(.x)),
paste(.x, collapse = ", ")
))
}
list_ce <- list(
.cost = cost,
.effect = effect
)
ce <- c(
lazy_dots(),
list_ce
)
strategy_names <- names(uneval_strategy_list)
if (is.null(strategy_names)) {
message("No named model -> generating names.")
strategy_names <- as.character(utils::as.roman(seq_along(uneval_strategy_list)))
names(uneval_strategy_list) <- strategy_names
}
if (any(strategy_names == "")) {
warning("Not all models are named -> generating names.")
strategy_names <- as.character(utils::as.roman(seq_along(uneval_strategy_list)))
names(uneval_strategy_list) <- strategy_names
}
if (! list_all_same(lapply(uneval_strategy_list,
function(x) sort(get_state_names(x))))) {
stop("State names differ between models.")
}
if (! list_all_same(lapply(uneval_strategy_list,
function(x) sort(get_state_value_names(x))))) {
stop("State value names differ between models.")
}
# Check state groups
check_state_groups(state_groups, get_state_names(uneval_strategy_list[[1]]))
state_time_limit <- complete_stl(
state_time_limit,
state_names = get_state_names(uneval_strategy_list[[1]]),
strategy_names = names(uneval_strategy_list),
cycles = cycles,
state_groups = state_groups
)
#eval_strategy_list <- list()
if (cores > 1 && parallel && .Platform$OS.type == "unix") {
eval_strategy_list <- mclapply(seq_len(length(uneval_strategy_list)), function(i) {
try(eval_strategy(
strategy = uneval_strategy_list[[i]],
parameters = parameters,
init = init,
cycles = cycles,
method = method,
expand_limit = state_time_limit[[i]],
inflow = inflow,
strategy_name = names(uneval_strategy_list)[i],
aux_params = aux_params,
disc_method = disc_method,
progress_reporter = create_progress_reporter(),
state_groups = state_groups,
individual_level = individual_level
))
}, mc.cores = cores)
plyr::l_ply(
eval_strategy_list,
function(x) {
if ("try-error" %in% class(x)) {
err_string <- as.character(attr(x,'condition'))
stop(clean_err_msg(err_string), call. = F)
}
}
)
} else {
eval_strategy_list <- lapply(seq_len(length(uneval_strategy_list)), function(i) {
eval_strategy(
strategy = uneval_strategy_list[[i]],
parameters = parameters,
init = init,
cycles = cycles,
method = method,
expand_limit = state_time_limit[[i]],
inflow = inflow,
strategy_name = names(uneval_strategy_list)[i],
aux_params = aux_params,
disc_method = disc_method,
progress_reporter = progress_reporter,
state_groups = state_groups,
individual_level = individual_level
)
})
}
names(eval_strategy_list) <- names(uneval_strategy_list)
list_res <- lapply(eval_strategy_list, get_total_state_values)
for (n in strategy_names) {
list_res[[n]]$.strategy_names <- n
}
res <-
bind_rows(list_res) %>%
mutate(!!!lazy_eval(ce, data = .))
root_strategy <- get_root_strategy(res)
noncomparable_strategy <- get_noncomparable_strategy(res)
if (is.null(central_strategy)) {
central_strategy <- get_central_strategy(res)
} else {
stopifnot(
length(central_strategy) == 1,
central_strategy %in% strategy_names
)
}
structure(
list(
run_model = res,
eval_strategy_list = eval_strategy_list,
uneval_strategy_list = uneval_strategy_list,
parameters = parameters,
aux_params = aux_params,
init = init,
inflow = inflow,
cycles = cycles,
method = method,
ce = ce,
root_strategy = root_strategy,
central_strategy = central_strategy,
noncomparable_strategy = noncomparable_strategy,
state_time_limit = state_time_limit,
frontier = if (! is.null(root_strategy)) get_frontier(res),
disc_method = disc_method,
create_progress_reporter = create_progress_reporter,
progress_reporter = progress_reporter,
state_groups = state_groups,
individual_level = individual_level,
cores = cores
),
class = c("run_model", class(res))
)
}
get_strategy_names <- function(x) {
get_model_results(x)$.strategy_names
}
get_strategy_count <- function(x) {
nrow(get_model_results(x))
}
get_state_value_names.run_model <- function(x) {
get_state_value_names(x$uneval_strategy_list[[1]])
}
get_total_state_values <- function(x) {
# faster than as.data.frame or as_data_frame
res <- as.list(colSums((x$values)[- 1]))
class(res) <- "data.frame"
attr(res, "row.names") <- c(NA, -1)
res$.n_indiv <- get_n_indiv(x)
res
}
get_inflow <- function(x) {
x$inflow
}
get_root_strategy <- function(x, ...) {
UseMethod("get_root_strategy")
}
get_root_strategy.default <- function(x, ...) {
if (! all(c(".cost", ".effect") %in% names(x))) {
warning("No cost and/or effect defined, cannot find root strategy.")
return(invisible(NULL))
}
(x %>%
arrange(.cost, desc(.effect)))$.strategy_names[1]
}
get_root_strategy.run_model <- function(x, ...) {
x$root_strategy
}
get_noncomparable_strategy <- function(x, ...) {
UseMethod("get_noncomparable_strategy")
}
get_noncomparable_strategy.default <- function(x, ...) {
if (! ".effect" %in% names(x)) {
warning("No effect defined, cannot find noncomparable strategy.")
return(invisible(NULL))
}
(x %>%
arrange(.effect) %>%
slice(1))$.strategy_names
}
get_noncomparable_strategy.run_model <- function(x, ...) {
x$noncomparable_strategy
}
get_central_strategy <- function(x, ...) {
UseMethod("get_central_strategy")
}
get_central_strategy.default <- function(x, ...) {
get_root_strategy(x)
}
get_central_strategy.run_model <- function(x, ...) {
x$central_strategy
}
get_effect <- function(x) {
get_model_results(x)$.effect
}
get_n_indiv.default <- function(x) {
get_model_results(x)$.n_indiv
}
get_n_indiv.combined_model <- function(x) {
get_n_indiv.default(x)
}
#' Get Strategy Values
#'
#' Given a result from [run_model()], return
#' cost and effect values for a specific strategy.
#'
#' @param x Result from [run_model()].
#' @param ... further arguments passed to or from other
#' methods.
#'
#' @return A data frame of values per state.
#' @export
get_values <- function(x, ...) {
UseMethod("get_values")
}
#' @rdname get_values
#' @export
get_values.run_model <- function(x, ...) {
res <- do.call(
bind_rows,
lapply(
get_strategy_names(x),
function(.n) {
get_values(x$eval_strategy_list[[.n]]) %>%
mutate(.strategy_names = .n)
}
)
)
reshape_long(
data = res,
key_col = "value_names",
value_col = "value",
gather_cols = names(res)[! names(res) %in%
c("markov_cycle",
".strategy_names")]
)
}
#' @rdname get_values
#' @export
get_values.eval_strategy <- function(x, ...) {
x$values
}
#' @rdname get_values
#' @export
get_values.list <- function(x, ...) {
x$values
}
#' Get State Membership Counts
#'
#' Given a result from [run_model()], return
#' state membership counts for a specific strategy.
#'
#' @param x Result from [run_model()].
#' @param ... further arguments passed to or from other
#' methods.
#'
#' @return A data frame of counts per state.
#' @export
get_counts <- function(x, ...) {
UseMethod("get_counts")
}
#' @rdname get_counts
#' @export
get_counts.run_model <- function(x, ...) {
res <- do.call(
bind_rows,
lapply(
get_strategy_names(x),
function(.n) {
get_counts(x$eval_strategy_list[[.n]]) %>%
mutate(
.strategy_names = .n,
markov_cycle = row_number())
}
)
)
reshape_long(
data = res,
key_col = "state_names",
value_col = "count",
gather_cols = names(res)[! names(res) %in%
c("markov_cycle",
".strategy_names")]
)
}
#' @rdname get_counts
#' @export
get_counts.eval_strategy <- function(x, ...) {
x$counts
}
#' @rdname get_counts
#' @export
get_counts.list <- function(x, ...) {
x$counts
}
get_uneval_init <- function(x) {
UseMethod("get_uneval_init")
}
get_uneval_init.default <- function(x) {
x$init
}
get_uneval_inflow <- function(x) {
UseMethod("get_uneval_inflow")
}
get_uneval_inflow.default <- function(x) {
x$inflow
}
get_ce <- function(x) {
x$ce
}
get_ce_cost <- function(x) {
get_ce(x)[[1]]
}
get_ce_effect <- function(x) {
get_ce(x)[[2]]
}
get_parameters <- function(x) {
x$parameters
}
get_cycles <- function(x) {
UseMethod("get_cycles")
}
get_cycles.run_model <- function(x) {
x$cycles
}
get_method <- function(x) {
UseMethod("get_method")
}
get_method.run_model <- function(x) {
x$method
}
get_state_names.run_model <- function(x, ...) {
get_state_names(get_states(x$uneval_strategy_list[[1]]))
}
get_expand_limit <- function(x, strategy) {
strategy <- check_strategy_index(x, strategy)
get_eval_strategy_list(x)[[strategy]]$expand_limit
}
get_eval_strategy_list <- function(x) {
x$eval_strategy_list
}
get_parameter_values <- function(x, ...) {
UseMethod("get_parameter_values")
}
get_parameter_values.updated_model <- function(x, ...) {
get_parameter_values(get_model(x), ...)
}
get_parameter_values.run_model <- function(x, parameter_names,
cycles = rep(1, length(parameter_names)),
strategy = 1) {
strategy <- check_strategy_index(x, strategy)
stopifnot(
cycles > 0,
cycles <= get_cycles(x),
all(parameter_names %in% get_parameter_names(x))
)
stats::setNames(
lapply(
seq_along(parameter_names),
function(i) {
as.vector(unlist(
get_eval_strategy_list(x)[[strategy]]$parameters[cycles[i], parameter_names[i]]
))
}),
parameter_names)
}
get_uneval_strategy_list <- function(x) {
x$uneval_strategy_list
}
get_state_time_limit <- function(x) {
x$state_time_limit
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.