Modelling the Cost-Effectiveness, Budget and Equity Impacts of Malaria Vaccine Introduction alongside other Malaria Prevention Strategies
Kenya, Uganda, and Burkina Faso
James Oguta- Jalorche Consulting Inc.
2026-06-16
Executive Summary
This analysis presents a Markov model for evaluating malaria intervention strategies in children under 5 years across Kenya, Uganda, and Burkina Faso. The model compares four intervention strategies across urban and rural settings.
Willingness-to-Pay (WTP) Thresholds Following WHO guidelines for cost-effectiveness thresholds in low- and middle-income countries, we use 0.5 × GDP per capita as the cost-effectiveness threshold for each country:
Country GDP per capita (USD) WTP Threshold (0.5× GDP) Kenya $2,300 $1,150 per DALY averted Uganda $1,100 $550 per DALY averted Burkina Faso $900 $450 per DALY averted
Interventions Being Compared The analysis evaluates four intervention strategies (all layered on top of ITNs as baseline prevention):
Strategy Target Population Efficacy Key Features Baseline (ITNs only) Universal 45% against transmission Current standard of care SMC + ITNs Children 2-5 years 85% against clinical malaria 4 cycles during rainy season Vaccine + ITNs Children 0-24 months 30% infection, 70% severe RTS,S/AS01 vaccine Combined Strategy All ages (layered) Synergistic effect Maximal intervention
Model Overview
Time horizon: 10 years
Cohort size: 10,000 children (scaled from birth to age 5+)
Discount rate: 3% per annum for both costs and health outcomes (DALYs)
Perspective: Healthcare system (direct medical costs + intervention delivery)
Outcomes: Clinical cases, severe cases, deaths, DALYs (discounted), costs (discounted)
Equity stratifiers: Residence (Rural/Urban) and Wealth quintile
# Clear environment to prevent conflicts from previous runs
rm(list = ls())
# Force garbage collection to free memory
gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 620141 33.2 1418608 75.8 NA 710035 38
## Vcells 1187332 9.1 8388608 64.0 49152 1963483 15# Load required packages
library(ggplot2)
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(readxl)
library(reshape2)
library(gridExtra)
library(scales)
# Set seed for reproducibility
set.seed(2026)
# Helper function for formatted numbers
format_number <- function(x) {
result <- ifelse(is.na(x) | is.infinite(x), "NA",
format(round(x, 0), big.mark = ",", scientific = FALSE))
return(result)
}
# Country-specific WTP thresholds (0.5 × GDP per capita)
country_wtp <- data.frame(
Country = c("Kenya", "Uganda", "Burkina Faso"),
GDP_per_capita = c(2300, 1100, 900),
WTP_threshold = c(1150, 550, 450),
WTP_very_CE = c(2300, 1100, 900) # 1× GDP for "very cost-effective"
)
cat("Environment setup complete\n")## Environment setup complete##
## ## Willingness-to-Pay Thresholds by Countrykable(country_wtp, caption = "Cost-Effectiveness Thresholds (0.5 × GDP per capita)",
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Country | GDP_per_capita | WTP_threshold | WTP_very_CE |
|---|---|---|---|
| Kenya | 2300 | 1150 | 2300 |
| Uganda | 1100 | 550 | 1100 |
| Burkina Faso | 900 | 450 | 900 |
## Environment setup complete## ✓ Environment setup complete## ✓ All packages loaded successfullyModel Structure
# ===================================================================
# MARKOV MODEL DIAGRAM
# ===================================================================
# Create Markov model diagram
if(!require(diagram)) {
install.packages("diagram")
library(diagram)
}
# Create a new plot
plot.new()
par(mar = c(2, 2, 3, 2), bg = "white")
# Define state positions
pos <- rbind(
c(0.5, 0.85), # Susceptible
c(0.5, 0.60), # Clinical
c(0.30, 0.35), # Severe
c(0.70, 0.35) # Death
)
# Define colors
state_colors <- c("Susceptible" = "#A8E6A3", "Clinical" = "#FFD3B6",
"Severe" = "#FF8B94", "Death" = "#AAAAAA")
transition_colors <- c("progression" = "#E41A1C", "mortality" = "#000000",
"recovery" = "#377EB8")
states <- c("Susceptible\n(Uninfected)", "Clinical Malaria\n(Non-severe)",
"Severe Malaria", "Death")
# Draw transitions
straightarrow(from = pos[1,], to = pos[2,], lwd = 2.5, arr.pos = 0.6,
arr.type = "triangle", lcol = transition_colors["progression"],
arr.col = transition_colors["progression"])
textrect(mid = c(0.55, 0.72), radx = 0.10, rady = 0.025,
lab = expression(bold(paste("Infection x", delta[clinical]))),
cex = 0.85, shadow.size = 0, box.col = "white",
col = transition_colors["progression"])
straightarrow(from = pos[2,], to = pos[3,], lwd = 2.5, arr.pos = 0.6,
arr.type = "triangle", lcol = transition_colors["progression"],
arr.col = transition_colors["progression"])
textrect(mid = c(0.35, 0.47), radx = 0.08, rady = 0.025,
lab = expression(bold(delta[severe])), cex = 0.85,
shadow.size = 0, box.col = "white",
col = transition_colors["progression"])
straightarrow(from = pos[3,], to = pos[4,], lwd = 2.5, arr.pos = 0.6,
arr.type = "triangle", lcol = transition_colors["mortality"],
arr.col = transition_colors["mortality"])
textrect(mid = c(0.5, 0.25), radx = 0.10, rady = 0.025,
lab = expression(bold(paste("Case fatality = CFR"))), cex = 0.85,
shadow.size = 0, box.col = "white",
col = transition_colors["mortality"])
curvedarrow(from = pos[2,], to = pos[1,], lwd = 2.5, arr.pos = 0.6,
curve = -0.4, lcol = transition_colors["recovery"],
arr.col = transition_colors["recovery"], arr.type = "triangle")
textrect(mid = c(0.30, 0.78), radx = 0.12, rady = 0.025,
lab = "Recovery +\nImmunity", cex = 0.80, shadow.size = 0,
box.col = "#A8E6A3", col = transition_colors["recovery"], font = 2)
curvedarrow(from = pos[3,], to = pos[1,], lwd = 2.5, arr.pos = 0.6,
curve = -0.5, lcol = transition_colors["recovery"],
arr.col = transition_colors["recovery"], arr.type = "triangle")
textrect(mid = c(0.12, 0.60), radx = 0.10, rady = 0.025,
lab = "Recovery", cex = 0.85, shadow.size = 0,
box.col = "#A8E6A3", col = transition_colors["recovery"], font = 2)
# Draw states
for(i in 1:4) {
textrect(mid = pos[i,], radx = 0.15, rady = 0.07, lab = states[i],
cex = 1.0, shadow.size = 0.01, box.col = state_colors[i],
col = "white", lwd = 3, font = 2)
}
title(main = "Markov Model Structure for Malaria Progression",
cex.main = 1.4, line = 0.8, font.main = 2)
text(x = 0.5, y = 0.96,
labels = "Monthly transitions with time-dependent discounting (3% annually)",
cex = 0.85, font = 3, col = "#555555")
Part 1: Country-Specific Parameters Kenya Parameters
# Kenya-specific parameters
kenya_params <- list(
country_name = "Kenya",
wtp_threshold = 1150, # 0.5 × $2,300 GDP per capita
geo_params = list(
Urban = list(
itn_coverage = 0.40,
vaccine_coverage = 0.79,
smc_coverage = 0.71,
eir_multiplier = 0.4,
cfr_multiplier = 1.0,
access_delay_days = 0.5,
care_seeking = 0.60,
malaria_positive = 0.13
),
Rural = list(
itn_coverage = 0.57,
vaccine_coverage = 0.62,
smc_coverage = 0.71,
eir_multiplier = 2.0,
cfr_multiplier = 1.6,
access_delay_days = 4.5,
care_seeking = 0.65,
malaria_positive = 0.28
)
),
age_params = data.frame(
Age_Group = c("0-6m", "6-24m", "2-5y", "5+"),
Population_Weight = c(0.13, 0.32, 0.55, 0.00),
Clinical_Risk = c(0.15, 0.16, 0.31, 0.04),
Severe_Risk = c(0.08, 0.06, 0.02, 0.005),
CFR_Base = c(0.037, 0.037, 0.037, 0.037)
),
base_params = list(
base_eir = 25,
smc_efficacy = 0.71,
vaccine_ve_infection = 0.38,
vaccine_ve_severe = 0.32,
itn_efficacy = 0.25,
cost_vaccine_full = 37.39,
cost_smc_per_child_per_year = 6.00,
cost_act = 1.50,
cost_outpatient = 8.00,
cost_hospital = 38.00,
cfr_base = 0.037,
dw_clinical = 0.006,
dw_severe = 0.133,
life_expectancy = 65,
discount_rate = 0.03
),
geo_pop_weights = c(Urban = 0.30, Rural = 0.70)
)
# Display Kenya parameters
kable(kenya_params$age_params, caption = "Kenya: Age-Specific Parameters",
digits = 3, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Age_Group | Population_Weight | Clinical_Risk | Severe_Risk | CFR_Base |
|---|---|---|---|---|
| 0-6m | 0.13 | 0.15 | 0.080 | 0.037 |
| 6-24m | 0.32 | 0.16 | 0.060 | 0.037 |
| 2-5y | 0.55 | 0.31 | 0.020 | 0.037 |
| 5+ | 0.00 | 0.04 | 0.005 | 0.037 |
Uganda Parameters
# Uganda-specific parameters
uganda_params <- list(
country_name = "Uganda",
wtp_threshold = 550, # 0.5 × $1,100 GDP per capita
geo_params = list(
Urban = list(
itn_coverage = 0.74,
vaccine_coverage = 0.78,
smc_coverage = 0.92,
eir_multiplier = 0.4,
cfr_multiplier = 1.0,
access_delay_days = 0.5,
care_seeking = 0.86,
malaria_positive = 0.51
),
Rural = list(
itn_coverage = 0.71,
vaccine_coverage = 0.78,
smc_coverage = 0.92,
eir_multiplier = 2.0,
cfr_multiplier = 1.6,
access_delay_days = 4.5,
care_seeking = 0.83,
malaria_positive = 0.54
)
),
age_params = data.frame(
Age_Group = c("0-6m", "6-24m", "2-5y", "5+"),
Population_Weight = c(0.14, 0.30, 0.56, 0.00),
Clinical_Risk = c(0.40, 0.57, 0.53, 0.10),
Severe_Risk = c(0.08, 0.06, 0.02, 0.005),
CFR_Base = c(0.027, 0.027, 0.027, 0.027)
),
base_params = list(
base_eir = 25,
smc_efficacy = 0.94,
vaccine_ve_infection = 0.36,
vaccine_ve_severe = 0.32,
itn_efficacy = 0.60,
cost_vaccine_full = 40.20,
cost_smc_per_child_per_year = 6.00,
cost_act = 1.20,
cost_outpatient = 8.00,
cost_hospital = 19.60,
cfr_base = 0.027,
dw_clinical = 0.006,
dw_severe = 0.133,
life_expectancy = 65,
discount_rate = 0.03
),
geo_pop_weights = c(Urban = 0.25, Rural = 0.75)
)
# Display Uganda parameters
kable(uganda_params$age_params, caption = "Uganda: Age-Specific Parameters",
digits = 3, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Age_Group | Population_Weight | Clinical_Risk | Severe_Risk | CFR_Base |
|---|---|---|---|---|
| 0-6m | 0.14 | 0.40 | 0.080 | 0.027 |
| 6-24m | 0.30 | 0.57 | 0.060 | 0.027 |
| 2-5y | 0.56 | 0.53 | 0.020 | 0.027 |
| 5+ | 0.00 | 0.10 | 0.005 | 0.027 |
Burkina Faso Parameters
# Burkina Faso-specific parameters
burkina_params <- list(
country_name = "Burkina Faso",
wtp_threshold = 450, # 0.5 × $900 GDP per capita
geo_params = list(
Urban = list(
itn_coverage = 0.82,
vaccine_coverage = 0.65,
smc_coverage = 0.36,
eir_multiplier = 0.4,
cfr_multiplier = 1.0,
access_delay_days = 0.5,
care_seeking = 0.77,
malaria_positive = 0.35
),
Rural = list(
itn_coverage = 0.80,
vaccine_coverage = 0.55,
smc_coverage = 0.49,
eir_multiplier = 2.0,
cfr_multiplier = 1.6,
access_delay_days = 4.5,
care_seeking = 0.73,
malaria_positive = 0.40
)
),
age_params = data.frame(
Age_Group = c("0-6m", "6-24m", "2-5y", "5+"),
Population_Weight = c(0.13, 0.30, 0.57, 0.00),
Clinical_Risk = c(0.20, 0.35, 0.40, 0.06),
Severe_Risk = c(0.08, 0.06, 0.02, 0.005),
CFR_Base = c(0.013, 0.013, 0.013, 0.013)
),
base_params = list(
base_eir = 25,
smc_efficacy = 0.69,
vaccine_ve_infection = 0.39,
vaccine_ve_severe = 0.29,
itn_efficacy = 0.23,
cost_vaccine_full = 26.08,
cost_smc_per_child_per_year = 6.00,
cost_act = 1.50,
cost_outpatient = 8.00,
cost_hospital = 25.00,
cfr_base = 0.013,
dw_clinical = 0.006,
dw_severe = 0.133,
life_expectancy = 65,
discount_rate = 0.03
),
geo_pop_weights = c(Urban = 0.28, Rural = 0.72)
)
# Display Burkina Faso parameters
kable(burkina_params$age_params, caption = "Burkina Faso: Age-Specific Parameters",
digits = 3, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Age_Group | Population_Weight | Clinical_Risk | Severe_Risk | CFR_Base |
|---|---|---|---|---|
| 0-6m | 0.13 | 0.20 | 0.080 | 0.013 |
| 6-24m | 0.30 | 0.35 | 0.060 | 0.013 |
| 2-5y | 0.57 | 0.40 | 0.020 | 0.013 |
| 5+ | 0.00 | 0.06 | 0.005 | 0.013 |
Part 2: Markov Model Implementation
# Main model function
run_markov_model <- function(params, intervention = "baseline",
time_horizon_years = 10, cohort_size = 10000) {
# Initialize accumulators
total_clinical <- 0
total_severe <- 0
total_deaths <- 0
total_dalys <- 0
total_treatment_cost <- 0
total_intervention_cost <- 0
results <- data.frame()
# Loop over geographic settings
for(geo in c("Urban", "Rural")) {
geo_pop <- cohort_size * params$geo_pop_weights[geo]
gp <- params$geo_params[[geo]]
bp <- params$base_params
# Calculate effective transmission
effective_eir <- bp$base_eir * gp$eir_multiplier
monthly_bites <- effective_eir / 12
monthly_infection_prob <- 1 - exp(-monthly_bites * 0.5)
# Apply ITN effect
itn_effect <- 1 - (bp$itn_efficacy * gp$itn_coverage)
monthly_infection_prob <- monthly_infection_prob * itn_effect
# Apply intervention-specific effects
if(intervention == "smc") {
smc_effect <- 1 - (bp$smc_efficacy * gp$smc_coverage)
monthly_infection_prob <- monthly_infection_prob * smc_effect
intervention_cost <- geo_pop * bp$cost_smc_per_child_per_year * time_horizon_years
} else if(intervention == "vaccine") {
vaccine_effect <- 1 - (bp$vaccine_ve_infection * gp$vaccine_coverage)
monthly_infection_prob <- monthly_infection_prob * vaccine_effect
intervention_cost <- geo_pop * bp$cost_vaccine_full * gp$vaccine_coverage
} else if(intervention == "combined") {
smc_effect <- 1 - (bp$smc_efficacy * gp$smc_coverage)
vaccine_effect <- 1 - (bp$vaccine_ve_infection * gp$vaccine_coverage)
monthly_infection_prob <- monthly_infection_prob * smc_effect * vaccine_effect
intervention_cost <- (geo_pop * bp$cost_smc_per_child_per_year * time_horizon_years) +
(geo_pop * bp$cost_vaccine_full * gp$vaccine_coverage)
} else {
intervention_cost <- 0
}
# Loop over age groups
for(age_idx in 1:nrow(params$age_params)) {
age_pop <- geo_pop * params$age_params$Population_Weight[age_idx]
if(age_pop == 0) next
clinical_risk <- params$age_params$Clinical_Risk[age_idx]
severe_risk <- params$age_params$Severe_Risk[age_idx]
cfr <- params$age_params$CFR_Base[age_idx] * gp$cfr_multiplier
susceptible <- age_pop
clinical_cases <- 0
severe_cases <- 0
deaths <- 0
discounted_treatment_cost <- 0
discounted_dalys <- 0
# Monthly simulation
for(month in 1:(time_horizon_years * 12)) {
discount_factor <- 1 / (1 + bp$discount_rate)^(month/12)
new_infections <- susceptible * monthly_infection_prob * clinical_risk
new_clinical <- new_infections
new_severe <- new_clinical * severe_risk
new_deaths <- new_severe * cfr
clinical_cases <- clinical_cases + new_clinical
severe_cases <- severe_cases + new_severe
deaths <- deaths + new_deaths
# Treatment costs (discounted)
treatment_cost <- (new_clinical * (bp$cost_outpatient + bp$cost_act) +
new_severe * (bp$cost_hospital + bp$cost_act)) * discount_factor
discounted_treatment_cost <- discounted_treatment_cost + treatment_cost
# DALYs (discounted)
yll <- new_deaths * bp$life_expectancy * discount_factor
yld <- (new_clinical * bp$dw_clinical * (14/365) +
new_severe * bp$dw_severe * (21/365)) * discount_factor
discounted_dalys <- discounted_dalys + yll + yld
# Update susceptible population
susceptible <- susceptible - new_infections
if(susceptible < age_pop * 0.1) susceptible <- age_pop * 0.1
}
# Accumulate totals
total_clinical <- total_clinical + clinical_cases
total_severe <- total_severe + severe_cases
total_deaths <- total_deaths + deaths
total_dalys <- total_dalys + discounted_dalys
total_treatment_cost <- total_treatment_cost + discounted_treatment_cost
# Store results
results <- rbind(results, data.frame(
geography = geo,
age_group = params$age_params$Age_Group[age_idx],
clinical_cases = clinical_cases,
severe_cases = severe_cases,
deaths = deaths,
dalys_discounted = discounted_dalys,
treatment_cost_discounted = discounted_treatment_cost,
stringsAsFactors = FALSE
))
}
total_intervention_cost <- total_intervention_cost + intervention_cost
}
total_costs <- total_treatment_cost + total_intervention_cost
return(list(
results = results,
total_clinical = total_clinical,
total_severe = total_severe,
total_deaths = total_deaths,
total_dalys = total_dalys,
total_costs = total_costs,
intervention_cost = total_intervention_cost,
treatment_cost = total_treatment_cost
))
}Model Results
Part 3: Run Model for All Countries
strategies <- c("baseline", "smc", "vaccine", "combined")
strategy_names <- c("Baseline (ITNs)", "SMC + ITNs", "Vaccine + ITNs", "Combined")
countries <- list(Kenya = kenya_params, Uganda = uganda_params, `Burkina Faso` = burkina_params)
# Store all results
all_results <- list()
comparison_table <- data.frame()
for(country_name in names(countries)) {
cat("\n")
cat("============================================================\n")
cat("RUNNING MODEL FOR:", country_name, "\n")
cat("============================================================\n\n")
params <- countries[[country_name]]
country_results <- list()
for(i in 1:length(strategies)) {
cat(sprintf(" Running %s...\n", strategy_names[i]))
result <- run_markov_model(
params = params,
intervention = strategies[i],
time_horizon_years = 10,
cohort_size = 10000
)
country_results[[i]] <- result
cat(sprintf(" Clinical cases: %s\n", format_number(result$total_clinical)))
cat(sprintf(" Severe cases: %s\n", format_number(result$total_severe)))
cat(sprintf(" Deaths: %s\n", format_number(result$total_deaths)))
cat(sprintf(" DALYs (discounted): %s\n", format_number(result$total_dalys)))
cat(sprintf(" Total costs (discounted): $%s\n", format_number(result$total_costs)))
cat("\n")
# Add to comparison table
comparison_table <- rbind(comparison_table, data.frame(
Country = country_name,
WTP_Threshold = params$wtp_threshold,
Strategy = strategy_names[i],
Clinical_Cases = result$total_clinical,
Severe_Cases = result$total_severe,
Deaths = result$total_deaths,
DALYs = result$total_dalys,
Costs = result$total_costs,
Intervention_Cost = result$intervention_cost,
Treatment_Cost = result$treatment_cost,
stringsAsFactors = FALSE
))
}
all_results[[country_name]] <- country_results
}##
## ============================================================
## RUNNING MODEL FOR: Kenya
## ============================================================
##
## Running Baseline (ITNs)...
## Clinical cases: 24,751
## Severe cases: 874
## Deaths: 48
## DALYs (discounted): 2,819
## Total costs (discounted): $242,418
##
## Running SMC + ITNs...
## Clinical cases: 15,648
## Severe cases: 570
## Deaths: 31
## DALYs (discounted): 1,848
## Total costs (discounted): $756,162
##
## Running Vaccine + ITNs...
## Clinical cases: 20,322
## Severe cases: 726
## Deaths: 40
## DALYs (discounted): 2,355
## Total costs (discounted): $451,456
##
## Running Combined...
## Clinical cases: 13,446
## Severe cases: 497
## Deaths: 27
## DALYs (discounted): 1,610
## Total costs (discounted): $985,729
##
##
## ============================================================
## RUNNING MODEL FOR: Uganda
## ============================================================
##
## Running Baseline (ITNs)...
## Clinical cases: 33,639
## Severe cases: 1,334
## Deaths: 55
## DALYs (discounted): 3,177
## Total costs (discounted): $300,459
##
## Running SMC + ITNs...
## Clinical cases: 10,255
## Severe cases: 410
## Deaths: 16
## DALYs (discounted): 986
## Total costs (discounted): $694,382
##
## Running Vaccine + ITNs...
## Clinical cases: 26,068
## Severe cases: 1,035
## Deaths: 42
## DALYs (discounted): 2,475
## Total costs (discounted): $548,097
##
## Running Combined...
## Clinical cases: 9,062
## Severe cases: 363
## Deaths: 15
## DALYs (discounted): 870
## Total costs (discounted): $996,568
##
##
## ============================================================
## RUNNING MODEL FOR: Burkina Faso
## ============================================================
##
## Running Baseline (ITNs)...
## Clinical cases: 32,451
## Severe cases: 1,195
## Deaths: 23
## DALYs (discounted): 1,368
## Total costs (discounted): $303,178
##
## Running SMC + ITNs...
## Clinical cases: 24,042
## Severe cases: 892
## Deaths: 17
## DALYs (discounted): 1,021
## Total costs (discounted): $826,775
##
## Running Vaccine + ITNs...
## Clinical cases: 26,793
## Severe cases: 991
## Deaths: 19
## DALYs (discounted): 1,139
## Total costs (discounted): $402,499
##
## Running Combined...
## Clinical cases: 20,223
## Severe cases: 754
## Deaths: 15
## DALYs (discounted): 867
## Total costs (discounted): $942,673Part 4: Cost-Effectiveness Results (0.5 × GDP Threshold)
# Calculate incremental metrics
ce_table <- comparison_table %>%
group_by(Country) %>%
mutate(
DALYs_Averted = first(DALYs) - DALYs,
Inc_Cost = Costs - first(Costs),
ICER = ifelse(DALYs_Averted > 0, Inc_Cost / DALYs_Averted, NA),
Is_Cost_Effective = ifelse(!is.na(ICER) & ICER < first(WTP_Threshold), "✓ Yes",
ifelse(!is.na(ICER), "✗ No", "Baseline")),
Is_Very_CE = ifelse(!is.na(ICER) & ICER < first(WTP_Threshold)*2, "✓ Yes",
ifelse(!is.na(ICER), "✗ No", "Baseline"))
) %>%
ungroup()
# Format for display
display_table <- ce_table %>%
mutate(
Clinical_Cases = format_number(Clinical_Cases),
Severe_Cases = format_number(Severe_Cases),
Deaths = format_number(Deaths),
DALYs = format_number(DALYs),
Costs = paste0("$", format_number(Costs)),
DALYs_Averted = format_number(DALYs_Averted),
ICER = ifelse(is.na(ICER), "Baseline", paste0("$", format_number(ICER)))
) %>%
select(Country, WTP_Threshold, Strategy, Clinical_Cases, Severe_Cases,
Deaths, DALYs, Costs, ICER, Is_Cost_Effective)
# Display results by country
for(country in unique(display_table$Country)) {
cat("\n")
cat("============================================================\n")
cat("COST-EFFECTIVENESS RESULTS:", country, "\n")
cat("============================================================\n\n")
wtp_value <- unique(display_table$WTP_Threshold[display_table$Country == country])
cat(sprintf("**WTP Threshold (0.5 × GDP):** $%s per DALY averted\n\n", format_number(wtp_value)))
country_table <- display_table %>% filter(Country == country) %>% select(-Country, -WTP_Threshold)
kable(country_table, caption = paste(country, "- Cost-Effectiveness Results"),
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
print()
}##
## ============================================================
## COST-EFFECTIVENESS RESULTS: Kenya
## ============================================================
##
## **WTP Threshold (0.5 × GDP):** $1,150 per DALY averted
##
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Kenya - Cost-Effectiveness Results</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> Clinical_Cases </th>
## <th style="text-align:left;"> Severe_Cases </th>
## <th style="text-align:left;"> Deaths </th>
## <th style="text-align:left;"> DALYs </th>
## <th style="text-align:left;"> Costs </th>
## <th style="text-align:left;"> ICER </th>
## <th style="text-align:left;"> Is_Cost_Effective </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Baseline (ITNs) </td>
## <td style="text-align:left;"> 24,751 </td>
## <td style="text-align:left;"> 874 </td>
## <td style="text-align:left;"> 48 </td>
## <td style="text-align:left;"> 2,819 </td>
## <td style="text-align:left;"> $242,418 </td>
## <td style="text-align:left;"> Baseline </td>
## <td style="text-align:left;"> Baseline </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 15,648 </td>
## <td style="text-align:left;"> 570 </td>
## <td style="text-align:left;"> 31 </td>
## <td style="text-align:left;"> 1,848 </td>
## <td style="text-align:left;"> $756,162 </td>
## <td style="text-align:left;"> $ 529 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 20,322 </td>
## <td style="text-align:left;"> 726 </td>
## <td style="text-align:left;"> 40 </td>
## <td style="text-align:left;"> 2,355 </td>
## <td style="text-align:left;"> $451,456 </td>
## <td style="text-align:left;"> $ 450 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 13,446 </td>
## <td style="text-align:left;"> 497 </td>
## <td style="text-align:left;"> 27 </td>
## <td style="text-align:left;"> 1,610 </td>
## <td style="text-align:left;"> $985,729 </td>
## <td style="text-align:left;"> $ 615 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## </tbody>
## </table>
## ============================================================
## COST-EFFECTIVENESS RESULTS: Uganda
## ============================================================
##
## **WTP Threshold (0.5 × GDP):** $550 per DALY averted
##
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Uganda - Cost-Effectiveness Results</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> Clinical_Cases </th>
## <th style="text-align:left;"> Severe_Cases </th>
## <th style="text-align:left;"> Deaths </th>
## <th style="text-align:left;"> DALYs </th>
## <th style="text-align:left;"> Costs </th>
## <th style="text-align:left;"> ICER </th>
## <th style="text-align:left;"> Is_Cost_Effective </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Baseline (ITNs) </td>
## <td style="text-align:left;"> 33,639 </td>
## <td style="text-align:left;"> 1,334 </td>
## <td style="text-align:left;"> 55 </td>
## <td style="text-align:left;"> 3,177 </td>
## <td style="text-align:left;"> $300,459 </td>
## <td style="text-align:left;"> Baseline </td>
## <td style="text-align:left;"> Baseline </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 10,255 </td>
## <td style="text-align:left;"> 410 </td>
## <td style="text-align:left;"> 16 </td>
## <td style="text-align:left;"> 986 </td>
## <td style="text-align:left;"> $694,382 </td>
## <td style="text-align:left;"> $ 180 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 26,068 </td>
## <td style="text-align:left;"> 1,035 </td>
## <td style="text-align:left;"> 42 </td>
## <td style="text-align:left;"> 2,475 </td>
## <td style="text-align:left;"> $548,097 </td>
## <td style="text-align:left;"> $ 353 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 9,062 </td>
## <td style="text-align:left;"> 363 </td>
## <td style="text-align:left;"> 15 </td>
## <td style="text-align:left;"> 870 </td>
## <td style="text-align:left;"> $996,568 </td>
## <td style="text-align:left;"> $ 302 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## </tbody>
## </table>
## ============================================================
## COST-EFFECTIVENESS RESULTS: Burkina Faso
## ============================================================
##
## **WTP Threshold (0.5 × GDP):** $450 per DALY averted
##
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Burkina Faso - Cost-Effectiveness Results</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> Clinical_Cases </th>
## <th style="text-align:left;"> Severe_Cases </th>
## <th style="text-align:left;"> Deaths </th>
## <th style="text-align:left;"> DALYs </th>
## <th style="text-align:left;"> Costs </th>
## <th style="text-align:left;"> ICER </th>
## <th style="text-align:left;"> Is_Cost_Effective </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Baseline (ITNs) </td>
## <td style="text-align:left;"> 32,451 </td>
## <td style="text-align:left;"> 1,195 </td>
## <td style="text-align:left;"> 23 </td>
## <td style="text-align:left;"> 1,368 </td>
## <td style="text-align:left;"> $303,178 </td>
## <td style="text-align:left;"> Baseline </td>
## <td style="text-align:left;"> Baseline </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 24,042 </td>
## <td style="text-align:left;"> 892 </td>
## <td style="text-align:left;"> 17 </td>
## <td style="text-align:left;"> 1,021 </td>
## <td style="text-align:left;"> $826,775 </td>
## <td style="text-align:left;"> $1,511 </td>
## <td style="text-align:left;"> ✗ No </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 26,793 </td>
## <td style="text-align:left;"> 991 </td>
## <td style="text-align:left;"> 19 </td>
## <td style="text-align:left;"> 1,139 </td>
## <td style="text-align:left;"> $402,499 </td>
## <td style="text-align:left;"> $ 435 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 20,223 </td>
## <td style="text-align:left;"> 754 </td>
## <td style="text-align:left;"> 15 </td>
## <td style="text-align:left;"> 867 </td>
## <td style="text-align:left;"> $942,673 </td>
## <td style="text-align:left;"> $1,278 </td>
## <td style="text-align:left;"> ✗ No </td>
## </tr>
## </tbody>
## </table># Summary interpretation
cat("\n\n## Interpretation by Country (0.5 × GDP per capita threshold):\n\n")##
##
## ## Interpretation by Country (0.5 × GDP per capita threshold):for(country in unique(ce_table$Country)) {
cat(sprintf("\n**%s (WTP = $%s per DALY):**\n",
country, format_number(unique(ce_table$WTP_Threshold[ce_table$Country == country]))))
country_data <- ce_table %>% filter(Country == country)
for(i in 2:nrow(country_data)) {
if(!is.na(country_data$ICER[i])) {
status <- ifelse(country_data$ICER[i] < country_data$WTP_Threshold[i],
"✓ COST-EFFECTIVE", "✗ NOT cost-effective")
cat(sprintf("- %s: %s (ICER = $%s)\n",
country_data$Strategy[i], status, format_number(country_data$ICER[i])))
}
}
}##
## **Kenya (WTP = $1,150 per DALY):**
## - SMC + ITNs: ✓ COST-EFFECTIVE (ICER = $529)
## - Vaccine + ITNs: ✓ COST-EFFECTIVE (ICER = $450)
## - Combined: ✓ COST-EFFECTIVE (ICER = $615)
##
## **Uganda (WTP = $550 per DALY):**
## - SMC + ITNs: ✓ COST-EFFECTIVE (ICER = $180)
## - Vaccine + ITNs: ✓ COST-EFFECTIVE (ICER = $353)
## - Combined: ✓ COST-EFFECTIVE (ICER = $302)
##
## **Burkina Faso (WTP = $450 per DALY):**
## - SMC + ITNs: ✗ NOT cost-effective (ICER = $1,511)
## - Vaccine + ITNs: ✓ COST-EFFECTIVE (ICER = $435)
## - Combined: ✗ NOT cost-effective (ICER = $1,278)Part 5: Geographic Subgroup Analysis
# Aggregate results by geography
geo_summary <- data.frame()
for(country_name in names(all_results)) {
params <- countries[[country_name]]
country_results <- all_results[[country_name]]
for(i in 1:length(strategies)) {
geo_temp <- country_results[[i]]$results %>%
group_by(geography) %>%
summarize(
Clinical_Cases = sum(clinical_cases),
Severe_Cases = sum(severe_cases),
Deaths = sum(deaths),
DALYs_Discounted = sum(dalys_discounted),
Treatment_Costs_Discounted = sum(treatment_cost_discounted),
.groups = 'drop'
) %>%
mutate(Country = country_name,
Strategy = strategy_names[i],
WTP_Threshold = params$wtp_threshold)
geo_summary <- rbind(geo_summary, geo_temp)
}
}
# Calculate geographic ICER matrix
geo_icer <- data.frame()
for(country_name in names(all_results)) {
params <- countries[[country_name]]
country_results <- all_results[[country_name]]
for(geo in c("Urban", "Rural")) {
# Get baseline results for this geography
baseline_results <- country_results[[1]]$results %>%
filter(geography == geo)
base_dalys <- sum(baseline_results$dalys_discounted)
base_cost <- sum(baseline_results$treatment_cost_discounted)
for(i in 2:length(strategies)) {
int_results <- country_results[[i]]$results %>%
filter(geography == geo)
int_dalys <- sum(int_results$dalys_discounted)
int_cost <- sum(int_results$treatment_cost_discounted)
# Get intervention cost for this geography
geo_weight <- params$geo_pop_weights[geo]
int_intervention_cost <- country_results[[i]]$intervention_cost * geo_weight
dalys_averted <- base_dalys - int_dalys
inc_cost <- (int_cost + int_intervention_cost) - base_cost
if(dalys_averted > 0) {
icer <- inc_cost / dalys_averted
is_ce <- icer < params$wtp_threshold
} else {
icer <- NA
is_ce <- FALSE
}
geo_icer <- rbind(geo_icer, data.frame(
Country = country_name,
Geography = geo,
WTP_Threshold = params$wtp_threshold,
Strategy = strategy_names[i],
DALYs_Averted = dalys_averted,
Inc_Cost = inc_cost,
ICER = icer,
Is_Cost_Effective = is_ce
))
}
}
}
# Display geographic ICER matrix
cat("\n\n## Geographic ICER Matrix by Country\n\n")##
##
## ## Geographic ICER Matrix by Countryfor(country in unique(geo_icer$Country)) {
cat(sprintf("\n### %s (WTP = $%s)\n",
country, format_number(unique(geo_icer$WTP_Threshold[geo_icer$Country == country]))))
geo_subset <- geo_icer %>%
filter(Country == country, !is.na(ICER)) %>%
mutate(
ICER_Display = paste0("$", format_number(ICER)),
DALYs_Averted = format_number(DALYs_Averted),
Status = ifelse(Is_Cost_Effective, "✓ CE", "✗ Not CE")
) %>%
select(Geography, Strategy, DALYs_Averted, ICER_Display, Status)
kable(geo_subset, caption = paste(country, "- Geographic ICER Results"),
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
print()
}##
## ### Kenya (WTP = $1,150)
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Kenya - Geographic ICER Results</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> </th>
## <th style="text-align:left;"> Geography </th>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> DALYs_Averted </th>
## <th style="text-align:left;"> ICER_Display </th>
## <th style="text-align:left;"> Status </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 100 </td>
## <td style="text-align:left;"> $1,673 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban1 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 59 </td>
## <td style="text-align:left;"> $1,153 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban2 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 133 </td>
## <td style="text-align:left;"> $1,794 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban3 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 872 </td>
## <td style="text-align:left;"> $ 398 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban4 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 406 </td>
## <td style="text-align:left;"> $ 348 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban5 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 1,077 </td>
## <td style="text-align:left;"> $ 469 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## </tbody>
## </table>
## ### Uganda (WTP = $550)
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Uganda - Geographic ICER Results</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> </th>
## <th style="text-align:left;"> Geography </th>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> DALYs_Averted </th>
## <th style="text-align:left;"> ICER_Display </th>
## <th style="text-align:left;"> Status </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Urban6 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 173 </td>
## <td style="text-align:left;"> $ 725 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban7 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 52 </td>
## <td style="text-align:left;"> $1,371 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban8 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 192 </td>
## <td style="text-align:left;"> $1,046 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban9 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 2,017 </td>
## <td style="text-align:left;"> $ 133 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban10 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 650 </td>
## <td style="text-align:left;"> $ 272 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban11 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 2,114 </td>
## <td style="text-align:left;"> $ 234 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## </tbody>
## </table>
## ### Burkina Faso (WTP = $450)
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Burkina Faso - Geographic ICER Results</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> </th>
## <th style="text-align:left;"> Geography </th>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> DALYs_Averted </th>
## <th style="text-align:left;"> ICER_Display </th>
## <th style="text-align:left;"> Status </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Urban12 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 23 </td>
## <td style="text-align:left;"> $7,031 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban13 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 23 </td>
## <td style="text-align:left;"> $1,477 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban14 </td>
## <td style="text-align:left;"> Urban </td>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 41 </td>
## <td style="text-align:left;"> $4,854 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban15 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> 324 </td>
## <td style="text-align:left;"> $1,122 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban16 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> 205 </td>
## <td style="text-align:left;"> $ 316 </td>
## <td style="text-align:left;"> ✓ CE </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Urban17 </td>
## <td style="text-align:left;"> Rural </td>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> 460 </td>
## <td style="text-align:left;"> $ 963 </td>
## <td style="text-align:left;"> ✗ Not CE </td>
## </tr>
## </tbody>
## </table>Part 6: Probabilistic Sensitivity Analysis (PSA)
n_simulations <- 1000
# Function to run PSA for a country
run_psa <- function(params, n_simulations = n_simulations) {
# Initialize results
psa_results <- data.frame(
sim_id = 1:n_simulations,
icer_smc = NA,
icer_vaccine = NA,
icer_combined = NA
)
# Progress bar
pb <- txtProgressBar(min = 0, max = n_simulations, style = 3)
for(sim in 1:n_simulations) {
# Sample parameters from uncertainty distributions
eir_sample <- rlnorm(1, meanlog = log(params$base_params$base_eir), sdlog = 0.25)
ve_infection_sample <- rbeta(1, shape1 = 35, shape2 = 82)
ve_severe_sample <- rbeta(1, shape1 = 75, shape2 = 32)
itn_efficacy_sample <- rbeta(1, shape1 = 25, shape2 = 75)
smc_efficacy_sample <- rbeta(1, shape1 = 70, shape2 = 30)
# Create temporary params with sampled values
temp_params <- params
temp_params$base_params$base_eir <- eir_sample
temp_params$base_params$vaccine_ve_infection <- ve_infection_sample
temp_params$base_params$vaccine_ve_severe <- ve_severe_sample
temp_params$base_params$itn_efficacy <- itn_efficacy_sample
temp_params$base_params$smc_efficacy <- smc_efficacy_sample
# Run model for each strategy
baseline_result <- run_markov_model(temp_params, "baseline")
smc_result <- run_markov_model(temp_params, "smc")
vaccine_result <- run_markov_model(temp_params, "vaccine")
combined_result <- run_markov_model(temp_params, "combined")
# Calculate ICERs
dalys_smc <- baseline_result$total_dalys - smc_result$total_dalys
cost_smc <- smc_result$total_costs - baseline_result$total_costs
psa_results$icer_smc[sim] <- ifelse(dalys_smc > 0, cost_smc / dalys_smc, NA)
dalys_vaccine <- baseline_result$total_dalys - vaccine_result$total_dalys
cost_vaccine <- vaccine_result$total_costs - baseline_result$total_costs
psa_results$icer_vaccine[sim] <- ifelse(dalys_vaccine > 0, cost_vaccine / dalys_vaccine, NA)
dalys_combined <- baseline_result$total_dalys - combined_result$total_dalys
cost_combined <- combined_result$total_costs - baseline_result$total_costs
psa_results$icer_combined[sim] <- ifelse(dalys_combined > 0, cost_combined / dalys_combined, NA)
setTxtProgressBar(pb, sim)
}
close(pb)
# Calculate summary statistics
psa_summary <- data.frame(
Strategy = c("SMC", "Vaccine", "Combined"),
WTP_Threshold = params$wtp_threshold,
ICER_Mean = c(mean(psa_results$icer_smc, na.rm = TRUE),
mean(psa_results$icer_vaccine, na.rm = TRUE),
mean(psa_results$icer_combined, na.rm = TRUE)),
ICER_Median = c(median(psa_results$icer_smc, na.rm = TRUE),
median(psa_results$icer_vaccine, na.rm = TRUE),
median(psa_results$icer_combined, na.rm = TRUE)),
ICER_Lower95 = c(quantile(psa_results$icer_smc, 0.025, na.rm = TRUE),
quantile(psa_results$icer_vaccine, 0.025, na.rm = TRUE),
quantile(psa_results$icer_combined, 0.025, na.rm = TRUE)),
ICER_Upper95 = c(quantile(psa_results$icer_smc, 0.975, na.rm = TRUE),
quantile(psa_results$icer_vaccine, 0.975, na.rm = TRUE),
quantile(psa_results$icer_combined, 0.975, na.rm = TRUE)),
Prob_CE = c(mean(psa_results$icer_smc <= params$wtp_threshold, na.rm = TRUE),
mean(psa_results$icer_vaccine <= params$wtp_threshold, na.rm = TRUE),
mean(psa_results$icer_combined <= params$wtp_threshold, na.rm = TRUE))
)
return(list(results = psa_results, summary = psa_summary))
}
# Run PSA for each country
cat("\n\n## Probabilistic Sensitivity Analysis Results\n\n")##
##
## ## Probabilistic Sensitivity Analysis Results## Using country-specific WTP thresholds (0.5 × GDP per capita)psa_summaries <- list()
for(country_name in names(countries)) {
cat(sprintf("\n### %s (WTP = $%s)\n",
country_name, format_number(countries[[country_name]]$wtp_threshold)))
cat("\nRunning PSA with", n_simulations,"simulations...\n")
psa_result <- run_psa(countries[[country_name]], n_simulations = n_simulations)
psa_summaries[[country_name]] <- psa_result
# Display summary
psa_summary_display <- psa_result$summary %>%
mutate(
ICER_Mean = paste0("$", format_number(ICER_Mean)),
ICER_Median = paste0("$", format_number(ICER_Median)),
ICER_Lower95 = paste0("$", format_number(ICER_Lower95)),
ICER_Upper95 = paste0("$", format_number(ICER_Upper95)),
Prob_CE = paste0(round(Prob_CE * 100, 1), "%")
) %>%
select(Strategy, ICER_Mean, ICER_Median, ICER_Lower95, ICER_Upper95, Prob_CE)
kable(psa_summary_display, caption = paste(country_name, "- PSA Summary"),
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
print()
}##
## ### Kenya (WTP = $1,150)
##
## Running PSA with 1000 simulations...
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Kenya - PSA Summary</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> ICER_Mean </th>
## <th style="text-align:left;"> ICER_Median </th>
## <th style="text-align:left;"> ICER_Lower95 </th>
## <th style="text-align:left;"> ICER_Upper95 </th>
## <th style="text-align:left;"> Prob_CE </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> SMC </td>
## <td style="text-align:left;"> $549 </td>
## <td style="text-align:left;"> $541 </td>
## <td style="text-align:left;"> $424 </td>
## <td style="text-align:left;"> $714 </td>
## <td style="text-align:left;"> 100% </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine </td>
## <td style="text-align:left;"> $620 </td>
## <td style="text-align:left;"> $602 </td>
## <td style="text-align:left;"> $424 </td>
## <td style="text-align:left;"> $900 </td>
## <td style="text-align:left;"> 100% </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> $665 </td>
## <td style="text-align:left;"> $657 </td>
## <td style="text-align:left;"> $531 </td>
## <td style="text-align:left;"> $836 </td>
## <td style="text-align:left;"> 100% </td>
## </tr>
## </tbody>
## </table>
## ### Uganda (WTP = $550)
##
## Running PSA with 1000 simulations...
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Uganda - PSA Summary</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> ICER_Mean </th>
## <th style="text-align:left;"> ICER_Median </th>
## <th style="text-align:left;"> ICER_Lower95 </th>
## <th style="text-align:left;"> ICER_Upper95 </th>
## <th style="text-align:left;"> Prob_CE </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> SMC </td>
## <td style="text-align:left;"> $170 </td>
## <td style="text-align:left;"> $167 </td>
## <td style="text-align:left;"> $120 </td>
## <td style="text-align:left;"> $244 </td>
## <td style="text-align:left;"> 100% </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine </td>
## <td style="text-align:left;"> $294 </td>
## <td style="text-align:left;"> $285 </td>
## <td style="text-align:left;"> $183 </td>
## <td style="text-align:left;"> $443 </td>
## <td style="text-align:left;"> 99.8% </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> $261 </td>
## <td style="text-align:left;"> $257 </td>
## <td style="text-align:left;"> $201 </td>
## <td style="text-align:left;"> $349 </td>
## <td style="text-align:left;"> 100% </td>
## </tr>
## </tbody>
## </table>
## ### Burkina Faso (WTP = $450)
##
## Running PSA with 1000 simulations...
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Burkina Faso - PSA Summary</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> ICER_Mean </th>
## <th style="text-align:left;"> ICER_Median </th>
## <th style="text-align:left;"> ICER_Lower95 </th>
## <th style="text-align:left;"> ICER_Upper95 </th>
## <th style="text-align:left;"> Prob_CE </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> SMC </td>
## <td style="text-align:left;"> $1,566 </td>
## <td style="text-align:left;"> $1,545 </td>
## <td style="text-align:left;"> $1,214 </td>
## <td style="text-align:left;"> $2,054 </td>
## <td style="text-align:left;"> 0% </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine </td>
## <td style="text-align:left;"> $ 689 </td>
## <td style="text-align:left;"> $ 670 </td>
## <td style="text-align:left;"> $ 430 </td>
## <td style="text-align:left;"> $1,082 </td>
## <td style="text-align:left;"> 4.4% </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> $1,451 </td>
## <td style="text-align:left;"> $1,436 </td>
## <td style="text-align:left;"> $1,132 </td>
## <td style="text-align:left;"> $1,888 </td>
## <td style="text-align:left;"> 0% </td>
## </tr>
## </tbody>
## </table>Part 6.1: Cost effectiveness Plane
# Define GLOBAL number of simulations (same for all countries and interventions)
N_SIMULATIONS <- 1000 # Change this value as needed
cat(sprintf("\n=== PSA CONFIGURATION ===\n"))##
## === PSA CONFIGURATION ===## Number of simulations per strategy: 1,000cat(sprintf("Total simulations per country: %s (3 strategies)\n",
format(N_SIMULATIONS * 3, big.mark = ",")))## Total simulations per country: 3,000 (3 strategies)cat(sprintf("Total simulations across all countries: %s\n",
format(N_SIMULATIONS * 3 * length(countries), big.mark = ",")))## Total simulations across all countries: 9,000## ===========================# Complete PSA function using the global N_SIMULATIONS
run_complete_psa <- function(params, n_simulations = N_SIMULATIONS) {
# Initialize results data frame to store all simulations
psa_results <- data.frame(
sim_id = 1:n_simulations,
# SMC
inc_effect_smc = NA,
inc_cost_smc = NA,
icer_smc = NA,
# Vaccine
inc_effect_vaccine = NA,
inc_cost_vaccine = NA,
icer_vaccine = NA,
# Combined
inc_effect_combined = NA,
inc_cost_combined = NA,
icer_combined = NA
)
# Progress bar
pb <- txtProgressBar(min = 0, max = n_simulations, style = 3)
for(sim in 1:n_simulations) {
# Sample parameters from uncertainty distributions
eir_sample <- rlnorm(1, meanlog = log(params$base_params$base_eir), sdlog = 0.25)
ve_infection_sample <- rbeta(1, shape1 = 35, shape2 = 82)
ve_severe_sample <- rbeta(1, shape1 = 75, shape2 = 32)
itn_efficacy_sample <- rbeta(1, shape1 = 25, shape2 = 75)
smc_efficacy_sample <- rbeta(1, shape1 = 70, shape2 = 30)
cost_vaccine_sample <- rlnorm(1, meanlog = log(params$base_params$cost_vaccine_full), sdlog = 0.1)
# Create temporary params with sampled values
temp_params <- params
temp_params$base_params$base_eir <- eir_sample
temp_params$base_params$vaccine_ve_infection <- ve_infection_sample
temp_params$base_params$vaccine_ve_severe <- ve_severe_sample
temp_params$base_params$itn_efficacy <- itn_efficacy_sample
temp_params$base_params$smc_efficacy <- smc_efficacy_sample
temp_params$base_params$cost_vaccine_full <- cost_vaccine_sample
# Run model for each strategy
baseline_result <- run_markov_model(temp_params, "baseline")
smc_result <- run_markov_model(temp_params, "smc")
vaccine_result <- run_markov_model(temp_params, "vaccine")
combined_result <- run_markov_model(temp_params, "combined")
# Calculate incremental effects and costs for SMC
psa_results$inc_effect_smc[sim] <- baseline_result$total_dalys - smc_result$total_dalys
psa_results$inc_cost_smc[sim] <- smc_result$total_costs - baseline_result$total_costs
if(psa_results$inc_effect_smc[sim] > 0 && !is.na(psa_results$inc_effect_smc[sim])) {
psa_results$icer_smc[sim] <- psa_results$inc_cost_smc[sim] / psa_results$inc_effect_smc[sim]
}
# Calculate incremental effects and costs for Vaccine
psa_results$inc_effect_vaccine[sim] <- baseline_result$total_dalys - vaccine_result$total_dalys
psa_results$inc_cost_vaccine[sim] <- vaccine_result$total_costs - baseline_result$total_costs
if(psa_results$inc_effect_vaccine[sim] > 0 && !is.na(psa_results$inc_effect_vaccine[sim])) {
psa_results$icer_vaccine[sim] <- psa_results$inc_cost_vaccine[sim] / psa_results$inc_effect_vaccine[sim]
}
# Calculate incremental effects and costs for Combined
psa_results$inc_effect_combined[sim] <- baseline_result$total_dalys - combined_result$total_dalys
psa_results$inc_cost_combined[sim] <- combined_result$total_costs - baseline_result$total_costs
if(psa_results$inc_effect_combined[sim] > 0 && !is.na(psa_results$inc_effect_combined[sim])) {
psa_results$icer_combined[sim] <- psa_results$inc_cost_combined[sim] / psa_results$inc_effect_combined[sim]
}
setTxtProgressBar(pb, sim)
}
close(pb)
return(psa_results)
}
# Run PSA for all countries using the global N_SIMULATIONS
cat("\n=== RUNNING PSA SIMULATIONS ===\n")##
## === RUNNING PSA SIMULATIONS ===cat(sprintf("Running %s simulations per strategy for each country...\n", format(N_SIMULATIONS, big.mark = ",")))## Running 1,000 simulations per strategy for each country...## This will take several minutes...psa_results_all <- list()
for(country_name in names(countries)) {
cat(sprintf("\n========================================\n"))
cat(sprintf("Running PSA for %s...\n", country_name))
cat(sprintf("========================================\n"))
psa_results_all[[country_name]] <- run_complete_psa(countries[[country_name]],
n_simulations = N_SIMULATIONS)
# Display summary
n_smc <- sum(!is.na(psa_results_all[[country_name]]$icer_smc))
n_vaccine <- sum(!is.na(psa_results_all[[country_name]]$icer_vaccine))
n_combined <- sum(!is.na(psa_results_all[[country_name]]$icer_combined))
cat(sprintf("\n ✓ SMC: %d of %d valid simulations (%.1f%%)\n",
n_smc, N_SIMULATIONS, n_smc/N_SIMULATIONS*100))
cat(sprintf(" ✓ Vaccine: %d of %d valid simulations (%.1f%%)\n",
n_vaccine, N_SIMULATIONS, n_vaccine/N_SIMULATIONS*100))
cat(sprintf(" ✓ Combined: %d of %d valid simulations (%.1f%%)\n",
n_combined, N_SIMULATIONS, n_combined/N_SIMULATIONS*100))
}##
## ========================================
## Running PSA for Kenya...
## ========================================
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
##
## ✓ SMC: 1000 of 1000 valid simulations (100.0%)
## ✓ Vaccine: 1000 of 1000 valid simulations (100.0%)
## ✓ Combined: 1000 of 1000 valid simulations (100.0%)
##
## ========================================
## Running PSA for Uganda...
## ========================================
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
##
## ✓ SMC: 1000 of 1000 valid simulations (100.0%)
## ✓ Vaccine: 1000 of 1000 valid simulations (100.0%)
## ✓ Combined: 1000 of 1000 valid simulations (100.0%)
##
## ========================================
## Running PSA for Burkina Faso...
## ========================================
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
##
## ✓ SMC: 1000 of 1000 valid simulations (100.0%)
## ✓ Vaccine: 1000 of 1000 valid simulations (100.0%)
## ✓ Combined: 1000 of 1000 valid simulations (100.0%)# Prepare data for plotting - include ALL valid PSA runs
psa_plot_data <- data.frame()
for(country_name in names(psa_results_all)) {
params <- countries[[country_name]]
psa_results <- psa_results_all[[country_name]]
# Add SMC points (all valid simulations)
valid_smc <- !is.na(psa_results$inc_effect_smc) &
!is.na(psa_results$inc_cost_smc) &
psa_results$inc_effect_smc > 0 &
is.finite(psa_results$inc_effect_smc) &
is.finite(psa_results$inc_cost_smc)
if(sum(valid_smc) > 0) {
psa_plot_data <- rbind(psa_plot_data, data.frame(
Country = country_name,
Strategy = "SMC + ITNs",
WTP_Threshold = params$wtp_threshold,
Incremental_DALYs = psa_results$inc_effect_smc[valid_smc],
Incremental_Costs = psa_results$inc_cost_smc[valid_smc]
))
}
# Add Vaccine points (all valid simulations)
valid_vaccine <- !is.na(psa_results$inc_effect_vaccine) &
!is.na(psa_results$inc_cost_vaccine) &
psa_results$inc_effect_vaccine > 0 &
is.finite(psa_results$inc_effect_vaccine) &
is.finite(psa_results$inc_cost_vaccine)
if(sum(valid_vaccine) > 0) {
psa_plot_data <- rbind(psa_plot_data, data.frame(
Country = country_name,
Strategy = "Vaccine + ITNs",
WTP_Threshold = params$wtp_threshold,
Incremental_DALYs = psa_results$inc_effect_vaccine[valid_vaccine],
Incremental_Costs = psa_results$inc_cost_vaccine[valid_vaccine]
))
}
# Add Combined points (all valid simulations)
valid_combined <- !is.na(psa_results$inc_effect_combined) &
!is.na(psa_results$inc_cost_combined) &
psa_results$inc_effect_combined > 0 &
is.finite(psa_results$inc_effect_combined) &
is.finite(psa_results$inc_cost_combined)
if(sum(valid_combined) > 0) {
psa_plot_data <- rbind(psa_plot_data, data.frame(
Country = country_name,
Strategy = "Combined",
WTP_Threshold = params$wtp_threshold,
Incremental_DALYs = psa_results$inc_effect_combined[valid_combined],
Incremental_Costs = psa_results$inc_cost_combined[valid_combined]
))
}
}
cat(sprintf("\n✓ Total points for plotting: %d\n", nrow(psa_plot_data)))##
## ✓ Total points for plotting: 9000cat(sprintf(" (Expected: %d simulations × 3 strategies × 3 countries = %d points)\n",
N_SIMULATIONS, N_SIMULATIONS * 3 * 3))## (Expected: 1000 simulations × 3 strategies × 3 countries = 9000 points)# Calculate mean points for each strategy-country combination
mean_points <- psa_plot_data %>%
group_by(Country, Strategy) %>%
summarise(
Mean_DALYs = mean(Incremental_DALYs, na.rm = TRUE),
Mean_Costs = mean(Incremental_Costs, na.rm = TRUE),
.groups = 'drop'
)
# Create a data frame for the floating label (showing global N_SIMULATIONS)
floating_labels <- data.frame(
Country = names(countries),
x_pos = Inf,
y_pos = Inf,
label = paste0("n = ", format(N_SIMULATIONS, big.mark = ","), " simulations per strategy")
)
# Create color palette for strategies
strategy_colors <- c(
"SMC + ITNs" = "#377EB8",
"Vaccine + ITNs" = "#4DAF4A",
"Combined" = "#E41A1C"
)
# Create the scatter plot with floating label showing global N_SIMULATIONS
p_psa_global <- ggplot() +
# Add all points with low opacity
geom_point(data = psa_plot_data,
aes(x = Incremental_DALYs, y = Incremental_Costs,
color = Strategy),
alpha = 0.15, size = 1.2) +
# Add 95% confidence ellipses
stat_ellipse(data = psa_plot_data,
aes(x = Incremental_DALYs, y = Incremental_Costs,
color = Strategy, fill = Strategy),
type = "norm", level = 0.95, alpha = 0.08,
geom = "polygon", show.legend = FALSE) +
# Add mean points (larger, more visible)
geom_point(data = mean_points,
aes(x = Mean_DALYs, y = Mean_Costs, color = Strategy),
size = 7, shape = 18, stroke = 2) +
# Add WTP threshold line for each country
geom_abline(data = unique(psa_plot_data[, c("Country", "WTP_Threshold")]),
aes(intercept = 0, slope = WTP_Threshold),
linetype = "dashed", color = "red", size = 1.2, alpha = 0.8) +
# Add quadrant reference lines
geom_hline(yintercept = 0, linetype = "solid", color = "gray50", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "solid", color = "gray50", alpha = 0.5) +
# Add floating label showing N_SIMULATIONS (top-right corner of each facet)
geom_text(data = floating_labels,
aes(x = x_pos, y = y_pos, label = label),
hjust = 1.1, vjust = 1.5, size = 4.5,
color = "gray20", fontface = "bold", alpha = 0.8) +
# Facet by country
facet_wrap(~Country, scales = "free", ncol = 1) +
# Labels and theme
labs(
title = "Probabilistic Sensitivity Analysis: Cost-Effectiveness Plane",
subtitle = paste0(
"Each point = one PSA simulation | Large diamonds = mean values | Ellipses = 95% CI\n",
"Kenya WTP: $1,150 | Uganda WTP: $550 | Burkina Faso WTP: $450 per DALY averted"
),
x = "Incremental DALYs Averted (higher is better)",
y = "Incremental Cost (USD, lower is better)"
) +
theme_bw() +
theme(
legend.position = "bottom",
legend.box = "vertical",
legend.title = element_text(size = 11, face = "bold"),
legend.text = element_text(size = 10),
plot.title = element_text(size = 16, face = "bold"),
plot.subtitle = element_text(size = 11, color = "gray40"),
strip.text = element_text(size = 13, face = "bold"),
strip.background = element_rect(fill = "gray95", color = "gray70"),
panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
panel.grid.major = element_line(linetype = "solid", color = "gray85"),
panel.border = element_rect(color = "gray70", linewidth = 0.8)
) +
scale_color_manual(name = "Strategy", values = strategy_colors) +
scale_x_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE)) +
scale_y_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE))
# Display the plot
print(p_psa_global)
# Save high-resolution version
ggsave("Figure_PSA_Global_N.png", p_psa_global, width = 14, height = 16, dpi = 300)# Load required packages
library(ggplot2)
library(dplyr)
library(patchwork)
library(gridExtra)
# Define global number of simulations
N_SIMULATIONS <- 1000
# Prepare PSA plot data (if not already done)
if(!exists("psa_plot_data") || nrow(psa_plot_data) == 0) {
psa_plot_data <- data.frame()
for(country_name in names(psa_results_all)) {
params <- countries[[country_name]]
psa_results <- psa_results_all[[country_name]]
# Add SMC points
valid_smc <- !is.na(psa_results$inc_effect_smc) &
!is.na(psa_results$inc_cost_smc) &
psa_results$inc_effect_smc > 0
if(sum(valid_smc) > 0) {
psa_plot_data <- rbind(psa_plot_data, data.frame(
Country = country_name,
Strategy = "SMC + ITNs",
WTP_Threshold = params$wtp_threshold,
Incremental_DALYs = psa_results$inc_effect_smc[valid_smc],
Incremental_Costs = psa_results$inc_cost_smc[valid_smc]
))
}
# Add Vaccine points
valid_vaccine <- !is.na(psa_results$inc_effect_vaccine) &
!is.na(psa_results$inc_cost_vaccine) &
psa_results$inc_effect_vaccine > 0
if(sum(valid_vaccine) > 0) {
psa_plot_data <- rbind(psa_plot_data, data.frame(
Country = country_name,
Strategy = "Vaccine + ITNs",
WTP_Threshold = params$wtp_threshold,
Incremental_DALYs = psa_results$inc_effect_vaccine[valid_vaccine],
Incremental_Costs = psa_results$inc_cost_vaccine[valid_vaccine]
))
}
# Add Combined points
valid_combined <- !is.na(psa_results$inc_effect_combined) &
!is.na(psa_results$inc_cost_combined) &
psa_results$inc_effect_combined > 0
if(sum(valid_combined) > 0) {
psa_plot_data <- rbind(psa_plot_data, data.frame(
Country = country_name,
Strategy = "Combined",
WTP_Threshold = params$wtp_threshold,
Incremental_DALYs = psa_results$inc_effect_combined[valid_combined],
Incremental_Costs = psa_results$inc_cost_combined[valid_combined]
))
}
}
}
# Calculate mean points
mean_points <- psa_plot_data %>%
group_by(Country, Strategy) %>%
summarise(
Mean_DALYs = mean(Incremental_DALYs, na.rm = TRUE),
Mean_Costs = mean(Incremental_Costs, na.rm = TRUE),
.groups = 'drop'
)
# Color palette
strategy_colors <- c(
"SMC + ITNs" = "#377EB8",
"Vaccine + ITNs" = "#4DAF4A",
"Combined" = "#E41A1C"
)
# Function to create country-specific plot
create_country_plot <- function(country_name, psa_plot_data, mean_points, strategy_colors, N_SIMULATIONS) {
country_data <- psa_plot_data %>% filter(Country == country_name)
country_mean <- mean_points %>% filter(Country == country_name)
wtp_value <- unique(country_data$WTP_Threshold)[1]
# Create the plot
p <- ggplot() +
# Add all points
geom_point(data = country_data,
aes(x = Incremental_DALYs, y = Incremental_Costs,
color = Strategy),
alpha = 0.2, size = 1.5) +
# Add 95% confidence ellipses
stat_ellipse(data = country_data,
aes(x = Incremental_DALYs, y = Incremental_Costs,
color = Strategy, fill = Strategy),
type = "norm", level = 0.95, alpha = 0.1,
geom = "polygon", show.legend = FALSE) +
# Add mean points
geom_point(data = country_mean,
aes(x = Mean_DALYs, y = Mean_Costs, color = Strategy),
size = 6, shape = 18, stroke = 2) +
# Add WTP threshold line
geom_abline(intercept = 0, slope = wtp_value,
linetype = "dashed", color = "red", size = 1.2, alpha = 0.8) +
# Add quadrant lines
geom_hline(yintercept = 0, linetype = "solid", color = "gray50", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "solid", color = "gray50", alpha = 0.5) +
# Labels
labs(
title = country_name,
subtitle = paste0("WTP = $", format(wtp_value, big.mark = ","), " per DALY"),
x = "Incremental DALYs Averted",
y = "Incremental Cost (USD)"
) +
theme_bw() +
theme(
legend.position = "none",
plot.title = element_text(size = 12, face = "bold"),
plot.subtitle = element_text(size = 9, color = "gray40"),
panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
panel.grid.major = element_line(linetype = "solid", color = "gray85"),
panel.border = element_rect(color = "gray70", linewidth = 0.5)
) +
scale_color_manual(values = strategy_colors) +
scale_x_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE)) +
scale_y_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE))
# Add annotation for number of simulations
n_sims_per_strategy <- nrow(country_data) / 3
p <- p + annotate("text", x = -Inf, y = Inf,
label = paste0("n = ", format(N_SIMULATIONS, big.mark = ","), " sims/strategy"),
hjust = -0.1, vjust = 1.5, size = 3.5,
color = "gray30", fontface = "italic")
return(p)
}
# Create individual country plots
p_kenya <- create_country_plot("Kenya", psa_plot_data, mean_points, strategy_colors, N_SIMULATIONS)
p_uganda <- create_country_plot("Uganda", psa_plot_data, mean_points, strategy_colors, N_SIMULATIONS)
p_burkina <- create_country_plot("Burkina Faso", psa_plot_data, mean_points, strategy_colors, N_SIMULATIONS)
# Create combined plot (all countries in one)
p_combined <- ggplot() +
# Add all points with different shapes for countries
geom_point(data = psa_plot_data,
aes(x = Incremental_DALYs, y = Incremental_Costs,
color = Strategy, shape = Country),
alpha = 0.15, size = 1.5) +
# Add 95% confidence ellipses for each strategy-country combination
stat_ellipse(data = psa_plot_data,
aes(x = Incremental_DALYs, y = Incremental_Costs,
color = Strategy, linetype = Country),
type = "norm", level = 0.95, alpha = 0,
size = 0.8, show.legend = TRUE) +
# Add mean points
geom_point(data = mean_points,
aes(x = Mean_DALYs, y = Mean_Costs,
color = Strategy, shape = Country),
size = 5, stroke = 1.5) +
# Add WTP threshold lines for each country
geom_abline(data = unique(psa_plot_data[, c("Country", "WTP_Threshold")]),
aes(intercept = 0, slope = WTP_Threshold, color = Country),
linetype = "dashed", size = 0.8, alpha = 0.5, show.legend = FALSE) +
# Add quadrant lines
geom_hline(yintercept = 0, linetype = "solid", color = "gray50", alpha = 0.5) +
geom_vline(xintercept = 0, linetype = "solid", color = "gray50", alpha = 0.5) +
# Labels
labs(
title = "COMBINED: PSA Cost-Effectiveness Plane - All Countries",
subtitle = paste0(
"Dashed lines = WTP thresholds (0.5× GDP) | Ellipses = 95% CI | Points = PSA simulations\n",
"Kenya WTP: $1,150 | Uganda WTP: $550 | Burkina Faso WTP: $450 per DALY"
),
x = "Incremental DALYs Averted (higher is better)",
y = "Incremental Cost (USD, lower is better)",
color = "Strategy",
shape = "Country",
linetype = "Country"
) +
theme_bw() +
theme(
legend.position = "bottom",
legend.box = "vertical",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 9),
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10, color = "gray40"),
panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
panel.grid.major = element_line(linetype = "solid", color = "gray85"),
panel.border = element_rect(color = "gray70", linewidth = 0.8)
) +
scale_color_manual(name = "Strategy", values = strategy_colors) +
scale_shape_manual(name = "Country",
values = c("Kenya" = 16, "Uganda" = 17, "Burkina Faso" = 18)) +
scale_linetype_manual(name = "Country",
values = c("Kenya" = "solid", "Uganda" = "dashed", "Burkina Faso" = "dotted")) +
scale_x_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE)) +
scale_y_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE)) +
# Add annotation for total simulations
annotate("text", x = -Inf, y = Inf,
label = paste0("Total simulations: ", format(N_SIMULATIONS, big.mark = ",")),
hjust = -0.1, vjust = 1.5, size = 4,
color = "gray30", fontface = "bold")
# Method 1: Using patchwork without the & operator
country_plots <- (p_kenya | p_uganda | p_burkina) +
plot_annotation(
title = "COUNTRY-SPECIFIC PSA Results",
subtitle = paste0("Each panel shows ", format(N_SIMULATIONS, big.mark = ","),
" simulations per strategy (", format(N_SIMULATIONS * 3, big.mark = ","), " total)")
)
# Add a single theme to the patchwork object using the & operator correctly
# Convert the patchwork object to a list and apply theme to each plot
country_plots$patches$plots <- lapply(country_plots$patches$plots, function(x) {
x + theme(legend.position = "none")
})
# Display the country plots
print(country_plots)
# Save the plots
ggsave("Figure_PSA_Country_Specific.png", country_plots, width = 18, height = 8, dpi = 300)
ggsave("Figure_PSA_Combined.png", p_combined, width = 14, height = 10, dpi = 300)Part 7: Cost-Effectiveness Acceptability Curves
# Create CEAC data for all countries with country-specific WTP ranges
ceac_data <- data.frame()
for(country_name in names(psa_summaries)) {
psa_results <- psa_summaries[[country_name]]$results
wtp_threshold <- countries[[country_name]]$wtp_threshold
wtp_range <- seq(0, wtp_threshold * 3, length.out = 100)
for(wtp in wtp_range) {
ceac_data <- rbind(ceac_data, data.frame(
Country = country_name,
WTP = wtp,
WTP_Threshold = wtp_threshold,
SMC = mean(psa_results$icer_smc <= wtp, na.rm = TRUE),
Vaccine = mean(psa_results$icer_vaccine <= wtp, na.rm = TRUE),
Combined = mean(psa_results$icer_combined <= wtp, na.rm = TRUE)
))
}
}
# Melt for plotting
ceac_melted <- ceac_data %>%
pivot_longer(cols = c(SMC, Vaccine, Combined),
names_to = "Strategy", values_to = "Probability")
# Create CEAC plot
p_ceac <- ggplot(ceac_melted, aes(x = WTP, y = Probability, color = Strategy, group = Strategy)) +
geom_line(size = 1.2) +
facet_wrap(~Country, ncol = 1, scales = "free_x") +
geom_vline(aes(xintercept = WTP_Threshold), linetype = "dashed", color = "red", size = 0.8) +
labs(title = "Cost-Effectiveness Acceptability Curves by Country",
subtitle = paste0("Red dashed line = WTP threshold (0.5 × GDP per capita)"),
x = "Willingness-to-pay threshold ($ per DALY averted)",
y = "Probability of cost-effectiveness") +
theme_minimal() +
scale_color_manual(values = c("SMC" = "#377EB8", "Vaccine" = "#4DAF4A", "Combined" = "#E41A1C")) +
scale_y_continuous(labels = percent, limits = c(0, 1)) +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold"))
print(p_ceac)
# Load required packages
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
# Define WTP range from 0 to 1000
wtp_range <- seq(0, 1000, by = 10)
# Create CEAC data for all countries with WTP range 0-1000
ceac_data <- data.frame()
for(country_name in names(psa_results_all)) {
params <- countries[[country_name]]
psa_results <- psa_results_all[[country_name]]
cat(sprintf("Processing CEAC for %s...\n", country_name))
# Calculate probability of cost-effectiveness at each WTP threshold
for(wtp in wtp_range) {
ceac_data <- rbind(ceac_data, data.frame(
Country = country_name,
WTP = wtp,
WTP_Threshold = params$wtp_threshold,
SMC = mean(psa_results$icer_smc <= wtp, na.rm = TRUE),
Vaccine = mean(psa_results$icer_vaccine <= wtp, na.rm = TRUE),
Combined = mean(psa_results$icer_combined <= wtp, na.rm = TRUE)
))
}
}## Processing CEAC for Kenya...
## Processing CEAC for Uganda...
## Processing CEAC for Burkina Faso...# Reshape data for plotting
ceac_long <- ceac_data %>%
pivot_longer(cols = c(SMC, Vaccine, Combined),
names_to = "Strategy",
values_to = "Probability")
# Strategy colors
strategy_colors <- c(
"SMC" = "#377EB8",
"Vaccine" = "#4DAF4A",
"Combined" = "#E41A1C"
)
# Create individual CEAC plots for each country
create_ceac_plot <- function(country_data, country_name, wtp_threshold) {
p <- ggplot(country_data, aes(x = WTP, y = Probability, color = Strategy)) +
geom_line(size = 1.2) +
geom_ribbon(aes(ymin = 0, ymax = Probability, fill = Strategy), alpha = 0.1) +
# Add vertical line at country-specific WTP threshold
geom_vline(xintercept = wtp_threshold,
linetype = "dashed", color = "red", size = 0.8) +
# Add annotation for WTP threshold
annotate("text", x = wtp_threshold, y = 0.95,
label = paste0("WTP = $", format(wtp_threshold, big.mark = ",")),
size = 3.5, angle = 90, hjust = -0.1, color = "red") +
# Labels
labs(
title = country_name,
subtitle = paste0("WTP threshold: $", format(wtp_threshold, big.mark = ","), " per DALY"),
x = "Willingness-to-pay threshold ($ per DALY averted)",
y = "Probability of cost-effectiveness"
) +
theme_minimal() +
theme(
legend.position = "none",
plot.title = element_text(size = 12, face = "bold"),
plot.subtitle = element_text(size = 9, color = "gray40"),
panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
panel.grid.major = element_line(linetype = "solid", color = "gray85"),
panel.border = element_rect(fill = NA, color = "gray70", linewidth = 0.5)
) +
scale_color_manual(values = strategy_colors) +
scale_fill_manual(values = strategy_colors) +
scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, by = 200)) +
scale_y_continuous(labels = scales::percent, limits = c(0, 1))
return(p)
}
# Create individual country CEAC plots
p_ceac_kenya <- create_ceac_plot(
ceac_long %>% filter(Country == "Kenya"),
"Kenya",
unique(ceac_data$WTP_Threshold[ceac_data$Country == "Kenya"])[1]
)
p_ceac_uganda <- create_ceac_plot(
ceac_long %>% filter(Country == "Uganda"),
"Uganda",
unique(ceac_data$WTP_Threshold[ceac_data$Country == "Uganda"])[1]
)
p_ceac_burkina <- create_ceac_plot(
ceac_long %>% filter(Country == "Burkina Faso"),
"Burkina Faso",
unique(ceac_data$WTP_Threshold[ceac_data$Country == "Burkina Faso"])[1]
)
# Create combined CEAC plot (all countries)
p_ceac_combined <- ggplot(ceac_long, aes(x = WTP, y = Probability, color = Strategy, linetype = Country)) +
geom_line(size = 1) +
# Add vertical lines for each country's WTP threshold
geom_vline(data = unique(ceac_data[, c("Country", "WTP_Threshold")]),
aes(xintercept = WTP_Threshold, color = Country),
linetype = "dashed", size = 0.6, alpha = 0.6) +
# Labels
labs(
title = "COMBINED: Cost-Effectiveness Acceptability Curves",
subtitle = "WTP range: $0-$1,000 per DALY averted | Dashed lines = Country-specific WTP thresholds (0.5× GDP)",
x = "Willingness-to-pay threshold ($ per DALY averted)",
y = "Probability of cost-effectiveness",
color = "Strategy",
linetype = "Country"
) +
theme_minimal() +
theme(
legend.position = "bottom",
legend.box = "vertical",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 9),
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10, color = "gray40"),
panel.grid.minor = element_line(linetype = "dotted", color = "gray90"),
panel.grid.major = element_line(linetype = "solid", color = "gray85"),
panel.border = element_rect(fill = NA, color = "gray70", linewidth = 0.5)
) +
scale_color_manual(values = c(strategy_colors,
"Kenya" = "#E69F00",
"Uganda" = "#56B4E9",
"Burkina Faso" = "#009E73")) +
scale_linetype_manual(values = c("Kenya" = "dashed",
"Uganda" = "dotted",
"Burkina Faso" = "dotdash")) +
scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, by = 200)) +
scale_y_continuous(labels = scales::percent, limits = c(0, 1))
# Arrange all CEAC plots in a single figure
ceac_country_plots <- (p_ceac_kenya | p_ceac_uganda | p_ceac_burkina) +
plot_annotation(
title = "Cost-Effectiveness Acceptability Curves (CEAC) by Country",
subtitle = paste0("WTP range: $0 - $1,000 per DALY averted | ",
"Each curve shows probability of cost-effectiveness"),
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, color = "gray40", hjust = 0.5)
)
)
# Add a legend for the country plots (using a separate plot)
legend_plot <- ggplot(ceac_long %>% filter(Country == "Kenya"),
aes(x = WTP, y = Probability, color = Strategy)) +
geom_line(size = 1.2) +
scale_color_manual(values = strategy_colors) +
theme_void() +
theme(legend.position = "bottom",
legend.title = element_text(size = 10, face = "bold"))
# Extract legend
library(cowplot)
ceac_legend <- get_legend(legend_plot)
# Combine country plots with legend
ceac_country_with_legend <- plot_grid(
ceac_country_plots,
ceac_legend,
ncol = 1,
rel_heights = c(1, 0.05)
)
# Display the CEAC plots
print(ceac_country_with_legend)
# Save the plots
ggsave("Figure_CEAC_Country_Specific.png", ceac_country_with_legend, width = 15, height = 6, dpi = 300)
ggsave("Figure_CEAC_Combined.png", p_ceac_combined, width = 12, height = 8, dpi = 300)Part 8: Budget Impact Analysis
# Function to calculate budget impact
calculate_budget_impact <- function(params, bia_horizon_years = 5,
national_u5_population = 5000000) {
# Scale cohort results to national population
scaling_factor <- national_u5_population / 10000
# Get discounted costs from model (10-year)
strategies_list <- c("baseline", "smc", "vaccine", "combined")
total_costs_10yr <- c()
for(i in 1:length(strategies_list)) {
result <- run_markov_model(params, strategies_list[i],
time_horizon_years = 10, cohort_size = 10000)
total_costs_10yr[i] <- result$total_costs
}
# Annualize (simple linear approximation)
annual_costs <- total_costs_10yr / 10 * scaling_factor
# Calculate 5-year costs (discounted)
discount_rate <- params$base_params$discount_rate
five_year_discounted <- sapply(annual_costs, function(cost) {
sum(cost / (1 + discount_rate)^(0:(bia_horizon_years-1)))
})
# Create results table
bia_results <- data.frame(
Strategy = c("Baseline (ITNs)", "SMC + ITNs", "Vaccine + ITNs", "Combined"),
Annual_Cost_M = annual_costs / 1e6,
Five_Year_Cost_M = five_year_discounted / 1e6,
Incremental_5yr_M = (five_year_discounted - five_year_discounted[1]) / 1e6
)
return(bia_results)
}
# Calculate budget impact for each country
cat("\n\n## Budget Impact Analysis (5-year horizon, National Scale)\n\n")##
##
## ## Budget Impact Analysis (5-year horizon, National Scale)## Assumption: 5 million children under 5 nationallybia_all <- data.frame()
for(country_name in names(countries)) {
cat(sprintf("\n### %s\n", country_name))
bia_results <- calculate_budget_impact(countries[[country_name]],
bia_horizon_years = 5,
national_u5_population = 5000000)
bia_results$Country <- country_name
bia_all <- rbind(bia_all, bia_results)
# Display
bia_display <- bia_results %>%
mutate(
Annual_Cost_M = paste0("$", format_number(Annual_Cost_M), "M"),
Five_Year_Cost_M = paste0("$", format_number(Five_Year_Cost_M), "M"),
Incremental_5yr_M = paste0("$", format_number(Incremental_5yr_M), "M")
) %>%
select(Strategy, Annual_Cost_M, Five_Year_Cost_M, Incremental_5yr_M)
kable(bia_display, caption = paste(country_name, "- Budget Impact (5-year)"),
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
print()
# Financial sustainability indicators
cat("\n**Financial Sustainability Indicators:**\n")
cat(sprintf("- Annual health budget (assumed): $500M\n"))
cat(sprintf("- Vaccine strategy share of budget: %.2f%%\n",
bia_results$Annual_Cost_M[3] / 500 * 100))
cat(sprintf("- Combined strategy share of budget: %.2f%%\n",
bia_results$Annual_Cost_M[4] / 500 * 100))
}##
## ### Kenya
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Kenya - Budget Impact (5-year)</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> Annual_Cost_M </th>
## <th style="text-align:left;"> Five_Year_Cost_M </th>
## <th style="text-align:left;"> Incremental_5yr_M </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Baseline (ITNs) </td>
## <td style="text-align:left;"> $12M </td>
## <td style="text-align:left;"> $ 57M </td>
## <td style="text-align:left;"> $ 0M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> $38M </td>
## <td style="text-align:left;"> $178M </td>
## <td style="text-align:left;"> $121M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> $23M </td>
## <td style="text-align:left;"> $106M </td>
## <td style="text-align:left;"> $ 49M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> $49M </td>
## <td style="text-align:left;"> $232M </td>
## <td style="text-align:left;"> $175M </td>
## </tr>
## </tbody>
## </table>
## **Financial Sustainability Indicators:**
## - Annual health budget (assumed): $500M
## - Vaccine strategy share of budget: 4.51%
## - Combined strategy share of budget: 9.86%
##
## ### Uganda
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Uganda - Budget Impact (5-year)</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> Annual_Cost_M </th>
## <th style="text-align:left;"> Five_Year_Cost_M </th>
## <th style="text-align:left;"> Incremental_5yr_M </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Baseline (ITNs) </td>
## <td style="text-align:left;"> $15M </td>
## <td style="text-align:left;"> $ 71M </td>
## <td style="text-align:left;"> $ 0M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> $35M </td>
## <td style="text-align:left;"> $164M </td>
## <td style="text-align:left;"> $ 93M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> $27M </td>
## <td style="text-align:left;"> $129M </td>
## <td style="text-align:left;"> $ 58M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> $50M </td>
## <td style="text-align:left;"> $235M </td>
## <td style="text-align:left;"> $164M </td>
## </tr>
## </tbody>
## </table>
## **Financial Sustainability Indicators:**
## - Annual health budget (assumed): $500M
## - Vaccine strategy share of budget: 5.48%
## - Combined strategy share of budget: 9.97%
##
## ### Burkina Faso
## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Burkina Faso - Budget Impact (5-year)</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Strategy </th>
## <th style="text-align:left;"> Annual_Cost_M </th>
## <th style="text-align:left;"> Five_Year_Cost_M </th>
## <th style="text-align:left;"> Incremental_5yr_M </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Baseline (ITNs) </td>
## <td style="text-align:left;"> $15M </td>
## <td style="text-align:left;"> $ 72M </td>
## <td style="text-align:left;"> $ 0M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> SMC + ITNs </td>
## <td style="text-align:left;"> $41M </td>
## <td style="text-align:left;"> $195M </td>
## <td style="text-align:left;"> $123M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Vaccine + ITNs </td>
## <td style="text-align:left;"> $20M </td>
## <td style="text-align:left;"> $ 95M </td>
## <td style="text-align:left;"> $ 23M </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Combined </td>
## <td style="text-align:left;"> $47M </td>
## <td style="text-align:left;"> $222M </td>
## <td style="text-align:left;"> $151M </td>
## </tr>
## </tbody>
## </table>
## **Financial Sustainability Indicators:**
## - Annual health budget (assumed): $500M
## - Vaccine strategy share of budget: 4.02%
## - Combined strategy share of budget: 9.43%# Budget impact plot
p_budget <- ggplot(bia_all, aes(x = Country, y = Five_Year_Cost_M, fill = Strategy)) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
geom_text(aes(label = paste0("$", format_number(Five_Year_Cost_M), "M")),
position = position_dodge(0.7), vjust = -0.5, size = 3) +
labs(title = "5-Year Budget Requirements by Country",
subtitle = "Discounted at 3% annually | National scale: 5 million children",
x = "", y = "Total Cost (Millions USD)") +
theme_minimal() +
scale_fill_manual(values = c("Baseline (ITNs)" = "gray50",
"Vaccine + ITNs" = "#4DAF4A",
"SMC + ITNs" = "#377EB8",
"Combined" = "#E41A1C")) +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold"))
print(p_budget)
Part 9: Wealth Quintile Equity Analysis
# Define wealth quintile parameters
wealth_params <- data.frame(
Quintile = c("Q1 (Poorest)", "Q2", "Q3", "Q4", "Q5 (Richest)"),
Population_Share = rep(0.20, 5),
itn_coverage = c(0.30, 0.45, 0.55, 0.65, 0.75),
vaccine_coverage = c(0.25, 0.40, 0.55, 0.65, 0.75),
smc_coverage = c(0.35, 0.45, 0.55, 0.65, 0.70),
eir_multiplier = c(2.5, 1.8, 1.2, 0.8, 0.5),
cfr_multiplier = c(2.0, 1.5, 1.2, 0.9, 0.7),
access_delay_days = c(7, 5, 3, 2, 1),
treatment_seeking = c(0.40, 0.55, 0.65, 0.75, 0.85),
malnutrition_rate = c(0.35, 0.25, 0.18, 0.12, 0.08)
)
# Function to run wealth-stratified analysis
run_wealth_analysis <- function(params, intervention = "vaccine",
time_horizon_years = 10, cohort_size = 10000) {
results <- data.frame()
total_dalys <- 0
total_costs <- 0
for(q in 1:nrow(wealth_params)) {
quintile <- wealth_params$Quintile[q]
quintile_pop <- cohort_size * wealth_params$Population_Share[q]
# Create modified params for this quintile
quintile_params <- params
quintile_params$geo_params$Urban$itn_coverage <- wealth_params$itn_coverage[q]
quintile_params$geo_params$Rural$itn_coverage <- wealth_params$itn_coverage[q]
quintile_params$geo_params$Urban$vaccine_coverage <- wealth_params$vaccine_coverage[q]
quintile_params$geo_params$Rural$vaccine_coverage <- wealth_params$vaccine_coverage[q]
quintile_params$geo_params$Urban$smc_coverage <- wealth_params$smc_coverage[q]
quintile_params$geo_params$Rural$smc_coverage <- wealth_params$smc_coverage[q]
# Adjust EIR and CFR
quintile_params$geo_params$Urban$eir_multiplier <- quintile_params$geo_params$Urban$eir_multiplier *
wealth_params$eir_multiplier[q]
quintile_params$geo_params$Rural$eir_multiplier <- quintile_params$geo_params$Rural$eir_multiplier *
wealth_params$eir_multiplier[q]
for(age_idx in 1:nrow(quintile_params$age_params)) {
quintile_params$age_params$CFR_Base[age_idx] <- quintile_params$age_params$CFR_Base[age_idx] *
wealth_params$cfr_multiplier[q]
}
# Run model for this quintile
if(intervention == "baseline") {
result <- run_markov_model(quintile_params, "baseline",
time_horizon_years, quintile_pop)
} else if(intervention == "smc") {
result <- run_markov_model(quintile_params, "smc",
time_horizon_years, quintile_pop)
} else if(intervention == "vaccine") {
result <- run_markov_model(quintile_params, "vaccine",
time_horizon_years, quintile_pop)
} else {
result <- run_markov_model(quintile_params, "combined",
time_horizon_years, quintile_pop)
}
results <- rbind(results, data.frame(
Quintile = quintile,
Quintile_Num = q,
DALYs = result$total_dalys,
Costs = result$total_costs,
Intervention_Cost = result$intervention_cost,
Treatment_Cost = result$treatment_cost
))
total_dalys <- total_dalys + result$total_dalys
total_costs <- total_costs + result$total_costs
}
return(list(results = results, total_dalys = total_dalys, total_costs = total_costs))
}
# Run wealth analysis for Kenya (as example)
cat("\n\n## Wealth Quintile Stratification Analysis (Kenya - Vaccine Strategy)\n\n")##
##
## ## Wealth Quintile Stratification Analysis (Kenya - Vaccine Strategy)cat(sprintf("**WTP Threshold for Kenya:** $%s per DALY\n\n", format_number(kenya_params$wtp_threshold)))## **WTP Threshold for Kenya:** $1,150 per DALY# Baseline wealth results
baseline_wealth <- run_wealth_analysis(kenya_params, "baseline")
vaccine_wealth <- run_wealth_analysis(kenya_params, "vaccine")
# Calculate quintile-specific results
wealth_results <- baseline_wealth$results %>%
rename(Baseline_DALYs = DALYs, Baseline_Costs = Costs) %>%
select(Quintile, Quintile_Num, Baseline_DALYs, Baseline_Costs) %>%
left_join(
vaccine_wealth$results %>%
rename(Vaccine_DALYs = DALYs, Vaccine_Costs = Costs) %>%
select(Quintile, Quintile_Num, Vaccine_DALYs, Vaccine_Costs),
by = c("Quintile", "Quintile_Num")
) %>%
mutate(
DALYs_Averted = Baseline_DALYs - Vaccine_DALYs,
Inc_Cost = (Vaccine_Costs + vaccine_wealth$results$Intervention_Cost) - Baseline_Costs,
ICER = ifelse(DALYs_Averted > 0, Inc_Cost / DALYs_Averted, NA),
Is_CE = ifelse(!is.na(ICER) & ICER < kenya_params$wtp_threshold, "✓ Yes", "✗ No")
)
# Display wealth equity results
wealth_display <- wealth_results %>%
mutate(
DALYs_Averted = format_number(DALYs_Averted),
ICER = ifelse(is.na(ICER), "N/A", paste0("$", format_number(ICER)))
) %>%
select(Quintile, DALYs_Averted, ICER, Is_CE)
kable(wealth_display, caption = "Wealth Quintile Equity Results (Kenya - Vaccine + ITNs)",
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
print()## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Wealth Quintile Equity Results (Kenya - Vaccine + ITNs)</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Quintile </th>
## <th style="text-align:left;"> DALYs_Averted </th>
## <th style="text-align:left;"> ICER </th>
## <th style="text-align:left;"> Is_CE </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Q1 (Poorest) </td>
## <td style="text-align:left;"> 94 </td>
## <td style="text-align:left;"> $ 352 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q2 </td>
## <td style="text-align:left;"> 104 </td>
## <td style="text-align:left;"> $ 516 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q3 </td>
## <td style="text-align:left;"> 102 </td>
## <td style="text-align:left;"> $ 734 </td>
## <td style="text-align:left;"> ✓ Yes </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q4 </td>
## <td style="text-align:left;"> 76 </td>
## <td style="text-align:left;"> $1,176 </td>
## <td style="text-align:left;"> ✗ No </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q5 (Richest) </td>
## <td style="text-align:left;"> 53 </td>
## <td style="text-align:left;"> $1,995 </td>
## <td style="text-align:left;"> ✗ No </td>
## </tr>
## </tbody>
## </table># Calculate concentration index
calculate_concentration_index <- function(health_gains, wealth_rank) {
order_idx <- order(wealth_rank)
health_gains_sorted <- health_gains[order_idx]
pop_shares <- rep(0.2, 5)
cum_pop <- cumsum(pop_shares)
cum_health <- cumsum(health_gains_sorted / sum(health_gains_sorted))
area <- 0
for(i in 1:5) {
area <- area + (cum_pop[i] - ifelse(i==1, 0, cum_pop[i-1])) *
(cum_health[i] + ifelse(i==1, 0, cum_health[i-1])) / 2
}
ci <- 2 * area - 1
return(ci)
}
ci_value <- calculate_concentration_index(wealth_results$DALYs_Averted,
wealth_results$Quintile_Num)
cat("\n**Equity Metrics:**\n")##
## **Equity Metrics:**## - Concentration Index: 0.102if(ci_value > 0.1) {
cat("- Interpretation: Strongly pro-rich distribution\n")
cat("- Recommendation: Increase targeted outreach to poorest quintiles\n")
} else if(ci_value > 0) {
cat("- Interpretation: Mildly pro-rich distribution\n")
cat("- Recommendation: Strengthen equity-focused implementation\n")
} else if(ci_value < -0.1) {
cat("- Interpretation: Strongly pro-poor distribution (desirable)\n")
cat("- Recommendation: Maintain current equity-focused approaches\n")
} else {
cat("- Interpretation: Equitable distribution\n")
cat("- Recommendation: Monitor equity indicators continuously\n")
}## - Interpretation: Strongly pro-rich distribution
## - Recommendation: Increase targeted outreach to poorest quintiles##
## **Coverage Gaps by Wealth Quintile:**coverage_gaps <- wealth_params %>%
mutate(
ITN_Gap = max(itn_coverage) - itn_coverage,
Vaccine_Gap = max(vaccine_coverage) - vaccine_coverage,
SMC_Gap = max(smc_coverage) - smc_coverage
) %>%
select(Quintile, ITN_Gap, Vaccine_Gap, SMC_Gap)
kable(coverage_gaps, caption = "Coverage Gaps Relative to Richest Quintile",
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
print()## <table class="table table-striped table-hover table-condensed" style="color: black; margin-left: auto; margin-right: auto;">
## <caption>Coverage Gaps Relative to Richest Quintile</caption>
## <thead>
## <tr>
## <th style="text-align:left;"> Quintile </th>
## <th style="text-align:right;"> ITN_Gap </th>
## <th style="text-align:right;"> Vaccine_Gap </th>
## <th style="text-align:right;"> SMC_Gap </th>
## </tr>
## </thead>
## <tbody>
## <tr>
## <td style="text-align:left;"> Q1 (Poorest) </td>
## <td style="text-align:right;"> 0.45 </td>
## <td style="text-align:right;"> 0.50 </td>
## <td style="text-align:right;"> 0.35 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q2 </td>
## <td style="text-align:right;"> 0.30 </td>
## <td style="text-align:right;"> 0.35 </td>
## <td style="text-align:right;"> 0.25 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q3 </td>
## <td style="text-align:right;"> 0.20 </td>
## <td style="text-align:right;"> 0.20 </td>
## <td style="text-align:right;"> 0.15 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q4 </td>
## <td style="text-align:right;"> 0.10 </td>
## <td style="text-align:right;"> 0.10 </td>
## <td style="text-align:right;"> 0.05 </td>
## </tr>
## <tr>
## <td style="text-align:left;"> Q5 (Richest) </td>
## <td style="text-align:right;"> 0.00 </td>
## <td style="text-align:right;"> 0.00 </td>
## <td style="text-align:right;"> 0.00 </td>
## </tr>
## </tbody>
## </table>Part 10: Visualizations DALYs Averted by Strategy and Country
p_dalys <- ce_table %>%
filter(Strategy != "Baseline (ITNs)") %>%
mutate(Strategy = factor(Strategy, levels = c("Vaccine + ITNs", "SMC + ITNs", "Combined"))) %>%
ggplot(aes(x = Country, y = DALYs_Averted, fill = Strategy)) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
geom_text(aes(label = format_number(DALYs_Averted)),
position = position_dodge(0.7), vjust = -0.5, size = 3.5) +
labs(title = "DALYs Averted by Country and Strategy",
subtitle = "10-year horizon, 10,000 children | Discounted at 3% per annum",
x = "", y = "DALYs Averted") +
theme_minimal() +
scale_fill_manual(values = c("Vaccine + ITNs" = "#4DAF4A",
"SMC + ITNs" = "#377EB8",
"Combined" = "#E41A1C")) +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold"))
print(p_dalys)
ICER Comparison Across Countries
p_icer <- ce_table %>%
filter(Strategy != "Baseline (ITNs)", !is.na(ICER)) %>%
mutate(Strategy = factor(Strategy, levels = c("Vaccine + ITNs", "SMC + ITNs", "Combined"))) %>%
ggplot(aes(x = Strategy, y = ICER, fill = Country)) +
geom_bar(stat = "identity", position = position_dodge(0.9), width = 0.8) +
geom_hline(yintercept = 800, linetype = "dashed", color = "red", size = 1) +
geom_text(aes(label = paste0("$", format_number(ICER))),
position = position_dodge(0.9), vjust = -0.5, size = 3) +
labs(title = "ICER Comparison Across Countries",
subtitle = paste0("WTP threshold: $800 per DALY averted (red dashed line)"),
x = "", y = "ICER ($ per DALY averted)") +
theme_minimal() +
scale_fill_manual(values = c("Kenya" = "#E69F00", "Uganda" = "#56B4E9",
"Burkina Faso" = "#009E73")) +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold")) +
coord_flip()
print(p_icer)
Clinical Cases by Strategy
p_cases <- comparison_table %>%
ggplot(aes(x = Country, y = Clinical_Cases, fill = Strategy)) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
labs(title = "Clinical Cases by Country and Strategy",
subtitle = "10-year horizon, 10,000 children (undiscounted)",
x = "", y = "Clinical Cases") +
theme_minimal() +
scale_fill_manual(values = c("Baseline (ITNs)" = "gray50",
"Vaccine + ITNs" = "#4DAF4A",
"SMC + ITNs" = "#377EB8",
"Combined" = "#E41A1C")) +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold")) +
scale_y_continuous(labels = comma)
print(p_cases)
Geographic Distribution of DALYs
p_geo <- geo_summary %>%
ggplot(aes(x = geography, y = DALYs_Discounted/1000, fill = Strategy)) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
facet_wrap(~Country, ncol = 1, scales = "free_y") +
labs(title = "Geographic Distribution of DALYs by Country",
subtitle = "Discounted at 3% per annum | Values in thousands",
x = "", y = "Discounted DALYs (thousands)") +
theme_minimal() +
scale_fill_manual(values = c("Baseline (ITNs)" = "gray50",
"Vaccine + ITNs" = "#4DAF4A",
"SMC + ITNs" = "#377EB8",
"Combined" = "#E41A1C")) +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold"))
print(p_geo)
Wealth Quintile Analysis Visualization
# DALYs averted by wealth quintile
p_wealth_dalys <- wealth_results %>%
ggplot(aes(x = Quintile, y = DALYs_Averted, fill = Quintile)) +
geom_bar(stat = "identity", width = 0.7) +
geom_text(aes(label = format_number(DALYs_Averted)), vjust = -0.5, size = 3.5) +
labs(title = "DALYs Averted by Wealth Quintile (Kenya - Vaccine Strategy)",
subtitle = "Poorest (Q1) to Richest (Q5) | 10-year horizon, 10,000 children",
x = "Wealth Quintile", y = "DALYs Averted") +
theme_minimal() +
scale_fill_brewer(palette = "Blues", name = "Quintile") +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold")) +
scale_y_continuous(labels = comma)
print(p_wealth_dalys)
# ICER by wealth quintile
p_wealth_icer <- wealth_results %>%
filter(!is.na(ICER)) %>%
ggplot(aes(x = Quintile, y = ICER, fill = Quintile)) +
geom_bar(stat = "identity", width = 0.7) +
geom_hline(yintercept = 800, linetype = "dashed", color = "red", size = 1) +
geom_text(aes(label = paste0("$", format_number(ICER))), vjust = -0.5, size = 3.5) +
labs(title = "ICER by Wealth Quintile (Kenya - Vaccine Strategy)",
subtitle = paste0("WTP threshold: $800 per DALY (red dashed line)"),
x = "Wealth Quintile", y = "ICER ($ per DALY averted)") +
theme_minimal() +
scale_fill_brewer(palette = "Greens", name = "Quintile") +
theme(legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold")) +
scale_y_continuous(labels = comma, limits = c(0, 250))
print(p_wealth_icer)
Part 11: Complete Results Table
##
##
## # Appendix: Complete Results# Full comparison table
full_table <- ce_table %>%
mutate(
Clinical_Cases = format_number(Clinical_Cases),
Severe_Cases = format_number(Severe_Cases),
Deaths = format_number(Deaths),
DALYs = format_number(DALYs),
Costs = paste0("$", format_number(Costs)),
Intervention_Cost = paste0("$", format_number(Intervention_Cost)),
Treatment_Cost = paste0("$", format_number(Treatment_Cost)),
DALYs_Averted = format_number(DALYs_Averted),
Inc_Cost = ifelse(is.na(Inc_Cost), "Baseline", paste0("$", format_number(Inc_Cost))),
ICER = ifelse(is.na(ICER), "Baseline", paste0("$", format_number(ICER)))
)
kable(full_table, caption = "Complete Results - All Countries and Strategies",
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
font_size = 12) %>%
scroll_box(width = "100%", height = "500px")| Country | WTP_Threshold | Strategy | Clinical_Cases | Severe_Cases | Deaths | DALYs | Costs | Intervention_Cost | Treatment_Cost | DALYs_Averted | Inc_Cost | ICER | Is_Cost_Effective | Is_Very_CE |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Kenya | 1150 | Baseline (ITNs) | 24,751 | 874 | 48 | 2,819 | $242,418 | $ 0 | $242,418 | 0 | $ 0 | Baseline | Baseline | Baseline |
| Kenya | 1150 | SMC + ITNs | 15,648 | 570 | 31 | 1,848 | $756,162 | $600,000 | $156,162 | 972 | $513,744 | $ 529 | ✓ Yes | ✓ Yes |
| Kenya | 1150 | Vaccine + ITNs | 20,322 | 726 | 40 | 2,355 | $451,456 | $250,887 | $200,569 | 465 | $209,039 | $ 450 | ✓ Yes | ✓ Yes |
| Kenya | 1150 | Combined | 13,446 | 497 | 27 | 1,610 | $985,729 | $850,887 | $134,842 | 1,209 | $743,311 | $ 615 | ✓ Yes | ✓ Yes |
| Uganda | 550 | Baseline (ITNs) | 33,639 | 1,334 | 55 | 3,177 | $300,459 | $ 0 | $300,459 | 0 | $ 0 | Baseline | Baseline | Baseline |
| Uganda | 550 | SMC + ITNs | 10,255 | 410 | 16 | 986 | $694,382 | $600,000 | $ 94,382 | 2,190 | $393,923 | $ 180 | ✓ Yes | ✓ Yes |
| Uganda | 550 | Vaccine + ITNs | 26,068 | 1,035 | 42 | 2,475 | $548,097 | $313,560 | $234,537 | 702 | $247,638 | $ 353 | ✓ Yes | ✓ Yes |
| Uganda | 550 | Combined | 9,062 | 363 | 15 | 870 | $996,568 | $913,560 | $ 83,008 | 2,306 | $696,109 | $ 302 | ✓ Yes | ✓ Yes |
| Burkina Faso | 450 | Baseline (ITNs) | 32,451 | 1,195 | 23 | 1,368 | $303,178 | $ 0 | $303,178 | 0 | $ 0 | Baseline | Baseline | Baseline |
| Burkina Faso | 450 | SMC + ITNs | 24,042 | 892 | 17 | 1,021 | $826,775 | $600,000 | $226,775 | 347 | $523,598 | $1,511 | ✗ No | ✗ No |
| Burkina Faso | 450 | Vaccine + ITNs | 26,793 | 991 | 19 | 1,139 | $402,499 | $150,742 | $251,757 | 229 | $ 99,321 | $ 435 | ✓ Yes | ✓ Yes |
| Burkina Faso | 450 | Combined | 20,223 | 754 | 15 | 867 | $942,673 | $750,742 | $191,931 | 500 | $639,496 | $1,278 | ✗ No | ✗ No |