## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 5,
  message = FALSE,
  warning = FALSE
)

## ----eval = FALSE-------------------------------------------------------------
# # install.packages("devtools")
# # devtools::install_github("username/shapleyHALE")

## ----eval = FALSE-------------------------------------------------------------
# library(shapleyHALE)
# library(data.table)

## ----eval = FALSE-------------------------------------------------------------
# data("demo_gbd")
# data("age_info_demo")
# 
# str(demo_gbd)
# head(demo_gbd)

## ----eval = FALSE-------------------------------------------------------------
# age_info_demo

## ----eval = FALSE-------------------------------------------------------------
# cause_ids <- unique(demo_gbd$cause_id)
# age_ids   <- age_info_demo$age_id
# 
# mx_1990   <- create_rate_matrix(demo_gbd[measure_id == 1], 1990, cause_ids, age_ids)
# mx_2023   <- create_rate_matrix(demo_gbd[measure_id == 1], 2023, cause_ids, age_ids)
# yld_1990  <- create_rate_matrix(demo_gbd[measure_id == 3], 1990, cause_ids, age_ids)
# yld_2023  <- create_rate_matrix(demo_gbd[measure_id == 3], 2023, cause_ids, age_ids)
# 
# # Quick check
# dim(mx_1990)
# head(mx_1990[, 1:5])

## ----eval = FALSE-------------------------------------------------------------
# total_mx_1990  <- rowSums(mx_1990)
# total_yld_1990 <- rowSums(yld_1990)
# total_mx_2023  <- rowSums(mx_2023)
# total_yld_2023 <- rowSums(yld_2023)
# 
# hale_1990 <- compute_hale(total_mx_1990, total_yld_1990, age_info_demo)
# hale_2023 <- compute_hale(total_mx_2023, total_yld_2023, age_info_demo)
# 
# c(HALE_1990 = round(hale_1990$HALE0, 2),
#   HALE_2023 = round(hale_2023$HALE0, 2),
#   change    = round(hale_2023$HALE0 - hale_1990$HALE0, 2))

## ----eval = FALSE-------------------------------------------------------------
# set.seed(2026) # for reproducibility
# decomp <- run_shapley_decomp(
#   base_mx    = mx_1990,
#   base_yld   = yld_1990,
#   target_mx  = mx_2023,
#   target_yld = yld_2023,
#   causes     = cause_ids,
#   ages       = age_info_demo,
#   n_perm     = 500  # In practice use 5000 or more
# )
# 
# # Total contributions by cause
# total_contrib <- data.table(
#   cause_id = names(decomp$total_effect),
#   total    = decomp$total_effect,
#   death    = decomp$death_effect,
#   disability = decomp$disability_effect
# )
# total_contrib[order(-total)][1:5]

## ----eval = FALSE-------------------------------------------------------------
# sum(decomp$total_effect)
# decomp$total_change

## ----eval = FALSE-------------------------------------------------------------
# inter_df <- compute_shapley_interactions(
#   base_mx    = mx_1990,
#   base_yld   = yld_1990,
#   target_mx  = mx_2023,
#   target_yld = yld_2023,
#   causes     = cause_ids,
#   ages       = age_info_demo,
#   top_n      = 10,
#   n_perm     = 200,       # small for demo, use 1000+ in practice
#   point_est  = decomp$total_effect
# )
# 
# # Display top interactions
# head(inter_df[order(-abs(interaction))], 10)

## ----eval = FALSE-------------------------------------------------------------
# boot_res <- bootstrap_shapley(
#   deaths_data  = demo_gbd[measure_id == 1],
#   ylds_data    = demo_gbd[measure_id == 3],
#   base_year    = 1990,
#   target_year  = 2023,
#   cause_ids    = cause_ids,
#   age_ids      = age_ids,
#   ages_info    = age_info_demo,
#   n_boot       = 200,        # use 1000+ in real analysis
#   n_perm_inner = 100,        # use 5000+ in real analysis
#   parallel     = FALSE
# )
# 
# # Total effect summary for top causes
# boot_res$boot_total[order(-abs(point_est))][1:5]
# The output contains point estimates and bootstrap‑based 95% uncertainty intervals (lower, upper) for total, death, and disability effects.