Probabilistic Sensitivity Analysis Results–Kenya CVD Microsimulation Model
James Oguta
March 17, 2026
- 0.1 Code Description
- 1 Clear the workspace
- 2 Import PSA Datasets
- 3 Rename datasets
- 4 Add strategy to the intervention datasets
- 5 Combine into one dataset for plotting
- 6 Reshape data to wide format
- 7 Summarise the mean costs and DALYs for each intervention
- 8 Plot CE Plane for costs and DALYs for all interventions
- 9 Calculate Incremental Cost and Effect
- 10 Cost-Effectiveness Plane (CE Plane) Plot
- 11 Draw CEAC and save as image
- 12 Add no intervention options
- 13 CE Plane for all the three interventions
- 14 CE-Plane for only Empower Health vs Usual Care with WTP line
- 15 CE-Plane for only Empower Health vs Usual Care with WTP line adding the mean ICER point in different color
- 16 The above CE-Plane with different WTP lines (1000, 2000, 6000 USD)
- 17 Cost-Effectiveness Acceptability Curve (All Interventions)
- 18 Summary of ICERs for all the intervention pairs
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
1.1 Load Libraries
The script begins by loading the necessary R libraries:
## 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##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smiths2 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
4 Add strategy to the intervention datasets
5 Combine into one dataset for plotting
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.4148097 Summarise the mean costs and DALYs for each intervention
## mean_Dcosts mean_Ddalys strategy sim
## Min. : 743.0 Min. :7.312 Length:1000 Min. : 1.0
## 1st Qu.: 920.1 1st Qu.:7.386 Class :character 1st Qu.: 250.8
## Median : 959.8 Median :7.402 Mode :character Median : 500.5
## Mean : 960.8 Mean :7.403 Mean : 500.5
## 3rd Qu.:1000.0 3rd Qu.:7.420 3rd Qu.: 750.2
## Max. :1197.0 Max. :7.501 Max. :1000.0## mean_Dcosts mean_Ddalys strategy sim
## Min. :4645 Min. :7.238 Length:1000 Min. : 1.0
## 1st Qu.:4812 1st Qu.:7.288 Class :character 1st Qu.: 250.8
## Median :4849 Median :7.299 Mode :character Median : 500.5
## Mean :4849 Mean :7.300 Mean : 500.5
## 3rd Qu.:4885 3rd Qu.:7.311 3rd Qu.: 750.2
## Max. :5062 Max. :7.368 Max. :1000.0## mean_Dcosts mean_Ddalys strategy sim
## Min. :4008 Min. :7.326 Length:1000 Min. : 1.0
## 1st Qu.:4229 1st Qu.:7.396 Class :character 1st Qu.: 250.8
## Median :4279 Median :7.411 Mode :character Median : 500.5
## Mean :4279 Mean :7.412 Mean : 500.5
## 3rd Qu.:4326 3rd Qu.:7.428 3rd Qu.: 750.2
## Max. :4566 Max. :7.504 Max. :1000.0# Print the mean costs and DALYs for each intervention
mean_costs_no_int <- mean(subset_psa_no_int$mean_Dcosts)
mean_dalys_no_int <- mean(subset_psa_no_int$mean_Ddalys)
mean_costs_emp_health <- mean(subset_psa_Emp_Health$mean_Dcosts)
mean_dalys_emp_health <- mean(subset_psa_Emp_Health$mean_Ddalys)
mean_costs_uc <- mean(subset_psa_UC$mean_Dcosts)
mean_dalys_uc <- mean(subset_psa_UC$mean_Ddalys)
cat("No Intervention - Mean Costs:", mean_costs_no_int, "Mean DALYs:", mean_dalys_no_int, "\n")## No Intervention - Mean Costs: 960.8445 Mean DALYs: 7.402547cat("Empower Health - Mean Costs:", mean_costs_emp_health, "Mean DALYs:", mean_dalys_emp_health, "\n")## Empower Health - Mean Costs: 4848.937 Mean DALYs: 7.299519## Usual Care - Mean Costs: 4278.653 Mean DALYs: 7.4117958 Plot CE Plane for costs and DALYs for all interventions
ggplot(psa_results, aes(x = mean_Ddalys, y = mean_Dcosts, color = strategy)) +
geom_point(alpha = 0.4, size = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
labs(
title = "Cost-Effectiveness Plane for All Interventions",
x = "DALYs Averted",
y = "Costs (USD)",
color = "Intervention"
) +
theme_minimal()
9 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.
## 496.1 558.3 570.5 570.3 582.6 636.5## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08807 0.10793 0.11227 0.11228 0.11649 0.13671## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4656 4998 5081 5087 5172 563310 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, 500, 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 = "black",
xlab = "Willingness-to-pay Threshold (USD per DALY averted)",
ylab = "Probability Cost-Effective",
main = "Cost-Effectiveness Acceptability Curve")
grid()
11 Draw CEAC and save as image
png(filename = "CEAC_Plot.png",
width = 800, height = 600)
plot(wtp_values, ceac,
type = "l",
lwd = 2,
col = "black",
xlab = "Willingness-to-pay Threshold (USD per DALY averted)",
ylab = "Probability Cost-Effective",
main = "Cost-Effectiveness Acceptability Curve (Empower Health vs Usual Care)")
grid()
dev.off()## quartz_off_screen
## 212 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)
13 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.
# CE- Plane for all interventions saved as image
png(filename = "CE_Plane_All_Interventions.png",
width = 800, height = 600)
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()
dev.off()## quartz_off_screen
## 214 CE-Plane for only Empower Health vs Usual Care with WTP line
png(filename = "CE_Plane_Empower_Health_vs_Usual_Care.png",
width = 800, height = 600)
ggplot(cep_long %>% filter(comparison == "Empower Health vs. Usual Care"),
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: Empower Health vs Usual Care",
x = "Incremental DALYs Averted",
y = "Incremental Costs (USD)",
color = "Comparison"
) +
theme_minimal()
dev.off()## quartz_off_screen
## 215 CE-Plane for only Empower Health vs Usual Care with WTP line adding the mean ICER point in different color
mean_delta_cost <- mean(psa_wide$delta_cost, na.rm = TRUE)
mean_delta_dalys <- mean(psa_wide$delta_dalys, na.rm = TRUE)
ce_plane_psa <- ggplot(cep_long %>% filter(comparison == "Empower Health vs. Usual Care"),
aes(x = delta_dalys, y = delta_cost, color = comparison)) +
geom_point(alpha = 0.4, size = 1.2) +
geom_point(aes(x = mean_delta_dalys, y = mean_delta_cost), color = "blue", size = 4, shape = 17) + # Mean ICER point
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 (Empower Health vs Usual Care) with Mean ICER Point",
x = "Incremental DALYs Averted",
y = "Incremental Costs (USD)",
color = "ICERs"
) +
theme_minimal()
ce_plane_psa## Warning in geom_point(aes(x = mean_delta_dalys, y = mean_delta_cost), color = "blue", : All aesthetics have length 1, but the data has 1000 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
# CE-Plane for only Empower Health vs Usual Care with WTP line adding
the mean ICER point in different color- Saved as pdf
# Calculate means from the same dataset used for plotting
plot_data <- cep_long %>% filter(comparison == "Empower Health vs. Usual Care")
mean_delta_cost <- mean(plot_data$delta_cost, na.rm = TRUE)
mean_delta_dalys <- mean(plot_data$delta_dalys, na.rm = TRUE)
max_x <- max(plot_data$delta_dalys, na.rm = TRUE)
# For WTP line angle calculation
wtp_value <- 1000 # WTP threshold
psa_ce_plane <- ggplot(plot_data,
aes(x = delta_dalys, y = delta_cost, color = comparison)) +
geom_point(alpha = 0.4, size = 1.2) +
# Mean ICER point
annotate("point",
x = mean_delta_dalys,
y = mean_delta_cost,
color = "blue", size = 4, shape = 17) +
annotate("text",
x = mean_delta_dalys,
y = mean_delta_cost,
label = "Mean ICER",
color = "blue",
vjust = -1, hjust = 0.5, size = 3) +
# Reference lines
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
# WTP line - consider if you want it through origin or through mean point
geom_abline(slope = wtp_value, intercept = 0,
linetype = "dotted", color = "black", size = 1) +
# Better WTP label placement
annotate("text",
x = max_x * 0.8,
y = wtp_value * max_x * 0.8,
label = paste0("WTP = $", wtp_value, "/DALY"),
hjust = 0, vjust = -0.5,
angle = atan(wtp_value * (max_x/mean_delta_cost)) * 180/pi,
size = 3.5) +
labs(
x = "Incremental DALYs Averted",
y = "Incremental Costs (USD)",
color = "ICERs for individual PSA runs"
) +
theme_minimal() +
theme(
legend.position = "bottom",
plot.title = element_text(hjust = 0.5)
)
psa_ce_plane
16 The above CE-Plane with different WTP lines (1000, 2000, 6000 USD)
max_x <- max(plot_data$delta_dalys, na.rm = TRUE)
max_y <- max(plot_data$delta_cost, na.rm = TRUE)
psa_ce_plane <- ggplot(plot_data,
aes(x = delta_dalys, y = delta_cost, color = comparison)) +
geom_point(alpha = 0.4, size = 1.2) +
# Mean ICER point
annotate("point",
x = mean_delta_dalys,
y = mean_delta_cost,
color = "blue", size = 4, shape = 17) +
annotate("text",
x = mean_delta_dalys,
y = mean_delta_cost,
label = "Mean ICER",
color = "blue",
vjust = -1, hjust = 0.5, size = 3) +
# Reference lines
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
# Multiple WTP lines - ALL DOTTED
geom_abline(slope = 1000, intercept = 0,
linetype = "dotted", color = "red", size = 0.8, alpha = 0.7) +
geom_abline(slope = 2000, intercept = 0,
linetype = "dotted", color = "black", size = 0.8, alpha = 0.7) +
geom_abline(slope = 6000, intercept = 0,
linetype = "dotted", color = "blue", size = 0.8, alpha = 0.7) +
# WTP labels - place them at the edge of the plot
annotate("text",
x = max_x * 0.9,
y = 1000 * max_x * 0.9,
label = "WTP = $1,000/DALY",
hjust = 1, vjust = -0.5,
color = "black",
angle = atan(1000 * (max_x/max_y)) * 180/pi,
size = 3) +
annotate("text",
x = max_x * 0.8,
y = 2000 * max_x * 0.8,
label = "WTP = $2,000/DALY",
hjust = 1, vjust = -0.5,
color = "black",
angle = atan(2000 * (max_x/max_y)) * 180/pi,
size = 3) +
annotate("text",
x = max_x * 0.7,
y = 6000 * max_x * 0.7,
label = "WTP = $6,000/DALY",
hjust = 1, vjust = -0.5,
color = "black",
angle = atan(6000 * (max_x/max_y)) * 180/pi,
size = 3) +
labs(
x = "Incremental DALYs Averted",
y = "Incremental Costs (USD)",
color = "ICERs for individual PSA runs"
) +
theme_minimal() +
theme(
legend.position = "bottom",
plot.title = element_text(hjust = 0.5)
)
psa_ce_plane
17 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)
# Draw CEAC for only Empower Health vs Usual Care and save as pdf
# CEAC for Empower Health vs Usual Care
ceac_emp_vs_uc <- data.frame(
wtp = wtp_values,
prob_cost_effective = numeric(length(wtp_values))
)
for(i in seq_along(wtp_values)) {
wtp <- wtp_values[i]
nmb_emp <- (max(psa_wide_dalys$Empower_Health) - psa_wide_dalys$Empower_Health) * wtp - psa_wide_costs$Empower_Health
nmb_uc <- (max(psa_wide_dalys$Usual_Care) - psa_wide_dalys$Usual_Care) * wtp - psa_wide_costs$Usual_Care
best_strategy <- ifelse(nmb_emp > nmb_uc, "Empower_Health", "Usual_Care")
ceac_emp_vs_uc$prob_cost_effective[i] <- mean(best_strategy == "Empower_Health")
}
# Plot CEAC
ceac_plot <- ggplot(ceac_emp_vs_uc, aes(x = wtp, y = prob_cost_effective)) +
geom_line(color = "blue", size = 1) +
labs(
title = "Cost-Effectiveness Acceptability Curve: Empower Health vs Usual Care",
x = "Willingness-to-pay Threshold (USD)",
y = "Probability Cost-Effective"
) +
theme_minimal()
ceac_plot
18 Summary of ICERs for all the intervention pairs
## Summary: Empower Health vs No Intervention## delta_cost_EmpNoInt delta_dalys_EmpNoInt icer_EmpNoInt
## Min. :3865 Min. :0.07339 Min. :29310
## 1st Qu.:3884 1st Qu.:0.09746 1st Qu.:35767
## Median :3888 Median :0.10310 Median :37707
## Mean :3888 Mean :0.10303 Mean :37957
## 3rd Qu.:3893 3rd Qu.:0.10884 3rd Qu.:39867
## Max. :3902 Max. :0.13311 Max. :52668##
## Summary: Usual Care vs No Intervention## delta_cost_UCNoInt delta_dalys_UCNoInt icer_UCNoInt
## Min. :3265 Min. :-0.014684 Min. :-907493
## 1st Qu.:3309 1st Qu.:-0.010685 1st Qu.:-422141
## Median :3318 Median :-0.009197 Median :-360416
## Mean :3318 Mean :-0.009248 Mean :-373142
## 3rd Qu.:3326 3rd Qu.:-0.007842 3rd Qu.:-311591
## Max. :3369 Max. :-0.003598 Max. :-229444## Summary: Empower Health vs Usual Care## delta_cost delta_dalys icer
## Min. :496.1 Min. :0.08807 Min. :4656
## 1st Qu.:558.3 1st Qu.:0.10793 1st Qu.:4998
## Median :570.5 Median :0.11227 Median :5081
## Mean :570.3 Mean :0.11228 Mean :5087
## 3rd Qu.:582.6 3rd Qu.:0.11649 3rd Qu.:5172
## Max. :636.5 Max. :0.13671 Max. :5633