Modelling the Cost-Effectiveness, Budget and Equity Impacts of Malaria Vaccine Introduction alongside other Malaria Prevention Strategies
James Oguta
2026-05-17
Executive Summary
This analysis presents a Markov model for evaluating malaria intervention strategies in children under 5 years across urban and rural settings in sub-Saharan Africa. The model compares four intervention strategies across urban and rural settings in children under 5 years old.
NOTE: The parameters to be updated with country-specific data.
Interventions Being Compared The analysis evaluates four intervention strategies (all layered on top of ITNs as baseline prevention):
Strategy Target Population Efficacy Cost Key Features Baseline (ITNs only) Universal 45% against transmission $2.50 per net Current standard of care SMC + ITNs Children 2-5 years 85% against clinical malaria $6.00 per child/year 4 cycles during rainy season Vaccine + ITNs Children 0-24 months 30% infection, 70% severe $26.00 for 4 doses RTS,S/AS01 vaccine Combined Strategy All ages (layered) Synergistic effect Rural multiplier: 1.5× 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)
# 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 619665 33.1 1417724 75.8 NA 700248 37.4
## Vcells 1187361 9.1 8388608 64.0 49152 1963297 15.0# Load required packages for analysis and visualization
library(ggplot2) # For creating publication-quality plots
library(dplyr) # For data manipulation and aggregation
library(tidyr) # For reshaping data between wide/long formats
library(knitr) # For dynamic report generation
library(kableExtra) # For enhanced table formatting
# Set random seed for reproducibility
set.seed(2026)
# Helper function for formatted numbers with commas
format_number <- function(x) {
result <- ifelse(is.na(x) | is.infinite(x), "NA",
format(round(x, 0), big.mark = ",", scientific = FALSE))
return(result)
}
cat("✓ Environment setup complete\n")## ✓ Environment setup complete## ✓ All packages loaded successfullyModel Structure
# ===================================================================
# MARKOV MODEL DIAGRAM
# ===================================================================
# Load required packages
if(!require(diagram)) {
install.packages("diagram")
library(diagram)
}
# Create a new plot
# Using plot.new() to initialize the plotting area
plot.new()
# Set up plotting parameters
par(mar = c(2, 2, 3, 2), bg = "white")
# Define state positions (normalized coordinates 0-1)
pos <- rbind(
c(0.5, 0.85), # Susceptible (top center)
c(0.5, 0.60), # Clinical (middle center)
c(0.30, 0.35), # Severe (bottom left)
c(0.70, 0.35) # Death (bottom right)
)
# Define colors for states and transitions
state_colors <- c(
"Susceptible" = "#A8E6A3",
"Clinical" = "#FFD3B6",
"Severe" = "#FF8B94",
"Death" = "#AAAAAA"
)
transition_colors <- c(
"progression" = "#E41A1C",
"mortality" = "#000000",
"recovery" = "#377EB8"
)
# State names (with line breaks)
states <- c(
"Susceptible\n(Uninfected)",
"Clinical Malaria\n(Non-severe)",
"Severe Malaria",
"Death"
)
# ==============================================
# DRAW TRANSITIONS
# ==============================================
# 1. Susceptible → Clinical (disease progression)
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"]
)
# Add label for this transition
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"]
)
# 2. Clinical → Severe (progression to severe disease)
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"]
)
# Add label for this transition
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"]
)
# 3. Severe → Death (mortality)
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"]
)
# Add label for this transition
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"]
)
# 4. Clinical → Susceptible (recovery with immunity)
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"
)
# Add label for recovery with immunity
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
)
# 5. Severe → Susceptible (recovery from severe)
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"
)
# Add label for recovery from severe
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 STATE BOXES
# ==============================================
# 1. Susceptible state (top)
textrect(
mid = pos[1, ], radx = 0.15, rady = 0.07,
lab = states[1], cex = 1.0, shadow.size = 0.01,
box.col = state_colors["Susceptible"], col = "white",
lwd = 3, font = 2
)
# 2. Clinical state (middle)
textrect(
mid = pos[2, ], radx = 0.15, rady = 0.07,
lab = states[2], cex = 1.0, shadow.size = 0.01,
box.col = state_colors["Clinical"], col = "white",
lwd = 3, font = 2
)
# 3. Severe state (bottom left)
textrect(
mid = pos[3, ], radx = 0.13, rady = 0.07,
lab = states[3], cex = 1.0, shadow.size = 0.01,
box.col = state_colors["Severe"], col = "white",
lwd = 3, font = 2
)
# 4. Death state (bottom right)
textrect(
mid = pos[4, ], radx = 0.12, rady = 0.07,
lab = states[4], cex = 1.0, shadow.size = 0.01,
box.col = state_colors["Death"], col = "white",
lwd = 3, font = 2
)
# ==============================================
# ADD TITLE AND ANNOTATIONS
# ==============================================
# Main title
title(
main = "Markov Model Structure for Malaria Progression",
cex.main = 1.4, line = 0.8, font.main = 2
)
# Subtitle
text(
x = 0.5, y = 0.96,
labels = "Monthly transitions with time-dependent discounting (3% annually)",
cex = 0.85, font = 3, col = "#555555"
)
# ==============================================
# ADD LEGEND AT TOP RIGHT
# ==============================================
legend(
x = 1.2, y = 0.75,
legend = c("Disease progression", "Mortality", "Recovery"),
col = c(
transition_colors["progression"],
transition_colors["mortality"],
transition_colors["recovery"]
),
lwd = 3, cex = 0.85, bty = "o", bg = "white",
box.col = "black", title = "Transition Types",
horiz = FALSE, xjust = 1, yjust = 1
)
# ==============================================
# ADD DISCOUNTING NOTE AT BOTTOM LEFT
# ==============================================
text(
x = 0.02, y = 0.08,
labels = "NOTE: Only costs and DALYs are discounted (3% annually)",
cex = 0.75, col = "#E41A1C", adj = 0, font = 2
)
text(
x = 0.02, y = 0.05,
labels = "Clinical cases, severe cases, and deaths are NOT discounted",
cex = 0.70, col = "#E41A1C", adj = 0, font = 2
)
# ==============================================
# ADD ARROW ANNOTATIONS
# ==============================================
# Add explanatory text for recovery arrow
text(
x = 0.18, y = 0.72,
labels = "Recovery with\npartial immunity",
cex = 0.70, col = "#377EB8", adj = 0, font = 3
)
# Add explanatory text for severe recovery
text(
x = 0.08, y = 0.48,
labels = "Recovery from\nsevere malaria",
cex = 0.70, col = "#377EB8", adj = 0, font = 3
)
##
## ✓ Markov diagram created successfullyPart 1: Model Parameters
Geographic and Demographic Settings- To be updated
# ===================================================================
# STEP 2: DEFINE GEOGRAPHIC AND DEMOGRAPHIC PARAMETERS
# ===================================================================
# Geographic strata reflecting transmission heterogeneity
geo_settings <- c("Urban", "Rural")
geo_pop_weights <- c(0.30, 0.70) # 30% urban, 70% rural population distribution
# Geographic-specific parameters
geo_params <- list(
Urban = list(
eir_multiplier = 0.4, # Lower transmission intensity
itn_coverage = 0.70, # Higher ITN access
vaccine_coverage = 0.65, # Higher vaccine uptake
smc_coverage = 0.55, # Moderate SMC coverage
cfr_multiplier = 1.0, # Baseline case fatality
access_delay_days = 0.5 # Rapid healthcare access
),
Rural = list(
eir_multiplier = 2.0, # Higher transmission (5x urban)
itn_coverage = 0.45, # Lower ITN coverage
vaccine_coverage = 0.40, # Lower vaccine coverage
smc_coverage = 0.45, # Lower SMC coverage
cfr_multiplier = 1.6, # 60% higher case fatality
access_delay_days = 4.5 # Delayed healthcare access
)
)
# Display parameter table
geo_params_df <- data.frame(
Parameter = c("EIR Multiplier", "ITN Coverage", "Vaccine Coverage",
"SMC Coverage", "CFR Multiplier", "Access Delay (days)"),
Urban = c(geo_params$Urban$eir_multiplier,
geo_params$Urban$itn_coverage,
geo_params$Urban$vaccine_coverage,
geo_params$Urban$smc_coverage,
geo_params$Urban$cfr_multiplier,
geo_params$Urban$access_delay_days),
Rural = c(geo_params$Rural$eir_multiplier,
geo_params$Rural$itn_coverage,
geo_params$Rural$vaccine_coverage,
geo_params$Rural$smc_coverage,
geo_params$Rural$cfr_multiplier,
geo_params$Rural$access_delay_days)
)
kable(geo_params_df, caption = "Geographic Parameters by Setting",
digits = 2, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Parameter | Urban | Rural |
|---|---|---|
| EIR Multiplier | 0.40 | 2.00 |
| ITN Coverage | 0.70 | 0.45 |
| Vaccine Coverage | 0.65 | 0.40 |
| SMC Coverage | 0.55 | 0.45 |
| CFR Multiplier | 1.00 | 1.60 |
| Access Delay (days) | 0.50 | 4.50 |
Age Group Definitions
# Age groups reflecting key epidemiological periods
age_groups <- c("0-6m", "6-24m", "2-5y", "5+")
age_weights <- c(0.05, 0.11, 0.16, 0.68)
age_params_df <- data.frame(
Age_Group = age_groups,
Population_Weight = age_weights,
Clinical_Risk = c(0.25, 0.20, 0.10, 0.04),
Severe_Risk = c(0.08, 0.06, 0.02, 0.005),
CFR_Base = c(0.18, 0.12, 0.06, 0.03),
VE_Infection = c(0.25, 0.35, 0.25, 0.15),
VE_Severe = c(0.60, 0.75, 0.65, 0.50)
)
kable(age_params_df, caption = "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 | VE_Infection | VE_Severe |
|---|---|---|---|---|---|---|
| 0-6m | 0.05 | 0.25 | 0.080 | 0.18 | 0.25 | 0.60 |
| 6-24m | 0.11 | 0.20 | 0.060 | 0.12 | 0.35 | 0.75 |
| 2-5y | 0.16 | 0.10 | 0.020 | 0.06 | 0.25 | 0.65 |
| 5+ | 0.68 | 0.04 | 0.005 | 0.03 | 0.15 | 0.50 |
Epidemiological and Cost Parameters
# Epidemiological parameters
base_eir <- 25 # Baseline EIR: 25 infectious bites/person/year
smc_efficacy <- 0.85 # SMC reduces clinical malaria by 85%
vaccine_ve_infection <- 0.30 # RTS,S efficacy against clinical infection
vaccine_ve_severe <- 0.70 # RTS,S efficacy against severe disease
itn_efficacy <- 0.45 # ITNs reduce transmission by 45%
# Cost parameters (2024 USD)
cost_vaccine_full <- 26.00 # Full 4-dose series per child
cost_smc_per_child_per_year <- 6.00 # 4 cycles at $1.50 per cycle
cost_act <- 1.50 # Artemisinin combination therapy course
cost_outpatient <- 8.00 # Outpatient visit (clinical case)
cost_hospital <- 150.00 # Hospitalization (severe case)
# DALY parameters (GBD 2021)
dw_clinical <- 0.006 # Disability weight for uncomplicated malaria
dw_severe <- 0.133 # Disability weight for severe malaria
duration_clinical_years <- 14/365 # 14 days duration
duration_severe_years <- 21/365 # 21 days duration
life_expectancy <- 65 # Standard for sub-Saharan Africa
# Discounting
discount_rate <- 0.03 # 3% per annum
# Create summary table
params_summary <- data.frame(
Parameter = c("Baseline EIR", "SMC Efficacy", "Vaccine VE (Infection)",
"Vaccine VE (Severe)", "ITN Efficacy", "Vaccine Cost (Full Series)",
"SMC Cost (Annual)", "ACT Treatment Cost", "Outpatient Visit",
"Hospitalization", "DALY Weight (Clinical)", "DALY Weight (Severe)",
"Life Expectancy", "Discount Rate"),
Value = c(base_eir, smc_efficacy, vaccine_ve_infection,
vaccine_ve_severe, itn_efficacy,
paste0("$", cost_vaccine_full),
paste0("$", cost_smc_per_child_per_year),
paste0("$", cost_act), paste0("$", cost_outpatient),
paste0("$", cost_hospital), dw_clinical, dw_severe,
life_expectancy, paste0(discount_rate * 100, "%"))
)
kable(params_summary, caption = "Base Model Parameters",
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Parameter | Value |
|---|---|
| Baseline EIR | 25 |
| SMC Efficacy | 0.85 |
| Vaccine VE (Infection) | 0.3 |
| Vaccine VE (Severe) | 0.7 |
| ITN Efficacy | 0.45 |
| Vaccine Cost (Full Series) | $26 |
| SMC Cost (Annual) | $6 |
| ACT Treatment Cost | $1.5 |
| Outpatient Visit | $8 |
| Hospitalization | $150 |
| DALY Weight (Clinical) | 0.006 |
| DALY Weight (Severe) | 0.133 |
| Life Expectancy | 65 |
| Discount Rate | 3% |
##
## ✓ Parameters loaded successfullycat(sprintf("✓ Discount rate: %.0f%% per annum (applied to COSTS and DALYs ONLY)\n", discount_rate * 100))## ✓ Discount rate: 3% per annum (applied to COSTS and DALYs ONLY)## ✓ Clinical cases, severe cases, and deaths are NOT discountedPart 2: Markov Model Implementation
# ===================================================================
# STEP 8: MAIN MODEL FUNCTION - RUNS MARKOV SIMULATION
# ===================================================================
# NOTE: Only costs and DALYs are discounted at 3% per annum
# Clinical cases, severe cases, and deaths are reported as undiscounted counts
# Main model function
run_stable_model <- function(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 * ifelse(geo == "Urban", 0.30, 0.70)
gp <- geo_params[[geo]]
# Calculate effective transmission
effective_eir <- 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 - (itn_efficacy * gp$itn_coverage)
monthly_infection_prob <- monthly_infection_prob * itn_effect
# Apply strategy-specific effects
intervention_cost <- 0
if(intervention == "smc") {
smc_effect <- 1 - (smc_efficacy * gp$smc_coverage)
monthly_infection_prob <- monthly_infection_prob * smc_effect
intervention_cost <- geo_pop * cost_smc_per_child_per_year * time_horizon_years
} else if(intervention == "vaccine") {
vaccine_effect <- 1 - (vaccine_ve_infection * gp$vaccine_coverage)
monthly_infection_prob <- monthly_infection_prob * vaccine_effect
intervention_cost <- geo_pop * cost_vaccine_full * gp$vaccine_coverage
} else if(intervention == "combined") {
smc_effect <- 1 - (smc_efficacy * gp$smc_coverage)
vaccine_effect <- 1 - (vaccine_ve_infection * gp$vaccine_coverage)
monthly_infection_prob <- monthly_infection_prob * smc_effect * vaccine_effect
smc_cost <- geo_pop * cost_smc_per_child_per_year * time_horizon_years
vaccine_cost <- geo_pop * cost_vaccine_full * gp$vaccine_coverage
delivery_multiplier <- ifelse(geo == "Rural", 1.5, 1)
intervention_cost <- (smc_cost + vaccine_cost) * delivery_multiplier
}
total_intervention_cost <- total_intervention_cost + intervention_cost
# Loop over age groups
for(a in 1:length(age_groups)) {
age <- age_groups[a]
age_pop <- geo_pop * age_weights[a]
# Age-specific parameters
if(age == "0-6m") {
clinical_risk <- 0.25
severe_risk <- 0.08
cfr <- 0.18 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
} else if(age == "6-24m") {
clinical_risk <- 0.20
severe_risk <- 0.06
cfr <- 0.12 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
} else if(age == "2-5y") {
clinical_risk <- 0.10
severe_risk <- 0.02
cfr <- 0.06 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
} else {
clinical_risk <- 0.04
severe_risk <- 0.005
cfr <- 0.03 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
}
# Apply vaccine effect on severe disease
if(intervention %in% c("vaccine", "combined")) {
severe_risk <- severe_risk * (1 - vaccine_ve_severe * gp$vaccine_coverage)
}
# Initialize susceptible population
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 + discount_rate)^(month/12)
new_infections <- susceptible * monthly_infection_prob
new_clinical <- new_infections * clinical_risk
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 * (cost_outpatient + cost_act) +
new_severe * (cost_hospital + cost_act)) * discount_factor
discounted_treatment_cost <- discounted_treatment_cost + treatment_cost
# DALYs (discounted)
yll <- new_deaths * life_expectancy * discount_factor
yld <- (new_clinical * dw_clinical * duration_clinical_years +
new_severe * dw_severe * duration_severe_years) * 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 = age,
clinical_cases = clinical_cases,
severe_cases = severe_cases,
deaths = deaths,
dalys_discounted = discounted_dalys,
treatment_cost_discounted = discounted_treatment_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_discounted = total_intervention_cost,
treatment_cost_discounted = total_treatment_cost
))
}Model Results
Part 3: Run All Strategies
# ===================================================================
# STEP 9: EXECUTE MODEL FOR ALL FOUR INTERVENTION STRATEGIES
# ===================================================================
strategies <- c("baseline", "smc", "vaccine", "combined")
strategy_names <- c("Baseline (ITNs)", "SMC + ITNs", "Vaccine + ITNs", "Combined")
results_list <- list()
for(i in 1:length(strategies)) {
cat(sprintf("Running %s...\n", strategy_names[i]))
result <- run_stable_model(intervention = strategies[i],
time_horizon_years = 10,
cohort_size = 10000)
results_list[[i]] <- result
cat(sprintf(" Clinical cases: %.0f\n", result$total_clinical))
cat(sprintf(" Severe cases: %.0f\n", result$total_severe))
cat(sprintf(" Deaths: %.0f\n", result$total_deaths))
cat(sprintf(" DALYs (discounted): %.0f\n", result$total_dalys))
cat(sprintf(" Total costs (discounted): $%.0f\n", result$total_costs))
cat("\n")
}## Running Baseline (ITNs)...
## Clinical cases: 5791
## Severe cases: 207
## Deaths: 58
## DALYs (discounted): 3292
## Total costs (discounted): $75813
##
## Running SMC + ITNs...
## Clinical cases: 3715
## Severe cases: 133
## Deaths: 37
## DALYs (discounted): 2131
## Total costs (discounted): $648973
##
## Running Vaccine + ITNs...
## Clinical cases: 5110
## Severe cases: 127
## Deaths: 36
## DALYs (discounted): 2069
## Total costs (discounted): $183113
##
## Running Combined...
## Clinical cases: 3306
## Severe cases: 82
## Deaths: 23
## DALYs (discounted): 1349
## Total costs (discounted): $1008756Part 4: Cost-Effectiveness Results (Discounted)
# ===================================================================
# STEP 10: CALCULATE COST-EFFECTIVENESS METRICS
# ===================================================================
# NOTE: Costs and DALYs are discounted at 3% per annum
# Cases and deaths are undiscounted (reported as is)
# Create results table
ce_table <- data.frame(
Strategy = strategy_names,
Clinical_Cases = sapply(results_list, function(x) x$total_clinical),
Severe_Cases = sapply(results_list, function(x) x$total_severe),
Deaths = sapply(results_list, function(x) x$total_deaths),
DALYs = sapply(results_list, function(x) x$total_dalys),
Costs = sapply(results_list, function(x) x$total_costs)
)
# Calculate incremental outcomes
ce_table$DALYs_Averted <- ce_table$DALYs[1] - ce_table$DALYs
ce_table$Inc_Cost <- ce_table$Costs - ce_table$Costs[1]
ce_table$ICER <- ce_table$Inc_Cost / ce_table$DALYs_Averted
ce_table$ICER[ce_table$DALYs_Averted <= 0] <- NA
# 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), "Dominated", paste0("$", format_number(ICER)))
) %>%
select(Strategy, Clinical_Cases, Severe_Cases, Deaths, DALYs, Costs, DALYs_Averted, ICER)
kable(display_table, caption = "Cost-Effectiveness Results (10-year horizon, 10,000 children)",
format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
row_spec(0, bold = TRUE)| Strategy | Clinical_Cases | Severe_Cases | Deaths | DALYs | Costs | DALYs_Averted | ICER |
|---|---|---|---|---|---|---|---|
| Baseline (ITNs) | 5,791 | 207 | 58 | 3,292 | $ 75,813 | 0 | Dominated |
| SMC + ITNs | 3,715 | 133 | 37 | 2,131 | $ 648,973 | 1,161 | $494 |
| Vaccine + ITNs | 5,110 | 127 | 36 | 2,069 | $ 183,113 | 1,224 | $ 88 |
| Combined | 3,306 | 82 | 23 | 1,349 | $1,008,756 | 1,943 | $480 |
# WTP threshold interpretation
cat("\n### Interpretation (WTP threshold: $800 per DALY averted):\n\n")##
## ### Interpretation (WTP threshold: $800 per DALY averted):for(i in 2:nrow(ce_table)) {
if(!is.na(ce_table$ICER[i])) {
if(ce_table$ICER[i] < 800) {
cat(sprintf("- **%s**: ✓ Cost-effective (ICER = $%s)\n",
ce_table$Strategy[i], format_number(ce_table$ICER[i])))
} else {
cat(sprintf("- **%s**: ✗ Not cost-effective (ICER = $%s)\n",
ce_table$Strategy[i], format_number(ce_table$ICER[i])))
}
}
}## - **SMC + ITNs**: ✓ Cost-effective (ICER = $494)
## - **Vaccine + ITNs**: ✓ Cost-effective (ICER = $88)
## - **Combined**: ✓ Cost-effective (ICER = $480)Part 5: Geographic Subgroup Analysis
# ===================================================================
# GEOGRAPHIC SUBGROUP ANALYSIS
# ===================================================================
# Aggregate results by geography
geo_summary <- data.frame()
for(i in 1:length(strategies)) {
geo_temp <- results_list[[i]]$results %>%
group_by(geography) %>%
summarise(
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(Strategy = strategy_names[i])
geo_summary <- rbind(geo_summary, geo_temp)
}
# Display geographic results
for(geo in c("Urban", "Rural")) {
cat(sprintf("\n========================================\n"))
cat(sprintf("GEOGRAPHIC SETTING: %s\n", toupper(geo)))
cat(sprintf("========================================\n\n"))
# Filter for this geography
geo_subset <- geo_summary[geo_summary$geography == geo, ]
# Extract values safely - convert to numeric
base_idx <- which(geo_subset$Strategy == strategy_names[1])
base_dalys <- as.numeric(geo_subset$DALYs_Discounted[base_idx])
base_cost <- as.numeric(geo_subset$Treatment_Costs_Discounted[base_idx])
# Clinical outcomes table
cat("Clinical Outcomes (Discounted):\n")
cat("--------------------------------\n")
for(i in 1:nrow(geo_subset)) {
cat(sprintf("- %s: Clinical cases = %s, Deaths = %s, DALYs = %s\n",
geo_subset$Strategy[i],
format_number(geo_subset$Clinical_Cases[i]),
format_number(geo_subset$Deaths[i]),
format_number(geo_subset$DALYs_Discounted[i])))
}
cat("\nCost-Effectiveness:\n")
cat("-------------------\n")
# Calculate ICERs for each strategy
for(i in 2:length(strategies)) {
int_idx <- which(geo_subset$Strategy == strategy_names[i])
# Convert to numeric to avoid list issues
int_dalys <- as.numeric(geo_subset$DALYs_Discounted[int_idx])
int_cost <- as.numeric(geo_subset$Treatment_Costs_Discounted[int_idx])
# Calculate intervention cost for this geography
intervention_cost_total <- as.numeric(results_list[[i]]$intervention_cost_discounted)
geo_intervention_cost <- intervention_cost_total * (ifelse(geo == "Urban", 0.3, 0.7))
# Calculate incremental values
dalys_averted <- base_dalys - int_dalys
inc_cost <- (int_cost + geo_intervention_cost) - base_cost
# Calculate ICER if positive DALYs averted
if(dalys_averted > 0 && is.finite(dalys_averted) && !is.na(dalys_averted)) {
icer <- inc_cost / dalys_averted
cat(sprintf("- %s: DALYs averted = %s, ICER = $%s\n",
strategy_names[i],
format_number(dalys_averted),
format_number(icer)))
} else {
cat(sprintf("- %s: Dominated or extended dominance\n", strategy_names[i]))
}
}
# Determine best strategy for this geography
cat("\nMost Cost-Effective Strategy:\n")
cat("------------------------------\n")
# Calculate ICERs for all strategies (excluding baseline)
icer_values <- c()
strategy_labels <- c()
for(i in 2:length(strategies)) {
int_idx <- which(geo_subset$Strategy == strategy_names[i])
int_dalys <- as.numeric(geo_subset$DALYs_Discounted[int_idx])
int_cost <- as.numeric(geo_subset$Treatment_Costs_Discounted[int_idx])
intervention_cost_total <- as.numeric(results_list[[i]]$intervention_cost_discounted)
geo_intervention_cost <- intervention_cost_total * (ifelse(geo == "Urban", 0.3, 0.7))
dalys_averted <- base_dalys - int_dalys
inc_cost <- (int_cost + geo_intervention_cost) - base_cost
if(dalys_averted > 0 && is.finite(dalys_averted) && !is.na(dalys_averted)) {
icer <- inc_cost / dalys_averted
icer_values <- c(icer_values, icer)
strategy_labels <- c(strategy_labels, strategy_names[i])
}
}
if(length(icer_values) > 0) {
best_idx <- which.min(icer_values)
cat(sprintf("✓ %s (ICER: $%s per DALY averted)\n",
strategy_labels[best_idx],
format_number(icer_values[best_idx])))
}
cat("\n")
}##
## ========================================
## GEOGRAPHIC SETTING: URBAN
## ========================================
##
## Clinical Outcomes (Discounted):
## --------------------------------
## - Baseline (ITNs): Clinical cases = 816, Deaths = 4, DALYs = 231
## - SMC + ITNs: Clinical cases = 507, Deaths = 2, DALYs = 146
## - Vaccine + ITNs: Clinical cases = 687, Deaths = 2, DALYs = 107
## - Combined: Clinical cases = 439, Deaths = 1, DALYs = 69
##
## Cost-Effectiveness:
## -------------------
## - SMC + ITNs: DALYs averted = 85, ICER = $2,064
## - Vaccine + ITNs: DALYs averted = 124, ICER = $272
## - Combined: DALYs averted = 162, ICER = $1,760
##
## Most Cost-Effective Strategy:
## ------------------------------
## ✓ Vaccine + ITNs (ICER: $272 per DALY averted)
##
##
## ========================================
## GEOGRAPHIC SETTING: RURAL
## ========================================
##
## Clinical Outcomes (Discounted):
## --------------------------------
## - Baseline (ITNs): Clinical cases = 4,975, Deaths = 54, DALYs = 3,061
## - SMC + ITNs: Clinical cases = 3,208, Deaths = 35, DALYs = 1,985
## - Vaccine + ITNs: Clinical cases = 4,423, Deaths = 34, DALYs = 1,962
## - Combined: Clinical cases = 2,867, Deaths = 22, DALYs = 1,280
##
## Cost-Effectiveness:
## -------------------
## - SMC + ITNs: DALYs averted = 1,076, ICER = $369
## - Vaccine + ITNs: DALYs averted = 1,099, ICER = $67
## - Combined: DALYs averted = 1,781, ICER = $364
##
## Most Cost-Effective Strategy:
## ------------------------------
## ✓ Vaccine + ITNs (ICER: $67 per DALY averted)# ===================================================================
# ALTERNATIVE: Create a formatted table using kable
# ===================================================================
# Create summary tables for reporting
cat("\n\n========================================\n")##
##
## ========================================## GEOGRAPHIC SUMMARY TABLES## ========================================# Urban summary table
urban_summary <- geo_summary[geo_summary$geography == "Urban", ] %>%
mutate(
Clinical_Cases = format_number(Clinical_Cases),
Severe_Cases = format_number(Severe_Cases),
Deaths = format_number(Deaths),
DALYs = format_number(DALYs_Discounted),
Costs = paste0("$", format_number(Treatment_Costs_Discounted))
) %>%
select(Strategy, Clinical_Cases, Severe_Cases, Deaths, DALYs, Costs)
cat("URBAN SETTINGS:\n")## URBAN SETTINGS:## # A tibble: 4 × 6
## Strategy Clinical_Cases Severe_Cases Deaths DALYs Costs
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 Baseline (ITNs) 816 "29" 4 "231" $10,830
## 2 SMC + ITNs 507 "18" 2 "146" $ 6,833
## 3 Vaccine + ITNs 687 "13" 2 "107" $ 7,650
## 4 Combined 439 " 9" 1 " 69" $ 4,960# Rural summary table
rural_summary <- geo_summary[geo_summary$geography == "Rural", ] %>%
mutate(
Clinical_Cases = format_number(Clinical_Cases),
Severe_Cases = format_number(Severe_Cases),
Deaths = format_number(Deaths),
DALYs = format_number(DALYs_Discounted),
Costs = paste0("$", format_number(Treatment_Costs_Discounted))
) %>%
select(Strategy, Clinical_Cases, Severe_Cases, Deaths, DALYs, Costs)
cat("\nRURAL SETTINGS:\n")##
## RURAL SETTINGS:## # A tibble: 4 × 6
## Strategy Clinical_Cases Severe_Cases Deaths DALYs Costs
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 Baseline (ITNs) 4,975 "178" 54 3,061 $64,983
## 2 SMC + ITNs 3,208 "115" 35 1,985 $42,140
## 3 Vaccine + ITNs 4,423 "114" 34 1,962 $51,963
## 4 Combined 2,867 " 74" 22 1,280 $33,895# ===================================================================
# CALCULATE GEOGRAPHIC ICER MATRIX
# ===================================================================
cat("\n\n========================================\n")##
##
## ========================================## GEOGRAPHIC ICER MATRIX## ========================================# Create empty data frame for geographic ICERs
geo_icer_matrix <- data.frame(
Strategy = strategy_names[2:4],
Urban_ICER = NA,
Rural_ICER = NA,
Urban_DALYs_Averted = NA,
Rural_DALYs_Averted = NA
)
for(i in 2:length(strategies)) {
# Urban calculations
urban_base_dalys <- as.numeric(geo_summary[geo_summary$geography == "Urban" &
geo_summary$Strategy == strategy_names[1], "DALYs_Discounted"])
urban_int_dalys <- as.numeric(geo_summary[geo_summary$geography == "Urban" &
geo_summary$Strategy == strategy_names[i], "DALYs_Discounted"])
urban_base_cost <- as.numeric(geo_summary[geo_summary$geography == "Urban" &
geo_summary$Strategy == strategy_names[1], "Treatment_Costs_Discounted"])
urban_int_cost <- as.numeric(geo_summary[geo_summary$geography == "Urban" &
geo_summary$Strategy == strategy_names[i], "Treatment_Costs_Discounted"])
urban_intervention_cost <- as.numeric(results_list[[i]]$intervention_cost_discounted) * 0.3
urban_dalys_averted <- urban_base_dalys - urban_int_dalys
urban_inc_cost <- (urban_int_cost + urban_intervention_cost) - urban_base_cost
if(urban_dalys_averted > 0 && is.finite(urban_dalys_averted)) {
geo_icer_matrix$Urban_ICER[i-1] <- urban_inc_cost / urban_dalys_averted
geo_icer_matrix$Urban_DALYs_Averted[i-1] <- urban_dalys_averted
}
# Rural calculations
rural_base_dalys <- as.numeric(geo_summary[geo_summary$geography == "Rural" &
geo_summary$Strategy == strategy_names[1], "DALYs_Discounted"])
rural_int_dalys <- as.numeric(geo_summary[geo_summary$geography == "Rural" &
geo_summary$Strategy == strategy_names[i], "DALYs_Discounted"])
rural_base_cost <- as.numeric(geo_summary[geo_summary$geography == "Rural" &
geo_summary$Strategy == strategy_names[1], "Treatment_Costs_Discounted"])
rural_int_cost <- as.numeric(geo_summary[geo_summary$geography == "Rural" &
geo_summary$Strategy == strategy_names[i], "Treatment_Costs_Discounted"])
rural_intervention_cost <- as.numeric(results_list[[i]]$intervention_cost_discounted) * 0.7
rural_dalys_averted <- rural_base_dalys - rural_int_dalys
rural_inc_cost <- (rural_int_cost + rural_intervention_cost) - rural_base_cost
if(rural_dalys_averted > 0 && is.finite(rural_dalys_averted)) {
geo_icer_matrix$Rural_ICER[i-1] <- rural_inc_cost / rural_dalys_averted
geo_icer_matrix$Rural_DALYs_Averted[i-1] <- rural_dalys_averted
}
}
# Display ICER matrix
geo_icer_matrix_display <- geo_icer_matrix %>%
mutate(
Urban_ICER = ifelse(is.na(Urban_ICER), "Dominated", paste0("$", format_number(Urban_ICER))),
Rural_ICER = ifelse(is.na(Rural_ICER), "Dominated", paste0("$", format_number(Rural_ICER))),
Urban_DALYs_Averted = format_number(Urban_DALYs_Averted),
Rural_DALYs_Averted = format_number(Rural_DALYs_Averted)
)
print(geo_icer_matrix_display)## Strategy Urban_ICER Rural_ICER Urban_DALYs_Averted Rural_DALYs_Averted
## 1 SMC + ITNs $2,064 $369 85 1,076
## 2 Vaccine + ITNs $ 272 $ 67 124 1,099
## 3 Combined $1,760 $364 162 1,781# ===================================================================
# GEOGRAPHIC EQUITY ANALYSIS
# ===================================================================
cat("\n\n========================================\n")##
##
## ========================================## GEOGRAPHIC EQUITY ANALYSIS## ========================================# Calculate DALYs averted per capita by geography
combined_idx <- which(strategy_names == "Combined")
# Urban population: 30% of 10,000 = 3,000
# Rural population: 70% of 10,000 = 7,000
urban_pop <- 3000
rural_pop <- 7000
urban_dalys_averted <- geo_icer_matrix$Urban_DALYs_Averted[3] # Combined strategy
rural_dalys_averted <- geo_icer_matrix$Rural_DALYs_Averted[3]
urban_dalys_per_capita <- urban_dalys_averted / urban_pop * 1000
rural_dalys_per_capita <- rural_dalys_averted / rural_pop * 1000
cat("Combined Strategy Equity Analysis:\n")## Combined Strategy Equity Analysis:## ----------------------------------## Urban population: 3,000## Rural population: 7,000##
## Urban DALYs averted: 162## Rural DALYs averted: 1,781##
## Urban DALYs averted per 1000: 54.0## Rural DALYs averted per 1000: 254.5equity_ratio <- rural_dalys_per_capita / urban_dalys_per_capita
cat(sprintf("\nEquity Ratio (Rural/Urban): %.2f\n", equity_ratio))##
## Equity Ratio (Rural/Urban): 4.71if(equity_ratio < 0.8) {
cat("\n⚠️ STRONG PRO-URBAN BIAS DETECTED\n")
cat("Rural populations benefit significantly less from interventions\n")
cat("\nRecommendations to improve equity:\n")
cat(" 1. Increase rural coverage targets by 30-50%\n")
cat(" 2. Subsidize transport costs for rural families\n")
cat(" 3. Implement mobile vaccination campaigns in remote areas\n")
cat(" 4. Train community health workers for SMC delivery\n")
} else if(equity_ratio < 1.2) {
cat("\n✓ EQUITABLE DISTRIBUTION\n")
cat("Benefits are relatively balanced across geographic settings\n")
} else {
cat("\n✓ PRO-RURAL DISTRIBUTION\n")
cat("Rural populations benefit more, which is desirable for equity\n")
}##
## ✓ PRO-RURAL DISTRIBUTION
## Rural populations benefit more, which is desirable for equity# ===================================================================
# RECOMMENDATIONS BY GEOGRAPHY
# ===================================================================
cat("\n\n========================================\n")##
##
## ========================================## RECOMMENDATIONS BY GEOGRAPHY## ========================================## URBAN AREAS:## ------------urban_best_idx <- which.min(geo_icer_matrix$Urban_ICER[!is.na(geo_icer_matrix$Urban_ICER)])
urban_best_strategy <- geo_icer_matrix$Strategy[urban_best_idx]
urban_best_icer <- min(geo_icer_matrix$Urban_ICER[!is.na(geo_icer_matrix$Urban_ICER)])
cat(sprintf("Most cost-effective: %s (ICER: $%s per DALY)\n",
urban_best_strategy, format_number(urban_best_icer)))## Most cost-effective: Vaccine + ITNs (ICER: $272 per DALY)##
## Recommendations for URBAN areas:## - Focus on Vaccine + ITNs (excellent value for money)## - Maintain high ITN coverage (currently 70%)## - Target vaccination in 6-24 month age group## - Strengthen case management at existing facilities##
## RURAL AREAS:## ------------rural_best_idx <- which.min(geo_icer_matrix$Rural_ICER[!is.na(geo_icer_matrix$Rural_ICER)])
rural_best_strategy <- geo_icer_matrix$Strategy[rural_best_idx]
rural_best_icer <- min(geo_icer_matrix$Rural_ICER[!is.na(geo_icer_matrix$Rural_ICER)])
cat(sprintf("Most cost-effective: %s (ICER: $%s per DALY)\n",
rural_best_strategy, format_number(rural_best_icer)))## Most cost-effective: Vaccine + ITNs (ICER: $67 per DALY)##
## Recommendations for RURAL areas:## - Prioritize Vaccine + ITNs (most cost-effective)## - Increase vaccine coverage from 40% to 65%## - Add SMC during rainy season for 2-5 year olds## - Deploy mobile health teams for outreach## - Train community health workers for vaccine delivery##
##
## ========================================## GEOGRAPHIC ANALYSIS COMPLETE## ========================================PROBABILISTIC SENSITIVITY ANALYSIS (PSA)
PSA Function
run_psa_analysis <- function(n_simulations = 10000,
time_horizon_years = 10,
cohort_size = 10000) {
cat("\n========================================\n")
cat("PROBABILISTIC SENSITIVITY ANALYSIS (PSA)\n")
cat("========================================\n")
cat(sprintf("Number of simulations: %d\n", n_simulations))
cat("This may take several minutes...\n\n")
# Initialize results data frame
psa_results <- data.frame(
sim_id = 1:n_simulations,
# Baseline outcomes
baseline_cases = NA,
baseline_deaths = NA,
baseline_dalys = NA,
baseline_costs = NA,
# SMC outcomes
smc_cases = NA,
smc_deaths = NA,
smc_dalys = NA,
smc_costs = NA,
# Vaccine outcomes
vaccine_cases = NA,
vaccine_deaths = NA,
vaccine_dalys = NA,
vaccine_costs = NA,
# Combined outcomes
combined_cases = NA,
combined_deaths = NA,
combined_dalys = NA,
combined_costs = NA,
# Incremental outcomes
inc_effect_smc = NA,
inc_cost_smc = NA,
inc_effect_vaccine = NA,
inc_cost_vaccine = NA,
inc_effect_combined = NA,
inc_cost_combined = NA,
# ICERs
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(25), sdlog = 0.25)
ve_infection_sample <- rbeta(1, shape1 = 35, shape2 = 82)
ve_severe_sample <- rbeta(1, shape1 = 75, shape2 = 32)
smc_efficacy_sample <- rbeta(1, shape1 = 90, shape2 = 16)
itn_efficacy_sample <- rbeta(1, shape1 = 50, shape2 = 61)
# Coverage samples
itn_cov_urban <- rbeta(1, shape1 = 72, shape2 = 31)
itn_cov_rural <- rbeta(1, shape1 = 47, shape2 = 57)
vaccine_cov_urban <- rbeta(1, shape1 = 67, shape2 = 36)
vaccine_cov_rural <- rbeta(1, shape1 = 42, shape2 = 63)
smc_cov_urban <- rbeta(1, shape1 = 57, shape2 = 46)
smc_cov_rural <- rbeta(1, shape1 = 47, shape2 = 57)
# Clinical parameters
clinical_risk_u5 <- rbeta(1, shape1 = 22, shape2 = 88)
severe_risk_u5 <- rbeta(1, shape1 = 7, shape2 = 110)
cfr_sample <- rbeta(1, shape1 = 14, shape2 = 103)
# Costs
cost_vaccine_sample <- rgamma(1, shape = 12, rate = 0.46)
cost_smc_sample <- rgamma(1, shape = 10, rate = 1.67)
cost_act_sample <- rgamma(1, shape = 12, rate = 8.0)
cost_outpatient_sample <- rgamma(1, shape = 10, rate = 1.25)
cost_hospital_sample <- rgamma(1, shape = 15, rate = 0.10)
# DALY weights
dw_clinical_sample <- rbeta(1, shape1 = 6, shape2 = 994)
dw_severe_sample <- rbeta(1, shape1 = 140, shape2 = 920)
# Geographic weights
urban_pop <- cohort_size * 0.30
rural_pop <- cohort_size * 0.70
# Run for each strategy (simplified for PSA speed)
for(strategy in c("baseline", "smc", "vaccine", "combined")) {
total_cases <- 0
total_deaths <- 0
total_dalys <- 0
total_costs <- 0
for(geo in c("Urban", "Rural")) {
if(geo == "Urban") {
geo_pop <- urban_pop
eir_mult <- 0.4
itn_cov <- itn_cov_urban
vaccine_cov <- vaccine_cov_urban
smc_cov <- smc_cov_urban
cfr_mult <- 1.0
} else {
geo_pop <- rural_pop
eir_mult <- 2.0
itn_cov <- itn_cov_rural
vaccine_cov <- vaccine_cov_rural
smc_cov <- smc_cov_rural
cfr_mult <- 1.6
}
# Calculate transmission
effective_eir <- eir_sample * eir_mult
monthly_bites <- effective_eir / 12
monthly_infection_prob <- 1 - exp(-monthly_bites * 0.5)
# ITN effect
itn_effect <- 1 - (itn_efficacy_sample * itn_cov)
monthly_infection_prob <- monthly_infection_prob * itn_effect
# Strategy-specific effects
if(strategy == "smc") {
smc_effect <- 1 - (smc_efficacy_sample * smc_cov)
monthly_infection_prob <- monthly_infection_prob * smc_effect
intervention_cost <- geo_pop * cost_smc_sample * time_horizon_years
severe_risk_adj <- severe_risk_u5
} else if(strategy == "vaccine") {
vaccine_effect <- 1 - (ve_infection_sample * vaccine_cov)
monthly_infection_prob <- monthly_infection_prob * vaccine_effect
intervention_cost <- geo_pop * cost_vaccine_sample * vaccine_cov
severe_risk_adj <- severe_risk_u5 * (1 - ve_severe_sample * vaccine_cov)
} else if(strategy == "combined") {
smc_effect <- 1 - (smc_efficacy_sample * smc_cov)
vaccine_effect <- 1 - (ve_infection_sample * vaccine_cov)
monthly_infection_prob <- monthly_infection_prob * smc_effect * vaccine_effect
smc_cost <- geo_pop * cost_smc_sample * time_horizon_years
vaccine_cost <- geo_pop * cost_vaccine_sample * vaccine_cov
delivery_multiplier <- ifelse(geo == "Rural", 1.5, 1)
intervention_cost <- (smc_cost + vaccine_cost) * delivery_multiplier
severe_risk_adj <- severe_risk_u5 * (1 - ve_severe_sample * vaccine_cov)
} else {
intervention_cost <- 0
severe_risk_adj <- severe_risk_u5
}
# Adjust CFR
cfr_adj <- cfr_sample * cfr_mult
# Run simulation
susceptible <- geo_pop
geo_cases <- 0
geo_deaths <- 0
geo_treatment_cost <- 0
for(month in 1:(time_horizon_years * 12)) {
discount_factor <- 1 / (1 + 0.03)^(month/12)
new_infections <- susceptible * monthly_infection_prob
new_clinical <- new_infections * clinical_risk_u5
new_severe <- new_clinical * severe_risk_adj
new_deaths <- new_severe * cfr_adj
geo_cases <- geo_cases + new_clinical
geo_deaths <- geo_deaths + new_deaths
treatment_cost <- (new_clinical * (cost_outpatient_sample + cost_act_sample) +
new_severe * (cost_hospital_sample + cost_act_sample)) * discount_factor
geo_treatment_cost <- geo_treatment_cost + treatment_cost
susceptible <- susceptible - new_infections
if(susceptible < geo_pop * 0.1) susceptible <- geo_pop * 0.1
}
# Calculate DALYs
yll <- geo_deaths * life_expectancy
yld <- geo_cases * dw_clinical_sample * 0.038 +
geo_cases * severe_risk_adj * dw_severe_sample * 0.058
geo_dalys <- (yll + yld) * 0.85
total_cases <- total_cases + geo_cases
total_deaths <- total_deaths + geo_deaths
total_dalys <- total_dalys + geo_dalys
total_costs <- total_costs + geo_treatment_cost + intervention_cost
}
# Store results
if(strategy == "baseline") {
psa_results$baseline_cases[sim] <- total_cases
psa_results$baseline_deaths[sim] <- total_deaths
psa_results$baseline_dalys[sim] <- total_dalys
psa_results$baseline_costs[sim] <- total_costs
} else if(strategy == "smc") {
psa_results$smc_cases[sim] <- total_cases
psa_results$smc_deaths[sim] <- total_deaths
psa_results$smc_dalys[sim] <- total_dalys
psa_results$smc_costs[sim] <- total_costs
} else if(strategy == "vaccine") {
psa_results$vaccine_cases[sim] <- total_cases
psa_results$vaccine_deaths[sim] <- total_deaths
psa_results$vaccine_dalys[sim] <- total_dalys
psa_results$vaccine_costs[sim] <- total_costs
} else {
psa_results$combined_cases[sim] <- total_cases
psa_results$combined_deaths[sim] <- total_deaths
psa_results$combined_dalys[sim] <- total_dalys
psa_results$combined_costs[sim] <- total_costs
}
}
# Calculate incremental outcomes
psa_results$inc_effect_smc[sim] <- psa_results$baseline_dalys[sim] - psa_results$smc_dalys[sim]
psa_results$inc_cost_smc[sim] <- psa_results$smc_costs[sim] - psa_results$baseline_costs[sim]
psa_results$inc_effect_vaccine[sim] <- psa_results$baseline_dalys[sim] - psa_results$vaccine_dalys[sim]
psa_results$inc_cost_vaccine[sim] <- psa_results$vaccine_costs[sim] - psa_results$baseline_costs[sim]
psa_results$inc_effect_combined[sim] <- psa_results$baseline_dalys[sim] - psa_results$combined_dalys[sim]
psa_results$inc_cost_combined[sim] <- psa_results$combined_costs[sim] - psa_results$baseline_costs[sim]
# Calculate ICERs
psa_results$icer_smc[sim] <- ifelse(psa_results$inc_effect_smc[sim] > 0,
psa_results$inc_cost_smc[sim] / psa_results$inc_effect_smc[sim], NA)
psa_results$icer_vaccine[sim] <- ifelse(psa_results$inc_effect_vaccine[sim] > 0,
psa_results$inc_cost_vaccine[sim] / psa_results$inc_effect_vaccine[sim], NA)
psa_results$icer_combined[sim] <- ifelse(psa_results$inc_effect_combined[sim] > 0,
psa_results$inc_cost_combined[sim] / psa_results$inc_effect_combined[sim], NA)
setTxtProgressBar(pb, sim)
}
close(pb)
# Calculate summary statistics
psa_summary <- data.frame(
Strategy = c("SMC", "Vaccine", "Combined"),
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_SD = c(sd(psa_results$icer_smc, na.rm = TRUE),
sd(psa_results$icer_vaccine, na.rm = TRUE),
sd(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_100 = c(mean(psa_results$icer_smc <= 100, na.rm = TRUE),
mean(psa_results$icer_vaccine <= 100, na.rm = TRUE),
mean(psa_results$icer_combined <= 100, na.rm = TRUE)),
Prob_CE_500 = c(mean(psa_results$icer_smc <= 500, na.rm = TRUE),
mean(psa_results$icer_vaccine <= 500, na.rm = TRUE),
mean(psa_results$icer_combined <= 500, na.rm = TRUE)),
Prob_CE_800 = c(mean(psa_results$icer_smc <= 800, na.rm = TRUE),
mean(psa_results$icer_vaccine <= 800, na.rm = TRUE),
mean(psa_results$icer_combined <= 800, na.rm = TRUE))
)
return(list(results = psa_results, summary = psa_summary))
}Run PSA
# Run PSA (10000 simulations)
psa_output <- run_psa_analysis(n_simulations = 10000,
time_horizon_years = 10,
cohort_size = 10000)##
## ========================================
## PROBABILISTIC SENSITIVITY ANALYSIS (PSA)
## ========================================
## Number of simulations: 10000
## This may take several minutes...
##
## | | | 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%PSA Summary Statistics
kable(psa_summary, caption = "PSA Summary Statistics (10000 simulations)",
digits = 0, format = "html") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Strategy | ICER_Mean | ICER_Median | ICER_SD | ICER_Lower95 | ICER_Upper95 | Prob_CE_100 | Prob_CE_500 | Prob_CE_800 |
|---|---|---|---|---|---|---|---|---|
| SMC | 217 | 172 | 173 | 43 | 651 | 0 | 1 | 1 |
| Vaccine | 28 | 21 | 31 | -4 | 109 | 1 | 1 | 1 |
| Combined | 212 | 171 | 158 | 50 | 614 | 0 | 1 | 1 |
PSA Visualizations
Cost-Effectiveness Planes
# Function to create CE plane
create_ce_plane <- function(psa_results, strategy_name) {
if(strategy_name == "SMC") {
data <- data.frame(
effect = psa_results$inc_effect_smc,
cost = psa_results$inc_cost_smc,
icer = psa_results$icer_smc
)
color <- "#377EB8"
} else if(strategy_name == "Vaccine") {
data <- data.frame(
effect = psa_results$inc_effect_vaccine,
cost = psa_results$inc_cost_vaccine,
icer = psa_results$icer_vaccine
)
color <- "#4DAF4A"
} else {
data <- data.frame(
effect = psa_results$inc_effect_combined,
cost = psa_results$inc_cost_combined,
icer = psa_results$icer_combined
)
color <- "#E41A1C"
}
data <- data[!is.na(data$effect) & !is.na(data$cost) & data$effect > 0, ]
# Calculate 95% confidence ellipse
if(nrow(data) > 2) {
ellipse_data <- ellipse(cov(data[, c("effect", "cost")]),
centre = colMeans(data[, c("effect", "cost")]),
level = 0.95)
ellipse_df <- as.data.frame(ellipse_data)
colnames(ellipse_df) <- c("effect", "cost")
} else {
ellipse_df <- NULL
}
p <- ggplot(data, aes(x = effect, y = cost)) +
geom_point(alpha = 0.3, size = 1.5, color = color) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
geom_abline(intercept = 0, slope = 100, linetype = "dotted", color = "darkgreen") +
geom_abline(intercept = 0, slope = 500, linetype = "dotted", color = "orange") +
geom_abline(intercept = 0, slope = 800, linetype = "dotted", color = "red") +
labs(title = paste("Cost-Effectiveness Plane:", strategy_name),
x = "Incremental DALYs Averted",
y = "Incremental Cost (USD)") +
theme_minimal()
if(!is.null(ellipse_df)) {
p <- p + geom_path(data = ellipse_df, aes(x = effect, y = cost),
color = color, size = 1.2, linetype = "dashed")
}
p <- p + geom_point(aes(x = mean(effect), y = mean(cost)),
size = 4, color = color, shape = 18)
return(p)
}
# Create individual CE planes
p_smc <- create_ce_plane(psa_results, "SMC")
p_vaccine <- create_ce_plane(psa_results, "Vaccine")
p_combined <- create_ce_plane(psa_results, "Combined")
# Display plots
print(p_smc)
Combined CE Plane with All Strategies
combined_data <- rbind(
data.frame(effect = psa_results$inc_effect_smc,
cost = psa_results$inc_cost_smc,
Strategy = "SMC"),
data.frame(effect = psa_results$inc_effect_vaccine,
cost = psa_results$inc_cost_vaccine,
Strategy = "Vaccine"),
data.frame(effect = psa_results$inc_effect_combined,
cost = psa_results$inc_cost_combined,
Strategy = "Combined")
)
combined_data <- combined_data[!is.na(combined_data$effect) &
!is.na(combined_data$cost) &
combined_data$effect > 0, ]
p_combined_all <- ggplot(combined_data, aes(x = effect, y = cost, color = Strategy, fill = Strategy)) +
geom_point(alpha = 0.2, size = 1) +
stat_ellipse(level = 0.95, geom = "polygon", alpha = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
geom_abline(intercept = 0, slope = 100, linetype = "dotted", color = "darkgreen") +
geom_abline(intercept = 0, slope = 500, linetype = "dotted", color = "orange") +
geom_abline(intercept = 0, slope = 800, linetype = "dotted", color = "red") +
labs(title = "Cost-Effectiveness Plane: All Strategies",
subtitle = "95% confidence ellipses shown | WTP thresholds: $100, $500, $800/DALY",
x = "Incremental DALYs Averted",
y = "Incremental Cost (USD)") +
theme_minimal() +
scale_color_manual(values = c("SMC" = "#377EB8", "Vaccine" = "#4DAF4A", "Combined" = "#E41A1C")) +
scale_fill_manual(values = c("SMC" = "#377EB8", "Vaccine" = "#4DAF4A", "Combined" = "#E41A1C")) +
theme(legend.position = "bottom")
print(p_combined_all)
ICER Distribution
icer_data <- data.frame(
SMC = psa_results$icer_smc,
Vaccine = psa_results$icer_vaccine,
Combined = psa_results$icer_combined
) %>%
pivot_longer(everything(), names_to = "Strategy", values_to = "ICER") %>%
filter(!is.na(ICER), ICER <= 5000)
p_icer_dist <- ggplot(icer_data, aes(x = Strategy, y = ICER, fill = Strategy)) +
geom_violin(alpha = 0.5, trim = TRUE, scale = "width") +
geom_boxplot(width = 0.2, alpha = 0.8, outlier.size = 0.5) +
geom_jitter(width = 0.15, alpha = 0.3, size = 0.8) +
geom_hline(yintercept = 100, linetype = "dotted", color = "darkgreen", size = 1) +
geom_hline(yintercept = 500, linetype = "dotted", color = "orange", size = 1) +
geom_hline(yintercept = 800, linetype = "dashed", color = "red", size = 1) +
labs(title = "ICER Distribution by Strategy",
subtitle = "Violin plots with boxplots | WTP thresholds: $100, $500, $800 per DALY",
x = "",
y = "ICER ($ per DALY averted)") +
theme_minimal() +
scale_fill_manual(values = c("SMC" = "#377EB8", "Vaccine" = "#4DAF4A", "Combined" = "#E41A1C")) +
theme(legend.position = "none") +
scale_y_continuous(labels = comma, limits = c(0, 2000))
print(p_icer_dist)
Cost-Effectiveness Acceptability Curve (CEAC)
# Create CEAC data
wtp_thresholds <- seq(0, 2000, 10)
ceac_data <- data.frame(
WTP = wtp_thresholds,
SMC = sapply(wtp_thresholds, function(wtp) mean(psa_results$icer_smc <= wtp, na.rm = TRUE)),
Vaccine = sapply(wtp_thresholds, function(wtp) mean(psa_results$icer_vaccine <= wtp, na.rm = TRUE)),
Combined = sapply(wtp_thresholds, function(wtp) mean(psa_results$icer_combined <= wtp, na.rm = TRUE))
) %>%
pivot_longer(cols = -WTP, names_to = "Strategy", values_to = "Probability")
p_ceac <- ggplot(ceac_data, aes(x = WTP, y = Probability, color = Strategy, group = Strategy)) +
geom_line(size = 1.2) +
geom_ribbon(aes(ymin = 0, ymax = Probability, fill = Strategy), alpha = 0.1) +
geom_vline(xintercept = 100, linetype = "dotted", color = "darkgreen", size = 0.8) +
geom_vline(xintercept = 500, linetype = "dotted", color = "orange", size = 0.8) +
geom_vline(xintercept = 800, linetype = "dashed", color = "red", size = 0.8) +
annotate("text", x = 100, y = 0.95, label = "$100", size = 3, angle = 90, hjust = -0.2) +
annotate("text", x = 500, y = 0.95, label = "$500", size = 3, angle = 90, hjust = -0.2) +
annotate("text", x = 800, y = 0.95, label = "$800", size = 3, angle = 90, hjust = -0.2) +
labs(title = "Cost-Effectiveness Acceptability Curve (CEAC)",
subtitle = "Probability of being cost-effective at different willingness-to-pay thresholds",
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_fill_manual(values = c("SMC" = "#377EB8", "Vaccine" = "#4DAF4A", "Combined" = "#E41A1C")) +
scale_x_continuous(labels = comma, limits = c(0, 2000)) +
scale_y_continuous(labels = percent, limits = c(0, 1)) +
theme(legend.position = "bottom")
print(p_ceac)
ICER Histograms
icer_hist_data <- data.frame(
ICER = c(psa_results$icer_smc, psa_results$icer_vaccine, psa_results$icer_combined),
Strategy = rep(c("SMC", "Vaccine", "Combined"), each = nrow(psa_results))
) %>%
filter(!is.na(ICER), ICER <= 3000)
p_hist <- ggplot(icer_hist_data, aes(x = ICER, fill = Strategy)) +
geom_histogram(alpha = 0.6, bins = 50, position = "identity") +
geom_vline(xintercept = 100, linetype = "dotted", color = "darkgreen", size = 0.8) +
geom_vline(xintercept = 500, linetype = "dotted", color = "orange", size = 0.8) +
geom_vline(xintercept = 800, linetype = "dashed", color = "red", size = 0.8) +
facet_wrap(~Strategy, ncol = 1, scales = "free_y") +
labs(title = "ICER Distribution Histograms",
subtitle = "Vertical lines: WTP thresholds ($100, $500, $800 per DALY)",
x = "ICER ($ per DALY averted)",
y = "Frequency") +
theme_minimal() +
scale_fill_manual(values = c("SMC" = "#377EB8", "Vaccine" = "#4DAF4A", "Combined" = "#E41A1C")) +
theme(legend.position = "none") +
scale_x_continuous(labels = comma, limits = c(0, 2000))
print(p_hist)
Figures Summary
# Clinical cases plot
clinical_data <- data.frame(
Strategy = strategy_names,
Clinical_Cases = sapply(results_list, function(x) x$total_clinical)
)
p1 <- ggplot(clinical_data, aes(x = reorder(Strategy, -Clinical_Cases),
y = Clinical_Cases, fill = Strategy)) +
geom_bar(stat = "identity", width = 0.7) +
geom_text(aes(label = format_number(Clinical_Cases)),
vjust = -0.5, size = 3.5) +
labs(title = "Clinical Cases by Strategy (Undiscounted)",
subtitle = "10-year horizon, 10,000 children",
x = "", y = "Clinical Cases") +
theme_minimal() +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_fill_brewer(palette = "Set2")
print(p1)
# DALYs averted plot
dalys_data <- ce_table[2:4, ]
p2 <- ggplot(dalys_data, aes(x = reorder(Strategy, -DALYs_Averted),
y = DALYs_Averted, fill = Strategy)) +
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 Strategy (Discounted at 3%)",
subtitle = "Health outcomes discounted at 3% per annum",
x = "", y = "DALYs Averted") +
theme_minimal() +
theme(legend.position = "none") +
scale_fill_brewer(palette = "Set2")
print(p2)
# ICER comparison plot
icer_data_plot <- ce_table[2:4, ]
icer_data_plot$CE_Threshold <- ifelse(icer_data_plot$ICER < 800, "Cost-Effective", "Not Cost-Effective")
p3 <- ggplot(icer_data_plot, aes(x = reorder(Strategy, ICER), y = ICER, fill = CE_Threshold)) +
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))),
hjust = -0.1, size = 3.5) +
labs(title = "Cost-Effectiveness Comparison",
subtitle = "WTP threshold: $800 per DALY averted (red dashed line)",
x = "", y = "ICER ($ per DALY averted)") +
theme_minimal() +
theme(legend.position = "bottom") +
scale_fill_manual(values = c("Cost-Effective" = "#4DAF4A", "Not Cost-Effective" = "#E41A1C")) +
coord_flip()
print(p3)
# Geographic distribution plot
p4 <- ggplot(geo_summary, aes(x = geography, y = DALYs_Discounted/1000, fill = Strategy)) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
labs(title = "Geographic Distribution of DALYs (Discounted)",
subtitle = "Discounted at 3% per annum",
x = "", y = "Discounted DALYs (thousands)") +
theme_minimal() +
scale_fill_brewer(palette = "Set1", name = "Intervention") +
theme(legend.position = "bottom")
print(p4)
BUDGET IMPACT ANALYSIS
# ===================================================================
# BUDGET IMPACT ANALYSIS (BIA) - CORRECTED VERSION
# ===================================================================
# Budget Impact Analysis for Malaria Interventions
# 5-year budget planning with national scale-up scenarios
cat("\n\n========================================\n")##
##
## ========================================## BUDGET IMPACT ANALYSIS## ========================================# ===================================================================
# PART 1: CREATE CE RESULTS
# ===================================================================
ce_results <- data.frame(
Strategy = strategy_names,
Clinical_Cases = sapply(results_list, function(x) x$total_clinical),
Severe_Cases = sapply(results_list, function(x) x$total_severe),
Deaths = sapply(results_list, function(x) x$total_deaths),
DALYs = sapply(results_list, function(x) x$total_dalys),
Costs = sapply(results_list, function(x) x$total_costs)
)
ce_results$DALYs_Averted <- ce_results$DALYs[1] - ce_results$DALYs
ce_results$Inc_Cost <- ce_results$Costs - ce_results$Costs[1]
ce_results$ICER <- ce_results$Inc_Cost / ce_results$DALYs_Averted
ce_results$ICER[ce_results$DALYs_Averted <= 0] <- NA
cat("✓ ce_results created successfully\n")## ✓ ce_results created successfully# ===================================================================
# PART 2: BUDGET IMPACT PARAMETERS
# ===================================================================
# National population parameters (example country: 5 million children under 5)
national_u5_population <- 5000000
cohort_size <- 10000
scale_factor <- national_u5_population / cohort_size
# Time horizon for budget impact (5 years)
bia_horizon_years <- 5
annual_discount_rate <- 0.03
# Target coverage targets by year (gradual scale-up)
coverage_targets <- data.frame(
Year = 1:5,
Vaccine_Urban = c(0.65, 0.70, 0.75, 0.80, 0.85),
Vaccine_Rural = c(0.40, 0.50, 0.60, 0.70, 0.75),
SMC_Urban = c(0.55, 0.60, 0.65, 0.70, 0.75),
SMC_Rural = c(0.45, 0.55, 0.65, 0.70, 0.75),
ITN_Urban = c(0.70, 0.75, 0.80, 0.85, 0.90),
ITN_Rural = c(0.45, 0.55, 0.65, 0.70, 0.75)
)
# Population distribution
urban_pop_u5 <- national_u5_population * 0.30 # 30% urban
rural_pop_u5 <- national_u5_population * 0.70 # 70% rural
# ===================================================================
# PART 3: CALCULATE STRATEGY COSTS (5-YEAR, UNDISCOUNTED)
# ===================================================================
# Function to calculate annual costs for a strategy
calculate_strategy_costs <- function(strategy, year) {
# Get coverage for this year
vaccine_cov_urban <- coverage_targets$Vaccine_Urban[year]
vaccine_cov_rural <- coverage_targets$Vaccine_Rural[year]
smc_cov_urban <- coverage_targets$SMC_Urban[year]
smc_cov_rural <- coverage_targets$SMC_Rural[year]
itn_cov_urban <- coverage_targets$ITN_Urban[year]
itn_cov_rural <- coverage_targets$ITN_Rural[year]
# ITN costs (annual distribution, 2-year lifespan)
itn_cost_per_net <- 2.50
itn_annual_cost_urban <- urban_pop_u5 * itn_cost_per_net * itn_cov_urban / 2
itn_annual_cost_rural <- rural_pop_u5 * itn_cost_per_net * itn_cov_rural / 2
itn_total <- itn_annual_cost_urban + itn_annual_cost_rural
if(strategy == "baseline") {
return(itn_total)
}
# SMC costs (annual, children 2-5 years = 27% of U5 population)
smc_eligible_u5_prop <- 0.27
smc_cost_per_child <- 6.00
smc_annual_cost_urban <- urban_pop_u5 * smc_eligible_u5_prop * smc_cost_per_child * smc_cov_urban
smc_annual_cost_rural <- rural_pop_u5 * smc_eligible_u5_prop * smc_cost_per_child * smc_cov_rural
smc_total <- smc_annual_cost_urban + smc_annual_cost_rural
# Vaccine costs (annual, children 0-24 months = 16% of U5 population)
vaccine_eligible_u5_prop <- 0.16
vaccine_cost_per_child <- 26.00
vaccine_annual_cost_urban <- urban_pop_u5 * vaccine_eligible_u5_prop * vaccine_cost_per_child * vaccine_cov_urban
vaccine_annual_cost_rural <- rural_pop_u5 * vaccine_eligible_u5_prop * vaccine_cost_per_child * vaccine_cov_rural
vaccine_total <- vaccine_annual_cost_urban + vaccine_annual_cost_rural
# Delivery cost multipliers
rural_delivery_multiplier <- 1.5
if(strategy == "smc") {
total <- itn_total + smc_total
} else if(strategy == "vaccine") {
total <- itn_total + vaccine_total
} else if(strategy == "combined") {
total <- itn_total + smc_total + (vaccine_total * rural_delivery_multiplier)
} else {
total <- itn_total
}
return(total)
}
# ===================================================================
# PART 4: CALCULATE ANNUAL AND CUMULATIVE COSTS
# ===================================================================
# Create budget impact data frame
bia_data <- data.frame(
Year = 1:bia_horizon_years,
Baseline_Cost = NA,
SMC_Cost = NA,
Vaccine_Cost = NA,
Combined_Cost = NA,
Incremental_SMC = NA,
Incremental_Vaccine = NA,
Incremental_Combined = NA
)
for(year in 1:bia_horizon_years) {
bia_data$Baseline_Cost[year] <- calculate_strategy_costs("baseline", year)
bia_data$SMC_Cost[year] <- calculate_strategy_costs("smc", year)
bia_data$Vaccine_Cost[year] <- calculate_strategy_costs("vaccine", year)
bia_data$Combined_Cost[year] <- calculate_strategy_costs("combined", year)
bia_data$Incremental_SMC[year] <- bia_data$SMC_Cost[year] - bia_data$Baseline_Cost[year]
bia_data$Incremental_Vaccine[year] <- bia_data$Vaccine_Cost[year] - bia_data$Baseline_Cost[year]
bia_data$Incremental_Combined[year] <- bia_data$Combined_Cost[year] - bia_data$Baseline_Cost[year]
}
# Calculate discounted costs (3% annually)
bia_data$Discount_Factor <- 1 / (1 + annual_discount_rate)^(bia_data$Year - 1)
for(col in c("Baseline_Cost", "SMC_Cost", "Vaccine_Cost", "Combined_Cost",
"Incremental_SMC", "Incremental_Vaccine", "Incremental_Combined")) {
bia_data[paste0(col, "_Discounted")] <- bia_data[[col]] * bia_data$Discount_Factor
}
# Calculate cumulative totals
cumulative_bia <- data.frame(
Strategy = c("Baseline (ITNs)", "SMC + ITNs", "Vaccine + ITNs", "Combined"),
Total_Undiscounted = c(
sum(bia_data$Baseline_Cost),
sum(bia_data$SMC_Cost),
sum(bia_data$Vaccine_Cost),
sum(bia_data$Combined_Cost)
),
Total_Discounted = c(
sum(bia_data$Baseline_Cost_Discounted),
sum(bia_data$SMC_Cost_Discounted),
sum(bia_data$Vaccine_Cost_Discounted),
sum(bia_data$Combined_Cost_Discounted)
)
)
cumulative_bia$Incremental_Undiscounted <- cumulative_bia$Total_Undiscounted - cumulative_bia$Total_Undiscounted[1]
cumulative_bia$Incremental_Discounted <- cumulative_bia$Total_Discounted - cumulative_bia$Total_Discounted[1]
# Format numbers for display
format_number <- function(x) {
format(round(x, 0), big.mark = ",", scientific = FALSE)
}
# ===================================================================
# PART 5: ANNUAL BUDGET TABLE
# ===================================================================
cat("\n\n========================================\n")##
##
## ========================================## BUDGET IMPACT ANALYSIS RESULTS## ========================================## Annual Budget Requirements (5-Year Scale-Up, National Level)## -----------------------------------------------------------# Format annual budget table
annual_budget_table <- bia_data %>%
mutate(
Year = Year,
Baseline = paste0("$", format_number(Baseline_Cost / 1e6), "M"),
SMC = paste0("$", format_number(SMC_Cost / 1e6), "M"),
Vaccine = paste0("$", format_number(Vaccine_Cost / 1e6), "M"),
Combined = paste0("$", format_number(Combined_Cost / 1e6), "M"),
Incremental_Vaccine = paste0("$", format_number(Incremental_Vaccine / 1e6), "M"),
Incremental_Combined = paste0("$", format_number(Incremental_Combined / 1e6), "M")
) %>%
select(Year, Baseline, SMC, Vaccine, Combined, Incremental_Vaccine, Incremental_Combined)
print(annual_budget_table)## Year Baseline SMC Vaccine Combined Incremental_Vaccine Incremental_Combined
## 1 1 $3M $ 7M $13M $22M $10M $19M
## 2 2 $4M $ 8M $15M $26M $12M $22M
## 3 3 $4M $10M $18M $30M $13M $25M
## 4 4 $5M $10M $20M $33M $15M $28M
## 5 5 $5M $11M $21M $35M $16M $30M# ===================================================================
# PART 6: CUMULATIVE BUDGET TABLE
# ===================================================================
cat("\n\nCumulative 5-Year Budget Requirements\n")##
##
## Cumulative 5-Year Budget Requirements## ---------------------------------------cumulative_display <- cumulative_bia %>%
mutate(
Total_Undiscounted = paste0("$", format_number(Total_Undiscounted / 1e6), "M"),
Total_Discounted = paste0("$", format_number(Total_Discounted / 1e6), "M"),
Incremental_Undiscounted = paste0("$", format_number(Incremental_Undiscounted / 1e6), "M"),
Incremental_Discounted = paste0("$", format_number(Incremental_Discounted / 1e6), "M")
)
print(cumulative_display)## Strategy Total_Undiscounted Total_Discounted Incremental_Undiscounted
## 1 Baseline (ITNs) $ 21M $ 20M $ 0M
## 2 SMC + ITNs $ 47M $ 44M $ 25M
## 3 Vaccine + ITNs $ 87M $ 82M $ 66M
## 4 Combined $146M $137M $125M
## Incremental_Discounted
## 1 $ 0M
## 2 $ 24M
## 3 $ 62M
## 4 $117M# ===================================================================
# PART 7: COST PER DALY BY BUDGET SCENARIO (USING CE_RESULTS)
# ===================================================================
# Safely extract DALYs averted
if(exists("ce_results") && nrow(ce_results) >= 4) {
dalys_averted_vaccine <- ce_results$DALYs_Averted[3]
dalys_averted_combined <- ce_results$DALYs_Averted[4]
# Scale to national level
dalys_averted_national_vaccine <- dalys_averted_vaccine * scale_factor
dalys_averted_national_combined <- dalys_averted_combined * scale_factor
# Calculate cost per DALY at national level
cost_per_daly_vaccine <- cumulative_bia$Total_Discounted[3] / dalys_averted_national_vaccine
cost_per_daly_combined <- cumulative_bia$Total_Discounted[4] / dalys_averted_national_combined
} else {
# Default values if ce_results not available
dalys_averted_national_vaccine <- 13780 * scale_factor / 10000
dalys_averted_national_combined <- 21950 * scale_factor / 10000
cost_per_daly_vaccine <- cumulative_bia$Total_Discounted[3] / dalys_averted_national_vaccine
cost_per_daly_combined <- cumulative_bia$Total_Discounted[4] / dalys_averted_national_combined
}
budget_scenarios <- data.frame(
Scenario = c("Base Case", "Optimistic (20% lower costs)",
"Pessimistic (20% higher costs)", "High Coverage (90% target)"),
Vaccine_Cost_Per_DALY = c(
cost_per_daly_vaccine,
cost_per_daly_vaccine * 0.8,
cost_per_daly_vaccine * 1.2,
cost_per_daly_vaccine * 0.85
),
Combined_Cost_Per_DALY = c(
cost_per_daly_combined,
cost_per_daly_combined * 0.8,
cost_per_daly_combined * 1.2,
cost_per_daly_combined * 0.85
),
Vaccine_Budget_M = c(
cumulative_bia$Total_Discounted[3] / 1e6,
cumulative_bia$Total_Discounted[3] / 1e6 * 0.8,
cumulative_bia$Total_Discounted[3] / 1e6 * 1.2,
cumulative_bia$Total_Discounted[3] / 1e6 * 1.1
),
Combined_Budget_M = c(
cumulative_bia$Total_Discounted[4] / 1e6,
cumulative_bia$Total_Discounted[4] / 1e6 * 0.8,
cumulative_bia$Total_Discounted[4] / 1e6 * 1.2,
cumulative_bia$Total_Discounted[4] / 1e6 * 1.1
)
)
cat("\n\nBudget Scenario Analysis\n")##
##
## Budget Scenario Analysis## ========================## Scenario Vaccine_Cost_Per_DALY Combined_Cost_Per_DALY
## 1 Base Case 133.8743 140.8530
## 2 Optimistic (20% lower costs) 107.0994 112.6824
## 3 Pessimistic (20% higher costs) 160.6491 169.0236
## 4 High Coverage (90% target) 113.7931 119.7250
## Vaccine_Budget_M Combined_Budget_M
## 1 81.89784 136.8507
## 2 65.51827 109.4805
## 3 98.27740 164.2208
## 4 90.08762 150.5358# ===================================================================
# PART 8: FINANCIAL SUSTAINABILITY ANALYSIS
# ===================================================================
# Calculate percentage of health budget
assumed_health_budget_million <- 5000 # $5 billion annual health budget
health_budget_share <- cumulative_bia$Total_Discounted[4] / 1e6 / (assumed_health_budget_million * bia_horizon_years) * 100
# Calculate cost per child per year
cost_per_child_vaccine <- cumulative_bia$Total_Discounted[3] / national_u5_population / bia_horizon_years
cost_per_child_combined <- cumulative_bia$Total_Discounted[4] / national_u5_population / bia_horizon_years
cat("\n\nFinancial Sustainability Indicators\n")##
##
## Financial Sustainability Indicators## =====================================## Cost per child per year (Vaccine strategy): $3.28## Cost per child per year (Combined strategy): $5.47## Share of annual health budget (Combined): 0.55%if(exists("dalys_averted_national_vaccine")) {
cat(sprintf("Total DALYs averted nationally (Vaccine): %s\n", format_number(dalys_averted_national_vaccine)))
cat(sprintf("Total DALYs averted nationally (Combined): %s\n", format_number(dalys_averted_national_combined)))
}## Total DALYs averted nationally (Vaccine): 611,752
## Total DALYs averted nationally (Combined): 971,585## Cost per DALY averted (Vaccine): $134## Cost per DALY averted (Combined): $141# ===================================================================
# PART 9: FUNDING SCENARIOS
# ===================================================================
funding_scenarios <- data.frame(
Scenario = c("Current Funding Only", "Gavi Support Added", "Domestic + Gavi", "Full Domestic"),
Annual_Gap_Vaccine_M = c(
max(0, cumulative_bia$Total_Discounted[3]/bia_horizon_years/1e6 - 10),
max(0, cumulative_bia$Total_Discounted[3]/bia_horizon_years/1e6 - 25),
max(0, cumulative_bia$Total_Discounted[3]/bia_horizon_years/1e6 - 40),
cumulative_bia$Total_Discounted[3]/bia_horizon_years/1e6
),
Annual_Gap_Combined_M = c(
max(0, cumulative_bia$Total_Discounted[4]/bia_horizon_years/1e6 - 10),
max(0, cumulative_bia$Total_Discounted[4]/bia_horizon_years/1e6 - 25),
max(0, cumulative_bia$Total_Discounted[4]/bia_horizon_years/1e6 - 40),
cumulative_bia$Total_Discounted[4]/bia_horizon_years/1e6
)
)
cat("\n\nFunding Gap Analysis (Annual)\n")##
##
## Funding Gap Analysis (Annual)## ==============================## Scenario Annual_Gap_Vaccine_M Annual_Gap_Combined_M
## 1 Current Funding Only 6.379567 17.370137
## 2 Gavi Support Added 0.000000 2.370137
## 3 Domestic + Gavi 0.000000 0.000000
## 4 Full Domestic 16.379567 27.370137# ===================================================================
# PART 10: BUDGET IMPACT VISUALIZATIONS
# ===================================================================
# Plot 1: Annual budget comparison by strategy
p_budget_annual <- ggplot(bia_data, aes(x = Year)) +
geom_line(aes(y = Baseline_Cost / 1e6, color = "Baseline (ITNs)"), size = 1.2) +
geom_line(aes(y = SMC_Cost / 1e6, color = "SMC + ITNs"), size = 1.2) +
geom_line(aes(y = Vaccine_Cost / 1e6, color = "Vaccine + ITNs"), size = 1.2) +
geom_line(aes(y = Combined_Cost / 1e6, color = "Combined"), size = 1.2) +
labs(title = "Annual Budget Requirements by Strategy (Undiscounted)",
subtitle = "5-year gradual scale-up to target coverage | National: 5M children",
x = "Year", y = "Annual Cost (Millions USD)") +
theme_minimal() +
scale_color_manual(values = c("Baseline (ITNs)" = "gray50",
"SMC + ITNs" = "#377EB8",
"Vaccine + ITNs" = "#4DAF4A",
"Combined" = "#E41A1C"),
name = "Strategy") +
theme(legend.position = "bottom")
# Plot 2: Incremental budget over baseline
p_budget_incremental <- ggplot(bia_data, aes(x = Year)) +
geom_area(aes(y = Incremental_Vaccine / 1e6, fill = "Vaccine + ITNs"), alpha = 0.6) +
geom_area(aes(y = Incremental_Combined / 1e6, fill = "Combined"), alpha = 0.6) +
geom_line(aes(y = Incremental_Vaccine / 1e6, color = "Vaccine + ITNs"), size = 1.2) +
geom_line(aes(y = Incremental_Combined / 1e6, color = "Combined"), size = 1.2) +
labs(title = "Incremental Budget Requirements vs Baseline",
subtitle = "Additional funding needed above current ITN spending",
x = "Year", y = "Incremental Annual Cost (Millions USD)") +
theme_minimal() +
scale_fill_manual(values = c("Vaccine + ITNs" = "#4DAF4A", "Combined" = "#E41A1C"),
name = "Strategy") +
scale_color_manual(values = c("Vaccine + ITNs" = "#4DAF4A", "Combined" = "#E41A1C"),
name = "Strategy", guide = "none") +
theme(legend.position = "bottom")
# Plot 3: Cumulative budget bar chart
p_budget_cumulative <- ggplot(cumulative_bia, aes(x = Strategy, y = Total_Discounted / 1e6, fill = Strategy)) +
geom_bar(stat = "identity", width = 0.7) +
geom_text(aes(label = paste0("$", format_number(Total_Discounted / 1e6), "M")),
vjust = -0.5, size = 3.5) +
labs(title = "Cumulative 5-Year Budget Requirements (Discounted at 3%)",
subtitle = "National scale: 5 million children under 5",
x = "", y = "Total Cost (Millions USD)") +
theme_minimal() +
scale_fill_manual(values = c("Baseline (ITNs)" = "gray50",
"SMC + ITNs" = "#377EB8",
"Vaccine + ITNs" = "#4DAF4A",
"Combined" = "#E41A1C"),
name = "Strategy") +
theme(legend.position = "none")
# Plot 4: Coverage scale-up trajectories
coverage_long <- coverage_targets %>%
pivot_longer(cols = -Year, names_to = "Setting", values_to = "Coverage") %>%
separate(Setting, into = c("Intervention", "Geography"), sep = "_")
p_coverage_scaleup <- ggplot(coverage_long, aes(x = Year, y = Coverage, color = Intervention, linetype = Geography)) +
geom_line(size = 1.2) +
geom_point(size = 3) +
labs(title = "Coverage Scale-Up Trajectories",
subtitle = "Gradual increase to target coverage over 5 years",
x = "Year", y = "Coverage (%)") +
theme_minimal() +
scale_color_manual(values = c("Vaccine" = "#4DAF4A", "SMC" = "#377EB8", "ITN" = "gray50"),
name = "Intervention") +
scale_linetype_manual(values = c("Urban" = "solid", "Rural" = "dashed"),
name = "Geography") +
scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(0, 1)) +
theme(legend.position = "bottom")
# Save budget impact plots (only if not in R Markdown rendering)
if(!interactive() || !exists("knitr")) {
ggsave("BIA_Annual_Budget.png", p_budget_annual, width = 10, height = 6)
ggsave("BIA_Incremental_Budget.png", p_budget_incremental, width = 10, height = 6)
ggsave("BIA_Cumulative_Budget.png", p_budget_cumulative, width = 10, height = 7)
ggsave("BIA_Coverage_ScaleUp.png", p_coverage_scaleup, width = 10, height = 6)
}
# Display plots
print(p_budget_annual)
# ===================================================================
# PART 11: PRINT SUMMARY
# ===================================================================
cat("\n\n========================================\n")##
##
## ========================================## BUDGET IMPACT ANALYSIS COMPLETE## ========================================##
## Key Budget Findings:## --------------------cat(sprintf("• Vaccine + ITNs 5-year cost: $%.1f million (discounted)\n", cumulative_bia$Total_Discounted[3] / 1e6))## • Vaccine + ITNs 5-year cost: $81.9 million (discounted)cat(sprintf("• Combined strategy 5-year cost: $%.1f million (discounted)\n", cumulative_bia$Total_Discounted[4] / 1e6))## • Combined strategy 5-year cost: $136.9 million (discounted)cat(sprintf("• Additional funding needed for vaccine: $%.1f million\n", cumulative_bia$Incremental_Discounted[3] / 1e6))## • Additional funding needed for vaccine: $62.1 millioncat(sprintf("• Additional funding needed for combined: $%.1f million\n", cumulative_bia$Incremental_Discounted[4] / 1e6))## • Additional funding needed for combined: $117.1 million## • Cost per child per year (vaccine): $3.28## • Cost per DALY averted (vaccine): $134##
## Policy Implications:## --------------------if(cost_per_daly_vaccine < 800) {
cat("✓ Vaccine strategy is cost-effective and financially feasible\n")
}## ✓ Vaccine strategy is cost-effective and financially feasibleif(cumulative_bia$Incremental_Discounted[4] / 1e6 > 1000) {
cat("⚠️ Combined strategy requires substantial additional funding (>$1B)\n")
cat(" Consider phased implementation: start with vaccine, add SMC later\n")
}
cat("\nGenerated budget impact plots:\n")##
## Generated budget impact plots:## - BIA_Annual_Budget.png## - BIA_Incremental_Budget.png## - BIA_Cumulative_Budget.png## - BIA_Coverage_ScaleUp.png