0.1 Code Description

The following R script processes the results from 1000 PSA runs. It takes in the PSA results(res_no_intervention_parallel, res_Empower_Health_parallel and res_Usual_Care_parallel) dataframes and performs subsequent analyses. The script also generates a plot of the simulation results.

1 Clear the workspace

rm(list = ls())

1.1 Load Libraries

The script begins by loading the necessary R libraries:

library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## âś” dplyr     1.1.4     âś” readr     2.1.6
## âś” forcats   1.0.1     âś” stringr   1.6.0
## âś” ggplot2   4.0.2     âś” tibble    3.3.0
## âś” lubridate 1.9.4     âś” tidyr     1.3.1
## âś” purrr     1.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## âś– dplyr::filter() masks stats::filter()
## âś– dplyr::lag()    masks stats::lag()
## â„ą Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(reshape2) # For reshaping data
## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(ggplot2) # For plotting
library(dampack) # For processing PSA results 
library(scales) # For formatting axes in plots
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor

2 Import PSA Datasets

# Import PSA datasets
load('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/psa_results_no_trt_3.rda')
load('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/psa_results_Emp_Health_3.rda')
load('/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/psa_results_Usual_Care_3.rda')

3 Rename datasets

subset_psa_no_int <- psa_results_no_trt_3[, c("mean_Dcosts", "mean_Ddalys")]
subset_psa_Emp_Health <- psa_results_Emp_Health_3[, c("mean_Dcosts", "mean_Ddalys")]
subset_psa_UC <- psa_results_Usual_Care_3[, c("mean_Dcosts", "mean_Ddalys")]

4 Add strategy to the intervention datasets

subset_psa_no_int$strategy <- "No_Intervention"
subset_psa_Emp_Health$strategy <- "Empower_Health"
subset_psa_UC$strategy <- "Usual_Care"
subset_psa_no_int$sim <- 1:nrow(subset_psa_no_int)
subset_psa_Emp_Health$sim <- 1:nrow(subset_psa_Emp_Health)
subset_psa_UC$sim <- 1:nrow(subset_psa_UC)

5 Combine into one dataset for plotting

psa_results <- bind_rows(subset_psa_no_int, subset_psa_Emp_Health, subset_psa_UC)

6 Reshape data to wide format

psa_wide_latest <- reshape(psa_results, 
                timevar = "strategy", 
                idvar = "sim", 
                direction = "wide")
head(psa_wide_latest)
##   sim mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
## 1   1                    1811.939                    7.491575
## 2   2                    1291.123                    7.551356
## 3   3                    1865.120                    7.527824
## 4   4                    1137.146                    7.479628
## 5   5                    1196.374                    7.492590
## 6   6                    1478.981                    7.503381
##   mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1                   5444.877                   7.374506               5077.466
## 2                   5068.829                   7.407125               4539.085
## 3                   5470.953                   7.398133               5089.859
## 4                   4982.714                   7.363545               4449.003
## 5                   5015.383                   7.368091               4501.201
## 6                   5212.053                   7.379397               4770.692
##   mean_Ddalys.Usual_Care
## 1               7.487483
## 2               7.541026
## 3               7.516364
## 4               7.473171
## 5               7.484676
## 6               7.494333
# Save the psa_wide dataset as an .rda file
save(psa_wide_latest, file = "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/psa_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/psa_results_latest.csv", row.names = FALSE)

7 Load CSV File

data <- read.csv("/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/psa_results_latest.csv")
# Set Kenya GDP per capita
gdp_per_capita <- 2000

# Define GDP-based WTP thresholds
wtp_thresholds_gdp <- c(
  gdp_per_capita * 0.5,   # 0.5x GDP = $1,000
  gdp_per_capita * 1,     # 1x GDP = $2,000
  gdp_per_capita * 3      # 3x GDP = $6,000
)

# Define maximum WTP for plots (up to $20,000 as requested)
max_wtp <- 10000

# Define WTP grid for analysis (0 to max_wtp)
wtp_grid <- seq(0, max_wtp, by = 500)

# View structure
str(data)
## 'data.frame':    1000 obs. of  7 variables:
##  $ sim                        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ mean_Dcosts.No_Intervention: num  1812 1291 1865 1137 1196 ...
##  $ mean_Ddalys.No_Intervention: num  7.49 7.55 7.53 7.48 7.49 ...
##  $ mean_Dcosts.Empower_Health : num  5445 5069 5471 4983 5015 ...
##  $ mean_Ddalys.Empower_Health : num  7.37 7.41 7.4 7.36 7.37 ...
##  $ mean_Dcosts.Usual_Care     : num  5077 4539 5090 4449 4501 ...
##  $ mean_Ddalys.Usual_Care     : num  7.49 7.54 7.52 7.47 7.48 ...

8 Calculate the means and standard deviations for the two interventions

# Calculate means and standard deviations for the two interventions
summary_stats <- data %>%
  summarise(
    # Means
    mean_cost_empower = mean(mean_Dcosts.Empower_Health),
    mean_cost_usual = mean(mean_Dcosts.Usual_Care),
    mean_dalys_empower = mean(mean_Ddalys.Empower_Health),
    mean_dalys_usual = mean(mean_Ddalys.Usual_Care),

    # Standard deviations
    sd_cost_empower = sd(mean_Dcosts.Empower_Health),
    sd_cost_usual = sd(mean_Dcosts.Usual_Care),
    sd_dalys_empower = sd(mean_Ddalys.Empower_Health),
    sd_dalys_usual = sd(mean_Ddalys.Usual_Care)
  )

print("Summary Statistics:")
## [1] "Summary Statistics:"
print(summary_stats)
##   mean_cost_empower mean_cost_usual mean_dalys_empower mean_dalys_usual
## 1          5181.166        4720.514           7.382205         7.498045
##   sd_cost_empower sd_cost_usual sd_dalys_empower sd_dalys_usual
## 1        272.3464       369.956       0.01988259     0.02753667
print(paste("Kenya GDP per capita: $", gdp_per_capita))
## [1] "Kenya GDP per capita: $ 2000"
print(paste("GDP-based thresholds: 0.5x = $", wtp_thresholds_gdp[1],
            ", 1x = $", wtp_thresholds_gdp[2],
            ", 3x = $", wtp_thresholds_gdp[3]))
## [1] "GDP-based thresholds: 0.5x = $ 1000 , 1x = $ 2000 , 3x = $ 6000"

9 Calculate incremental costs and DALYs

# Calculate incremental differences
incremental <- data %>%
  mutate(
    # Cost difference (Empower Health - Usual Care)
    inc_cost = mean_Dcosts.Empower_Health - mean_Dcosts.Usual_Care,

    # DALY difference (negative means DALYs averted)
    inc_dalys = mean_Ddalys.Usual_Care-mean_Ddalys.Empower_Health,

    # DALYs averted (positive means better outcome)
    dalys_averted = mean_Ddalys.Usual_Care - mean_Ddalys.Empower_Health,

    # ICER
    icer = inc_cost / dalys_averted
  )

