Designing a stepped-wedge trial

library(stepwedgepower)

A stepped-wedge trial assigns clusters to sequences. Every sequence starts under control; sequences then cross over to the intervention at different periods, and once a sequence has crossed it stays crossed. By the end of the trial every cluster has received the intervention.

The package encodes exactly that and nothing more. A cluster can be a hospital ward, a clinic, a school, a general practice, a village, or an individual clinician; a sequence can be a region, a wave, a specialty, or an arbitrary randomisation group. The software does not care.

Design and assumptions are kept apart, because they answer different questions:

[sw_design()] – who crosses over, and when.
[sw_assumptions()] – what the outcome is expected to do.

The design

The most common design gives one baseline period and then crosses one sequence per period.

design <- sw_design(
  clusters_per_sequence = c(8, 8, 8, 8),
  crossover_period = c(2, 3, 4, 5),
  n_periods = 5
)
design
#> <sw_design>
#>   4 sequences, 5 periods, 32 clusters total
#>   clusters per sequence: 8, 8, 8, 8 
#>   treatment schedule (rows = sequences, cols = periods):
#>            Period 1 Period 2 Period 3 Period 4 Period 5
#> Sequence 1        0        1        1        1        1
#> Sequence 2        0        0        1        1        1
#> Sequence 3        0        0        0        1        1
#> Sequence 4        0        0        0        0        1

Designs are not restricted to that shape. Sequences may hold unequal numbers of clusters, several sequences may cross at once, and there may be more than one baseline period. Anything expressible as a 0/1 schedule is expressible here:

schedule <- matrix(
  c(0, 0, 1, 1, 1,
    0, 0, 0, 1, 1,
    0, 0, 0, 0, 1),
  nrow = 3, byrow = TRUE
)

custom <- sw_design(
  clusters_per_sequence = c(8, 10, 12),  # unequal
  treatment = schedule                   # two baseline periods
)
custom
#> <sw_design>
#>   3 sequences, 5 periods, 30 clusters total
#>   clusters per sequence: 8, 10, 12 
#>   treatment schedule (rows = sequences, cols = periods):
#>            Period 1 Period 2 Period 3 Period 4 Period 5
#> Sequence 1        0        0        1        1        1
#> Sequence 2        0        0        0        1        1
#> Sequence 3        0        0        0        0        1

The assumptions

assumptions <- sw_assumptions(
  baseline_prob = 0.05,     # control-arm risk
  treatment_or = 1.6,       # target effect
  icc = 0.05,               # cluster correlation
  period_effects = c(0, 0.05, 0.10, 0.15, 0.20),  # secular trend, log-odds
  n_per_cluster_period = 20
)
assumptions
#> <sw_assumptions>
#>   outcome: binary 
#>   treatment: OR = 1.600 (log-odds 0.470)
#>   cluster_sd = 0.416 (icc = 0.050)
#>   baseline (prob scale): 0.050 
#>   period effects (log-odds): 0.00, 0.05, 0.10, 0.15, 0.20

Three points worth drawing out.

Specify the ICC, not the random-effect SD. Most people can defend an ICC and almost nobody can defend a random-intercept standard deviation. Supply icc and the conversion is handled; supply cluster_sd if you prefer. Supplying both is an error unless they agree.

icc_to_cluster_sd(0.05)
#> [1] 0.4161141
cluster_sd_to_icc(0.4161)
#> [1] 0.04999678

Model the secular trend explicitly. Outcomes often drift over time for reasons unrelated to the intervention. In a stepped-wedge design time and treatment are partially confounded by construction, so a trend you do not simulate is a trend you have not accounted for. period_effects puts it in the generating model, where the analysis has to contend with it.

Cluster sizes need not be equal. Supply a single number, a sequence-by-period matrix, or a function:

variable_size <- sw_assumptions(
  baseline_prob = 0.05, treatment_or = 1.6, icc = 0.05,
  n_per_cluster_period = function(n) rpois(n, lambda = 20) + 1
)
trial <- simulate_swcrt(design, variable_size, seed = 1)
range(trial$n)
#> [1] 10 31

Simulating and analysing one trial

trial <- simulate_swcrt(design, assumptions, seed = 123)
head(trial)
#>   cluster_id   sequence sequence_idx period intervention  n events  true_prob
#> 1          1 Sequence 1            1      1            0 20      2 0.04001509
#> 2          1 Sequence 1            1      2            1 20      1 0.06551863
#> 3          1 Sequence 1            1      3            1 20      2 0.06864724
#> 4          1 Sequence 1            1      4            1 20      2 0.07191375
#> 5          1 Sequence 1            1      5            1 20      2 0.07532312
#> 6          2 Sequence 1            1      1            0 20      1 0.04564163

The columns are deliberately neutral: cluster_id, sequence, period, intervention, n, events. Any aggregated cluster-period dataset in this shape can be fed to the analysis function, whether it came from simulate_swcrt() or from a real trial.

fit <- analyze_swcrt(trial)
fit$estimate    # log-odds
#> [1] 0.5249781
fit$p_value
#> [1] 0.01214173

The default analysis is the standard Hussey-Hughes model: treatment, a categorical period effect, and a cluster random intercept. You can rename any column, or hand over a formula outright:

fit_stepwedge_model(
  trial,
  formula = cbind(events, n - events) ~ intervention + factor(period) +
    (1 | cluster_id)
)

Power

power <- power_swcrt(design, assumptions, nsim = 1000, seed = 1, n_cores = 4)
summary(power)
plot(power)

The returned object reports more than a single number, and the extras are not decoration. The Monte Carlo interval tells you how much of the estimate is simulation noise. The failure and singular-fit rates tell you whether the analysis model is actually estimable at this design size – a “power” computed from replicates that mostly failed to converge is not a power. Bias and coverage tell you whether the estimator is behaving.

The interval for power is exact (Clopper-Pearson) rather than Wald. A Wald interval collapses to zero width when the estimate hits 0 or 1, which is exactly the high-power regime a power analysis is trying to land in.

Sensitivity

Power depends on assumptions you cannot verify in advance, so vary them.

grid <- power_grid(
  design, assumptions,
  vary = list(
    treatment_or = c(1.3, 1.5, 1.75, 2.0),
    icc = c(0.02, 0.05, 0.10),
    clusters_per_sequence = list(c(6, 6, 6, 6), c(8, 8, 8, 8), c(10, 10, 10, 10))
  ),
  nsim = 500, n_cores = 4, seed = 1
)
summary(grid)
plot(grid)

Simulations are reproducible under a given seed regardless of n_cores: each replicate draws from its own random-number stream, so adding cores makes the run faster without moving the answer.

Checking the engine

validate_engine() re-runs the power calculation with an independently written reference implementation that shares no code with the main engine, and reports whether the two agree within Monte Carlo error.

validate_engine(design, assumptions, nsim = 500, seed = 1)

compare_with_swdpwr() will also line the simulation up against a closed-form calculation, if you have that package installed. Treat the result as information rather than as a target: closed-form stepped-wedge power is asymptotic in the number of clusters, and with the small cluster counts these trials usually run it can report power well above what a simulated trial achieves. Where the two disagree, the simulation is the one that reflects the model you will actually fit.