Introduction This report presents the results of a one-way sensitivity analysis comparing Empower Health versus Usual Care. The analysis evaluates how variations in key input parameters affect:

Costs for both interventions

Health outcomes (DALYs - Disability-Adjusted Life Years)

Incremental Cost-Effectiveness Ratios (ICERs)

The base case represents the reference scenario, and all parameters are varied individually to assess their impact on model outputs.

Clear the workspace

rm(list = ls())

Load Required Libraries

The script begins by loading the necessary R libraries:

library(tidyverse)
library(ggplot2)
library(scales)
library(gridExtra)
library(patchwork)
library(knitr)
library(DT)
library(kableExtra)
library(cowplot)

# Set theme for consistent plotting
theme_set(theme_minimal())

Import OWSA Input Parameters

# Import OWSA input parameters from csv file
owsa_input_data <- read_csv('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/owsa_data_12_params.csv')

# View the OWSA input parameters
head(owsa_input_data)
## # A tibble: 6 × 15
##   ...1  parameter_varied level rr_ST_D rr_PS_D rr_MI_D rr_PM_D dw_ST dw_PS dw_MI
##   <chr> <chr>            <chr>   <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 Base… base             base     7.43    2.07    5.84    2.21 0.659 0.019 0.432
## 2 rr_S… rr_ST_D          low      6.5     2.07    5.84    2.21 0.659 0.019 0.432
## 3 rr_S… rr_ST_D          high     8.5     2.07    5.84    2.21 0.659 0.019 0.432
## 4 rr_P… rr_PS_D          low      7.43    1.3     5.84    2.21 0.659 0.019 0.432
## 5 rr_P… rr_PS_D          high     7.43    3.32    5.84    2.21 0.659 0.019 0.432
## 6 rr_M… rr_MI_D          low      7.43    2.07    3.72    2.21 0.659 0.019 0.432
## # ℹ 5 more variables: dw_PM <dbl>, c_ST <dbl>, c_PS <dbl>, c_MI <dbl>,
## #   c_PM <dbl>

Import OWSA Datasets

# Import PSA datasets
load('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_results_no_trt_3.rda')
load('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_results_Emp_Health_3.rda')
load('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_results_Usual_Care_3.rda')

Rename datasets

subset_owsa_no_int <- owsa_results_no_trt_3[, c("mean_Dcosts", "mean_Ddalys")]
subset_owsa_Emp_Health <- owsa_results_Emp_Health_3[, c("mean_Dcosts", "mean_Ddalys")]
subset_owsa_UC <- owsa_results_Usual_Care_3[, c("mean_Dcosts", "mean_Ddalys")]

Add strategy to the intervention datasets

subset_owsa_no_int$strategy <- "No_Intervention"
subset_owsa_Emp_Health$strategy <- "Empower_Health"
subset_owsa_UC$strategy <- "Usual_Care"
subset_owsa_no_int$sim <- 1:nrow(subset_owsa_no_int)
subset_owsa_Emp_Health$sim <- 1:nrow(subset_owsa_Emp_Health)
subset_owsa_UC$sim <- 1:nrow(subset_owsa_UC)

Combine into one dataset for plotting

owsa_results <- bind_rows(subset_owsa_no_int, subset_owsa_Emp_Health, subset_owsa_UC)

Reshape data to wide format

owsa_wide_latest <- reshape(owsa_results, 
                timevar = "strategy", 
                idvar = "sim", 
                direction = "wide")
head(owsa_wide_latest)
##   sim mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
## 1   1                    1427.196                    7.510546
## 2   2                    1444.969                    7.501338
## 3   3                    1407.505                    7.520783
## 4   4                    1506.145                    7.474602
## 5   5                    1352.529                    7.547547
## 6   6                    1465.298                    7.487675
##   mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1                   5170.826                   7.384656               4704.749
## 2                   5187.349                   7.378609               4729.082
## 3                   5154.920                   7.390712               4684.278
## 4                   5245.941                   7.360518               4804.647
## 5                   5097.188                   7.411406               4608.070
## 6                   5212.731                   7.368221               4758.170
##   mean_Ddalys.Usual_Care
## 1               7.501453
## 2               7.492091
## 3               7.509498
## 4               7.468264
## 5               7.538061
## 6               7.479856
# # Save the psa_wide dataset as an .rda file
# save(owsa_wide_latest, file = "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_wide_latest.rda")
# # Save the psa_wide dataset as a .csv file
# write.csv(psa_wide_latest, file = "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_results_latest.csv", row.names = FALSE)

