Overview

♦️ Standard PERMANOVA (via adonis2) is the workhorse for microbiome association testing. But it has two well-documented blind spots that require careful consideration:

Problem What goes wrong What permaBrioche does
Invalid permutations under repeated-measures designs Excess false positives (small p-values) for subject-level covariates Confounding-aware block permutations
Inflated R² under null, even without repeated-measures R² > 0 even under the null (no true association) Null-centered R² or Hájek effect size

♦️ This demo walks through both problems live, using HMP2 microbiome data bundled with the package.

Installation

devtools::install_github("biobakery/permabrioche", force = TRUE)

GitHub: https://github.com/biobakery/permabrioche
Questions / issues: https://github.com/biobakery/permabrioche/issues

Setup

library(permabrioche)   # our package
library(vegan)          # for vegdist; adonis2 comparison
library(ggplot2)        # for plotting
set.seed(1234)

♦️ We’ll use data derived from the Human Microbiome Project 2 (HMP2), a longitudinal study following individuals over time with repeated stool sampling. It’s already saved in the package - we just have to load it:

d_name    <- load(system.file("extdata", "d_hmp.Rda",    package = "permabrioche"))
bugs_name <- load(system.file("extdata", "bugs_hmp.Rda", package = "permabrioche"))

meta <- get(d_name)
bugs <- get(bugs_name)

# Align row names
rownames(meta) <- meta$sample
stopifnot(identical(rownames(meta), rownames(bugs)))

# Downsample to 25 subjects for live demo
keep_subjects <- sample(unique(meta$subject), 25)
meta <- meta[meta$subject %in% keep_subjects, ]
bugs <- bugs[rownames(meta), ]

♦️ We have 343 samples from 25 subjects, each containing abundances of 100 taxa :

cat("Samples:", nrow(meta), "\n")
#> Samples: 343
cat("Subjects:", length(unique(meta$subject)), "\n")
#> Subjects: 25
cat("Taxa:", ncol(bugs), "\n")
#> Taxa: 100

Data structure

♦️ Let’s use the ‘head’ function to double check that each subject appears multiple times (longitudinal samples). We will show that using standard adonis2 on this type of data can be problematic (and what you should do instead):

knitr::kable(head(meta[, c("sample", "subject")], 8),
             caption = "First few rows of sample-level metadata")
First few rows of sample-level metadata
sample subject
CSM5FZ4E_P CSM5FZ4E_P C3003
CSM5FZ4G_P CSM5FZ4G_P C3003
CSM5FZ4K_P CSM5FZ4K_P C3003
CSM5FZ4M CSM5FZ4M C3003
CSM5LLGB_P CSM5LLGB_P M2014
CSM5MCV1_P CSM5MCV1_P C3007
CSM5MCV5_P CSM5MCV5_P C3007
CSM5MCVB_P CSM5MCVB_P C3007 
colnames(bugs) <- paste0("bug", seq_len(ncol(bugs)))   # tidy display
knitr::kable(head(bugs[, 1:6]), caption = "First few rows of the abundance data")
First few rows of the abundance data
bug1 bug2 bug3 bug4 bug5 bug6
CSM5FZ4E_P 3.15531 0.09501 0.15372 39.43736 16.52267 0.00000
CSM5FZ4G_P 14.81362 0.40161 0.00446 23.65032 11.47776 0.00000
CSM5FZ4K_P 9.92632 0.86028 0.00504 25.66730 7.92792 0.00000
CSM5FZ4M 5.78900 0.22401 0.00000 27.08071 12.07513 0.00000
CSM5LLGB_P 2.85017 4.62317 55.70357 0.11023 0.01096 5.65554
CSM5MCV1_P 0.00000 0.00000 0.00845 5.52779 0.00000 0.80608

Create the distance matrix

♦️ One of the main ingredients in a PERMANOVA is the distance matrix. Just like normal, we can use ‘vegdist’ from the ‘vegan’ package to get our distance matrix:

dist_bc <- vegdist(bugs, method = "bray")
Visualized example of microbial abundances across two groups

Visualized example of microbial abundances across two groups

♦️ PERMANOVA answers the question of “are the abundances of these two groups actually different or is it just due to chance?”

