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.

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 (Moore et al. 2022). 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) (Moore et al. 2022). The experiment was conducted with real TANF recipients in Washington, D.C.

data("tanf")
head(tanf)
#> # A tibble: 6 × 21
#>   ic_case_id   service_center appt_date  condition recert_month letter_sent_date
#>   <chr>        <chr>          <date>     <fct>     <fct>        <date>          
#> 1 42c43695d7d… Anacostia      2017-06-26 no_letter July         2017-06-19      
#> 2 faff78c811a… Anacostia      2017-06-26 no_letter July         2017-06-19      
#> 3 31ba0f81073… Congress_Heig… 2017-06-26 no_letter July         2017-06-19      
#> 4 dc669332bb9… Congress_Heig… 2017-06-26 no_letter July         2017-06-19      
#> 5 a8c8de78fff… Congress_Heig… 2017-06-26 no_letter July         2017-06-19      
#> 6 e4c12ca3447… Congress_Heig… 2017-06-26 no_letter July         2017-06-19      
#> # ℹ 15 more variables: recert_id <dbl>, return_to_sender <chr>,
#> #   pdc_status <chr>, renewal_date <date>, notice_date.x <date>,
#> #   days_betwn_notice_and_recert_due <dbl>, cert_period_start <date>,
#> #   cert_period_end <date>, recert_status <chr>, denial_reason <chr>,
#> #   recert_month_year <chr>, notice_date.y <date>, recert_status_dcas <chr>,
#> #   date_of_recert <date>, success <dbl>
glimpse(tanf)
#> Rows: 3,517
#> Columns: 21
#> $ ic_case_id                       <chr> "42c43695d7d06af7f66fd072c3ee2593", "…
#> $ service_center                   <chr> "Anacostia", "Anacostia", "Congress_H…
#> $ appt_date                        <date> 2017-06-26, 2017-06-26, 2017-06-26, …
#> $ condition                        <fct> no_letter, no_letter, no_letter, no_l…
#> $ recert_month                     <fct> July, July, July, July, July, July, J…
#> $ letter_sent_date                 <date> 2017-06-19, 2017-06-19, 2017-06-19, …
#> $ recert_id                        <dbl> 151516, 153505, 151887, 155288, 15033…
#> $ return_to_sender                 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
#> $ pdc_status                       <chr> "Closed", "Closed", "Pending closure"…
#> $ renewal_date                     <date> 2017-07-31, 2017-07-31, 2017-07-31, …
#> $ notice_date.x                    <date> 2017-05-29, 2017-05-29, 2017-05-29, …
#> $ days_betwn_notice_and_recert_due <dbl> 63, 63, 63, 63, 63, 63, 63, 63, 63, 6…
#> $ cert_period_start                <date> 2016-11-01, 2017-03-07, 2016-11-01, …
#> $ cert_period_end                  <date> 2017-07-31, 2017-07-31, 2017-07-31, …
#> $ recert_status                    <chr> "Notice Sent", "Notice Sent", "Notice…
#> $ denial_reason                    <chr> "No response from client", "No respon…
#> $ recert_month_year                <chr> "July 2017", "July 2017", "July 2017"…
#> $ notice_date.y                    <date> 2017-05-29, 2017-05-29, 2017-05-29, …
#> $ recert_status_dcas               <chr> "Denied due to failure of action", "D…
#> $ date_of_recert                   <date> NA, NA, NA, 2017-06-27, NA, NA, NA, …
#> $ success                          <dbl> 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0…

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-Westort, Coppock, and Green 2021) 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. 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.

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 Moore et al. (2022) 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.

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.

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)
#>         [,1]
#> control 0.20
#> T1      0.35
#> T2      0.50
#> T3      0.27
#> T4      0.33

set.seed(123)
sim_from_scratch <- simulate_mab(
  n = 2000,
  t = 20,
  p = p,
  algorithm = "thompson",
  random_assign_prop = 0.3,
  contrasts = "all"
)
#> Registered S3 method overwritten by 'clubSandwich':
#>   method    from    
#>   bread.mlm sandwich

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:

class(sim_from_scratch)
#> [1] "single_param_mab" "single_mab"       ".mab"             "list"
class(first_sim)
#> [1] "single_rct_mab" "single_mab"     ".mab"           "list"
class(tanf_simulations)
#> [1] "multi_rct_mab" "multi_mab"     ".mab"          "list"

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

