# =============================================================================
# Machakos County dog rabies vaccination: Monte Carlo cost model

#
# Reproduces all results in the manuscript: base case, donated-vaccine
# scenario, optimised scenarios (90, 60 and 360 days), PRCC sensitivity
# analysis, and the cost component breakdown. Writes every result table to
# ../output/ as CSV.
#
# Reporting conventions
#   1. Headline point estimates are medians (cost distributions are right-skewed).
#   2. All KES to USD conversions use a fixed rate (usd_rate <- 129).
#   3. Uncertainty intervals are the 5th and 95th percentiles (90% UI).
#
# Requirements: R >= 4.1; packages mc2d, sensitivity, dplyr, ggplot2.
# Run time: about 2 minutes on a standard laptop.
#
# Usage: set the working directory to the code/ folder, then source this file.
#   setwd(".../machakos_rabies_cost_model/code"); source("machakos_cost_model.R")
# =============================================================================

rm(list = ls())

library(mc2d)
library(ggplot2)
library(dplyr)
library(sensitivity)

out_dir <- "../output"
dir.create(out_dir, showWarnings = FALSE)

n_sim    <- 50000
usd_rate <- 129
set.seed(123)

rtriang <- function(n, min, mode, max) mc2d::rtriang(n, min = min, mode = mode, max = max)

# =============================================================================
# 1. Base-case inputs (Table 2)
# =============================================================================

dog_population    <- rtriang(n_sim, 408702, 421079, 424950)   # P
coverage          <- rtriang(n_sim, 0.5, 0.7, 0.8)            # C
throughput        <- rtriang(n_sim, 20, 40, 100)              # VT, dogs/team/day
vaccine_wastage   <- rtriang(n_sim, 0.02, 0.16, 0.20)         # W, proportion
campaign_duration <- rtriang(n_sim, 60, 90, 360)              # CD, days

# This draw is NOT used in any result. It is retained only so that the random
# number stream is identical to the run reported in the manuscript; deleting
# it would change every subsequent draw and therefore the reported figures.
exchange_rate_unused <- runif(n_sim, 100, 150)

D     <- dog_population * coverage                       # dogs to vaccinate
teams <- ceiling((D / campaign_duration) / throughput)   # T

coordination_staff_salaries <- rep(2244000, n_sim)       # a, 100% FTE
admin_m_e_cost              <- rep(230000, n_sim)        # b
fixed_cost_kes <- coordination_staff_salaries + admin_m_e_cost

vaccine_cost             <- rtriang(n_sim, 50, 60, 70)   # c
certificate_cost         <- rtriang(n_sim, 20, 23, 30)   # d
core_delivery_cost       <- rtriang(n_sim, 5, 6, 10)     # e
non_core_medicament_cost <- rtriang(n_sim, 10, 14, 15)   # f

vc_total_kes <- (vaccine_cost * (1 + vaccine_wastage) + certificate_cost +
                 core_delivery_cost * (1 + vaccine_wastage) + non_core_medicament_cost) * D

personnel_cost      <- rtriang(n_sim, 3000, 4000, 5000)     # h, per team per day
transport_cost      <- rtriang(n_sim, 3000, 5000, 10000)    # i, per vehicle per day
equipment_cost_team <- rtriang(n_sim, 5000, 7000, 10000)    # g, per team
mobilization_cost   <- rtriang(n_sim, 20000, 35000, 35000)  # j, per team

personnel_total    <- personnel_cost * campaign_duration
transport_total    <- transport_cost * campaign_duration * 0.5   # one vehicle per two teams
equipment_total    <- equipment_cost_team
mobilization_total <- mobilization_cost

svc_total_kes <- (personnel_total + transport_total + equipment_total + mobilization_total) * teams

base_cost_kes    <- fixed_cost_kes + vc_total_kes + svc_total_kes
base_cost_usd    <- base_cost_kes / usd_rate
cost_per_dog_kes <- base_cost_kes / D
cost_per_dog_usd <- cost_per_dog_kes / usd_rate