Problem 1 — Repeat Measures → Invalid Permutations → False Positives

♦️ Now we will show that if you have repeated measures in your data, the p-values from adonis2 will be inflated with false positives. This means that these p-values are misleadingly low and are not statistically valid.

❇️ Question: What is a “repeated measure”?

💠 GOAL: Imagine we want to know whether delivery mode is associated with gut microbiome composition.

♦️ We simulate delivery mode here by assigning a random binary label at the subject level — so by construction there is no true association:

set.seed(1234)
# Assign a random binary exposure CONSTANT within subject — true effect = 0
subjects    <- unique(meta$subject)
subj_labels <- setNames(
  sample(c(0, 1), length(subjects), replace = TRUE),
  subjects
)
meta$exposure_invariant <- subj_labels[meta$subject]

cat("Exposure is invariant within subject:\n")
#> Exposure is invariant within subject:
table(meta$subject, meta$exposure_invariant)[1:5, ]
#>        
#>          0  1
#>   C3003  0 10
#>   C3007  0  3
#>   C3013  0 22
#>   C3016 20  0
#>   C3024  1  0

♦️From the following illustration, you can see that our data has repeated measures. Within the two groups, there are multiple samples from the same subjects.

What standard adonis2 does

♦️ Standard adonis2 permutes sample labels without regard to which samples belong to the same subject. Because subjects cluster strongly in microbiome space, the null distribution is artificially compressed — incorrectly yielding far too many small p-values.

♦️ Let’s see that in action:

set.seed(1234)
# ---- adonis2: unadjusted (incorrect for repeated measures) ----
res_adonis_unadj <- adonis2(
  dist_bc ~ exposure_invariant,
  data         = meta,
  permutations = 999
)
cat("=== adonis2 (unadjusted) ===\n")
#> === adonis2 (unadjusted) ===
print(res_adonis_unadj)
#> Permutation test for adonis under reduced model
#> Permutation: free
#> Number of permutations: 999
#> 
#> adonis2(formula = dist_bc ~ exposure_invariant, data = meta, permutations = 999)
#>           Df SumOfSqs      R2      F Pr(>F)    
#> Model      1    4.021 0.04356 15.532  0.001 ***
#> Residual 341   88.286 0.95644                  
#> Total    342   92.307 1.00000                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(1234)
# ---- adonis2: adding subject as a covariate (still wrong) ----
res_adonis_varsubj <- adonis2(
  dist_bc ~ subject + exposure_invariant,
  data         = meta,
  permutations = 999,
  by = "terms"
)
cat("=== adonis2 (exposure + subject) ===\n")
#> === adonis2 (exposure + subject) ===
print(res_adonis_varsubj)
#> Permutation test for adonis under reduced model
#> Terms added sequentially (first to last)
#> Permutation: free
#> Number of permutations: 999
#> 
#> adonis2(formula = dist_bc ~ subject + exposure_invariant, data = meta, permutations = 999, by = "terms")
#>           Df SumOfSqs      R2      F Pr(>F)    
#> subject   24   63.108 0.68368 28.638  0.001 ***
#> Residual 318   29.199 0.31632                  
#> Total    342   92.307 1.00000                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Key point: Even when you add subject to the formula, adonis2 still permutes labels across subjects and actually drops invariant variables entirely.

♦️ This p-vaue calculation is an issue for both invariant covariates (ones that never change across a subject’s repeated samples) and variant covariates (ones that do change across a subject’s repeated samples).

Proposed Solution

What permaBrioche does to correct this:

♦️PERMANOVA_repeat_measures() permutes invariant covariates across subjects (not across samples), so the permutation scheme respects the repeated-measures structure. For variant covariates, the permutation is done within subjects. You do not have to specify which one exists in your data - the function automatically handles it for you 🙂.

set.seed(1234)
res_perm <- PERMANOVA_repeat_measures(
  formula           = dist_bc ~ exposure_invariant,
  data              = meta,
  sample_id         = "sample",
  blocking_variable = "subject",
  permutations      = 999
)