head(first_sim$new_data)
#> # A tibble: 6 × 10
#>   condition success period_number mab_condition mab_success impute_req
#>   <fct>       <dbl>         <int> <chr>               <dbl>      <dbl>
#> 1 no_letter       0             1 no_letter               0          0
#> 2 no_letter       0             1 no_letter               0          0
#> 3 no_letter       0             1 no_letter               0          0
#> 4 no_letter       1             1 no_letter               1          0
#> 5 no_letter       0             1 no_letter               0          0
#> 6 no_letter       0             1 no_letter               0          0
#> # ℹ 4 more variables: impute_block <chr>, treatment_block <chr>,
#> #   mab_assign_prob <dbl>, ipw_weights <dbl>

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 Moore et al. (2022) MAB replication.

head(sim_from_scratch$bandit$statistic)
#> # A tibble: 6 × 6
#>        control       t1    t2       t3      t4 period_number
#>          <dbl>    <dbl> <dbl>    <dbl>   <dbl>         <int>
#> 1 0.000848     0.00338  0.958 0.0105   0.0272              1
#> 2 0.000647     0.000647 0.846 0.0173   0.135               2
#> 3 0.0000195    0.00496  0.888 0.00700  0.100               3
#> 4 0.00000804   0.00428  0.973 0.00216  0.0207              4
#> 5 0.000000140  0.00985  0.987 0.000856 0.00227             5
#> 6 0.0000000467 0.00109  0.995 0.000289 0.00389             6
head(sim_from_scratch$bandit$assignment_prob)
#> # A tibble: 6 × 6
#>   control     t1    t2     t3     t4 period_number
#>     <dbl>  <dbl> <dbl>  <dbl>  <dbl>         <int>
#> 1  0.2    0.2    0.2   0.2    0.2                1
#> 2  0.0606 0.0624 0.731 0.0674 0.0791             2
#> 3  0.0605 0.0605 0.652 0.0721 0.155              3
#> 4  0.0600 0.0635 0.681 0.0649 0.130              4
#> 5  0.0600 0.0630 0.741 0.0615 0.0745             5
#> 6  0.0600 0.0669 0.751 0.0606 0.0616             6
sim_from_scratch$bandit$assignment_quant
#> control      t1      t2      t3      t4 
#>     134     136    1440     136     154

# 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 (Hadad et al. 2021), the IPW estimator (Offer-Westort, Coppock, and Green 2021), and a traditional OLS estimator. A quick overview of each is:

• AW-AWIPW is unbiased and asymptotically normal (Hadad et al. 2021), 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-Westort, Coppock, and Green (2021) 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.

head(sim_from_scratch$means)
#> # A tibble: 6 × 5
#>    mean     se mab_condition estimator    df
#>   <dbl>  <dbl> <chr>         <chr>     <dbl>
#> 1 0.169 0.0334 control       AW-AIPW      NA
#> 2 0.360 0.0431 t1            AW-AIPW      NA
#> 3 0.496 0.0136 t2            AW-AIPW      NA
#> 4 0.285 0.0396 t3            AW-AIPW      NA
#> 5 0.342 0.0395 t4            AW-AIPW      NA
#> 6 0.172 0.0342 control       IPW        1995
# 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)
#> # A tibble: 10 × 6
#>    arm1  arm2        est     se z_stat  p_value
#>    <chr> <chr>     <dbl>  <dbl>  <dbl>    <dbl>
#>  1 t1    control  0.191  0.0545  3.50  4.58e- 4
#>  2 t2    control  0.327  0.0361  9.07  1.14e-19
#>  3 t3    control  0.116  0.0518  2.24  2.51e- 2
#>  4 t4    control  0.173  0.0517  3.35  8.06e- 4
#>  5 t2    t1       0.136  0.0452  3.01  2.62e- 3
#>  6 t3    t1      -0.0751 0.0585 -1.28  1.99e- 1
#>  7 t4    t1      -0.0179 0.0585 -0.306 7.59e- 1
#>  8 t3    t2      -0.211  0.0419 -5.04  4.56e- 7
#>  9 t4    t2      -0.154  0.0417 -3.69  2.25e- 4
#> 10 t4    t3       0.0572 0.0559  1.02  3.06e- 1

Note on Causal Inference and Variance

