---
title: "Introduction to `{whatifbandit}`"
bibliography: "REFERENCES.bib"
vignette: >
  %\VignetteIndexEntry{Introduction to `{whatifbandit}`}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
knitr:
  opts_chunk:
    collapse: true
    comment: '#>'
message: FALSE
warnings: FALSE
---

# Overview

Traditional randomized controlled trials (RCTs) are the workhorse of causal inference in the social
sciences. Randomization of treatments provides unconfoundedness, and the fixed, equal probability of treatment assignment means causal effects can be validly estimated even with the simplest, most intuitive estimators. However, these designs come with two distinct costs:

- The number of treatments to be tested must be limited or pruned in advance by researchers to maintain the experiment's power given a fixed sample size.
- Participants continue to be assigned weaker treatments, even after incoming data reveal it may be more optimal for them to receive a different one.

In a two- or three-armed, one-shot experiment, these costs are negligible, but academics and policy experimenters frequently run trials over long periods of time in which participants' welfare is directly at stake, making these costs non-negligible.

Response-adaptive designs address these costs directly by tying assignment probabilities to the
evidence a trial has already produced rather than fixing them in advance. As the data accumulates, the
probability of assignment shifts toward treatments that appear to be working, and away from ones
that aren't. This gives researchers room to test more treatments up front, since underperforming
arms will be phased out or dropped. At the same time, participants become concentrated on the treatments that
actually work as the trial unfolds. For simplicity, response-adaptive designs will be referred to as
simply "adaptive designs" or "adaptive trials" for the remainder of the vignette, even though
response adaptation is only one form of adaptive experimental design more broadly.

Although there are many ways to conduct a response-adaptive design, `{whatifbandit}` uses multi-armed bandit (MAB) algorithms, since they provide a clear framework for navigating the exploration-versus-exploitation trade-off in treatment assignment.

The central question `{whatifbandit}` is ***What if my experiment had been run as a
bandit trial instead?*** The package name reflects this directly: *bandit* for
multi-armed bandit, and *whatif* for the counterfactual question the package
lets you explore. This vignette walks through the package's two main modes of use, its supporting
design features, and how to interpret the estimates it produces.

# Two Ways to Explore a Bandit Design

`{whatifbandit}` supports two distinct workflows, depending on whether you
already have trial data in hand or are planning a trial that hasn't occurred
yet.

## 1. Resimulation

If you have data from a trial that already ran as a standard RCT,
`{whatifbandit}` can resimulate that trial as though a bandit algorithm had
governed assignment instead. Using the outcomes you actually observed, the
package imputes new outcomes for participants under the counterfactual adaptive
assignment path, reconstructing a plausible "what-if" version of your
experiment. This is useful for retrospectively asking whether an adaptive design (in the same
experimental conditions) would have delivered better outcomes for participants, more precise
estimates, or both.


## 2. Simulation from Scratch

If you are designing a trial that hasn't been run yet, or doing other forms of testing, `{whatifbandit}` can
simulate a bandit trial directly from a set of population parameters you supply
(for example, assumed true success probabilities for each arm). This is useful
for:

- power analysis for a planned adaptive trial,
- comparing how different bandit algorithms would behave under your assumed
 conditions, and
- exploring design choices (batch size, enforced exploration, number of arms)
 before collecting any real data.

## Scope and Limitations

`{whatifbandit}` is currently built specifically for experiments with **binary
outcomes** (success/failure), without covariates (except for discrete blocking variables). Support for other outcome types may be added in
future releases, but any design considered in this vignette assumes a binary
response.

```{r}
#| output: FALSE
library(whatifbandit)
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(forcats)
options(scipen = 9999)
```

# Example Data: `tanf`

Examples throughout this vignette use the `tanf` dataset bundled with the
package. It contains anonymized Temporary Assistance for Needy Families (TANF) recertification
records from a field experiment which tested the impact of specific notification letters on
recertification results [@moore2022]. These new notification letters were designed to reduce the
cognitive load required to understand the instructions, and were touted as a solution to reduce the
program's notorious churn (where welfare recipients let their benefits lapse and reapply to the
program, as opposed to recertifying their benefits) [@moore2022]. The experiment was conducted with
real TANF recipients in Washington, D.C.

```{r}
data("tanf")
head(tanf)
glimpse(tanf)
```

# Core Package Features

To customize the simulation, `{whatifbandit}` provides this core set of simulation features:

- **Two bandit algorithms**: probability matching via Thompson Sampling (TS), and
UCB1 (Upper Confidence Bound).
- **Flexible assignment periods**: individual-level, batch-based, or date-based
 assignment cutoffs for bandit updates.
- **A configurable lookback window and discounting**: complete control over the number of periods
used for bandit updates, and how to weight them for the computation.
- **Perfect and imperfect information**: you can simulate the real-time trials
 where outcomes are not known immediately (`delayed_feedback`), so the
 algorithm must make decisions with incomplete information, just as you would in the field.
