randomizr 2.0.1 was written by Alexander Coppock working with Claude
(Anthropic), across design, implementation, tests and documentation,
with contributions to balanced_ra() and cube-on-X by
Macartan Humphreys working with Cursor (various models).
While we have reviewed the code base for 2.0.1, we did not write it.
The guarantee offered here then is not that every line of code has been
vouched for. It is that the package can be shown to do what it says: the
table below states each promise in words and then measures it, live, in
this document. Every function reproduces randomizr 1.0.1’s draws for a
given seed as well, with one documented exception covering
strata_rs() and strata_and_cluster_rs(). That
is checked separately, in
tests/testthat/test_stream_compat.R, against output
recorded from an installed 1.0.1. The package contains many additional
tests in tests/testthat.
A random assignment function makes two kinds of promise and we check these in different ways.
An exact promise holds on every draw or it is
broken. For example, that exactly m units are
treated, that no cluster is split across conditions, that every block
gets its stated count. There is no tolerance to set: the check runs the
design a few thousand times and checks that the property holds every
time.
A probability promise is about a long-run average, so it can
only be checked against sampling error. That each unit is
treated with probability 0.25 does not mean any particular set of draws
will show exactly 0.25. With sims draws, a single unit’s
estimated probability has a standard error of
sqrt(p(1-p)/sims), and the check asks whether every unit
sits within five of those of its target. Five is chosen so that a
correct package essentially never fails the check by chance, while a
wrong marginal probability, which would be off by a systematic amount
rather than a Monte Carlo one, is still caught.
library(randomizr)
set.seed(20260824)
sims <- 2000
# An exact promise: the property either held on every draw or it did not.
exact <- function(draw, property) {
ok <- replicate(sims, isTRUE(property(draw())))
list(result = sprintf("holds on %d of %d draws", sum(ok), sims),
holds = all(ok))
}
# A probability promise: no unit's realized rate may sit more than five Monte
# Carlo standard errors from the probability the package claims for it.
probability <- function(draw, target, condition = 1) {
realized <- rowMeans(replicate(sims, draw()) == condition)
gap <- max(abs(realized - target))
tolerance <- 5 * sqrt(max(target * (1 - target)) / sims)
list(result = sprintf("largest gap %.3f, tolerance %.3f", gap, tolerance),
holds = gap < tolerance)
}One hundred units throughout: five blocks of twenty for the blocked
designs, twenty clusters of five for the clustered ones, and for
balanced_ra() a probability that differs for every
unit.
N <- 100
blocks <- rep(c("a", "b", "c", "d", "e"), each = 20)
clusters <- rep(1:20, each = 5)
regions <- rep(c("north", "south"), each = 50)
p_unit <- seq(0.05, 0.95, length.out = N)
x <- rnorm(N)
# helpers the properties below are written in terms of
by_group <- function(z, g) tapply(z, g, sum)
uniform_within <- function(z, g) all(tapply(z, g, function(x) length(unique(x))) == 1)checks <- list(
list("simple_ra()", "each unit is treated with probability 0.25",
probability(function() simple_ra(N, prob = 0.25), 0.25)),
list("simple_ra()", "the number treated is random, not fixed",
local({
m <- replicate(sims, sum(simple_ra(N, prob = 0.25)))
list(result = sprintf("sd %.2f against a binomial %.2f, %d distinct counts",
sd(m), sqrt(N * 0.25 * 0.75), length(unique(m))),
holds = abs(sd(m) / sqrt(N * 0.25 * 0.75) - 1) < 0.1)
})),
list("complete_ra()", "each unit is treated with probability 25/100",
probability(function() complete_ra(N, m = 25), 0.25)),
list("complete_ra()", "exactly 25 units are treated",
exact(function() complete_ra(N, m = 25), function(z) sum(z) == 25)),
list("complete_ra(num_arms = 3)", "the three arms get 33, 33 and 34 units",
exact(function() complete_ra(N, num_arms = 3),
function(z) all(sort(as.vector(table(z))) == c(33, 33, 34)))),
list("complete_ra(num_arms = 3)", "each unit reaches arm T1 with probability 1/3",
probability(function() complete_ra(N, num_arms = 3), 1 / 3, condition = "T1")),
list("block_ra()", "exactly half of every block of 20 is treated",
exact(function() block_ra(blocks = blocks, prob = 0.5),
function(z) all(by_group(z, blocks) == 10))),
list("block_ra()", "each unit is treated with its block's probability, 0.5",
probability(function() block_ra(blocks = blocks, prob = 0.5), 0.5)),
list("cluster_ra()", "no cluster is ever split across conditions",
exact(function() cluster_ra(clusters = clusters),
function(z) uniform_within(z, clusters))),
list("cluster_ra()", "each of the 20 clusters is treated with probability 0.5",
probability(function() tapply(cluster_ra(clusters = clusters), clusters, `[`, 1), 0.5)),
list("block_and_cluster_ra()", "clusters stay intact and each region treats 25 units",
exact(function() block_and_cluster_ra(clusters = clusters, blocks = regions),
function(z) uniform_within(z, clusters) && all(by_group(z, regions) == 25))),
list("balanced_ra()", "each unit is treated with its own probability",
probability(function() balanced_ra(prob_unit = p_unit), p_unit)),
list("balanced_ra()", "the number treated is the floor or ceiling of its target",
exact(function() balanced_ra(prob_unit = p_unit),
function(z) sum(z) %in% c(floor(sum(p_unit)), ceiling(sum(p_unit))))),
list("balanced_ra(formula = ~ x)", "each unit is treated with its own probability",
probability(function() balanced_ra(prob_unit = p_unit, formula = ~ x), p_unit)),
list("balanced_ra(formula = ~ x)", "the number treated is the floor or ceiling of its target",
exact(function() balanced_ra(prob_unit = p_unit, formula = ~ x),
function(z) sum(z) %in% c(floor(sum(p_unit)), ceiling(sum(p_unit))))),
list("declare_ra() + conduct_ra()", "realized rates match the declared probabilities",
local({
declaration <- declare_ra(
blocks = blocks,
block_m_each = rbind(c(15, 5), c(14, 6), c(13, 7), c(12, 8), c(11, 9)))
probability(function() conduct_ra(declaration),
declaration$probabilities_matrix[, 2])
})),
list("complete_rs()", "exactly 25 units are sampled",
exact(function() complete_rs(N, n = 25), function(s) sum(s) == 25)),
list("strata_rs()", "exactly half of every stratum is sampled",
exact(function() strata_rs(strata = blocks, prob = 0.5),
function(s) all(by_group(s, blocks) == 10)))
)
guarantees <- data.frame(
Function = vapply(checks, function(x) x[[1]], character(1)),
Promise = vapply(checks, function(x) x[[2]], character(1)),
Measured = vapply(checks, function(x) x[[3]]$result, character(1)),
Holds = vapply(checks, function(x) x[[3]]$holds, logical(1))
)| Function | Promise | Measured | Holds |
|---|---|---|---|
| simple_ra() | each unit is treated with probability 0.25 | largest gap 0.025, tolerance 0.048 | TRUE |
| simple_ra() | the number treated is random, not fixed | sd 4.26 against a binomial 4.33, 30 distinct counts | TRUE |
| complete_ra() | each unit is treated with probability 25/100 | largest gap 0.029, tolerance 0.048 | TRUE |
| complete_ra() | exactly 25 units are treated | holds on 2000 of 2000 draws | TRUE |
| complete_ra(num_arms = 3) | the three arms get 33, 33 and 34 units | holds on 2000 of 2000 draws | TRUE |
| complete_ra(num_arms = 3) | each unit reaches arm T1 with probability 1/3 | largest gap 0.026, tolerance 0.053 | TRUE |
| block_ra() | exactly half of every block of 20 is treated | holds on 2000 of 2000 draws | TRUE |
| block_ra() | each unit is treated with its block’s probability, 0.5 | largest gap 0.029, tolerance 0.056 | TRUE |
| cluster_ra() | no cluster is ever split across conditions | holds on 2000 of 2000 draws | TRUE |
| cluster_ra() | each of the 20 clusters is treated with probability 0.5 | largest gap 0.025, tolerance 0.056 | TRUE |
| block_and_cluster_ra() | clusters stay intact and each region treats 25 units | holds on 2000 of 2000 draws | TRUE |
| balanced_ra() | each unit is treated with its own probability | largest gap 0.043, tolerance 0.056 | TRUE |
| balanced_ra() | the number treated is the floor or ceiling of its target | holds on 2000 of 2000 draws | TRUE |
| balanced_ra(formula = ~ x) | each unit is treated with its own probability | largest gap 0.033, tolerance 0.056 | TRUE |
| balanced_ra(formula = ~ x) | the number treated is the floor or ceiling of its target | holds on 2000 of 2000 draws | TRUE |
| declare_ra() + conduct_ra() | realized rates match the declared probabilities | largest gap 0.026, tolerance 0.056 | TRUE |
| complete_rs() | exactly 25 units are sampled | holds on 2000 of 2000 draws | TRUE |
| strata_rs() | exactly half of every stratum is sampled | holds on 2000 of 2000 draws | TRUE |
Every promise in the table is met, on 2000 draws of each design, at the seed set at the top of this document. The numbers are computed when the vignette is built, so they are the numbers your installed copy produces rather than numbers recorded from a run somewhere else.
The table is a demonstration, not a proof, and it is deliberately a small one.
It checks one configuration per promise. The test suite checks many: it runs the same properties across argument forms, arm counts, block and cluster structures, and the edge cases that have caused bugs, which is where a guarantee is actually enforced. What this document adds is that the promises are stated in words a reader can disagree with, and measured where a reader can see them.
It says nothing about designs it does not draw.
balanced_ra() with three or more arms and
blocks is the case worth naming, because there the
within-block counts stay tight and the overall count does not: it is the
sum of independent block landings and can move by more than one. That is
a documented limit of the implementation rather than an oversight, and
the “Introduction to balanced_ra” vignette gives the reason. The count
promise also carries one documented exception, a rare floating-point
path these draws did not exercise; Section 6 of that vignette describes
it.
It does not check the distribution of assignments beyond its first
moment, and neither does the test suite. That each unit is treated with
the right probability leaves open whether whole assignment vectors are
drawn with the right joint distribution. What the package offers on that
is enumeration rather than measurement:
obtain_permutation_matrix() lists every assignment a
declaration can produce and
obtain_permutation_probabilities() gives each one its
probability, and the suite confirms that the enumeration is complete,
that no assignment is listed twice, and that the probabilities sum to
one. It does not draw a design many times and compare the realized
frequencies of those assignments against those probabilities. For
randomization inference that gap does not bite, because
conduct_ra() and the permutation machinery are the same
code path. It is still the part of this document that is enumerated
rather than shown.