R/run_model_define.R

In heemod: Markov Models for Health Economic Evaluations

#' 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. See details.
#' @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.
#'   
#' @details Counting method represents where the transition should occur,
#' based on https://journals.sagepub.com/doi/10.1177/0272989X09340585:
#' "beginning" overestimates costs and "end" underestimates costs.
#' @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 = c("life-table", "beginning", "end"),
                      cost = NULL, effect = NULL,
                      state_time_limit = NULL,
                      central_strategy = NULL,
                      inflow = rep(0L, get_state_number(get_states(list(...)[[1]])))) {
  
  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 = new_quosure(substitute(cost), env = parent.frame()),
    effect = new_quosure(substitute(effect), env = parent.frame()),
    state_time_limit = state_time_limit,
    central_strategy = central_strategy,
    inflow = inflow
  )
}

#' @export
#' @rdname run_model
run_model_ <- function(uneval_strategy_list,
                       parameters,
                       init,
                       cycles,
                       method,
                       cost, effect,
                       state_time_limit,
                       central_strategy,
                       inflow) {
  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 = ", ")
    ))
  }
  
  ce <- list(
    .cost = cost,
    .effect = effect
  )
  
  
  strategy_names <- names(uneval_strategy_list)
  
  if (is.null(strategy_names)) {
    if (!identical(Sys.getenv("TESTTHAT"), "true"))
      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.")
  }
  
  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
  )
  
  eval_strategy_list <- list()
  
  for (n in names(uneval_strategy_list)) {
    eval_strategy_list[[n]] <- eval_strategy(
      strategy = uneval_strategy_list[[n]], 
      parameters = parameters,
      init = init, 
      cycles = cycles,
      method = method,
      expand_limit = state_time_limit[[n]],
      inflow = inflow,
      strategy_name = n
    )
  }
  
  list_res <- lapply(eval_strategy_list, get_total_state_values)
  
  for (n in strategy_names) {
    list_res[[n]]$.strategy_names <- n
  }
  
  p <- dplyr::bind_rows(list_res) 
  nr <- nrow(p)
  
  tab_res <- lapply(ce, function(x){
    res <- rlang::eval_tidy(x, data = p)
    if (length(res) == 1){
      return(rep(res, nr))
    }
    res
  }) 

    res <- dplyr::bind_cols(p, tab_res)
  
  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,
      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)
    ),
    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 dplyr::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 %>% 
      dplyr::arrange(.data$.cost, desc(.data$.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 %>% 
      dplyr::arrange(.data$.effect) %>% 
      dplyr::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]]) %>% 
          dplyr::mutate(.strategy_names = .n)
      }
    )
  )
  
  reshape_long(
    data = res,
    key_col = "value_names",
    value_col = "value",
    gather_cols = names(res)[! names(res) %in% 
                               c("model_time",
                                 ".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]]) %>% 
          dplyr::mutate(
            .strategy_names = .n,
            model_time = row_number())
      }
    )
  )
  
  reshape_long(
    data = res,
    key_col = "state_names",
    value_col = "count",
    gather_cols = names(res)[! names(res) %in% 
                               c("model_time",
                                 ".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), ...)
}

#' @export
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]]$complete_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
}