# Summary of incremental differences and ICER
icer_summary <- incremental %>%
  summarise(
    mean_inc_cost = mean(inc_cost),
    sd_inc_cost = sd(inc_cost),
    ci_inc_cost_lower = quantile(inc_cost, 0.025),
    ci_inc_cost_upper = quantile(inc_cost, 0.975),

    mean_inc_dalys = mean(inc_dalys),
    sd_inc_dalys = sd(inc_dalys),
    ci_inc_dalys_lower = quantile(inc_dalys, 0.025),
    ci_inc_dalys_upper = quantile(inc_dalys, 0.975),

    mean_dalys_averted = mean(dalys_averted),
    sd_dalys_averted = sd(dalys_averted),
    ci_dalys_averted_lower = quantile(dalys_averted, 0.025),
    ci_dalys_averted_upper = quantile(dalys_averted, 0.975),

    mean_icer = (mean_inc_cost/mean_inc_dalys),
    median_icer = median(icer, na.rm = TRUE),
    icer_lower = quantile(icer, 0.025, na.rm = TRUE),
    icer_upper = quantile(icer, 0.975, na.rm = TRUE)
  )

print("ICER Summary (Empower Health vs Usual Care):")
## [1] "ICER Summary (Empower Health vs Usual Care):"
print(icer_summary)
##   mean_inc_cost sd_inc_cost ci_inc_cost_lower ci_inc_cost_upper mean_inc_dalys
## 1      460.6522    99.47481          228.1888          628.5122      0.1158396
##   sd_inc_dalys ci_inc_dalys_lower ci_inc_dalys_upper mean_dalys_averted
## 1  0.008999507         0.09938518          0.1338674          0.1158396
##   sd_dalys_averted ci_dalys_averted_lower ci_dalys_averted_upper mean_icer
## 1      0.008999507             0.09938518              0.1338674  3976.638
##   median_icer icer_lower icer_upper
## 1    4044.724   1933.614   5590.016

10 Test statistical differences in costs and effects

# Paired t-tests
ttest_cost <- t.test(data$mean_Dcosts.Empower_Health,
                     data$mean_Dcosts.Usual_Care, paired = TRUE)

ttest_dalys <- t.test(data$mean_Ddalys.Empower_Health,
                      data$mean_Ddalys.Usual_Care, paired = TRUE)

print("Paired t-test for Costs:")
## [1] "Paired t-test for Costs:"
print(ttest_cost)
## 
##  Paired t-test
## 
## data:  data$mean_Dcosts.Empower_Health and data$mean_Dcosts.Usual_Care
## t = 146.44, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  454.4793 466.8251
## sample estimates:
## mean difference 
##        460.6522
print("Paired t-test for DALYs:")
## [1] "Paired t-test for DALYs:"
print(ttest_dalys)
## 
##  Paired t-test
## 
## data:  data$mean_Ddalys.Empower_Health and data$mean_Ddalys.Usual_Care
## t = -407.04, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.1163981 -0.1152811
## sample estimates:
## mean difference 
##      -0.1158396

11 Plot the distribution of mean costs and DALYs for each intervention

# Use boxplots to visualize the distribution of costs and DALYs for each intervention
cost_dalys_long <- data %>%
  pivot_longer(cols = c(mean_Dcosts.Empower_Health, mean_Dcosts.Usual_Care, mean_Ddalys.Empower_Health, mean_Ddalys.Usual_Care),
               names_to = "variable", values_to = "value") %>%
  mutate(
    intervention = case_when(
      variable %in% c("mean_Dcosts.Empower_Health", "mean_Ddalys.Empower_Health") ~ "Empower Health",
      variable %in% c("mean_Dcosts.Usual_Care", "mean_Ddalys.Usual_Care") ~ "Usual Care"
    ),
    measure = case_when(
      variable %in% c("mean_Dcosts.Empower_Health", "mean_Dcosts.Usual_Care") ~ "Cost",
      variable %in% c("mean_Ddalys.Empower_Health", "mean_Ddalys.Usual_Care") ~ "DALYs"
    )
  )
results_boxplot <- ggplot(cost_dalys_long, aes(x = intervention, y = value, fill = intervention)) +
  geom_boxplot() +
  facet_wrap(~ measure, scales = "free_y") +
  labs(title = "Distribution of Mean Costs and DALYs by Intervention",
       x = "Intervention",
       y = "Value") +
  theme_minimal() +
  scale_fill_manual(values = c("Empower Health" = "blue", "Usual Care" = "green")) +
  theme(legend.position = "none")

# View Plot
print(results_boxplot)

# Plot cost-effectiveness plane for mean costs and DALYs

# Mean values
mean_uc_dalys <- mean(data$mean_Ddalys.Usual_Care, na.rm = TRUE)
mean_uc_costs <- mean(data$mean_Dcosts.Usual_Care, na.rm = TRUE)

mean_eh_dalys <- mean(data$mean_Ddalys.Empower_Health, na.rm = TRUE)
mean_eh_costs <- mean(data$mean_Dcosts.Empower_Health, na.rm = TRUE)



# Cost-effectiveness scatter plot with legend
ce_scatter <- ggplot() +
  geom_point(
    data = data,
    aes(
      x = mean_Ddalys.Usual_Care,
      y = mean_Dcosts.Usual_Care,
      color = "Usual Care"
    ),
    alpha = 0.3, size = 1
  ) +
  geom_point(
    data = data,
    aes(
      x = mean_Ddalys.Empower_Health,
      y = mean_Dcosts.Empower_Health,
      color = "Empower Health"
    ),
    alpha = 0.3, size = 1
  ) +
  geom_point(
    aes(
      x = mean_uc_dalys,
      y = mean_uc_costs,
      color = "Usual Care"
    ),
    size = 4, shape = 18
  ) +
  geom_point(
    aes(
      x = mean_eh_dalys,
      y = mean_eh_costs,
      color = "Empower Health"
    ),
    size = 4, shape = 18
  ) +
  labs(
    title = "Cost-Effectiveness Scatter Plot (Costs vs DALYs)",
    x = "DALYs",
    y = "Costs (USD)",
    color = "Intervention"
  ) +
  scale_color_manual(
    values = c(
      "Usual Care" = "green",
      "Empower Health" = "blue"
    )
  ) +
  scale_x_continuous(labels = number_format(accuracy = 0.001)) +
  scale_y_continuous(labels = dollar_format()) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(ce_scatter)

# Plot CE Plane for ICERs

