全球卫生支出 2000–2023：公平、韧性、未来

# =============================================================================
# 程序/07_plot_static.R
# -----------------------------------------------------------------------------
# 静态图工厂函数：输入 tidy tibble，输出 ggplot 对象。
# 所有函数签名：plot_*(data, ..., title, subtitle, caption)
# =============================================================================

if (!exists("proj_root", mode = "function")) {
  source(file.path("程序", "00_utils.R"))
  source(file.path("程序", "06_plot_theme.R"))
}

# =============================================================================
# A. 趋势 / 时序
# =============================================================================

#' 全球三源（政府/私人/外援）占比时序堆叠面积
#'
#' @param master 宽表 (含 gghed_che, pvtd_che, ext_che, year, pop)
plot_source_area <- function(master,
                             title = "全球卫生支出来源结构演化 Global Health Financing Sources",
                             subtitle = NULL) {
  ensure_pkgs(c("dplyr", "tidyr", "ggplot2", "scales"))
  has_pop <- "pop" %in% names(master) && any(is.finite(master$pop))
  safe_wmean <- function(x, w) {
    ok <- is.finite(x) & is.finite(w) & w > 0
    if (!any(ok)) return(mean(x, na.rm = TRUE))
    stats::weighted.mean(x[ok], w[ok], na.rm = TRUE)
  }
  if (is.null(subtitle)) {
    subtitle <- if (has_pop) {
      "2000–2023 · 按人口加权均值，占总开支 % of CHE"
    } else {
      "2000–2023 · 未加权均值，占总开支 % of CHE"
    }
  }
  d <- master |>
    dplyr::group_by(.data$year) |>
    dplyr::summarise(
      gghed = if (has_pop) {
        safe_wmean(.data$gghed_che, .data$pop)
      } else mean(.data$gghed_che, na.rm = TRUE),
      pvtd  = if (has_pop) {
        safe_wmean(.data$pvtd_che, .data$pop)
      } else mean(.data$pvtd_che, na.rm = TRUE),
      ext   = if (has_pop) {
        safe_wmean(.data$ext_che, .data$pop)
      } else mean(.data$ext_che, na.rm = TRUE),
      .groups = "drop"
    ) |>
    tidyr::pivot_longer(c("gghed", "pvtd", "ext"),
                        names_to = "source", values_to = "pct") |>
    dplyr::mutate(source = factor(.data$source, levels = c("gghed", "pvtd", "ext")))
  ggplot2::ggplot(d, ggplot2::aes(.data$year, .data$pct, fill = .data$source)) +
    ggplot2::geom_area(alpha = 0.9, position = "stack") +
    scale_fill_ghs_source(
      labels = c(gghed = "政府 GGHE-D", pvtd = "私人 PVT-D", ext = "外援 EXT")
    ) +
    ggplot2::scale_y_continuous(labels = scales::label_percent(scale = 1)) +
    labs_ghs(

library(dplyr)
rank_oops <- master |>
  dplyr::filter(year == max(year, na.rm = TRUE)) |>
  dplyr::select(country_name, continent, hf3_che, che_pc_usd2023) |>
  dplyr::arrange(dplyr::desc(hf3_che))
utils::head(rank_oops, 10)

gini_weighted <- function(x, w = NULL) {
  if (is.null(w)) w <- rep(1, length(x))
  ok <- is.finite(x) & is.finite(w) & x >= 0 & w > 0
  x <- x[ok]; w <- w[ok]
  if (length(x) < 2) return(NA_real_)
  ord <- order(x)
  x <- x[ord]; w <- w[ord]
  cum_w <- cumsum(w)
  cum_xw <- cumsum(x * w)
  total_w <- sum(w); total_xw <- sum(x * w)
  if (total_xw == 0) return(0)
  g <- 1 - 2 * sum((cum_xw - x * w / 2) * w) / (total_w * total_xw)
  g
}

#' Theil-T 指数（人口加权）
theil_t <- function(x, w = NULL) {
  if (is.null(w)) w <- rep(1, length(x))
  ok <- is.finite(x) & is.finite(w) & x > 0 & w > 0
  x <- x[ok]; w <- w[ok]
  if (length(x) < 2) return(NA_real_)
  mean_x <- stats::weighted.mean(x, w)
  if (mean_x == 0) return(0)
  sum(w * x / sum(w * x) * log(x / mean_x))
}

#' Atkinson 指数 (ε)
atkinson <- function(x, w = NULL, eps = 0.5) {
  if (is.null(w)) w <- rep(1, length(x))
  ok <- is.finite(x) & is.finite(w) & x > 0 & w > 0
  x <- x[ok]; w <- w[ok]
  if (length(x) < 2) return(NA_real_)
  mean_x <- stats::weighted.mean(x, w)
  if (eps == 1) {
    ede <- exp(stats::weighted.mean(log(x), w))
  } else {
    ede <- (stats::weighted.mean(x^(1 - eps), w))^(1 / (1 - eps))
  }
  1 - ede / mean_x
}

# ---- 2. 时间序列的不平等 ---------------------------------------------------
#' 按 year 计算全球人均 CHE 的不平等指数面板
inequality_by_year <- function(master, value_col = "che_pc_usd2023",
                               weight_col = "pop") {
  ensure_pkgs(c("dplyr", "tidyr", "tibble"))
  val <- rlang::sym(value_col)
  wgt <- rlang::sym(weight_col)
  master |>
    dplyr::filter(
      is.finite(!!val), !!val > 0,
      is.finite(!!wgt) | TRUE
    ) |>
    dplyr::group_by(.data$year) |>
    dplyr::summarise(
      n            = dplyr::n(),
      gini_eq      = gini_weighted(!!val, NULL),
      gini_pop     = if (weight_col %in% names(master)) {
        gini_weighted(!!val, !!wgt)
      } else NA_real_,
      theil_pop    = if (weight_col %in% names(master)) {
        theil_t(!!val, !!wgt)
      } else NA_real_,
      atk05        = atkinson(!!val, NULL, eps = 0.5),
      atk1         = atkinson(!!val, NULL, eps = 1),
      mean_val     = mean(!!val, na.rm = TRUE),
      median_val   = median(!!val, na.rm = TRUE),
      .groups      = "drop"
    )
}

covid_shock <- function(master, base_year = 2019,
                        shock_years = 2020:2022) {
  ensure_pkgs("dplyr")
  master |>
    dplyr::filter(.data$year %in% c(base_year, shock_years)) |>
    dplyr::mutate(
      phase = dplyr::if_else(.data$year == base_year, "base", "shock")
    ) |>
    dplyr::group_by(.data$iso3_code, .data$country_name) |>
    dplyr::summarise(
      base_che    = mean(.data$che_usd2023[.data$phase == "base"], na.rm = TRUE),
      shock_che   = mean(.data$che_usd2023[.data$phase == "shock"], na.rm = TRUE),
      base_oops   = mean(.data$hf3_che[.data$phase == "base"], na.rm = TRUE),
      shock_oops  = mean(.data$hf3_che[.data$phase == "shock"], na.rm = TRUE),
      base_gghed  = mean(.data$gghed_che[.data$phase == "base"], na.rm = TRUE),
      shock_gghed = mean(.data$gghed_che[.data$phase == "shock"], na.rm = TRUE),
      .groups     = "drop"
    ) |>
    dplyr::mutate(
      che_delta_pct   = (.data$shock_che  - .data$base_che)  / .data$base_che  * 100,
      oops_delta_pp   =  .data$shock_oops - .data$base_oops,
      gghed_delta_pp  =  .data$shock_gghed - .data$base_gghed
    )
}

#' log(人均 CHE) 增速 ~ log(起始水平) + 固定效应
#' 输入：master 宽表（含 iso3_code, year, che_pc_usd2023, continent）
fit_beta_convergence <- function(master, start_year = 2000, end_year = 2023) {
  ensure_pkgs(c("dplyr"))
  if (!requireNamespace("fixest", quietly = TRUE)) {
    warning("fixest not installed; beta-convergence skipped")
    return(NULL)
  }
  panel <- master |>
    dplyr::filter(.data$year %in% c(start_year, end_year),
                  is.finite(.data$che_pc_usd2023),
                  .data$che_pc_usd2023 > 0) |>
    dplyr::select("iso3_code", "continent", "year", "che_pc_usd2023") |>
    tidyr::pivot_wider(names_from = "year", values_from = "che_pc_usd2023",
                       names_prefix = "y_") |>
    tidyr::drop_na() |>
    dplyr::mutate(
      log_start  = log(!!rlang::sym(paste0("y_", start_year))),
      log_end    = log(!!rlang::sym(paste0("y_", end_year))),
      growth     = (.data$log_end - .data$log_start) / (end_year - start_year)
    )
  if (nrow(panel) < 10) return(NULL)
  mod <- tryCatch(
    fixest::feols(growth ~ log_start | continent, data = panel),
    error = function(e) lm(growth ~ log_start + continent, data = panel)
  )
  list(
    panel = panel,
    model = mod,
    beta  = stats::coef(mod)[["log_start"]],
    half_life = if ("log_start" %in% names(stats::coef(mod))) {
      -log(2) / stats::coef(mod)[["log_start"]]
    } else NA_real_
  )
}

#' log(人均 CHE) ~ log(GDP/cap) + income_group | country + year
fit_panel_fe <- function(master) {
  ensure_pkgs("dplyr")
  if (!requireNamespace("fixest", quietly = TRUE)) {
    warning("fixest not installed; returning NULL")
    return(NULL)
  }
  dat <- master |>
    dplyr::filter(is.finite(.data$che_pc_usd2023), .data$che_pc_usd2023 > 0,
                  is.finite(.data$gdp_pc_usd),      .data$gdp_pc_usd > 0) |>
    dplyr::mutate(
      log_che_pc = log(.data$che_pc_usd2023),
      log_gdp_pc = log(.data$gdp_pc_usd)
    )
  if (nrow(dat) < 100) return(NULL)
  tryCatch(
    fixest::feols(log_che_pc ~ log_gdp_pc | iso3_code + year, data = dat,
                  cluster = ~iso3_code),
    error = function(e) NULL
  )
}

fit_pca <- function(df, vars, scale = TRUE) {
  ensure_pkgs(c("dplyr"))
  keep <- df |>
    dplyr::select(dplyr::all_of(c("iso3_code", "country_name", vars))) |>
    tidyr::drop_na()
  if (nrow(keep) < 5) {
    warning("PCA: fewer than 5 complete observations; returning NULL")
    return(NULL)
  }
  X <- as.matrix(keep[, vars, drop = FALSE])
  rownames(X) <- keep$iso3_code
  pca <- stats::prcomp(X, center = TRUE, scale. = scale)
  list(
    pca    = pca,
    scores = as.data.frame(pca$x) |>
      tibble::rownames_to_column("iso3_code") |>
      dplyr::left_join(keep[, c("iso3_code", "country_name")], by = "iso3_code"),
    loadings = as.data.frame(pca$rotation) |> tibble::rownames_to_column("variable"),
    var_explained = summary(pca)$importance["Proportion of Variance", ]
  )
}

# ---- 2. k-means 聚类 --------------------------------------------------------
#' 在 PCA 前两主成分上做 k-means
fit_cluster <- function(pca_obj, k = 4, seed = 1L) {
  if (is.null(pca_obj)) return(NULL)
  set.seed(seed)
  X <- pca_obj$scores[, grep("^PC", names(pca_obj$scores))]
  if (ncol(X) < 2) {
    warning("cluster: need at least 2 PCs")
    return(NULL)
  }
  km <- stats::kmeans(X[, 1:min(ncol(X), 5)], centers = k, nstart = 25)
  scores <- pca_obj$scores
  scores$cluster <- factor(km$cluster)
  list(km = km, scores = scores, k = k)
}

#' 用 fable 或 forecast 对一国一指标做 h 期预测
fit_forecast <- function(ts_vec, years, h = 5, method = c("auto", "prophet")) {
  method <- match.arg(method)
  ts_vec <- as.numeric(ts_vec)
  ok <- is.finite(ts_vec)
  if (sum(ok) < 5) return(NULL)
  ts_vec <- ts_vec[ok]; years <- years[ok]
  if (method == "auto" && requireNamespace("forecast", quietly = TRUE)) {
    ts_obj <- stats::ts(ts_vec, start = min(years), frequency = 1)
    fit <- tryCatch(forecast::auto.arima(ts_obj), error = function(e) NULL)
    if (is.null(fit)) return(NULL)
    fc <- forecast::forecast(fit, h = h)
    return(tibble::tibble(
      year  = (max(years) + 1):(max(years) + h),
      point = as.numeric(fc$mean),
      lo_80 = as.numeric(fc$lower[, 1]),
      hi_80 = as.numeric(fc$upper[, 1]),
      lo_95 = as.numeric(fc$lower[, 2]),
      hi_95 = as.numeric(fc$upper[, 2]),
      method = "ARIMA"
    ))
  }
  if (method == "prophet" && requireNamespace("prophet", quietly = TRUE)) {
    dfp <- data.frame(
      ds = as.Date(paste0(years, "-01-01")),
      y  = ts_vec
    )
    m <- tryCatch(
      prophet::prophet(dfp, yearly.seasonality = FALSE, weekly.seasonality = FALSE,
                       daily.seasonality = FALSE, verbose = FALSE),
      error = function(e) NULL
    )
    if (is.null(m)) return(NULL)
    future <- prophet::make_future_dataframe(m, periods = h, freq = "year")
    pred <- stats::predict(m, future)
    tail_n <- utils::tail(pred, h)
    return(tibble::tibble(
      year  = as.integer(format(tail_n$ds, "%Y")),
      point = tail_n$yhat,
      lo_80 = tail_n$yhat_lower,
      hi_80 = tail_n$yhat_upper,
      lo_95 = tail_n$yhat_lower,
      hi_95 = tail_n$yhat_upper,
      method = "Prophet"
    ))
  }
  NULL
}

library(dplyr)
aid_rank <- master |>
  dplyr::filter(year == max(year, na.rm = TRUE)) |>
  dplyr::select(country_name, continent, ext_che, gghed_che, hf3_che, che_pc_usd2023) |>
  dplyr::arrange(dplyr::desc(ext_che))
head(aid_rank, 10)

library(dplyr)
fit <- master |>
  dplyr::filter(year == max(year, na.rm = TRUE),
                is.finite(che_pc_usd2023), is.finite(life_exp)) |>
  lm(life_exp ~ log(che_pc_usd2023), data = _)
summary(fit)

library(dplyr)
rk <- master |>
  dplyr::filter(year %in% c(2000, 2010, 2023),
                is.finite(che_pc_usd2023)) |>
  dplyr::group_by(year) |>
  dplyr::mutate(rank = rank(-che_pc_usd2023)) |>
  dplyr::ungroup() |>
  dplyr::filter(rank <= 25)

fit <- lm(life_exp ~ log(che_pc_usd2023),
          data = subset(master, year == max(year, na.rm = TRUE) &
                        is.finite(che_pc_usd2023) & is.finite(life_exp)))
summary(fit)