Take a simulation whose outcome at time \(t\) depends on \(n\) factors, \(Y_t = f_t(x_1, \ldots, x_n)\). Each factor has an observed trajectory \(x_i^1\), the one the calibrated model uses, and a counterfactual trajectory \(x_i^0\). The goal is to split the total difference
\[ \Delta Y_t = f_t(x_1^1, \ldots, x_n^1) - f_t(x_1^0, \ldots, x_n^0) \]
into contributions \(\delta_i\) with \(\sum_i \delta_i = \Delta Y_t\) exactly.
For \(n\) factors there are \(2^n\) on/off combinations.
make_run_table() enumerates them, and
shapley_decompose() runs the model once per row.
Following Shapley (1953), treat the set \(S\) of factors held at their observed values as a coalition, with value \(v(S) = f_t(x_S^1, x_{N \setminus S}^0)\). The Shapley value of factor \(i\) is its marginal effect averaged over every coalition of the other factors:
\[ \phi_i = \sum_{S \subseteq N \setminus \{i\}} \frac{|S|!\,(n - |S| - 1)!}{n!}\,\bigl[v(S \cup \{i\}) - v(S)\bigr]. \]
It is the only allocation that is efficient (the shares sum to \(\Delta Y_t\)), symmetric (interchangeable factors get equal shares), additive, and gives zero to a factor that never changes the outcome.
Interactions are split evenly. In the example below, the outcome is
\(10 + 5a + 3b + 2ab\). Factor
a gets its main effect plus half of the interaction, and
factor b gets the same treatment:
rt <- make_run_table(c("a", "b"))
rt[, y := 10 + 5 * a + 3 * b + 2 * a * b]
decomp_algebra(rt, list(summary_vars = "y", decomp_vars = c("a", "b"),
decomp_method = "shapley"))
#> measure a b
#> <char> <num> <num>
#> 1: y 6 4Das Gupta’s (1993) symmetric standardization, the classical demographic tool, weights each marginal effect by \(1 / \bigl[n \binom{n-1}{k}\bigr]\), where \(k = |S|\). Expanding the binomial gives
\[ \frac{1}{n \binom{n-1}{k}} = \frac{k!\,(n-1-k)!}{n \cdot (n-1)!} = \frac{k!\,(n-k-1)!}{n!}, \]
which is the Shapley weight. For binary factors the two methods are
identical, so the package implements one weight formula.
method = "das_gupta" is accepted as a deprecated alias for
"shapley".
What \(x_i^0\) is determines the question being answered.
cf_type |
\(x_i^0(t)\) | Question |
|---|---|---|
"total_impact" |
\(0\) (no programme) | How much did this factor avert, in total? |
"change_since" |
\(x_i(t_0)\), held at
baseline_year |
How much of the change since \(t_0\) is due to this factor? |
Your get_baseline_value() and
build_counterfactual() methods implement these definitions
for your model’s parameters. For example, “zero” for a coverage factor
might mean no coverage, while for a rate it might mean a reference
rate.
Every factor is either observed or counterfactual, but the
counterfactual does not have to be zero. To ask what coverage
above some reference achieved, such as half of observed
coverage or a target level, define that reference as the counterfactual
in your model’s get_baseline_value() and
build_counterfactual(). The decomposition is unchanged and
still needs \(2^n\) runs per step. It
attributes the difference between the observed world and the reference
world.
Two related questions are different methods, not settings of this one. Splitting credit along a continuous path from off to observed is the Aumann-Shapley value, which differs from the Shapley value only when factors interact non-linearly. Asking what would happen at 50% coverage is a scenario projection, not an attribution.
Simulations carry state forward. Infections averted in 2005 change
the susceptible pool in 2010. A single decomposition at a late year
would credit every past policy at once, so
shapley_decompose() works year by year, using the
sequential Shapley construction of Carter (2025).
At each decomposition year \(y\):
The Shapley value at step \(y\) is
therefore the marginal impact of the year-\(y\) policy, conditional on the past and
undiluted by the future. Because the model is simulated forward from
\(y\), extended_decomp
records that step’s effect at every later year \(T \ge y\), written \(\phi_i^{(y)}(T)\).
sequential_totals() sums the steps:
\[ \text{total}_i(T) = \sum_{y \le T} \phi_i^{(y)}(T). \]
Summed over factors, this telescopes to observed minus all-off
from the first decomposition year at every \(T\), so the totals stay exact. That only
holds when decomp_years is contiguous, and
sequential_totals() checks it.
Exact Shapley values need all \(2^n\) coalitions at every step, but two kinds of coalition repeat a run that already exists:
shapley_decompose() simulates each distinct scenario
once. It confirms that the parameters are identical before reusing a
run, so a reused run always equals the one it replaces.
run_summary in the result reports, for each step, how many
runs were simulated and which factors were null.
Shapley LS (1953). A value for n-person games. In Kuhn HW, Tucker AW (eds), Contributions to the Theory of Games II, 307–317. Princeton University Press. https://doi.org/10.1515/9781400881970-018
Das Gupta P (1993). Standardization and Decomposition of Rates: A User’s Manual. Current Population Reports P23-186. U.S. Bureau of the Census. https://www2.census.gov/library/publications/1993/demographics/p23-186.pdf
Carter A (2025). Quantifying averted burden from infectious disease control policies. PhD dissertation, University of Washington. https://hdl.handle.net/1773/53310
Carter A. A model-agnostic framework for factorial decomposition of discrete-time simulation outcomes, with a proof of equivalence between Das Gupta and Shapley value decomposition. Manuscript in preparation.