- **Block-randomized and cluster-randomized designs**
- **Control augmentation and random assignment proportions**: these exploration parameters enforce a
minimum level of continued exploration for the algorithms, allowing for valid treatment effect estimation.
- **Joint hypothesis testing** (***Experimental***): bootstrap or randomization inference [@offer-westort2021] to test
 whether *any* treatment effect exists across arms.
- **pairwise contrasts**: pairwise contrasts can be automatically computed versus the control, versus the current best arm,
 or all pairwise comparisons among treatments.

The package is also built to scale effectively. When conducting multiple
replications of the same design, parallel processing can be used via
[future](https://future.futureverse.org/). [data.table](https://rdatatable.gitlab.io/data.table/) is
also supported for potentially large datasets.

I encourage a look at the documentation for `mab_from_rct()`, `simulate_mab()`, and `joint_test()`
for a more detailed look at all these features and their defaults.

# Resimulating an RCT as a Bandit Trial

## A Single Resimulation

To resimulate an RCT as a bandit trial, use `mab_from_rct()`. It requires a formula relating the
outcome to the treatment, with any potential blocking or clustering included in `block()` and
`cluster()` wrappers, and how the cutoff points are determined.


```{r}
set.seed(543645)
first_sim <- mab_from_rct(
  success ~ condition,
  data = tanf,
  algorithm = "ucb1",
  period_method = "batch",
  period_length = 1000,
  whole_experiment = TRUE
)
```

Here, `period_method = "batch"` with `period_length = 1000` tells the algorithm
to recompute assignment probabilities every 1,000 participants.
`whole_experiment = TRUE` means the probabilities used to impute new outcomes are estimated from
the entire dataset (instead of the data received up to that point).

## Blocking and Date-Based Periods with Delayed Feedback

The original experiment in @moore2022 randomized within blocks defined by the closest service
center. There is also a natural delay from when notification letters were assigned (and sent)
compared to when researchers received the recertification outcome. The example below accounts for these features to make the
experiment a more realistic resimulation.

```{r}
#| include: FALSE
load("tanf_sims.RData")
```

```{r}
#| eval: FALSE
future::plan("multisession", workers = 2)
set.seed(53254)
tanf_simulations <- mab_from_rct(
  success ~ condition + block(service_center),
  data = tanf,
  algorithm = "thompson",
  period_method = "date",
  time_unit = "month",
  period_length = 1,
  delayed_feedback = TRUE,
  assignment_date_col = letter_sent_date,
  success_date_col = date_of_recert,
  date_col = appt_date,
  month_col = recert_month,
  whole_experiment = FALSE,
  random_assign_prop = 0.3,
  r = 100,
  keep_data = TRUE,
  seed = TRUE
)
future::plan("sequential")
```

Some things to highlight are:

- `block(service_center)` inside the formula tells `{whatifbandit}` to treat
 `service_center` as a blocking variable, so treatment assignment and outcome
 imputation respect the original block structure.
- `period_method = "date"` together with `time_unit = "month"` and
 `period_length = 1` recomputes assignment probabilities once per calendar
month, rather than after a fixed number of participants. This mimics the assignment schedule of the
original experiment.
- `delayed_feedback = TRUE`, along with the `assignment_date_col`, `success_date_col`, `date_col`,
and `month_col` arguments, tells the algorithm which outcomes were actually observable by the time
each new assignment decision was made, capturing the time lag.
- `random_assign_prop = 0.3` enforces that 30% of assignment probability mass is
 held out for random (non-adaptive) assignment across arms in every period,
 guaranteeing a minimum level of continued exploration rather than letting the
 algorithm fully commit to the arm it currently believes is best.
- ` r = 100` requests 100 independent resimulations of the trial since any single resimulation is just one draw of
 assignments and imputations.
- `keep_data = TRUE` forces the `new_data` object (explained below) to be saved for each iteration. By default, they are not saved to limit the memory footprint of the final output.
- Running `future::plan("multisession")` before the call means the 100 resimulations are conducted in
parallel across multiple workers. `seed` is an argument passed to `future::future()` for parallelization, setting it to `TRUE`
ensures that, given the same seed in the previous `set.seed()`, the same results will be returned
(essentially normal R behavior for RNG with seeds). Users have complete control over the
parallelization, and can take advantage of `...` to pass arguments to `furrr::furrr_options()`
which control the parallel setup.

# Simulating a New Adaptive Trial from Scratch

When no trial data exists yet, `simulate_mab()` lets you specify assumed true
success probabilities for each arm and simulate how a bandit algorithm would
behave under those assumptions.

```{r}
p <- matrix(
  c(0.20, 0.35, 0.5, 0.27, 0.33),
  ncol = 1,
  dimnames = list(
    c("control", "T1", "T2", "T3", "T4"),
    NULL
  )
)
print(p)

set.seed(123)
sim_from_scratch <- simulate_mab(
  n = 2000,
  t = 20,
  p = p,
  algorithm = "thompson",
  random_assign_prop = 0.3,
  contrasts = "all"
)
```

In this example:

- `p` encodes the researcher's assumed true success probability for each arm.
- `n = 2000` sets the total number of simulated participants, and `t = 20` sets
 the number of periods over which assignment probabilities are updated.
- `contrasts = "all"` will compute all $\binom{5}{2}$ pairwise effect contrasts between the treatment arms.

# Understanding The Output

## Class Structure

Both `mab_from_rct()` and `simulage_mab()` return an S3 `.mab` class with the same fields.
Additionally, each object has a second class identifier, marking which function produced it and
whether multiple trials were run: `single_rct_mab`, `multi_rct_mab`, `single_param_mab`, and
`multi_param_mab`. Using the computed examples above:

```{r}
class(sim_from_scratch)
class(first_sim)
class(tanf_simulations)
```

This class structure at the moment has no implemented generics, but plans in the future are to write
generics which expedite common analyses on the final output objects.

Regardless of function or number of replications conducted, every `.mab` object contains the same
fields.

## Output Data

First, the `new_data` from the resulting MAB procedure is stored, allowing users to interact
with the results directly.^[By default when multiple simulations are conducted `new_data` is not
returned to limit the memory required for the output object, but can be requested by setting
`keep_data = TRUE`.]
```{r}
head(first_sim$new_data)
```

## MAB Results

Second, `bandit` contains the diagnostic information of the bandit procedure, which can be used to
examine the algorithm's behavior (`statistic`) and the resulting assignment_probabilities
(`assignment_prob`) over the course of the experiment. The final assignment counts to each treatment
are included as well (`assignment_quant`). Below is a glimpse of each object, along with an example
plot of the assignment probabilities over time for the @moore2022 MAB replication.

```{r}
head(sim_from_scratch$bandit$statistic)
head(sim_from_scratch$bandit$assignment_prob)
sim_from_scratch$bandit$assignment_quant

# Assignment probabilities over time
tanf_simulations$bandit$assignment_prob |>
  pivot_longer(
    cols = c(no_letter, open_appt, specific_appt),
    values_to = "p",
    names_to = "mab_condition"
  ) |>
  summarize(
    avg_p = mean(p),
    sd = sd(p),
    .by = c(mab_condition, period_number)
  ) |>
  mutate(
    lo = avg_p - 2 * sd,
    hi = avg_p + 2 * sd,
    mab_condition = str_to_title(str_replace(mab_condition, "_", " "))
  ) |>
  ggplot(aes(
    x = period_number,
    y = avg_p,
    color = mab_condition,
    fill = mab_condition
  )) +
  geom_ribbon(aes(ymax = hi, ymin = lo), alpha = 0.15, color = NA) +
  geom_line(linewidth = 0.8) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Period",
    y = "Probability of Assignment",
    title = "Assignment Probability to Treatment
Arms Over Time",
    subtitle = "With 2 Standard Deviation Uncertainty Ribbon",
    color = "Treatment Arm",
    fill = "Treatment Arm"
  ) +
  theme_minimal() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(face = "italic")
  )
```

## Conditional Expectations and Treatment Effects

Third, `means` and `contrasts` contain the final conditional expectation estimates for each treatment
arm (`means`) and treatment effect estimates (`contrasts`) for the requested configuration, along
with their standard errors.
Treatment effects can be requested via the `contrasts` argument for the following configurations:

- `"best"`: effect relative to the best treatment arm, as determined by the maximum UCB1 value or
TS probabilities.
- `"control"`: effect relative to the labelled control arm.
- `"both"`: both of the above
- `"all"`: all $\binom{n}{k}$ pairwise differences among arms

Estimates are reported for the AW-AIPW estimator [@hadad2021a], the IPW estimator
[@offer-westort2021], and a traditional OLS estimator. A quick overview of each is:

- AW-AWIPW is unbiased and asymptotically normal [@hadad2021a], so it can be used for hypothesis
tests with the standard normal distribution.
- IPW is unbiased with a non-normal asymptotic sampling distribution, but @offer-westort2021
that in extremely large samples t-tests may still provide proper coverage.
- OLS is biased and completely invalid to use for inference, but is provided for comparison
nonetheless.

Below is a glimpse of each object along with examples of a confidence interval plot for the AW-AIPW
mean estimates, and a p-value computation for the AW-AIPW treatment effects.

```{r}
head(sim_from_scratch$means)
# Plotting 95% CI using Normal distribution for AW-AIPW estimates
sim_from_scratch$means |>
  filter(estimator == "AW-AIPW") |>
  mutate(
    lo = mean + qnorm(0.025) * se,
    hi = mean + qnorm(0.975) * se,
    mab_condition = str_to_title(mab_condition),
    mab_condition = fct_reorder(mab_condition, mean)
  ) |>
  ggplot(aes(x = mean, y = mab_condition)) +
  geom_errorbar(
    aes(xmin = lo, xmax = hi),
    orientation = "y",
    width = 0.15,
    color = "grey40",
    linewidth = 0.6
  ) +
  geom_point(size = 2.5, color = "#2c3e50") +
  theme_minimal(base_size = 12) +
  labs(
    y = NULL,
    x = "Probability of Success",
    title = "Estimated Probability of Success By Arm",
    subtitle = "AW-AIPW Point Estimates with 95% Confidence Interval"
  ) +
  theme(
    panel.grid.major.y = element_blank(),
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(face = "italic")
  )
# Simple test of no treatment effect between each pair
sim_from_scratch$contrasts |>
  filter(estimator == "AW-AIPW") |>
  mutate(
    z_stat = est / se,
    p_value = pnorm(abs(z_stat), lower.tail = FALSE) * 2
  ) |>
  select(arm1, arm2, est, se, z_stat, p_value)
```

### Note on Causal Inference and Variance

Adaptive assignment probabilities do not strictly violate the assumptions required for causal
inference [@rubin1974; @rubin1977; @rubin1978; @rubin1990; @rubin2005; @holland1986]. Potential
outcomes remain well-defined, and more complex
assumptions like SUTVA derive from the design and administration, not the fact
that assignment probabilities are constant. Most importantly, however, is that treatment assignment is
unconfounded conditional on the history of previous assignments and outcomes, which is why the
simple IPW estimator is unbiased under an adaptive trial. Similar to how IPW uses a known or estimated
propensity score when treatment assignment is unconfounded only conditional on a set of
observed covariates to remove bias, the IPW in the adaptive case uses the exact probabilities of assignment
produced by the MAB algorithm to account for the bias. Additionally, just as
researchers model pre-treatment covariates in RCTs purely to improve efficiency, the AW-AIPW's
adaptive weighting serves the same role, controlling the variance compared to the IPW without
impacting unbiasedness. However, the adaptive weights play an even more crucial role than covariate
adjustments, since they create the conditions necessary to prove a central limit theorem
[@hadad2021a], allowing for valid hypothesis testing using the asymptotically normal sampling distribution.

The only assumption that you need to be careful about is positivity [@rubin1978; @holland1986], since a MAB
algorithm can drop a treatment arm permanently from the trial. This assumption is not only required
for a causal analysis but also for the asymptotic normality of the AW-AIPW to hold because of the
assumptions imposed on the adaptive weights required for the central limit theorem [@hadad2021a].
`{whatifbandit}` can ensure this is the case with`control_augment` and `random_assign_prop`, which
prevents assignment probabilities of 0 to the control arm or to all the arms jointly. If these
settings were unused and arms were dropped from the trial, then their estimated mean or treatment
effect would be unreliable, and any hypothesis test conducted would be invalid. However,
heuristically it could be inferred that the dropped treatment arm is worse than the rest, since it
is not deemed worthy of any of the assignment probability mass by the adaptive algorithm.

This result would be the intended behavior, and one of the potential benefits to running a MAB
experiment. In an experiment with many treatment arms, a traditional RCT, with equal assignment
probabilities, would produce inefficient estimates, since the share of participants assigned to the
subsets of best and worst arms is the same, conditional on the subset size. However, an
adaptive design, using an optimal MAB algorithm to generate assignment probabilities, would quickly
recognize which arms are worse and reallocate assignment mass to the subset of best treatment arms.
Therefore, the subset of best treatment arms, by having more observations, will have more precise
final estimates, while the subset of worst treatment arms, by being dropped, will have less precise
or invalid estimates, as opposed to an RCT where all arms are estimated with the same degree of
precision, creating an important variance reallocation trade-off that can be controlled with `random_assign_prop`
and `control_augment`.

An example of this, using the realistic resimulation of @moore2022,
is below:

```{r}
# Selecting a random trial from the 100 resimulations
set.seed(12345)
i <- sample(100, 1)
bandit <- tanf_simulations$means |> filter(trial == i & estimator == "AW-AIPW")
rct <- estimatr::lm_robust(
  success ~ condition - 1,
  data = tanf,
  se_type = "HC2"
)

bandit
tibble(
  treatment = str_replace(names(rct$coefficients), "condition", ""),
  coefs = rct$coefficients,
  se = rct$std.error
)

bandit_assignments <- tanf_simulations$new_data |>
  unnest(data) |>
  filter(
    trial == i &
      mab_condition == "specific_appt"
  ) |>
  nrow()
regular_assignments <- tanf |> filter(condition == "specific_appt") |> nrow()

coef_diff <- abs(
  bandit$mean[bandit$mab_condition == "specific_appt"] -
    coefficients(rct)["conditionspecific_appt"]
)
se_diff <- abs(
  bandit$se[bandit$mab_condition == "specific_appt"] -
    rct$std.error["conditionspecific_appt"]
)
```

As you can see above, the standard errors for the open-appointment treatment and the no-letter
control group are larger than the corresponding standard errors from the RCT. However, the specific
appointment letter features a standard error `r se_diff` lower than the RCT, and final estimate that
differs by `r coef_diff`, suggesting that the `r bandit_assignments - regular_assignments`
additional assignments both marginally tightened the estimate and potentially moved it closer to the true underlying
probability of success for the arm. In fact, even with increased uncertainty on the other arms'
estimates, this does not preclude us from detecting treatment effects with respect to the no-letter
group for both arms as the original experiment did [@moore2022], though this may not be true in all
100 of the resimulations. Additionally, the relative ranking of the treatments is the same as in the
original trial, preserving the conclusions in @moore2022 with the added benefits described above.

```{r}
bandit |>
  bind_cols(
    bandit |>
      filter(mab_condition == "no_letter") |>
      select(control = mean, control_se = se)
  ) |>
  mutate(
    z_se = sqrt(se^2 + control_se^2),
    z_stat = (mean - control) / z_se,
    p_value = pnorm(abs(z_stat), lower.tail = FALSE) * 2
  ) |>
  filter(mab_condition != "no_letter") |>
  select(mab_condition, z = z_stat, se = z_se, p_value)
```

## Other Returned Objects

The other objects returned are the F-statistics from the IPW and OLS regressions (`f-stats`), the
`lm_robust` objects from the IPW and OLS regressions (`models`),^[Only provided if clustering is
used so that arbitrary contrasts can be estimated later with the appropriate CR2 covariance matrix
estimator using `clubSandwich::linear_contrast()`. By default, with multiple simulations, these are
not returned since each is large memory-wise but can be requested by setting `keep_models = TRUE`.]
and a `config` object containing the arguments and original `call` object for easy replication since
there are so many customizable settings.

# Joint Hypothesis Tests

`{whatifbandit}` also implements MAB-aware joint hypothesis tests either via a randomization inference procedure in
@offer-westort2021 or a bootstrap test, testing whether any treatment effect exists across arms.
Both tests are experimental, unverified, and noted to have low power, but are included for
completeness nonetheless. Both tests use the F-statistic from the IPW regression as its test
statistic, and simulate its sampling distribution under the appropriate $H_0$ to compute a p-value.

```{r}
# joint_test parallelizes through the same options as the original sim
sim_from_scratch$config$parallel <- furrr::furrr_options(seed = TRUE)
set.seed(123)
boot <- joint_test(sim_from_scratch, method = "bootstrap", r = 100)
rand <- joint_test(sim_from_scratch, method = "randomization", r = 100)
```

Both return a list with the observed F-statistic (`f-stat`), the null distribution generated
under $H_0$ (`null_distribution`), the resulting `p_value` (The proportion of null F-statistics
greater than the observed one), the `method` used, the number of simulations to perform the test (`
r`), and the number of those simulations with a non-`NA` F-statistic (`effective_r`; Some iterations
can produce extremely skewed treatment assignment, resulting in an undefined F-statistic)

```{r}
boot$f_stat
boot$p_value
boot$method
```


# Estimating Regret and Welfare Benefits

Up to this point, only estimator efficiency has been discussed as a benefit using response-adaptive
trials. Everything else discussed: means, treatment effects, and joint tests can all be computed
with data from a valid RCT, in a much more straightforward fashion. What these inferential tools
don't tell us is whether adaptivity has
produced any tangible benefits that would be worth incurring the additional
complexity costs of adaptive trials in both administration and estimation. This section focuses on
tangible gains, which could be used to persuade a stakeholder or policymaker to approve of an
adaptive experiment as opposed to a statistician. To explore these
potential benefits and drawbacks, I'll walk through a regret analysis of an adaptive trial with the
`tanf` dataset from @moore2022.

## Visualizing Regret

Regret is a commonly used metric in the MAB literature [@lai1985; @agrawal2012; @agrawal2017;
@auer2002],^[Not an exhaustive list of the MAB or the regret literature, simply a small sample of
the papers referring to `{whatifbandit}`'s implemented algorithms.] to quantify and compare decision algorithms. For a single observation in the
trial, regret is simply the difference between the outcome under the optimal treatment and the
outcome under the selected treatment. Thus it quantifies how much you **regret** your selected
treatment, and when you select the optimal treatment the regret is **0**. Similar to the fundamental
problem of causal inference, both of these outcomes for each observation cannot be directly observed,
so expected regret is used instead. The computation is the same, but outcomes are just replaced with
the means for each treatment arm. Added over time cumulatively, regret characterizes the asymptotic
performance of MAB algorithms. Algorithms are traditionally assessed based on the growth rate of
their cumulative expected regret, which for optimal algorithms is bounded below by
$\mathcal{O}(\ln(T))$, as proven by @lai1985.^[This result is extremely important for MAB
algorithms, as it shows optimal algorithms, even on infinite time horizons, must always continue to sample
from sub-optimal arms. Regret cannot be eliminated, nor can the growth stop since $\lim_{T \to
\infty} \ln(T)$ diverges to $\infty$. Therefore, if an algorithm were to eliminate regret,
it would only be doing so trivially (such as only selecting one arm the whole trial, which just happens
to be the best) since the theorem would bind an optimal algorithm.]

Although `{whatifbandit}` does not have a dedicated function to compute cumulative expected regret,^[Feature planned for a future release.] it can easily be done using the available data, so long as the
outcomes are stored (remember to set `keep_data = TRUE` for multiple
simulations). To illustrate regret clearly, I use a pure bandit trial, with no enforced exploration,
under sequential assignment (1-observation batches) to match the theoretical MAB
problem.

```{r}
#| include: FALSE
load("pure_sims.RData")
```

```{r}
#| eval: FALSE
set.seed(53245)
pure_ucb_sims <- mab_from_rct(
  success ~ condition,
  data = tanf,
  algorithm = "ucb1",
  period_method = "individual",
  keep_data = TRUE,
  seed = TRUE,
  r = 100,
  whole_experiment = TRUE
)
pure_ts_sims <- update(pure_ucb_sims$config, algorithm = "thompson")
```

It is important to note that under sequential assignment, the ordering of our data can significantly
impact the results, because early observations (under no enforced exploration) set the initial
trajectory of the future assignments. While this is true for all adaptive trials, those with larger
period sizes are less susceptible since the algorithm receives more representative and stable
within-period estimates of each treatment arm before and during the adaptive procedure.

The `tanf` data is inherently ordered by the corresponding appointment dates, even though the
original trial assigned by month. However, for an RCT without a natural ordering of observations,
such as one where all the data was collected in a single period without any time delays, it may be
beneficial to randomly permute the data in each iteration of the trial as a robustness
check.^[Feature planned for a future release.]

Regret is computed with respect to the original RCT mean estimates, since they are perfectly valid
in their own right, and agnostic of the adaptive environment. However, the AW-AIPW and IPW estimators
are unbiased for the true mean, so the corresponding estimates could be used as well. For an original
MAB trial, the posited population parameters should be used instead, since they are the truth.

Regret is commonly visualized in a plot, like the one below, where the relative rates of cumulative
growth can be examined to evaluate each adaptive algorithm. As a benchmark, the same cumulative
regret measure is also computed under the original RCT assignment scheme under the same
ordering.^[Even though the batch-size of the original RCT is changing, a resimulation is not
required. This is because assignments are not dependent on previous observations, so treatment
assignments can be made all at once before the trial, no matter how they will be administered over
time]. For this non-adaptive benchmark, regret will accrue linearly, contrasting with the
expected
sub-linear regret of the bandit algorithms.

```{r}
# RCT Mean Estimates
tanf_estimates <- coef(lm(success ~ condition - 1, tanf))
names(tanf_estimates) <- str_replace(names(tanf_estimates), "condition", "")
best <- names(tanf_estimates)[which.max(tanf_estimates)]

# Original RCT Data. Adding new columns to mimic output of MAB experimental data.
rct_regret <- tanf |>
  mutate(period_number = row_number(), mab_condition = condition)

# Selecting a random trial of the 100 previous simulations for regret analysis
set.seed(0934)
i <- sample.int(100, 1)

bandit_regrets <- lapply(list(pure_ucb_sims, pure_ts_sims), \(mab) {
  mab$new_data |> unnest(data) |> filter(trial == i)
})

lapply(
  list(
    "Static" = rct_regret,
    "UCB1" = bandit_regrets[[1]],
    "Thompson Sampling" = bandit_regrets[[2]]
  ),
  \(df) {
    df |>
      mutate(
        regret = tanf_estimates[[best]] - tanf_estimates[mab_condition]
      ) |>
      arrange(period_number) |>
      mutate(c_regret = cumsum(regret))
  }
) |>
  bind_rows(.id = "Algorithm") |>
  select(Algorithm, c_regret) |>
  mutate(x = row_number(), .by = Algorithm) |>
  ggplot(aes(x = x, y = c_regret, color = Algorithm)) +
  geom_line(linewidth = 0.8) +
  theme_minimal(base_size = 12) +
  labs(
    y = "Cumulative Regret",
    x = "Participants Assigned",
    title = "Cumulative Expected Regret By Adaptive Algorithm",
    subtitle = paste0(
      "Static Refers To Original RCT; Simulation ",
      i,
      " Selected"
    )
  ) +
  theme(
    panel.grid.major.y = element_blank(),
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(face = "italic")
  )
```

As expected, the original RCT has linear cumulative regret, while both TS and UCB1 exhibit
sub-linear regret growth patterns, showing they make more optimal assignments, in terms of
regret, than the block-randomized design used by @moore2022.

## Quantifying Welfare

Regret and welfare are two sides of the same coin. While a regret analysis asks "*What is the
expected loss from a sub-optimal decision?*", a welfare analysis asks "*What is the benefit from an
additional optimal decision?*" For a traditional MAB analysis, the analog
to regret is reward, and the cumulative expected reward can be calculated, where each reward
is the expected outcome of the treatment arm selected. However, for experiments in the social
sciences, this simple metric is uninformative.

Firstly, for experiments where the outcome is binary,
like @moore2022, the cumulative reward is simply the sum of the estimated probabilities of success
corresponding to each assigned outcome. While this could be used to compare between algorithms,
since a higher cumulative reward corresponds to assigning more participants to the better treatments,
there is no way to assess the magnitude of the gain. The real-world payoff of assigning more people
to the most effective treatments, for @moore2022, is the monetary value of the TANF benefits that
would not have been retained under the traditional RCT because of the sub-optimal treatment
assignments. In fact, to go even further, it would not just be the monetary value, but how the
additional benefits improve the standard of living for recipients. This, however, is a much less defined
measure, requiring a more normative analysis of the TANF program and whether more money received is
better, going far beyond the scope of this vignette.

Secondly, even for experiments with
continuous outcomes, the measured outcome of interest is not always what researchers truly care
about. Take, for example, an educational experiment where the outcome variable of the study is
standardized test scores. Here the expected cumulative reward has magnitude, as you can clearly
identify the additional points scored due to the assignment algorithm. However, researchers may only
care about using test scores to proxy wage differentials or college acceptances, so once again the
outcome variable must be translated into a downstream quantity that can be used to accurately
gauge the additional value the adaptive trial is creating due to more optimal assignments.

I refer to the results of these outcome translations as **welfare**, an experiment-context-specific
measure of how the experiment's outcome variable maps onto a concrete gain from receiving a more
effective treatment, which can be used to identify the tangible benefits of adaptive trials against
traditional RCTs. The complexity of **welfare** depends on the experimental context, and certainly
the same experiment can have different **welfare** measures depending on the benefit of interest, and the
level of plausibility desired by the researcher.

As explained above, in the @moore2022 experiment, the estimated **welfare** comes from how the
additional TANF recertifications map to the additional monetary value of TANF benefits which would
not have been received under the initial trial. For simplicity's sake, even though the experiment
was conducted in 2017 [@moore2022], assume it is conducted today, and subject to the current TANF
benefit structure. Below are the monthly TANF benefits in Washington, D.C. effective from October 1,
2025 to September 30, 2026, for each listed household size:

```{r}
tanf_benefits <- data.frame(
  hh_size = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
  benefit = c(504, 629, 803, 983, 1136, 1335, 1532, 1691, 1863, 2024)
)
```

I assume that all observations in the data have the same household size, and
to avoid picking a single household size, I present the calculations for all 10, providing a clear
lower and upper bound on the potential welfare. For a more realistic point estimate, a weighted
average of the benefits, based on the distribution of household sizes for TANF recipients in
Washington, D.C. could be used.

Once again `{whatifbandit}` does not have a dedicated function for computing welfare, because of the
experiment-specific computations required.^[Though it is planned for a future release.]
Below is the welfare computed using the same UCB1 and
TS pure bandit trials from the regret section above, using the average number of additional
recertifications compared to the RCT across the 100 trials multiplied by each benefit amount. Since
I kept the output data (`keep_data = TRUE`) from each of the 100 trials, I can use the observed number of successful
recertifications from each MAB trial to compute the average number of additional successes when
compared to the original RCT.

```{r}
true_success_n <- sum(tanf$success)
data_method <- lapply(
  list("Thompson Sampling" = pure_ts_sims, "UCB1" = pure_ucb_sims),
  \(mab) {
    mab$new_data |>
      unnest(data) |>
      summarize(
        add_success_n = sum(mab_success) - true_success_n,
        .by = trial
      ) |>
      summarize(mean_add_success_n = mean(add_success_n))
  }
) |>
  bind_rows(.id = "algorithm") |>
  pivot_wider(names_from = algorithm, values_from = mean_add_success_n) |>
  cbind(tanf_benefits) |>
  mutate(
    Welfare_TS = `Thompson Sampling` * benefit,
    Welfare_UCB1 = UCB1 * benefit
  ) |>
  select(hh_size, benefit, Welfare_TS, Welfare_UCB1)

data_method
lo <- min(data_method[1, c("Welfare_TS", "Welfare_UCB1")])
hi <- max(data_method[10, c("Welfare_TS", "Welfare_UCB1")])
```

In the case that `keep_data = FALSE`, the outcomes from each
simulated trial are unknown, so the above approach cannot be used directly. Instead, another method
is required, relying on what is available in the standard output. The assignment quantities per
treatment per trial are always reported, so using the probabilities of successful recertification
from the original RCT, the expected number of successes for each trial can be estimated. From here,
the computation is the same: calculate the average additional number of successes compared to
the RCT, using the expected successful recertifications per trial, then multiply by the benefit
amount.^[In the future, the
number of successes may be reported in standard output to avoid this issue.]

```{r}
expectation_method <- lapply(
  list("Thompson Sampling" = pure_ts_sims, "UCB1" = pure_ucb_sims),
  \(mab) {
    quants <- mab$bandit$assignment_quant[, -1] |> as.matrix()
    ests <- tanf_estimates[colnames(quants)]
    return(mean(quants %*% ests - true_success_n))
  }
) |>
  bind_rows(.id = "algorithm") |>
  cbind(tanf_benefits) |>
  mutate(
    Welfare_TS = `Thompson Sampling` * benefit,
    Welfare_UCB1 = UCB1 * benefit
  ) |>
  select(hh_size, benefit, Welfare_TS, Welfare_UCB1)
expectation_method
```

Both methods produce similar welfare estimates. In general, the welfare is positive,
indicating that for the @moore2022 experiment, a response-adaptive trial could have produced anywhere
between \$`r format(lo, big.mark = ",")` and \$`r format(hi, big.mark = ",")` per month in TANF benefits from the additional recertifications in
the experiment.^[The lower bound is the minimum of the welfare under a household size of 1 (UCB1). The upper bound is the
maximum welfare under a household size of 10 (TS). Values are from the `keep_data = TRUE` method.
This range is only representative of an adaptive trial administered under the same settings as
the simulation used (pure bandit, no enforced exploration) and would likely be reduced in a more
realistic simulation where exploration is enforced to ensure valid causal inference.] For how
many months to extend the value for is uncertain and recipient dependent, since nothing prevents a
TANF recipient that failed to
recertify from reapplying to the program later, since they are still eligible under the income limits.^[This program churn is the exact problem @moore2022
attempts to tackle by creating better notification letters.]


# Concluding Remarks

The power of adaptive experiments, with enforced exploration, lies in how they straddle the line between
pure bandits and RCTs. Compared to a traditional RCT,
response adaptive trials produce optimal assignment patterns, resulting in additional real-world
gains from running the experiment. Although it is presumed the most effective treatment identified
by an RCT would be implemented permanently afterwards, the value of more optimal assignments
*during* the trial is non-negligible. This is especially the case when the experiment has a large
sample size, or takes place over a long period of time, since a longer trial delays when the most
effective treatment can be permanently rolled out. It is also the case for experiments involving
real, unknowing participants, since the effect of receiving the better treatment can meaningfully
improve their lives, as in @moore2022.

In the @moore2022 experiment, all the participants were D.C. TANF
recipients up for recertification during the timeframe of the study; they did not opt into the trial
and had no way of knowing they were participating until a letter was received. Receiving
one of the treatment letters, however, would significantly increase a recipient's chance of
recertifying, thus keeping their TANF benefits without any lapse. Understanding this, an optimal
method of treatment assignment should assign as many participants as possible to one of the
treatment letters or better yet the best treatment letter, to the extent that treatment effects
could still be reliably detected by hypothesis testing. This is precisely what occurred in the adaptive resimulations in this
vignette: reallocating participants to the better-performing treatment letters did not prevent us from
reaching the same conclusion as @moore2022 did, while also receiving the welfare gains associated
from the additional recertifications over the course of the trial.

Response adaptive designs, and their associated welfare gains that come from more optimal
assignment, on their own, do not violate the assumptions required for valid causal inference.
Treatment effects remain estimable, but simply require more complicated estimators than
traditional RCTs. What adaptive designs do require, though, is a clear understanding of the
precision tradeoffs they impose, since these tradeoffs can still impact whether an estimated effect
is determined to be statistically significant. `{whatifbandit}` is designed to help researchers
navigate this tradeoff, making it possible to explore, resimulate, and evaluate adaptive designs
before committing to them in the field.

# Future Plans
Plans for this package are ambitious and expansive:

- Support for Contextual Bandits (covariate-aware adaptation).
- Continuous outcomes.
- Post-simulation helper functions to accelerate analysis.
- Fixed confidence settings and early stopping for fixed-budget settings.
- Adaptive treatment pool, as in switching treatments in and out of the potential assignment pool
over time. This mimics an infinitely running experiment where new treatments are added in real time.
- Exporting modular internal functions for more fine-tuned user control, and to help others run
their own adaptive experiments in the field.

# Getting Help and Contributions

For more complete details on individual function arguments and behavior, consult
the full package documentation (`?mab_from_rct`, `?simulate_mab`, `?joint_test`).

If you have specific questions about the package, feel free to reach out to me by
email at <no9857a@american.edu>. If you encounter a bug, please open an issue on
[GitHub](https://github.com/Noch05/whatifbandit/issues) with a reproducible
example.

# References

:::{#refs}
:::