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/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_no_trt_2.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_Emp_Health_2.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_Usual_Care_2.rda')

3 Rename datasets

subset_psa_no_int <- psa_results_no_trt_2[, c("mean_Dcosts", "mean_Ddalys")]
subset_psa_Emp_Health <- psa_results_Emp_Health_2[, c("mean_Dcosts", "mean_Ddalys")]
subset_psa_UC <- psa_results_Usual_Care_2[, 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 <- reshape(psa_results, 
                timevar = "strategy", 
                idvar = "sim", 
                direction = "wide")
head(psa_wide)
##   sim mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
## 1   1                    967.3705                    7.399176
## 2   2                    884.9262                    7.435233
## 3   3                    907.1375                    7.424413
## 4   4                   1007.9599                    7.382568
## 5   5                    973.2876                    7.397506
## 6   6                    951.3212                    7.405980
##   mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1                   4851.712                   7.298338               4283.573
## 2                   4783.454                   7.320417               4188.935
## 3                   4797.463                   7.315476               4213.083
## 4                   4887.848                   7.287305               4333.161
## 5                   4865.832                   7.293878               4297.776
## 6                   4842.831                   7.301330               4268.467
##   mean_Ddalys.Usual_Care
## 1               7.409644
## 2               7.441969
## 3               7.432341
## 4               7.394076
## 5               7.405577
## 6               7.414809
# Save the psa_wide dataset as an .rda file
save(psa_wide, file = "/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_wide.rda")
# Save the psa_wide dataset as a .csv file
write.csv(psa_wide, file = "/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_1st.csv", row.names = FALSE)

7 Load CSV File

data <- read.csv("/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_1st.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  967 885 907 1008 973 ...
##  $ mean_Ddalys.No_Intervention: num  7.4 7.44 7.42 7.38 7.4 ...
##  $ mean_Dcosts.Empower_Health : num  4852 4783 4797 4888 4866 ...
##  $ mean_Ddalys.Empower_Health : num  7.3 7.32 7.32 7.29 7.29 ...
##  $ mean_Dcosts.Usual_Care     : num  4284 4189 4213 4333 4298 ...
##  $ mean_Ddalys.Usual_Care     : num  7.41 7.44 7.43 7.39 7.41 ...

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          4848.937        4278.653           7.299519         7.411795
##   sd_cost_empower sd_cost_usual sd_dalys_empower sd_dalys_usual
## 1        53.87085      71.59616       0.01686957     0.02313703
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.Empower_Health - mean_Ddalys.Usual_Care,

    # 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(icer, na.rm = TRUE),
    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       570.284    17.83834          534.1923          602.8531     -0.1122763
##   sd_inc_dalys ci_inc_dalys_lower ci_inc_dalys_upper mean_dalys_averted
## 1  0.006324055         -0.1243319        -0.09978764          0.1122763
##   sd_dalys_averted ci_dalys_averted_lower ci_dalys_averted_upper mean_icer
## 1      0.006324055             0.09978764              0.1243319  5086.519
##   median_icer icer_lower icer_upper
## 1    5080.567   4849.276   5353.795

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 = 1011, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  569.1770 571.3909
## sample estimates:
## mean difference 
##         570.284
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 = -561.43, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.1126687 -0.1118838
## sample estimates:
## mean difference 
##      -0.1122763

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"
    )
  )
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")

# 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" = "blue",
      "Empower Health" = "green"
    )
  ) +
  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)

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"
  ) +

  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"
    )
  )

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)

# 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%   5079.28    5079.209 5071.376 5087.157                               0
##      prob_cost_effective_at_1x_GDP prob_cost_effective_at_3x_GDP
## 2.5%                             0                             1

12 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 / incremental$dalys_averted, na.rm = TRUE)



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("Mean 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       $4,848.94
## 6                Mean Cost - Usual Care       $4,278.65
## 7           Mean DALYs - Empower Health          7.2995
## 8               Mean DALYs - Usual Care          7.4118
## 9                 Mean Incremental Cost         $570.28
## 10                   Mean DALYs Averted          0.1123
## 11                            Mean ICER          $5,087
## 12                          Median ICER          $5,081
## 13              95% CI ICER (Bootstrap) $5,071 - $5,087
## 14 Probability ICER < 0.5x GDP ($1,000)            0.0%
## 15   Probability ICER < 1x GDP ($2,000)            0.0%
## 16   Probability ICER < 3x GDP ($6,000)          100.0%
## 17            EVPI at 0.5x GDP ($1,000)              $0
## 18              EVPI at 1x GDP ($2,000)              $0
## 19              EVPI at 3x GDP ($6,000)              $0

13 Save all plots and summary table to files

# Save plots
# ggsave("cost_effectiveness_scatter.png", plot = ce_scatter, width = 8, height = 6)
# ggsave("icer_scatter.png", plot = icer_scatter, width = 8, height = 6)
# ggsave("ceac_plot.png", plot = ceac_plot, width = 8, height = 6)
# ggsave("nmb_plot.png", plot = nmb_plot, width = 8, height = 6)
# ggsave("evpi_plot.png", plot = evpi_plot, width = 8, height = 6)
# ggsave("ce_plane.png", plot = ce_plane, width = 8, height = 6)