---
title: "Get started"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Get started}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
fig.height = 4.2, dpi = 110)
```
shapleydecomp answers one question about a simulation model: **how much of an
outcome is due to each of its inputs?** You give it a calibrated model and a
list of factors (treatment coverage, a prevention programme, a demographic
trend). It runs the model under every combination of observed and
counterfactual factor values, and splits the outcome among the factors with
Shapley values. The attributions always add up exactly to the total.
The package knows nothing about your model. It talks to the model through
four S3 methods, so any discrete-time simulation, in R or wrapped from C++,
can be decomposed once those methods exist.
```{r setup, message = FALSE}
library(shapleydecomp)
```
## 1. A model
This vignette uses a small SIR epidemic with two interventions:
- `beta_reduction`, a prevention programme that cuts transmission by up to 50%;
- `treatment_effect`, treatment that speeds recovery.
Both scale up linearly from 1990 to 2020. The model object holds the fixed
structure, and the parameters hold the inputs the decomposition will vary.
That split matters, because shapleydecomp only modifies `params`.
```{r model}
create_sir_model <- function(population = 1e5, years = 1990:2020) {
structure(
list(population = population, years = years,
S0 = population - 10, I0 = 10, R0 = 0),
class = "sir_model"
)
}
create_sir_params <- function(years = 1990:2020) {
yrs <- as.character(years)
list(
years = years,
base_beta = 0.3,
gamma = 0.1,
beta_reduction = setNames(seq(0, 0.5, length.out = length(years)), yrs),
treatment_effect = setNames(seq(0, 0.15, length.out = length(years)), yrs)
)
}
```
## 2. The interface
Four methods are required.
| Method | Job |
|---|---|
| `run_simulation()` | Run the model on a parameter set and return its raw output. |
| `summarize_results()` | Turn raw output into a `data.table` with a `year` column and one column per outcome. |
| `get_baseline_value()` | Say what "counterfactual" means for a factor. This is called once per factor. |
| `build_counterfactual()` | Return `params` with one factor switched to its counterfactual from a given year onward. |
```{r methods}
run_simulation.sir_model <- function(model, params, ...) {
n <- length(model$years)
S <- I <- R <- inf <- deaths <- numeric(n)
S[1] <- model$S0; I[1] <- model$I0; R[1] <- model$R0
for (t in 2:n) {
yr <- as.character(model$years[t])
beta <- params$base_beta * (1 - params$beta_reduction[[yr]])
gamma <- params$gamma + params$treatment_effect[[yr]]
inf[t] <- beta * S[t - 1] * I[t - 1] / model$population
deaths[t] <- 0.01 * I[t - 1]
rec <- gamma * I[t - 1]
S[t] <- S[t - 1] - inf[t]
I[t] <- I[t - 1] + inf[t] - rec - deaths[t]
R[t] <- R[t - 1] + rec
}
data.table::data.table(year = model$years, infections = inf, deaths = deaths)
}
summarize_results.sir_model <- function(model, sim_output, ...) sim_output
get_baseline_value.sir_model <- function(model, params, factor_name, cf_type,
baseline_year, ...) {
if (cf_type == CF_TOTAL_IMPACT) return(list(type = "zero"))
list(type = "change_since",
value = params[[factor_name]][[as.character(baseline_year)]])
}
build_counterfactual.sir_model <- function(model, params, factor_name,
decomp_year, cf_type,
baseline_values, ...) {
idx <- params$years >= decomp_year
params[[factor_name]][idx] <-
if (cf_type == CF_TOTAL_IMPACT) 0 else baseline_values[[factor_name]]$value
params
}
```
### Registering the methods
In a script run at the top level, defining `run_simulation.sir_model` is
enough, because S3 dispatch finds it in the global environment. When the
methods are defined somewhere else, such as inside a function, a vignette or a
test, register them so dispatch can find them:
```{r register}
.S3method("run_simulation", "sir_model", run_simulation.sir_model)
.S3method("summarize_results", "sir_model", summarize_results.sir_model)
.S3method("get_baseline_value", "sir_model", get_baseline_value.sir_model)
.S3method("build_counterfactual", "sir_model", build_counterfactual.sir_model)
```
If your model lives in a package, put `@export` on each method and roxygen
writes the `S3method()` lines into `NAMESPACE` for you. That is the
recommended home for a production model.
## 3. Validate
`validate_model_interface()` checks that the methods exist. Given parameters,
it also runs each method once and checks the shape of what comes back. Run it
before a long decomposition.
```{r validate}
model <- create_sir_model()
params <- create_sir_params()
validate_model_interface(
model, params,
factors = c("beta_reduction", "treatment_effect"),
summary_vars = c("infections", "deaths"),
baseline_year = 2000
)
```
## 4. Decompose
```{r decompose}
res <- shapley_decompose(
model, params,
factors = c("beta_reduction", "treatment_effect"),
summary_vars = c("infections", "deaths"),
decomp_years = 2000:2020,
cf_type = "total_impact",
verbose = FALSE
)
names(res)
head(res$decomp)
```
Each row of `res$decomp` is one decomposition year and one outcome. Each
factor column holds the Shapley value of *that year's* change in the factor:
the outcome with the factor on, minus the outcome with it at its
counterfactual. A negative value means the factor lowered the outcome.
The `deaths` rows are zero because deaths in this model follow infections
with a one-year lag, so a policy enacted in year $y$ cannot change deaths in
year $y$. Its effect on later years is recorded in `res$extended_decomp`.
`sequential_totals()` adds the steps up into each factor's total impact at
every year:
```{r totals}
tot <- sequential_totals(res)
tot[measure == "infections" & year %in% c(2005, 2010, 2020)]
```
```{r plot-totals, fig.alt = "Infections averted by each intervention per year, stacked, 2000 to 2020."}
inf <- tot[measure == "infections"]
op <- par(mar = c(3, 4.5, 1, 1))
averted <- -t(as.matrix(inf[, .(beta_reduction, treatment_effect)]))
colnames(averted) <- inf$year
barplot(averted, col = c("#3b6ea8", "#d08c2b"), border = NA, las = 1,
ylab = "Infections averted", ylim = c(0, max(colSums(averted)) * 1.1),
legend.text = c("Prevention (beta_reduction)", "Treatment"),
args.legend = list(x = "topleft", bty = "n"))
par(op)
```
The attributions are **exact**: at every year, the factor totals add up to
the observed outcome minus the outcome with both interventions off from 2000.
```{r exact}
off <- res$step_scenarios[step_year == 1999, .(year, infections)]
obs <- run_simulation(model, params)[, .(year, observed = infections)]
chk <- merge(merge(inf, obs, by = "year"), off, by = "year")
chk[, gap := beta_reduction + treatment_effect - (observed - infections)]
max(abs(chk$gap))
```
## 5. Two questions, two counterfactuals
`cf_type` decides what "the factor is off" means, and so which question the
decomposition answers.
- `"total_impact"` sets a factor to zero. It asks *how much did this factor
avert in total?*
- `"change_since"` holds a factor at its `baseline_year` value. It asks *how
much of the change since the baseline year is due to this factor?*
```{r change-since}
since <- shapley_decompose(
model, params,
factors = c("beta_reduction", "treatment_effect"),
summary_vars = "infections",
decomp_years = 2000:2020,
cf_type = "change_since",
baseline_year = 2000,
verbose = FALSE
)
sequential_totals(since)[year == 2020]
```
The `"change_since"` effects are smaller here because both programmes were
already running in 2000. Only their scale-up after 2000 is attributed.
## 6. Going further
- **Speed.** Implement `supports_warm_start()`, `get_warm_start_state()` and
`run_from_warm_start()`, and each counterfactual restarts from the observed
state the year before the decomposition year instead of from the start of
the simulation. On Unix, `n_cores` runs all the counterfactual
simulations, across every decomposition year, in parallel.
- **One-time setup.** `prepare_decomp()` runs once before the loop. Use it to
convert inputs (for example counts to coverage) or to cache the observed
simulation.
- **The method.** `vignette("methods")` covers the Shapley weights, the
equivalence with Das Gupta's standardization, how the sequential
decomposition is built, and which runs are skipped.