cat("=== permaBrioche: PERMANOVA_repeat_measures ===\n")
#> === permaBrioche: PERMANOVA_repeat_measures ===
print(res_perm)
#> Permutation test for adonis under reduced model
#> Overall (omnibus) test of all terms jointly
#> Permutation: blocked by subject
#> Number of permutations: 999
#> PERMANOVA_repeat_measures(formula = dist_bc ~ exposure_invariant, 
#>     data = meta, sample_id = "sample", blocking_variable = "subject", 
#>     permutations = 999)
#>           Df SumOfSqs      R2      F Pr(>F)
#> Model      1    3.693 0.04001 14.212   0.32
#> Residual 341   88.614 0.95999              
#> Total    342   92.307 1.00000

Under the null (no true effect), a valid method should give p ≈ Uniform(0,1). permaBrioche achieves this while a standard run of adonis2 does not.

Side-by-side comparison

Comparison under the null (true effect = 0)
Method p-value Valid
adonis2 (unadjusted) 0.001 ❌ No
adonis2 (+ subject covariate) 0.001 ❌ No
permaBrioche: PERMANOVA_repeat_measures 0.320 ✅ Yes
♦️ PERMANOVA_repeat_measures() can also handle formulas with multiple covariates of the same type on the right-hand side. For example, if all covariates are invariant within subjects (e.g., sex, treatment arm, genotype), the function performs the appropriate between-subject permutations for each term. Likewise, if all covariates vary within subjects (e.g., diet, medication use, time-varying biomarkers), it automatically uses within-subject permutations. No additional arguments are required. However it does not support formulas with multiple covariates of differing types (e.g. sex and medication use).

Problem 2 — Non-centered R²

Why is it common to see large R² values in microbiome research - even when the true effect is minimal?

♦️ Even when there is actually no association between a covariate and the microbiome (no-association null hypothesis), PERMANOVA R² follows a Beta distribution centered above zero:

\[E_0[R^2] = \frac{g - 1}{N - 1}\]

where g = number of levels of the covariate of interest and N = number of samples. With small N this can be substantial.

♦️ So as an example, if your data has 30 samples and your covariate has 4 levels, even though there is no true effect, your R² will be centered around 0.103! 😯:

N <- 30 # 30 samples
g <- 4    # exposure with four levels
expected_null_R2 <- (g - 1) / (N - 1)
cat(sprintf(
  "For example, with N=%d samples and g=%d groups, E[R² | null] = %.3f\n",
  N, g, expected_null_R2
))
#> For example, with N=30 samples and g=4 groups, E[R² | null] = 0.103

❇️ Question: In microbiome data, we often have small sample sizes. As sample size (N) gets smaller, does this deviation increase or decrease?



E[R² | null] = (g−1)/(N−1) = 0.103

♦️ In our data, with 343 samples and 2 levels in our covariate, the deviation from 0 is quite small but still meaningful. For small sample sizes and covariates with multiple levels, this deviation increases.

set.seed(1234)
# Simulate 200 null R² values using random label permutations
null_r2 <- replicate(200, {
  shuffled_exposure <- sample(meta$exposure_invariant)
  meta$shuffled_exposure <- shuffled_exposure
  fit <- adonis2(dist_bc ~ shuffled_exposure, data = meta, permutations = 0)
  fit["Model", "R2"]
})

# In our case, we can calculate the expected R^2 under the null as
g <- 2
N <- nrow(meta)
expected_null_R2 <- (g - 1) / (N - 1)

ggplot(data.frame(R2 = null_r2), aes(x = R2)) +
  geom_histogram(bins = 30, fill = "#4E79A7", colour = "white", alpha = 0.85) +
  geom_vline(xintercept = expected_null_R2, colour = "#E15759", linewidth = 1.2,
             linetype = "dashed") +
  geom_vline(xintercept = 0, colour = "black", linewidth = 0.8) +
  annotate("text", x = expected_null_R2 + 0.003, y = Inf,
           label = sprintf("E[R²|null] = %.3f", expected_null_R2),
           hjust = 0, vjust = 2, colour = "#E15759", size = 4) +
  labs(
    title    = "Simulated null distribution of PERMANOVA R²",
    subtitle = "Labels randomly permuted 200 times — true effect is zero",
    x        = "R²", y        = "Count"
  ) +
  theme_minimal(base_size = 13)

