Shapley decomposition of HLE variation

shapley.R

load("Data/jdm_estimates.Rda")
library(tidyverse)
library(DemoTools)

calc_Lx <- function(lmx){
  n <- length(lmx)
  ax <- c(.5,rep(2.5,n-1))
  qx <- lt_id_ma_q(nMx = exp(lmx), 
                   nax = ax, 
                   AgeInt = rep(5,n))
  lx <- c(1,cumprod(1-qx))
  dx <- lx[1:n] * qx
  Lx <- 5 * lx[1:n] - (5 - ax) * dx
  Lx
}

# make this easier for me to remember colnames:
data_area <- data.out |> 
  select(area = clusterid, sex = sex, 
         year = year, age = age.start, 
         lmx = ETAhat.mort, prev = phat, 
         exposure = exposures)

# create Basque reference schedules for Lx, prev. These vary by
# year and sex.
ref_sullivan <-
  data_area |> 
  mutate(Dx_hat = exp(lmx) * exposure,
         cases_hat = prev * exposure) |> 
  group_by(sex, year, age) |> 
  summarize(Dx_hat = sum(Dx_hat),
            cases_hat = sum(cases_hat),
            exposure = sum(exposure), .groups = "drop") |> 
  group_by(sex, year) |> 
  mutate(mx_0 = Dx_hat / exposure,
         prev_0 = cases_hat / exposure,
         Lx_0 = calc_Lx(log(mx_0)))

# Join the basque reference to each area
dat_cf <- data_area |>
  group_by(area,sex,year) |> 
  mutate(Lx = calc_Lx(lmx)) |> 
  left_join(
    ref_sullivan |> select(sex, year, age, Lx_0, prev_0),
    by = c("sex", "year", "age")
  )

# counterfactuals:
hle_components <- dat_cf |>
  mutate(
    healthy      = 1 - prev,
    healthy_0    = 1 - prev_0,
    
    # Observed HLE: area mortality * area prevalence
    hle_obs_x    = Lx   * healthy,
    
    # Mortality-only counterfactual: area mortality * Basque prevalence
    hle_L_only_x = Lx   * healthy_0,
    
    # Prevalence-only counterfactual: Basque mortality * area prevalence
    hle_p_only_x = Lx_0 * healthy,
    
    # Fully reference schedule: Basque mortality * Basque prevalence
    hle_ref_x    = Lx_0 * healthy_0
  ) |>
  group_by(area, sex, year) |>
  summarise(
    HLE_obs    = sum(hle_obs_x, na.rm = TRUE),
    HLE_L_only = sum(hle_L_only_x, na.rm = TRUE),
    HLE_p_only = sum(hle_p_only_x, na.rm = TRUE),
    HLE_ref    = sum(hle_ref_x, na.rm = TRUE),
    .groups = "drop"
  )

shapley_var <- hle_components |>
  group_by(sex, year) |>
  summarise(
    V_obs    = var(HLE_obs, na.rm = TRUE),
    V_Lonly  = var(HLE_L_only, na.rm = TRUE),
    V_ponly  = var(HLE_p_only, na.rm = TRUE),
    # this is the shapley part
    V_prev = 0.5 * (V_ponly + V_obs - V_Lonly),
    V_LE   = 0.5 * (V_Lonly + V_obs - V_ponly),
    
    share_prev = V_prev / V_obs,
    share_LE   = V_LE / V_obs,
    
    check_sum = share_prev + share_LE,
    
    .groups = "drop"
  )