---
title: "Fitting species abundance distributions in R"
description: "Fit and compare rank-abundance models in R with vegan radfit: pre-emption, lognormal, Zipf and Zipf-Mandelbrot, with AIC selection and Preston octave plots."
date: "2026-05-05 09:00"
date-modified: "2026-10-07"
categories: [R, community ecology, biodiversity, vegan, ecology tutorial]
image: thumbnail.png
image-alt: "Rank-abundance curves for two ecological communities with fitted SAD models on a log abundance axis"
---
*Corrected 7 October 2026: the introduction said the two example communities have similar species lists, though they are unnamed synthetic counts that differ widely in richness, and now describes them as one dominated by a few species and one more even; the octave section said the forest peak sits exactly where a lognormal puts its mode, and now says close to where the fitted lognormal puts it.*
*Corrected 30 September 2026: the octave plot section described Preston's octaves as the classes used in the plots; it now says that Preston's class limits double and that he split species sitting on a limit between neighbouring octaves, while the plots here use simpler classes without that split.*
*Updated 27 September 2026: a new section, The best model changes with the area pooled, fits the five models to the Barro Colorado Island tree census at 1, 5, 25 and 50 hectares and counts how often each one wins, and the practical note on unequal effort now says rarefy rather than rescale.*
A diversity index turns a whole community into one number. Shannon, Simpson and evenness are useful, but they throw away the shape of the abundance distribution: how steeply the common species tower over the rare ones, and whether there is a mode in the middle. That shape is the species abundance distribution (SAD), one of the oldest descriptions in community ecology (Fisher et al. 1943; Preston 1948).
This post fits parametric SAD models to ranked counts with `vegan::radfit`, compares them with AIC, and shows how two communities, one dominated by a few species and one more even, end up under different models. The point is not that a fitted model reveals a mechanism. It is that the rank-abundance curve carries structure that a single index hides, and that model selection gives a compact, comparable summary of that structure.
## Two synthetic communities
The data are illustrative. One community has strong dominance, like an early-successional site where one or two species pre-empt most of the available space. The other is more even, with many species of intermediate abundance and a few rare ones, closer to a mature stand.
```{r}
#| label: data
#| message: false
#| warning: false
library(vegan)
library(ggplot2)
library(dplyr)
library(tidyr)
# Quarry: early-successional, niche pre-emption (geometric dominance).
# Each species takes a fixed fraction k of the resource left by the previous one.
set.seed(19)
S1 <- 15; k <- 0.38; Ntot1 <- 2600
rem <- 1; props <- numeric(S1)
for (i in 1:S1) { props[i] <- k * rem; rem <- rem - props[i] }
props <- props / sum(props)
quarry <- rpois(S1, props * Ntot1); quarry <- quarry[quarry > 0]
# Forest: mature, lognormal SAD (more even, with an internal mode in log abundance).
set.seed(204)
raw <- rlnorm(55, meanlog = 1.7, sdlog = 1.05)
forest <- round(raw * 4200 / sum(raw)); forest <- forest[forest > 0]
c(quarry_S = length(quarry), quarry_N = sum(quarry),
forest_S = length(forest), forest_N = sum(forest))
```
The quarry holds `r length(quarry)` species in `r sum(quarry)` individuals; the forest holds `r length(forest)` species in `r sum(forest)`. Sorted from most to least abundant, the quarry falls away fast (951 down to a single individual), while the forest descends more gently.
## What the indices say, and what they hide
```{r}
#| label: indices
idx <- function(x) {
H <- diversity(x, "shannon")
data.frame(
S = sum(x > 0), N = sum(x),
Shannon = round(H, 2),
Simpson = round(diversity(x, "simpson"), 3),
Pielou_J = round(H / log(sum(x > 0)), 2),
Hill_q1 = round(exp(H), 1))
}
knitr::kable(rbind(Quarry = idx(quarry), Forest = idx(forest)))
```
The forest is more diverse and far more even (Pielou's J of `r sprintf("%.2f", idx(forest)$Pielou_J)` against `r sprintf("%.2f", idx(quarry)$Pielou_J)`). These numbers tell you the forest is flatter, but they do not tell you the form of the flattening: whether the abundances decline geometrically, follow a lognormal hump, or fall on a power law. For that you need the curve itself.
## Rank-abundance curves and the radfit models
A rank-abundance (dominance/diversity) plot puts species rank on the x-axis and abundance on a log y-axis (Whittaker 1965). `radfit` fits five models to these ranked counts as Poisson generalised linear models, choosing parameters by minimising deviance (Wilson 1991):
- **Null** (broken stick): no free parameters, a neutral partition of the resource axis.
- **Pre-emption** (geometric series): one parameter, a constant fraction taken at each rank. A straight line on the log axis.
- **Lognormal**: two parameters, a hump in log abundance.
- **Zipf**: two parameters, a power law with a heavy tail of rare species.
- **Zipf-Mandelbrot**: three parameters, a more flexible version of Zipf.
```{r}
#| label: radfit
#| message: false
#| warning: false
rq <- radfit(quarry)
rf <- radfit(forest)
rq
rf
```
```{r}
#| label: fig-rad
#| fig-width: 8.6
#| fig-height: 4.2
#| fig-cap: "Rank-abundance curves with three fitted SAD models. The quarry follows a straight pre-emption line on the log axis; the forest bends away from it toward a lognormal."
#| fig-alt: "Rank-abundance curves for two communities on a log axis with fitted SAD models: the quarry follows a straight pre-emption line while the forest bends toward a lognormal."
ink <- "#16241d"; forest_col <- "#275139"; gold <- "#cda23f"; red <- "#b5534e"
faint <- "#5d6b61"; line_col <- "#dad9ca"; paper <- "#f5f4ee"
fitcurve <- function(r, model, site) {
f <- r$models[[model]]$fitted.values
data.frame(rank = seq_along(f), fit = f, model = model, site = site)
}
mods <- c("Preemption", "Lognormal", "Zipf")
fit <- bind_rows(
lapply(mods, fitcurve, r = rq, site = "Quarry (pre-emption)"),
lapply(mods, fitcurve, r = rf, site = "Forest (lognormal)"))
obs <- bind_rows(
data.frame(rank = seq_along(quarry), abund = sort(quarry, decreasing = TRUE),
site = "Quarry (pre-emption)"),
data.frame(rank = seq_along(forest), abund = sort(forest, decreasing = TRUE),
site = "Forest (lognormal)"))
lev <- c("Quarry (pre-emption)", "Forest (lognormal)")
fit$site <- factor(fit$site, lev); obs$site <- factor(obs$site, lev)
fit$model <- factor(fit$model, levels = mods)
te <- theme_minimal(base_size = 12) + theme(
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = line_col, linewidth = 0.3),
plot.background = element_rect(fill = paper, color = NA),
panel.background = element_rect(fill = paper, color = NA),
strip.text = element_text(face = "bold", color = ink, hjust = 0),
axis.title = element_text(color = "#46604a"),
legend.position = "top", legend.title = element_blank())
ggplot() +
geom_line(data = fit, aes(rank, fit, color = model), linewidth = 0.9) +
geom_point(data = obs, aes(rank, abund), color = ink, size = 1.6) +
facet_wrap(~site, scales = "free_x") +
scale_y_log10(labels = scales::label_number()) +
scale_color_manual(values = c(Preemption = forest_col, Lognormal = gold, Zipf = red)) +
labs(x = "Species rank (most to least abundant)", y = "Abundance (log scale)") + te
```
For the quarry, pre-emption wins outright: an AIC of `r sprintf("%.1f", AIC(rq$models$Preemption))` against `r sprintf("%.1f", AIC(rq$models$Lognormal))` for the lognormal and `r sprintf("%.1f", AIC(rq$models$Null))` for the null. The fitted fraction is about `r sprintf("%.2f", coef(rq$models$Preemption))`, and that is the share of the resource still unclaimed that each species takes, so each species in rank order holds roughly `r sprintf("%.2f", 1 - coef(rq$models$Preemption))` of what the one above it held. On the log axis that is a straight line, and the points sit on it. For the forest the order reverses: the lognormal wins (AIC `r sprintf("%.1f", AIC(rf$models$Lognormal))`), the pre-emption model is the worst fit (`r sprintf("%.1f", AIC(rf$models$Preemption))`), and the points trace the gentle S-curve that a geometric series cannot follow.
## The Preston octave plot
Preston (1948) binned species into octaves, abundance classes whose limits double (1, 2, 4, 8 and so on), and split the species whose count falls exactly on a limit equally between the two neighbouring octaves. The plots here use the simpler classes 1, 2 to 3, 4 to 7 and so on, with no splitting, which keeps all the singletons in the lowest octave. A lognormal community shows an internal mode: a peak at some intermediate octave, with fewer species both rarer and commoner. A geometric community has no such peak; species are spread roughly evenly across octaves.
```{r}
#| label: fig-octave
#| fig-width: 8.6
#| fig-height: 3.6
#| fig-cap: "Preston octave plots. The forest has a clear internal mode near octave 6, the signature of a lognormal; the quarry is flat across octaves."
#| fig-alt: "Preston octave plots of species counts: the forest shows an internal peak near octave six while the quarry stays flat across octaves."
octdf <- function(x, site) {
cl <- floor(log2(sort(x, decreasing = TRUE)))
tb <- tabulate(cl + 1, nbins = max(cl) + 1)
data.frame(octave = seq_along(tb) - 1, nsp = tb, site = site)
}
oct <- bind_rows(octdf(quarry, "Quarry (pre-emption)"),
octdf(forest, "Forest (lognormal)"))
oct$site <- factor(oct$site, lev)
ggplot(oct, aes(octave, nsp)) +
geom_col(fill = forest_col, width = 0.85) +
facet_wrap(~site) +
scale_x_continuous(breaks = 0:9) +
labs(x = expression(paste("Abundance octave (", log[2], " class)")),
y = "Number of species") + te
```
The forest piles up near octave 6, close to where the fitted lognormal puts its mode. The quarry has no peak, which is why the lognormal fits it poorly and the geometric line fits it well.
## Model selection across all five models
Ranking by AIC makes the contrast explicit. The best model has a delta of zero; everything else is measured against it.
```{r}
#| label: fig-aic
#| fig-width: 8.6
#| fig-height: 3.4
#| fig-cap: "Delta-AIC by model. Pre-emption is best for the quarry, lognormal for the forest. Mandelbrot is runner-up in both, about four AIC units behind in the quarry (there the gap is the parameter penalty alone) but far behind in the forest, where it also fits worse."
#| fig-alt: "Bar chart of delta-AIC by model, with pre-emption lowest for the quarry and lognormal lowest for the forest, and Mandelbrot runner-up in both, close behind in the quarry and far behind in the forest."
getaic <- function(r) sapply(r$models, function(m) tryCatch(AIC(m), error = function(e) NA))
aicdf <- bind_rows(
data.frame(model = names(getaic(rq)), aic = getaic(rq), site = "Quarry (pre-emption)"),
data.frame(model = names(getaic(rf)), aic = getaic(rf), site = "Forest (lognormal)")) %>%
group_by(site) %>% mutate(dAIC = aic - min(aic, na.rm = TRUE), best = dAIC == 0) %>% ungroup()
aicdf$site <- factor(aicdf$site, lev)
aicdf$model <- factor(aicdf$model,
levels = rev(c("Null", "Preemption", "Lognormal", "Zipf", "Mandelbrot")))
ggplot(aicdf, aes(dAIC, model)) +
geom_segment(aes(x = 0, xend = dAIC, yend = model), color = line_col, linewidth = 1) +
geom_point(aes(color = best), size = 3) +
facet_wrap(~site, scales = "free_x") +
scale_color_manual(values = c(`FALSE` = faint, `TRUE` = forest_col), guide = "none") +
labs(x = expression(paste(Delta, "AIC (0 = best model)")), y = NULL) + te
```
One detail rewards attention. For the quarry, Mandelbrot fits the data exactly as well as pre-emption: both reach a deviance of about `r sprintf("%.1f", rq$models$Preemption$deviance)`. But Mandelbrot carries three parameters against the one of pre-emption, so its AIC is `r sprintf("%.1f", AIC(rq$models$Mandelbrot))` against `r sprintf("%.1f", AIC(rq$models$Preemption))`, a delta of 4: two parameters at two AIC units each. Equal fit, worse score. This is parsimony in action: when two models describe the data equally well, the simpler one is preferred, and AIC does the bookkeeping. The extra flexibility of Mandelbrot only pays off in the forest, where the lognormal still beats it.
## The best model changes with the area pooled
Each community above is one sample, fitted once. A real survey can pool its plots in more than one way, and the winner depends on how many individuals go into the fit. The Barro Colorado Island tree data that ship with vegan (`BCI`) show this on real counts: 50 one-hectare quadrats that together make a 50-hectare plot, with every tree of at least 10 cm diameter at breast height tallied in each. The chunk fits the five models to every quadrat on its own, to 40 random sets of 5 quadrats, to 40 random sets of 25, and to the whole plot, and keeps the model with the lowest AIC, the criterion used above (`AIC()` on a `radfit` object returns the same five values as `getaic()`, which the chunk checks). The pool sizes, the 40 pools and the seed were fixed before anything was run.
```{r}
#| label: bci-pools
data("BCI", package = "vegan")
stopifnot(dim(BCI) == c(50, 225), sum(BCI) == 21457, all(colSums(BCI) > 0))
# fit the five models to one count vector; warnings are counted and muffled,
# and the final Mandelbrot model is refitted to check that it converged
bci_fit <- function(x) {
x <- x[x > 0]; msg <- character(0); n_other <- 0
r <- withCallingHandlers(radfit(x), warning = function(w) {
in_zm <- any(vapply(sys.calls(), function(cl) identical(cl[[1]], quote(rad.zipfbrot)), TRUE))
msg <<- c(msg, conditionMessage(w)); n_other <<- n_other + !in_zm
invokeRestart("muffleWarning") })
zm <- r$models$Mandelbrot; zm_rank <- seq_along(zm$y)
zm_glm <- glm(zm$y ~ log(zm_rank + zm$coefficients[["beta"]]), family = poisson,
offset = rep(log(sum(zm$y)), length(zm$y)))
stopifnot(n_other == 0, zm_glm$converged, isTRUE(all.equal(AIC(r), getaic(r))),
abs(AIC(zm) - (AIC(zm_glm) + 2)) < 1e-6)
bci_oct <- octdf(x, "")
data.frame(S = length(x), winner = names(which.min(AIC(r))), warned = length(msg) > 0,
mode_oct = bci_oct$octave[which.max(bci_oct$nsp)], msg = paste(unique(msg), collapse = "; "))
}
set.seed(1948)
bci_pools <- list(`1 ha` = lapply(1:50, function(i) unlist(BCI[i, ])),
`5 ha` = replicate(40, colSums(BCI[sample(50, 5), ]), simplify = FALSE),
`25 ha` = replicate(40, colSums(BCI[sample(50, 25), ]), simplify = FALSE),
`50 ha` = list(colSums(BCI)))
bci_res <- bind_rows(lapply(bci_pools, function(p) bind_rows(lapply(p, bci_fit))), .id = "area")
bci_res$area <- factor(bci_res$area, names(bci_pools))
bci_tab <- bci_res %>% group_by(area) %>%
summarise(fits = n(), median_S = median(S), Lognormal = sum(winner == "Lognormal"),
Mandelbrot = sum(winner == "Mandelbrot"), Zipf = sum(winner == "Zipf"),
other = sum(!winner %in% c("Lognormal", "Mandelbrot", "Zipf")),
mode_above_0 = sum(mode_oct > 0), hump_LN = sum(mode_oct > 0 & winner == "Lognormal"),
hump_ZM = sum(mode_oct > 0 & winner == "Mandelbrot"), warned = sum(warned))
knitr::kable(bci_tab)
bt <- function(a, col) bci_tab[[col]][bci_tab$area == a]
stopifnot(bt("1 ha", "Lognormal") == 0, bt("50 ha", "Lognormal") == 1)
bci_msg <- unique(bci_res$msg[bci_res$warned])
stopifnot(sum(bci_tab$warned) > 0, length(bci_msg) == 1)
```
On a single hectare Mandelbrot wins in `r bt("1 ha", "Mandelbrot")` of the 50 quadrats and Zipf in `r bt("1 ha", "Zipf")`; the lognormal wins none. Pool 5 quadrats and the two split the 40 pools almost evenly, `r bt("5 ha", "Lognormal")` for the lognormal and `r bt("5 ha", "Mandelbrot")` for Mandelbrot; pool 25 and the lognormal wins `r bt("25 ha", "Lognormal")` of 40, and it wins for the whole plot too. The forest is the same in every row; what changes is the number of trees read from it, and with it the number of species, from a median of `r bt("1 ha", "median_S")` in a quadrat to `r bt("50 ha", "median_S")` in the whole plot.
Preston's (1948) veil is the textbook reading of this drift. A sample small enough to miss many of the rare species cuts the octave plot on its left: the species it does see pile up in the lowest octaves, the curve falls from the first octave, and a heavy tail of rare species is what Zipf and Mandelbrot are built to fit. As the sample grows every count grows with it, the curve moves right past the veil line, and the hump of the lognormal comes into view. The `mode_above_0` column of the table follows it: the modal octave lies above octave 0 in `r bt("1 ha", "mode_above_0")` of 50 single quadrats, `r bt("5 ha", "mode_above_0")` of 40 five-quadrat pools and `r bt("25 ha", "mode_above_0")` of 40 pools of 25. The column moves with the totals, not pool by pool: among the five-quadrat pools whose mode lies above octave 0 the lognormal wins `r bt("5 ha", "hump_LN")` and Mandelbrot `r bt("5 ha", "hump_ZM")` (the `hump_LN` and `hump_ZM` columns), so the veil reading fits the trend across areas, while the winner of a single pool of a given size turns on more than where its mode sits.
```{r}
#| label: fig-bci-octaves
#| fig-width: 8.6
#| fig-height: 3.6
#| fig-cap: "Preston octave plots of the Barro Colorado Island trees, averaged over the pools at each area (the whole plot is a single pool). A single hectare peaks in the lowest octave; as more quadrats are pooled the peak moves to higher octaves and a hump appears."
#| fig-alt: "Four bar charts of the mean number of species per abundance octave, for 1, 5, 25 and 50 hectares: at 1 hectare the bars fall from a peak at octave 0, at 5 hectares the peak sits at octave 1, and at 25 and 50 hectares the bars rise to a hump near octaves 3 and 4 before falling away."
bci_oct_all <- bind_rows(lapply(names(bci_pools), function(a)
bind_rows(lapply(bci_pools[[a]], function(x) octdf(x[x > 0], a))) %>%
complete(octave = 0:max(octave), site, fill = list(nsp = 0)) %>%
group_by(site, octave) %>% summarise(nsp = sum(nsp) / length(bci_pools[[a]]), .groups = "drop")))
bci_oct_all$site <- factor(bci_oct_all$site, names(bci_pools))
ggplot(bci_oct_all, aes(octave, nsp)) +
geom_col(fill = forest_col, width = 0.85) +
facet_wrap(~site, nrow = 1) +
scale_x_continuous(breaks = seq(0, 10, 2)) +
labs(x = expression(paste("Abundance octave (", log[2], " class)")),
y = "Mean number of species") + te
```
`radfit` warned in `r sum(bci_tab$warned)` of the `r sum(bci_tab$fits)` fits, with the message "`r paste(bci_msg, collapse = "; ")`". Every warning came from inside the Mandelbrot search over its third parameter (the chunk stops if one comes from another model), and the chunk refits each final Mandelbrot model and stops unless it converged, so the warnings came from trial values during the search and the counts use every fit. The counts are for these trees, one census of stems of at least 10 cm, and these pools: two sets of 25 quadrats drawn from 50 share half their quadrats on average (25 x 25 / 50), so the 40 pools are not independent, and a random set of quadrats is not a contiguous block of forest. None of the five models is the right one for Barro Colorado Island. The winner describes the sample as much as the forest, so compare winning models only between samples with similar numbers of individuals, the same trap as [Pielou's evenness under unequal sampling effort](../pielou-evenness-and-sampling-effort/).
## How far to read these fits
A fitted SAD model is a description, not a diagnosis. Many processes can generate the same curve, and the SAD on its own rarely discriminates among them (McGill et al. 2007). A geometric fit is consistent with sequential niche pre-emption, but it does not prove it; a lognormal can arise from the central limit theorem acting on multiplicative growth, with no special ecology required. Treat the winning model as a compact summary of community structure that lets you compare sites, not as evidence for a mechanism.
A few practical notes. The models are fitted to ranked counts, so they describe relative abundance, not absolute density; to compare communities sampled at different efforts, rarefy them to a common number of individuals (`vegan::rrarefy()`) rather than rescaling, because a rescaled count is no longer a count and does not bring back the species the smaller sample missed (the section above shows what sample size does to the winner). With strongly geometric data the three-parameter Mandelbrot model can fail to converge and `radfit` returns `NA` for it, which is informative rather than alarming: the data do not support the extra parameters. And `radfit` works on a single community vector here, but it also accepts a site-by-species matrix and fits every row at once, which is the usual way to compare many plots.
## References
Fisher RA, Corbet AS, Williams CB 1943 Journal of Animal Ecology 12(1):42-58 (10.2307/1411)
Preston FW 1948 Ecology 29(3):254-283 (10.2307/1930989)
Whittaker RH 1965 Science 147(3655):250-260 (10.1126/science.147.3655.250)
Wilson JB 1991 Journal of Vegetation Science 2(1):35-46 (10.2307/3235896)
McGill BJ, Etienne RS, Gray JS, et al. 2007 Ecology Letters 10(10):995-1015 (10.1111/j.1461-0248.2007.01094.x)
Oksanen J, Simpson G, Blanchet F, et al. 2022 vegan: Community Ecology Package, R package version 2.6-4
## Related tutorials
- [Rarefaction and accumulation curves](../rarefaction-accumulation-curves/)
- [Functional diversity](../functional-diversity-fd/)
- [When not to use Shannon](../when-not-to-use-shannon/)
- [The exact Ewens test: is your SAD neutral?](../testing-neutrality-with-the-ewens-test/)
