library(ggplot2)
library(MASS)
library(vegan)
te_paper <- "#f5f4ee"
te_ink <- "#16241d"
te_body <- "#2c3a31"
te_forest <- "#275139"
te_rust <- "#b5534e"
te_gold <- "#c9b458"
te_line <- "#dad9ca"
theme_datasheet <- function() {
theme_minimal(base_size = 12) +
theme(plot.background = element_rect(fill = te_paper, colour = NA),
panel.background = element_rect(fill = te_paper, colour = NA),
panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
panel.grid.minor = element_blank(),
text = element_text(colour = te_body),
plot.title = element_text(colour = te_ink, face = "bold"),
plot.subtitle = element_text(colour = te_body),
axis.text = element_text(colour = te_body))
}Pseudocounts turn read depth into a result
R
vegan
PERMANOVA
microbial ecology
compositional data
ecology tutorial
A fixed 0.5 pseudocount makes a sparse taxon’s log-ratio a function of library size, so unequal sequencing depth is reported as a group effect. Worked in R.
Corrected 14 October 2026: the post put the excess of the Wilcoxon test on proportions down to ties, when it comes from the shallow libraries missing the taxon more often; it put the excess of the two-part test down to its presence part and said the combination worked against it, when most of it comes from the Welch part on the detected samples and in the middle of the sweep the combined test rejects more often than either part, and the hurdle model it suggests now needs depth in both of its parts; it said a threefold depth difference was enough for every result above, when only the single-taxon sections ran at that ratio; and it said the ordination result depended entirely on sparsity, which its three settings cannot separate from richness and spread.
Corrected 7 October 2026: the post called adding a half and renormalising the Bayesian-multiplicative replacement, which instead imputes only the zeros and is not covered by the argument; it said an absent taxon’s log-ratio gap is exactly minus the log depth ratio, when the taxon’s own term in the centring mean makes it slightly smaller; and it said the negative binomial model held its rate, when its test runs slightly above nominal on the best-observed taxon at either depth, although depth does not mislead it.
Corrected 1 October 2026: the depth shift of a sparse taxon was computed by putting the expected count in place of the count, which overstated it, and is now its expectation under Poisson sampling. The summaries of the two-part test, the depth-scaled pseudocount, the zero-free community analysis, the ordination of real tables and the dispersion test no longer claim more than the results or Anderson and Walsh show, the case McMurdie and Holmes made against rarefying is no longer reduced to power, the Monte Carlo standard errors now cover rates far from the nominal one, and a second summary of the zeros post now matches the first.
Corrected 27 September 2026: the summary of the zeros post now describes a zero read as a sampling zero or a value below a detection limit, not as a detection-limit value in every case.
A soil survey sequences the 16S gene in two batches. The first batch went to the machine in March and came back with a median of about seventeen thousand reads per sample; the second went in June, after the facility changed its loading protocol, and came back with about fifty-two thousand. The two batches happen to correspond to the two treatments, because that is how the field season ran. Every analyst in the room knows this is a confound. The question is what to do about it, and the usual answer is a compositional one: divide the depth out on the log scale with a centred log-ratio, which is scale invariant and does not throw away reads. Zeros cannot be logged, so add 0.5 to every count first.
That last sentence is where the trouble is. The detection-limit post showed that a constant substituted at a moving limit turns the laboratory into a trend, and the zeros post said a zero read is a sampling zero or a value below a detection limit, not necessarily a true absence, and in a sequencing library either one is a non-detect at a limit set by depth; this post is where the two meet the advice of the metabarcoding check: the log-ratio that check ranked at the top reports the sequencer on any taxon sparse enough to have zeros, and the rarefaction it called wasteful is the repair.
The three posts being set against each other are worth naming exactly. Checking a metabarcoding analysis ran this same confounded design, with libraries of seventeen thousand against fifty-two thousand reads, and found that raw counts rejected a true null in 0.973 of runs while rarefied counts, proportions and the centred log-ratio held at 0.060, 0.062 and 0.055, within a Monte Carlo standard error of one another. It recommended dividing by the total, which keeps 0.743 power against the rarefied 0.543, and it concluded that rarefying “is wasteful, and the waste is measurable in power”. Its focal taxon had a share of 0.0006 of the reads, which at those depths is about ten reads a sample and almost never a zero. Compositional analysis and the zero problem showed that the log-ratio of a taxon that was zero tracks whatever value is imputed for it, and named the zero as a sampling zero or a value below a detection limit. Values below the detection limit in R showed what happens when that limit moves during a programme: a flat population acquires a decline.
Put those together and the prediction is specific. The detection limit for a sequencing library is of order one read in N, so it moves with depth, and the pseudocount is a constant substituted at it. That should make the log-ratio of any sparse taxon a function of log N, and in a design where depth is confounded with group it should be reported as a group effect. The metabarcoding check did not see this because its focal taxon was not sparse. Most of a real amplicon table is.
This post measures four things. First the size of the bias, as a function of a taxon’s expected count, which is arithmetic and not simulation. Then the false positive rate of seven analyses on a null taxon as its zero share is swept from about 0.03 to about 0.9. Then the same bias acting on a whole table at once, through a PERMANOVA on the Aitchison distance under a null grouping, in the regime Weiss and colleagues warned about in 2017, with a measurement of which taxa carry it. Then whether an ordination would have shown the analyst any of this, which turns out to depend on how sparse the table is.
The bias has a closed form, and it is about small counts
Write the centred log-ratio of taxon k in a sample with library size N as the log of its padded count minus the mean of the padded logs. If most taxa in the sample are abundant enough that their counts swamp the pseudocount, that mean is close to log N plus a term that does not depend on N. So the centred log-ratio of one taxon with true share p, at library size N, is approximately the log of its count plus the pseudocount, minus log N, minus a constant, where the count has expected value pN.
Take the difference between two groups whose only difference is depth. The constant cancels, and what is left is the difference in the mean of the log of a padded count. That is not the log of the padded expected count, because the count is random and the log is curved; under Poisson sampling it is a sum over the possible counts, which is arithmetic with no simulation in it:
depth_lo <- 17000
depth_hi <- 52000
pseudo <- 0.5
ratio_d <- depth_hi / depth_lo
log_gap <- log(ratio_d)
clr_shift <- function(mu_lo, ratio, pc) {
mean_log <- function(m) { # mean of log(count + pc) for a Poisson count
k <- 0:ceiling(m + 12 * sqrt(m) + 30)
sum(dpois(k, m) * log(k + pc))
}
vapply(mu_lo, function(m) mean_log(ratio * m) - mean_log(m), numeric(1)) -
log(ratio)
}
mu_show <- c(0, 0.5, 1, 2, 5, 20)
shift_tab <- data.frame(expected_count = mu_show,
shift = clr_shift(mu_show, ratio_d, pseudo))
half_at <- uniroot(function(m) clr_shift(m, ratio_d, pseudo) + 0.5 * log_gap,
c(1e-8, 100))$root
tenth_at <- uniroot(function(m) clr_shift(m, ratio_d, pseudo) + 0.1 * log_gap,
c(1e-8, 1e4))$root
print(round(shift_tab, 3), row.names = FALSE) expected_count shift
0.0 -1.118
0.5 -0.421
1.0 -0.167
2.0 -0.028
5.0 0.002
20.0 0.000A taxon that is absent from every sample has a padded count of 0.5 in both groups, so under the approximation above its centred log-ratio differs between the groups by minus the log depth ratio: -1.118 at 17000 against 52000 reads. The actual gap is a little smaller, by about one part in the number of taxa, because the taxon’s own padded log is itself one of the terms in the centring mean. That endpoint is closed form and owes nothing to any simulation below.
The interesting part is not the endpoint but the approach to it. At an expected count of one read in the shallow library the shift is already -0.167, at two reads -0.028, and at five reads 0.002. Half of the full bias is present at an expected count of 0.35 reads, and a tenth of it survives out to 1.2 reads. The biological noise in the simulations below has a log standard deviation of 0.4, so the shift at one expected read is a substantial fraction of one sample’s variation, and it is a shift of the mean, which a test on fifteen samples per group can find.
A zero is the limiting case of a small count, and an expected count of a read or two still carries part of the same depth-dependent shrinkage. The curve is steep, though: by an expected count of five reads most of the bias has gone. So the taxa that matter are the ones whose expected count is of order one, and in a real table those are the taxa that turn up as zeros somewhere. The community section below tests that directly, by deleting them.
One arithmetic point disposes of an arm before it is simulated. Adding 0.5 and renormalising to N plus 0.5 times the number of taxa, which is the posterior mean of the proportions under a Jeffreys prior, differs from adding 0.5 alone by a per-sample scale factor, and the centred log-ratio removes per-sample scale factors:
set.seed(4491)
demo_counts <- matrix(rpois(30 * 12, 4), 30, 12)
clr_rows <- function(mat) t(apply(log(mat), 1, function(z) z - mean(z)))
plain_clr <- clr_rows(demo_counts + pseudo)
bayes_clr <- clr_rows((demo_counts + pseudo) /
(rowSums(demo_counts) + pseudo * ncol(demo_counts)))
max_gap <- max(abs(plain_clr - bayes_clr))The two agree to within floating-point rounding. Any replacement that is a constant plus a per-sample rescaling is the fixed pseudocount wearing a different name, and it inherits every property measured below. The Bayesian-multiplicative replacement is not of that kind: it imputes only the zeros and rescales the non-zero counts by a common factor, which keeps their ratios, so this argument does not cover it and it is not measured here.
mu_grid <- exp(seq(log(0.001), log(100), length.out = 300))
shift_df <- rbind(
data.frame(mu = mu_grid, shift = clr_shift(mu_grid, ratio_d, pseudo),
ratio = sprintf("depth ratio %.1f (17k vs 52k)", ratio_d)),
data.frame(mu = mu_grid, shift = clr_shift(mu_grid, 10, pseudo),
ratio = "depth ratio 10 (2k vs 20k)"))
ggplot(shift_df, aes(mu, shift, colour = ratio)) +
geom_hline(yintercept = c(-log_gap, -log(10)), linetype = "dashed",
colour = te_body, linewidth = 0.4) +
geom_vline(xintercept = 1, linetype = "dashed",
colour = te_body, linewidth = 0.5) +
geom_line(linewidth = 1) +
scale_x_log10(breaks = c(0.001, 0.01, 0.1, 1, 10, 100),
labels = c("0.001", "0.01", "0.1", "1", "10", "100")) +
scale_colour_manual(values = c(te_rust, te_forest), name = NULL) +
labs(x = "expected count in the shallow library (reads)",
y = "shift in the centred log-ratio",
title = "The pseudocount bites at counts of order one",
subtitle = "Poisson expectation, no simulation; dashed vertical line: one expected read") +
theme_datasheet() +
theme(legend.position = "bottom")
A null taxon, swept from common to nearly absent
The design follows the metabarcoding check so that the comparison is a comparison and not a change of subject: two groups of fifteen samples, multinomial reads, lognormal biological noise on every taxon, and the same two depth regimes, one confounded at 17000 reads against 52000 and one balanced at thirty thousand against thirty thousand. The difference is the focal taxon. Its true share is identical in the two groups in every cell, so every rejection is a false one, and the share is swept across five values so that its zero share runs from almost none to almost all.
n_taxa <- 40L
n_each <- 15L
sd_bio <- 0.4
n_rep <- 600L
grp_num <- rep(0:1, each = n_each)
grp_fac <- factor(grp_num)
focal <- 1L
share_grid <- c(2e-4, 5e-5, 2e-5, 8e-6, 3e-6)
depth_flat <- 30000
mc_se_rate <- sqrt(0.05 * 0.95 / n_rep)
set.seed(20260920)
base_share <- exp(rnorm(n_taxa, 0, 1.2))
base_share <- base_share / sum(base_share)
cat("taxa", n_taxa, " samples per group", n_each,
" biological log sd", sd_bio, "\n")taxa 40 samples per group 15 biological log sd 0.4 cat("focal shares swept:", format(share_grid, scientific = TRUE), "\n")focal shares swept: 2e-04 5e-05 2e-05 8e-06 3e-06 cat("replicates per cell", n_rep,
" Monte Carlo se at the nominal rate", round(mc_se_rate, 4), "\n")replicates per cell 600 Monte Carlo se at the nominal rate 0.0089 Every constant above was fixed before the first cell ran. The reference community is drawn once, outside the replication loop, so that the five cells differ in the focal taxon’s share and in nothing else; the focal share is set explicitly and the other thirty-nine taxa enter the log-ratio only through the mean of their logs, which is stable from draw to draw. The community section later redraws its base composition in every replicate, because there the composition is the object of study rather than a fixed background.
Seven analyses of that same focal taxon are compared. The first six are applied to the identical table in each replicate. The first is the arm the metabarcoding check ranked at the top under the confound: a centred log-ratio with a fixed 0.5 added, tested with a Welch t test. The second scales the pseudocount with the library, 0.5 times N over the median N, which is the obvious repair and the one a compositional reviewer usually asks for. The third rarefies every library down to the smallest in the run and then takes the same log-ratio; rarefaction here is multinomial resampling with replacement, the usual approximation, rather than subsampling without it. The fourth is a Wilcoxon test on the raw proportions, which needs no pseudocount at all. The fifth is a negative binomial generalised linear model with the log library size as an offset, which puts depth in the model rather than in the data. The sixth is the log of a padded proportion, which is the arm the metabarcoding check actually recommended, written with the same 0.5 and the same denominator convention. The seventh is the two-part test a microbiome statistician asks for when a log-ratio looks bimodal; it runs on its own draws from the same generator, and it shares the next section with a rank test on the same log-ratio values.
gen_table <- function(share, d_lo, d_hi) {
lib <- c(rep(d_lo, n_each), rep(d_hi, n_each))
p_vec <- base_share
p_vec[focal] <- share
p_vec <- p_vec / sum(p_vec)
cnt <- matrix(0L, 2 * n_each, n_taxa)
for (i in seq_len(2 * n_each)) {
q_i <- p_vec * exp(rnorm(n_taxa, 0, sd_bio))
cnt[i, ] <- as.integer(rmultinom(1, lib[i], q_i / sum(q_i)))
}
list(cnt = cnt, lib = lib)
}
clr_focal <- function(cnt, pc) {
lg <- log(cnt + pc)
lg[, focal] - rowMeans(lg)
}
welch_p <- function(y) t.test(y[grp_num == 0], y[grp_num == 1])$p.value
rank_p <- function(y) { # normal approximation whenever values tie, in every R version
a <- y[grp_num == 0]; b <- y[grp_num == 1]
wilcox.test(a, b, exact = if (anyDuplicated(c(a, b))) FALSE else NULL)$p.value
}
rarefy_rows <- function(cnt, depth) {
t(apply(cnt, 1, function(r_i) rmultinom(1, depth, r_i / sum(r_i))[, 1]))
}
one_rep <- function(share, d_lo, d_hi) {
tab <- gen_table(share, d_lo, d_hi)
cnt <- tab$cnt
lib <- tab$lib
prop <- cnt[, focal] / lib
cl_fix <- clr_focal(cnt, pseudo)
c(zero_share = mean(cnt[, focal] == 0),
clr_fixed = welch_p(cl_fix),
clr_scaled = welch_p(clr_focal(cnt, pseudo * lib / median(lib))),
clr_rarefy = welch_p(clr_focal(rarefy_rows(cnt, min(lib)), pseudo)),
wilcox_prop = suppressWarnings(rank_p(prop)),
nb_offset = tryCatch(suppressWarnings(
summary(glm.nb(cnt[, focal] ~ grp_num + offset(log(lib))))
)$coefficients["grp_num", 4], error = function(e) NA_real_),
clr_rank = suppressWarnings(rank_p(cl_fix)),
prop_pad = tryCatch(welch_p(log((cnt[, focal] + pseudo) /
(lib + pseudo * n_taxa))),
error = function(e) NA_real_))
}
run_cell <- function(share, d_lo, d_hi) {
mat <- t(replicate(n_rep, one_rep(share, d_lo, d_hi)))
out <- c(zero_share = mean(mat[, "zero_share"]),
colMeans(mat[, -1] < 0.05, na.rm = TRUE))
c(out, usable_wilcox = sum(!is.na(mat[, "wilcox_prop"])),
usable_nb = sum(!is.na(mat[, "nb_offset"])),
usable_pad = sum(!is.na(mat[, "prop_pad"])))
}
set.seed(7731)
res_conf <- t(vapply(share_grid, function(s) run_cell(s, depth_lo, depth_hi),
numeric(11)))
res_flat <- t(vapply(share_grid, function(s) run_cell(s, depth_flat, depth_flat),
numeric(11)))
rate_conf <- as.data.frame(res_conf)
rate_flat <- as.data.frame(res_flat)
rate_conf$share <- share_grid
rate_flat$share <- share_grid
z_nominal <- function(r) (r - 0.05) / mc_se_rate
show_cols <- c("clr_fixed", "clr_scaled", "clr_rarefy", "wilcox_prop",
"nb_offset")
hold_max <- max(c(rate_conf$clr_rarefy, rate_conf$nb_offset))
nb_min <- min(rate_conf$nb_offset)
rar_min <- min(rate_conf$clr_rarefy)
scaled_max <- max(rate_conf$clr_scaled)
flat_max <- max(rate_flat[, show_cols])
pad_gap <- max(abs(rate_conf$prop_pad - rate_conf$clr_fixed))
wilcox_unusable <- n_rep - rate_conf$usable_wilcox[5]
nb_unusable <- n_rep - rate_conf$usable_nb[5]
pad_unusable <- n_rep - rate_conf$usable_pad[5]
print(round(rate_conf[, c("zero_share", show_cols, "prop_pad")], 3),
row.names = FALSE) zero_share clr_fixed clr_scaled clr_rarefy wilcox_prop nb_offset prop_pad
0.025 0.062 0.047 0.055 0.047 0.063 0.060
0.259 0.102 0.140 0.048 0.070 0.057 0.092
0.521 0.588 0.135 0.040 0.142 0.038 0.597
0.755 0.980 0.127 0.028 0.193 0.032 0.987
0.894 1.000 0.050 0.013 0.075 0.007 1.000print(round(rate_flat[, c("zero_share", show_cols)], 3), row.names = FALSE) zero_share clr_fixed clr_scaled clr_rarefy wilcox_prop nb_offset
0.009 0.047 0.047 0.050 0.045 0.075
0.233 0.038 0.038 0.040 0.037 0.060
0.538 0.055 0.055 0.043 0.052 0.047
0.774 0.040 0.040 0.027 0.047 0.012
0.905 0.028 0.028 0.015 0.019 0.000cat("negative binomial fits that returned an estimate, per cell, of", n_rep, ":",
rate_conf$usable_nb, "\n")negative binomial fits that returned an estimate, per cell, of 600 : 600 600 600 600 583 The five cells of the confounded regime put the focal taxon’s zero share at 0.03, 0.26, 0.52, 0.76, 0.89. The fixed pseudocount log-ratio rejects a true null in 0.062 of tables at the first of those, which reproduces the metabarcoding check’s verdict, and then in 0.102, 0.588, 0.980 and 1.000 as the taxon gets sparser, against a Monte Carlo standard error of 0.0089 at the nominal rate. At the sparsest cell the test rejects every table it is given.
At balanced depth the same analysis on the same generator rejects in 0.047, 0.038, 0.055, 0.040, 0.028 across the same five focal shares. The design constant across the two panels is the share, not the zero share: at equal depth the focal taxon is absent in 0.01, 0.23, 0.54, 0.77, 0.90 of samples, which differs from the confounded panel most, in relative terms, at the first cell. The pseudocount is not the problem on its own; the pseudocount with a depth difference is.
The repairs separate cleanly, and they miss in the other direction. Rarefying to the smallest library first and then taking the identical log-ratio gives 0.055, 0.048, 0.040, 0.028, 0.013 in the confounded regime, and the negative binomial model with a log depth offset gives 0.063, 0.057, 0.038, 0.032, 0.007. Neither picks up the depth contrast: across those ten confounded-regime rates the highest either of them reaches is 0.063, 1.5 Monte Carlo standard errors above the nominal rate, and the count model’s best-observed cell, 0.063, is matched at equal depth by 0.075, so its slight excess there is not a depth effect. Both go conservative as the taxon empties, and the count model goes furthest: its lowest cell is 0.007, 4.8 standard errors below nominal, against the rarefied 0.013 at 4.1 below. At the sparsest cell the negative binomial fit also failed in 17 of 600 replicates, the tables in which the taxon was absent from every sample, and that cell’s rate is computed on the rest.
Scaling the pseudocount with the library is a different matter: it peaks at 0.140 at a zero share of 0.26, 10.1 standard errors above nominal, and it is above nominal in the middle cells of the sweep. The reason is visible in the design: scaling the pseudocount equalises what is added to a zero relative to depth, but the two groups do not have the same number of zeros to begin with, so the correction is applied to different fractions of each group.
The Wilcoxon test on proportions, which never touches a pseudocount, reaches 0.193 at its worst cell. Its failure is a different one, and not the ties as such: at equal depth, where the zeros tie about as often, the same test reaches only 0.052. A rank test asks whether one group tends to sit above the other, and at unequal depth the shallow libraries miss the taxon more often, so more of their values are zeros at the bottom of the ranking, while their non-zero proportions are coarser and tend to be larger; the ranks carry depth although the mean proportion is the same in both groups. At the sparsest cell the test is undefined in 17 of 600 replicates, where the taxon was absent from every sample in the table. Where values tie, the test here uses the normal approximation with a tie-corrected variance, which is what a default wilcox.test() did before R 4.6.0; the call asks for it explicitly so that the rates do not depend on the version of R, and What a rank test actually tests sets that approximation against the exact conditional test that R 4.6.0 computes instead.
The arm the metabarcoding check recommended does not escape either, once a zero forces a pseudocount back into it. The log of a padded proportion, the log of x plus 0.5 over N plus 0.5 times the number of taxa, tested with the same Welch test, rejects in 0.060, 0.092, 0.597, 0.987, 1.000. The largest difference between that curve and the centred log-ratio curve anywhere in the sweep is 0.010, or 1.1 Monte Carlo standard errors, so the two are the same curve for practical purposes and the padded proportion is left out of the figure below to avoid drawing one line on top of another. Like the Wilcoxon test, it is undefined in 17 of 600 replicates at the sparsest cell: where the taxon is absent from the whole table the padded proportion takes one constant value in each group, and the two constants differ by very nearly the log depth ratio, which is the closed form of the first section arriving without a test to report it. Centring is not what breaks. The constant added under the log is.
meth_key <- c(clr_fixed = "CLR, fixed 0.5",
clr_scaled = "CLR, depth-scaled 0.5",
clr_rarefy = "rarefy, then CLR",
wilcox_prop = "Wilcoxon on proportions",
nb_offset = "NB GLM, log-depth offset")
long_rates <- function(tab, lab) {
do.call(rbind, lapply(names(meth_key), function(m) {
data.frame(zero_share = tab$zero_share, rate = tab[[m]],
method = meth_key[[m]], regime = lab)
}))
}
rates_df <- rbind(
long_rates(rate_conf, sprintf("%d vs %d reads", depth_lo, depth_hi)),
long_rates(rate_flat, sprintf("%d vs %d reads", depth_flat, depth_flat)))
rates_df$method <- factor(rates_df$method, levels = unname(meth_key))
ggplot(rates_df, aes(zero_share, rate, colour = method)) +
geom_hline(yintercept = 0.05, linetype = "dashed",
colour = te_body, linewidth = 0.5) +
geom_line(linewidth = 0.9) +
geom_point(size = 1.8) +
facet_wrap(~regime) +
scale_colour_manual(values = c(te_rust, te_gold, te_forest,
"#8d8b78", te_ink), name = NULL) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "zero share of the focal taxon", y = "false positive rate at 0.05",
title = "Only rarefying and the count model resist the depth contrast",
subtitle = "the focal taxon's true share is identical in both groups in every cell") +
theme_datasheet() +
theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 2))
The test is not the problem, and a zero-aware test is no cure
A reasonable objection is that the failure is the Welch test rather than the pseudocount. A centred log-ratio on a sparse taxon is bimodal: the zero samples pile up at the padded value and the rest sit above it, and a t test on such a thing has no business being trusted. Two checks settle this. The first applies a Wilcoxon rank sum test to the identical log-ratio values. The second is the two-part test that is usually proposed for zero-inflated compositions: a Fisher exact test on presence and absence by group, a Welch test on the log-ratio among the samples where the taxon was seen, and Fisher’s method to combine the two.
two_part_p <- function(cnt, lib) {
x_focal <- cnt[, focal]
present <- x_focal > 0
p_pres <- tryCatch(
fisher.test(table(factor(present, c(FALSE, TRUE)), grp_num))$p.value,
error = function(e) NA_real_)
cl_val <- clr_focal(cnt, pseudo)
keep <- present & is.finite(cl_val)
p_pos <- if (sum(keep & grp_num == 0) >= 2 && sum(keep & grp_num == 1) >= 2) {
tryCatch(t.test(cl_val[keep & grp_num == 0],
cl_val[keep & grp_num == 1])$p.value,
error = function(e) NA_real_)
} else NA_real_
# shallow minus deep mean log-ratio among the detected samples
pos_gap <- if (is.na(p_pos)) NA_real_ else
mean(cl_val[keep & grp_num == 0]) - mean(cl_val[keep & grp_num == 1])
p_both <- c(p_pres, p_pos)
p_both <- p_both[!is.na(p_both)]
p_comb <- if (!length(p_both)) NA_real_ else
pchisq(-2 * sum(log(p_both)), df = 2 * length(p_both), lower.tail = FALSE)
c(comb = p_comb, pres = p_pres, pos = p_pos, pos_gap = pos_gap)
}
set.seed(5518)
part_p <- lapply(share_grid, function(s) {
t(replicate(n_rep, {
tab <- gen_table(s, depth_lo, depth_hi)
two_part_p(tab$cnt, tab$lib)
}))
})
two_part_rate <- vapply(part_p, function(m) mean(m[, "comb"] < 0.05, na.rm = TRUE),
numeric(1))
pres_rate <- vapply(part_p, function(m) mean(m[, "pres"] < 0.05, na.rm = TRUE),
numeric(1))
# a table where the Welch part cannot run counts as no rejection
pos_rate <- vapply(part_p, function(m) mean(!is.na(m[, "pos"]) & m[, "pos"] < 0.05),
numeric(1))
pos_missing <- vapply(part_p, function(m) mean(is.na(m[, "pos"])), numeric(1))
pos_gap <- vapply(part_p, function(m) mean(m[, "pos_gap"], na.rm = TRUE), numeric(1))
gap_cell <- which.max(two_part_rate - rate_conf$clr_fixed)
form_tab <- data.frame(zero_share = round(rate_conf$zero_share, 2),
clr_welch = rate_conf$clr_fixed,
clr_wilcoxon = rate_conf$clr_rank,
two_part = two_part_rate,
presence_only = pres_rate,
detected_welch_only = pos_rate,
no_welch_part = pos_missing)
print(round(form_tab, 3), row.names = FALSE) zero_share clr_welch clr_wilcoxon two_part presence_only detected_welch_only
0.03 0.062 0.057 0.023 0.000 0.065
0.26 0.102 0.088 0.915 0.422 0.895
0.52 0.588 0.375 0.977 0.300 0.952
0.76 0.980 0.912 0.647 0.143 0.598
0.89 1.000 0.998 0.092 0.025 0.108
no_welch_part
0.000
0.000
0.027
0.333
0.862The rank test on the same log-ratio values rejects in 0.057, 0.088, 0.375, 0.912, 0.998 across the sweep. It is better than the Welch test in the middle of the range, though still far above nominal, and just as bad at the end. The failure is in the values, not in the test applied to them.
In the middle of the sweep the two-part test is worse than either. On independently drawn tables from the same generator it rejects in 0.023, 0.915, 0.977, 0.647, 0.092, reaching 0.977 at a zero share of 0.52. The widest gap is at a zero share of 0.26, where the two-part test rejects 0.915 of tables and the plain Welch test on the same log-ratio rejects 0.102 on its own tables. Making the test zero-aware makes it depth-aware in the wrong direction, and mostly not through the presence part. A deeper library does see a rare taxon more often, but the Fisher test on presence alone rejects in 0.000, 0.422, 0.300, 0.143, 0.025 across the sweep, while the Welch test on the detected samples alone rejects in 0.065, 0.895, 0.952, 0.598, 0.108, counting a table where it cannot run as no rejection. A sample counts as detected only when it has at least one read, which is a larger share of a shallow library, so among the detected samples the shallow group’s log-ratios sit higher: by 0.60 of a log unit on average at a zero share of 0.26. Fisher’s method then adds the evidence of the two parts together, and although a Fisher exact p on a two by two presence table is discrete and conservative, in the three middle cells the combined test rejects more often than either part. The rate falls again over the two sparsest cells, below both of the other tests: there the taxon is seen in few samples, so both parts have less to work with, and the Welch part cannot run at all in 33 and 86 per cent of those tables, which leaves the presence table on its own.
The same bias moves a whole community
The single-taxon result is the easy one to argue away: analysts do not test one taxon, they run a PERMANOVA. Weiss and colleagues warned in 2017 that normalisations other than rarefying are vulnerable to artefacts from library size, and reported that for groups differing about tenfold in average library size rarefying lowers the false discovery rate of a differential abundance analysis. They did not test a log-ratio distance under a null grouping, and their normalisation set was rarefying against DESeq-type and cumulative-sum scaling rather than the centred log-ratio. This section is an extension of that warning to the arrangement an ecologist is most likely to run, in the simplest form that shows it. There is no group effect at all in the generator: the two groups are drawn from the same community, and the only thing that differs is the library size.
The base composition is redrawn in every replicate, so the rejection rate is an average over communities rather than a property of one lucky draw. Five analyses run on each table: the Aitchison distance, which is Euclidean distance on the centred log-ratio with a fixed 0.5, Bray-Curtis on proportions, both of those again after rarefying to the smallest library, and the Aitchison distance once more on only those taxa that have no zero anywhere in the table. The last one is the arm that says where the effect lives, and the code also keeps the number of taxa it had to delete.
n_taxa2 <- 60L
n_rep2 <- 120L
n_perm <- 199L
depth_pairs <- list(c(20000, 20000), c(5000, 20000), c(2000, 20000),
c(1000, 20000), c(500, 20000))
gen_table2 <- function(d_lo, d_hi) {
b_vec <- exp(rnorm(n_taxa2, 0, 1.5))
b_vec <- b_vec / sum(b_vec)
lib <- c(rep(d_lo, n_each), rep(d_hi, n_each))
cnt <- matrix(0L, 2 * n_each, n_taxa2)
for (i in seq_len(2 * n_each)) {
q_i <- b_vec * exp(rnorm(n_taxa2, 0, sd_bio))
cnt[i, ] <- as.integer(rmultinom(1, lib[i], q_i / sum(q_i)))
}
list(cnt = cnt, lib = lib)
}
aitchison <- function(cnt) dist(clr_rows(cnt + pseudo))
one_rep2 <- function(d_lo, d_hi) {
tab <- gen_table2(d_lo, d_hi)
cnt <- tab$cnt
cnt_r <- rarefy_rows(cnt, min(tab$lib))
d_ait <- aitchison(cnt)
complete <- apply(cnt, 2, min) > 0
c(aitchison = adonis2(d_ait ~ grp_fac, permutations = n_perm)$`Pr(>F)`[1],
bray = adonis2(vegdist(cnt / rowSums(cnt), "bray") ~ grp_fac,
permutations = n_perm)$`Pr(>F)`[1],
aitchison_rar = adonis2(aitchison(cnt_r) ~ grp_fac,
permutations = n_perm)$`Pr(>F)`[1],
bray_rar = adonis2(vegdist(cnt_r / rowSums(cnt_r), "bray") ~ grp_fac,
permutations = n_perm)$`Pr(>F)`[1],
dispersion = permutest(betadisper(d_ait, grp_fac),
permutations = n_perm)$tab$`Pr(>F)`[1],
aitchison_nz = if (sum(complete) >= 3) {
adonis2(aitchison(cnt[, complete, drop = FALSE]) ~ grp_fac,
permutations = n_perm)$`Pr(>F)`[1]
} else NA_real_,
taxa_with_zero = sum(!complete),
min_kept = if (any(complete)) {
median(apply(cnt[, complete, drop = FALSE], 2, min))
} else NA_real_,
zero_share = mean(cnt == 0))
}
set.seed(9081)
act_two <- as.data.frame(t(vapply(depth_pairs, function(dp) {
mat <- t(replicate(n_rep2, one_rep2(dp[1], dp[2])))
c(depth_lo = dp[1], ratio = dp[2] / dp[1],
zero_share = mean(mat[, "zero_share"]),
taxa_with_zero = median(mat[, "taxa_with_zero"]),
min_kept = median(mat[, "min_kept"], na.rm = TRUE),
colMeans(mat[, 1:5] < 0.05),
aitchison_nz = mean(mat[, "aitchison_nz"] < 0.05, na.rm = TRUE))
}, numeric(11))))
mc_se2 <- sqrt(0.25 / n_rep2)
mc_se2_nom <- sqrt(0.05 * 0.95 / n_rep2)
print(round(act_two, 3), row.names = FALSE) depth_lo ratio zero_share taxa_with_zero min_kept aitchison bray
20000 1 0.002 1 48.5 0.075 0.058
5000 4 0.011 6 14.0 0.167 0.033
2000 10 0.036 15 7.0 0.633 0.067
1000 20 0.062 24 4.5 0.900 0.042
500 40 0.112 35 3.0 1.000 0.150
aitchison_rar bray_rar dispersion aitchison_nz
0.058 0.067 0.075 0.058
0.083 0.058 0.867 0.083
0.042 0.058 1.000 0.125
0.025 0.033 1.000 0.100
0.033 0.033 1.000 0.117With both groups at twenty thousand reads the Aitchison PERMANOVA rejects the null grouping in 0.075 of tables, which is consistent with the nominal rate. Shorten one group’s libraries and it climbs: 0.167 at a fourfold depth difference, 0.633 at tenfold, 0.900 at twentyfold, and 1.000 at fortyfold, with a Monte Carlo standard error of at most 0.046. Bray-Curtis on proportions is far better behaved but not clean: it runs from 0.033 to 0.150, and its worst cell is 5.0 standard errors above nominal on the standard error of 0.020 that applies at the nominal rate. Rarefying to the smallest library first holds both distances down: the highest rate reached by the rarefied Aitchison distance anywhere on the grid is 0.083 and by rarefied Bray-Curtis 0.067.
Which taxa carry it is a separate question from how many cells are zero, and the two counts give opposite impressions. At the tenfold contrast, where the Aitchison test rejects 63 per cent of null groupings, only 3.6 per cent of the cells in the table are zero, which sounds like far too few to matter. The zeros are not spread evenly over the table. They sit in a median of 15 of the 60 taxa, 25 per cent of them. Delete exactly those taxa, keep the ones observed in every sample, and the Aitchison test on what is left rejects in 0.125 of tables instead of 0.633. The same deletion takes the fortyfold contrast from 1.000 to 0.117.
So most of it is carried by the taxa that touch the pseudocount. Not all of it: at the fortyfold contrast the zero-free analysis is still 3.4 standard errors above nominal, so the taxa observed in every sample carry a smaller part. The survivors are not an abundant subset: at the tenfold contrast the median retained taxon’s smallest count anywhere in the table is 7 reads. The Poisson calculation of the first section puts the mean depth shift of a taxon with that expected count at 0.001 of a log unit, against -2.303 for a cell that is zero. The share of cells is the wrong denominator for the question, because Euclidean distance adds up whole taxa, and a taxon that is zero in some of the shallow libraries is displaced along its own axis in every sample it has.
The dispersion test tells the same story in the vocabulary of four common PERMANOVA mistakes. Betadisper on the Aitchison distance rejects in 0.075 of runs at equal depth and 1.000 at tenfold, so the shallow group is genuinely more spread out on this distance as well as displaced. Anderson and Walsh set out in 2013 why a PERMANOVA on groups with unequal dispersion is not always testing what its user thinks, and found it largely unaffected by unequal dispersion in balanced designs like this one; here the unequal dispersion is manufactured by the pseudocount and the depth together, and a reader who checks betadisper and finds it significant will at least know that the groups differ in spread as well as in location.
perm_key <- c(aitchison = "Aitchison (CLR, fixed 0.5)",
aitchison_nz = "Aitchison, zero-free taxa only",
bray = "Bray-Curtis on proportions",
aitchison_rar = "rarefy, then Aitchison",
dispersion = "betadisper on Aitchison")
perm_df <- do.call(rbind, lapply(names(perm_key), function(m) {
data.frame(ratio = act_two$ratio, rate = act_two[[m]], test = perm_key[[m]])
}))
perm_df$test <- factor(perm_df$test, levels = unname(perm_key))
ggplot(perm_df, aes(ratio, rate, colour = test)) +
geom_hline(yintercept = 0.05, linetype = "dashed",
colour = te_body, linewidth = 0.5) +
geom_line(linewidth = 0.9) +
geom_point(size = 2) +
scale_x_log10(breaks = act_two$ratio,
labels = sprintf("%.0f", act_two$ratio)) +
scale_y_continuous(limits = c(0, 1)) +
scale_colour_manual(values = c(te_rust, te_ink, te_forest, "#8d8b78",
te_gold), name = NULL) +
labs(x = "ratio of the two groups' library sizes",
y = "rejection rate under a null grouping",
title = "A depth contrast read as a community difference",
subtitle = "no group effect exists in the generator at any point on these curves") +
theme_datasheet() +
theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 3))
Whether an ordination shows it depends on the sparsity
A natural defence is that this would be visible in an ordination: run a principal component analysis on the centred log-ratio matrix and the depth gradient will be the first axis, so anyone who looks at the picture will see the problem. That is a claim about a population of tables and not about one, so one table cannot settle it. The arrangement here is forty samples with no groups at all, library sizes spread smoothly from one thousand to forty thousand reads, and three settings that change how sparse the table is: sixty taxa with a moderate spread of abundances, two hundred taxa with the same spread, and two hundred with a wider one. Each draw gets a fresh community. For each draw the code regresses each of the first six principal components on log library size and keeps the r squared.
n_grad <- 200L
n_grad_ad <- 120L
n_samp3 <- 40L
lib3 <- round(exp(seq(log(1000), log(40000), length.out = n_samp3)))
log_lib <- log(lib3)
grad_set <- list(c(60, 1.5), c(200, 1.5), c(200, 2.5))
grad_lab <- c("60 taxa, spread 1.5", "200 taxa, spread 1.5",
"200 taxa, spread 2.5")
gen_grad <- function(n_tx, sd_base) {
b_vec <- exp(rnorm(n_tx, 0, sd_base))
b_vec <- b_vec / sum(b_vec)
cnt <- matrix(0L, n_samp3, n_tx)
for (i in seq_len(n_samp3)) {
q_i <- b_vec * exp(rnorm(n_tx, 0, sd_bio))
cnt[i, ] <- as.integer(rmultinom(1, lib3[i], q_i / sum(q_i)))
}
cnt
}
depth_r2 <- function(cnt) {
sc <- prcomp(clr_rows(cnt + pseudo))$x[, 1:6]
apply(sc, 2, function(z) summary(lm(z ~ log_lib))$r.squared)
}
set.seed(771)
grad_raw <- lapply(grad_set, function(g) t(replicate(n_grad, {
cnt <- gen_grad(as.integer(g[1]), g[2])
c(depth_r2(cnt), zero = mean(cnt == 0))
})))
names(grad_raw) <- grad_lab
grad_sum <- data.frame(
zeros = vapply(grad_raw, function(m) median(m[, "zero"]), numeric(1)),
pc1_med = vapply(grad_raw, function(m) median(m[, 1]), numeric(1)),
pc1_q25 = vapply(grad_raw, function(m) unname(quantile(m[, 1], 0.25)),
numeric(1)),
pc1_q75 = vapply(grad_raw, function(m) unname(quantile(m[, 1], 0.75)),
numeric(1)),
pc1_best = vapply(grad_raw,
function(m) mean(apply(m[, 1:6], 1, which.max) == 1),
numeric(1)))
print(round(grad_sum, 3)) zeros pc1_med pc1_q25 pc1_q75 pc1_best
60 taxa, spread 1.5 0.030 0.159 0.047 0.339 0.6
200 taxa, spread 1.5 0.101 0.833 0.786 0.865 1.0
200 taxa, spread 2.5 0.361 0.978 0.971 0.983 1.0set.seed(3324)
grad_ad <- t(replicate(n_grad_ad, {
cnt <- gen_grad(60L, 1.5)
cnt_r <- rarefy_rows(cnt, min(lib3))
c(aitchison = adonis2(aitchison(cnt) ~ log_lib,
permutations = n_perm)$`Pr(>F)`[1],
bray = adonis2(vegdist(cnt / rowSums(cnt), "bray") ~ log_lib,
permutations = n_perm)$`Pr(>F)`[1],
rarefied = adonis2(aitchison(cnt_r) ~ log_lib,
permutations = n_perm)$`Pr(>F)`[1])
}))
grad_rate <- colMeans(grad_ad < 0.05)
mc_se_ad <- sqrt(0.05 * 0.95 / n_grad_ad)
cat("gradient PERMANOVA rejection at 0.05:", round(grad_rate, 3), "\n")gradient PERMANOVA rejection at 0.05: 0.683 0.075 0.1 The answer depends on how sparse the table is, and not in the direction that would be comfortable. On the sixty-taxon setting, the one the community section used and the least sparse of the three at 3.0 per cent zero cells, log library size explains a median of 0.159 of the variance in PC1, with the middle half of the draws between 0.047 and 0.339, and PC1 is the component most strongly related to depth in 120 of the 200 draws. Raise the richness to two hundred taxa and the zero share goes to 10.1 per cent, log depth explains a median 0.833 of PC1, and PC1 is the depth axis in 200 of 200 draws. Widen the abundance spread as well, to 36.1 per cent zeros, and the median rises to 0.978.
A real 16S table has hundreds to thousands of amplicon variants and is mostly zeros, so on this measure it sits at the sparse end of that range, and there, when library sizes vary as widely as they do here and no real community structure competes with them, the first axis of a log-ratio ordination can be largely a picture of library size. An analyst who reads that axis as community structure without checking it against depth may be reading the sequencer. The comfortable version of the claim, that the depth gradient would announce itself on PC1 and so could not be missed, is the one that fails: it fails exactly where the table is nearly complete, and there the artefact does not go away, it stops dominating the first axis.
grad_long <- do.call(rbind, lapply(seq_along(grad_raw), function(j) {
m <- grad_raw[[j]]
data.frame(
regime = sprintf("%s\n%.0f per cent zeros", grad_lab[j],
100 * median(m[, "zero"])),
pc = factor(rep(1:6, each = nrow(m))),
r2 = as.vector(m[, 1:6]))
}))
grad_long$regime <- factor(grad_long$regime, levels = unique(grad_long$regime))
grad_long$first <- ifelse(grad_long$pc == "1", "PC1", "later")
ggplot(grad_long, aes(pc, r2, fill = first)) +
geom_boxplot(colour = te_ink, linewidth = 0.35, outlier.size = 0.5,
outlier.colour = te_body) +
facet_wrap(~regime, nrow = 1) +
scale_fill_manual(values = c(PC1 = te_rust, later = te_line),
guide = "none") +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "principal component of the CLR matrix",
y = "r squared against log library size",
title = "The sparser the table, the more depth lands on the first axis",
subtitle = sprintf("%d draws per panel; forty samples, no groups, depth from 1k to 40k reads", n_grad)) +
theme_datasheet()
The test that names depth as a predictor finds it in most draws even where the picture is weakest. On the sixty-taxon setting, a PERMANOVA with log library size as a continuous predictor rejects in 0.683 of draws on the Aitchison distance, against 0.075 on Bray-Curtis on proportions and 0.100 after rarefying, over 120 draws with a standard error of 0.020 at the nominal rate. That is the regime where the picture would not have warned anyone: in half the draws log depth accounts for less than 0.159 of PC1, and the p value finds the shift anyway. So an ordination is evidence when it does show a depth gradient and is not evidence when it does not.
What to report
Give the zero share of every taxon you make a claim about, next to the claim. The false positive rate of the fixed-pseudocount log-ratio was 0.062 at a zero share of 0.03 and 0.980 at a zero share of 0.76, on the same generator, the same depths and the same test. One number separates a trustworthy result from a certainty of the wrong kind, and it costs nothing to print.
Report the library sizes by group, and report them as a distribution, not a mean. A threefold difference in median depth between two treatments is enough for the single-taxon results above, and a fourfold one was already enough for the Aitchison PERMANOVA to reject a null grouping in 0.167 of tables. If the depths are balanced, say so: at equal depth the highest rate any of the five analyses in the figure reached, at any sparsity, was 0.075.
Say which pseudocount was used and whether it is a constant. If it is, the claim about any sparse taxon is partly a claim about that constant, which is the point the zero problem makes in general. Adding 0.5 and renormalising is the same thing under a log-ratio, up to floating-point rounding, so calling it a posterior mean does not buy anything here; the Bayesian-multiplicative replacement, which imputes only the zeros, is a different method and was not tested.
For a differential abundance question on a sparse taxon under unequal depth, report the count model, and report where it is conservative. The negative binomial fit with an offset of log library size never picked up the depth contrast: 0.063, 0.057, 0.038, 0.032, 0.007 across the sweep. Its highest rate is 1.5 Monte Carlo standard errors above nominal, and on the best-observed taxon it is 0.063 here and 0.075 at equal depth, so its test is slightly liberal there whatever the depth; its lowest is 4.8 below, so on the sparsest taxa it is conservative rather than valid, and it loses power there accordingly. It does not discard reads, and depth belongs in the model. If a distance-based community test is the question instead, rarefy: on the Aitchison distance at a tenfold depth contrast the null grouping was rejected in 0.633 of tables without rarefying and 0.042 with it.
Run betadisper whenever the distance is an Aitchison distance and the depths are unequal. It rejected in 1.000 of the tenfold-contrast tables here, and a significant dispersion test is the cheapest available warning that a location claim is not safe.
Honest limits
The repair is not symmetric with the problem. Rarefying protects the error rate and costs power, exactly as checking a metabarcoding analysis measured on a taxon that had no zeros, and nothing above measures that cost at these sparsities. Power is one part of the case McMurdie and Holmes made against rarefying in 2014, on the ground that discarding reads is statistically inefficient; they also reported a high false positive rate for rarefied counts in their own simulations of differential abundance, a design different from the one here. This post does not dispute the power argument; it bounds where it applies. Rarefaction is implemented here as multinomial resampling with replacement rather than subsampling without it, which is the usual approximation and adds a little variance to the shallow group’s own counts. The honest reading of the two posts together is that the ranking depends on the zero share: where a taxon is well observed the log-ratio is both valid and more powerful, and where it is sparse the log-ratio is invalid and the comparison of power does not arise. No single rule covers a whole table, and a table has taxa of both kinds in it.
The library sizes within a group are fixed here, not drawn. That keeps the closed-form endpoint simple and the simulation clean, and it is not how a run comes off a machine: the metabarcoding check used a lognormal library size with a log standard deviation of 0.55 around each group’s median, which is more realistic. How much within-group spread in depth dilutes the contrast at a given median ratio is not measured here, and should not be assumed large.
Whether a threefold median depth difference between batches is typical is a question about laboratories, not about statistics, and this post cannot answer it. It is the value the earlier post used and it was kept so that the two are comparable. The community-level PERMANOVA curve is the better guide to the general shape: the rejection rate under a null grouping was 0.167 at a fourfold ratio and 1.000 at fortyfold, so the artefact grows smoothly with the depth ratio rather than appearing suddenly at one value.
Every table here is multinomial with lognormal biological noise and no overdispersion beyond that, no compositional group effect, and no correlation structure among taxa. Real amplicon data have all three. The negative binomial model was fitted to data that a negative binomial does not exactly describe, and depth still did not mislead it, though its test ran slightly above nominal on the best-observed taxon at either depth; that is reassuring rather than conclusive, and a model-based analysis should be checked against its own residuals rather than trusted because it worked in a simulation with a different generator.
The two-part test measured here is one specific construction, Fisher exact on presence plus a Welch test on the observed log-ratios, combined by Fisher’s method. A hurdle model that puts log depth into the count part for the detected samples and into the presence part as well should not fail the way this one does, because the depth dependence of both parts would then be modelled rather than tested, though that is not measured here. The result above should be read as a warning about the naive two-part test, not as a claim that no zero-aware method can work.
The community generator has no group effect at all, so everything reported there is a false positive rate and nothing is a power comparison. A rarefied analysis that holds its error rate is not thereby the better analysis; it is the safer one. Deciding between them needs the effect size that matters to the study, which is a design question and not a simulation result.
The ordination section varies richness and the spread of abundances and nothing else. Its depth gradient is smooth and even and spans a factor of forty, the biological noise has the same log standard deviation everywhere, and there is no community structure at all for the first axis to compete with. A real survey has none of that tidiness, and its first axis usually has real ecology on it. What the three settings establish is the direction of the dependence, not a threshold zero share at which depth takes over PC1.
Finally, 600 replicates per cell in the single-taxon sections and 120 in the community section fix the resolution of every rate quoted: a single-taxon rate has a standard error of 0.0089 at the nominal rate and of at most 0.020 anywhere, and a community rate one of at most 0.046. Differences smaller than a few of those standard errors are not differences, and the replication was chosen before the runs and not revised.
References
Weiss S, Xu ZZ, Peddada S, Amir A, Bittinger K, Gonzalez A, Lozupone C, Zaneveld JR, Vazquez-Baeza Y, Birmingham A, Hyde ER, Knight R 2017 Microbiome 5(1):27 (10.1186/s40168-017-0237-y)
McMurdie PJ, Holmes S 2014 PLoS Computational Biology 10(4):e1003531 (10.1371/journal.pcbi.1003531)
Anderson MJ, Walsh DCI 2013 Ecological Monographs 83(4):557-574 (10.1890/12-2010.1)