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)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── 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

2 Import PSA Datasets

# Import PSA datasets
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/res_no_intervention_parallel.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/res_Empower_Health_parallel.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/res_Usual_Care_parallel.rda')

3 Rename datasets

subset_psa_no_int <- res_no_intervention_parallel[, c("mean_Dcosts", "mean_Ddalys")]
subset_psa_Emp_Health <- res_Empower_Health_parallel[, c("mean_Dcosts", "mean_Ddalys")]
subset_psa_UC <- res_Usual_Care_parallel[, 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                    2209.915                    4.930132
## 2   2                    2127.673                    6.324136
## 3   3                    2090.577                    6.306282
## 4   4                    2099.303                    6.335162
## 5   5                    2002.332                    6.967782
## 6   6                    2103.297                    6.386979
##   mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1                   6058.030                   4.555382               6059.929
## 2                   5722.752                   6.002305               5742.818
## 3                   5702.600                   5.971203               5712.601
## 4                   5702.067                   6.008426               5715.055
## 5                   5507.058                   6.644298               5517.995
## 6                   5695.474                   6.060224               5711.839
##   mean_Ddalys.Usual_Care
## 1               4.681237
## 2               6.109798
## 3               6.084559
## 4               6.117239
## 5               6.750024
## 6               6.166955

7 Calculate Incremental Cost and Effect

# Calculate incremental cost and incremental DALYs (Usual_Care as comparator)
psa_wide$delta_cost <- psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.Usual_Care
psa_wide$delta_dalys <- psa_wide$mean_Ddalys.Usual_Care - psa_wide$mean_Ddalys.Empower_Health
# Avoid division by zero in ICER
psa_wide$icer <- with(psa_wide, ifelse(delta_dalys == 0, NA, delta_cost / delta_dalys))

# Summary of incremental costs and DALYs
summary(psa_wide$delta_cost)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -34.235 -17.443 -12.882 -13.496  -8.948   6.163
summary(psa_wide$delta_dalys)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1006  0.1065  0.1098  0.1112  0.1164  0.1263
# Summarize ICER
summary(psa_wide$icer)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -340.19 -161.65 -116.28 -123.97  -78.28   49.95

8 Cost-Effectiveness Plane (CE Plane) Plot

# Clean NA from icers for plotting points
valid <- !is.na(psa_wide$delta_dalys) & !is.na(psa_wide$delta_cost)

plot(psa_wide$delta_dalys[valid], psa_wide$delta_cost[valid],
     xlab = "Incremental DALYs Averted (Usual Care vs Empower Health)",
     ylab = "Incremental Cost (USD)",
     main = "Cost-Effectiveness Plane",
     pch = 19, col = rgb(0, 0, 1, 0.3))

abline(h = 0, col = "gray", lty = 2)  # horizontal zero cost line
abline(v = 0, col = "gray", lty = 2)  # vertical zero effect line

# Histogram of ICER

hist(psa_wide$icer,
     breaks = 50,
     main = "Distribution of ICERs (Empower Health vs Usual Care)",
     xlab = "ICER (USD per DALY averted)",
     col = "lightblue",
     xlim = c(min(psa_wide$icer, na.rm=TRUE), quantile(psa_wide$icer, 0.95, na.rm=TRUE)))

# Cost-Effectiveness Acceptability Curve (CEAC)

# Define range of WTP thresholds, e.g., 0 to 5000 USD per DALY averted
wtp_values <- seq(0, 5000, by = 50)

# Initialize vector to store probability cost-effective at each WTP
ceac <- numeric(length(wtp_values))

for(i in seq_along(wtp_values)) {
  wtp <- wtp_values[i]
  # Cost-effective if ICER <= WTP
  ceac[i] <- mean(psa_wide$icer <= wtp, na.rm = TRUE)
}

# Plot CEAC
plot(wtp_values, ceac,
     type = "l",
     lwd = 2,
     col = "darkgreen",
     xlab = "Willingness-to-pay Threshold (USD per DALY averted)",
     ylab = "Probability Cost-Effective",
     main = "Cost-Effectiveness Acceptability Curve")
grid()

# Add no intervention options

# Empower_Health vs No_Intervention
psa_wide$delta_cost_EmpNoInt <- psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.No_Intervention
psa_wide$delta_dalys_EmpNoInt <- psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Empower_Health
psa_wide$icer_EmpNoInt <- psa_wide$delta_cost_EmpNoInt / psa_wide$delta_dalys_EmpNoInt

# Usual_Care vs No_Intervention
psa_wide$delta_cost_UCNoInt <- psa_wide$mean_Dcosts.Usual_Care - psa_wide$mean_Dcosts.No_Intervention
psa_wide$delta_dalys_UCNoInt <- psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Usual_Care
psa_wide$icer_UCNoInt <- psa_wide$delta_cost_UCNoInt / psa_wide$delta_dalys_UCNoInt

par(mfrow = c(1,2))

# Empower Health vs No Intervention
plot(psa_wide$delta_dalys_EmpNoInt, psa_wide$delta_cost_EmpNoInt,
     xlab = "Incremental DALYs Averted (No Intervention vs Empower Health)",
     ylab = "Incremental Cost (USD)",
     main = "CE Plane: Empower Health vs No Intervention",
     pch = 19, col = rgb(1, 0, 0, 0.3))
abline(h = 0, v = 0, col = "gray", lty = 2)

# Usual Care vs No Intervention
plot(psa_wide$delta_dalys_UCNoInt, psa_wide$delta_cost_UCNoInt,
     xlab = "Incremental DALYs Averted (No Intervention vs Usual Care)",
     ylab = "Incremental Cost (USD)",
     main = "CE Plane: Usual Care vs No Intervention",
     pch = 19, col = rgb(0, 0, 1, 0.3))
abline(h = 0, v = 0, col = "gray", lty = 2)

par(mfrow = c(1,1))

9 CE Plane for all the three interventions

# Calculate incremental values including Empower vs Usual
cep_data <- data.frame(
  sim = psa_wide$sim,
  delta_cost_emp = psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.No_Intervention,
  delta_dalys_emp = psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Empower_Health,
  
  delta_cost_uc = psa_wide$mean_Dcosts.Usual_Care - psa_wide$mean_Dcosts.No_Intervention,
  delta_dalys_uc = psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Usual_Care,
  
  delta_cost_emp_vs_uc = psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.Usual_Care,
  delta_dalys_emp_vs_uc = psa_wide$mean_Ddalys.Usual_Care - psa_wide$mean_Ddalys.Empower_Health
)

# Convert to long format for ggplot
cep_long <- bind_rows(
  data.frame(
    sim = cep_data$sim,
    delta_cost = cep_data$delta_cost_emp,
    delta_dalys = cep_data$delta_dalys_emp,
    comparison = "Empower Health vs. No Intervention"
  ),
  data.frame(
    sim = cep_data$sim,
    delta_cost = cep_data$delta_cost_uc,
    delta_dalys = cep_data$delta_dalys_uc,
    comparison = "Usual Care vs. No Intervention"
  ),
  data.frame(
    sim = cep_data$sim,
    delta_cost = cep_data$delta_cost_emp_vs_uc,
    delta_dalys = cep_data$delta_dalys_emp_vs_uc,
    comparison = "Empower Health vs. Usual Care"
  )
)


# Plot CEP with WTP threshold of $1000/DALY
ggplot(cep_long, aes(x = delta_dalys, y = delta_cost, color = comparison)) +
  geom_point(alpha = 0.4, size = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  # Add WTP line
  geom_abline(slope = 1000, intercept = 0, linetype = "dotted", color = "black", size = 1) +
  annotate("text", x = max(cep_long$delta_dalys, na.rm = TRUE) * 0.9,
           y = 1000 * max(cep_long$delta_dalys, na.rm = TRUE) * 0.9,
           label = "WTP = $1000/DALY", hjust = 1, vjust = -0.5, angle = atan(1000)*180/pi, size = 3.5) +
  labs(
    title = "Cost-Effectiveness Plane at US$1000 WTP Threshold",
    x = "Incremental DALYs Averted",
    y = "Incremental Costs (USD)",
    color = "Comparison"
  ) +
  theme_minimal()
## 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.

10 Cost-Effectiveness Acceptability Curve (All Interventions)

# Step 1: Define WTP thresholds
wtp_values <- seq(0, 20000, by = 100)

# Step 2: Pivot long to wide
library(reshape2)
psa_wide_costs <- dcast(psa_results, sim ~ strategy, value.var = "mean_Dcosts")
psa_wide_dalys <- dcast(psa_results, sim ~ strategy, value.var = "mean_Ddalys")

# Step 3: Calculate NMB for each strategy at each WTP
strategies <- c("No_Intervention", "Usual_Care", "Empower_Health")
ceac_matrix <- matrix(NA, nrow = length(wtp_values), ncol = length(strategies))
colnames(ceac_matrix) <- strategies

for (i in seq_along(wtp_values)) {
  wtp <- wtp_values[i]
  
  # Calculate NMB: NMB = DALYs averted × WTP − Cost
  nmb_df <- data.frame(
    No_Intervention = (max(psa_wide_dalys) - psa_wide_dalys$No_Intervention) * wtp - psa_wide_costs$No_Intervention,
    Usual_Care      = (max(psa_wide_dalys) - psa_wide_dalys$Usual_Care) * wtp - psa_wide_costs$Usual_Care,
    Empower_Health  = (max(psa_wide_dalys) - psa_wide_dalys$Empower_Health) * wtp - psa_wide_costs$Empower_Health
  )
  
  # Find which strategy has the max NMB per simulation
  best_strategy <- apply(nmb_df, 1, function(x) colnames(nmb_df)[which.max(x)])
  
  # Compute CEAC (probability each strategy is cost-effective)
  ceac_matrix[i, ] <- table(factor(best_strategy, levels = strategies)) / nrow(nmb_df)
}

# Step 4: Plot CEAC using matplot (base R)
matplot(wtp_values, ceac_matrix, type = "l", lty = 1, lwd = 2,
        col = c("black", "red", "blue"),
        xlab = "Willingness-to-pay Threshold (USD)",
        ylab = "Probability Cost-Effective",
        main = "Cost-Effectiveness Acceptability Curve")

legend("right", legend = strategies, col = c("black", "red", "blue"), lty = 1, lwd = 2)

# Summary of ICERs for all the intervention pairs

cat("Summary: Empower Health vs No Intervention\n")
## Summary: Empower Health vs No Intervention
summary(psa_wide[, c("delta_cost_EmpNoInt", "delta_dalys_EmpNoInt", "icer_EmpNoInt")])
##  delta_cost_EmpNoInt delta_dalys_EmpNoInt icer_EmpNoInt  
##  Min.   :3420        Min.   :0.3012       Min.   :10084  
##  1st Qu.:3540        1st Qu.:0.3219       1st Qu.:10671  
##  Median :3650        Median :0.3354       Median :10866  
##  Mean   :3655        Mean   :0.3377       Mean   :10834  
##  3rd Qu.:3771        3rd Qu.:0.3533       3rd Qu.:10998  
##  Max.   :3911        Max.   :0.3879       Max.   :11489
cat("\nSummary: Usual Care vs No Intervention\n")
## 
## Summary: Usual Care vs No Intervention
summary(psa_wide[, c("delta_cost_UCNoInt", "delta_dalys_UCNoInt", "icer_UCNoInt")])
##  delta_cost_UCNoInt delta_dalys_UCNoInt  icer_UCNoInt  
##  Min.   :3450       Min.   :0.1996      Min.   :14933  
##  1st Qu.:3557       1st Qu.:0.2156      1st Qu.:15881  
##  Median :3664       Median :0.2256      Median :16208  
##  Mean   :3668       Mean   :0.2266      Mean   :16214  
##  3rd Qu.:3777       3rd Qu.:0.2374      3rd Qu.:16512  
##  Max.   :3907       Max.   :0.2616      Max.   :17426
cat("Summary: Empower Health vs Usual Care\n")
## Summary: Empower Health vs Usual Care
summary(psa_wide[, c("delta_cost", "delta_dalys", "icer")])
##    delta_cost       delta_dalys          icer        
##  Min.   :-34.235   Min.   :0.1006   Min.   :-340.19  
##  1st Qu.:-17.443   1st Qu.:0.1065   1st Qu.:-161.65  
##  Median :-12.882   Median :0.1098   Median :-116.28  
##  Mean   :-13.496   Mean   :0.1112   Mean   :-123.97  
##  3rd Qu.: -8.948   3rd Qu.:0.1164   3rd Qu.: -78.28  
##  Max.   :  6.163   Max.   :0.1263   Max.   :  49.95