The formula interface handles multi-arm trials directly. When the
treatment variable takes more than two values, conduct_ri
produces one result per pairwise comparison against the baseline
condition (the first condition in sorted order). The null distribution
for each comparison is constructed by re-randomizing only the units
assigned to those two conditions while leaving the remaining arms at
their observed values. Each arm therefore gets its own permutation
draws, reflecting the conditional design.
The summary and plot return one row (or panel) per comparison. Here
T2 has a detectable effect while T3 does
not.
Why conditioning on the two-arm comparison matters
The pairwise null distribution for T2 vs T1 is constructed by
re-randomizing only the T1 and T2 units. A natural but wrong alternative
is to permute all three arms and then extract the T2-vs-T1 difference
from each permutation. The problem: T3-assigned units carry their T3
potential outcomes. When a full permutation shuffles them into T1 or T2,
their outcomes pollute the null distribution, since those outcomes
reflect the T3 effect rather than anything about T2 vs T1.
The contamination is most visible when T3 has a large effect. Here T3
has an effect of 3, while T2 has an effect of 0.4:
The two approaches differ only in how they build a simulated
assignment. The naive one draws a fresh assignment for all three arms.
The correct one shuffles the labels of the two arms being compared among
the units that actually received them, leaving every other unit at its
observed arm. permute_within does the shuffling:
permute_within <- function(z, arms) {
shuffled <- z %in% arms
z[shuffled] <- sample(z[shuffled])
z
}
Rather than simulate one assignment at a time, we stack every
simulation into one long data frame: sims copies of the
experiment, each with its own simulated assignment, then a single
grouped summary at the end. arm_dims takes such a frame and
returns the simulated difference-in-means for one arm against the T1
baseline:
arm_dims <- function(draws, arm) {
draws |>
filter(Z_sim %in% c("T1", arm)) |>
summarize(mean_Y = mean(Y), .by = c(approach, sim, Z_sim)) |>
pivot_wider(names_from = Z_sim, values_from = mean_Y) |>
mutate(est_sim = .data[[arm]] - T1)
}
set.seed(7)
declaration_contam <- declare_ra(N = 90, num_arms = 3)
dat_contam <-
tibble(Z = conduct_ra(declaration_contam)) |>
mutate(Y = 0.4 * (Z == "T2") + 3.0 * (Z == "T3") + rnorm(n()))
obs_dim <-
dat_contam |>
mutate(approach = "observed", sim = 1, Z_sim = Z) |>
arm_dims("T2") |>
pull(est_sim)
null_dims <-
bind_rows(
"Naive: full permutation" =
expand_grid(sim = 1:1000, dat_contam) |>
mutate(Z_sim = conduct_ra(declaration_contam), .by = sim),
"Correct: conditional permutation" =
expand_grid(sim = 1:1000, dat_contam) |>
mutate(Z_sim = permute_within(Z, c("T1", "T2")), .by = sim),
.id = "approach"
) |>
arm_dims("T2")
The naive null distribution is much wider because units with \(Y \approx 3\) (the T3 units) randomly land
in T1 or T2 across permutations, adding variance that has nothing to do
with the T2 vs T1 comparison:
gg_df <-
null_dims |>
mutate(approach = fct_relevel(approach, "Naive: full permutation"))
ggplot(gg_df, aes(x = est_sim)) +
geom_histogram(bins = 40, fill = "grey80", colour = "white") +
geom_vline(xintercept = obs_dim, colour = "steelblue", linewidth = 1) +
facet_wrap(~ approach) +
labs(x = "Simulated T2 vs T1 difference-in-means", y = NULL) +
theme_bw() +
theme(legend.position = "none")
null_dims |>
summarize(sd_null = sd(est_sim),
p_value = mean(abs(est_sim) >= abs(obs_dim)),
.by = approach)
## # A tibble: 2 × 3
## approach sd_null p_value
## <chr> <dbl> <dbl>
## 1 Naive: full permutation 0.455 0.235
## 2 Correct: conditional permutation 0.244 0.022
The naive null’s inflated variance swamps the T2 signal: its p-value
is too large. The correct null’s narrower spread gives a p-value that
reflects the actual T2 effect. conduct_ri() always uses the
conditional approach.
The direction of the distortion, however, is not fixed. It depends on
whether T3 units add or remove variance when mixed into the T2/T1
comparison. The previous example had high-variance T3 units, which
widened the naive null and made it conservative. A different but equally
plausible design produces the opposite: a leveling
treatment, where T3 brings all units to a similar outcome
level, eliminating baseline heterogeneity. T3 units then have very low
outcome variance. When a full permutation shuffles them into the T2 or
T1 group, they anchor the group means, narrowing the naive null
distribution for the T2 comparison. The same effect that made naive
conservative before now makes it anti-conservative.
At the same time, high-variance T1/T2 units contaminate the T3
comparison when they shuffle in, widening that null and reducing power
precisely where we have a real effect to detect.
The same two helpers do the work here. One replication draws an
experiment, builds both null distributions for both comparisons, and
returns four p-values:
one_replication <- function(declaration, sims = 300) {
dat_lev <-
tibble(Z = conduct_ra(declaration)) |>
mutate(Y = if_else(Z == "T3",
rnorm(n(), mean = 1.5, sd = 0.1),
rnorm(n(), mean = 0, sd = 6)))
draws <- function(arm) {
bind_rows(
naive =
expand_grid(sim = 1:sims, dat_lev) |>
mutate(Z_sim = conduct_ra(declaration), .by = sim),
correct =
expand_grid(sim = 1:sims, dat_lev) |>
mutate(Z_sim = permute_within(Z, c("T1", arm)), .by = sim),
.id = "approach"
) |>
arm_dims(arm) |>
mutate(arm = arm)
}
observed <-
map(c("T2", "T3"), \(arm)
dat_lev |>
mutate(approach = "observed", sim = 1, Z_sim = Z) |>
arm_dims(arm) |>
mutate(arm = arm)) |>
list_rbind() |>
select(arm, est_obs = est_sim)
map(c("T2", "T3"), draws) |>
list_rbind() |>
left_join(observed, by = "arm") |>
summarize(p_value = mean(abs(est_sim) >= abs(est_obs)), .by = c(arm, approach))
}
Running it 300 times gives the rejection rate of each approach. T2’s
sharp null is true, so its rejection rate should sit at 0.05. T3 has a
real effect, so a higher rate there is power:
set.seed(42)
declaration_lev <- declare_ra(N = 90, num_arms = 3)
leveling_results <-
tibble(replication = 1:300) |>
mutate(p_values = map(replication, \(r) one_replication(declaration_lev))) |>
unnest(p_values) |>
summarize(rejection_rate = mean(p_value < 0.05), .by = c(arm, approach))
leveling_results
## # A tibble: 4 × 3
## arm approach rejection_rate
## <chr> <chr> <dbl>
## 1 T2 naive 0.0967
## 2 T2 correct 0.0533
## 3 T3 naive 0.17
## 4 T3 correct 0.253
The naive T2 test rejects at roughly twice the nominal rate: the
low-variance T3 units narrow its null distribution, making modest
observed T2 differences look extreme. The naive T3 test has
substantially lower power: high-variance T1/T2 units inflate its null
distribution, making the real T3 effect harder to detect. The correct
conditional approach maintains the nominal type I error rate for T2 and
the correct power for T3 in both scenarios. The conditional permutation
is not a technicality: it is what makes the test valid.
