Power for a community study

vegan

PERMANOVA

community ecology

study design

ecology tutorial

How many plots does a PERMANOVA need? Simulation shows the answer depends on the dissimilarity index, which decides what counts as a community difference.

Author

Tidy Ecology

Published

2026-08-05

Sample size questions have a standard answer for a t-test and no answer at all for a community study. The literature is full of designs with five or eight plots per group, an adonis2 result, and no statement of what the study could have detected. The gap is not that power analysis is hard here; it is that the first input is missing. Before you can ask how many plots, you have to say what a difference between communities means, and that decision is made by the dissimilarity index, not by the test.

This post simulates a PERMANOVA design end to end: a known difference, a sample size grid, three common indices, and a null case as a calibration check. The result is that the power of the same design, on the same data, over the same range of sample sizes, spans the full interval from nothing to certainty depending on which index you picked in the first line of the script.

A community and two ways of changing it

Thirty species with log-normal mean abundances and negative binomial counts, so the mean-variance relationship is realistic rather than Poisson. Two treatment scenarios, both plausible field outcomes.

library(vegan)
library(ggplot2)

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"),
          strip.text       = element_text(colour = te_ink, face = "bold"),
          axis.text        = element_text(colour = te_body))
}

n_sp <- 30
set.seed(20260805)
mu_base <- exp(rnorm(n_sp, mean = 1.6, sd = 0.9))

## scenario A: every taxon multiplied by 1.5, so proportions are untouched
mu_A <- mu_base * 1.5

## scenario B: turnover, with the community total held at the reference total
odd  <- seq_len(n_sp) %% 2 == 1
mu_B <- mu_base
mu_B[odd]  <- mu_base[odd] * 1.6
mu_B[!odd] <- mu_base[!odd] *
  (sum(mu_base) - sum(mu_base[odd] * 1.6)) / sum(mu_base[!odd])

a_mult    <- unique(round(mu_A / mu_base, 6))
b_up      <- max(mu_B / mu_base)
total_ref <- sum(mu_base)
total_A   <- sum(mu_A)
total_B   <- sum(mu_B)
b_down    <- min(mu_B / mu_base)

Scenario A is more of everything: every taxon is 1.5 times as abundant in the treated plots, so the expected community total rises from 337 to 506 individuals while relative composition is untouched. Scenario B is turnover: half the species rise by 60 per cent, the other half fall to 48 per cent of their reference means, and the expected total stays at 337. A field ecologist would call both a treatment effect. Two of the three indices below disagree.

scen <- data.frame(
  species  = rep(order(order(mu_base)), 3),
  mu       = c(mu_base, mu_A, mu_B),
  scenario = rep(c("reference", "A: all taxa x 1.5", "B: turnover, total held"),
                 each = n_sp))
scen$scenario <- factor(scen$scenario,
                        levels = c("reference", "A: all taxa x 1.5", "B: turnover, total held"))

ggplot(scen, aes(x = species, y = mu, colour = scenario, shape = scenario)) +
  geom_point(size = 2.1, alpha = 0.9) +
  scale_y_log10() +
  scale_colour_manual(values = c("reference" = te_body,
                                 "A: all taxa x 1.5" = te_forest,
                                 "B: turnover, total held" = te_rust)) +
  scale_shape_manual(values = c(16, 17, 15)) +
  labs(x = "species, ordered by reference abundance",
       y = "expected abundance (log scale)", colour = NULL, shape = NULL,
       title = "Two treatment scenarios a field study would both call an effect") +
  theme_datasheet() +
  theme(legend.position = "top")

Three series of points across thirty species ordered by reference abundance. The proportional scenario sits uniformly above the reference; the turnover scenario alternates above and below it. — Figure 1: Expected abundance per species under the reference community and the two treatment scenarios.

Three indices, one design