# =============================================================================
# 2. Donated-vaccine scenario
# =============================================================================

vc_donated_kes <- (certificate_cost + core_delivery_cost * (1 + vaccine_wastage) +
                   non_core_medicament_cost) * D
donated_cost_kes <- fixed_cost_kes + vc_donated_kes + svc_total_kes
donated_cost_usd <- donated_cost_kes / usd_rate

# =============================================================================
# 3. Optimised scenarios (Table 3): 90, 60 and 360 days
#    Throughput ~ triangular(40, 80, 100); wastage fixed at 0.06;
#    coordination salaries at 50% FTE; duration fixed.
# =============================================================================

set.seed(123)
throughput_opt      <- rtriang(n_sim, 40, 80, 100)
vaccine_wastage_opt <- rep(0.06, n_sim)
fixed_cost_kes_opt  <- coordination_staff_salaries * 0.5 + admin_m_e_cost
D_opt               <- D

vc_total_kes_opt <- (vaccine_cost * (1 + vaccine_wastage_opt) + certificate_cost +
                     core_delivery_cost * (1 + vaccine_wastage_opt) + non_core_medicament_cost) * D_opt

run_optimised <- function(duration) {
  cd     <- rep(duration, n_sim)
  t_opt  <- ceiling((D_opt / cd) / throughput_opt)
  svc    <- (personnel_cost * cd + transport_cost * cd * 0.5 +
             equipment_cost_team + mobilization_cost) * t_opt
  total  <- fixed_cost_kes_opt + vc_total_kes_opt + svc
  list(teams = t_opt, svc = svc, total_kes = total, total_usd = total / usd_rate,
       cpd_kes = total / D_opt, cpd_usd = total / D_opt / usd_rate, cd = cd)
}

opt90  <- run_optimised(90)
opt60  <- run_optimised(60)
opt360 <- run_optimised(360)

# =============================================================================
# 4. Per-iteration outputs
# =============================================================================

write.csv(data.frame(
  iteration = seq_len(n_sim), dog_population, coverage, dogs_to_vaccinate = D,
  throughput, campaign_duration, vaccine_wastage, teams,
  fixed_cost_kes, variable_cost_kes = vc_total_kes, semivariable_cost_kes = svc_total_kes,
  total_cost_kes = base_cost_kes, total_cost_usd = base_cost_usd,
  cost_per_dog_kes, cost_per_dog_usd),
  file.path(out_dir, "simulation_results_base.csv"), row.names = FALSE)

write.csv(data.frame(
  iteration = seq_len(n_sim), dogs_to_vaccinate = D, teams,
  total_cost_kes = donated_cost_kes, total_cost_usd = donated_cost_usd,
  cost_per_dog_kes = donated_cost_kes / D, cost_per_dog_usd = donated_cost_kes / D / usd_rate),
  file.path(out_dir, "simulation_results_donated.csv"), row.names = FALSE)

for (nm in c("opt90", "opt60", "opt360")) {
  s <- get(nm)
  write.csv(data.frame(
    iteration = seq_len(n_sim), dogs_to_vaccinate = D_opt, throughput = throughput_opt,
    campaign_duration = s$cd, teams = s$teams,
    fixed_cost_kes = fixed_cost_kes_opt, variable_cost_kes = vc_total_kes_opt,
    semivariable_cost_kes = s$svc, total_cost_kes = s$total_kes, total_cost_usd = s$total_usd,
    cost_per_dog_kes = s$cpd_kes, cost_per_dog_usd = s$cpd_usd),
    file.path(out_dir, paste0("simulation_results_", nm, ".csv")), row.names = FALSE)
}

# =============================================================================
# 5. Summary statistics, all scenarios
# =============================================================================

