Compares two independent groups by referring the observed statistic to the distribution of that same statistic over reassignments of the observed scores to the two groups. That reference distribution, not a normal or a t distribution, supplies the p-value, so the test needs no assumption about the shape of the population. Alongside the test, the function reports the effect sizes that answer the question the test only screens: the mean difference and its randomization-based interval, the standardized mean difference with a noncentral t interval, the common language effect size, and Cliff's delta.
randomization_test(
x = NULL,
group = NULL,
data = NULL,
group_1 = NULL,
group_2 = NULL,
statistic = c("mean", "t"),
alternative = c("two_sided", "less", "greater"),
exact = NULL,
n_resamples = 10000L,
seed = NULL,
conf_level = 0.95,
shift_ci = TRUE
)A data.frame with a term column and a numeric
value column, in three blocks.
The test: mean_difference (first group minus second),
statistic (the statistic actually referred to the reference
distribution), p_value, and p_value_se (the Monte Carlo
standard error of the p-value, NA under exact
enumeration, which has no Monte Carlo error).
The intervals and effect sizes: shift_lower_limit and
shift_upper_limit (the randomization interval for the shift,
obtained by inverting the test); normal_theory_lower_limit and
normal_theory_upper_limit (Welch's t interval on the same
mean difference, reported for contrast); smd with
smd_lower_limit and smd_upper_limit from
ci_smd; cles with cles_lower_limit and
cles_upper_limit from cles; and
cliff_delta with cliff_delta_lower_limit and
cliff_delta_upper_limit from cliff_delta.
The design: n_1, n_2, N, n_evaluated (how
many reassignments were actually used), and exact (1 if every
reassignment was enumerated, 0 if they were sampled).
Non-numeric information travels on attributes rather than in the
value column: statistic_name, method
(“exact enumeration” or “Monte Carlo”),
alternative, group_labels, response_name,
group_name, seed, observed_statistic, and
reference_distribution, the vector of statistics over the
reassignments that plot_randomization_test draws.
A formula of the form y ~ group, or a numeric response
vector to be paired with group. Leave NULL when the two
samples are supplied through group_1 and group_2.
A grouping variable the same length as x, with
exactly two levels after unused levels are dropped. Ignored when
x is a formula.
An optional data.frame in which to find the variables
named in the formula.
The two samples supplied directly as numeric vectors, an alternative to the formula and response-plus-grouping interfaces. Lengths need not be equal.
One of "mean" (default) or "t".
"mean" uses the difference in group means. "t" uses the
studentized (Welch) statistic, which divides that difference by its
separate-variances standard error and is the better choice when the
groups may differ in variance (see Details).
One of "two_sided" (default; the base-R
spelling "two.sided" is accepted as an alias), "less",
or "greater". The direction refers to the first group minus the
second, the same orientation t.test uses.
Logical. If NULL (default), the test enumerates
every reassignment when choose(N, n_1) is at most 50,000 and
samples reassignments otherwise. If TRUE, enumeration is forced
(refused above 1,000,000 reassignments). If FALSE, Monte Carlo
is forced.
Number of randomly drawn reassignments when
enumeration is not used. Default 10000L.
Optional integer seed for the Monte Carlo branch. Default
NULL, which leaves the user's current RNG state intact; supply
an integer for reproducibility. When a seed is supplied the RNG state
in place before the call is restored on exit.
Confidence level for every interval reported, the
inverted randomization interval included. Default 0.95.
Logical. Compute the randomization-based interval for
the shift by inverting the test? Default TRUE. Setting it to
FALSE reports the two endpoints as NA and skips the
inversion, which is the expensive part of the call.
Ken Kelley [email protected]
What the randomization distribution is. Suppose N participants were randomly assigned, \(n_1\) to one condition and \(n_2\) to the other. Under the null hypothesis that the condition a participant received made no difference to that participant's score, each score would have been the same number no matter which group the participant landed in. The assignment actually used was one draw from the \(\binom{N}{n_1}\) assignments the randomization could equally well have produced, so every one of those assignments was equally likely, and each of them yields a value of the test statistic. Those values are the randomization distribution. The p-value is the proportion of them at least as extreme as the value the experiment actually produced.
Why there is no normality assumption. Nothing in that argument mentions a population, a normal curve, or a sampling model. The probability comes from the coin flips the experimenter performed, which are known exactly because the experimenter performed them. This is the inferential logic Fisher (1935) used to introduce experimental design, and it is where Chapter 1 of Maxwell, Delaney, and Kelley (2027) starts, for the same reason: the validity of the test rests on the randomization rather than on assumptions a data analyst cannot check.
What the test does and does not license. A small p-value says the observed separation between the groups would rarely arise from reassignment alone, which is evidence that the assignment mattered. It does not say how much it mattered, and with a large N an uninteresting difference will produce a small p-value. It also does not, by itself, license generalization beyond the participants at hand: randomization licenses a causal claim about these units, while generalization to a population is a separate argument that rests on how the units were recruited. That is why this function reports effect sizes with intervals rather than a p-value alone.
Why the studentized statistic. With \(n_1 = n_2\) and equal
population variances the two statistics give the same p-value to
within the discreteness of the reference distribution, because the
denominator of the studentized statistic is then nearly constant across
reassignments. When the variances differ and the groups are unbalanced
they part company. Reassigning scores between groups of unequal size
mixes the two variances in proportions that the observed assignment does
not have, so the reference distribution for the raw mean difference is
built under a null that is false in a second way, and the test's actual
Type I error rate drifts away from the nominal level. The studentized
statistic rescales each reassignment by its own separate-variances
standard error, which removes most of that drift and remains
asymptotically valid under heteroscedasticity (Janssen, 1997; Neuhaus,
1993). Use statistic = "t" whenever unequal variances are
plausible, which for unbalanced designs is nearly always.
Exact or Monte Carlo. When choose(N, n_1) is at most
50,000 every reassignment is enumerated and the p-value is exact:
it is a count divided by a known total, with no approximation anywhere.
Above that threshold n_resamples reassignments are drawn at
random and the p-value is
\((r + 1) / (m + 1)\), where r counts the sampled reassignments
at least as extreme as the observed one and m is
n_resamples. Adding one to each part counts the observed
assignment, which is itself a legitimate reassignment; without it a
p-value of exactly zero could be reported for a hypothesis the
data cannot rule out, and the test would be anticonservative (Phipson &
Smyth, 2010). The reported p_value_se is
\(\sqrt{\hat p (1 - \hat p) / m}\), the standard error of the
resampling itself. It describes how much the p-value would move
if the reassignments were drawn again, not how much it would move in a
new experiment. Raising n_resamples shrinks it at the usual
\(1/\sqrt{m}\) rate.
The randomization interval, and how it differs from the normal
theory one. Suppose the treatment adds a constant \(\delta\) to every
score it touches. Subtracting \(\delta\) from each first-group score
should then leave scores that are exchangeable across groups, so the
randomization test applied to the subtracted data is a test of
\(H_0\!: \mathrm{shift} = \delta\). The set of \(\delta\) for which
that test does not reject at level \(1 - \) conf_level is a
confidence interval for the shift, and it is reported as
shift_lower_limit and shift_upper_limit. Inverting a test
this way is the general recipe (Ernst, 2004); the endpoints are located
by bisection on the p-value, using the same reassignments
throughout so the interval and the test agree.
The contrast with normal_theory_lower_limit and
normal_theory_upper_limit, which are Welch's t limits on
the same mean difference, is worth reading whenever both are printed.
The randomization interval is exactly the set of shifts the test being
run does not reject, so the test and the interval can never disagree.
The normal theory interval instead assumes the sampling distribution of
the mean difference has a known shape; it is smooth, symmetric about the
point estimate, and can extend past the range the data can support.
The randomization interval is discrete, need not be symmetric, and in a
very small design is unbounded, a correct statement of how little
information the design carries rather than a defect:
with three observations per group the smallest attainable two-sided
p-value is 2/20 = 0.10, so no shift can be rejected at the 5%
level and the 95% interval is the whole real line. The randomization
interval also inherits the shift model, so it answers a narrower
question than the test does: the test needs only exchangeability, while
the interval needs the treatment to move every score by the same amount.
Effect sizes. Every effect size reported here comes from the
package function that owns it, so the numbers match a direct call.
smd and ci_smd supply the standardized mean
difference and its noncentral t interval; cles
supplies the common language effect size, the probability that a
randomly drawn score from the first group exceeds one from the second,
by transforming those limits through \(\Phi(\cdot/\sqrt 2)\); and
cliff_delta supplies Cliff's delta with its consistent
interval. The standardized mean difference and the common language
effect size are normal theory quantities, so their intervals lean on the
assumption the test itself avoids. Cliff's delta does not: it is a
function of the ordering of the observations alone, which makes it the
natural effect size companion to a randomization test. Reporting all
three lets a reader see whether the distribution-free and normal theory
summaries tell the same story.
Ranks give the Wilcoxon test. Replacing the scores by their
ranks and running this test with statistic = "mean" reproduces
the exact Wilcoxon rank sum test, since the rank sum is a monotone
function of the difference in mean ranks. That equivalence is a useful
check and a reminder of what the rank test is: a randomization test on
transformed data.
Edgington, E. S., & Onghena, P. (2007). Randomization tests (4th ed.). Chapman & Hall/CRC.
Ernst, M. D. (2004). Permutation methods: A basis for exact inference. Statistical Science, 19(4), 676--685. tools:::Rd_expr_doi("10.1214/088342304000000396")
Fisher, R. A. (1935). The design of experiments. Oliver & Boyd.
Janssen, A. (1997). Studentized permutation tests for non-i.i.d. hypotheses and the generalized Behrens-Fisher problem. Statistics & Probability Letters, 36(1), 9--21. tools:::Rd_expr_doi("10.1016/S0167-7152(97)00043-6")
Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 1 on the logic of randomization and the randomization test.)
Neuhaus, G. (1993). Conditional rank tests for the two-sample problem under random censorship. The Annals of Statistics, 21(4), 1760--1779. tools:::Rd_expr_doi("10.1214/aos/1176349396")
Phipson, B., & Smyth, G. K. (2010). Permutation p-values should never be zero: Calculating exact p-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology, 9(1), Article 39. tools:::Rd_expr_doi("10.2202/1544-6115.1585")
Pitman, E. J. G. (1937). Significance tests which may be applied to samples from any populations. Supplement to the Journal of the Royal Statistical Society, 4(1), 119--130.
plot_randomization_test for the figure that shows
the reference distribution, the observed statistic, and the rejection
region; randomization_test_paired for the sign-flip
sibling used with paired observations;
t.test and wilcox.test for
the parametric and rank-based alternatives; smd,
ci_smd, cles, and
cliff_delta for the effect sizes reported here.
Other hypothesis tests:
adjusted_means(),
ancova(),
anova_within(),
ci_dunnett(),
ci_scheffe(),
ci_tukey_kramer(),
compare_cov_structures(),
contrast_test(),
correlations_test(),
equivalence_r(),
equivalence_smd(),
factorial_anova(),
manova_split_plot(),
mauchly_test(),
mixed_anova(),
obrien_test(),
pairwise_within(),
randomization_test_paired(),
regions_of_significance(),
simple_effects_AB(),
summary_t_test(),
welch_t()
# 1. Ten observations, so every one of the choose(10, 5) = 252
# reassignments is enumerated and the p-value is exact.
treatment <- c(80, 84, 79, 88, 83)
control <- c(72, 75, 68, 81, 74)
randomization_test(group_1 = treatment, group_2 = control)
# 2. The studentized statistic, preferable when the groups may differ
# in variance.
randomization_test(group_1 = treatment, group_2 = control,
statistic = "t")
# 3. Formula interface: weekly drinking in the two comparable arms of
# the drinks_trial data, a right-skewed outcome, which is exactly
# where a distribution-free test earns its keep. With 37 and 32
# participants there are far too many reassignments to enumerate, so
# 10,000 are drawn and the p-value carries a Monte Carlo standard
# error.
cra <- droplevels(subset(drinks_trial, treatment != "CRA + Disulfiram"))
set.seed(113)
randomization_test(drinks_per_week ~ treatment, data = cra, seed = 113)
# 4. A one-sided test, and the one-sided interval that goes with it.
randomization_test(group_1 = treatment, group_2 = control,
alternative = "greater")
# 5. On ranks, the test is the exact Wilcoxon rank sum test.
y <- c(treatment, control)
g <- rep(c("treatment", "control"), each = 5)
res <- randomization_test(rank(y), g)
res$value[res$term == "p_value"]
wilcox.test(treatment, control, exact = TRUE)$p.value
Run the code above in your browser using DataLab