Bray-Curtis on raw counts, binary Jaccard on presence and absence, and Euclidean distance on Hellinger-transformed counts. All three are standard, all three are defensible, and they answer different questions. Bray-Curtis on raw counts responds to total abundance because it is not closed. The Hellinger transformation divides each plot by its own total, so it discards abundance and keeps proportions. Binary Jaccard discards abundance entirely.

draw <- function(mu, n, size = 5) {
  matrix(rnbinom(n * n_sp, mu = rep(mu, each = n), size = size), n, n_sp)
}

dist_of <- function(m, index) {
  if (index == "Bray-Curtis") vegdist(m, "bray")
  else if (index == "Jaccard (binary)") vegdist(m, "jaccard", binary = TRUE)
  else dist(decostand(m, "hellinger"))
}

one_run <- function(mu_treat, n, index, base = mu_base, size = 5, perms = 199) {
  m <- rbind(draw(base, n, size), draw(mu_treat, n, size))
  g <- factor(rep(c("reference", "treated"), each = n))
  a <- adonis2(dist_of(m, index) ~ g, permutations = perms)
  c(p = a$`Pr(>F)`[1], r2 = a$R2[1])
}

n_grid  <- c(5, 8, 12, 20)
indices <- c("Bray-Curtis", "Jaccard (binary)", "Hellinger")
reps    <- 200

grid <- expand.grid(n = n_grid, index = indices,
                    scenario = c("A: all taxa x 1.5", "B: turnover, total held"),
                    stringsAsFactors = FALSE)
grid$power <- NA_real_
grid$r2    <- NA_real_

set.seed(11)
for (k in seq_len(nrow(grid))) {
  mu_t <- if (grid$scenario[k] == "A: all taxa x 1.5") mu_A else mu_B
  out <- replicate(reps, one_run(mu_t, grid$n[k], grid$index[k]))
  grid$power[k] <- mean(out["p", ] <= 0.05)
  grid$r2[k]    <- median(out["r2", ])
}

pick <- function(sc, ix, nn) {
  grid$power[grid$scenario == sc & grid$index == ix & grid$n == nn]
}
sc_A <- "A: all taxa x 1.5"
sc_B <- "B: turnover, total held"
mc_se <- sqrt(0.05 * 0.95 / reps)

At twenty plots per group, scenario A is detected by Bray-Curtis in 1.000 of runs, by binary Jaccard in 0.065 and by Hellinger in 0.050. The last two numbers are not low power. They are the nominal error rate: those indices see no difference at all, because the difference was removed before the test ran. Eight plots per group would have been ample for the first index and useless for the other two.

Scenario B reverses part of the ranking. Bray-Curtis reaches 0.985 at five plots per group, Hellinger 0.995, and binary Jaccard climbs slowly from 0.120 to 0.720 over the whole grid.

offsets <- c("Bray-Curtis" = -0.45, "Jaccard (binary)" = 0, "Hellinger" = 0.45)
grid$n_off <- grid$n + offsets[grid$index]

ggplot(grid, aes(x = n_off, y = power, colour = index, shape = index)) +
  geom_hline(yintercept = 0.05, linetype = "dotted", colour = te_ink) +
  annotate("text", x = 12, y = 0.115, label = "nominal 0.05",
           colour = te_ink, size = 3.3) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 2.6) +
  facet_wrap(~ scenario) +
  scale_x_continuous(breaks = n_grid) +
  scale_y_continuous(limits = c(0, 1.02)) +
  scale_colour_manual(values = c("Bray-Curtis" = te_forest,
                                 "Jaccard (binary)" = te_rust,
                                 "Hellinger" = te_gold)) +
  scale_shape_manual(values = c(16, 17, 15)) +
  labs(x = "plots per group", y = "power at alpha = 0.05", colour = NULL, shape = NULL,
       title = "The same design, the same effect, three answers") +
  theme_datasheet() +
  theme(legend.position = "top")