summarise_scenario <- function(name, total_kes, cpd_kes, teams) {
  q <- function(x, p) unname(quantile(x, p))
  data.frame(
    scenario         = name,
    median_total_kes = median(total_kes), mean_total_kes = mean(total_kes),
    p5_total_kes     = q(total_kes, .05), p95_total_kes  = q(total_kes, .95),
    median_total_usd = median(total_kes) / usd_rate, mean_total_usd = mean(total_kes) / usd_rate,
    p5_total_usd     = q(total_kes, .05) / usd_rate, p95_total_usd = q(total_kes, .95) / usd_rate,
    median_cpd_kes   = median(cpd_kes), mean_cpd_kes = mean(cpd_kes),
    p5_cpd_kes       = q(cpd_kes, .05), p95_cpd_kes = q(cpd_kes, .95),
    median_cpd_usd   = median(cpd_kes) / usd_rate, mean_cpd_usd = mean(cpd_kes) / usd_rate,
    p5_cpd_usd       = q(cpd_kes, .05) / usd_rate, p95_cpd_usd = q(cpd_kes, .95) / usd_rate,
    median_teams     = median(teams), p5_teams = q(teams, .05), p95_teams = q(teams, .95))
}

summary_table <- rbind(
  summarise_scenario("Base case",            base_cost_kes,    cost_per_dog_kes,      teams),
  summarise_scenario("Donated vaccines",     donated_cost_kes, donated_cost_kes / D,  teams),
  summarise_scenario("Optimised (90-day)",   opt90$total_kes,  opt90$cpd_kes,  opt90$teams),
  summarise_scenario("Optimised (60-day)",   opt60$total_kes,  opt60$cpd_kes,  opt60$teams),
  summarise_scenario("Optimised (360-day)",  opt360$total_kes, opt360$cpd_kes, opt360$teams))

base_med <- summary_table$median_total_usd[1]
summary_table$pct_change_vs_base_median <- round((base_med - summary_table$median_total_usd) / base_med * 100, 1)

write.csv(summary_table, file.path(out_dir, "summary_statistics.csv"), row.names = FALSE)
print(summary_table)

# =============================================================================
# 6. Cost component breakdown (median KES per component, Figure 7)
# =============================================================================

component_medians <- function(name, wast, Dn, pers, trans, equip, mob, tm, fixed) {
  data.frame(scenario = name,
    component = c("Vaccine", "Certificate", "Core delivery", "Non-core medicaments",
                  "Personnel", "Transport", "Equipment", "Mobilisation", "Fixed (coordination and admin)"),
    median_kes = c(
      median(vaccine_cost * (1 + wast) * Dn), median(certificate_cost * Dn),
      median(core_delivery_cost * (1 + wast) * Dn), median(non_core_medicament_cost * Dn),
      median(pers * tm), median(trans * tm), median(equip * tm), median(mob * tm), median(fixed)))
}

component_table <- rbind(
  component_medians("Base case", vaccine_wastage, D, personnel_total, transport_total,
                    equipment_total, mobilization_total, teams, fixed_cost_kes),
  component_medians("Optimised (90-day)", vaccine_wastage_opt, D_opt, personnel_cost * 90,
                    transport_cost * 90 * 0.5, equipment_cost_team, mobilization_cost, opt90$teams, fixed_cost_kes_opt),
  component_medians("Optimised (60-day)", vaccine_wastage_opt, D_opt, personnel_cost * 60,
                    transport_cost * 60 * 0.5, equipment_cost_team, mobilization_cost, opt60$teams, fixed_cost_kes_opt),
  component_medians("Optimised (360-day)", vaccine_wastage_opt, D_opt, personnel_cost * 360,
                    transport_cost * 360 * 0.5, equipment_cost_team, mobilization_cost, opt360$teams, fixed_cost_kes_opt))

write.csv(component_table, file.path(out_dir, "cost_components_by_scenario.csv"), row.names = FALSE)

# =============================================================================
# 7. PRCC sensitivity analysis (base case, outcome = cost per dog, USD)
# =============================================================================

