R/acceptability_curve.R

In heemod: Markov Models for Health Economic Evaluations

#' Acceptability Curve from Probabilistic Analysis
#' 
#' @param x Result from [run_psa()].
#' @param wtp_thresholds willingness to pay thresholds
#'   
#' @return A data frame with columns `.ceac` (the 
#'   cost-effectiveness acceptability threshold), 
#'   `.model` (treatments or models), `.n` (the 
#'   number of cases in which the treatment was most 
#'   cost-effective), and `.p` (the proportion of cases
#'   where the treatment was most effective).
#'   
#' @keywords internal
acceptability_curve <- function(x, wtp_thresholds) {
  x %>% 
    dplyr::mutate(.key = 1) %>% 
    dplyr::left_join(
      tibble::tibble(
        .ceac = wtp_thresholds,
        .key = 1
      ),
      by = ".key",
      relationship = "many-to-many"
    ) %>% 
    dplyr::group_by(.data$.index, .data$.ceac) %>% 
    dplyr::mutate(
      .nmb = .data$.effect * .data$.ceac - .data$.cost,
      .top_strategy = .data$.nmb == max(.data$.nmb),
      .top_strategy = .data$.top_strategy & cumsum(.data$.top_strategy) == 1
      # in case 2 nmb are identical, pick first
    ) %>% 
    dplyr::group_by(.data$.ceac, .data$.strategy_names) %>% 
    dplyr::summarise(.n = sum(.data$.top_strategy)) %>% 
    dplyr::group_by(.data$.ceac) %>% 
    dplyr::mutate(.p = .data$.n / sum(.data$.n))
}

generate_wtp <- function(max_wtp,
                         min_wtp = max_wtp / 1000,
                         n, log_scale) {
  stopifnot(
    max_wtp > 0
  )
  if (log_scale) {
    res <- seq(from = log(min_wtp, base = 10),
               to = log(max_wtp, base = 10),
               length.out = n)
    10 ^ res
  } else {
    seq(from = 0, to = max_wtp, length.out = n)
  }
}


globalVariables(c(     ".ceac", ".cost", ".dcost", ".deffect", ".denom", ".dref", ".effect", ".evpi", ".hjust", ".icer",
                       ".index", ".mod", ".n_indiv", ".p", ".par_names", ".par_value", ".prop", ".strategy_names",
                       ".type", ".weights", "Model", "Value", "count", "model_time", "res", "value", "maxtime", "selector"))