Extract first column from the OWSA input parameters dataset

# Extract the first column from the OWSA input parameters dataset-name it scenario name
owsa_input_data <- owsa_input_data %>%
  rename(scenario_name = 1)
# View the updated OWSA input parameters dataset
head(owsa_input_data)
## # A tibble: 6 × 15
##   scenario_name parameter_varied level rr_ST_D rr_PS_D rr_MI_D rr_PM_D dw_ST
##   <chr>         <chr>            <chr>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>
## 1 Base_Case     base             base     7.43    2.07    5.84    2.21 0.659
## 2 rr_ST_D_low   rr_ST_D          low      6.5     2.07    5.84    2.21 0.659
## 3 rr_ST_D_high  rr_ST_D          high     8.5     2.07    5.84    2.21 0.659
## 4 rr_PS_D_low   rr_PS_D          low      7.43    1.3     5.84    2.21 0.659
## 5 rr_PS_D_high  rr_PS_D          high     7.43    3.32    5.84    2.21 0.659
## 6 rr_MI_D_low   rr_MI_D          low      7.43    2.07    3.72    2.21 0.659
## # ℹ 7 more variables: dw_PS <dbl>, dw_MI <dbl>, dw_PM <dbl>, c_ST <dbl>,
## #   c_PS <dbl>, c_MI <dbl>, c_PM <dbl>

Extract scenario name and append to the OWSA results dataset

# Extract scenario name and append to the OWSA results dataset-assume rows are in the same order-move scenario name to the first column
owsa_wide_latest <- owsa_wide_latest %>%
  mutate(scenario_name = owsa_input_data$scenario_name) %>%
  select(scenario_name, everything())
# View the updated OWSA results dataset
head(owsa_wide_latest)
##   scenario_name sim mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
## 1     Base_Case   1                    1427.196                    7.510546
## 2   rr_ST_D_low   2                    1444.969                    7.501338
## 3  rr_ST_D_high   3                    1407.505                    7.520783
## 4   rr_PS_D_low   4                    1506.145                    7.474602
## 5  rr_PS_D_high   5                    1352.529                    7.547547
## 6   rr_MI_D_low   6                    1465.298                    7.487675
##   mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1                   5170.826                   7.384656               4704.749
## 2                   5187.349                   7.378609               4729.082
## 3                   5154.920                   7.390712               4684.278
## 4                   5245.941                   7.360518               4804.647
## 5                   5097.188                   7.411406               4608.070
## 6                   5212.731                   7.368221               4758.170
##   mean_Ddalys.Usual_Care
## 1               7.501453
## 2               7.492091
## 3               7.509498
## 4               7.468264
## 5               7.538061
## 6               7.479856

Save the updated OWSA results dataset as an .rda file

save(owsa_wide_latest, file = "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_wide_latest.rda")
# Save the updated OWSA results dataset as a .csv file
write.csv(owsa_wide_latest, file = "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_results_latest.csv", row.names = FALSE)

Load CSV File

df <- read.csv("/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/owsa_results_latest.csv")

# View structure
str(df)
## 'data.frame':    25 obs. of  8 variables:
##  $ scenario_name              : chr  "Base_Case" "rr_ST_D_low" "rr_ST_D_high" "rr_PS_D_low" ...
##  $ sim                        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ mean_Dcosts.No_Intervention: num  1427 1445 1408 1506 1353 ...
##  $ mean_Ddalys.No_Intervention: num  7.51 7.5 7.52 7.47 7.55 ...
##  $ mean_Dcosts.Empower_Health : num  5171 5187 5155 5246 5097 ...
##  $ mean_Ddalys.Empower_Health : num  7.38 7.38 7.39 7.36 7.41 ...
##  $ mean_Dcosts.Usual_Care     : num  4705 4729 4684 4805 4608 ...
##  $ mean_Ddalys.Usual_Care     : num  7.5 7.49 7.51 7.47 7.54 ...