X <- data.frame(vaccine_cost, certificate_cost, core_delivery_cost, non_core_medicament_cost,
                personnel_cost, transport_cost, equipment_cost_team, mobilization_cost,
                coverage, throughput, campaign_duration, dog_population, vaccine_wastage)

set.seed(123)
prcc_result <- pcc(X, cost_per_dog_usd, rank = TRUE, nboot = 200)

prcc_df <- data.frame(
  variable = rownames(prcc_result$PRCC),
  prcc     = prcc_result$PRCC[, "original"],
  ci_lower = prcc_result$PRCC[, "min. c.i."],
  ci_upper = prcc_result$PRCC[, "max. c.i."]) %>%
  arrange(desc(abs(prcc)))

write.csv(prcc_df, file.path(out_dir, "prcc_results.csv"), row.names = FALSE)
print(prcc_df)

# =============================================================================
# 8. Figures
# =============================================================================

fig4 <- ggplot(data.frame(x = base_cost_usd / 1e6), aes(x)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.05, fill = "#AEC7E8", colour = "black", alpha = 0.7) +
  geom_density(colour = "#1F77B4", linewidth = 1.2) +
  geom_vline(xintercept = median(base_cost_usd) / 1e6, linetype = "dashed") +
  geom_vline(xintercept = quantile(base_cost_usd, c(.05, .95)) / 1e6, linetype = "dotted", colour = "grey40") +
  labs(x = "Total campaign cost (million USD)", y = "Density") +
  theme_minimal(base_size = 14)
ggsave(file.path(out_dir, "Figure4_total_cost_distribution.png"), fig4, width = 8, height = 5, dpi = 300, bg = "white")

fig5 <- ggplot(prcc_df, aes(x = reorder(gsub("_", " ", variable), abs(prcc)), y = prcc)) +
  geom_col(fill = "#1F77B4") +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  coord_flip() +
  labs(x = "Input variable", y = "Partial rank correlation coefficient") +
  theme_minimal(base_size = 14)
ggsave(file.path(out_dir, "Figure5_prcc.png"), fig5, width = 8, height = 6, dpi = 300, bg = "white")

fig6 <- ggplot(summary_table[1:2, ], aes(x = scenario, y = median_total_usd / 1e6, fill = scenario)) +
  geom_col(width = 0.5, colour = "black") +
  geom_errorbar(aes(ymin = p5_total_usd / 1e6, ymax = p95_total_usd / 1e6), width = 0.15) +
  scale_fill_manual(values = c("#E69F00", "#56B4E9")) +
  labs(x = NULL, y = "Total campaign cost (million USD)") +
  theme_minimal(base_size = 14) + theme(legend.position = "none")
ggsave(file.path(out_dir, "Figure6_base_vs_donated.png"), fig6, width = 6, height = 5, dpi = 300, bg = "white")

component_table$scenario  <- factor(component_table$scenario, levels = unique(component_table$scenario))
component_table$component <- factor(component_table$component, levels = unique(component_table$component))
fig7 <- ggplot(component_table, aes(x = scenario, y = median_kes / 1e6, fill = component)) +
  geom_col(width = 0.6, colour = "white", linewidth = 0.3) +
  scale_fill_manual(values = c("#08306B", "#2171B5", "#4292C6", "#6BAED6", "#9ECAE1",
                               "#C6DBEF", "#969696", "#BDBDBD", "#D9D9D9")) +
  labs(x = NULL, y = "Median total campaign cost (million KES)", fill = "Cost component") +
  theme_minimal(base_size = 13)
ggsave(file.path(out_dir, "Figure7_cost_components.png"), fig7, width = 9, height = 6.5, dpi = 300, bg = "white")

# =============================================================================
# 9. Session information
# =============================================================================

writeLines(capture.output(sessionInfo()), file.path(out_dir, "sessionInfo.txt"))
cat("Done. Outputs written to", normalizePath(out_dir), "\n")