| Title: | Shapley Value Decomposition for Discrete-Time Simulation Models |
|---|---|
| Description: | Attributes changes in the outcomes of a discrete-time simulation model to its input factors using Shapley values (Shapley 1953) <doi:10.1515/9781400881970-018>, which for binary factors coincide with the symmetric standardization of Das Gupta (1993) <https://www2.census.gov/library/publications/1993/demographics/p23-186.pdf>. Any calibrated model can be decomposed by implementing four S3 methods that run the model, summarize its output, and construct counterfactual inputs. Counterfactuals either remove a factor entirely or hold it at its value in a baseline year, and a sequential decomposition attributes each year's policy conditional on the past. |
| Authors: | Austin Carter [aut, cre] |
| Maintainer: | Austin Carter <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.0 |
| Built: | 2026-11-03 23:50:47 UTC |
| Source: | https://github.com/cran/shapleydecomp |
Given the observed parameters, produce a modified parameter set where the
specified factor is set to its counterfactual value from
decomp_year onward.
build_counterfactual( model, params, factor_name, decomp_year, cf_type, baseline_values, ... )build_counterfactual( model, params, factor_name, decomp_year, cf_type, baseline_values, ... )
model |
An S3 model object |
params |
Observed model parameters |
factor_name |
Character, name of the factor to modify |
decomp_year |
Integer, the year from which the counterfactual applies |
cf_type |
Character, one of |
baseline_values |
Named list of baseline values for each factor
(pre-computed by |
... |
Additional arguments |
Modified model parameters with the factor set to its counterfactual value.
Other model interface:
get_baseline_value(),
prepare_decomp(),
run_simulation(),
summarize_results(),
validate_model_interface(),
warm_start
Two types of counterfactual analysis are supported:
CF_CHANGE_SINCE: The counterfactual value is the factor's
value at the baseline year (typically the first decomposition year).
This measures the impact of changes in the factor since that
year.
CF_TOTAL_IMPACT: The counterfactual value is zero (or the
model's notion of "no intervention"). This measures the total
impact of the factor, not just changes since a reference point.
CF_CHANGE_SINCE CF_TOTAL_IMPACTCF_CHANGE_SINCE CF_TOTAL_IMPACT
An object of class character of length 1.
An object of class character of length 1.
A length-one character string: "change_since" or
"total_impact", the values accepted by the cf_type
argument of shapley_decompose.
CF_CHANGE_SINCE CF_TOTAL_IMPACTCF_CHANGE_SINCE CF_TOTAL_IMPACT
Decompose all summary variables for a single group
decomp_algebra(data, params)decomp_algebra(data, params)
data |
A |
params |
Parameter list (see |
A data.table with one row per summary variable: a
measure column and one column per factor.
# Outcomes for every on/off combination of two factors, where the # factors interact: y = 10 + 5a + 3b + 2ab rt <- make_run_table(c("a", "b")) rt[, y := 10 + 5 * a + 3 * b + 2 * a * b] # The interaction term is split evenly: a gets 5 + 1, b gets 3 + 1 decomp_algebra(rt, list(summary_vars = "y", decomp_vars = c("a", "b"), decomp_method = "shapley"))# Outcomes for every on/off combination of two factors, where the # factors interact: y = 10 + 5a + 3b + 2ab rt <- make_run_table(c("a", "b")) rt[, y := 10 + 5 * a + 3 * b + 2 * a * b] # The interaction term is split evenly: a gets 5 + 1, b gets 3 + 1 decomp_algebra(rt, list(summary_vars = "y", decomp_vars = c("a", "b"), decomp_method = "shapley"))
For each factor, computes the weighted average marginal effect of switching the factor from 0 (counterfactual) to 1 (observed), averaging over all combinations of the other factors.
decomp_pair(dt, params, summary_var)decomp_pair(dt, params, summary_var)
dt |
A |
params |
Parameter list with |
summary_var |
Character, the name of the summary variable column. |
A one-row data.table with a measure column and one
column per factor.
rt <- make_run_table(c("a", "b")) rt[, y := 10 + 5 * a + 3 * b + 2 * a * b] decomp_pair(rt, list(decomp_vars = c("a", "b"), decomp_method = "shapley"), "y")rt <- make_run_table(c("a", "b")) rt[, y := 10 + 5 * a + 3 * b + 2 * a * b] decomp_pair(rt, list(decomp_vars = c("a", "b"), decomp_method = "shapley"), "y")
Returns the value of a factor at the baseline year (for
CF_CHANGE_SINCE) or the zero/null value (for
CF_TOTAL_IMPACT). Called once per factor before the decomposition
loop begins.
get_baseline_value(model, params, factor_name, cf_type, baseline_year, ...)get_baseline_value(model, params, factor_name, cf_type, baseline_year, ...)
model |
An S3 model object |
params |
Observed model parameters |
factor_name |
Character, name of the factor |
cf_type |
Character, |
baseline_year |
Integer, the reference year for |
... |
Additional arguments |
The baseline value for the factor (format is model-specific).
Other model interface:
build_counterfactual(),
prepare_decomp(),
run_simulation(),
summarize_results(),
validate_model_interface(),
warm_start
Creates a data.table with one row per counterfactual combination.
Each factor column contains 0 (counterfactual) or 1 (observed), encoding
whether that factor is "on" or "off" for that run.
make_run_table(factors)make_run_table(factors)
factors |
Character vector of factor names. |
A data.table with 2^n rows: run_id (integer), and one
column per factor with values 0 or 1.
make_run_table(c("treatment", "prevention", "demography"))make_run_table(c("treatment", "prevention", "demography"))
One-time preparation step called before the decomposition loop. Use this for expensive setup operations like running an initial simulation to extract denominators, converting counts to percentages, etc.
prepare_decomp(model, params, ...)prepare_decomp(model, params, ...)
model |
An S3 model object |
params |
Model parameters |
... |
Additional arguments |
The default method returns the model and params unchanged.
A list with components:
model |
Possibly modified model object |
params |
Possibly modified parameters |
full_sim |
Optional: output of a full (observed) simulation, used for warm starts and as the "all-on" reference |
full_summary |
Optional: summary of the full simulation |
Other model interface:
build_counterfactual(),
get_baseline_value(),
run_simulation(),
summarize_results(),
validate_model_interface(),
warm_start
Applies the Shapley decomposition to each group × decomp-year in the data.
run_decomp(data, params)run_decomp(data, params)
data |
A |
params |
A list with at least:
|
A data.table with one row per group × decomp-year ×
measure, and one column per factor giving its decomposition effect.
# Outcomes for every on/off combination of two factors, where the # factors interact: y = 10 + 5a + 3b + 2ab rt <- make_run_table(c("a", "b")) rt[, y := 10 + 5 * a + 3 * b + 2 * a * b] rt[, decomp_year := 2000L] run_decomp(rt, list(summary_vars = "y", decomp_vars = c("a", "b"), decomp_method = "shapley"))# Outcomes for every on/off combination of two factors, where the # factors interact: y = 10 + 5a + 3b + 2ab rt <- make_run_table(c("a", "b")) rt[, y := 10 + 5 * a + 3 * b + 2 * a * b] rt[, decomp_year := 2000L] run_decomp(rt, list(summary_vars = "y", decomp_vars = c("a", "b"), decomp_method = "shapley"))
Execute the discrete-time simulation model and return raw simulation output. This is called once for each counterfactual combination at each decomposition year.
run_simulation(model, params, ...)run_simulation(model, params, ...)
model |
An S3 model object (class determines dispatch) |
params |
Model parameters (possibly modified for counterfactual) |
... |
Additional arguments passed to the method |
Raw simulation output (model-specific format). Will be passed to
summarize_results().
Other model interface:
build_counterfactual(),
get_baseline_value(),
prepare_decomp(),
summarize_results(),
validate_model_interface(),
warm_start
Sums the step attributions of a shapley_decompose() run across
decomposition steps: the total impact of factor i at year T is
, where is the
Shapley value at year T of the year-y policy (extended_decomp). Summed
over factors this telescopes to (observed) - (all policies off from the first
decomp year) at every T, so it is only meaningful when decomp_years
is contiguous.
sequential_totals(result)sequential_totals(result)
result |
The list returned by |
data.table with one row per group x year x measure and one column per factor.
shapley_decompose, whose examples call this function.
This is the main user-facing function. Given a model object that implements the shapleydecomp S3 interface, it performs a complete factorial decomposition of the specified factors across all decomposition years.
shapley_decompose( model, params, factors, summary_vars, decomp_years, method = "shapley", cf_type = CF_CHANGE_SINCE, baseline_year = min(decomp_years), extended = TRUE, step_scenarios = TRUE, verbose = TRUE, n_cores = 1L, group_vars = character(), ... )shapley_decompose( model, params, factors, summary_vars, decomp_years, method = "shapley", cf_type = CF_CHANGE_SINCE, baseline_year = min(decomp_years), extended = TRUE, step_scenarios = TRUE, verbose = TRUE, n_cores = 1L, group_vars = character(), ... )
model |
An S3 model object. Must have methods for
|
params |
Model parameters (model-specific list/object). |
factors |
Character vector of factor names to decompose. |
summary_vars |
Character vector of summary variable names (columns
in the output of |
decomp_years |
Integer vector of years at which to perform
decomposition (e.g., |
method |
Character, decomposition method. Only
|
cf_type |
Character, counterfactual type:
|
baseline_year |
Integer, the reference year for |
extended |
Logical, whether to compute the extended factor decomposition (effects at all output years >= decomp year). Default TRUE. |
step_scenarios |
Logical, whether to compute step scenarios (all-on time series per decomp year). Default TRUE. |
verbose |
Logical, print progress messages. Default TRUE. |
n_cores |
Integer, number of cores for parallel counterfactual
simulation runs. Default 1 (sequential). Values > 1 run every
counterfactual simulation, across all decomp years, in one
|
group_vars |
Character vector of columns in the output of
|
... |
Additional arguments passed to model methods. |
The decomposition is sequential (Carter 2025, Quantifying averted
burden from infectious disease control policies, PhD dissertation,
University of Washington, https://hdl.handle.net/1773/53310). At
each decomp year y, every policy before y is at its observed value, every
policy after y is at its counterfactual, and only the year-y policy varies
across the 2^n coalitions. The attribution at step y is therefore the
marginal impact of the year-y policy, conditional on the past and undiluted
by the future. extended_decomp records that impact at every later
year, and sequential_totals() sums the steps into each factor's
total impact at every year.
A list with components:
data.table — standard year-by-year decomposition (one row per decomp year × measure, factor effects in columns)
data.table — all-on time series per decomp year
(columns: step_year, year, plus summary vars). NULL if
step_scenarios = FALSE.
data.table — extended decomposition (factor
effects at all years >= decomp year). NULL if extended = FALSE.
data.table — the factorial design table used.
data.table — per decomp year: the number of coalition runs, how many were simulated, how many were reused from the previous step, and the factors that were null players (their counterfactual left the params unchanged, so their effect is 0).
List of decomposition parameters used.
vignette("shapleydecomp") for a worked example, and
validate_model_interface to check a model first.
# A two-factor SIR model shipped with the package sir <- new.env() sys.source(system.file("examples", "sir_model.R", package = "shapleydecomp"), envir = sir) model <- sir$create_sir_model(years = 1990:2010) params <- sir$create_sir_params(years = 1990:2010) res <- shapley_decompose( model, params, factors = c("beta_reduction", "treatment_effect"), summary_vars = c("infections", "deaths"), decomp_years = 2000:2010, cf_type = "total_impact", verbose = FALSE ) res$decomp # Each factor's total impact at every year sequential_totals(res)# A two-factor SIR model shipped with the package sir <- new.env() sys.source(system.file("examples", "sir_model.R", package = "shapleydecomp"), envir = sir) model <- sir$create_sir_model(years = 1990:2010) params <- sir$create_sir_params(years = 1990:2010) res <- shapley_decompose( model, params, factors = c("beta_reduction", "treatment_effect"), summary_vars = c("infections", "deaths"), decomp_years = 2000:2010, cf_type = "total_impact", verbose = FALSE ) res$decomp # Each factor's total impact at every year sequential_totals(res)
Extract summary measures (e.g., deaths, infections) from raw simulation
output. Must return a data.table with at minimum a year
column and one column per summary variable specified in the decomposition.
summarize_results(model, sim_output, ...)summarize_results(model, sim_output, ...)
model |
An S3 model object |
sim_output |
Raw simulation output from |
... |
Additional arguments |
A data.table with columns: year, plus one column per
summary variable (e.g., deaths, infections).
Other model interface:
build_counterfactual(),
get_baseline_value(),
prepare_decomp(),
run_simulation(),
validate_model_interface(),
warm_start
Checks that the required S3 methods exist for the given model class and runs a minimal smoke test if parameters and factors are provided.
validate_model_interface( model, params = NULL, factors = NULL, summary_vars = NULL, baseline_year = NULL, verbose = TRUE )validate_model_interface( model, params = NULL, factors = NULL, summary_vars = NULL, baseline_year = NULL, verbose = TRUE )
model |
An S3 model object |
params |
Optional: model parameters for smoke test |
factors |
Optional: character vector of factor names |
summary_vars |
Optional: character vector of summary variable names |
baseline_year |
Optional: integer baseline year |
verbose |
Logical, print diagnostic messages |
Invisible TRUE if all checks pass; stops with an error otherwise.
Other model interface:
build_counterfactual(),
get_baseline_value(),
prepare_decomp(),
run_simulation(),
summarize_results(),
warm_start
sir <- new.env() sys.source(system.file("examples", "sir_model.R", package = "shapleydecomp"), envir = sir) validate_model_interface( sir$create_sir_model(), sir$create_sir_params(), factors = c("beta_reduction", "treatment_effect"), summary_vars = c("infections", "deaths"), baseline_year = 2000 )sir <- new.env() sys.source(system.file("examples", "sir_model.R", package = "shapleydecomp"), envir = sir) validate_model_interface( sir$create_sir_model(), sir$create_sir_params(), factors = c("beta_reduction", "treatment_effect"), summary_vars = c("infections", "deaths"), baseline_year = 2000 )
Optional methods that let shapley_decompose() skip re-simulating
the years before each decomposition year. Up to year every
coalition run is identical to the observed run, so a model that can
restart from saved state only needs to simulate from onward.
supports_warm_start(model, ...) get_warm_start_state(model, sim_output, year, ...) run_from_warm_start(model, warm_state, params, output_years, ...)supports_warm_start(model, ...) get_warm_start_state(model, sim_output, year, ...) run_from_warm_start(model, warm_state, params, output_years, ...)
model |
An S3 model object |
... |
Additional arguments. |
sim_output |
Observed simulation output from |
year |
Integer, the year whose end-of-year state to extract |
warm_state |
Observed simulation output to restart from (the
|
params |
Model parameters (possibly counterfactual) |
output_years |
Integer vector of years the output must cover |
When supports_warm_start() returns TRUE,
shapley_decompose() calls
run_from_warm_start(model, warm_state, params, output_years,
warm_start_year = y - 1, ...) for every counterfactual run instead of
run_simulation(). warm_state is the observed simulation
output (the full_sim returned by prepare_decomp()).
run_from_warm_start() extracts the state at the end of
warm_start_year from it, typically by calling
get_warm_start_state(model, warm_state, warm_start_year), and
simulates forward with the counterfactual params.
shapley_decompose() itself never calls get_warm_start_state().
It is the conventional helper for that extraction step.
The default supports_warm_start() returns FALSE. The default
get_warm_start_state() and run_from_warm_start() methods
throw an error.
supports_warm_start(): logical.
get_warm_start_state(): a model state object in a model-specific
format. run_from_warm_start(): raw simulation output in the same
format as run_simulation(), covering at least
output_years.
Other model interface:
build_counterfactual(),
get_baseline_value(),
prepare_decomp(),
run_simulation(),
summarize_results(),
validate_model_interface()
# A stock that decays by 20% a year and receives input x: # s_t = 0.8 * s_{t-1} + x_t run_simulation.stock_model <- function(model, params, ...) { s <- as.numeric(stats::filter(params$x, 0.8, method = "recursive")) data.table::data.table(year = params$years, stock = s) } summarize_results.stock_model <- function(model, sim_output, ...) sim_output get_baseline_value.stock_model <- function(model, params, factor_name, cf_type, baseline_year, ...) { list(type = "zero") } build_counterfactual.stock_model <- function(model, params, factor_name, decomp_year, cf_type, baseline_values, ...) { params$x[params$years >= decomp_year] <- 0 params } # Warm start: restart from the observed stock at the end of # warm_start_year and simulate only the remaining years supports_warm_start.stock_model <- function(model, ...) TRUE get_warm_start_state.stock_model <- function(model, sim_output, year, ...) { s <- sim_output$stock[sim_output$year == year] # Before the first simulated year there is no stock yet if (length(s) == 0) 0 else s } run_from_warm_start.stock_model <- function(model, warm_state, params, output_years, warm_start_year, ...) { s0 <- get_warm_start_state(model, warm_state, warm_start_year) keep <- params$years > warm_start_year s <- as.numeric(stats::filter(params$x[keep], 0.8, method = "recursive", init = s0)) data.table::data.table(year = params$years[keep], stock = s) } for (g in c("run_simulation", "summarize_results", "get_baseline_value", "build_counterfactual", "supports_warm_start", "get_warm_start_state", "run_from_warm_start")) { .S3method(g, "stock_model", get(paste0(g, ".stock_model"))) } params <- list(years = 2000:2010, x = rep(1, 11)) warm <- shapley_decompose(structure(list(), class = "stock_model"), params, factors = "x", summary_vars = "stock", decomp_years = 2000:2010, cf_type = "total_impact", verbose = FALSE) # The same model without warm starts gives the same answer .S3method("supports_warm_start", "stock_cold", function(model, ...) FALSE) cold <- shapley_decompose(structure(list(), class = c("stock_cold", "stock_model")), params, factors = "x", summary_vars = "stock", decomp_years = 2000:2010, cf_type = "total_impact", verbose = FALSE) all.equal(warm$decomp, cold$decomp)# A stock that decays by 20% a year and receives input x: # s_t = 0.8 * s_{t-1} + x_t run_simulation.stock_model <- function(model, params, ...) { s <- as.numeric(stats::filter(params$x, 0.8, method = "recursive")) data.table::data.table(year = params$years, stock = s) } summarize_results.stock_model <- function(model, sim_output, ...) sim_output get_baseline_value.stock_model <- function(model, params, factor_name, cf_type, baseline_year, ...) { list(type = "zero") } build_counterfactual.stock_model <- function(model, params, factor_name, decomp_year, cf_type, baseline_values, ...) { params$x[params$years >= decomp_year] <- 0 params } # Warm start: restart from the observed stock at the end of # warm_start_year and simulate only the remaining years supports_warm_start.stock_model <- function(model, ...) TRUE get_warm_start_state.stock_model <- function(model, sim_output, year, ...) { s <- sim_output$stock[sim_output$year == year] # Before the first simulated year there is no stock yet if (length(s) == 0) 0 else s } run_from_warm_start.stock_model <- function(model, warm_state, params, output_years, warm_start_year, ...) { s0 <- get_warm_start_state(model, warm_state, warm_start_year) keep <- params$years > warm_start_year s <- as.numeric(stats::filter(params$x[keep], 0.8, method = "recursive", init = s0)) data.table::data.table(year = params$years[keep], stock = s) } for (g in c("run_simulation", "summarize_results", "get_baseline_value", "build_counterfactual", "supports_warm_start", "get_warm_start_state", "run_from_warm_start")) { .S3method(g, "stock_model", get(paste0(g, ".stock_model"))) } params <- list(years = 2000:2010, x = rep(1, 11)) warm <- shapley_decompose(structure(list(), class = "stock_model"), params, factors = "x", summary_vars = "stock", decomp_years = 2000:2010, cf_type = "total_impact", verbose = FALSE) # The same model without warm starts gives the same answer .S3method("supports_warm_start", "stock_cold", function(model, ...) FALSE) cold <- shapley_decompose(structure(list(), class = c("stock_cold", "stock_model")), params, factors = "x", summary_vars = "stock", decomp_years = 2000:2010, cf_type = "total_impact", verbose = FALSE) all.equal(warm$decomp, cold$decomp)