Provides a unified framework for designing, simulating, and analyzing implementation rollout trials, including stepped wedge, sequential rollout, head-to-head, multi-condition, and rollout implementation optimization designs. The package enables users to flexibly specify rollout schedules, incorporate site-level and nested data structures, generate outcomes under rich hierarchical models, and evaluate analytic strategies through simulation-based power analysis. By separating data generation from model fitting, the tools support assessment of bias, Type I error, and robustness to model misspecification. The workflow integrates with standard mixed-effects modeling approaches and the tidyverse ecosystem, offering transparent and reproducible tools for implementation scientists and applied statisticians.
Rollout trials, including stepped wedge and sequential rollout designs, allow staggered implementation of interventions across sites while collecting outcome data over time. These designs are increasingly used in implementation science because they balance practical, ethical, and statistical considerations when rolling out interventions in real-world settings.
The rollout R package supports the design and analysis of rollout
trials using simulation-based methods. This tutorial demonstrates how to
simulate, analyze, and estimate power for a stepped wedge trial with
multiple sites per cohort. The workflow aligns with the following
structure:
By following this pipeline, you can efficiently evaluate and optimize your rollout trial design before conducting your study.
rollout PackageYou can install the release version directly from CRAN:
install.packages("rollout")
Or, you can install the development version of rollout directly from
GitHub using the remotes package:
# Install remotes package if needed
install.packages("remotes")
# Install rollout from GitHub
remotes::install_github("iancero/rollout")
We begin by loading the required packages and setting a seed for
reproducibility. Throughout this tutorial, we will use the tidyverse
for data manipulation and plotting, and the rollout package for
design, simulation, and evaluation of rollout trials.
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.5.2
#> Warning: package 'ggplot2' was built under R version 4.5.2
#> Warning: package 'purrr' was built under R version 4.5.2
#> Warning: package 'stringr' was built under R version 4.5.2
library(rollout)
set.seed(1234) # for reproducibility
We will simulate a stepped wedge rollout trial with:
"intv") at a different
time while data are collected continuously.We use a tribble() for clarity, then convert to long format using
pivot_schedule_longer().
# Create a stepped wedge schedule with 8 sites, 4 cohorts, 8 timepoints
rollout_schedule <- tribble(
~cohort, ~site, ~t1, ~t2, ~t3, ~t4, ~t5, ~t6, ~t7, ~t8,
1, "A", "ctrl", "ctrl", "intv", "intv", "intv", "intv", "intv", "intv",
1, "B", "ctrl", "ctrl", "intv", "intv", "intv", "intv", "intv", "intv",
2, "C", "ctrl", "ctrl", "ctrl", "intv", "intv", "intv", "intv", "intv",
2, "D", "ctrl", "ctrl", "ctrl", "intv", "intv", "intv", "intv", "intv",
3, "E", "ctrl", "ctrl", "ctrl", "ctrl", "intv", "intv", "intv", "intv",
3, "F", "ctrl", "ctrl", "ctrl", "ctrl", "intv", "intv", "intv", "intv",
4, "G", "ctrl", "ctrl", "ctrl", "ctrl", "ctrl", "ctrl", "intv", "intv",
4, "H", "ctrl", "ctrl", "ctrl", "ctrl", "ctrl", "ctrl", "intv", "intv"
)
# Transform to long format for simulation
long_schedule <- rollout_schedule |>
pivot_schedule_longer(time_cols = starts_with("t"))
# Preview structure
long_schedule |>
head()
#> # A tibble: 6 × 5
#> cohort site chron_time condition local_time
#> <dbl> <chr> <dbl> <fct> <int>
#> 1 1 A 1 ctrl 0
#> 2 1 A 2 ctrl 1
#> 3 1 A 3 intv 0
#> 4 1 A 4 intv 1
#> 5 1 A 5 intv 2
#> 6 1 A 6 intv 3
Next, we specify the unit-level information for each site, which is sometimes available in advance of a trial. For this example, each site will have a slightly different number of units, reflecting real-world differences in site sizes.
We will then join this information to our long-format rollout schedule
using join_info(). This will expand the dataset to include a row for
each unit at each timepoint.
# Create unit-level information with variable site sizes
unit_info <- tribble(
~site, ~n_units,
"A", 18,
"B", 22,
"C", 20,
"D", 19,
"E", 21,
"F", 20,
"G", 23,
"H", 17
)
# Join unit-level info to the long-format schedule
design_df <- join_info(
long_schedule,
unit_info = unit_info,
by = "site",
uncount_vars = "n_units",
.ids = "unit_id"
)
# Preview the expanded design
design_df |>
head()
#> # A tibble: 6 × 6
#> cohort site chron_time condition local_time unit_id
#> <dbl> <chr> <dbl> <fct> <int> <int>
#> 1 1 A 1 ctrl 0 1
#> 2 1 A 1 ctrl 0 2
#> 3 1 A 1 ctrl 0 3
#> 4 1 A 1 ctrl 0 4
#> 5 1 A 1 ctrl 0 5
#> 6 1 A 1 ctrl 0 6
We will now add simulation parameters to our design. A key advantage
of rollout is the ability to expand a parameter grid, enabling the
evaluation of multiple scenarios within the same simulation pipeline.
In this example: - We will evaluate two intervention effect sizes
(b_intv = 0.2 for a small effect, b_intv = 0.5 for a moderate
effect). - We will specify a constant site-level standard deviation
(sigma_site = 0.5) and unit-level residual standard deviation
(sigma_unit = 1).
design_df <- design_df |>
add_parameter(
b_intv = c(0.2, 0.5),
sigma_site = 0.5,
sigma_unit = 1
)
# Preview to confirm parameter expansion
design_df |>
head()
#> # A tibble: 6 × 9
#> cohort site chron_time condition local_time unit_id b_intv sigma_site
#> <dbl> <chr> <dbl> <fct> <int> <int> <dbl> <dbl>
#> 1 1 A 1 ctrl 0 1 0.2 0.5
#> 2 1 A 1 ctrl 0 1 0.5 0.5
#> 3 1 A 1 ctrl 0 2 0.2 0.5
#> 4 1 A 1 ctrl 0 2 0.5 0.5
#> 5 1 A 1 ctrl 0 3 0.2 0.5
#> 6 1 A 1 ctrl 0 3 0.5 0.5
#> # ℹ 1 more variable: sigma_unit <dbl>
After adding simulation parameters, we expand our design to include replicates. This enables us to simulate multiple trials under the same design to estimate power and evaluate bias efficiently.
Here, we use 10 replicates to keep the example lightweight. In practice,
the number of replicates should be chosen based on the precision you
want for your power estimate. A common guideline is:
n = p*(1-p)/(se^2), where se is the desired standard error of the
power estimate and p is the target power (for example, p = 0.80).
design_df <- design_df |>
initialize_replicates(n = 10)
# Preview structure with replicate IDs
design_df |>
head()
#> # A tibble: 6 × 10
#> sample_id cohort site chron_time condition local_time unit_id b_intv
#> <int> <dbl> <chr> <dbl> <fct> <int> <int> <dbl>
#> 1 1 1 A 1 ctrl 0 1 0.2
#> 2 2 1 A 1 ctrl 0 1 0.2
#> 3 3 1 A 1 ctrl 0 1 0.2
#> 4 4 1 A 1 ctrl 0 1 0.2
#> 5 5 1 A 1 ctrl 0 1 0.2
#> 6 6 1 A 1 ctrl 0 1 0.2
#> # ℹ 2 more variables: sigma_site <dbl>, sigma_unit <dbl>
We now simulate the outcome data according to our specified design. This prepares the dataset for model fitting, power estimation, and bias evaluation.
First, we add a small negative time trend (b_time = -0.05) to
represent gradual improvement over time, independent of intervention.
design_df <- design_df |>
add_parameter(b_time = -0.05)
Next, we add fixed effects for: - The intervention effect, applied when
condition == "intv". - The time trend, scaled by chron_time.
design_df <- design_df |>
add_fixed_effect(
intv_effect = b_intv * as.numeric(condition == "intv")) |>
add_fixed_effect(
time_trend = b_time * chron_time)
We add a random intercept at the site level to simulate between-site variability.
design_df <- design_df |>
add_random_effect(
site_intercept = rnorm(1, mean = 0, sd = sigma_site),
.nesting = "site"
)
We add unit-level residual error with variance sigma_unit^2.
design_df <- design_df |>
add_error(variance = sigma_unit^2)
Finally, we generate a linear outcome (y_linear) by summing all fixed, random, and error components.
design_df <- design_df |>
add_linear_outcome(output_col = "y_linear")
# Preview simulated data
design_df |>
head()
#> # A tibble: 6 × 16
#> sample_id cohort site chron_time condition local_time unit_id b_intv
#> <int> <dbl> <chr> <dbl> <fct> <int> <int> <dbl>
#> 1 1 1 A 1 ctrl 0 1 0.2
#> 2 2 1 A 1 ctrl 0 1 0.2
#> 3 3 1 A 1 ctrl 0 1 0.2
#> 4 4 1 A 1 ctrl 0 1 0.2
#> 5 5 1 A 1 ctrl 0 1 0.2
#> 6 6 1 A 1 ctrl 0 1 0.2
#> # ℹ 8 more variables: sigma_site <dbl>, sigma_unit <dbl>, b_time <dbl>,
#> # .intv_effect <dbl>, .time_trend <dbl>, .site_intercept <dbl>, .error <dbl>,
#> # y_linear <dbl>
With our simulated outcomes generated, we now fit our planned analysis model across all replicates and parameter scenarios.
This process is automatically conducted in parallel, but you’ll need to
load any packages for that parallel process explicitly (see line
packages = "lmerTest").
In this example: - We use linear mixed-effects models to account for
site-level clustering. - Our model includes: - condition (intervention
vs control), - chron_time (continuous time trend), - A random
intercept for site. - We group by sample_id (replicates) and
b_intv (effect size scenario).
This approach allows us to evaluate bias and power across both replicates and effect size conditions.
fitted_models <- design_df |>
group_by(sample_id, b_intv) |>
nest() |>
ungroup() |>
fit_models(
.x = data,
.f = ~ lmer(
formula = y_linear ~ condition + chron_time + (1 | site),
data = .x),
packages = "lmerTest") # load lmerTest in each parallel process
fitted_models |>
head()
#> # A tibble: 6 × 4
#> Loading required namespace: lmerTest
#> sample_id b_intv data model
#> <int> <dbl> <list> <list>
#> 1 1 0.2 <tibble [1,280 × 14]> <lmrMdLmT>
#> 2 2 0.2 <tibble [1,280 × 14]> <lmrMdLmT>
#> 3 3 0.2 <tibble [1,280 × 14]> <lmrMdLmT>
#> 4 4 0.2 <tibble [1,280 × 14]> <lmrMdLmT>
#> 5 5 0.2 <tibble [1,280 × 14]> <lmrMdLmT>
#> 6 6 0.2 <tibble [1,280 × 14]> <lmrMdLmT>
After fitting the models, we now extract key results and evaluate bias and power for our intervention effect.
We will: - Extract model estimates using
extract_model_results(). - Summarize power by calculating the
proportion of p-values below .05 (done by default in
evaluate_model_results). - **Summarize bias** of the intervention effect usingeval_bias(). - **Estimate a critical value** by calculating the observed value of b-estimates above the 97.5th percentile, usingeval_quantile()`.
# Extract tidy model results
model_results <- fitted_models |>
extract_model_results()
# Summarize bias and power
summary_results <- model_results |>
group_by(b_intv, term) |>
evaluate_model_results(
bias = eval_bias(estimate, term = c("(Intercept)" = 0, conditionintv = b_intv)),
critical_val = eval_quantile(estimate, term = c("conditionintv" = 0.975))
)
# View summary of bias and power
summary_results
#> # A tibble: 10 × 8
#> # Groups: b_intv [2]
#> b_intv term n_models mean_estimate mean_std.error power bias critical_val
#> <dbl> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.2 (Inte… 10 0.333 0.110 0.45 -0.155 NA
#> 2 0.2 chron… 10 0.333 0.110 0.45 NA NA
#> 3 0.2 condi… 10 0.333 0.110 0.45 NA 0.328
#> 4 0.2 sd__(… 10 0.333 0.110 0.45 NA NA
#> 5 0.2 sd__O… 10 0.333 0.110 0.45 NA NA
#> 6 0.5 (Inte… 10 0.333 0.110 0.45 -0.192 NA
#> 7 0.5 chron… 10 0.333 0.110 0.45 NA NA
#> 8 0.5 condi… 10 0.333 0.110 0.45 NA 0.590
#> 9 0.5 sd__(… 10 0.333 0.110 0.45 NA NA
#> 10 0.5 sd__O… 10 0.333 0.110 0.45 NA NA