# Optional helpers for dynamic annotation placement
x_pos <- max(incremental$dalys_averted, na.rm = TRUE) * 0.7
y_max <- max(incremental$inc_cost, na.rm = TRUE)
x_max <- max(incremental$dalys_averted, na.rm = TRUE) * 1.2
xmax <- x_max
mean_icer <- mean(incremental$inc_cost) / mean(incremental$dalys_averted)
icer_scatter <- ggplot(incremental, aes(x = dalys_averted, y = inc_cost)) +
  geom_point(alpha = 0.3, color = "darkgreen", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50", linewidth = 0.6) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50", linewidth = 0.6) +

  # GDP-based WTP lines
  geom_abline(
    intercept = 0, slope = wtp_thresholds_gdp[1],
    linetype = "dotted", color = "orange", alpha = 0.7, linewidth = 1
  ) +
  geom_abline(
    intercept = 0, slope = wtp_thresholds_gdp[2],
    linetype = "dotted", color = "blue", alpha = 0.7, linewidth = 1
  ) +
  geom_abline(
    intercept = 0, slope = wtp_thresholds_gdp[3],
    linetype = "dotted", color = "red", alpha = 0.7, linewidth = 1
  ) +

  # Mean point
  geom_point(
    x = mean(incremental$dalys_averted, na.rm = TRUE),
    y = mean(incremental$inc_cost, na.rm = TRUE),
    color = "black", size = 4, shape = 18
  ) +

  # Annotations
  annotate(
    "text", x = x_pos, y = wtp_thresholds_gdp[1] * x_pos,
    label = paste0("0.5x GDP = $", format(wtp_thresholds_gdp[1], big.mark = ",")),
    angle = 20, vjust = -0.5, size = 3, color = "orange"
  ) +
  annotate(
    "text", x = x_pos, y = wtp_thresholds_gdp[2] * x_pos,
    label = paste0("1x GDP = $", format(wtp_thresholds_gdp[2], big.mark = ",")),
    angle = 20, vjust = -0.5, size = 3, color = "blue"
  ) +
  annotate(
    "text", x = x_pos, y = wtp_thresholds_gdp[3] * x_pos,
    label = paste0("3x GDP = $", format(wtp_thresholds_gdp[3], big.mark = ",")),
    angle = 20, vjust = -0.5, size = 3, color = "red"
  ) +
  # Mean ICER label
  annotate("text",
           x = xmax * 0.55,
           y = mean_icer * xmax * 0.55,
           label = paste0("ICER = ", dollar(mean_icer), "/DALY averted"),
           color = "black", size = 3.5, hjust = 0, vjust = 0, fontface = "bold") +

  labs(
    title = "Incremental Cost-Effectiveness Scatter Plot",
    subtitle = paste0(
      "Empower Health vs Usual Care | Kenya GDP thresholds: 0.5x = $",
      format(wtp_thresholds_gdp[1], big.mark = ","),
      ", 1x = $", format(wtp_thresholds_gdp[2], big.mark = ","),
      ", 3x = $", format(wtp_thresholds_gdp[3], big.mark = ",")
    ),
    x = "Incremental DALYs Averted (Empower Health better →)",
    y = "Incremental Cost (USD) (Empower Health more expensive →)"
  ) +
  scale_x_continuous(labels = number_format(accuracy = 0.001)) +
  scale_y_continuous(labels = dollar_format()) +
  theme_minimal()

print(icer_scatter)

# Plot CEAC for the two interventions

# 3. Cost-Effectiveness Acceptability Curve with GDP thresholds
# Function to calculate Net Monetary Benefit (NMB)
calculate_nmb <- function(cost, dalys, wtp) {
  return(wtp * (7.8 - dalys) - cost)  # Using 7.8 as reference (max DALYs in dataset)
}

# Calculate NMB at different WTP thresholds
nmb_data <- data %>%
  crossing(wtp = wtp_grid) %>%
  mutate(
    nmb_usual = calculate_nmb(mean_Dcosts.Usual_Care,
                              mean_Ddalys.Usual_Care, wtp),
    nmb_empower = calculate_nmb(mean_Dcosts.Empower_Health,
                                mean_Ddalys.Empower_Health, wtp),
    nmb_diff = nmb_empower - nmb_usual
  )

# Calculate probability of being cost-effective at each threshold
ce_acceptability <- nmb_data %>%
  group_by(wtp) %>%
  summarise(
    prob_usual = mean(nmb_usual > nmb_empower),
    prob_empower = mean(nmb_empower > nmb_usual),
    prob_equal = mean(nmb_empower == nmb_usual)
  ) %>%
  pivot_longer(cols = c(prob_usual, prob_empower),
               names_to = "intervention",
               values_to = "probability") %>%
  mutate(
    intervention = case_when(
      intervention == "prob_usual" ~ "Usual Care",
      intervention == "prob_empower" ~ "Empower Health"
    )
  )
