R/run_model_summary.R

In PolicyAnalysisInc/heRoMod: Health Economic Modeling

#**************************************************************************
#* 
#* Original work Copyright (C) 2016  Antoine Pierucci
#* Original work Copyright (C) 2017  Matt Wiener
#*
#* 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/>.
#**************************************************************************


#' @export
print.run_model <- function(x, ...) {
  print(summary(x, ...))
}

#' Summarise Markov Model Results
#' 
#' @param object Output from [run_model()].
#' @param threshold ICER threshold (possibly several) for
#'   net monetary benefit computation.
#' @param ... additional arguments affecting the summary 
#'   produced.
#'   
#' @return A `summary_run_model` object.
#' @export
summary.run_model <- function(object, threshold = NULL, strategy_order = NULL, ...) {
  if (! all(c(".cost", ".effect") %in% names(get_model_results(object)))) {
    warning("No cost and/or effect defined, model summary unavailable.")
    return(invisible(NULL))
  }
  
  res_values <- get_model_results(object) %>% 
    select(
      - .cost,
      - .effect
    ) %>% 
    as.data.frame()
  
  if (! is.null(threshold)) {
    stopifnot(
      all(threshold >= 0),
      ! any(duplicated(threshold))
    )
    
    res_nmb <- object %>% 
      scale.run_model(center = FALSE) %>% 
      select(.strategy_names, .cost, .effect) %>% 
      as.data.frame()
    
    res_nmb_strat <- character()
    
    for (tr in threshold) {
      res_nmb[format(tr, digits = 2)] <- res_nmb$.effect * tr -
        res_nmb$.cost
      
      res_nmb_strat <- c(
        res_nmb_strat,
        res_nmb$.strategy_names[res_nmb[format(tr, digits = 2)] ==
                                  max(res_nmb[format(tr, digits = 2)])][1]
      )
    }
    
  } else {
    res_nmb <- NULL
    res_nmb_strat <- NULL
  }
  
  res_comp <- object %>%
    scale.run_model(...) %>% 
    compute_icer(strategy_order = strategy_order) %>% 
    as.data.frame()
  
  structure(
    list(
      res_values = res_values,
      res_comp = res_comp,
      res_nmb = res_nmb,
      res_nmb_strat = res_nmb_strat,
      cycles = get_cycles(object),
      init = get_uneval_init(object),
      method = get_method(object),
      frontier = get_frontier(object)
    ),
    class = "summary_run_model"
  )
}

get_init.summary_run_model <- function(x) {
  x$init
}

get_model_results <- function(x) {
  UseMethod("get_model_results")
}

get_model_results.run_model <- function(x) {
  x$run_model
}

get_cost <- function(x) {
  get_model_results(x)$.cost
}

get_effect <- function(x) {
  get_model_results(x)$.effect
}

#' Normalize Cost and Effect
#' 
#' Normalize cost and effect values taking base model as a 
#' reference.
#' @name heRomod_scale
#' @param x Result of [run_model()] or 
#'   [run_psa()].
#' @param center Center results around base model?
#' @param scale Scale results to individual values?
#'   
#' @return Input with normalized `.cost` and 
#'   `.effect`, ordered by `.effect`.
#'   
#' @keywords internal
NULL

#' @rdname heRomod_scale
scale.run_model <- function(x, center = FALSE, scale = TRUE) {
  bm <- get_central_strategy(x)
  res <- tibble::tibble(
    .strategy_names = get_strategy_names(x),
    .cost = get_cost(x),
    .effect = get_effect(x),
    .n_indiv = get_n_indiv(x)
  )
  
  if (center) {
    res$.cost <- res$.cost -
      res$.cost[res$.strategy_names == bm]
    res$.effect <- res$.effect -
      res$.effect[res$.strategy_names == bm]
  }
  
  if (scale) {
    res$.cost = res$.cost / res$.n_indiv
    res$.effect = res$.effect / res$.n_indiv
  }
  res
}