Base Case Results The base case represents the reference scenario against which all parameter variations are compared.

# Extract base case
base_case <- df %>% filter(scenario_name == "Base_Case")

base_costs_empower <- base_case$mean_Dcosts.Empower_Health
base_dalys_empower <- base_case$mean_Ddalys.Empower_Health
base_costs_usual <- base_case$mean_Dcosts.Usual_Care
base_dalys_usual <- base_case$mean_Ddalys.Usual_Care

base_inc_cost <- base_costs_empower - base_costs_usual
base_inc_dalys <- base_dalys_empower - base_dalys_usual
base_icer <- base_inc_cost / (-base_inc_dalys)

# Create base case summary table
base_case_table <- data.frame(
  Metric = c("Empower Health Costs", 
             "Usual Care Costs",
             "Incremental Cost (Empower - Usual)",
             "Empower Health DALYs",
             "Usual Care DALYs",
             "Incremental DALYs (Empower - Usual)",
             "ICER (Cost per DALY averted)"),
  Value = c(paste0("$", round(base_costs_empower, 2)),
            paste0("$", round(base_costs_usual, 2)),
            paste0("$", round(base_inc_cost, 2)),
            round(base_dalys_empower, 4),
            round(base_dalys_usual, 4),
            round(base_inc_dalys, 4),
            paste0("$", round(base_icer, 2)))
)

kable(base_case_table, caption = "Base Case Results", booktabs = TRUE)
Base Case Results
Metric Value
Empower Health Costs $5170.83
Usual Care Costs $4704.75
Incremental Cost (Empower - Usual) $466.08
Empower Health DALYs 7.3847
Usual Care DALYs 7.5015
Incremental DALYs (Empower - Usual) -0.1168
ICER (Cost per DALY averted) $3990.47

Prepare Sensitivity Analysis Data

# Filter out base case for sensitivity analysis
df_sens <- df %>% filter(scenario_name != "Base_Case")

icer_results <- df_sens %>%
  mutate(
    # Incremental costs and DALYs
    inc_cost = mean_Dcosts.Empower_Health - mean_Dcosts.Usual_Care,
    inc_dalys = mean_Ddalys.Empower_Health - mean_Ddalys.Usual_Care,
    ICER = inc_cost / (-inc_dalys),
    
    # Extract parameter name
    parameter = case_when(
      grepl("^rr_ST_D", scenario_name) ~ "rr_ST_D",
      grepl("^rr_PS_D", scenario_name) ~ "rr_PS_D",
      grepl("^rr_MI_D", scenario_name) ~ "rr_MI_D",
      grepl("^rr_PM_D", scenario_name) ~ "rr_PM_D",
      grepl("^dw_ST", scenario_name) ~ "dw_ST",
      grepl("^dw_PS", scenario_name) ~ "dw_PS",
      grepl("^dw_MI", scenario_name) ~ "dw_MI",
      grepl("^dw_PM", scenario_name) ~ "dw_PM",
      grepl("^c_ST", scenario_name) ~ "c_ST",
      grepl("^c_PS", scenario_name) ~ "c_PS",
      grepl("^c_MI", scenario_name) ~ "c_MI",
      grepl("^c_PM", scenario_name) ~ "c_PM"
    ),
    
    # Extract level (low/high)
    level = ifelse(grepl("_low$", scenario_name), "low", "high")
  ) %>%
  select(scenario_name, parameter, level, 
         mean_Dcosts.Empower_Health, mean_Ddalys.Empower_Health,
         mean_Dcosts.Usual_Care, mean_Ddalys.Usual_Care,
         inc_cost, inc_dalys, ICER)

# Display ICER results
icer_display <- icer_results %>%
  select(parameter, level, ICER, inc_cost, inc_dalys) %>%
  mutate(ICER = round(ICER, 2),
         inc_cost = round(inc_cost, 2),
         inc_dalys = round(inc_dalys, 4))