The distribution is entirely above zero even though there is no real effect. Raw R² can be quite misleading as an effect size.

Proposed Solution 1 — Null-centered R²

♦️PERMANOVA_repeat_measures() with center_R2 = TRUE subtracts the mean null R² estimated from the permutation distribution, giving an effect size centered at zero under the null:

\[\tilde{R}^2 = R^2 - \frac{g-1}{N-1}\]

set.seed(1234)
res_perm_centered <- PERMANOVA_repeat_measures(
  formula           = dist_bc ~ exposure_invariant,
  data              = meta,
  sample_id         = "sample",
  blocking_variable = "subject",
  permutations      = 999,
  center_R2         = TRUE
)

cat("=== permaBrioche: null-centered R² ===\n")
#> === permaBrioche: null-centered R² ===
print(res_perm_centered)
#> Permutation test for adonis under reduced model
#> Overall (omnibus) test of all terms jointly
#> Permutation: blocked by subject
#> Number of permutations: 999
#> PERMANOVA_repeat_measures(formula = dist_bc ~ exposure_invariant, 
#>     data = meta, sample_id = "sample", blocking_variable = "subject", 
#>     permutations = 999, center_R2 = TRUE)
#>           Df SumOfSqs      R2      F Pr(>F) R2_centered
#> Model      1    3.693 0.04001 14.212   0.32   0.0031862
#> Residual 341   88.614 0.95999                          
#> Total    342   92.307 1.00000

Interpretation: R2_centered ≈ 0 confirms no real effect beyond what chance partitioning would produce. Importantly, the p-value is unchanged.

Proposed Solution 2 — Hájek Distance-Based Effect Size

A more interpretable effect size

💠 Using the null-centered R² helps us get a more honest value, but the interpretation is still clunky: “The covariate explains roughly X percentage points more of the variance in community dissimilarity than would be expected by chance.”

There is nothing wrong with using the null-centered R², but it could be nice to have a more concretely interpretable effect size. This is why we introduce the “Hájek estimator”.

♦️The Hájek estimator answers:

On average, how far does exposure move samples away from the control centroid (center), in units of the chosen dissimilarity metric?

♦️This is directly interpretable — e.g. “exposure increases Bray-Curtis distance from the control centroid by 0.05 units” — unlike the abstract variance-explained framing of R².

♦️ To use this estimator, we can simply call the following function:

hajek_repeat_measures()

set.seed(1234)
res_hajek <- hajek_repeat_measures(
  formula           = dist_bc ~ exposure_invariant,
  data              = meta,
  bugs              = bugs,
  sample_id         = "sample",
  blocking_variable = "subject",
  covariate_name    = "exposure_invariant",
  permutations      = 999,
  method            = "bray"
)

cat("=== Hájek effect size ===\n")
#> === Hájek effect size ===
cat("Observed τ:", round(res_hajek$observed, 4), "\n")
#> Observed τ: 0.0218
cat("p-value:   ", round(res_hajek$pval, 4), "\n")
#> p-value:    0.8188

Interpretation: - τ ≈ 0 → exposure does not move samples away from the control centroid - τ > 0 → exposure moves samples farther from the control group - τ < 0 → exposure moves samples closer to the control group

Because our exposure was random, we expect τ ≈ 0 and p ≈ non-significant.

Location–Dispersion Decomposition

Separating where from how spread out

♦️ One interesting feature about this estimand is that in the Euclidean distance setting, the Hájek effect decomposes as:

\[\tau = \tau_{\text{location}} + \tau_{\text{dispersion}}\]

  • τ_location: did the group centroid shift? (systematic change in mean composition)
  • τ_dispersion: did within-group variability change? (increased heterogeneity)

♦️ This is useful because two studies can have the same total τ for completely different biological reasons.

hajek_repeat_measures_loc_and_disp()

♦️ The final function in permaBrioche allows us to calculate this decomposition:

This function requires data at two levels — sample level and subject level:

set.seed(1234)
# --- Sample-level blocking vector ---
blocks <- meta$subject
names(blocks) <- meta$sample
blocks <- blocks[rownames(bugs)]