Two panels of power curves against plots per group, with the three series offset slightly so they do not overlap and a labelled dotted line at the nominal rate. In the proportional panel the Bray-Curtis curve climbs to the top of the scale while the Jaccard and Hellinger curves stay flat on the nominal line. In the turnover panel Bray-Curtis and Hellinger both sit at the top of the scale across the grid while Jaccard rises steadily from low to high. — Figure 2: Power to detect each scenario at the five per cent level, over 200 simulated studies per cell.

The null case, which is the part that gets skipped

A power curve without a null run is a curve with no zero on it. Here both groups are drawn from the reference community, so every rejection is a false one.

set.seed(303)
null_rate <- sapply(indices, function(ix) {
  mean(replicate(reps, one_run(mu_base, 12, ix)["p"]) <= 0.05)
})
null_rate

     Bray-Curtis Jaccard (binary)        Hellinger 
           0.075            0.045            0.060

At twelve plots per group the three indices reject in 0.075, 0.045 and 0.060 of runs. With 200 replicates the Monte Carlo standard error on a rate near the nominal level is 0.015, so all three sit within sampling error of 0.05. That matters for reading the power table: the flat Jaccard and Hellinger curves in scenario A are exactly this rate, not a weak signal.

null_df <- data.frame(index = factor(names(null_rate), levels = indices),
                      rate = as.numeric(null_rate))
null_df$lo <- pmax(0, null_df$rate - 1.96 * mc_se)
null_df$hi <- null_df$rate + 1.96 * mc_se

ggplot(null_df, aes(x = index, y = rate, fill = index)) +
  geom_col(width = 0.5) +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.12, colour = te_ink) +
  geom_hline(yintercept = 0.05, linetype = "dotted", colour = te_ink) +
  annotate("text", x = 0.62, y = 0.062, label = "nominal", colour = te_ink, size = 3.4) +
  scale_fill_manual(values = c("Bray-Curtis" = te_forest,
                               "Jaccard (binary)" = te_rust,
                               "Hellinger" = te_gold), guide = "none") +
  labs(x = NULL, y = "false positive rate", title = "Calibration check: no difference in the data") +
  theme_datasheet()

Three bars close to the five per cent nominal line, each with a Monte Carlo error bar spanning it. — Figure 3: False positive rate at twelve plots per group when the two groups are drawn from the same community.

Effect size on each scale

The reason the three curves separate is visible in the variance explained. The same twenty-plot study gives a median R-squared of 0.150 for Bray-Curtis under scenario A, against 0.026 for Hellinger. Under scenario B those become 0.256 and 0.239. The index is not filtering noise; it is defining the signal.

Would binary Jaccard do better in a sparser community?

The matrix above is species-rich and nearly filled, so presence and absence carries almost nothing. That is a property of the simulation, not of Jaccard, and it is worth measuring rather than assuming. The same means divided by twelve give a matrix with many structural zeros.

mu_s  <- mu_base / 12
mu_sA <- mu_s * 1.5
mu_sB <- mu_s
mu_sB[odd]  <- mu_s[odd] * 1.6
mu_sB[!odd] <- mu_s[!odd] * (sum(mu_s) - sum(mu_s[odd] * 1.6)) / sum(mu_s[!odd])

set.seed(202)
occupancy <- mean(draw(mu_s, 500) > 0)
sparse <- expand.grid(index = indices, scenario = c(sc_A, sc_B), stringsAsFactors = FALSE)
sparse$power <- NA_real_
for (k in seq_len(nrow(sparse))) {
  mu_t <- if (sparse$scenario[k] == sc_A) mu_sA else mu_sB
  sparse$power[k] <- mean(replicate(reps, one_run(mu_t, 12, sparse$index[k],
                                                  base = mu_s)["p"]) <= 0.05)
}
sparse_jac_A <- sparse$power[sparse$index == "Jaccard (binary)" & sparse$scenario == sc_A]
rich_jac_A   <- pick(sc_A, "Jaccard (binary)", 12)

With 44 per cent of cells occupied, binary Jaccard detects scenario A in 0.205 of runs at twelve plots per group, against 0.075 in the rich matrix. Adding individuals to a sparse community also adds detections, so a presence index stops being blind to abundance. The lesson survives in a weaker form: the index still decides what the study can see, and now sparsity decides how much.