datatable(icer_display, 
          options = list(pageLength = 12, scrollX = TRUE),
          caption = "ICER Results for All Scenarios")

Parameter Impact Analysis This section quantifies how much each parameter affects the model outputs when varied from low to high values.

calculate_impact <- function(df, outcome_col, outcome_name, base_value) {
  df %>%
    group_by(parameter) %>%
    summarise(
      low_value = .data[[outcome_col]][level == "low"],
      high_value = .data[[outcome_col]][level == "high"],
      abs_change = high_value - low_value,
      rel_change_pct = (high_value - low_value) / low_value * 100,
      low_vs_base = (.data[[outcome_col]][level == "low"] - base_value) / base_value * 100,
      high_vs_base = (.data[[outcome_col]][level == "high"] - base_value) / base_value * 100,
      .groups = "drop"
    ) %>%
    mutate(
      abs_impact = abs(abs_change),
      outcome = outcome_name,
      base_value = base_value
    )
}

# Calculate impacts for all outcomes
outcomes <- list(
  list(col = "mean_Dcosts.Empower_Health", name = "Costs_Empower_Health", base = base_costs_empower),
  list(col = "mean_Ddalys.Empower_Health", name = "DALYs_Empower_Health", base = base_dalys_empower),
  list(col = "mean_Dcosts.Usual_Care", name = "Costs_Usual_Care", base = base_costs_usual),
  list(col = "mean_Ddalys.Usual_Care", name = "DALYs_Usual_Care", base = base_dalys_usual),
  list(col = "ICER", name = "ICER", base = base_icer)
)

all_impacts <- bind_rows(lapply(outcomes, function(out) {
  calculate_impact(icer_results, out$col, out$name, out$base)
}))

# Most influential parameters
most_influential <- all_impacts %>%
  group_by(outcome) %>%
  arrange(desc(abs_impact)) %>%
  slice_head(n = 3) %>%
  select(outcome, parameter, abs_change, rel_change_pct) %>%
  mutate(rel_change_pct = round(rel_change_pct, 1),
         abs_change = round(abs_change, 2))

kable(most_influential, 
      caption = "Most Influential Parameters for Each Outcome",
      col.names = c("Outcome", "Parameter", "Absolute Change", "Relative Change (%)"),
      booktabs = TRUE)
Most Influential Parameters for Each Outcome
Outcome Parameter Absolute Change Relative Change (%)
Costs_Empower_Health c_MI 152.82 3.0
Costs_Empower_Health rr_PS_D -148.75 -2.8
Costs_Empower_Health c_ST 141.93 2.8
Costs_Usual_Care c_ST 199.04 4.3
Costs_Usual_Care c_MI 197.59 4.3
Costs_Usual_Care rr_PS_D -196.58 -4.1
DALYs_Empower_Health rr_PS_D 0.05 0.7
DALYs_Empower_Health rr_PM_D 0.04 0.6
DALYs_Empower_Health rr_MI_D 0.03 0.4
DALYs_Usual_Care rr_PS_D 0.07 0.9
DALYs_Usual_Care rr_PM_D 0.05 0.7
DALYs_Usual_Care rr_MI_D 0.04 0.5
ICER c_ST -488.96 -11.5
ICER c_MI -383.34 -9.2
ICER dw_ST -277.89 -6.7

Combined Tornado Diagrams The following figure combines all five tornado diagrams into a single comprehensive visualization, showing the sensitivity of each outcome to parameter variations. Parameters are ordered with the largest impact at the top (descending order) within each panel.

# Function to prepare tornado plot data (largest impact at TOP)
prepare_tornado_data <- function(data, outcome_col, base_value) {
  plot_data <- data %>%
    select(parameter, level, !!sym(outcome_col)) %>%
    pivot_wider(names_from = level, values_from = !!sym(outcome_col)) %>%
    mutate(
      low = as.numeric(low),
      high = as.numeric(high),
      range = abs(high - low)
    ) %>%
    filter(!is.na(range)) %>%
    arrange(desc(range)) %>%  # largest first
    mutate(
      parameter = factor(parameter, levels = parameter),
      base_value = base_value
    )
  
  return(plot_data)
}

