From 9512782994a018d1b040b42963fb686448b8c989 Mon Sep 17 00:00:00 2001 From: "Matthew B. Kaufmann" <57465744+mbkauf@users.noreply.github.com> Date: Thu, 30 Jul 2026 12:57:52 -0500 Subject: [PATCH 1/6] Updates to the owsa_tornado function. This includes planned changes, formatting, and package names. --- R/owsa.R | 109 +++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 82 insertions(+), 27 deletions(-) diff --git a/R/owsa.R b/R/owsa.R index 11bd642..234c125 100644 --- a/R/owsa.R +++ b/R/owsa.R @@ -156,10 +156,23 @@ plot.owsa <- function(x, txtsize = 12, #' Tornado plot of a one-way sensitivity analysis #' #' @param owsa an owsa object +#' @param return choose whether to return the plot, underlying data, or both +#' @param txtsize text size used for the plot; input for add_common_aes function #' @param min_rel_diff this function only plots #' parameters that lead to a relative change in the outcome greater than or equal #' to \code{min_rel_diff}, which must be between 0 and 1. The default (0) is that #' no strategies are filtered. +#' @param col color used for the plot; input for add_common_aes function +#' @param n_y_ticks number of y-axis ticks; input for add_common_aes function +#' @param ylim vector of y-axis limits, or NULL, which sets limits +#' automatically; input for add_common_aes function +#' @param ybreaks vector of y-axis breaks, will override \code{n_y_ticks} if +#' provided; input for add_common_aes function +#' @param select_str vector of strategy names to examine (max of 2). If none are +#' input, it will default to choose only the first of the strategies in the owsa +#' object. +#' @param outcome_name The name of the outcome to show by strategy, either IMB, +#' NMB, or ICER #' @inheritParams add_common_aes #' @inheritParams owsa_opt_strat #' @return If \code{return == "plot"}, a \code{ggplot2} tornado plot derived from the \code{owsa} @@ -172,9 +185,53 @@ plot.owsa <- function(x, txtsize = 12, owsa_tornado <- function(owsa, return = c("plot", "data"), txtsize = 12, min_rel_diff = 0, col = c("full", "bw"), - n_y_ticks = 8, ylim = NULL, ybreaks = NULL) { + n_y_ticks = 8, ylim = NULL, ybreaks = NULL, + select_str = NULL, outcome_name = NULL) { + + if (length(select_str) > 2) { + stop("Please select a max of 2 strategies") + } + + if (is.null(outcome_name)) { + outcome_name <- "Outcome" + } + + if (is.null(select_str)) { + select_str <- owsa$strategy[1] + owsa <- owsa %>% + dplyr::filter(strategy == select_str[1]) + y_label <- paste0(outcome_name, " (", select_str[1], ")") + avg <- median(owsa$outcome_val) # MK: Replacement for res_out + + } else { + if (length(select_str) == 2) { + owsa_base <- owsa %>% + dplyr::filter(strategy == select_str[1]) + owsa_comp <- owsa %>% + dplyr::filter(strategy == select_str[2]) + owsa_join <- dplyr::inner_join(owsa_base, owsa_comp, + by = c("parameter", "param_val")) + owsa_join$outcome_val <- owsa_join$outcome_val.y - owsa_join$outcome_val.x + owsa_join$strategy <- owsa_join$strategy.x + owsa <- owsa_join %>% + dplyr::select(parameter, strategy, param_val, outcome_val) + y_label <- paste0(outcome_name, " (", select_str[2], + " - ", select_str[1], ")") + # MK: Replacement for res_out + avg <- median(owsa_comp$outcome_val) - median(owsa_base$outcome_val) + + } else { + owsa <- owsa %>% + dplyr::filter(strategy == select_str) + y_label <- paste0(outcome_name, " (", select_str, ")") + avg <- median(owsa$outcome_val) # MK: Replacement for res_out + } + } + + # MK: Removed rel_diff parameter <- param_val <- outcome_val <- strategy <- outcome_val.low <- - outcome_val.high <- abs_diff <- rel_diff <- NULL + outcome_val.high <- abs_diff <- NULL + # check that is owsa object if (!is_owsa(owsa)) { stop("must provide an owsa object created with owsa()") @@ -186,40 +243,38 @@ owsa_tornado <- function(owsa, return = c("plot", "data"), } owsa_filt <- owsa %>% - group_by(parameter, param_val) %>% - arrange(outcome_val) %>% - slice(n()) %>% - select(-strategy) %>% - ungroup() + dplyr::group_by(parameter, param_val) %>% + dplyr::arrange(outcome_val) %>% + dplyr::slice(n()) %>% + dplyr::select(-strategy) %>% + dplyr::ungroup() # group by parameter and strategy mins <- owsa_filt %>% - group_by(parameter) %>% - filter(param_val == min(param_val)) + dplyr::group_by(parameter) %>% + dplyr::filter(param_val == min(param_val)) maxes <- owsa_filt %>% - group_by(parameter) %>% - filter(param_val == max(param_val)) - - avg <- median(owsa_filt$outcome_val) + dplyr::group_by(parameter) %>% + dplyr::filter(param_val == max(param_val)) - min_max <- inner_join(mins, maxes, by = c("parameter"), - suffix = c(".low", ".high")) %>% - mutate(abs_diff = abs(outcome_val.high - outcome_val.low), - rel_diff = abs_diff / outcome_val.low) %>% - filter(rel_diff >= min_rel_diff) %>% - arrange(-abs_diff) + min_max <- dplyr::inner_join(mins, maxes, by = c("parameter"), + suffix = c(".low", ".high")) %>% + dplyr::mutate(abs_diff = abs(outcome_val.high - outcome_val.low), + rel_diff = abs_diff / outcome_val.low) %>% + dplyr::arrange(-abs_diff) # return either plot or data ret <- match.arg(return) if (ret == "plot") { - g <- ggplot(min_max, aes(x = reorder(min_max$parameter, min_max$abs_diff))) + - geom_bar(aes(y = outcome_val.low, fill = "Low"), - stat = "identity") + - geom_bar(aes(y = outcome_val.high, fill = "High"), - stat = "identity") + - labs(x = "Parameter", y = "Outcome") + - coord_flip() + g <- ggplot2::ggplot(min_max, aes(x = reorder(min_max$parameter, + min_max$abs_diff))) + + ggplot2::geom_bar(aes(y = outcome_val.low, fill = "Low"), + stat = "identity") + + ggplot2::geom_bar(aes(y = outcome_val.high, fill = "High"), + stat = "identity") + + ggplot2::labs(x = "Parameter", y = y_label) + + ggplot2::coord_flip() col <- match.arg(col) g <- add_common_aes(g, txtsize, col = col, col_aes = "fill", @@ -229,7 +284,7 @@ owsa_tornado <- function(owsa, return = c("plot", "data"), n_y_ticks = n_y_ticks, ybreaks = ybreaks, ylim = ylim) + - geom_hline(yintercept = avg, linetype = 3) + ggplot2::geom_hline(yintercept = avg, linetype = 3) return(g) } else { From a0e9b7420115c87225041e629c1e5f6383f4d35d Mon Sep 17 00:00:00 2001 From: "Matthew B. Kaufmann" <57465744+mbkauf@users.noreply.github.com> Date: Thu, 30 Jul 2026 15:53:47 -0500 Subject: [PATCH 2/6] correct reference for dsa vignette. It now references the correct one --- vignettes/psa_analysis.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/psa_analysis.Rmd b/vignettes/psa_analysis.Rmd index af10f46..f6c4284 100644 --- a/vignettes/psa_analysis.Rmd +++ b/vignettes/psa_analysis.Rmd @@ -154,7 +154,7 @@ plot(el, A one-way sensitivity analysis (OWSA) illustrates how the expected outcome of a decision model changes as function of a single input parameter in the PSA. Knowing that the parameters in a decision model can interact in complex and nonlinear ways to influence outcomes, the characterization of the role of a single parameter is not straightforward. To disentangle the effect of the parameter of interest, a regression model, or 'metamodel', is fit using the entire PSA that treats the specified outcome (Net-monetary benefit, QALYs, etc) as the dependent variable and the parameter as the independent variable. By omitting the other input parameters in the PSA in the regression model, the variation in the outcome that cannot be directly attributed to variation in the parameter of interest is treated as random statistical error (even if the model is purely deterministic). Consequently, any predicted outcome value derived from the metamodel is essentially conditional on the average value of all the other PSA parameters. Separate models are produced for each strategy included in the `strategies` argument of the `owsa()` function, and these models are used to predict outcome values over a range values of the parameter of interest. -The first input of the `owsa()` function is `sa_obj`, which stands for "sensitivity analysis object", because the owsa() function can also be used in the context of a determinstic sensitivity analysis on a `dsa` object. For an illustration of this functionality, enter `vignette("dsa_visualization", package = "dampack")` in the console. The `params` argument allows the user to execute a series of separate one-way sensitivity analyses with a single function call by specifying a vector of different parameters of interest. By default, the regression model used to quantify the effect of the parameter of interest on the expected outcome value is a second order polynomial linear model. The order of the polynomial in the regression model can be altered by changing the argument `poly.order`. A higher value of `poly.order` will yield a more flexible metamodel, which can potentially lead to overfitting. +The first input of the `owsa()` function is `sa_obj`, which stands for "sensitivity analysis object", because the owsa() function can also be used in the context of a determinstic sensitivity analysis on a `dsa` object. For an illustration of this functionality, enter `vignette("dsa_generation", package = "dampack")` in the console. The `params` argument allows the user to execute a series of separate one-way sensitivity analyses with a single function call by specifying a vector of different parameters of interest. By default, the regression model used to quantify the effect of the parameter of interest on the expected outcome value is a second order polynomial linear model. The order of the polynomial in the regression model can be altered by changing the argument `poly.order`. A higher value of `poly.order` will yield a more flexible metamodel, which can potentially lead to overfitting. ```{r} o <- owsa(psa_obj) From 12d3dc9f3b89c664b99243d127a3e5ee3068e774 Mon Sep 17 00:00:00 2001 From: "Matthew B. Kaufmann" <57465744+mbkauf@users.noreply.github.com> Date: Fri, 21 Aug 2026 10:02:46 -0500 Subject: [PATCH 3/6] addition of `format.icers` to get nicely formatted ICER tables --- R/icers.R | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/R/icers.R b/R/icers.R index f2df168..242fc79 100644 --- a/R/icers.R +++ b/R/icers.R @@ -396,3 +396,43 @@ plot.icers <- function(x, } return(icer_plot) } + + +#' Format CEA table +#' +#' \code{format_table_cea} formats the CEA table. +#' +#' @param x Object of class \code{icers} +#' @param currency string. with currency used in the cost-effectiveness analysis (CEA). +#' @param dig_cost integer. number of digits to round costs to +#' @param dig_eff integer. number of digits to round effects to +#' @param dig_icer integer. number of digits to round ICER to +#' @return a dataframe object - formatted CEA table +#' @export +format.icers <- function(x, + currency = "$", + dig_cost = 0, + dig_eff = 2, + dig_icer = 0) { + + table_cea <- x + table_cea$Cost <- scales::comma(round(table_cea$Cost, dig_cost)) + table_cea$Inc_Cost <- scales::comma(round(table_cea$Inc_Cost, dig_cost)) + table_cea$Effect <- round(table_cea$Effect, dig_eff) + table_cea$Inc_Effect <- round(table_cea$Inc_Effect, dig_eff) + table_cea$ICER <- scales::comma(round(table_cea$ICER, dig_icer)) + + colnames(table_cea)[colnames(table_cea) + %in% c("Cost", + "Effect", + "Inc_Cost", + "Inc_Effect", + "ICER")] <- + c(paste0("Costs (", currency, ")"), + "QALYs", + paste0("Incremental Costs (", currency, ")"), + "Incremental QALYs", + paste0("ICER (", currency, "/QALY)")) + + return(table_cea) +} From 404d5f81862d40fb4bb72c85a00be262f7a5c35f Mon Sep 17 00:00:00 2001 From: "Matthew B. Kaufmann" <57465744+mbkauf@users.noreply.github.com> Date: Fri, 21 Aug 2026 10:07:59 -0500 Subject: [PATCH 4/6] New y-axis label --- R/ceac.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/ceac.R b/R/ceac.R index 39509ea..075624d 100644 --- a/R/ceac.R +++ b/R/ceac.R @@ -215,7 +215,7 @@ plot.ceac <- function(x, color = !!sym(strat_name))) + geom_line() + xlab(paste("Willingness to Pay (Thousand ", currency, " / QALY)", sep = "")) + - ylab("Pr Cost-Effective") + ylab("Pr Preferred") if (points) { p <- p + geom_point(aes(color = !!sym(strat_name))) From 475bf3098e08c52bedcdd00e29745cba8cb5729b Mon Sep 17 00:00:00 2001 From: "Matthew B. Kaufmann" <57465744+mbkauf@users.noreply.github.com> Date: Fri, 21 Aug 2026 11:14:35 -0500 Subject: [PATCH 5/6] change label to "Most Cost-Effective" --- R/ceac.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/ceac.R b/R/ceac.R index 075624d..aa114dc 100644 --- a/R/ceac.R +++ b/R/ceac.R @@ -215,7 +215,7 @@ plot.ceac <- function(x, color = !!sym(strat_name))) + geom_line() + xlab(paste("Willingness to Pay (Thousand ", currency, " / QALY)", sep = "")) + - ylab("Pr Preferred") + ylab("Pr Most Cost-Effective") if (points) { p <- p + geom_point(aes(color = !!sym(strat_name))) From 4402a2fc4149e11dd48ba0c486777e002c68bc39 Mon Sep 17 00:00:00 2001 From: "Matthew B. Kaufmann" <57465744+mbkauf@users.noreply.github.com> Date: Tue, 25 Aug 2026 16:08:51 -0500 Subject: [PATCH 6/6] Updated function to format the ICER table, including removing dependency on scales package, and a switch to print() from format(). --- R/icers.R | 44 ++++++++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/R/icers.R b/R/icers.R index 242fc79..6b17142 100644 --- a/R/icers.R +++ b/R/icers.R @@ -398,9 +398,7 @@ plot.icers <- function(x, } -#' Format CEA table -#' -#' \code{format_table_cea} formats the CEA table. +#' Print the formatted CEA table #' #' @param x Object of class \code{icers} #' @param currency string. with currency used in the cost-effectiveness analysis (CEA). @@ -409,18 +407,34 @@ plot.icers <- function(x, #' @param dig_icer integer. number of digits to round ICER to #' @return a dataframe object - formatted CEA table #' @export -format.icers <- function(x, - currency = "$", - dig_cost = 0, - dig_eff = 2, - dig_icer = 0) { +print.icers <- function(x, + currency = "$", + dig_cost = 0, + dig_eff = 2, + dig_icer = 0) { + # check if icers object + if (!inherits(x, "icers")) { + stop("x must be an object of class icers") + } + + format_numbers <- function(col, n_digits) { + new_col <- rep("-", length(col)) + non_na <- !is.na(col) + if(n_digits > 0) { + new_col[non_na] <- format(col[non_na], big.mark = ",", nsmall = n_digits, digits = n_digits, scientific = FALSE) + } else{ + new_col[non_na] <- format(round(col[non_na], n_digits), big.mark = ",", nsmall = n_digits, scientific = FALSE) + } + + return(new_col) + } table_cea <- x - table_cea$Cost <- scales::comma(round(table_cea$Cost, dig_cost)) - table_cea$Inc_Cost <- scales::comma(round(table_cea$Inc_Cost, dig_cost)) - table_cea$Effect <- round(table_cea$Effect, dig_eff) - table_cea$Inc_Effect <- round(table_cea$Inc_Effect, dig_eff) - table_cea$ICER <- scales::comma(round(table_cea$ICER, dig_icer)) + table_cea$Cost <- format_numbers(table_cea$Cost, n_digits = dig_cost) + table_cea$Inc_Cost <- format_numbers(table_cea$Inc_Cost, n_digits = dig_cost) + table_cea$Effect <- format_numbers(table_cea$Effect, n_digits = dig_eff) + table_cea$Inc_Effect <- format_numbers(table_cea$Inc_Effect, n_digits = dig_eff) + table_cea$ICER <- format_numbers(table_cea$ICER, n_digits = dig_icer) colnames(table_cea)[colnames(table_cea) %in% c("Cost", @@ -434,5 +448,7 @@ format.icers <- function(x, "Incremental QALYs", paste0("ICER (", currency, "/QALY)")) - return(table_cea) + class(table_cea) <- "data.frame" + print(table_cea) + invisible(table_cea) }