comp <- rbind(
  data.frame(index = grid$index[grid$n == 12], scenario = grid$scenario[grid$n == 12],
             power = grid$power[grid$n == 12], matrix = "rich (few zeros)"),
  data.frame(index = sparse$index, scenario = sparse$scenario,
             power = sparse$power, matrix = "sparse (many zeros)"))
comp$index <- factor(comp$index, levels = indices)

ggplot(comp, aes(x = index, y = power, fill = matrix)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.62) +
  facet_wrap(~ scenario) +
  scale_fill_manual(values = c("rich (few zeros)" = te_forest,
                               "sparse (many zeros)" = te_gold)) +
  scale_y_continuous(limits = c(0, 1.05)) +
  labs(x = NULL, y = "power at twelve plots per group", fill = NULL,
       title = "Sparsity changes what a presence index can detect") +
  theme_datasheet() +
  theme(legend.position = "top", axis.text.x = element_text(angle = 20, hjust = 1))

Grouped bar chart comparing three indices under two scenarios in a rich and a sparse matrix. The Jaccard bars are the only ones that rise substantially when the matrix becomes sparse. — Figure 4: Power at twelve plots per group in the species-rich matrix and in a sparse version of the same community.

Turning this into a design

The machinery is short enough to adapt in an afternoon. Write down the community you expect, write down the change you would not want to miss, choose the index that matches the question you are actually asking, run the grid, and run the null case as well. The choice of index deserves the same paragraph in the methods as the sample size, because the two are not separable: reporting that a study had eighty per cent power says nothing until the index is named.

If the question is about total abundance, an index that normalises by plot totals will not answer it. If the question is about relative composition, an index that responds to total abundance will answer a different question and may answer it very convincingly. The dissimilarity index tutorial sets out which index tracks which gradient; this post prices those properties in plots per group.

Honest limits

Every number here belongs to this generator: thirty species, this abundance distribution, this dispersion parameter, balanced groups, and a treatment applied uniformly across taxa. Real communities have patchy responses, correlated species, and unbalanced designs, all of which move power. The rates carry a Monte Carlo standard error of about 0.035 at the steep part of a curve, so differences of a few percentage points between neighbouring cells are noise.

The deeper limit is that a simulated power analysis needs an assumed effect, and if you knew the effect you would not be running the study. The output is not a prediction; it is a statement of what the design could have found had the assumption been right, which is exactly what a reader needs when the result is not significant.

References

Anderson MJ 2001 Austral Ecology 26(1):32-46 (10.1111/j.1442-9993.2001.01070.pp.x)

Legendre P, Gallagher ED 2001 Oecologia 129(2):271-280 (10.1007/s004420100716)

Warton DI, Wright ST, Wang Y 2012 Methods in Ecology and Evolution 3(1):89-101 (10.1111/j.2041-210X.2011.00127.x)