# --- Sample-level within-subject permutation frame (empty = nothing varies within subject) ---
permute_within <- data.frame(row.names = rownames(bugs))

# --- Subject-level data ---
block_data <- aggregate(
  exposure_invariant ~ subject,
  data = meta,
  FUN  = function(x) x[1]
)
rownames(block_data) <- block_data$subject
block_data$subject   <- NULL

♦️ And we put these ingredients into the function in the following manner:

set.seed(1234)
res_ld <- hajek_repeat_measures_loc_and_disp(
  D              = dist_bc,
  permute_within = permute_within,
  blocks         = blocks,
  block_data     = block_data,
  bugs           = bugs,
  covariate_name = "exposure_invariant"
)

cat("=== Location–Dispersion Decomposition ===\n")
#> === Location–Dispersion Decomposition ===
print(res_ld)
#> $observed_tau
#> [1] -1368.944
#> 
#> $observed_loc
#> [1] -987.3894
#> 
#> $observed_disp
#> [1] -409.7602

♦️ We can use ggplot to nicely format and display this decomposition:

set.seed(1234)
# Extract components for a simple visual summary
components <- data.frame(
  Component = c("Total (τ)", "Location (τ_loc)", "Dispersion (τ_disp)"),
  Value     = c(
    res_ld$observed_tau,
    res_ld$observed_loc,
    res_ld$observed_disp
  )
)

ggplot(components, aes(x = Component, y = Value, fill = Component)) +
  geom_col(width = 0.5, show.legend = FALSE) +
  geom_hline(yintercept = 0, linewidth = 0.8) +
  scale_fill_manual(values = c("#4E79A7", "#59A14F", "#F28E2B")) +
  labs(
    title    = "Hájek effect size decomposition",
    subtitle = "Location = centroid shift · Dispersion = variability change",
    y        = "Effect size (Euclidean units)", x        = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(size = 11))

❇️ Question: Which of the following two outputs would indicate a systematic compositional shift rather than just indicating that one group is more heterogeneous?

Large location, small dispersion

Large dispersion, small location

Complete Function Summary

Function What it does Key argument
PERMANOVA_repeat_measures() Valid hypothesis test blocking_variable
PERMANOVA_repeat_measures(..., center_R2 = TRUE) + Null-centered R² center_R2 = TRUE
hajek_repeat_measures() Interpretable distance effect size covariate_name, method
hajek_repeat_measures_loc_and_disp() Decomposes effect into location + dispersion block_data, blocks

When to Use Which

Do you have longitudinal / clustered microbiome data?
│
├─ Is your covariate INVARIANT within subject? (the function will figure 
this out for you automatically)
│   (e.g. sex, delivery mode, genotype, treatment arm)
│   └─→ PERMANOVA_repeat_measures()   ← between-block permutations
│
└─ Is your covariate VARIANT within subject? (the function will figure 
this out for you automatically)
    (e.g. diet, medication use)
    └─→ PERMANOVA_repeat_measures()   ← within-block permutations

Do you want an interpretable effect size?
└─→ hajek_repeat_measures()

Do you want to understand HOW the microbiome changed?
└─→ hajek_repeat_measures_loc_and_disp()   (Euclidean distance)

Do you want a null-corrected variance summary?
└─→ PERMANOVA_repeat_measures(..., center_R2 = TRUE)

Take-Home Messages

🥇 1. adonis2 alone is not valid for repeated-measures data. Even adding subject as a covariate does not fix the permutation scheme for invariant covariates — false positive rates can exceed 70% on HMP2 data.

🥈 2. PERMANOVA R² is upwardly inflated under the null. A significant positive R² does not necessarily imply a real effect. Use null-centered R² or the Hájek estimator instead.

🥉️ 3. The Hájek estimator gives you a geometrically interpretable effect size in the units of your chosen dissimilarity metric — something R² does not provide, and the location–dispersion decomposition tells you whether a covariate shifts the centroid, increases heterogeneity, or both.

💠 Overall: When running PERMANOVA on any microbiome dataset, make sure the permutation scheme respects its structure, and consider pairing your R² with a null-centered effect size. 😄