Adaptive assignment probabilities do not strictly violate the assumptions required for causal inference (Donald B. Rubin 1974, 1977, 1978, 1990; Donald B. Rubin 2005; Holland 1986). 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 (Hadad et al. 2021), allowing for valid hypothesis testing using the asymptotically normal sampling distribution.

The only assumption that you need to be careful about is positivity (Donald B. Rubin 1978; Holland 1986), 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 (Hadad et al. 2021). {whatifbandit} can ensure this is the case withcontrol_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 Moore et al. (2022), is below:

# 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
#> # A tibble: 3 × 6
#>   trial  mean     se mab_condition estimator    df
#>   <dbl> <dbl>  <dbl> <chr>         <chr>     <dbl>
#> 1    14 0.375 0.0183 no_letter     AW-AIPW      NA
#> 2    14 0.443 0.0154 open_appt     AW-AIPW      NA
#> 3    14 0.424 0.0145 specific_appt AW-AIPW      NA
tibble(
  treatment = str_replace(names(rct$coefficients), "condition", ""),
  coefs = rct$coefficients,
  se = rct$std.error
)
#> # A tibble: 3 × 3
#>   treatment     coefs     se
#>   <chr>         <dbl>  <dbl>
#> 1 no_letter     0.401 0.0143
#> 2 open_appt     0.464 0.0145
#> 3 specific_appt 0.447 0.0145

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 0.0000336 lower than the RCT, and final estimate that differs by 0.0238212, suggesting that the 218 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 (Moore et al. 2022), 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 Moore et al. (2022) with the added benefits described above.

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)
#> # A tibble: 2 × 4
#>   mab_condition     z     se p_value
#>   <chr>         <dbl>  <dbl>   <dbl>
#> 1 open_appt      2.83 0.0240 0.00472
#> 2 specific_appt  2.08 0.0234 0.0374

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),² 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-Westort, Coppock, and Green (2021) 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₀ to compute a p-value.

# 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₀ (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)

boot$f_stat
#> [1] 303.482
boot$p_value
#> [1] 0.08
boot$method
#> [1] "bootstrap"

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 Moore et al. (2022).

Visualizing Regret

Regret is a commonly used metric in the MAB literature (Lai and Robbins 1985; Agrawal and Goyal 2012, 2017; Auer, Cesa-Bianchi, and Fischer 2002),³ 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 𝒪(ln (T)), as proven by Lai and Robbins (1985).⁴

Although {whatifbandit} does not have a dedicated function to compute cumulative expected regret,⁵ 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.

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

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.⁷. For this non-adaptive benchmark, regret will accrue linearly, contrasting with the expected sub-linear regret of the bandit algorithms.

# 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 Moore et al. (2022).

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 Moore et al. (2022), 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 Moore et al. (2022), 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 Moore et al. (2022) 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 (Moore et al. 2022), 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:

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

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
#>    hh_size benefit Welfare_TS Welfare_UCB1
#> 1        1     504   30950.64     15719.76
#> 2        2     629   38626.89     19618.51
#> 3        3     803   49312.23     25045.57
#> 4        4     983   60366.03     30659.77
#> 5        5    1136   69761.76     35431.84
#> 6        6    1335   81982.35     41638.65
#> 7        7    1532   94080.12     47783.08
#> 8        8    1691  103844.31     52742.29
#> 9        9    1863  114406.83     58106.97
#> 10      10    2024  124293.84     63128.56
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.⁹

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
#>    hh_size benefit Welfare_TS Welfare_UCB1
#> 1        1     504   30399.21     18954.09
#> 2        2     629   37938.70     23655.01
#> 3        3     803   48433.66     30198.68
#> 4        4     983   59290.52     36968.00
#> 5        5    1136   68518.85     42721.92
#> 6        6    1335   80521.72     50205.78
#> 7        7    1532   92403.95     57614.43
#> 8        8    1691  101994.17     63593.99
#> 9        9    1863  112368.51     70062.45
#> 10      10    2024  122079.36     76117.23

Both methods produce similar welfare estimates. In general, the welfare is positive, indicating that for the Moore et al. (2022) experiment, a response-adaptive trial could have produced anywhere between $15,719.76 and $124,293.8 per month in TANF benefits from the additional recertifications in the experiment.¹⁰ 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.¹¹

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 Moore et al. (2022).

In the Moore et al. (2022) 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 Moore et al. (2022) 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 with a reproducible example.

References