#' Compute ICER
#' 
#' Compute ICER for Markov models.
#' 
#' Models are ordered by effectiveness and ICER are computed
#' sequentially.
#' 
#' @param x Result of [run_model()].
#' @param strategy_order Order in which the strategies 
#'   should be sorted. Default: by increasing effect.
#' @param threshold ICER threshold for net monetary benefit
#'   computation.
#'   
#' @return A `data.frame` with computed ICER.
#'   
#' @keywords internal
compute_icer <- function(x, strategy_order = order(x$.effect),
                         threshold = 3e4) {
  
  if (is.null(strategy_order)) strategy_order <- order(x$.effect)
  
  stopifnot(length(threshold) == 1)
  
  ef <- get_frontier(x)
  tab <- x[strategy_order, ]
  
  tab$.icer <- NA
  tab$.dcost <- NA
  tab$.deffect <- NA
  tab$.dref <- NA
  tab$.nmb <- NA
  tab$.dnmb <- NA
  
  for (i in seq_len(nrow(tab))) {
    if (i == 1) {
      tab$.icer[i] <- NA
      tab$.nmb[i] <- tab$.effect[i] * threshold - tab$.cost[i]
      ref_nmb <- tab$.nmb[i]
      ref_cost <- tab$.cost[i]
      ref_effect <- tab$.effect[i]
      ref_name <- tab$.strategy_names[i]
      
    } else {
      tab$.nmb[i] <- tab$.effect[i] * threshold - tab$.cost[i]
      tab$.dnmb[i] <- tab$.nmb[i] - ref_nmb
      tab$.dcost[i] <- tab$.cost[i] - ref_cost
      tab$.deffect[i] <- tab$.effect[i] - ref_effect
      tab$.icer[i] <- tab$.dcost[i] / tab$.deffect[i]
      tab$.dref[i] <- ref_name
      
      if (tab$.strategy_names[i] %in% ef) {
        ref_nmb <- tab$.nmb[i]
        ref_cost <- tab$.cost[i]
        ref_effect <- tab$.effect[i]
        ref_name <- tab$.strategy_names[i]
      }
    }
  }
  tab
}

compute_pw_icer <- function(deffect, dcost) {
  
  # Compare costs & outcomes
  eq_cost <- dcost == 0
  lt_cost <- dcost < 0
  lte_cost <- dcost <= 0
  gt_cost <- dcost > 0
  gte_cost <- dcost >=0 
  eq_effect <- deffect == 0
  lt_effect <- deffect < 0
  lte_effect <- deffect <= 0
  gt_effect <- deffect > 0
  gte_effect <- deffect >= 0
  
  # First Quadrant (ICER of A vs B)
  icer <- dcost / deffect
  # Origin (Identical)
  icer[eq_cost & eq_effect] <- NaN
  # Second Quadrant (Dominanted)
  icer[(lt_effect & gte_cost) | (lte_effect & gt_cost)] <- Inf
  # Fourth Quadrant (Dominant)
  icer[(gt_effect & lte_cost) | (gte_effect & lt_cost)] <- -Inf
  # Third Quadrant (ICER of B vs A)
  icer[lt_effect & lt_cost] <- -icer[lt_effect & lt_cost]
  
  return(icer)
}

format_icer <- function(icer) {
  is_identical <- is.nan(icer)
  is_dominated <- icer == Inf
  is_dominant <- icer == -Inf
  is_le_lc <- is.finite(icer)  & icer > 0
  is_me_mc <- is.finite(icer) & icer < 0
  formatted_icer <- character(length(icer))
  formatted_icer[is_identical] <- 'Identical'
  formatted_icer[is_dominated] <- 'Dominated'
  formatted_icer[is_dominant] <- 'Dominant'
  formatted_icer[is_le_lc] <- as.character(round(icer[is_le_lc],0))
  formatted_icer[is_me_mc] <- paste0(round(-icer[is_me_mc], 0), '*')
  
  return(formatted_icer)
}

#' @export
print.summary_run_model <- function(x, ...) {
  cat(sprintf(
    "%i strateg%s run for %i cycle%s.\n\n",
    nrow(x$res_values),
    plur_y(nrow(x$res_values)),
    x$cycles,
    plur(x$cycles)
  ))
  cat("Initial state counts:\n\n")
  cat(paste(
    to_text_dots(get_uneval_init(x)),
    collapse = "\n"
  ))
  cat(sprintf(
    "\n\nCounting method: '%s'.\n\n", x$method
  ))
  
  print_results(x$res_values, x$res_comp, x$res_nmb)
}

print_results <- function(res_values, res_comp, res_nmb) {
  cat("Values:\n\n")
  rownames(res_values) <- res_values$.strategy_names
  res_values <- select(res_values,
                               - .strategy_names,
                               - .n_indiv)
  print(res_values)
  
  if (! is.null(res_nmb)) {
    cat("\nNet monetary benefit difference:\n\n")
    .strategy_names <- res_nmb$.strategy_names
    f <- function(x) x - min(x)
    res_nmb <- res_nmb %>% 
      select(
        - .strategy_names,
        - .cost,
        - .effect
      ) %>% 
      mutate_all(list(f))
    
    rownames(res_nmb) <- .strategy_names
    print(res_nmb)
  }
  
  if (nrow(res_values) > 1) {
    cat("\nEfficiency frontier:\n\n")
    cat(paste(get_frontier(res_comp), collapse = " -> "))
    cat("\n\nDifferences:\n\n")
    
    rownames(res_comp) <- res_comp$.strategy_names
    res_comp <- res_comp %>% 
      select(.dcost, .deffect,
                     .icer, .dref)
    res_comp$.icer <- format(res_comp$.icer)
    res_comp$.icer[res_comp$.icer == "NA"] <- "-"
    res_comp <- res_comp[-1, ]
    res_comp <- pretty_names(res_comp)
    
    print(res_comp)
  }
}