ce_acceptability
## # A tibble: 42 Ă— 4
##      wtp prob_equal intervention   probability
##    <dbl>      <dbl> <chr>                <dbl>
##  1     0          0 Usual Care           1    
##  2     0          0 Empower Health       0    
##  3   500          0 Usual Care           1    
##  4   500          0 Empower Health       0    
##  5  1000          0 Usual Care           0.995
##  6  1000          0 Empower Health       0.005
##  7  1500          0 Usual Care           0.99 
##  8  1500          0 Empower Health       0.01 
##  9  2000          0 Usual Care           0.974
## 10  2000          0 Empower Health       0.026
## # â„ą 32 more rows
ceac_plot <- ggplot(ce_acceptability, aes(x = wtp, y = probability, color = intervention)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = 0, ymax = probability, fill = intervention),
              alpha = 0.1, show.legend = FALSE) +
  # Add vertical lines for GDP thresholds
  geom_vline(xintercept = wtp_thresholds_gdp[1], linetype = "dashed",
             color = "orange", alpha = 0.7, size = 0.8) +
  geom_vline(xintercept = wtp_thresholds_gdp[2], linetype = "dashed",
             color = "blue", alpha = 0.7, size = 0.8) +
  geom_vline(xintercept = wtp_thresholds_gdp[3], linetype = "dashed",
             color = "red", alpha = 0.7, size = 0.8) +
  # Add labels for GDP thresholds
  annotate("text", x = wtp_thresholds_gdp[1] + 300, y = 0.95,
           label = paste0("0.5x GDP\n$", format(wtp_thresholds_gdp[1], big.mark = ",")),
           size = 3, color = "orange") +
  annotate("text", x = wtp_thresholds_gdp[2] + 300, y = 0.95,
           label = paste0("1x GDP\n$", format(wtp_thresholds_gdp[2], big.mark = ",")),
           size = 3, color = "blue") +
  annotate("text", x = wtp_thresholds_gdp[3] + 300, y = 0.95,
           label = paste0("3x GDP\n$", format(wtp_thresholds_gdp[3], big.mark = ",")),
           size = 3, color = "red") +
  labs(
    title = "Cost-Effectiveness Acceptability Curve",
    subtitle = paste0("Kenya GDP per capita = $", gdp_per_capita,
                      " | Thresholds: 0.5x, 1x, and 3x GDP"),
    x = "Willingness-to-Pay per DALY Averted (USD)",
    y = "Probability of Cost-Effectiveness",
    color = "Intervention"
  ) +
  scale_x_continuous(labels = scales::dollar_format(),
                     limits = c(0, max_wtp),
                     breaks = seq(0, max_wtp, by = 2000)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  theme_minimal() +
  theme(legend.position = "bottom")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## â„ą Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(ceac_plot)

# Combine relevant plots into a single figure for comparison

# Combine the state occupancy, CVD events, costs, and effectiveness plots for both interventions
combined_plot_psa <- cowplot::plot_grid(
  results_boxplot, ce_scatter, icer_scatter, ceac_plot, #first_instance_plot,  
  labels = c("A", "B", "C", "D", "E"),
  ncol = 2,
  align = "hv"
)
## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.
## Warning: Graphs cannot be horizontally aligned unless the axis parameter is
## set. Placing graphs unaligned.
# View the combined plot
print(combined_plot_psa)

# Save the combined plot
ggsave("combined_comparison_plot_psa.png", plot = combined_plot_psa, width = 15, height = 10)
ggsave("combined_comparison_plot_psa.pdf", plot = combined_plot_psa, width = 15, height = 10)
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## for 'Incremental DALYs Averted (Empower Health better →)' in 'mbcsToSbcs': ->
## substituted for → (U+2192)
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## for 'Incremental Cost (USD) (Empower Health more expensive →)' in 'mbcsToSbcs':
## -> substituted for → (U+2192)

12 NMB Plot with GDP thresholds

# 4. Net Monetary Benefit Plot with GDP thresholds
nmb_summary <- nmb_data %>%
  group_by(wtp) %>%
  summarise(
    nmb_usual_mean = mean(nmb_usual),
    nmb_usual_lower = quantile(nmb_usual, 0.025),
    nmb_usual_upper = quantile(nmb_usual, 0.975),
    nmb_empower_mean = mean(nmb_empower),
    nmb_empower_lower = quantile(nmb_empower, 0.025),
    nmb_empower_upper = quantile(nmb_empower, 0.975)
  ) %>%
  pivot_longer(cols = -wtp,
               names_to = c("intervention", "statistic"),
               names_pattern = "nmb_(.*)_(.*)") %>%
  pivot_wider(names_from = statistic, values_from = value)

nmb_plot <- ggplot(nmb_summary, aes(x = wtp, y = mean, color = intervention, fill = intervention)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
  # Add vertical lines for GDP thresholds
  geom_vline(xintercept = wtp_thresholds_gdp[1], linetype = "dashed",
             color = "orange", alpha = 0.7, size = 0.8) +
  geom_vline(xintercept = wtp_thresholds_gdp[2], linetype = "dashed",
             color = "blue", alpha = 0.7, size = 0.8) +
  geom_vline(xintercept = wtp_thresholds_gdp[3], linetype = "dashed",
             color = "red", alpha = 0.7, size = 0.8) +
  labs(
    title = "Net Monetary Benefit by Willingness-to-Pay Threshold",
    subtitle = paste0("Empower Health vs Usual Care | Kenya GDP thresholds (0.5x, 1x, 3x = $",
                      wtp_thresholds_gdp[1], ", $", wtp_thresholds_gdp[2], ", $",
                      wtp_thresholds_gdp[3], ")"),
    x = "Willingness-to-Pay per DALY Averted (USD)",
    y = "Net Monetary Benefit (USD)",
    color = "Intervention",
    fill = "Intervention"
  ) +
  scale_x_continuous(labels = scales::dollar_format(),
                     limits = c(0, max_wtp),
                     breaks = seq(0, max_wtp, by = 2000)) +
  scale_y_continuous(labels = scales::dollar_format()) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50")

print(nmb_plot)

# Bootstrap confidence intervals for ICER

# 6. Sensitivity Analysis: Bootstrap confidence intervals for ICER
set.seed(123)
n_bootstrap <- 5000

bootstrap_icers <- function(data, n_bootstrap) {
  boot_icers <- numeric(n_bootstrap)

  for (i in 1:n_bootstrap) {
    boot_sample <- data[sample(1:nrow(data), nrow(data), replace = TRUE), ]

    # Calculate ICER for bootstrap sample
    dalys_averted <- mean(boot_sample$mean_Ddalys.Usual_Care - boot_sample$mean_Ddalys.Empower_Health)
    cost_diff <- mean(boot_sample$mean_Dcosts.Empower_Health - boot_sample$mean_Dcosts.Usual_Care)
    boot_icers[i] <- cost_diff / dalys_averted
  }

  return(boot_icers)
}

boot_icers <- bootstrap_icers(data, n_bootstrap)

# Summary of bootstrap ICERs
boot_icer_summary <- data.frame(
  mean_icer = mean(boot_icers, na.rm = TRUE),
  median_icer = median(boot_icers, na.rm = TRUE),
  lower_95 = quantile(boot_icers, 0.025, na.rm = TRUE),
  upper_95 = quantile(boot_icers, 0.975, na.rm = TRUE),
  prob_cost_effective_at_0.5x_GDP = mean(boot_icers < wtp_thresholds_gdp[1], na.rm = TRUE),
  prob_cost_effective_at_1x_GDP = mean(boot_icers < wtp_thresholds_gdp[2], na.rm = TRUE),
  prob_cost_effective_at_3x_GDP = mean(boot_icers < wtp_thresholds_gdp[3], na.rm = TRUE)
)

print("Bootstrap ICER Confidence Intervals:")
## [1] "Bootstrap ICER Confidence Intervals:"
print(boot_icer_summary)
##      mean_icer median_icer lower_95 upper_95 prob_cost_effective_at_0.5x_GDP
## 2.5%  3976.661    3976.792 3921.322 4029.796                               0
##      prob_cost_effective_at_1x_GDP prob_cost_effective_at_3x_GDP
## 2.5%                             0                             1

13 Cost-Effectiveness Plane with Quadrants and GDP thresholds

mean_dalys <- mean(incremental$dalys_averted, na.rm = TRUE)
mean_cost <- mean(incremental$inc_cost, na.rm = TRUE)
xmax <- max(abs(incremental$dalys_averted), na.rm = TRUE) * 1.5
ymax <- max(abs(incremental$inc_cost), na.rm = TRUE) * 1.5
mean_icer <- mean(incremental$inc_cost) / mean(incremental$dalys_averted)



ce_plane <- ggplot(incremental,
                   aes(x = dalys_averted, y = inc_cost)) +
  geom_point(alpha = 0.3, color = "darkgreen", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  # Add GDP thresholds
  geom_abline(intercept = 0, slope = wtp_thresholds_gdp[1],
              linetype = "dotted", color = "orange", alpha = 0.7, size = 1) +
  geom_abline(intercept = 0, slope = wtp_thresholds_gdp[2],
              linetype = "dotted", color = "blue", alpha = 0.7, size = 1) +
  geom_abline(intercept = 0, slope = wtp_thresholds_gdp[3],
              linetype = "dotted", color = "red", alpha = 0.7, size = 1) +
  # Mean incremental point
  geom_point(x = mean_dalys, y = mean_cost,
             color = "black", size = 4, shape = 18) +

  # Mean ICER label
  annotate("text",
           x = xmax * 0.55,
           y = mean_icer * xmax * 0.55,
           label = paste0("ICER = ", dollar(mean_icer), "/DALY averted"),
           color = "black", size = 3, hjust = 0, vjust = -0.5) +
  # Add quadrants
  annotate("rect", xmin = 0, xmax = Inf, ymin = 0, ymax = Inf,
           alpha = 0.1, fill = "lightgreen") +
  annotate("rect", xmin = -Inf, xmax = 0, ymin = 0, ymax = Inf,
           alpha = 0.1, fill = "lightyellow") +
  annotate("rect", xmin = -Inf, xmax = 0, ymin = -Inf, ymax = 0,
           alpha = 0.1, fill = "lightcoral") +
  annotate("rect", xmin = 0, xmax = Inf, ymin = -Inf, ymax = 0,
           alpha = 0.1, fill = "lightblue") +
  # Add labels
  annotate("text", x = 0.04, y = 400, label = "Quadrant I\nCostly, Effective",
           size = 3, color = "darkgreen") +
  annotate("text", x = -0.02, y = 400, label = "Quadrant II\nCostly, Ineffective",
           size = 3, color = "darkred") +
  annotate("text", x = -0.02, y = -400, label = "Quadrant III\nCost-saving, Ineffective",
           size = 3, color = "darkred") +
  annotate("text", x = 0.04, y = -400, label = "Quadrant IV\nCost-saving, Effective",
           size = 3, color = "darkgreen") +
  labs(
    title = "Cost-Effectiveness Plane with Quadrants",
    subtitle = paste0("Empower Health vs Usual Care | Kenya GDP thresholds (0.5x, 1x, 3x = $",
                      wtp_thresholds_gdp[1], ", $", wtp_thresholds_gdp[2], ", $",
                      wtp_thresholds_gdp[3], ")"),
    x = "Incremental DALYs Averted (Empower Health better →)",
    y = "Incremental Cost (USD) (Empower Health more expensive →)"
  ) +
  scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
  scale_y_continuous(labels = scales::dollar_format()) +
  theme_minimal()

print(ce_plane)

# EVPI

# 7. Expected Value of Perfect Information (EVPI)
calculate_evpi <- function(data, wtp_grid) {
  evpi_results <- data.frame(wtp = wtp_grid, evpi = NA)

  for (i in seq_along(wtp_grid)) {
    wtp <- wtp_grid[i]

    # Calculate NMB for each intervention at this WTP
    nmb_data <- data %>%
      mutate(
        nmb_usual = wtp * (7.8 - mean_Ddalys.Usual_Care) - mean_Dcosts.Usual_Care,
        nmb_empower = wtp * (7.8 - mean_Ddalys.Empower_Health) - mean_Dcosts.Empower_Health
      ) %>%
      rowwise() %>%
      mutate(
        max_nmb = max(nmb_usual, nmb_empower)
      )

    # EVPI = Expected value with perfect information - Expected value with current information
    current_expected_value <- max(mean(nmb_data$nmb_usual), mean(nmb_data$nmb_empower))
    expected_value_perfect_info <- mean(nmb_data$max_nmb)

    evpi_results$evpi[i] <- expected_value_perfect_info - current_expected_value
  }

  return(evpi_results)
}

# Calculate EVPI using the same grid
evpi_results <- calculate_evpi(data, wtp_grid)

# EVPI Plot with GDP thresholds
evpi_plot <- ggplot(evpi_results, aes(x = wtp, y = evpi)) +
  geom_line(size = 1.2, color = "darkblue") +
  geom_ribbon(aes(ymin = 0, ymax = evpi), alpha = 0.3, fill = "darkblue") +
  # Add vertical lines for GDP thresholds
  geom_vline(xintercept = wtp_thresholds_gdp[1], linetype = "dashed",
             color = "orange", alpha = 0.7, size = 0.8) +
  geom_vline(xintercept = wtp_thresholds_gdp[2], linetype = "dashed",
             color = "blue", alpha = 0.7, size = 0.8) +
  geom_vline(xintercept = wtp_thresholds_gdp[3], linetype = "dashed",
             color = "red", alpha = 0.7, size = 0.8) +
  # Add labels
  annotate("text", x = wtp_thresholds_gdp[1] + 300, y = max(evpi_results$evpi) * 0.9,
           label = paste0("0.5x GDP\n$", format(wtp_thresholds_gdp[1], big.mark = ",")),
           size = 3, color = "orange") +
  annotate("text", x = wtp_thresholds_gdp[2] + 300, y = max(evpi_results$evpi) * 0.8,
           label = paste0("1x GDP\n$", format(wtp_thresholds_gdp[2], big.mark = ",")),
           size = 3, color = "blue") +
  annotate("text", x = wtp_thresholds_gdp[3] + 300, y = max(evpi_results$evpi) * 0.7,
           label = paste0("3x GDP\n$", format(wtp_thresholds_gdp[3], big.mark = ",")),
           size = 3, color = "red") +
  labs(
    title = "Expected Value of Perfect Information (EVPI)",
    subtitle = paste0("Kenya GDP per capita = $", gdp_per_capita,
                      " | Thresholds: 0.5x, 1x, and 3x GDP"),
    x = "Willingness-to-Pay per DALY Averted (USD)",
    y = "EVPI (USD per patient)"
  ) +
  scale_x_continuous(labels = scales::dollar_format(),
                     limits = c(0, max_wtp),
                     breaks = seq(0, max_wtp, by = 2000)) +
  scale_y_continuous(labels = scales::dollar_format()) +
  theme_minimal() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50")

print(evpi_plot)

# Summary table with GDP thresholds

summary_table <- data.frame(
  Metric = c("Kenya GDP per capita",
             "0.5x GDP Threshold", "1x GDP Threshold", "3x GDP Threshold",
             "Mean Cost - Empower Health", "Mean Cost - Usual Care",
             "Mean DALYs - Empower Health", "Mean DALYs - Usual Care",
             "Mean Incremental Cost", "Mean DALYs Averted",
             "Mean ICER", "Median ICER",
             # "95% CI ICER (Bootstrap)",
             paste0("Probability ICER < 0.5x GDP ($", format(wtp_thresholds_gdp[1], big.mark = ","), ")"),
             paste0("Probability ICER < 1x GDP ($", format(wtp_thresholds_gdp[2], big.mark = ","), ")"),
             paste0("Probability ICER < 3x GDP ($", format(wtp_thresholds_gdp[3], big.mark = ","), ")"),
             paste0("EVPI at 0.5x GDP ($", format(wtp_thresholds_gdp[1], big.mark = ","), ")"),
             paste0("EVPI at 1x GDP ($", format(wtp_thresholds_gdp[2], big.mark = ","), ")"),
             paste0("EVPI at 3x GDP ($", format(wtp_thresholds_gdp[3], big.mark = ","), ")")),
  Value = c(
    sprintf("$%s", format(gdp_per_capita, big.mark = ",")),
    sprintf("$%s", format(wtp_thresholds_gdp[1], big.mark = ",")),
    sprintf("$%s", format(wtp_thresholds_gdp[2], big.mark = ",")),
    sprintf("$%s", format(wtp_thresholds_gdp[3], big.mark = ",")),
    sprintf("$%s", format(round(mean(data$mean_Dcosts.Empower_Health), 2), big.mark = ",")),
    sprintf("$%s", format(round(mean(data$mean_Dcosts.Usual_Care), 2), big.mark = ",")),
    round(mean(data$mean_Ddalys.Empower_Health), 4),
    round(mean(data$mean_Ddalys.Usual_Care), 4),
    sprintf("$%s", format(round(icer_summary$mean_inc_cost, 2), big.mark = ",")),
    round(icer_summary$mean_dalys_averted, 4),
    sprintf("$%s", format(round(icer_summary$mean_icer, 0), big.mark = ",")),
    sprintf("$%s", format(round(icer_summary$median_icer, 0), big.mark = ",")),
    # sprintf("$%s - $%s",
    #         format(round(boot_icer_summary$lower_95, 0), big.mark = ","),
    #         format(round(boot_icer_summary$upper_95, 0), big.mark = ",")),
    sprintf("%.1f%%", boot_icer_summary$prob_cost_effective_at_0.5x_GDP * 100),
    sprintf("%.1f%%", boot_icer_summary$prob_cost_effective_at_1x_GDP * 100),
    sprintf("%.1f%%", boot_icer_summary$prob_cost_effective_at_3x_GDP * 100),
    sprintf("$%s", format(round(evpi_results$evpi[evpi_results$wtp == wtp_thresholds_gdp[1]], 2), big.mark = ",")),
    sprintf("$%s", format(round(evpi_results$evpi[evpi_results$wtp == wtp_thresholds_gdp[2]], 2), big.mark = ",")),
    sprintf("$%s", format(round(evpi_results$evpi[evpi_results$wtp == wtp_thresholds_gdp[3]], 2), big.mark = ","))
  )
)

print("Summary Table - Empower Health vs Usual Care with Kenya GDP Thresholds:")
## [1] "Summary Table - Empower Health vs Usual Care with Kenya GDP Thresholds:"
print(summary_table)
##                                  Metric     Value
## 1                  Kenya GDP per capita    $2,000
## 2                    0.5x GDP Threshold    $1,000
## 3                      1x GDP Threshold    $2,000
## 4                      3x GDP Threshold    $6,000
## 5            Mean Cost - Empower Health $5,181.17
## 6                Mean Cost - Usual Care $4,720.51
## 7           Mean DALYs - Empower Health    7.3822
## 8               Mean DALYs - Usual Care     7.498
## 9                 Mean Incremental Cost   $460.65
## 10                   Mean DALYs Averted    0.1158
## 11                            Mean ICER    $3,977
## 12                          Median ICER    $4,045
## 13 Probability ICER < 0.5x GDP ($1,000)      0.0%
## 14   Probability ICER < 1x GDP ($2,000)      0.0%
## 15   Probability ICER < 3x GDP ($6,000)    100.0%
## 16            EVPI at 0.5x GDP ($1,000)     $0.16
## 17              EVPI at 1x GDP ($2,000)     $1.63
## 18              EVPI at 3x GDP ($6,000)     $0.12

14 Plotting ICER Stabilization curves to diagnose the number of PSA iterations needed

# View the data structure
head(incremental)
##   sim mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
## 1   1                    1811.939                    7.491575
## 2   2                    1291.123                    7.551356
## 3   3                    1865.120                    7.527824
## 4   4                    1137.146                    7.479628
## 5   5                    1196.374                    7.492590
## 6   6                    1478.981                    7.503381
##   mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1                   5444.877                   7.374506               5077.466
## 2                   5068.829                   7.407125               4539.085
## 3                   5470.953                   7.398133               5089.859
## 4                   4982.714                   7.363545               4449.003
## 5                   5015.383                   7.368091               4501.201
## 6                   5212.053                   7.379397               4770.692
##   mean_Ddalys.Usual_Care inc_cost inc_dalys dalys_averted     icer
## 1               7.487483 367.4107 0.1129769     0.1129769 3252.088
## 2               7.541026 529.7438 0.1339012     0.1339012 3956.229
## 3               7.516364 381.0940 0.1182311     0.1182311 3223.296
## 4               7.473171 533.7110 0.1096266     0.1096266 4868.445
## 5               7.484676 514.1812 0.1165855     0.1165855 4410.334
## 6               7.494333 441.3610 0.1149356     0.1149356 3840.071
summary(incremental)
##       sim         mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
##  Min.   :   1.0   Min.   : 554.8              Min.   :7.410              
##  1st Qu.: 250.8   1st Qu.:1168.0              1st Qu.:7.487              
##  Median : 500.5   Median :1405.2              Median :7.506              
##  Mean   : 500.5   Mean   :1442.0              Mean   :7.506              
##  3rd Qu.: 750.2   3rd Qu.:1672.7              3rd Qu.:7.526              
##  Max.   :1000.0   Max.   :3002.2              Max.   :7.633              
##  mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
##  Min.   :4567               Min.   :7.318              Min.   :3880          
##  1st Qu.:4989               1st Qu.:7.369              1st Qu.:4463          
##  Median :5152               Median :7.382              Median :4684          
##  Mean   :5181               Mean   :7.382              Mean   :4721          
##  3rd Qu.:5344               3rd Qu.:7.395              3rd Qu.:4946          
##  Max.   :6219               Max.   :7.478              Max.   :6154          
##  mean_Ddalys.Usual_Care    inc_cost        inc_dalys       dalys_averted    
##  Min.   :7.410          Min.   : 65.02   Min.   :0.08889   Min.   :0.08889  
##  1st Qu.:7.480          1st Qu.:400.93   1st Qu.:0.10965   1st Qu.:0.10965  
##  Median :7.498          Median :472.93   Median :0.11567   Median :0.11567  
##  Mean   :7.498          Mean   :460.65   Mean   :0.11584   Mean   :0.11584  
##  3rd Qu.:7.516          3rd Qu.:531.92   3rd Qu.:0.12199   3rd Qu.:0.12199  
##  Max.   :7.616          Max.   :686.46   Max.   :0.14346   Max.   :0.14346  
##       icer     
##  Min.   : 513  
##  1st Qu.:3423  
##  Median :4045  
##  Mean   :3997  
##  3rd Qu.:4644  
##  Max.   :6374
# Calculate basic statistics
total_runs <- nrow(incremental)
final_icer <- mean(incremental$inc_cost) / mean(incremental$dalys_averted)
final_cost <- mean(incremental$inc_cost)
final_effect <- mean(incremental$dalys_averted)

cat("=== PSA DATA SUMMARY ===\n")
## === PSA DATA SUMMARY ===
cat("Total PSA runs:", total_runs, "\n")
## Total PSA runs: 1000
cat("Mean incremental cost:", dollar(final_cost), "\n")
## Mean incremental cost: $460.65
cat("Mean DALYs averted:", round(final_effect, 4), "\n")
## Mean DALYs averted: 0.1158
cat("Final ICER (all runs):", dollar(final_icer), "/DALY\n\n")
## Final ICER (all runs): $3,976.64 /DALY
# Define sample sizes for cumulative calculation
sample_sizes <- c(1, 5, 10, 25, 50, 75, 100, 150, 200, 250, 300, 350, 400, 
                  450, 500, 600, 700, 800, 900, total_runs)
sample_sizes <- unique(sample_sizes[sample_sizes <= total_runs])

# Calculate cumulative means
cumulative_data <- data.frame(
  n_runs = sample_sizes,
  cumulative_icer = sapply(sample_sizes, function(n) {
    mean(incremental$inc_cost[1:n]) / mean(incremental$dalys_averted[1:n])
  }),
  cumulative_cost = sapply(sample_sizes, function(n) {
    mean(incremental$inc_cost[1:n])
  }),
  cumulative_effect = sapply(sample_sizes, function(n) {
    mean(incremental$dalys_averted[1:n])
  })
)

# Calculate stabilization bands
band_5pct <- abs(final_icer * 0.05)
band_10pct <- abs(final_icer * 0.10)

# --------------------------------
# BEST OPTION: Single Faceted Plot (No patchwork needed)
# --------------------------------

# Reshape data for faceted plot
faceted_data <- cumulative_data %>%
  select(n_runs, cumulative_icer, cumulative_cost, cumulative_effect) %>%
  pivot_longer(cols = c(cumulative_icer, cumulative_cost, cumulative_effect),
               names_to = "metric",
               values_to = "value") %>%
  mutate(
    metric = recode(metric,
                    "cumulative_icer" = "ICER ($ / DALY averted)",
                    "cumulative_cost" = "Incremental Cost (USD)",
                    "cumulative_effect" = "Incremental DALYs Averted"),
    final_value = case_when(
      metric == "ICER ($ / DALY averted)" ~ final_icer,
      metric == "Incremental Cost (USD)" ~ final_cost,
      metric == "Incremental DALYs Averted" ~ final_effect
    )
  )

# Create faceted plot
faceted_plot <- ggplot(faceted_data, aes(x = n_runs, y = value)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(color = "steelblue", size = 1.5, alpha = 0.6) +
  geom_hline(aes(yintercept = final_value), 
             linetype = "dashed", color = "red", linewidth = 0.8) +
  facet_wrap(~metric, scales = "free_y", ncol = 1) +
  labs(
    title = "Cumulative ICER Stabilization Analysis",
    subtitle = paste0("Cumulative means at each sample size | Final values shown as red dashed lines | N = ", total_runs, " PSA runs"),
    x = "Number of PSA Runs (cumulative)",
    y = "Cumulative Mean Value"
  ) +
  scale_x_continuous(labels = comma, breaks = seq(0, total_runs, by = 200)) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 10, color = "gray40"),
    strip.text = element_text(face = "bold", size = 11),
    panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
    panel.grid.major = element_line(color = "gray85"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Display the faceted plot
print(faceted_plot)

# Save the plot
ggsave("icer_stabilization_faceted_plot.png", plot = faceted_plot, width = 10, height = 12)
ggsave("icer_stabilization_faceted_plot.pdf", plot = faceted_plot, width = 10, height = 12)

15 Summary Table of ICER Stabilization at Key Milestones

# --------------------------------
# Summary Table of Key Milestones
# --------------------------------

milestones <- c(1, 10, 25, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, total_runs)
milestones <- milestones[milestones <= total_runs]

summary_table <- data.frame(
  PSA_Runs = milestones,
  Cumulative_ICER = sapply(milestones, function(n) 
    dollar(round(mean(incremental$inc_cost[1:n]) / mean(incremental$dalys_averted[1:n]), 0))),
  Cumulative_Cost = sapply(milestones, function(n) 
    dollar(round(mean(incremental$inc_cost[1:n]), 0))),
  Cumulative_DALYs = sapply(milestones, function(n) 
    round(mean(incremental$dalys_averted[1:n]), 4))
)

cat("\n\n========================================\n")
## 
## 
## ========================================
cat("CUMULATIVE VALUES AT KEY MILESTONES\n")
## CUMULATIVE VALUES AT KEY MILESTONES
cat("========================================\n")
## ========================================
print(summary_table)
##    PSA_Runs Cumulative_ICER Cumulative_Cost Cumulative_DALYs
## 1         1          $3,252            $367           0.1130
## 2        10          $3,917            $452           0.1154
## 3        25          $3,921            $456           0.1163
## 4        50          $4,017            $463           0.1153
## 5       100          $3,914            $456           0.1166
## 6       200          $3,871            $451           0.1164
## 7       300          $3,904            $454           0.1163
## 8       400          $3,949            $459           0.1161
## 9       500          $3,935            $457           0.1161
## 10      600          $3,960            $459           0.1159
## 11      700          $3,960            $459           0.1160
## 12      800          $3,941            $457           0.1159
## 13      900          $3,964            $460           0.1160
## 14     1000          $3,977            $461           0.1158
# Convergence assessment
last_20_start <- floor(total_runs * 0.8)
last_20_icer <- sapply(seq(last_20_start, total_runs, by = 10), function(n) {
  mean(incremental$inc_cost[1:n]) / mean(incremental$dalys_averted[1:n])
})
relative_range <- (max(last_20_icer) - min(last_20_icer)) / final_icer * 100

cat("\n========================================\n")
## 
## ========================================
cat("CONVERGENCE ASSESSMENT\n")
## CONVERGENCE ASSESSMENT
cat("========================================\n")
## ========================================
cat("Relative range in last 20% of runs:", sprintf("%.2f%%", relative_range), "\n")
## Relative range in last 20% of runs: 0.89%
if(relative_range < 5) {
  cat("âś“ EXCELLENT: ICER well-stabilized\n")
} else if(relative_range < 10) {
  cat("âś“ GOOD: ICER reasonably stabilized\n")
} else {
  cat("âš  CAUTION: Consider more PSA runs\n")
}
## âś“ EXCELLENT: ICER well-stabilized

16 Assess the point of convergence

# Based on your data, let's analyze how many runs are sufficient
total_runs <- nrow(incremental)

cat("========================================\n")
## ========================================
cat("PSA SAMPLE SIZE DETERMINATION\n")
## PSA SAMPLE SIZE DETERMINATION
cat("========================================\n\n")
## ========================================
# --------------------------------
# 1. Analyze Convergence of Your Current Data
# --------------------------------

# Calculate cumulative ICER at different intervals
intervals <- c(50, 100, 200, 300, 400, 500, 600, 700, 800, 900, total_runs)
cumulative_icers <- sapply(intervals, function(n) {
  mean(incremental$inc_cost[1:n]) / mean(incremental$dalys_averted[1:n])
})

# Calculate stabilization metrics
final_icer <- tail(cumulative_icers, 1)
relative_diff <- abs(cumulative_icers - final_icer) / final_icer * 100

# Create convergence data
convergence_check <- data.frame(
  n_runs = intervals,
  cumulative_icer = cumulative_icers,
  diff_from_final_pct = relative_diff
)

cat("CONVERGENCE ANALYSIS:\n")
## CONVERGENCE ANALYSIS:
print(convergence_check)
##    n_runs cumulative_icer diff_from_final_pct
## 1      50        4017.026           1.0156278
## 2     100        3913.542           1.5866633
## 3     200        3871.459           2.6449246
## 4     300        3904.176           1.8221824
## 5     400        3948.753           0.7012224
## 6     500        3935.398           1.0370605
## 7     600        3959.528           0.4302540
## 8     700        3959.538           0.4300028
## 9     800        3941.065           0.8945502
## 10    900        3964.330           0.3095017
## 11   1000        3976.638           0.0000000
cat("\n")
# Determine when ICER stabilized (within 5% of final)
stabilized_at <- convergence_check$n_runs[which(convergence_check$diff_from_final_pct < 5)][1]
if(!is.na(stabilized_at)) {
  cat("âś“ ICER stabilized within 5% after approximately", stabilized_at, "runs\n")
} else {
  cat("âš  ICER has not yet stabilized within 5% after", total_runs, "runs\n")
}
## âś“ ICER stabilized within 5% after approximately 50 runs

17 Use the coefficient of variation to assess stability

# 2. Coefficient of Variation (CV) Analysis
# Calculate incremental costs and effects for all runs
incremental_costs <- incremental$inc_cost
incremental_effects <- incremental$dalys_averted
# Calculate ICER for each run
icers <- incremental_costs / incremental_effects
# Calculate CV for incremental costs, effects, and ICERs
cv_costs <- sd(incremental_costs, na.rm = TRUE) / mean(incremental_costs, na.rm = TRUE)
cv_effects <- sd(incremental_effects, na.rm = TRUE) / mean(incremental_effects, na.rm = TRUE)
cv_icers <- sd(icers, na.rm = TRUE) / mean(icers, na.rm = TRUE)
cat("COEFFICIENT OF VARIATION (CV) ANALYSIS:\n")
## COEFFICIENT OF VARIATION (CV) ANALYSIS:
cat("CV of Incremental Costs:", sprintf("%.2f%%", cv_costs * 100), "\n")
## CV of Incremental Costs: 21.59%
cat("CV of Incremental Effects:", sprintf("%.2f%%", cv_effects * 100), "\n")
## CV of Incremental Effects: 7.77%
cat("CV of ICERs:", sprintf("%.2f%%", cv_icers * 100), "\n")
## CV of ICERs: 22.73%
if(cv_icers < 0.5) {
  cat("âś“ LOW VARIABILITY: ICER estimates are stable\n")
} else if(cv_icers < 1) {
  cat("âś“ MODERATE VARIABILITY: Consider more runs for stability\n")
} else {
  cat("âš  HIGH VARIABILITY: More PSA runs likely needed\n")
}
## âś“ LOW VARIABILITY: ICER estimates are stable

Plot the cumulative ICER with stabilization bands

# 3. Plot Cumulative ICER with Stabilization Bands
cumulative_data <- data.frame(
  n_runs = seq(1, total_runs),
  cumulative_icer = sapply(seq(1, total_runs), function(n) {
    mean(incremental$inc_cost[1:n]) / mean(incremental$dalys_averted[1:n])
  })
)    
final_icer <- tail(cumulative_data$cumulative_icer, 1)                            
band_5pct <- abs(final_icer * 0.05)                            
band_10pct <- abs(final_icer * 0.10)                         
cumulative_plot <- ggplot(cumulative_data, aes(x = n_runs, y = cumulative_icer)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(color = "steelblue", size = 1.5, alpha = 0.6) +
  geom_hline(yintercept = final_icer, linetype = "dashed", color = "red", linewidth = 0.8) +
  geom_hline(yintercept = final_icer + band_5pct, linetype = "dotted", color = "orange", linewidth = 0.8) +
  geom_hline(yintercept = final_icer - band_5pct, linetype = "dotted", color = "orange", linewidth = 0.8) +
  geom_hline(yintercept = final_icer + band_10pct, linetype = "dotted", color = "blue", linewidth = 0.8) +
  geom_hline(yintercept = final_icer - band_10pct, linetype = "dotted", color = "blue", linewidth = 0.8) +
  labs(
    title = "Cumulative ICER Stabilization Analysis",
    subtitle = paste0("Final ICER: ", dollar(final_icer), "/DALY averted | Bands at ±5% (orange) and ±10% (blue)"),
    x = "Number of PSA Runs (cumulative)",
    y = "Cumulative Mean ICER ($ / DALY averted)"
  ) +
  scale_x_continuous(labels = comma, breaks = seq(0, total_runs, by = 200)) +
  scale_y_continuous(labels = dollar_format()) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 10, color = "gray40"),
    panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
    panel.grid.major = element_line(color = "gray85"),
    axis.text.x = element_text(angle = 45, hjust = 1)                               
  )                            
print(cumulative_plot)                             

# Save the plot                             
#ggsave("cumulative_icer_stabilization_plot.png", plot = cumulative_plot, width = 10, height = 6)                            
#ggsave("cumulative_icer_stabilization_plot.pdf", plot = cumulative_plot, width = 10, height = 6)

Plot the cumulative incremental cost and effect with stabilization bands-Combine with the ICER plot

# 4. Plot Cumulative Incremental Cost and Effect with Stabilization Bands
cumulative_data <- data.frame(
  n_runs = seq(1, total_runs),
  cumulative_icer = sapply(seq(1, total_runs), function(n) {
    mean(incremental$inc_cost[1:n]) / mean(incremental$dalys_averted[1:n])
  }),
  cumulative_cost = sapply(seq(1, total_runs), function(n) {
    mean(incremental$inc_cost[1:n])
  }),
  cumulative_effect = sapply(seq(1, total_runs), function(n) {
    mean(incremental$dalys_averted[1:n])
  })
)

final_icer <- tail(cumulative_data$cumulative_icer, 1)
final_cost <- tail(cumulative_data$cumulative_cost, 1)
final_effect <- tail(cumulative_data$cumulative_effect, 1)

band_5pct_icer <- abs(final_icer * 0.05)
band_5pct_cost <- abs(final_cost * 0.05)
band_5pct_effect <- abs(final_effect * 0.05)

# Reshape data for faceted plot
faceted_data <- cumulative_data %>%
  select(n_runs, cumulative_icer, cumulative_cost, cumulative_effect) %>%
  pivot_longer(cols = c(cumulative_icer, cumulative_cost, cumulative_effect),
               names_to = "metric",
               values_to = "value") %>%
  mutate(
    metric = recode(metric,
                    "cumulative_icer" = "ICER ($ / DALY averted)",
                    "cumulative_cost" = "Incremental Cost (USD)",
                    "cumulative_effect" = "Incremental DALYs Averted"),
    final_value = case_when(
      metric == "ICER ($ / DALY averted)" ~ final_icer,
      metric == "Incremental Cost (USD)" ~ final_cost,
      metric == "Incremental DALYs Averted" ~ final_effect
    ),
    band_5pct = case_when(
      metric == "ICER ($ / DALY averted)" ~ band_5pct_icer,
      metric == "Incremental Cost (USD)" ~ band_5pct_cost,
      metric == "Incremental DALYs Averted" ~ band_5pct_effect
    )
  )

# Create faceted plot
faceted_plot <- ggplot(faceted_data, aes(x = n_runs, y = value)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(color = "steelblue", size = 1.5, alpha = 0.6) +
  geom_hline(aes(yintercept = final_value), linetype = "dashed", color = "red", linewidth = 0.8) +
  geom_hline(aes(yintercept = final_value + band_5pct), linetype = "dotted", color = "blue", linewidth = 0.8) +
  geom_hline(aes(yintercept = final_value - band_5pct), linetype = "dotted", color = "blue", linewidth = 0.8) +
  labs(
    title = "Cumulative ICER, Incremental Cost, and Incremental Effect Stabilization Analysis",
    subtitle = paste0("Final values (At 1000 PSA runs) shown as red dashed lines | 5% stabilization bands shown in blue"),
    x = "Number of PSA Runs (cumulative)",
    y = "Cumulative Mean Value"
  ) +
  facet_wrap(~metric, scales = "free_y", ncol = 1) +
  scale_x_continuous(labels = comma, breaks = seq(0, total_runs, by = 200)) +
  scale_y_continuous(labels = dollar_format()) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 10, color = "gray40"),
    strip.text = element_text(face = "bold", size = 11),
    panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
    panel.grid.major = element_line(color = "gray85"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

print(faceted_plot)

# Save the plot
ggsave("cumulative_icer_cost_effect_stabilization_plot_psa.png", plot = faceted_plot, width = 10, height = 12)
ggsave("cumulative_icer_cost_effect_stabilization_plot_psa.pdf", plot = faceted_plot, width = 10, height = 12)