# Function to create a single tornado plot
create_tornado_single <- function(plot_data, title, x_label, x_scale = "identity") {
  
  p <- ggplot(plot_data) +
    geom_segment(aes(x = pmin(low, high),
                     xend = pmax(low, high),
                     y = parameter,
                     yend = parameter),
                 size = 2, color = "steelblue") +
    
    geom_point(aes(x = low, y = parameter),
               size = 2.5, color = "darkred") +
    
    geom_point(aes(x = high, y = parameter),
               size = 2.5, color = "darkgreen") +
    
    geom_vline(aes(xintercept = base_value),
               linetype = "dashed",
               color = "black",
               size = 0.6) +
    
    scale_y_discrete(limits = rev(levels(plot_data$parameter))) +  # ✅ THIS LINE FIXES IT
    
    labs(title = title,
         x = x_label,
         y = "") +
    
    theme_minimal() +
    theme(
      axis.text.y = element_text(size = 8),
      plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
      axis.title.x = element_text(size = 9),
      axis.text.x = element_text(size = 8),
      plot.margin = margin(5, 5, 5, 5)
    )
  
  if (x_scale == "comma") {
    p <- p + scale_x_continuous(labels = comma)
  }
  
  return(p)
}

# Prepare datasets
data_costs_empower <- prepare_tornado_data(
  icer_results, "mean_Dcosts.Empower_Health", base_costs_empower
)

data_dalys_empower <- prepare_tornado_data(
  icer_results, "mean_Ddalys.Empower_Health", base_dalys_empower
)

data_costs_usual <- prepare_tornado_data(
  icer_results, "mean_Dcosts.Usual_Care", base_costs_usual
)

data_dalys_usual <- prepare_tornado_data(
  icer_results, "mean_Ddalys.Usual_Care", base_dalys_usual
)

data_icer <- prepare_tornado_data(
  icer_results, "ICER", base_icer
)

# Create plots
p1 <- create_tornado_single(
  data_costs_empower,
  "A) Costs - Empower Health",
  "Costs (US $)",
  "comma"
)

p2 <- create_tornado_single(
  data_dalys_empower,
  "B) DALYs - Empower Health",
  "DALYs (disability-adjusted life years)"
)

p3 <- create_tornado_single(
  data_costs_usual,
  "C) Costs - Usual Care",
  "Costs (US $)",
  "comma"
)

p4 <- create_tornado_single(
  data_dalys_usual,
  "D) DALYs - Usual Care",
  "DALYs (disability-adjusted life years)"
)

p5 <- create_tornado_single(
  data_icer,
  "E) ICER (Empower Health vs Usual Care)",
  "ICER (US$ per DALY averted)",
  "comma"
)

# Combine plots
combined_tornado <- plot_grid(
  p1, p2, p3, p4, p5,
  ncol = 2,
  nrow = 3,
  align = "v",
  rel_heights = c(1, 1, 0.8)
)

# Title
title <- ggdraw() +
  draw_label(
    "Tornado Diagrams: Parameter Impact on Key Outcomes",
    fontface = "bold",
    size = 14,
    hjust = 0.5
  ) +
  theme(plot.margin = margin(0, 0, 5, 0))

# Subtitle
subtitle <- ggdraw() +
  draw_label(
    "Parameters ordered by impact (LARGEST at TOP) | Red = Low value, Green = High value, Dashed line = Base case",
    fontface = "italic",
    size = 10,
    hjust = 0.5,
    color = "gray30"
  ) +
  theme(plot.margin = margin(0, 0, 10, 0))

# Final plot
combined_tonardo_plot <- plot_grid(
  title,
  subtitle,
  combined_tornado,
  ncol = 1,
  rel_heights = c(0.05, 0.05, 1)
)

# Display
print(combined_tonardo_plot)

