---
title: "Shapley Decomposition of Health‑Adjusted Life Expectancy with shapleyHALE"
author: "Jun‑Yan Xi; Sun Yat-sen University; Email: xijy3@mail2.sysu.edu.cn; ORCID: https://orcid.org/0000-0001-7473-7783"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Shapley Decomposition of Health‑Adjusted Life Expectancy with shapleyHALE}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## 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:

```{r, eval = FALSE}
# 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).

```{r, eval = FALSE}
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:

```{r, eval = FALSE}
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:

```{r, eval = FALSE}
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.

```{r, 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])
```

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.

```{r, 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))
```

## 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.

```{r, 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]
```

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

```{r, eval = FALSE}
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.

```{r, 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)
```

## 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.

```{r, 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.
```

## 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.