Code
set.seed(1)
x <- c(rnorm(15), rnorm(15, mean = 0.6))
g <- rep(c("a", "b"), each = 15)
T_obs <- diff(tapply(x, g, mean))
B <- 1999
T_perm <- replicate(B, diff(tapply(x, sample(g), mean)))
(sum(abs(T_perm) >= abs(T_obs)) + 1) / (B + 1)[1] 0.106Multivariate comparisons: permutation and distance-based tests
Exchangeability, what PERMANOVA actually tests, and permuting at the right level.
Before reading: Basic testing and the independent unit in the study design.
1 State the multivariate question
A profile contains many measured features, but “the groups differ” can mean different things: a mean-vector contrast on a chosen scale, a dispersion difference, or a broader distributional difference. A distance-based statistic is a summary chosen to detect some alternatives; it is not itself an effect-size estimand with a universal biological meaning.
Specify the population, features, preprocessing, distance, and sampling unit. A whole-profile test and thousands of feature-wise tests answer different questions. A significant omnibus result does not identify which features differ or establish that group centroids alone changed.
Write a person’s profile as \mathbf Y_i=(Y_{i1},\ldots,Y_{ip}). A mean-vector question compares E[\mathbf Y\mid A=1] with E[\mathbf Y\mid A=0]. A dispersion question compares spread around the group centers. Equality of the whole distributions is stronger than either: it excludes changes in means, spread, tails, and dependence among features. Keep that distinction visible when choosing the null that will justify a permutation reference.
2 The permutation idea
You have a statistic T that measures how different two groups look. If the null and design make the allowed label assignments equally likely, you can hold the observed profiles fixed and ask where the observed statistic falls among those assignments. This is a conditional reference distribution for this dataset. The exchangeability or randomization argument is what turns that conditional calculation into a test with a repeated-sampling error guarantee.
Two things to get right from the start.
Exchangeability is the justification for the relabeling, not a synonym for every scientific null. The permutation argument requires that, under the null, the joint distribution of the data is unchanged by relabeling. That is true when the two groups have the same distribution. It is not true when they have the same mean but different variances, and a permutation test of a mean difference can then reject at well above its nominal rate (the Behrens-Fisher problem does not go away by permuting). For the simple unrestricted label-permutation test here, equal group distributions justify exchangeability. More general randomization, restricted, or suitably studentized permutation procedures can address other nulls under their own conditions; identify the justification for the procedure actually used.
The +1 in the p-value. With B random permutations, the p-value (b + 1)/(B + 1) includes the observed arrangement among the permutations and gives a valid, possibly conservative, Monte Carlo permutation p-value under the required exchangeability and sampling scheme (Phipson and Smyth 2010). Its smallest possible value is 1/(B + 1), so with 199 permutations you cannot see a p-value below 0.005, and a simulation study of level needs B large enough that the rejection threshold is well resolved. Reject when p \le \alpha, not p < \alpha: with B=99, the grid includes 0.05, whereas p<0.05 excludes it. Exact rejection probabilities depend on the valid randomization scheme, discreteness, and ties; the grid alone does not prove exact nominal size.
If all M allowed assignments are enumerated, including the observed one, the exact p-value is the number at least as extreme divided by M. The plus-one formula is for a random sample of assignments; do not add another observed arrangement to an already complete enumeration. With 12 people assigned six to each group, there are \binom{12}{6}=924 distinct label assignments. If each person has three visits, there are still 924 person-level assignments, rather than assignments of 36 independent visit labels.
3 What is being permuted
The permutation argument applies to units that are exchangeable under the null. Observations within the same subject are not exchangeable with observations from other subjects when subjects differ (an intraclass correlation above zero). If the factor you are testing is a subject-level factor, such as disease group, the exchangeable units are subjects, and the labels must be permuted across subjects, carrying each subject’s observations with them.
Permuting at the observation level treats N_{\text{obs}} correlated measurements as if they were N_{\text{obs}} independent ones. The test then rejects a true null far more often than it should.
Code
set.seed(2)
# Subjects nested in groups; `visits` observations per subject; intraclass correlation `icc`.
SimulateVisits <- function(n_subj = 20, visits = c(2, 6), icc = 0.5,
design = c("independent", "group-linked")) {
design <- match.arg(design); J <- 2 * n_subj
grp <- sample(rep(c("ctrl", "pat"), each = n_subj))
if (length(visits) == 1) {
M <- rep(visits, J)
} else if (design == "independent") {
M <- sample(rep(visits, length.out = J)) # schedule independent of assignment
} else {
M <- ifelse(grp == "ctrl", visits[1], visits[2]) # breaks record exchangeability
}
subj <- rep(seq_len(J), times = M)
u <- rnorm(J, sd = sqrt(icc))
y <- u[subj] + rnorm(length(subj), sd = sqrt(1 - icc))
data.frame(subj = subj, group = grp[subj], y = y)
}
# One-way ANOVA F, written out so the permutation loop is fast
AnovaF <- function(y, g) {
g <- as.integer(factor(g)); n <- length(y); k <- max(g)
y <- y - mean(y)
totals <- as.numeric(rowsum(y, g)); counts <- tabulate(g)
ss_a <- sum(totals^2 / counts); ss_w <- sum(y^2) - ss_a
(ss_a / (k - 1)) / (ss_w / (n - k))
}
PermTest <- function(d, B = 199, unit = c("observation", "subject")) {
unit <- match.arg(unit)
F_obs <- AnovaF(d$y, d$group)
subj_group <- tapply(d$group, d$subj, `[`, 1) # one label per subject, named by subject
F_perm <- replicate(B, {
g_perm <- if (unit == "observation") {
sample(d$group)
} else {
perm <- setNames(sample(subj_group), names(subj_group)) # shuffle labels across subjects
perm[as.character(d$subj)] # carry each subject's visits along
}
AnovaF(d$y, g_perm)
})
(sum(F_perm >= F_obs) + 1) / (B + 1)
}
# Prespecified rejection fractions; the group-linked case lacks the same null justification.
Level <- function(visits, unit, design = "independent", R = 200) {
p <- replicate(R, PermTest(SimulateVisits(n_subj = 20, visits = visits,
icc = 0.5, design = design), B = 99, unit = unit))
rate <- mean(p <= 0.05)
c(rejection = rate, mcse = sqrt(rate * (1 - rate) / R))
}
designs <- list(`4 visits each` = list(visits = 4, design = "independent"),
`2 or 6; assignment independent` = list(visits = c(2, 6), design = "independent"),
`2 vs 6; group-linked records` = list(visits = c(2, 6), design = "group-linked"))
grid <- expand.grid(design = names(designs), unit = c("observation", "subject"),
stringsAsFactors = FALSE)
res <- t(vapply(seq_len(nrow(grid)), function(i) {
cfg <- designs[[grid$design[i]]]
Level(cfg$visits, grid$unit[i], cfg$design)
}, numeric(2)))
comparison <- cbind(grid, as.data.frame(res))
comparison design unit rejection mcse
1 4 visits each observation 0.215 0.02904953
2 2 or 6; assignment independent observation 0.245 0.03041176
3 2 vs 6; group-linked records observation 0.180 0.02716616
4 4 visits each subject 0.035 0.01299519
5 2 or 6; assignment independent subject 0.040 0.01385641
6 2 vs 6; group-linked records subject 0.030 0.01206234For the first two designs, person assignment is independent of outcomes and visit schedules, and the generated outcomes have no treatment effect. Relabeling people reproduces the person-level randomization under its sharp null, even with unequal prespecified visit counts. Shuffling individual profiles does not reproduce that assignment and can yield a high rejection fraction under positive intraclass correlation.
The third design is deliberately different: record length depends on original group. Equal outcome means do not make those observed records exchangeable. Neither whole-person relabeling nor a conservative-looking empirical rejection fraction supplies the missing justification. The distinction is the assignment/observation process, not merely balance of the final row counts.
In the code, grp assigns one label per person and M gives that person’s record length. Under the independent schedule, a six-visit record could have received either group label. Under the group-linked schedule, a six-visit record identifies the original group. Relabeling it as a two-visit-group record while keeping all six visits does not recreate the observation process under a new assignment. This is why preserving clusters is necessary but does not, by itself, prove exchangeability.
Code
par(mar = c(5, 4, 1, 1))
plot(NA, xlim = c(0.6, 3.4), ylim = c(0, max(0.35, comparison$rejection + 2*comparison$mcse)),
xaxt = "n", xlab = "Visit schedule", ylab = "Rejection fraction at nominal 0.05")
axis(1, at = 1:3, labels = c("Balanced", "Unequal;\nindependent", "Group-linked\nrecords"), cex.axis = 0.8)
for (j in 1:2) {
unit <- c("observation", "subject")[j]
a <- comparison[comparison$unit == unit, ]
x <- 1:3 + c(-0.07, 0.07)[j]
color <- c("#D55E00", "#2F80C8")[j]
points(x, a$rejection, pch = 15 + j, col = color)
arrows(x, pmax(0, a$rejection - 2*a$mcse), x, a$rejection + 2*a$mcse,
angle = 90, code = 3, length = 0.04, col = color)
}
abline(h = 0.05, lty = 2)
legend("topright", c("Profile shuffle", "Person shuffle"), col = c("#D55E00", "#2F80C8"),
pch = 16:17, bty = "n", cex = 0.8)
The accompanying table retains the numerical results; colors distinguish observation (orange) and person (blue) schemes. The Monte Carlo uncertainty is substantial at 200 replicates, so do not attach a universal numerical rejection rate to a design label.
For a within-person factor, within-person swaps need their own randomization or symmetry argument; chronological visits are not freely exchangeable merely because they belong to one person. Blocks can restrict a justified assignment scheme. The permute documentation distinguishes within-plot, whole-plot, and blocking controls. None substitutes for stating the scientific null.
4 Distance-based tests: PERMANOVA
Now suppose each observation is a high-dimensional vector (a gene expression profile, a microbiome composition) and you want to test whether groups differ. One route is to choose a distance between observations, compute the N \times N distance matrix, and ask whether within-group distances are smaller than between-group distances.
PERMANOVA (Anderson 2001) does this with an ANOVA-style decomposition that needs only distances. Here N is the number of profiles, G is the number of groups, and n_k is group k’s size; G does not denote the number of genes. With d_{ij} the distance between profiles i and j,
\begin{aligned} SS_T&=N^{-1}\sum_{i<j}d_{ij}^2,\\ SS_W&=\sum_k n_k^{-1}\sum_{i<j:\,i,j\in k}d_{ij}^2,\\ SS_A&=SS_T-SS_W, \end{aligned}
and the pseudo-F is F = \frac{SS_A/(G-1)}{SS_W/(N-G)}. For Euclidean distance, the coordinate-wise sums of squares add before the ratio is formed: this is not the sum of univariate F-statistics. Other dissimilarities define another geometry; negative principal-coordinate eigenvalues can complicate the centroid interpretation. Permutation calibration still needs a justified null and relabeling scheme, regardless of distance.
The Euclidean identity behind the calculation is
\frac1N\sum_{i<j}\|\mathbf Y_i-\mathbf Y_j\|^2 =\sum_i\|\mathbf Y_i-\bar{\mathbf Y}\|^2.
Thus SS_T measures total spread around the overall center, SS_W adds spread around each group’s center, and their difference SS_A measures the separation of group centers in this geometry. Dividing by G-1 and N-G gives an ANOVA-like ratio. The familiar shape of that ratio does not establish an ordinary parametric F reference; the example calibrates it by relabeling.
Code
PermanovaF <- function(D, g) {
D <- as.matrix(D); N <- nrow(D); groups <- unique(g); G <- length(groups)
ss <- function(idx) sum(D[idx, idx][lower.tri(D[idx, idx])]^2) / length(idx)
SS_T <- ss(seq_len(N))
SS_W <- sum(sapply(groups, function(k) ss(which(g == k))))
((SS_T - SS_W) / (G - 1)) / (SS_W / (N - G))
}
set.seed(3)
X <- matrix(rnorm(40 * 5), 40, 5); g <- rep(c("a", "b"), each = 20)
X[g == "b", ] <- X[g == "b", ] + 0.5
D <- dist(X) # Euclidean
F_obs <- PermanovaF(D, g)
F_perm <- replicate(999, PermanovaF(D, sample(g)))
c(F = F_obs, p = (sum(F_perm >= F_obs) + 1) / 1000) F p
4.938003 0.001000 The following executes the same one-factor analysis in vegan, using an explicit permutation matrix from permute. It also checks the coordinate sum-of-squares identity. These packages are required for this and the later package examples.
Code
meta <- data.frame(g = factor(g))
perms <- permute::shuffleSet(nrow(X), nset = 199)
ad <- vegan::adonis2(D ~ g, data = meta, permutations = perms)
ss_total <- sum(sweep(X, 2, colMeans(X))^2)
ss_within <- sum(vapply(split(seq_len(nrow(X)), g), function(i) {
sum(sweep(X[i, , drop = FALSE], 2, colMeans(X[i, , drop = FALSE]))^2)
}, numeric(1)))
f_coordinates <- ((ss_total - ss_within)/(nlevels(meta$g)-1)) /
(ss_within/(nrow(X)-nlevels(meta$g)))
stopifnot(isTRUE(all.equal(F_obs, f_coordinates)),
isTRUE(all.equal(F_obs, as.numeric(ad[1, "F"]))))
adPermutation test for adonis under reduced model
Permutation: free
Number of permutations: 199
vegan::adonis2(formula = D ~ g, data = meta, permutations = perms)
Df SumOfSqs R2 F Pr(>F)
Model 1 24.172 0.115 4.938 0.005 **
Residual 38 186.016 0.885
Total 39 210.189 1.000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The F-statistics agree. The earlier hand-coded run used 999 random relabelings and this package call uses a fresh set of 199, so their Monte Carlo p-values need not match and have different minimum grid values. A p-value comparison must use the same justified reference and sampling scheme.
For models with several terms, by=NULL, sequential terms, and marginal tests answer different questions. Adding covariates does not make every raw-label permutation exact. Record both the model and the permutation method; see the official adonis2 documentation.
4.1 What PERMANOVA actually tests
Equality of group distributions is a sufficient exchangeability null for the simple independent-sample label test above. In a randomized experiment a relevant sharp no-effect null can instead justify the assignment reference. Neither should be silently replaced by the weaker assertion that population centroids are equal. The pseudo-F is sensitive to differences in location (the groups’ centroids differ) and also to differences in dispersion (one group is more spread out). A significant PERMANOVA does not tell you which. With unbalanced groups the problem is worse: Anderson and Walsh (2013) show that rejection relative to an equal-centroid scientific null can be distorted by unequal dispersions, especially with unequal group sizes. Those populations do not satisfy the equal-distribution null, so a rejection is not automatically a false positive for that stronger null. Warton, Wright, and Wang (2012) make the related point that distance-based methods confound location and dispersion because the mean-variance relationship in count data makes the more abundant group look more variable.
PERMDISP (Anderson 2006) assesses group dispersions through distances to group centers. In vegan::betadisper, the default center is a spatial median; centroid and bias-adjustment choices should be recorded. It is a useful companion diagnostic/test with its own conditions. A non-significant PERMDISP does not establish equal dispersion or prove that a PERMANOVA result is solely a location difference. Interpret both tests with their power, the distance, the balance of the design, and plots or estimates of the relevant quantities.
Code
dispersion <- vegan::betadisper(D, meta$g, type = "median", bias.adjust = TRUE)
vegan::permutest(dispersion, permutations = 199)$tab Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.06956812 0.06956812 0.1564657 199 0.705
Residuals 38 16.89563867 0.44462207 NA NA NACode
set.seed(7)
visual <- list(
`Population location shift` = rbind(matrix(rnorm(60), 30, 2),
sweep(matrix(rnorm(60), 30, 2), 2, c(1.5, 0), "+")),
`Population dispersion change` = rbind(matrix(rnorm(60), 30, 2),
matrix(rnorm(60, sd = 2), 30, 2)))
limits <- range(unlist(visual)); labels <- rep(0:1, each = 30)
par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
for (name in names(visual)) {
a <- visual[[name]]
plot(a, pch = 16, cex = 0.6, col = c("#D55E0088", "#2F80C888")[labels+1],
xlim = limits, ylim = limits, xlab = "Feature 1", ylab = "Feature 2", main = name, asp = 1)
centers <- rbind(colMeans(a[labels == 0, ]), colMeans(a[labels == 1, ]))
points(centers, pch = 4, cex = 1.5, lwd = 2, col = c("#D55E00", "#2F80C8"))
}
4.2 Choice of distance
The distance defines what “different” means, and the choice is part of the hypothesis.
- Euclidean on log-transformed expression: differences in every gene count equally; dominated by high-variance genes unless scaled.
- Bray-Curtis on compositions or counts: ignores joint absences, bounded in [0, 1], standard in microbiome work; it is not a metric (fails the triangle inequality), which matters for some downstream methods.
- Correlation distance 1 - r between profiles: ignores overall level, sensitive to shape.
- Aitchison (Euclidean after centered log-ratio): a compositional geometry for strictly positive components; zeros require an explicit measurement-aware treatment that can affect the analysis.
Any of these distances can be used within a permutation statistic, but validity still depends on the null and allowed relabelings. What differs is the alternative the test is powered against and how location and dispersion are confounded.
For example, two nonconstant profiles have correlation distance zero if one is a positive affine transformation of the other, even if their overall levels and magnitudes differ greatly. Euclidean distance retains such changes. Scaling each feature to unit variance before Euclidean analysis changes the relative contribution of features and therefore changes the comparison. Choose the geometry from the scientific question, rather than treating distance as an interchangeable software option.
If preprocessing uses group labels, such as selecting genes for a large observed group contrast, repeat that selection within every permutation when it is part of the tested procedure. Holding a label-selected feature set fixed gives the permuted analyses an advantage the observed analysis did not have to meet. A fixed, label-invariant transformation can have a different justification. Finally, an omnibus rejection followed by feature-wise claims creates a further testing family; one significant global test does not automatically validate every follow-up claim.
5 Repeated measures in distance-based tests
A person-level factor needs an assignment or exchangeability argument at the person level. strata=subject_id restricts row permutations within people; it does not move whole people between treatment groups. With a constant group label within person it is not the generic between-person test desired here.
For balanced records, Plots with free whole-plot permutations and Within(type="none") can preserve visit order while moving whole people. This executable example has 12 people, three visits each, and randomly assigned person-level groups. Under the appropriate sharp no-effect null, the whole-person reference reproduces the assignment.
Code
set.seed(8)
subject_id <- rep(seq_len(12), each = 3); visit <- rep(1:3, 12)
assigned <- sample(rep(0:1, each = 6))
repeated_meta <- data.frame(subject_id, visit, treatment = assigned[subject_id])
u <- matrix(rnorm(12*3), 12, 3)
profiles <- u[subject_id, ] + matrix(rnorm(36*3), 36, 3)
control <- permute::how(plots = permute::Plots(strata = factor(subject_id), type = "free"),
within = permute::Within(type = "none"))
whole_people <- permute::shuffleSet(36, nset = 199, control = control)
cbind(original_person = subject_id[1:9],
permuted_person = subject_id[whole_people[1, 1:9]],
original_visit = visit[1:9], permuted_visit = visit[whole_people[1, 1:9]]) original_person permuted_person original_visit permuted_visit
[1,] 1 1 1 1
[2,] 1 1 2 2
[3,] 1 1 3 3
[4,] 2 10 1 1
[5,] 2 10 2 2
[6,] 2 10 3 3
[7,] 3 9 1 1
[8,] 3 9 2 2
[9,] 3 9 3 3Code
vegan::adonis2(dist(profiles) ~ treatment, data = repeated_meta, permutations = whole_people)Permutation test for adonis under reduced model
Plots: factor(subject_id), plot permutation: free
Permutation: none
Number of permutations: 199
vegan::adonis2(formula = dist(profiles) ~ treatment, data = repeated_meta, permutations = whole_people)
Df SumOfSqs R2 F Pr(>F)
Model 1 9.223 0.04475 1.5927 0.57
Residual 34 196.875 0.95525
Total 35 206.097 1.00000 Whole-plot row permutations require compatible record structures. With unequal visit counts, do not force these balanced row blocks onto the data. A custom reassignment of person-level labels with statistic recomputation, as in PermTest(), can reproduce a known assignment when visits were prespecified independently of treatment. Informative follow-up or an observational disease grouping needs additional assumptions or a different analysis. See Units of analysis, replication, and study design.
6 Exercises
Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.
Checkpoint · MV-01 — Match the permutation to the design and null
Task type: Diagnosis and interpretation.
Twelve people are independently enrolled in an experiment and randomly assigned, six per treatment. Each contributes three repeated multivariate profiles. You want a treatment comparison. Describe a design-valid reassignment of treatment labels, explain why arbitrary profile-level relabeling is inappropriate, and say whether a significant PERMANOVA plus non-significant PERMDISP proves a pure centroid difference.
Ready to continue when: you can identify the randomization unit, preserve its dependent measurements, and limit the interpretation of an omnibus result.
Transfer: Why would subject-level exchangeability require additional justification in an observational disease comparison?
Hint
Reproduce the treatment assignment at the person level. Failure to reject a dispersion test does not establish equal dispersion.
Solution outline — open after trying
Under the appropriate randomization null, reassign the six treatment labels among people and carry each person’s profiles together. Arbitrary relabeling of the 36 profiles does not reproduce the experiment and breaks the within-person structure. Define the statistic and null before calibrating it; exact randomization inference uses the relevant sharp null, while weaker nulls need their own justification. A non-significant PERMDISP may reflect low power and does not isolate a centroid-only effect. The distance, statistic, design, and distributional alternatives still matter.
Additional practice (optional)
These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.
- Modify the first example so that the groups have equal means but standard deviations 1 and 3, and estimate the level of the permutation test of the mean difference over 1000 simulated datasets. Then repeat with unequal group sizes (10 and 40), both ways round. Relate what you see to exchangeability.
- Run the clustered-data example with
icc = 0and withicc = 0.9. How does the inflation of the observation-level test depend on the intraclass correlation and on the number of visits? - Verify numerically that
PermanovaF()with Euclidean distance equals the sum over variables of the univariate ANOVA sums of squares, decomposed the same way. - Simulate two groups with identical centroids and different dispersions (multiply one group’s deviations by 2), unbalanced 15 versus 45, and estimate PERMANOVA’s rejection rate both ways round. This is the Anderson and Walsh result; see it yourself.
7 Reading
The official adonis2, betadisper, and permute vignette document the executable package examples.
- Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26(1), 32 to 46.
- Anderson, M. J. and ter Braak, C. J. F. (2003). Permutation tests for multi-factorial analysis of variance. Journal of Statistical Computation and Simulation, 73(2), 85 to 113. The reference for which units to permute.
- Anderson, M. J. and Walsh, D. C. I. (2013). PERMANOVA, ANOSIM, and the Mantel test in the face of heterogeneous dispersions: what null hypothesis are you testing? Ecological Monographs, 83(4), 557 to 574.
- Warton, D. I., Wright, S. T., and Wang, Y. (2012). Distance-based multivariate analyses confound location and dispersion effects. Methods in Ecology and Evolution, 3(1), 89 to 101.
- Phipson, B. and Smyth, G. K. (2010). Permutation p-values should never be zero. Statistical Applications in Genetics and Molecular Biology, 9(1), Article 39.
- Ernst, M. D. (2004). Permutation methods: a basis for exact inference. Statistical Science, 19(4), 676 to 685. A short, careful general introduction.