Shapley Decomposition of Health‑Adjusted Life Expectancy with shapleyHALE

Introduction

Health‑adjusted life expectancy (HALE) integrates mortality and morbidity into a single summary measure of population health. Conventional decomposition methods assume that diseases act independently, which can obscure important interactions in an era of rising multimorbidity. The shapleyHALE package provides a Shapley‑value‑based framework that fully accounts for higher‑order interactions among causes, separates mortality from disability effects, and quantifies pairwise synergies or antagonisms. This vignette demonstrates the package’s core functionality using simulated Global Burden of Disease (GBD) data.

Installation

You can install the development version from GitHub (or locally) with:

# install.packages("devtools")
# devtools::install_github("username/shapleyHALE")

After installation, load the package along with data.table (which is used internally and recommended for data handling).

library(shapleyHALE)
library(data.table)

Data structure

The package expects cause‑specific mortality rates and years lived with disability (YLD) rates, typically obtained from the GBD study. The built‑in demo dataset demo_gbd mimics such data:

data("demo_gbd")
data("age_info_demo")

str(demo_gbd)
head(demo_gbd)

measure_id: 1 = mortality, 3 = YLD
year: 1990 or 2023
cause_id: character cause identifier (e.g., c1, c2)
age_id: numeric age group ID, matching age_info_demo
val, lower, upper: point estimate and 95% uncertainty interval of the rate

The dataset age_info_demo provides age group boundaries and widths:

age_info_demo

Step 1: Creating rate matrices

The function create_rate_matrix() extracts and reshapes cause‑age rates into a matrix (rows = ages, columns = causes). Here we build the baseline (1990) and target (2023) matrices for mortality and YLD.

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])

Step 2: HALE computation (Sullivan method) compute_hale() implements the standard Sullivan method. It requires age‑specific total mortality and YLD rates, plus the age information table.

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))

Step 3: Shapley decomposition

run_shapley_decomp() decomposes the total HALE change into additive contributions of individual causes, further separating mortality‑driven and disability‑driven components.

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]

The sum of all Shapley values equals the total HALE change (up to Monte Carlo error):

sum(decomp$total_effect)
decomp$total_change

Step 4: Second‑order interaction indices

compute_shapley_interactions() quantifies pairwise synergies (positive values) and antagonisms (negative values) among the leading causes.

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)

Step 5: Uncertainty propagation via parametric bootstrap

bootstrap_shapley() uses a parametric bootstrap (triangular distribution from GBD uncertainty intervals) to obtain 95% confidence intervals for all Shapley values.

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.

Conclusion

This vignette illustrated the core workflow of the shapleyHALE package:

building rate matrices from GBD‑style data,
computing HALE via the Sullivan method,
decomposing changes with Shapley values,
quantifying pairwise interactions,
and propagating uncertainty with a parametric bootstrap.

The framework is fully portable and can be applied to any population for which age‑ and cause‑specific rates are available. For further details on the methodology, please refer to the accompanying paper.