Equity Analysis of Empower Health vs Usual Care - Kenya CVD Model_Latest
James Oguta
2026-05-04
- Introduction
This analysis evaluates the distributional impact of Empower Health intervention compared to usual care across wealth quintiles, regions, and sex in Kenya. The analysis includes:
Health benefits: DALYs averted (disability-adjusted life years) - discounted at 3%
Life years gained: Years of life saved (mortality benefit only) - NOT discounted (undiscounted)
Cost burden: Incremental costs by equity stratifiers - discounted at 3%
Concentration curves for health benefits, life years, and costs
Cost-effectiveness ICERs by equity metric
Kenya GDP per capita = $2,200 (WTP threshold)
Important Note on Concentration Curve Interpretation: Curve below diagonal = Positive C-index (+)
For Health Benefits/Life Years: PRO-RICH (richer gain more benefit) - INEQUITABLE
For Costs: PROGRESSIVE (richer pay more) - EQUITABLE
Curve above diagonal = Negative C-index (-)
For Health Benefits/Life Years: PRO-POOR (poorer gain more benefit) - EQUITABLE
For Costs: REGRESSIVE (poorer pay more) - INEQUITABLE
Clear Workspace
- Load Required Libraries
- Define Parameters and Load Data
# Discount rate (3% per year for costs and DALYs only)
DISCOUNT_RATE <- 0.03
# Kenya GDP per capita (WTP threshold)
WTP_THRESHOLD <- 2090
# Health state definitions
death_states <- c("D")
morbidity_states <- c("MI", "ST", "PS", "PM")
alive_states <- c("NE", "MI", "ST", "PS", "PM") # States where individual is alive
# Define file paths
base_path <- "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/data/"
file_path_Emp <- "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/sim_results/Empower Health/sim_1_data.rds"
file_path_UC <- "/Users/jamesoguta/Documents/James Oguta/PhD Thesis_James Oguta/James Oguta Model/Kenya CVD Model/KenyaCVDModelFinal/sim_results/Usual Care/sim_1_data.rds"
# Load simulation data
cat("Loading simulation data...\n")## Loading simulation data...sim_data_Emp <- readRDS(file_path_Emp)
sim_data_UC <- readRDS(file_path_UC)
# Extract matrices for Empower Health
m_States_Emp <- sim_data_Emp$m_States
m_Costs_Emp <- sim_data_Emp$m_Costs
m_Effs_Emp <- sim_data_Emp$m_Effs
# Extract matrices for Usual Care
m_States_UC <- sim_data_UC$m_States
m_Costs_UC <- sim_data_UC$m_Costs
m_Effs_UC <- sim_data_UC$m_Effs
# Load individual characteristics
load(paste0(base_path, "m_indi_features_3.rda"))
cat("✓ All datasets loaded successfully\n")## ✓ All datasets loaded successfully## Individuals: 100000## Cycles: 61- Apply Discounting (Costs and DALYs only - NOT Life Years)
# Create discount factors
n_cycles <- ncol(m_States_UC)
cycles <- 0:(n_cycles - 1)
discount_factors <- 1 / (1 + DISCOUNT_RATE)^cycles
# Discounting function
discount_matrix <- function(matrix_data, discount_factors) {
sweep(matrix_data, 2, discount_factors, FUN = "*")
}
# Apply discounting to costs and DALYs ONLY
# Life years will remain undiscounted
cat("Applying discounting to costs and DALYs (3% per year)...\n")## Applying discounting to costs and DALYs (3% per year)...## Note: Life years are NOT discountedm_Costs_Emp_disc <- discount_matrix(m_Costs_Emp, discount_factors)
m_Effs_Emp_disc <- discount_matrix(m_Effs_Emp, discount_factors)
m_Costs_UC_disc <- discount_matrix(m_Costs_UC, discount_factors)
m_Effs_UC_disc <- discount_matrix(m_Effs_UC, discount_factors)
cat("✓ Discounting applied\n")## ✓ Discounting applied- Process Individual Characteristics
features_df <- as.data.frame(m_indi_features_3, stringsAsFactors = FALSE)
# Clean wealth variable
wealth_clean <- gsub(".*?([0-9]+).*", "\\1", features_df$wealth)
features_df$wealth_quintile_num <- as.numeric(wealth_clean)
features_df$wealth_quintile <- factor(features_df$wealth_quintile_num,
levels = 1:5,
labels = c("Q1_Poorest", "Q2", "Q3", "Q4", "Q5_Richest"))
# Process other variables
features_df$region <- as.factor(features_df$region)
features_df$female <- as.numeric(features_df$female)
features_df$male <- as.numeric(features_df$male)
features_df$sex <- ifelse(features_df$female == 1, "Female",
ifelse(features_df$male == 1, "Male", "Unknown"))
features_df$sex <- as.factor(features_df$sex)
features_df$individual_id <- rownames(m_States_UC)
# Display summary
cat("\n=== INDIVIDUAL CHARACTERISTICS ===\n")##
## === INDIVIDUAL CHARACTERISTICS ===## Wealth distribution:##
## Q1_Poorest Q2 Q3 Q4 Q5_Richest
## 18251 23373 23267 22430 12679##
## Sex distribution:##
## Female Male
## 47597 52403##
## Regional distribution:##
## Central Coast Eastern Nairobi North Eastern
## 15468 4444 19148 9419 2842
## Nyanza Rift Valley Western
## 10214 24660 13805- Process Intervention Arms (Life Years NOT Discounted)
process_arm <- function(states_matrix, costs_matrix_disc, dalys_matrix_disc, arm_name) {
n_ind <- nrow(states_matrix)
n_cycles <- ncol(states_matrix)
cat(sprintf("Processing %s: %d individuals\n", arm_name, n_ind))
# Calculate total discounted costs and DALYs
results <- data.frame(
individual_id = rownames(states_matrix),
arm = arm_name,
total_costs = rowSums(costs_matrix_disc, na.rm = TRUE),
total_dalys = rowSums(dalys_matrix_disc, na.rm = TRUE)
)
# Calculate LIFE YEARS (UNDISCOUNTED - simple count of alive cycles)
is_alive_matrix <- matrix(states_matrix %in% alive_states, nrow = n_ind, ncol = n_cycles)
life_years <- rowSums(is_alive_matrix, na.rm = TRUE)
# Detect death cycle
is_death_matrix <- matrix(states_matrix %in% death_states, nrow = n_ind, ncol = n_cycles)
death_cycle <- apply(is_death_matrix, 1, function(x) {
idx <- which(x)[1]
ifelse(is.na(idx), n_cycles + 1, idx)
})
# Detect morbidity
is_morbidity_matrix <- matrix(states_matrix %in% morbidity_states, nrow = n_ind, ncol = n_cycles)
morbidity_count <- rowSums(is_morbidity_matrix, na.rm = TRUE)
# Store results
results$death_cycle <- death_cycle
results$cycles_survived <- pmin(death_cycle, n_cycles)
results$life_years <- life_years
results$alive_at_end <- death_cycle > n_cycles
results$had_morbidity <- morbidity_count > 0
return(results)
}
# Process both arms
cat("\n=== PROCESSING INTERVENTION ARMS ===\n")##
## === PROCESSING INTERVENTION ARMS ===## Note: Life years are undiscounted, costs and DALYs are discounted at 3%## Processing Empower: 100000 individuals## Processing UsualCare: 100000 individuals- Create Comparison Dataset
df_compare <- results_Emp %>%
select(individual_id,
total_costs_Emp = total_costs,
total_dalys_Emp = total_dalys,
life_years_Emp = life_years,
cycles_survived_Emp = cycles_survived,
alive_at_end_Emp = alive_at_end,
had_morbidity_Emp = had_morbidity) %>%
left_join(
results_UC %>% select(individual_id,
total_costs_UC = total_costs,
total_dalys_UC = total_dalys,
life_years_UC = life_years,
cycles_survived_UC = cycles_survived,
alive_at_end_UC = alive_at_end,
had_morbidity_UC = had_morbidity),
by = "individual_id"
) %>%
left_join(features_df, by = "individual_id") %>%
mutate(
delta_costs = total_costs_Emp - total_costs_UC,
delta_dalys = total_dalys_Emp - total_dalys_UC,
dalys_averted = -delta_dalys,
delta_life_years = life_years_Emp - life_years_UC,
life_years_gained = delta_life_years,
delta_survival = cycles_survived_Emp - cycles_survived_UC,
delta_mortality = (1 - alive_at_end_Emp) - (1 - alive_at_end_UC),
delta_morbidity = had_morbidity_Emp - had_morbidity_UC
)
cat(sprintf("Comparison dataframe: %d rows\n", nrow(df_compare)))## Comparison dataframe: 100000 rows- Subgroup Analysis Function (with Life Years)
subgroup_analysis <- function(data, group_var) {
groups <- unique(data[[group_var]])
groups <- groups[!is.na(groups)]
results <- list()
for (g in groups) {
sub_data <- data[data[[group_var]] == g, ]
n <- nrow(sub_data)
if (n == 0) next
# Mean outcomes
mean_delta_costs <- mean(sub_data$delta_costs, na.rm = TRUE)
mean_dalys_averted <- mean(sub_data$dalys_averted, na.rm = TRUE)
mean_life_years_gained <- mean(sub_data$life_years_gained, na.rm = TRUE)
mean_delta_mortality <- mean(sub_data$delta_mortality, na.rm = TRUE)
# Standard errors
se_costs <- sd(sub_data$delta_costs, na.rm = TRUE) / sqrt(n)
se_dalys <- sd(sub_data$dalys_averted, na.rm = TRUE) / sqrt(n)
se_life <- sd(sub_data$life_years_gained, na.rm = TRUE) / sqrt(n)
# ICERs
if (mean_dalys_averted > 0) {
icer_daly <- mean_delta_costs / mean_dalys_averted
icer_daly_text <- sprintf("$%.0f", icer_daly)
cost_effective_daly <- icer_daly < WTP_THRESHOLD
} else {
icer_daly <- NA
icer_daly_text <- "NA"
cost_effective_daly <- FALSE
}
if (mean_life_years_gained > 0) {
icer_life <- mean_delta_costs / mean_life_years_gained
icer_life_text <- sprintf("$%.0f", icer_life)
cost_effective_life <- icer_life < WTP_THRESHOLD
} else {
icer_life <- NA
icer_life_text <- "NA"
cost_effective_life <- FALSE
}
# Confidence intervals
costs_ci <- c(mean_delta_costs - 1.96 * se_costs, mean_delta_costs + 1.96 * se_costs)
dalys_ci <- c(mean_dalys_averted - 1.96 * se_dalys, mean_dalys_averted + 1.96 * se_dalys)
life_ci <- c(mean_life_years_gained - 1.96 * se_life, mean_life_years_gained + 1.96 * se_life)
results[[as.character(g)]] <- data.frame(
subgroup = as.character(g),
n = n,
delta_costs = mean_delta_costs,
costs_ci_lower = costs_ci[1],
costs_ci_upper = costs_ci[2],
dalys_averted = mean_dalys_averted,
dalys_ci_lower = dalys_ci[1],
dalys_ci_upper = dalys_ci[2],
icer_daly = icer_daly,
icer_daly_text = icer_daly_text,
cost_effective_daly = cost_effective_daly,
life_years_gained = mean_life_years_gained,
life_ci_lower = life_ci[1],
life_ci_upper = life_ci[2],
icer_life = icer_life,
icer_life_text = icer_life_text,
cost_effective_life = cost_effective_life,
delta_mortality = mean_delta_mortality
)
}
return(do.call(rbind, results))
}- Run Subgroup Analyses
# Wealth quintile analysis
wealth_results <- subgroup_analysis(df_compare, "wealth_quintile")
wealth_results$subgroup <- factor(wealth_results$subgroup,
levels = c("Q1_Poorest", "Q2", "Q3", "Q4", "Q5_Richest"))
wealth_results <- wealth_results %>% arrange(subgroup)
cat("\n=== WEALTH QUINTILE RESULTS ===\n")##
## === WEALTH QUINTILE RESULTS ===print(wealth_results %>% select(subgroup, n, dalys_averted, life_years_gained, delta_costs,
icer_daly_text, icer_life_text))## subgroup n dalys_averted life_years_gained delta_costs
## Q1_Poorest Q1_Poorest 18251 0.09988889 0.3297901 384.5282
## Q2 Q2 23373 0.11187154 0.3587045 432.2140
## Q3 Q3 23267 0.11779849 0.3881893 480.9004
## Q4 Q4 22430 0.12649013 0.4216228 499.6531
## Q5_Richest Q5_Richest 12679 0.13123331 0.4737755 559.2831
## icer_daly_text icer_life_text
## Q1_Poorest $3850 $1166
## Q2 $3863 $1205
## Q3 $4082 $1239
## Q4 $3950 $1185
## Q5_Richest $4262 $1180##
## === SEX RESULTS ===print(sex_results %>% select(subgroup, n, dalys_averted, life_years_gained, delta_costs,
icer_daly_text, icer_life_text))## subgroup n dalys_averted life_years_gained delta_costs
## Female Female 47597 0.1150564 0.3968527 524.7948
## Male Male 52403 0.1183788 0.3818484 412.7432
## icer_daly_text icer_life_text
## Female $4561 $1322
## Male $3487 $1081# Regional analysis
region_counts <- table(df_compare$region)
region_percent <- prop.table(region_counts) * 100
major_regions <- names(region_counts)[region_percent >= 1]
df_regions_major <- df_compare %>% filter(region %in% major_regions)
region_results <- subgroup_analysis(df_regions_major, "region")
region_results <- region_results %>% arrange(desc(dalys_averted))
cat("\n=== REGIONAL RESULTS (Top 10) ===\n")##
## === REGIONAL RESULTS (Top 10) ===print(head(region_results %>% select(subgroup, n, dalys_averted, life_years_gained,
icer_daly_text), 10))## subgroup n dalys_averted life_years_gained
## Nairobi Nairobi 9419 0.13666225 0.4847648
## Coast Coast 4444 0.13272881 0.4520702
## Eastern Eastern 19148 0.13168074 0.4253708
## Central Central 15468 0.11966531 0.4016680
## Western Western 13805 0.11803885 0.3712423
## Nyanza Nyanza 10214 0.10238659 0.3379675
## Rift Valley Rift Valley 24660 0.10155691 0.3489457
## North Eastern North Eastern 2842 0.08816788 0.2758621
## icer_daly_text
## Nairobi $4260
## Coast $3663
## Eastern $3548
## Central $4274
## Western $3330
## Nyanza $4714
## Rift Valley $4232
## North Eastern $4602- Overall Results
# Overall means
mean_cost_Emp <- mean(df_compare$total_costs_Emp, na.rm = TRUE)
mean_cost_UC <- mean(df_compare$total_costs_UC, na.rm = TRUE)
mean_daly_Emp <- mean(df_compare$total_dalys_Emp, na.rm = TRUE)
mean_daly_UC <- mean(df_compare$total_dalys_UC, na.rm = TRUE)
mean_life_Emp <- mean(df_compare$life_years_Emp, na.rm = TRUE)
mean_life_UC <- mean(df_compare$life_years_UC, na.rm = TRUE)
delta_costs <- mean_cost_Emp - mean_cost_UC
delta_dalys <- mean_daly_Emp - mean_daly_UC
dalys_averted <- -delta_dalys
delta_life <- mean_life_Emp - mean_life_UC
life_years_gained <- delta_life
icer_daly <- delta_costs / dalys_averted
icer_life <- delta_costs / life_years_gained
cat("\n=== OVERALL RESULTS ===\n")##
## === OVERALL RESULTS ===cat(sprintf("Empower: Cost = $%.2f (discounted), DALYs = %.6f, Life Years = %.4f\n",
mean_cost_Emp, mean_daly_Emp, mean_life_Emp))## Empower: Cost = $5170.83 (discounted), DALYs = 7.384656, Life Years = 21.9610cat(sprintf("Usual Care: Cost = $%.2f, DALYs = %.6f, Life Years = %.4f\n",
mean_cost_UC, mean_daly_UC, mean_life_UC))## Usual Care: Cost = $4704.75, DALYs = 7.501453, Life Years = 21.5720##
## ICER (per DALY): $3990## ICER (per Life Year): $1198- Concentration Indices (DALYs, Life Years, and Costs)
# Function to calculate concentration curve with analytical standard error
concentration_curve_with_ci <- function(data, var_name, wealth_var = "wealth_quintile_num") {
valid <- complete.cases(data[[var_name]], data[[wealth_var]])
values <- data[[var_name]][valid]
wealth <- data[[wealth_var]][valid]
if (length(values) == 0) return(NULL)
# Sort by wealth (poorest to richest)
order_idx <- order(wealth)
values_sorted <- values[order_idx]
# Calculate cumulative shares
n <- length(values_sorted)
cum_pop <- (1:n) / n
cum_values <- cumsum(values_sorted) / sum(values_sorted)
# Calculate concentration index
mu <- mean(values_sorted)
if (mu == 0) return(NULL)
C <- (2 / (n * mu)) * sum(cum_pop * values_sorted) - 1
# Calculate analytical standard error for concentration index
# Using Kakwani et al. (1997) method
# Variance formula: var(C) = (1/n) * (sum((y/mu - 1)^2) - C^2)
y <- values_sorted
variance <- (1/n) * (sum((y/mu - 1)^2, na.rm = TRUE) - (C^2))
se_C <- sqrt(variance / n)
# 95% confidence interval
ci_lower <- C - 1.96 * se_C
ci_upper <- C + 1.96 * se_C
# Calculate standard errors for cumulative curve points
# Using the delta method for each point
cum_sum <- cumsum(y)
cum_sum_total <- cum_sum[n]
# Standard error for each point on the curve
se_curve <- numeric(n)
for (i in 1:n) {
# Approximate SE for cumulative proportion
p <- cum_pop[i]
cum_y <- cum_sum[i]
# Using formula: SE = sqrt(p*(1-p)/n) * (cum_y/cum_sum_total)
se_curve[i] <- sqrt(p * (1-p) / n) * (cum_y / cum_sum_total)
}
# Confidence bands for curve
curve_lower <- pmax(0, cum_values - 1.96 * se_curve)
curve_upper <- pmin(1, cum_values + 1.96 * se_curve)
return(list(
curve = data.frame(
Population = cum_pop,
Cumulative = cum_values,
Lower_CI = curve_lower,
Upper_CI = curve_upper
),
C_index = C,
C_ci_lower = ci_lower,
C_ci_upper = ci_upper,
se_C = se_C
))
}
# Generate concentration curves with CIs
health_curve <- concentration_curve_with_ci(df_compare, "dalys_averted", "wealth_quintile_num")
life_curve <- concentration_curve_with_ci(df_compare, "life_years_gained", "wealth_quintile_num")
costs_curve <- concentration_curve_with_ci(df_compare, "delta_costs", "wealth_quintile_num")
# Display results with CIs
cat("\n=== CONCENTRATION CURVE RESULTS WITH 95% CI ===\n")##
## === CONCENTRATION CURVE RESULTS WITH 95% CI ===if (!is.null(health_curve)) {
cat(sprintf("\nDALYs C-index: %.4f (95%% CI: %.4f - %.4f)\n",
health_curve$C_index, health_curve$C_ci_lower, health_curve$C_ci_upper))
cat(ifelse(health_curve$C_index > 0,
" → PRO-RICH (richer gain more DALY benefit)\n",
" → PRO-POOR (poorer gain more DALY benefit)\n"))
}##
## DALYs C-index: 0.0451 (95% CI: 0.0123 - 0.0780)
## → PRO-RICH (richer gain more DALY benefit)if (!is.null(life_curve)) {
cat(sprintf("\nLife Years C-index: %.4f (95%% CI: %.4f - %.4f)\n",
life_curve$C_index, life_curve$C_ci_lower, life_curve$C_ci_upper))
cat(ifelse(life_curve$C_index > 0,
" → PRO-RICH (richer gain more life years)\n",
" → PRO-POOR (poorer gain more life years)\n"))
}##
## Life Years C-index: 0.0622 (95% CI: 0.0278 - 0.0966)
## → PRO-RICH (richer gain more life years)if (!is.null(costs_curve)) {
cat(sprintf("\nCosts C-index: %.4f (95%% CI: %.4f - %.4f)\n",
costs_curve$C_index, costs_curve$C_ci_lower, costs_curve$C_ci_upper))
cat(ifelse(costs_curve$C_index > 0,
" → PROGRESSIVE (richer pay more) - EQUITABLE\n",
" → REGRESSIVE (poorer pay more) - INEQUITABLE\n"))
}##
## Costs C-index: 0.0646 (95% CI: 0.0437 - 0.0855)
## → PROGRESSIVE (richer pay more) - EQUITABLE- Figure 1: DALYs and Life Years Analysis with 95% CI Bands
# Overall means for reference lines
overall_dalys <- mean(df_compare$dalys_averted, na.rm = TRUE)
overall_life <- mean(df_compare$life_years_gained, na.rm = TRUE)
# Panel A: DALYs averted by wealth
panel_A <- ggplot(wealth_results, aes(x = dalys_averted, y = subgroup)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_vline(xintercept = overall_dalys, linetype = "dotted", color = "blue", size = 0.8) +
geom_point(size = 4, color = "steelblue") +
geom_errorbarh(aes(xmin = dalys_ci_lower, xmax = dalys_ci_upper), height = 0.2, size = 0.8, color = "steelblue") +
annotate("text", x = overall_dalys + 0.01, y = 5.2, label = sprintf("Mean = %.4f", overall_dalys), size = 3, color = "blue") +
labs(x = "\nDALYs Averted (discounted)", y = "Wealth Quintile\n",
title = "A) DALYs Averted by Wealth Quintile (3% Discount)",
subtitle = "Positive = beneficial | Blue dotted line = population mean") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
# Panel B: Life years gained by wealth
panel_B <- ggplot(wealth_results, aes(x = life_years_gained, y = subgroup)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_vline(xintercept = overall_life, linetype = "dotted", color = "blue", size = 0.8) +
geom_point(size = 4, color = "darkgreen") +
geom_errorbarh(aes(xmin = life_ci_lower, xmax = life_ci_upper), height = 0.2, size = 0.8, color = "darkgreen") +
annotate("text", x = overall_life + 0.01, y = 5.2, label = sprintf("Mean = %.4f", overall_life), size = 3, color = "blue") +
labs(x = "\nLife Years Gained (undiscounted)", y = "Wealth Quintile\n",
title = "B) Life Years Gained by Wealth Quintile (Undiscounted)",
subtitle = "Positive = beneficial | Blue dotted line = population mean") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
# Panel C: DALYs concentration curve with 95% CI band
if (!is.null(health_curve)) {
panel_C <- ggplot(health_curve$curve, aes(x = Population, y = Cumulative)) +
geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), alpha = 0.2, fill = "steelblue") +
geom_line(size = 1.2, color = "steelblue") +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50", size = 0.8) +
annotate("text", x = 0.7, y = 0.25,
label = sprintf("C-index = %.4f\n95%% CI [%.4f, %.4f]",
health_curve$C_index, health_curve$C_ci_lower, health_curve$C_ci_upper),
size = 3.5, fontface = "bold", color = "steelblue") +
annotate("text", x = 0.7, y = 0.12,
label = ifelse(health_curve$C_index > 0,
"Interpretation: PRO-RICH\n(Curve below diagonal = richer gain more DALY benefit)\n→ INEQUITABLE",
"Interpretation: PRO-POOR\n(Curve above diagonal = poorer gain more DALY benefit)\n→ EQUITABLE"),
size = 2.8, color = "steelblue", lineheight = 1.2) +
labs(x = "\nCumulative population (poorest to richest)",
y = "Cumulative DALYs averted (discounted)\n",
title = "C) DALYs Concentration Curve (3% Discount)",
subtitle = "Shaded area = 95% CI | Curve below diagonal (C>0) = PRO-RICH") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
} else {
panel_C <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
# Panel D: Life years concentration curve with 95% CI band
if (!is.null(life_curve)) {
panel_D <- ggplot(life_curve$curve, aes(x = Population, y = Cumulative)) +
geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), alpha = 0.2, fill = "darkgreen") +
geom_line(size = 1.2, color = "darkgreen") +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50", size = 0.8) +
annotate("text", x = 0.7, y = 0.25,
label = sprintf("C-index = %.4f\n95%% CI [%.4f, %.4f]",
life_curve$C_index, life_curve$C_ci_lower, life_curve$C_ci_upper),
size = 3.5, fontface = "bold", color = "darkgreen") +
annotate("text", x = 0.7, y = 0.12,
label = ifelse(life_curve$C_index > 0,
"Interpretation: PRO-RICH\n(Curve below diagonal = richer gain more life years)\n→ INEQUITABLE",
"Interpretation: PRO-POOR\n(Curve above diagonal = poorer gain more life years)\n→ EQUITABLE"),
size = 2.8, color = "darkgreen", lineheight = 1.2) +
labs(x = "\nCumulative population (poorest to richest)",
y = "Cumulative life years gained (undiscounted)\n",
title = "D) Life Years Concentration Curve (Undiscounted)",
subtitle = "Shaded area = 95% CI | Curve below diagonal (C>0) = PRO-RICH") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
} else {
panel_D <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
# Arrange all panels
grid.arrange(panel_A, panel_B, panel_C, panel_D, ncol = 2, nrow = 2,
top = grid::textGrob("Figure 1: DALYs (Discounted) and Life Years (Undiscounted) Equity Analysis",
gp = grid::gpar(fontsize = 14, fontface = "bold")))
13. Figure 2: Sex and Regional Analyses (4 Panels)
# Panel E: DALYs by sex
panel_E <- ggplot(sex_results, aes(x = subgroup, y = dalys_averted, fill = subgroup)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.6) +
geom_errorbar(aes(ymin = dalys_ci_lower, ymax = dalys_ci_upper), width = 0.2, size = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_hline(yintercept = overall_dalys, linetype = "dotted", color = "blue", size = 0.8) +
annotate("text", x = 2.5, y = overall_dalys + 0.01, label = sprintf("Mean = %.4f", overall_dalys), size = 3, color = "blue") +
labs(x = "\nSex", y = "DALYs Averted (discounted)\n",
title = "E) DALYs Averted by Sex (3% Discount)",
subtitle = "Blue dotted line = population mean") +
scale_fill_manual(values = c("Female" = "#FF6B6B", "Male" = "#4ECDC4")) +
theme_minimal(base_size = 11) +
theme(legend.position = "none", plot.title = element_text(face = "bold"))
# Panel F: Life years by sex
panel_F <- ggplot(sex_results, aes(x = subgroup, y = life_years_gained, fill = subgroup)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.6) +
geom_errorbar(aes(ymin = life_ci_lower, ymax = life_ci_upper), width = 0.2, size = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_hline(yintercept = overall_life, linetype = "dotted", color = "blue", size = 0.8) +
annotate("text", x = 2.5, y = overall_life + 0.01, label = sprintf("Mean = %.4f", overall_life), size = 3, color = "blue") +
labs(x = "\nSex", y = "Life Years Gained (undiscounted)\n",
title = "F) Life Years Gained by Sex (Undiscounted)",
subtitle = "Blue dotted line = population mean") +
scale_fill_manual(values = c("Female" = "#FF6B6B", "Male" = "#4ECDC4")) +
theme_minimal(base_size = 11) +
theme(legend.position = "none", plot.title = element_text(face = "bold"))
# Panel G: DALYs by region
if (nrow(region_results) >= 3) {
panel_G <- ggplot(region_results, aes(x = dalys_averted, y = reorder(subgroup, dalys_averted))) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_vline(xintercept = overall_dalys, linetype = "dotted", color = "blue", size = 0.8) +
geom_point(size = 3, color = "darkorange") +
geom_errorbarh(aes(xmin = dalys_ci_lower, xmax = dalys_ci_upper), height = 0.2, size = 0.8, color = "darkorange") +
annotate("text", x = overall_dalys + 0.01, y = nrow(region_results) + 0.5, label = sprintf("Mean = %.4f", overall_dalys), size = 3, color = "blue") +
labs(x = "\nDALYs Averted (discounted)", y = "Region\n",
title = "G) DALYs Averted by Region (3% Discount)",
subtitle = "Blue dotted line = population mean") +
theme_minimal(base_size = 10) +
theme(axis.text.y = element_text(size = 8), plot.title = element_text(face = "bold"))
} else {
panel_G <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
# Panel H: Life years by region
if (nrow(region_results) >= 3) {
panel_H <- ggplot(region_results, aes(x = life_years_gained, y = reorder(subgroup, dalys_averted))) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_vline(xintercept = overall_life, linetype = "dotted", color = "blue", size = 0.8) +
geom_point(size = 3, color = "darkorange") +
geom_errorbarh(aes(xmin = life_ci_lower, xmax = life_ci_upper), height = 0.2, size = 0.8, color = "darkorange") +
annotate("text", x = overall_life + 0.01, y = nrow(region_results) + 0.5, label = sprintf("Mean = %.4f", overall_life), size = 3, color = "blue") +
labs(x = "\nLife Years Gained (undiscounted)", y = "Region\n",
title = "H) Life Years Gained by Region (Undiscounted)",
subtitle = "Blue dotted line = population mean") +
theme_minimal(base_size = 10) +
theme(axis.text.y = element_text(size = 8), plot.title = element_text(face = "bold"))
} else {
panel_H <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
grid.arrange(panel_E, panel_F, panel_G, panel_H, ncol = 2, nrow = 2,
top = grid::textGrob("Figure 2: Sex and Regional Equity Analysis - DALYs & Life Years",
gp = grid::gpar(fontsize = 14, fontface = "bold")))
14. Figure 3: Cost Burden Analysis (4 Panels)
overall_costs <- mean(df_compare$delta_costs, na.rm = TRUE)
# Panel I: Cost burden by wealth
panel_I <- ggplot(wealth_results, aes(x = delta_costs, y = subgroup)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_vline(xintercept = overall_costs, linetype = "dotted", color = "blue", size = 0.8) +
geom_point(size = 4, color = "darkred") +
geom_errorbarh(aes(xmin = costs_ci_lower, xmax = costs_ci_upper), height = 0.2, size = 0.8, color = "darkred") +
annotate("text", x = overall_costs + 10, y = 5.2, label = sprintf("Mean = $%.0f", overall_costs), size = 3, color = "blue") +
labs(x = "\nIncremental Cost ($, discounted)", y = "Wealth Quintile\n",
title = "I) Cost Burden by Wealth Quintile",
subtitle = "Positive = higher cost for Empower | Blue dotted line = population mean") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
# Panel J: Costs concentration curve with 95% CI band
if (!is.null(costs_curve)) {
panel_J <- ggplot(costs_curve$curve, aes(x = Population, y = Cumulative)) +
geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), alpha = 0.2, fill = "darkred") +
geom_line(size = 1.2, color = "darkred") +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50", size = 0.8) +
annotate("text", x = 0.7, y = 0.25,
label = sprintf("C-index = %.4f\n95%% CI [%.4f, %.4f]",
costs_curve$C_index, costs_curve$C_ci_lower, costs_curve$C_ci_upper),
size = 3.5, fontface = "bold", color = "darkred") +
annotate("text", x = 0.7, y = 0.12,
label = ifelse(costs_curve$C_index > 0,
"Interpretation: PROGRESSIVE\n(Curve below diagonal = richer pay more)\n→ EQUITABLE (fair financing)",
"Interpretation: REGRESSIVE\n(Curve above diagonal = poorer pay more)\n→ INEQUITABLE (unfair burden)"),
size = 2.8, color = "darkred", lineheight = 1.2) +
labs(x = "\nCumulative population (poorest to richest)",
y = "Cumulative incremental costs (discounted)\n",
title = "J) Cost Burden Concentration Curve",
subtitle = "Shaded area = 95% CI | Curve below diagonal (C>0) = PROGRESSIVE") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold"))
} else {
panel_J <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
# Panel K: Cost burden by sex
panel_K <- ggplot(sex_results, aes(x = subgroup, y = delta_costs, fill = subgroup)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.6) +
geom_errorbar(aes(ymin = costs_ci_lower, ymax = costs_ci_upper), width = 0.2, size = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_hline(yintercept = overall_costs, linetype = "dotted", color = "blue", size = 0.8) +
annotate("text", x = 2.5, y = overall_costs + 10, label = sprintf("Mean = $%.0f", overall_costs), size = 3, color = "blue") +
labs(x = "\nSex", y = "Incremental Cost ($, discounted)\n",
title = "K) Cost Burden by Sex",
subtitle = "Blue dotted line = population mean") +
scale_fill_manual(values = c("Female" = "#FF6B6B", "Male" = "#4ECDC4")) +
theme_minimal(base_size = 11) +
theme(legend.position = "none", plot.title = element_text(face = "bold"))
# Panel L: Cost burden by region
if (nrow(region_results) >= 3) {
region_costs <- region_results %>%
arrange(desc(delta_costs)) %>%
mutate(cost_status = ifelse(delta_costs > 0, "Higher cost", "Lower cost"))
panel_L <- ggplot(region_costs, aes(x = delta_costs, y = reorder(subgroup, delta_costs), fill = cost_status)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.7) +
geom_errorbarh(aes(xmin = costs_ci_lower, xmax = costs_ci_upper), height = 0.2, size = 0.8) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 0.8) +
geom_vline(xintercept = overall_costs, linetype = "dotted", color = "blue", size = 0.8) +
annotate("text", x = overall_costs + 10, y = nrow(region_costs) + 0.5,
label = sprintf("Mean = $%.0f", overall_costs), size = 3, color = "blue") +
labs(x = "\nIncremental Cost ($, discounted)", y = "Region\n",
title = "L) Cost Burden by Region",
subtitle = "Blue dotted line = population mean | Red line = zero") +
scale_fill_manual(values = c("Higher cost" = "#D32F2F", "Lower cost" = "#2E7D32"),
name = "Cost Status") +
theme_minimal(base_size = 10) +
theme(axis.text.y = element_text(size = 8),
plot.title = element_text(face = "bold"),
legend.position = "bottom")
} else {
panel_L <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
grid.arrange(panel_I, panel_J, panel_K, panel_L, ncol = 2, nrow = 2,
top = grid::textGrob("Figure 3: Cost Burden Analysis - Wealth, Sex, and Regional Equity",
gp = grid::gpar(fontsize = 14, fontface = "bold")))
Figure 3: Cost Burden Analysis - Wealth, Sex, and Regional Equity
- Figure 4: ICERs by Equity Metric (4 Panels)
# Prepare ICER data
icer_wealth <- wealth_results %>%
filter(!is.na(icer_daly) & is.finite(icer_daly)) %>%
mutate(category = "Wealth Quintile",
label = as.character(subgroup),
ce_status = ifelse(icer_daly < WTP_THRESHOLD, "Cost-effective", "Not cost-effective"))
icer_sex <- sex_results %>%
filter(!is.na(icer_daly) & is.finite(icer_daly)) %>%
mutate(category = "Sex",
label = as.character(subgroup),
ce_status = ifelse(icer_daly < WTP_THRESHOLD, "Cost-effective", "Not cost-effective"))
icer_region <- region_results %>%
filter(!is.na(icer_daly) & is.finite(icer_daly)) %>%
mutate(category = "Region",
label = as.character(subgroup),
ce_status = ifelse(icer_daly < WTP_THRESHOLD, "Cost-effective", "Not cost-effective")) %>%
head(8)
# Panel M: ICERs by wealth
panel_M <- ggplot(icer_wealth, aes(x = icer_daly, y = reorder(label, -as.numeric(subgroup)), fill = ce_status)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.7) +
geom_vline(xintercept = icer_daly, linetype = "dotted", color = "gray50", size = 0.8) +
geom_text(aes(label = sprintf("$%.0f", icer_daly)), hjust = -0.2, size = 3) +
scale_fill_manual(values = c("Cost-effective" = "green", "Not cost-effective" = "red")) +
annotate("text", x = icer_daly + 50, y = 0.5, label = sprintf("Overall ICER = $%.0f", icer_daly), size = 3, color = "gray50") +
labs(x = "ICER ($ per DALY averted)", y = "",
title = "M) ICERs by Wealth Quintile") +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom", plot.title = element_text(face = "bold"))
# Panel N: ICERs by sex
panel_N <- ggplot(icer_sex, aes(x = icer_daly, y = label, fill = ce_status)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.6) +
geom_vline(xintercept = icer_daly, linetype = "dotted", color = "gray50", size = 0.8) +
geom_text(aes(label = sprintf("$%.0f", icer_daly)), hjust = -0.2, size = 3) +
scale_fill_manual(values = c("Cost-effective" = "green", "Not cost-effective" = "red")) +
annotate("text", x = icer_daly + 50, y = 2.2, label = sprintf("Overall ICER = $%.0f", icer_daly), size = 3, color = "gray50") +
labs(x = "ICER ($ per DALY averted)", y = "",
title = "N) ICERs by Sex") +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom", plot.title = element_text(face = "bold"))
# Panel O: ICERs by region
if (nrow(icer_region) >= 2) {
panel_O <- ggplot(icer_region, aes(x = icer_daly, y = reorder(label, icer_daly), fill = ce_status)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.6) +
geom_vline(xintercept = icer_daly, linetype = "dotted", color = "gray50", size = 0.8) +
geom_text(aes(label = sprintf("$%.0f", icer_daly)), hjust = -0.2, size = 2.5) +
scale_fill_manual(values = c("Cost-effective" = "green", "Not cost-effective" = "red")) +
annotate("text", x = icer_daly + 50, y = 8.5, label = sprintf("Overall ICER = $%.0f", icer_daly), size = 2.5, color = "gray50") +
labs(x = "ICER ($ per DALY averted)", y = "",
title = "O) ICERs by Region") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom", plot.title = element_text(face = "bold"))
} else {
panel_O <- ggplot() + annotate("text", x=0.5, y=0.5, label="Insufficient data") + theme_void()
}
# Panel P: Overall ICER comparison
icer_summary <- data.frame(
Metric = c("Overall", "Poorest (Q1)", "Richest (Q5)", "Male", "Female"),
ICER = c(icer_daly,
wealth_results$icer_daly[wealth_results$subgroup == "Q1_Poorest"],
wealth_results$icer_daly[wealth_results$subgroup == "Q5_Richest"],
sex_results$icer_daly[sex_results$subgroup == "Male"],
sex_results$icer_daly[sex_results$subgroup == "Female"])
)
icer_summary <- icer_summary %>%
mutate(ce_status = ifelse(ICER < WTP_THRESHOLD, "Cost-effective", "Not cost-effective"))
panel_P <- ggplot(icer_summary, aes(x = reorder(Metric, ICER), y = ICER, fill = ce_status)) +
geom_bar(stat = "identity", alpha = 0.8, width = 0.7) +
geom_text(aes(label = sprintf("$%.0f", ICER)), vjust = -0.5, size = 3) +
geom_hline(yintercept = WTP_THRESHOLD, linetype = "dashed", color = "blue", size = 0.8) +
scale_fill_manual(values = c("Cost-effective" = "green", "Not cost-effective" = "red")) +
annotate("text", x = 5.5, y = WTP_THRESHOLD + 50, label = sprintf("Kenya GDP = $%.0f", WTP_THRESHOLD), size = 2.5, color = "blue") +
labs(x = "", y = "ICER ($ per DALY averted)\n",
title = "P) ICER Comparison Across Subgroups") +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom", plot.title = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1))
grid.arrange(panel_M, panel_N, panel_O, panel_P, ncol = 2, nrow = 2,
top = grid::textGrob("Figure 4: Cost-Effectiveness ICERs by Equity Metric",
gp = grid::gpar(fontsize = 14, fontface = "bold")))
15. Export Results
# Save data
write.csv(wealth_results, "kenya_wealth_results.csv", row.names = FALSE)
write.csv(sex_results, "kenya_sex_results.csv", row.names = FALSE)
write.csv(region_results, "kenya_region_results.csv", row.names = FALSE)
# Save concentration curve data with CIs
if (!is.null(health_curve)) {
write.csv(health_curve$curve, "kenya_concentration_curve_dalys.csv", row.names = FALSE)
write.csv(data.frame(C_index = health_curve$C_index,
CI_lower = health_curve$C_ci_lower,
CI_upper = health_curve$C_ci_upper),
"kenya_concentration_curve_dalys_summary.csv", row.names = FALSE)
}
if (!is.null(life_curve)) {
write.csv(life_curve$curve, "kenya_concentration_curve_life_years.csv", row.names = FALSE)
write.csv(data.frame(C_index = life_curve$C_index,
CI_lower = life_curve$C_ci_lower,
CI_upper = life_curve$C_ci_upper),
"kenya_concentration_curve_life_summary.csv", row.names = FALSE)
}
if (!is.null(costs_curve)) {
write.csv(costs_curve$curve, "kenya_concentration_curve_costs.csv", row.names = FALSE)
write.csv(data.frame(C_index = costs_curve$C_index,
CI_lower = costs_curve$C_ci_lower,
CI_upper = costs_curve$C_ci_upper),
"kenya_concentration_curve_costs_summary.csv", row.names = FALSE)
}
# Policy summary with CIs
policy_summary <- data.frame(
Metric = c(
"Discount rate (costs & DALYs)",
"Life years discounting",
"Overall ICER (per DALY)",
"Kenya GDP per capita",
"Cost-effective?",
"DALYs C-index (95% CI)",
"DALYs distribution",
"Life Years C-index (95% CI)",
"Life Years distribution",
"Costs C-index (95% CI)",
"Costs distribution"
),
Value = c(
sprintf("%.1f%%", DISCOUNT_RATE * 100),
"Not discounted",
sprintf("$%.0f per DALY", icer_daly),
sprintf("$%.0f", WTP_THRESHOLD),
ifelse(icer_daly < WTP_THRESHOLD, "YES", "NO"),
ifelse(!is.null(health_curve),
sprintf("%.4f (%.4f-%.4f)", health_curve$C_index, health_curve$C_ci_lower, health_curve$C_ci_upper),
"NA"),
ifelse(!is.null(health_curve) && health_curve$C_index > 0, "PRO-RICH (INEQUITABLE)", "PRO-POOR (EQUITABLE)"),
ifelse(!is.null(life_curve),
sprintf("%.4f (%.4f-%.4f)", life_curve$C_index, life_curve$C_ci_lower, life_curve$C_ci_upper),
"NA"),
ifelse(!is.null(life_curve) && life_curve$C_index > 0, "PRO-RICH (INEQUITABLE)", "PRO-POOR (EQUITABLE)"),
ifelse(!is.null(costs_curve),
sprintf("%.4f (%.4f-%.4f)", costs_curve$C_index, costs_curve$C_ci_lower, costs_curve$C_ci_upper),
"NA"),
ifelse(!is.null(costs_curve) && costs_curve$C_index > 0, "PROGRESSIVE (EQUITABLE)", "REGRESSIVE (INEQUITABLE)")
)
)
write.csv(policy_summary, "kenya_equity_policy_summary.csv", row.names = FALSE)
saveRDS(df_compare, "kenya_equity_analysis_data.rds")
cat("\n✓ All results exported successfully\n")##
## ✓ All results exported successfully- Summary of Findings
## ========================================## EQUITY ANALYSIS COMPLETE## ========================================## Overall ICER: $3990 per DALY##
## --- CONCENTRATION INDICES WITH 95% CI ---if (!is.null(health_curve)) {
cat(sprintf("DALYs: %.4f (95%% CI: %.4f - %.4f) - %s\n",
health_curve$C_index, health_curve$C_ci_lower, health_curve$C_ci_upper,
ifelse(health_curve$C_index > 0, "PRO-RICH", "PRO-POOR")))
}## DALYs: 0.0451 (95% CI: 0.0123 - 0.0780) - PRO-RICHif (!is.null(life_curve)) {
cat(sprintf("Life Years: %.4f (95%% CI: %.4f - %.4f) - %s\n",
life_curve$C_index, life_curve$C_ci_lower, life_curve$C_ci_upper,
ifelse(life_curve$C_index > 0, "PRO-RICH", "PRO-POOR")))
}## Life Years: 0.0622 (95% CI: 0.0278 - 0.0966) - PRO-RICHif (!is.null(costs_curve)) {
cat(sprintf("Costs: %.4f (95%% CI: %.4f - %.4f) - %s\n",
costs_curve$C_index, costs_curve$C_ci_lower, costs_curve$C_ci_upper,
ifelse(costs_curve$C_index > 0, "PROGRESSIVE", "REGRESSIVE")))
}## Costs: 0.0646 (95% CI: 0.0437 - 0.0855) - PROGRESSIVE##
## ========================================