--- title: "Power for a community study" description: "How many plots does a PERMANOVA need? Simulation shows the answer depends on the dissimilarity index, which decides what counts as a community difference." date: "2026-08-05 11:00" categories: [R, vegan, PERMANOVA, community ecology, study design, ecology tutorial] image: thumbnail.png image-alt: "Power curves against sample size for three dissimilarity indices, showing one index reaching full power while the other two stay flat at the nominal rate." --- Sample size questions have a standard answer for a t-test and no answer at all for a community study. The literature is full of designs with five or eight plots per group, an `adonis2` result, and no statement of what the study could have detected. The gap is not that power analysis is hard here; it is that the first input is missing. Before you can ask how many plots, you have to say what a difference between communities means, and that decision is made by the dissimilarity index, not by the test. This post simulates a PERMANOVA design end to end: a known difference, a sample size grid, three common indices, and a null case as a calibration check. The result is that the power of the same design, on the same data, over the same range of sample sizes, spans the full interval from nothing to certainty depending on which index you picked in the first line of the script. ## A community and two ways of changing it Thirty species with log-normal mean abundances and negative binomial counts, so the mean-variance relationship is realistic rather than Poisson. Two treatment scenarios, both plausible field outcomes. ```{r setup} #| message: false #| warning: false library(vegan) library(ggplot2) 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"), strip.text = element_text(colour = te_ink, face = "bold"), axis.text = element_text(colour = te_body)) } ``` ```{r scenarios} n_sp <- 30 set.seed(20260805) mu_base <- exp(rnorm(n_sp, mean = 1.6, sd = 0.9)) ## scenario A: every taxon multiplied by 1.5, so proportions are untouched mu_A <- mu_base * 1.5 ## scenario B: turnover, with the community total held at the reference total odd <- seq_len(n_sp) %% 2 == 1 mu_B <- mu_base mu_B[odd] <- mu_base[odd] * 1.6 mu_B[!odd] <- mu_base[!odd] * (sum(mu_base) - sum(mu_base[odd] * 1.6)) / sum(mu_base[!odd]) a_mult <- unique(round(mu_A / mu_base, 6)) b_up <- max(mu_B / mu_base) total_ref <- sum(mu_base) total_A <- sum(mu_A) total_B <- sum(mu_B) b_down <- min(mu_B / mu_base) ``` Scenario A is more of everything: every taxon is `r sprintf("%.1f", a_mult)` times as abundant in the treated plots, so the expected community total rises from `r sprintf("%.0f", total_ref)` to `r sprintf("%.0f", total_A)` individuals while relative composition is untouched. Scenario B is turnover: half the species rise by `r sprintf("%.0f", 100 * (b_up - 1))` per cent, the other half fall to `r sprintf("%.0f", 100 * b_down)` per cent of their reference means, and the expected total stays at `r sprintf("%.0f", total_B)`. A field ecologist would call both a treatment effect. Two of the three indices below disagree. ```{r fig-scenarios} #| fig-cap: "Expected abundance per species under the reference community and the two treatment scenarios." #| fig-alt: "Three series of points across thirty species ordered by reference abundance. The proportional scenario sits uniformly above the reference; the turnover scenario alternates above and below it." scen <- data.frame( species = rep(order(order(mu_base)), 3), mu = c(mu_base, mu_A, mu_B), scenario = rep(c("reference", "A: all taxa x 1.5", "B: turnover, total held"), each = n_sp)) scen$scenario <- factor(scen$scenario, levels = c("reference", "A: all taxa x 1.5", "B: turnover, total held")) ggplot(scen, aes(x = species, y = mu, colour = scenario, shape = scenario)) + geom_point(size = 2.1, alpha = 0.9) + scale_y_log10() + scale_colour_manual(values = c("reference" = te_body, "A: all taxa x 1.5" = te_forest, "B: turnover, total held" = te_rust)) + scale_shape_manual(values = c(16, 17, 15)) + labs(x = "species, ordered by reference abundance", y = "expected abundance (log scale)", colour = NULL, shape = NULL, title = "Two treatment scenarios a field study would both call an effect") + theme_datasheet() + theme(legend.position = "top") ``` ## Three indices, one design Bray-Curtis on raw counts, binary Jaccard on presence and absence, and Euclidean distance on Hellinger-transformed counts. All three are standard, all three are defensible, and they answer different questions. Bray-Curtis on raw counts responds to total abundance because it is not closed. The Hellinger transformation divides each plot by its own total, so it discards abundance and keeps proportions. Binary Jaccard discards abundance entirely. ```{r machinery} draw <- function(mu, n, size = 5) { matrix(rnbinom(n * n_sp, mu = rep(mu, each = n), size = size), n, n_sp) } dist_of <- function(m, index) { if (index == "Bray-Curtis") vegdist(m, "bray") else if (index == "Jaccard (binary)") vegdist(m, "jaccard", binary = TRUE) else dist(decostand(m, "hellinger")) } one_run <- function(mu_treat, n, index, base = mu_base, size = 5, perms = 199) { m <- rbind(draw(base, n, size), draw(mu_treat, n, size)) g <- factor(rep(c("reference", "treated"), each = n)) a <- adonis2(dist_of(m, index) ~ g, permutations = perms) c(p = a$`Pr(>F)`[1], r2 = a$R2[1]) } ``` ```{r power-grid} n_grid <- c(5, 8, 12, 20) indices <- c("Bray-Curtis", "Jaccard (binary)", "Hellinger") reps <- 200 grid <- expand.grid(n = n_grid, index = indices, scenario = c("A: all taxa x 1.5", "B: turnover, total held"), stringsAsFactors = FALSE) grid$power <- NA_real_ grid$r2 <- NA_real_ set.seed(11) for (k in seq_len(nrow(grid))) { mu_t <- if (grid$scenario[k] == "A: all taxa x 1.5") mu_A else mu_B out <- replicate(reps, one_run(mu_t, grid$n[k], grid$index[k])) grid$power[k] <- mean(out["p", ] <= 0.05) grid$r2[k] <- median(out["r2", ]) } pick <- function(sc, ix, nn) { grid$power[grid$scenario == sc & grid$index == ix & grid$n == nn] } sc_A <- "A: all taxa x 1.5" sc_B <- "B: turnover, total held" mc_se <- sqrt(0.05 * 0.95 / reps) ``` At twenty plots per group, scenario A is detected by Bray-Curtis in `r sprintf("%.3f", pick(sc_A, "Bray-Curtis", 20))` of runs, by binary Jaccard in `r sprintf("%.3f", pick(sc_A, "Jaccard (binary)", 20))` and by Hellinger in `r sprintf("%.3f", pick(sc_A, "Hellinger", 20))`. The last two numbers are not low power. They are the nominal error rate: those indices see no difference at all, because the difference was removed before the test ran. Eight plots per group would have been ample for the first index and useless for the other two. Scenario B reverses part of the ranking. Bray-Curtis reaches `r sprintf("%.3f", pick(sc_B, "Bray-Curtis", 5))` at five plots per group, Hellinger `r sprintf("%.3f", pick(sc_B, "Hellinger", 5))`, and binary Jaccard climbs slowly from `r sprintf("%.3f", pick(sc_B, "Jaccard (binary)", 5))` to `r sprintf("%.3f", pick(sc_B, "Jaccard (binary)", 20))` over the whole grid. ```{r fig-power} #| fig-cap: "Power to detect each scenario at the five per cent level, over 200 simulated studies per cell." #| fig-alt: "Two panels of power curves against plots per group, with the three series offset slightly so they do not overlap and a labelled dotted line at the nominal rate. In the proportional panel the Bray-Curtis curve climbs to the top of the scale while the Jaccard and Hellinger curves stay flat on the nominal line. In the turnover panel Bray-Curtis and Hellinger both sit at the top of the scale across the grid while Jaccard rises steadily from low to high." offsets <- c("Bray-Curtis" = -0.45, "Jaccard (binary)" = 0, "Hellinger" = 0.45) grid$n_off <- grid$n + offsets[grid$index] ggplot(grid, aes(x = n_off, y = power, colour = index, shape = index)) + geom_hline(yintercept = 0.05, linetype = "dotted", colour = te_ink) + annotate("text", x = 12, y = 0.115, label = "nominal 0.05", colour = te_ink, size = 3.3) + geom_line(linewidth = 0.8) + geom_point(size = 2.6) + facet_wrap(~ scenario) + scale_x_continuous(breaks = n_grid) + scale_y_continuous(limits = c(0, 1.02)) + scale_colour_manual(values = c("Bray-Curtis" = te_forest, "Jaccard (binary)" = te_rust, "Hellinger" = te_gold)) + scale_shape_manual(values = c(16, 17, 15)) + labs(x = "plots per group", y = "power at alpha = 0.05", colour = NULL, shape = NULL, title = "The same design, the same effect, three answers") + theme_datasheet() + theme(legend.position = "top") ``` ## The null case, which is the part that gets skipped A power curve without a null run is a curve with no zero on it. Here both groups are drawn from the reference community, so every rejection is a false one. ```{r null-calibration} set.seed(303) null_rate <- sapply(indices, function(ix) { mean(replicate(reps, one_run(mu_base, 12, ix)["p"]) <= 0.05) }) null_rate ``` At twelve plots per group the three indices reject in `r sprintf("%.3f", null_rate[["Bray-Curtis"]])`, `r sprintf("%.3f", null_rate[["Jaccard (binary)"]])` and `r sprintf("%.3f", null_rate[["Hellinger"]])` of runs. With `r reps` replicates the Monte Carlo standard error on a rate near the nominal level is `r sprintf("%.3f", mc_se)`, so all three sit within sampling error of `r sprintf("%.2f", 0.05)`. That matters for reading the power table: the flat Jaccard and Hellinger curves in scenario A are exactly this rate, not a weak signal. ```{r fig-null} #| fig-cap: "False positive rate at twelve plots per group when the two groups are drawn from the same community." #| fig-alt: "Three bars close to the five per cent nominal line, each with a Monte Carlo error bar spanning it." null_df <- data.frame(index = factor(names(null_rate), levels = indices), rate = as.numeric(null_rate)) null_df$lo <- pmax(0, null_df$rate - 1.96 * mc_se) null_df$hi <- null_df$rate + 1.96 * mc_se ggplot(null_df, aes(x = index, y = rate, fill = index)) + geom_col(width = 0.5) + geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.12, colour = te_ink) + geom_hline(yintercept = 0.05, linetype = "dotted", colour = te_ink) + annotate("text", x = 0.62, y = 0.062, label = "nominal", colour = te_ink, size = 3.4) + scale_fill_manual(values = c("Bray-Curtis" = te_forest, "Jaccard (binary)" = te_rust, "Hellinger" = te_gold), guide = "none") + labs(x = NULL, y = "false positive rate", title = "Calibration check: no difference in the data") + theme_datasheet() ``` ## Effect size on each scale The reason the three curves separate is visible in the variance explained. The same twenty-plot study gives a median R-squared of `r sprintf("%.3f", grid$r2[grid$scenario == sc_A & grid$index == "Bray-Curtis" & grid$n == 20])` for Bray-Curtis under scenario A, against `r sprintf("%.3f", grid$r2[grid$scenario == sc_A & grid$index == "Hellinger" & grid$n == 20])` for Hellinger. Under scenario B those become `r sprintf("%.3f", grid$r2[grid$scenario == sc_B & grid$index == "Bray-Curtis" & grid$n == 20])` and `r sprintf("%.3f", grid$r2[grid$scenario == sc_B & grid$index == "Hellinger" & grid$n == 20])`. The index is not filtering noise; it is defining the signal. ## Would binary Jaccard do better in a sparser community? The matrix above is species-rich and nearly filled, so presence and absence carries almost nothing. That is a property of the simulation, not of Jaccard, and it is worth measuring rather than assuming. The same means divided by twelve give a matrix with many structural zeros. ```{r sparse} mu_s <- mu_base / 12 mu_sA <- mu_s * 1.5 mu_sB <- mu_s mu_sB[odd] <- mu_s[odd] * 1.6 mu_sB[!odd] <- mu_s[!odd] * (sum(mu_s) - sum(mu_s[odd] * 1.6)) / sum(mu_s[!odd]) set.seed(202) occupancy <- mean(draw(mu_s, 500) > 0) sparse <- expand.grid(index = indices, scenario = c(sc_A, sc_B), stringsAsFactors = FALSE) sparse$power <- NA_real_ for (k in seq_len(nrow(sparse))) { mu_t <- if (sparse$scenario[k] == sc_A) mu_sA else mu_sB sparse$power[k] <- mean(replicate(reps, one_run(mu_t, 12, sparse$index[k], base = mu_s)["p"]) <= 0.05) } sparse_jac_A <- sparse$power[sparse$index == "Jaccard (binary)" & sparse$scenario == sc_A] rich_jac_A <- pick(sc_A, "Jaccard (binary)", 12) ``` With `r sprintf("%.0f", 100 * occupancy)` per cent of cells occupied, binary Jaccard detects scenario A in `r sprintf("%.3f", sparse_jac_A)` of runs at twelve plots per group, against `r sprintf("%.3f", rich_jac_A)` in the rich matrix. Adding individuals to a sparse community also adds detections, so a presence index stops being blind to abundance. The lesson survives in a weaker form: the index still decides what the study can see, and now sparsity decides how much. ```{r fig-sparse} #| fig-cap: "Power at twelve plots per group in the species-rich matrix and in a sparse version of the same community." #| fig-alt: "Grouped bar chart comparing three indices under two scenarios in a rich and a sparse matrix. The Jaccard bars are the only ones that rise substantially when the matrix becomes sparse." comp <- rbind( data.frame(index = grid$index[grid$n == 12], scenario = grid$scenario[grid$n == 12], power = grid$power[grid$n == 12], matrix = "rich (few zeros)"), data.frame(index = sparse$index, scenario = sparse$scenario, power = sparse$power, matrix = "sparse (many zeros)")) comp$index <- factor(comp$index, levels = indices) ggplot(comp, aes(x = index, y = power, fill = matrix)) + geom_col(position = position_dodge(width = 0.7), width = 0.62) + facet_wrap(~ scenario) + scale_fill_manual(values = c("rich (few zeros)" = te_forest, "sparse (many zeros)" = te_gold)) + scale_y_continuous(limits = c(0, 1.05)) + labs(x = NULL, y = "power at twelve plots per group", fill = NULL, title = "Sparsity changes what a presence index can detect") + theme_datasheet() + theme(legend.position = "top", axis.text.x = element_text(angle = 20, hjust = 1)) ``` ## Turning this into a design The machinery is short enough to adapt in an afternoon. Write down the community you expect, write down the change you would not want to miss, choose the index that matches the question you are actually asking, run the grid, and run the null case as well. The choice of index deserves the same paragraph in the methods as the sample size, because the two are not separable: reporting that a study had eighty per cent power says nothing until the index is named. If the question is about total abundance, an index that normalises by plot totals will not answer it. If the question is about relative composition, an index that responds to total abundance will answer a different question and may answer it very convincingly. The [dissimilarity index tutorial](../choosing-a-dissimilarity-index/) sets out which index tracks which gradient; this post prices those properties in plots per group. ## Honest limits Every number here belongs to this generator: thirty species, this abundance distribution, this dispersion parameter, balanced groups, and a treatment applied uniformly across taxa. Real communities have patchy responses, correlated species, and unbalanced designs, all of which move power. The rates carry a Monte Carlo standard error of about `r sprintf("%.3f", sqrt(0.25 / reps))` at the steep part of a curve, so differences of a few percentage points between neighbouring cells are noise. The deeper limit is that a simulated power analysis needs an assumed effect, and if you knew the effect you would not be running the study. The output is not a prediction; it is a statement of what the design could have found had the assumption been right, which is exactly what a reader needs when the result is not significant. ## References Anderson MJ 2001 Austral Ecology 26(1):32-46 (10.1111/j.1442-9993.2001.01070.pp.x) Legendre P, Gallagher ED 2001 Oecologia 129(2):271-280 (10.1007/s004420100716) Warton DI, Wright ST, Wang Y 2012 Methods in Ecology and Evolution 3(1):89-101 (10.1111/j.2041-210X.2011.00127.x) ## Related tutorials - [Choosing a dissimilarity index](../choosing-a-dissimilarity-index/) - [Power analysis by simulation](../power-analysis-by-simulation/) - [Common PERMANOVA mistakes](../common-permanova-mistakes/) - [Pairwise PERMANOVA](../pairwise-permanova/)