R/run_model_define.R

In PolicyAnalysisInc/heRoMod: Health Economic Modeling

#**************************************************************************
#* 
#* 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
}