Probabilistic Sensitivity Analysis Results–Kenya CVD Microsimulation Model
James Oguta
December 1, 2025
- 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 Calculate Incremental Cost and Effect
- 8 Cost-Effectiveness Plane (CE Plane) Plot
- 9 Draw CEAC and save as image
- 10 Add no intervention options
- 11 CE Plane for all the three interventions
- 12 CE-Plane for only Empower Health vs Usual Care with WTP line
- 13 CE-Plane for only Empower Health vs Usual Care with WTP line adding the mean ICER point in different color
- 14 Cost-Effectiveness Acceptability Curve (All Interventions)
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:
## ── 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 4.0.0 ✔ 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##
## 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.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_Emp_Health.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/psa_results_Usual_Care.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 1206.473 4.743405
## 2 2 1228.786 6.133804
## 3 3 1187.193 6.112661
## 4 4 1198.016 6.148143
## 5 5 1129.835 6.776182
## 6 6 1203.293 6.188534
## mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1 7176.782 4.645208 6217.620
## 2 6771.298 6.168652 5904.517
## 3 6737.978 6.144695 5852.409
## 4 6741.760 6.175688 5866.236
## 5 6496.725 6.820126 5630.660
## 6 6726.379 6.225300 5850.450
## mean_Ddalys.Usual_Care
## 1 4.889149
## 2 6.375088
## 3 6.358847
## 4 6.386532
## 5 7.034411
## 6 6.4365487 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.
## 824.9 863.9 893.2 896.3 928.9 988.2## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1949 0.2090 0.2189 0.2203 0.2319 0.2535## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3899 4010 4073 4073 4137 42428 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()
9 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
## 210 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)
11 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
## 212 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
## 213 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
mean_delta_cost <- mean(psa_wide$delta_cost, na.rm = TRUE)
mean_delta_dalys <- mean(psa_wide$delta_dalys, na.rm = TRUE)
psa_ce_plane <- 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 = NULL,
x = "Incremental DALYs Averted",
y = "Incremental Costs (USD)",
color = "ICERs for each PSA run"
) +
theme_minimal()
psa_ce_plane## 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.
ggsave("CE_Plane_Empower_Health_vs_Usual_Care_Mean_ICER.pdf", plot = psa_ce_plane, width = 8, height = 6)## 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.14 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
## Summary: Empower Health vs No Intervention## delta_cost_EmpNoInt delta_dalys_EmpNoInt icer_EmpNoInt
## Min. :5272 Min. :-0.068524 Min. :-5153245
## 1st Qu.:5455 1st Qu.:-0.031808 1st Qu.: -175719
## Median :5630 Median :-0.005172 Median : -87703
## Mean :5647 Mean : 0.017057 Mean : -59057
## 3rd Qu.:5842 3rd Qu.: 0.069452 3rd Qu.: 78378
## Max. :6068 Max. : 0.123593 Max. :18722490##
## Summary: Usual Care vs No Intervention## delta_cost_UCNoInt delta_dalys_UCNoInt icer_UCNoInt
## Min. :4431 Min. :-0.2685 Min. :-40412
## 1st Qu.:4592 1st Qu.:-0.2418 1st Qu.:-29934
## Median :4738 Median :-0.2257 Median :-20973
## Mean :4751 Mean :-0.2032 Mean :-24820
## 3rd Qu.:4911 3rd Qu.:-0.1624 3rd Qu.:-18970
## Max. :5114 Max. :-0.1264 Max. :-16740## Summary: Empower Health vs Usual Care## delta_cost delta_dalys icer
## Min. :824.9 Min. :0.1949 Min. :3899
## 1st Qu.:863.9 1st Qu.:0.2090 1st Qu.:4010
## Median :893.2 Median :0.2189 Median :4073
## Mean :896.3 Mean :0.2203 Mean :4073
## 3rd Qu.:928.9 3rd Qu.:0.2319 3rd Qu.:4137
## Max. :988.2 Max. :0.2535 Max. :4242