# Save the combined tornado plot as a high-resolution image
ggsave("combined_tornado_plot.png", combined_tonardo_plot, width = 12, height = 15, dpi = 300)
ggsave("combined_tornado_plot.pdf", combined_tonardo_plot, width = 12, height = 15)

Detailed Individual Tornado Plots (Inverted - Largest at TOP) While the combined figure above provides a comprehensive overview, here are larger, individual tornado plots for closer inspection of each outcome. Costs - Empower Health

print(p1)

DALYs - Empower Health

print(p2)

Costs - Usual Care

print(p3)

DALYs - Usual Care

print(p4)

ICER (Cost per DALY Averted)

print(p5)

Parameter Ranking Summary Impact on ICER (Ranked by Magnitude)

parameter_icer_impact <- icer_results %>%
  group_by(parameter) %>%
  summarise(
    ICER_low = ICER[level == "low"],
    ICER_high = ICER[level == "high"],
    ICER_range = abs(ICER_high - ICER_low),
    ICER_pct_change = (ICER_high - ICER_low) / ICER_low * 100,
    more_cost_effective = ifelse(ICER_low < ICER_high, "low", "high"),
    .groups = "drop"
  ) %>%
  mutate(
    ICER_low = round(ICER_low, 2),
    ICER_high = round(ICER_high, 2),
    ICER_range = round(ICER_range, 2),
    ICER_pct_change = round(ICER_pct_change, 1)
  ) %>%
  arrange(desc(ICER_range))

kable(parameter_icer_impact, 
      caption = "Parameter Impact on ICER (Ranked by Impact)",
      col.names = c("Parameter", "ICER (Low)", "ICER (High)", 
                    "Range", "Change (%)", "More Cost-Effective Level"),
      booktabs = TRUE)
Parameter Impact on ICER (Ranked by Impact)
Parameter ICER (Low) ICER (High) Range Change (%) More Cost-Effective Level
c_ST 4234.95 3745.99 488.96 -11.5 high
c_MI 4182.14 3798.80 383.34 -9.2 high
dw_ST 4143.35 3865.46 277.89 -6.7 high
c_PS 4120.29 3860.64 259.65 -6.3 high
rr_PS_D 4095.69 3861.83 233.86 -5.7 high
dw_MI 4108.95 3876.36 232.59 -5.7 high
rr_PM_D 4021.45 3860.05 161.40 -4.0 high
rr_MI_D 4071.82 3925.63 146.18 -3.6 high
dw_PM 4036.97 3934.27 102.70 -2.5 high
rr_ST_D 4038.25 3962.09 76.16 -1.9 high
dw_PS 4016.54 3953.40 63.14 -1.6 high
c_PM 4018.11 3962.82 55.29 -1.4 high

Scenarios with Better Cost-Effectiveness Than Base

better_than_base <- icer_results %>%
  mutate(better = ICER < base_icer,
         improvement_pct = (base_icer - ICER) / base_icer * 100) %>%
  filter(better == TRUE) %>%
  select(parameter, level, ICER, improvement_pct) %>%
  mutate(ICER = round(ICER, 2),
         improvement_pct = round(improvement_pct, 1)) %>%
  arrange(improvement_pct)

if(nrow(better_than_base) > 0) {
  kable(better_than_base, 
        caption = "Scenarios with Better Cost-Effectiveness Than Base Case",
        col.names = c("Parameter", "Level", "ICER", "Improvement (%)"),
        booktabs = TRUE)
} else {
  cat("No scenarios had better cost-effectiveness than the base case.\n")
}
Scenarios with Better Cost-Effectiveness Than Base Case
Parameter Level ICER Improvement (%)
rr_ST_D high 3962.09 0.7
c_PM high 3962.82 0.7
dw_PS high 3953.40 0.9
dw_PM high 3934.27 1.4
rr_MI_D high 3925.63 1.6
dw_MI high 3876.36 2.9
dw_ST high 3865.46 3.1
rr_PS_D high 3861.83 3.2
rr_PM_D high 3860.05 3.3
c_PS high 3860.64 3.3
c_MI high 3798.80 4.8
c_ST high 3745.99 6.1