Fitting a body size spectrum

body size

power laws

statistics

ecology tutorial

Log bins shift the fitted size spectrum exponent by exactly one, and the likelihood never sees a bin. Fitting individual body mass data in R, no histogram.

Author

Tidy Ecology

Published

2026-08-08

Modified

2026-10-01

Corrected 1 October 2026: the post said the four binned quantities give four different slopes; counts per class and biomass per unit mass share one, so there are three, and What to report now says so.

Corrected 27 September 2026: Honest limits said the example retention curve reaches one well below the chosen limit; it now says the curve is still short of one at the limits the scan chooses. The generalised Pareto remark now says that alpha - 1, not alpha, is the reciprocal of the shape parameter.

Updated 25 September 2026: the post fitted the exponent but did not say how to test it. A new section, Testing the exponent when the lower limit is estimated, gives the exact test when the limit is the sample minimum and measures how often the usual test rejects a true exponent once the limit comes from a scan, with and without a bootstrap that re-runs the scan. What to report now takes the standard error from that bootstrap when the limit is scanned.

Weigh every individual in a benthic sample and the masses span four or five orders of magnitude, with thousands of tiny animals and a handful of large ones. The standard summary of that spread is the size spectrum: a claim that the number of individuals at mass m falls off as a power of m, and a single exponent that says how fast.

The exponent gets compared across sites, across years and against fishing pressure, so it matters that two people analysing the same sample arrive at the same number. They frequently do not, and the reason is rarely the ecology. It is that the exponent depends on whether the counts were divided by the width of the bin they fell in, on how many bins there were, and on where the smallest mass in the fit was set. None of those is a property of the animals.

This post separates the three. The bin width turns out to shift the fitted slope by exactly one, the number of bins turns out to bias it in a direction that can be measured, and the maximum likelihood estimator turns out not to need bins at all.

A sample of individual masses

Individual dry masses in milligrams, drawn from a bounded power law so that the answer is known. The density is proportional to m^b between a lower and an upper mass, and b is the quantity every method below is trying to recover.

library(ggplot2)
library(patchwork)

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))
}

rbounded_pareto <- function(n, b, lo, hi) {
  u <- runif(n)
  ((hi^(b + 1) - lo^(b + 1)) * u + lo^(b + 1))^(1 / (b + 1))
}

b_true <- -2.05
m_lo   <- 0.5
m_hi   <- 8000
n_ind  <- 1200

set.seed(11)
mass <- rbounded_pareto(n_ind, b_true, m_lo, m_hi)
mass_max <- max(mass)

The sample holds 1200 individuals between 0.50 and 285 milligrams. Nothing in the analysis below uses b_true again until the end.

The bin width is not a plotting choice

Bins that double in width are the convention, because equal bins on the arithmetic scale put almost every animal in the first one. The bookkeeping question is what goes on the vertical axis: the raw count in each doubling class, or that count divided by the width of the class.

brk <- 2^(floor(log2(m_lo)):ceiling(log2(mass_max)))
cls <- cut(mass, brk, right = FALSE)

bm_in <- as.numeric(tapply(mass, cls, sum))
bm_in[is.na(bm_in)] <- 0

sp <- data.frame(
  lower   = brk[-length(brk)],
  upper   = brk[-1],
  n_in    = as.numeric(table(cls)),
  biomass = bm_in
)
sp$width <- sp$upper - sp$lower
sp$mid   <- sqrt(sp$lower * sp$upper)
sp$dens  <- sp$n_in / sp$width
sp       <- sp[sp$n_in > 0, ]

fit_count <- lm(log(n_in) ~ log(mid), data = sp)
fit_dens  <- lm(log(dens) ~ log(mid), data = sp)

slope_count <- unname(coef(fit_count)[2])
slope_dens  <- unname(coef(fit_dens)[2])
slope_shift <- slope_count - slope_dens
n_bin_used  <- nrow(sp)

The two lines are fitted to the same 10 bins holding the same 1200 animals. Counting gives a slope of -0.990. Dividing by bin width gives -1.990. The gap between them is 1.0000, and it is not approximately one, it is one.

The reason is that a doubling class starting at mass a has width a, so its width is proportional to its own position on the axis. Dividing the count by a subtracts log(a) from the vertical coordinate while the horizontal coordinate stays at log(a), which tilts the whole line by exactly one unit of slope. Any binning scheme with widths proportional to mass does this; the factor is a property of the bins, not of the sample.

pl <- ggplot(sp, aes(mid, n_in)) +
  geom_point(colour = te_forest, size = 2.4) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              colour = te_rust, linewidth = 0.7) +
  scale_x_log10() + scale_y_log10() +
  labs(x = "body mass (mg)", y = "individuals in class",
       title = "Counts per class",
       subtitle = sprintf("fitted slope %.2f", slope_count)) +
  theme_datasheet()

pr <- ggplot(sp, aes(mid, dens)) +
  geom_point(colour = te_forest, size = 2.4) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              colour = te_rust, linewidth = 0.7) +
  scale_x_log10() + scale_y_log10(labels = scales::label_comma(drop0trailing = TRUE)) +
  labs(x = "body mass (mg)", y = "individuals per mg",
       title = "Counts per unit mass",
       subtitle = sprintf("fitted slope %.2f", slope_dens)) +
  theme_datasheet()

pl + pr + plot_annotation(theme = theme_datasheet())

Two log-log scatter panels side by side. The left panel plots raw counts per doubling class against class midpoint with a fitted line of slope near minus one. The right panel plots counts divided by class width against the same midpoints, with a visibly steeper fitted line of slope near minus two. — Figure 1: One sample of individual masses, binned once, plotted two ways.

This is the same bookkeeping that trips up dispersal kernels, where the density of distances carries a ring area that the density of positions does not, and the same failure to divide by the width of the measure produces a scale parameter that is wrong by a fixed factor.

The bookkeeping continues into biomass. Total mass in a doubling class rises one further power of m above the count, and dividing that by class width takes one back off, so four quantities from one sample give only three different slopes in expectation: biomass per unit mass shares the slope of counts per class.

slope_bm      <- unname(coef(lm(log(biomass) ~ log(mid), data = sp))[2])
slope_bm_dens <- unname(coef(lm(log(biomass / width) ~ log(mid), data = sp))[2])

Counts give -0.99, counts per unit mass give -1.99, biomass gives -0.02 and biomass per unit mass gives -1.02. A published slope near -1 and a published slope near -2 can describe an identical sample. The number is uninterpretable without the sentence that says which of the four it is, and that sentence is missing from a lot of papers.

The likelihood never sees a bin

If the masses really do follow a power law above some lower limit, the exponent has a closed form. Writing the density as proportional to m^b with b below minus one, and setting alpha = -b, the maximum likelihood estimate from n observed masses above a lower limit m0 is a single expression, and so is its standard error.

pareto_alpha <- function(x, m0) 1 + length(x) / sum(log(x / m0))

alpha_hat <- pareto_alpha(mass, m_lo)
b_hat     <- -alpha_hat
alpha_se  <- (alpha_hat - 1) / sqrt(length(mass))

nll_bounded <- function(b, x, lo, hi) {
  -sum(log((b + 1) * x^b / (hi^(b + 1) - lo^(b + 1))))
}
b_bounded <- optimize(nll_bounded, c(-4, -1.01),
                      x = mass, lo = m_lo, hi = m_hi)$minimum

The estimate is -2.049 with a standard error of 0.030, against a truth of -2.05. Allowing for the upper bound leaves it at -2.049: the largest animal in the tray sits far below the bound, so that refinement has nothing to correct here.

Two things are worth noticing about that expression. It contains no bins, no midpoints and no regression, so the choice of doubling class or decade class or any other class cannot influence it. And it is the same estimator that sits behind the generalised Pareto fit used for extremes: a Pareto tail is a generalised Pareto with the scale tied to the threshold, and alpha - 1 is the reciprocal of the generalised Pareto shape parameter.

What the bins cost

Comparing estimators needs repetition rather than one sample, so here are four hundred fresh samples at two sample sizes, each fitted by the closed form and by a log-binned regression at three bin counts.

n_rep <- 400

binned_b <- function(x, m0, n_bin) {
  bb <- exp(seq(log(m0), log(max(x) * 1.0001), length.out = n_bin + 1))
  ct <- as.numeric(table(cut(x, bb, right = FALSE)))
  wd <- diff(bb)
  md <- sqrt(bb[-length(bb)] * bb[-1])
  keep <- ct > 0
  unname(coef(lm(log(ct[keep] / wd[keep]) ~ log(md[keep])))[2])
}

n_small <- 200
bin_grid <- c(6, 10, 20)

set.seed(202)
sweep_one <- function(n_draw) {
  out <- t(replicate(n_rep, {
    x <- rbounded_pareto(n_draw, b_true, m_lo, m_hi)
    c(-pareto_alpha(x, m_lo), vapply(bin_grid, function(k) binned_b(x, m_lo, k), 0))
  }))
  colnames(out) <- c("likelihood", paste0("bins", bin_grid))
  data.frame(
    method = colnames(out),
    n_draw = n_draw,
    bias   = colMeans(out) - b_true,
    rmse   = sqrt(colMeans((out - b_true)^2)),
    mc_se  = apply(out, 2, sd) / sqrt(n_rep),
    row.names = NULL
  )
}
mc <- rbind(sweep_one(n_small), sweep_one(n_ind))

mle_rmse_small <- mc$rmse[mc$method == "likelihood" & mc$n_draw == n_small]
b20_rmse_small <- mc$rmse[mc$method == "bins20"     & mc$n_draw == n_small]
b20_bias_small <- mc$bias[mc$method == "bins20"     & mc$n_draw == n_small]
b06_bias_small <- mc$bias[mc$method == "bins6"      & mc$n_draw == n_small]
mle_bias_big   <- mc$bias[mc$method == "likelihood" & mc$n_draw == n_ind]
b20_bias_big   <- mc$bias[mc$method == "bins20"     & mc$n_draw == n_ind]
rmse_ratio     <- b20_rmse_small / mle_rmse_small
mc_se_max      <- max(mc$mc_se)
b06_bias_big   <- mc$bias[mc$method == "bins6"  & mc$n_draw == n_ind]
b10_bias_big   <- mc$bias[mc$method == "bins10" & mc$n_draw == n_ind]
bin_gap_big    <- b10_bias_big - b06_bias_big

The closed form is essentially unbiased at both sample sizes, with a bias of -0.002 at 1200 individuals. Every binned fit is biased towards a shallower spectrum, the bias grows with the number of bins rather than shrinking, and at 200 individuals the twenty bin version is off by +0.245 on average with a root mean squared error 3.5 times the likelihood’s.

mc$method <- factor(mc$method, levels = c("likelihood", paste0("bins", bin_grid)),
                    labels = c("likelihood", sprintf("%d bins", bin_grid)))
mc$label  <- factor(ifelse(mc$n_draw == n_small,
                           sprintf("%d individuals", n_small),
                           sprintf("%d individuals", n_ind)))

pa <- ggplot(mc, aes(method, bias, colour = label, group = label)) +
  geom_hline(yintercept = 0, colour = te_line, linewidth = 0.6) +
  geom_line(linewidth = 0.6) + geom_point(size = 2.6) +
  scale_colour_manual(values = c(te_forest, te_rust), name = NULL) +
  labs(x = NULL, y = "bias in b", title = "Binning biases the exponent upwards") +
  theme_datasheet() + theme(legend.position = "none")

pb <- ggplot(mc, aes(method, rmse, colour = label, group = label)) +
  geom_line(linewidth = 0.6) + geom_point(size = 2.6) +
  scale_colour_manual(values = c(te_forest, te_rust), name = NULL) +
  labs(x = NULL, y = "root mean squared error", title = "And it costs precision as well") +
  theme_datasheet() + theme(legend.position = "bottom")

pa / pb + plot_annotation(theme = theme_datasheet())

Two panels of points joined by lines. The upper panel shows bias against method with the likelihood estimator sitting on zero and the three binned regressions rising above it, more steeply at the smaller sample size. The lower panel shows root mean squared error with the same ordering and a wider gap at the smaller sample size. — Figure 2: Bias and error of four estimators of the spectrum exponent, over 400 samples.

The mechanism sits in the top bins. Under the truth, the expected count in the largest classes of a small sample is below one, and a bin cannot hold a fraction of an animal. It either holds none, in which case the logarithm is undefined and the bin is dropped, or it holds at least one, in which case it stands above where the line should pass.

expected_in <- function(bb, n_draw) {
  n_draw * (bb[-1]^(b_true + 1) - bb[-length(bb)]^(b_true + 1)) /
    (m_hi^(b_true + 1) - m_lo^(b_true + 1))
}

set.seed(21)
floor_stat <- replicate(n_rep, {
  x  <- rbounded_pareto(n_small, b_true, m_lo, m_hi)
  bb <- exp(seq(log(m_lo), log(max(x) * 1.0001), length.out = 11))
  ct <- as.numeric(table(cut(x, bb, right = FALSE)))
  ex <- expected_in(bb, n_small)
  thin <- ex < 1
  c(sum(thin), sum(ct[thin] > 0),
    if (any(thin & ct > 0)) mean(ct[thin & ct > 0] / ex[thin & ct > 0]) else NA_real_)
})

thin_bins  <- mean(floor_stat[1, ])
thin_kept  <- mean(floor_stat[2, ])
thin_ratio <- mean(floor_stat[3, ], na.rm = TRUE)

Across the four hundred small samples, an average of 1.7 of the ten bins expect fewer than one animal, 1.2 of them survive into the regression, and those survivors carry on average 5.0 times the count the truth calls for. They sit at the far right of the horizontal range, where a point pulls hardest on a slope, and they can only ever be too high. Adding bins adds more of them, which is why the bias grew rather than shrank.

Where the spectrum starts

The lower limit is the remaining knob, and it is the one that moves the answer furthest. Gear misses small animals: a sieve, a net or a sorter retains large individuals reliably and small ones erratically, so the observed distribution rolls over at the bottom even when the underlying one does not.

n_pool     <- 300000
m_pool_lo  <- 0.05
retain_50  <- 1.5
retain_sd  <- 0.4

set.seed(1)
pool     <- rbounded_pareto(n_pool, b_true, m_pool_lo, m_hi)
p_retain <- 1 / (1 + exp(-(log(pool) - log(retain_50)) / retain_sd))
caught   <- pool[runif(n_pool) < p_retain]

b_naive <- -pareto_alpha(caught, min(caught))

The retained sample holds 11424 individuals. Fitting the closed form from the smallest one in the tray gives -1.267 against a truth of -2.05, an error far larger than anything the binning produced. The estimator is doing exactly what it was asked; the request was wrong.

The usual repair is to fit above a lower limit chosen so that the part of the sample above it looks Pareto, measured by the Kolmogorov-Smirnov distance between the empirical distribution above the limit and the fitted one.

ks_pareto <- function(x, m0) {
  z <- sort(x[x >= m0])
  k <- length(z)
  a <- pareto_alpha(z, m0)
  theo <- 1 - (z / m0)^(-(a - 1))
  list(alpha = a, n_above = k,
       ks = max(abs(c(theo - (seq_len(k) - 1) / k, seq_len(k) / k - theo))))
}

n_grid  <- 60
min_fit <- 100
lo_q    <- 0.002
hi_q    <- 0.92
xmin_grid <- exp(seq(log(unname(quantile(caught, lo_q))),
                     log(unname(quantile(caught, hi_q))), length.out = n_grid))

scan_out <- do.call(rbind, lapply(xmin_grid, function(m0) {
  if (sum(caught >= m0) < min_fit) return(NULL)
  s <- ks_pareto(caught, m0)
  data.frame(xmin = m0, b = -s$alpha, ks = s$ks, n_above = s$n_above)
}))

pick     <- which.min(scan_out$ks)
xmin_ks  <- scan_out$xmin[pick]
b_ks     <- scan_out$b[pick]
n_above  <- scan_out$n_above[pick]
b_span   <- range(scan_out$b)

Minimising that distance picks a lower limit of 3.3 milligrams, keeps 3667 individuals and returns -2.054, which lands on the truth. The more useful output is the rest of the curve: across the whole grid of candidate limits the exponent runs from -1.29 to -2.10. That is most of a unit of exponent, from one tray of animals, produced by moving one knob that the methods section usually does not mention.

ggplot(scan_out, aes(xmin, b)) +
  geom_hline(yintercept = b_true, colour = te_gold, linetype = "dashed",
             linewidth = 0.7) +
  geom_vline(xintercept = xmin_ks, colour = te_rust, linetype = "dashed",
             linewidth = 0.7) +
  geom_line(colour = te_forest, linewidth = 0.9) +
  scale_x_log10() +
  labs(x = "lower mass limit of the fit (mg)", y = "fitted exponent b",
       title = "One tray of animals, most of a unit of exponent",
       subtitle = sprintf("dashed gold: truth %.2f; dashed red: KS limit %.1f mg",
                          b_true, xmin_ks)) +
  theme_datasheet()

A line on a logarithmic horizontal axis showing the fitted exponent rising in magnitude from about minus one point three at the smallest lower limits to about minus two at larger ones, then flattening. A vertical dashed line marks the limit chosen by the Kolmogorov-Smirnov criterion and a horizontal dashed line marks the true exponent, and the two meet on the curve. — Figure 3: Fitted exponent against the lower mass limit, with the limit chosen by minimum KS distance.

Choosing a lower limit this way is the same operation as choosing a threshold for a generalised Pareto fit to extremes, with the same character: the criterion is a rule for picking, not evidence that the pick is right. The curve above is the honest report, and a single number pulled off it is not.

Testing the exponent when the lower limit is estimated

A common next question is whether the exponent differs from a particular value, for instance b = -2, the value at which each doubling class of mass holds the same biomass. When the lower limit is taken as the smallest mass in the sample, which is also the maximum likelihood estimate of the Pareto scale, the likelihood ratio test has an exact null distribution and no chi-square approximation is needed. Write r0 = (alpha0 - 1) / (alpha_hat - 1) for the hypothesised value alpha0 (alpha - 1 is the Pareto shape parameter and the lower limit is its scale, so r0 is the ratio of hypothesised to estimated shape). Then -2 log(LR) = 2 n (r0 - 1 - log(r0)), and it depends on the data only through W = 2 n r0 = 2 (alpha0 - 1) sum(log(x / min(x))). Under the null the values of log(x / m), with m the true lower limit, are exponential with rate alpha0 - 1, and measuring them from the sample minimum instead leaves W distributed exactly as chi-square on 2(n - 1) degrees of freedom (the minimum’s own term is zero, and by the memoryless property of the exponential the other n - 1 values, measured from the minimum, are again independent exponentials with rate alpha0 - 1).

The ratio statistic is U-shaped in W, with its minimum at W = 2n, so an exact p-value adds the two tails of that chi-square beyond the two values of W that give the observed statistic.

lr_pareto <- function(alpha_hat, k, alpha0) {
  r0 <- (alpha0 - 1) / (alpha_hat - 1)
  2 * k * (r0 - 1 - log(r0))
}

# P(LR >= lr) when W is chi-square on df_w: LR = W - 2n - 2n log(W / 2n) is
# U-shaped in W with its minimum at W = 2n, so add the two outer tails of W
lr_tail <- function(lr, n, df_w) {
  h    <- function(v) v - 2 * n - 2 * n * log(v / (2 * n)) - lr
  w_lo <- uniroot(h, c(n * exp(-1 - lr / (2 * n)), 2 * n), tol = 1e-10)$root
  w_hi <- uniroot(h, c(2 * n, 2 * n * (2 + lr / n)), tol = 1e-10)$root
  pchisq(w_lo, df_w) + pchisq(w_hi, df_w, lower.tail = FALSE)
}

pareto_exact_test <- function(x, alpha0) {
  n  <- length(x)
  w  <- 2 * (alpha0 - 1) * sum(log(x / min(x)))
  lr <- lr_pareto(pareto_alpha(x, min(x)), n, alpha0)
  if (lr <= 0) return(c(lr = 0, w = w, p_exact = 1, p_chisq1 = 1))
  c(lr = lr, w = w,
    p_exact  = lr_tail(lr, n, 2 * (n - 1)),
    p_chisq1 = pchisq(lr, 1, lower.tail = FALSE))
}

tray_test <- pareto_exact_test(mass, 2)

n_exact <- c(10, 25, 100)
r_exact <- 4000
alpha0  <- -b_true

set.seed(3341)
exact_check <- do.call(rbind, lapply(n_exact, function(n_draw) {
  sims <- replicate(r_exact, pareto_exact_test(rbounded_pareto(n_draw, b_true, m_lo, Inf), alpha0))
  data.frame(n_draw = n_draw,
             size_chisq1 = mean(sims["p_chisq1", ] < 0.05),
             size_exact  = mean(sims["p_exact", ] < 0.05),
             w_mean_ratio = mean(sims["w", ]) / (2 * (n_draw - 1)))
}))
exact_mcse <- sqrt(0.05 * 0.95 / r_exact)
# the chi-square(1) test's own size is the same two-tail sum, taken at its critical value
exact_chisq1 <- vapply(n_exact, function(n_draw) lr_tail(qchisq(0.95, 1), n_draw, 2 * (n_draw - 1)), 0)

On the tray of 1200 animals from the start of the post, the test of b = -2 gives an exact p-value of 0.095 and a chi-square(1) p-value of 0.094, so the tray does not separate -2 from its true -2.05 at the 0.05 level. At that sample size the two p-values agree. They part company in small samples. At a nominal 0.05 the chi-square(1) reference has an exact size of 0.078 at n = 10, 0.060 at n = 25 and 0.053 at n = 100 (the same two-tail sum, taken at the chi-square(1) critical value). Over 4000 samples drawn with the null true it rejects in 0.075, 0.060 and 0.052 of them, while the exact version gives 0.047, 0.047 and 0.050, within Monte Carlo error of 0.05 at every size (the standard error is 0.0034). The simulated mean of W divided by 2(n - 1) is 1.001 to 1.007, as the chi-square null says it should be.

That settles the case the textbook covers. It is not the case in this post, where the lower limit came from the Kolmogorov-Smirnov scan, and the scan looked at the same animals the test then uses. Clauset, Shalizi and Newman (2009) point out that the standard error (alpha_hat - 1) / sqrt(n) treats the lower limit as known, and they carry the uncertainty in the limit through a nonparametric bootstrap that re-runs the choice of limit on every resample. The same objection applies to the likelihood ratio: nothing in the chi-square(1) reference knows that the limit was picked.

The size of that effect is measurable. The chunk below re-runs the post’s own scan (the same grid of 60 candidate limits, the same minimum of 100 animals above the limit) in a faster form, checks that it picks the same limit on the gear-filtered tray, and then applies the chi-square(1) test of the true exponent after the scan to two kinds of sample. One is a clean power law with no lower tail to cut away, where a limit fixed at the true m_lo is available as a reference. The other is fresh samples drawn through the same gear as in the selectivity section, where the retention curve makes the bottom of the distribution not a power law at all.

scan_ks <- function(x) {
  xs <- sort(x)
  n  <- length(xs)
  lx <- log(xs)
  ends <- log(quantile(xs, c(lo_q, hi_q), names = FALSE))
  g  <- exp(seq(ends[1], ends[2], length.out = n_grid))
  first <- findInterval(g, xs, left.open = TRUE) + 1   # first animal at or above each limit
  keep  <- n - first + 1 >= min_fit
  g <- g[keep]; first <- first[keep]; lg <- log(g)
  best <- Inf; i_best <- NA
  for (i in seq_along(g)) {
    lz <- lx[first[i]:n] - lg[i]
    k  <- length(lz)
    a  <- 1 + k / sum(lz)
    # fitted minus empirical CDF just after each point; the KS distance is its largest
    # excursion on either side of each step
    gap <- 1 - exp(-(a - 1) * lz) - seq_len(k) / k
    d <- max(max(gap) + 1 / k, -min(gap))
    if (d < best) { best <- d; i_best <- i }
  }
  z <- xs[first[i_best]:n]
  c(xmin = g[i_best], alpha = pareto_alpha(z, g[i_best]), n_above = length(z))
}
xmin_fast <- unname(scan_ks(caught)["xmin"])
scan_matches_post <- xmin_fast == xmin_ks

r_gear <- function(n_draw) {
  kept <- numeric(0)
  while (length(kept) < n_draw) {
    pool_r <- rbounded_pareto(40 * n_draw, b_true, m_pool_lo, m_hi)
    p_r    <- 1 / (1 + exp(-(log(pool_r) - log(retain_50)) / retain_sd))
    kept   <- c(kept, pool_r[runif(length(pool_r)) < p_r])
  }
  kept[seq_len(n_draw)]
}

n_scan <- c(300, n_ind)
r_scan <- 1500
crit1  <- qchisq(0.95, 1)
# with the limit known, W is chi-square on 2n, so the size of the chi-square(1) test can be computed exactly
known_size <- vapply(n_scan, function(n_draw) lr_tail(crit1, n_draw, 2 * n_draw), 0)

set.seed(3345)
scan_size <- do.call(rbind, lapply(n_scan, function(n_draw) {
  clean <- replicate(r_scan, {
    x <- rbounded_pareto(n_draw, b_true, m_lo, Inf)
    s <- scan_ks(x)
    c(fixed = lr_pareto(pareto_alpha(x, m_lo), n_draw, alpha0) > crit1,
      scan  = lr_pareto(s["alpha"], s["n_above"], alpha0) > crit1,
      alpha = s[["alpha"]], se = (s[["alpha"]] - 1) / sqrt(s[["n_above"]]))
  })
  gear <- replicate(r_scan, {
    s <- scan_ks(r_gear(n_draw))
    c(scan = lr_pareto(s["alpha"], s["n_above"], alpha0) > crit1,
      xmin = s["xmin"], b = -s["alpha"], se = (s[["alpha"]] - 1) / sqrt(s[["n_above"]]))
  })
  data.frame(n_draw = n_draw,
             clean_fixed = mean(clean[1, ]), clean_scan = mean(clean[2, ]),
             gear_scan = mean(gear[1, ]), gear_xmin = median(gear[2, ]),
             gear_b = mean(gear[3, ]),
             clean_sd_ratio = sd(clean[3, ]) / mean(clean[4, ]),
             gear_sd_ratio  = sd(-gear[3, ]) / mean(gear[4, ]))
}))
scan_size$gear_retain <- 1 / (1 + exp(-(log(scan_size$gear_xmin) - log(retain_50)) / retain_sd))
mcse <- function(p, reps) sqrt(p * (1 - p) / reps)
scan_mcse_max <- max(mcse(unlist(scan_size[, c("clean_fixed", "clean_scan", "gear_scan")]), r_scan))
# distance of the fixed-limit rates from their exact size, in Monte Carlo standard errors
fixed_z <- (scan_size$clean_fixed - known_size) / mcse(known_size, r_scan)

Run on the gear-filtered tray, the faster scan returns exactly the same limit of 3.3 milligrams. With the limit fixed at the truth, the chi-square(1) test rejects the true exponent in 0.060 of 1500 clean samples at n = 300 and 0.052 at n = 1200, against an exact size of 0.050 and 0.050 (computed from the chi-square on 2n degrees of freedom that W follows when the limit is known): 1.8 and 0.4 Monte Carlo standard errors away, so the simulation agrees with the algebra. After the scan, on the same samples, the rate is 0.089 and 0.084. That is a modest excess on a clean power law, and quadrupling the sample barely moves it.

On the gear-filtered samples the same test rejects the true exponent in 0.202 and 0.196 of samples, 4.0 and 3.9 times the nominal rate. The scan’s median limit is 2.0 and 2.8 milligrams, and at the limit itself the gear keeps only 0.67 and 0.82 of the animals it meets, so part of the rolled-over body is inside the fit and the mean exponent is -1.996 and -2.032 against -2.05. The Monte Carlo standard error of every rate in this paragraph and the last is at most 0.0104.

The bootstrap Clauset and colleagues recommend is the obvious repair: resample the animals with replacement, re-run the scan on each resample, and use the spread of the resampled exponents in place of the formula. Whether it restores the size can be measured the same way, at a cost, so the chunk below does it at the smaller sample size only. Two tests come from the same resamples: a Wald test with the bootstrap standard deviation as the standard error, and a percentile interval.

n_boot_draw <- 300
r_boot_clean <- 200
r_boot_gear  <- 500
b_boot <- 100

boot_one <- function(x) {
  s  <- scan_ks(x)
  ab <- replicate(b_boot, scan_ks(x[sample.int(length(x), replace = TRUE)])["alpha"])
  se_boot <- sd(ab)
  q_ab <- quantile(ab, c(0.025, 0.975), names = FALSE)
  c(wald  = abs(s[["alpha"]] - alpha0) / se_boot > qnorm(0.975),
    perc  = alpha0 < q_ab[1] | alpha0 > q_ab[2],
    ratio = se_boot / ((s[["alpha"]] - 1) / sqrt(s[["n_above"]])),
    alpha = s[["alpha"]], xmin = s[["xmin"]], se = se_boot, centre = mean(ab) - s[["alpha"]])
}

set.seed(3346)
boot_clean <- replicate(r_boot_clean, boot_one(rbounded_pareto(n_boot_draw, b_true, m_lo, Inf)))
boot_gear  <- replicate(r_boot_gear, boot_one(r_gear(n_boot_draw)))
boot_size  <- rbind(clean = rowMeans(boot_clean), gear = rowMeans(boot_gear))
boot_reps  <- c(clean = r_boot_clean, gear = r_boot_gear)
boot_se    <- mcse(boot_size[, c("wald", "perc")], boot_reps)
boot_mcse_max <- max(boot_se)
# share of the gear-filtered excess over 0.05 that the bootstrap leaves in place
boot_left <- (boot_size["gear", c("wald", "perc")] - 0.05) / (scan_size$gear_scan[1] - 0.05)
boot_left_se <- sqrt(boot_se["gear", ]^2 + boot_left^2 * mcse(scan_size$gear_scan[1], r_scan)^2) /
  (scan_size$gear_scan[1] - 0.05)
# clean samples: bootstrap rates below the chi-square(1) rate after the scan, in standard errors
clean_gap_z <- (scan_size$clean_scan[1] - boot_size["clean", c("wald", "perc")]) /
  sqrt(boot_se["clean", ]^2 + mcse(scan_size$clean_scan[1], r_scan)^2)

# where the gear-filtered Wald rejections come from
g_alpha <- boot_gear["alpha", ]
g_se    <- boot_gear["se", ]
g_rej   <- boot_gear["wald", ] == 1
gear_diag <- c(
  bias_sd      = (mean(g_alpha) - alpha0) / sd(g_alpha),
  wald_one_sd  = mean(abs(g_alpha - alpha0) / sd(g_alpha) > qnorm(0.975)),
  cor_se       = cor(g_alpha, g_se),
  cor_xmin     = cor(g_alpha, boot_gear["xmin", ]),
  cor_xmin_se  = cor(boot_gear["xmin", ], g_se),
  se_rejected  = mean(g_se[g_rej]) / mean(g_se),
  below_truth  = mean(g_alpha[g_rej] < alpha0),
  centre_sd    = mean(boot_gear["centre", ]) / sd(g_alpha))

With 200 clean and 500 gear-filtered samples and 100 resamples per sample, the answer differs between the two kinds of data. On the clean power law both bootstrap rates are consistent with nominal: the Wald version rejects the true exponent in 0.040 of the samples and the percentile version in 0.065, against 0.089 for the chi-square(1) test after the scan at the same sample size, and neither is more than 0.9 Monte Carlo standard errors from 0.05. The excess on clean data was small to begin with, though, and at this replication only the Wald rate is clearly below the chi-square(1) one (3.1 standard errors of the difference, against 1.3 for the percentile rate). On the gear-filtered samples the bootstrap helps but does not finish the job: 0.124 and 0.084 against 0.202, both still above 0.05 by at least 2.7 Monte Carlo standard errors, so 49 per cent of the excess over 0.05 is still there for the Wald version and 22 per cent for the percentile version (each share with a Monte Carlo standard error of at most 0.10). The Monte Carlo standard error of the four bootstrap rates is at most 0.017.

On average the bootstrap interval is wide enough. Over the 1500 samples of the previous chunk at n = 300, the real spread of the exponent after the scan is 1.19 times the average formula standard error on clean samples and 1.24 times on gear-filtered ones, and the bootstrap standard error is on average 1.28 and 1.30 times the formula’s, so it covers that gap: the extra width is the uncertainty in the limit that the formula leaves out. On the gear-filtered samples the trouble is which samples get the narrow intervals. Across the 500 of them the fitted exponent is on average 0.54 of its own standard deviation shallower than the truth, and a Wald test that used that one standard deviation for every sample would reject in 0.078. With each sample’s own bootstrap standard error it rejects in 0.124, because the error and the interval move together: when the scan settles on a lower limit, more of the rolled-over body enters the fit, the exponent comes out shallower, and the bootstrap spread comes out narrower. Across samples the steepness of the fitted exponent correlates at 0.61 with the chosen limit and at 0.59 with its bootstrap standard error, and the chosen limit at 0.48 with the bootstrap standard error; all 62 rejections are of an exponent shallower than the truth, and the rejected samples have a bootstrap standard error 0.74 times the average. Resampling reproduces the shift rather than correcting it: the net offset of the resampled exponents from the sample’s own estimate, averaged over the samples, is 0.02 standard deviations, against the 0.54 by which the estimates miss the truth, so a bootstrap estimate of the bias would find almost none of it.

size_tab <- rbind(
  data.frame(sample = "clean power law", n_draw = scan_size$n_draw,
             test = "limit fixed at truth", rate = scan_size$clean_fixed, reps = r_scan),
  data.frame(sample = "clean power law", n_draw = scan_size$n_draw,
             test = "KS scan, chi-square(1)", rate = scan_size$clean_scan, reps = r_scan),
  data.frame(sample = "gear-filtered", n_draw = scan_size$n_draw,
             test = "KS scan, chi-square(1)", rate = scan_size$gear_scan, reps = r_scan),
  data.frame(sample = c("clean power law", "gear-filtered"), n_draw = n_boot_draw,
             test = "KS scan, bootstrap SE", rate = boot_size[, "wald"], reps = boot_reps),
  data.frame(sample = c("clean power law", "gear-filtered"), n_draw = n_boot_draw,
             test = "KS scan, bootstrap percentile", rate = boot_size[, "perc"], reps = boot_reps)
)
size_tab$test <- factor(size_tab$test, levels = rev(c("limit fixed at truth", "KS scan, chi-square(1)",
                                                     "KS scan, bootstrap SE", "KS scan, bootstrap percentile")))
size_tab$n_lab <- factor(sprintf("n = %d", size_tab$n_draw),
                         levels = sprintf("n = %d", n_scan))
size_tab$se <- mcse(size_tab$rate, size_tab$reps)

ggplot(size_tab, aes(rate, test, colour = n_lab)) +
  geom_vline(xintercept = 0.05, colour = te_gold, linetype = "dashed", linewidth = 0.7) +
  geom_errorbar(aes(xmin = rate - 2 * se, xmax = rate + 2 * se), width = 0.2,
                orientation = "y", linewidth = 0.6,
                position = position_dodge(width = 0.5)) +
  geom_point(size = 2.6, position = position_dodge(width = 0.5)) +
  facet_wrap(~ sample, ncol = 1, scales = "free_y") +
  scale_colour_manual(values = c(te_forest, te_rust), name = NULL) +
  labs(x = "rejection rate of the true exponent at nominal 0.05", y = NULL,
       title = "Picking the limit from the data inflates the test",
       subtitle = "bars: two Monte Carlo standard errors; dashed gold: nominal 0.05") +
  theme_datasheet() + theme(legend.position = "bottom")

Two stacked panels of points with horizontal error bars of two Monte Carlo standard errors on a rejection-rate axis, with a dashed gold line at the nominal 0.05; green points are samples of 300 animals and red points samples of 1200. In the upper panel, for a clean power law, the test with the limit fixed at the truth sits on or just right of the line, the chi-square(1) test after the KS scan sits near 0.09 at both sample sizes, and the bootstrap standard-error and percentile tests at 300 animals sit near 0.04 and 0.065 with bars that cross the line. In the lower panel, for gear-filtered samples, the chi-square(1) test after the scan sits near 0.2 at both sample sizes, the bootstrap standard-error test near 0.12 with its bar well clear of the line, and the bootstrap percentile test near 0.08 with its bar just clear of it. — Figure 4: Rejection rate of the true exponent at nominal 0.05: lower limit fixed, picked by the KS scan, and with the scan re-run inside a bootstrap.

In practice, then: if the lower limit is the sample minimum, the exact test above is the answer and needs nothing else. If the limit came from a scan, a chi-square(1) p-value at the chosen limit rejects a true exponent too often, a little on a clean power law and a lot when the gear has left a rolled-over body behind. The resample-and-rescan bootstrap is the better default: on a clean power law both of its rates are within Monte Carlo error of 0.05, but on gear-filtered data it removes only part of the excess, because the scan shifts the exponent, each resample reproduces the shift, and the most shifted samples get the narrowest intervals. The sensitivity curve still belongs in the report, because a test at one chosen limit says nothing about the others. These rates belong to one exponent, one scan grid and one retention curve, and a coarser grid, a larger minimum tail or a steeper retention curve will move them. The bootstrap was measured at one sample size and with 100 resamples, against the 1000 or so Clauset and colleagues suggest; a noisier bootstrap standard error tends to make the Wald version reject more often, not less.

What to report

Say which spectrum the exponent belongs to. Counts per class, counts per unit mass and biomass per class differ by whole units of slope, biomass per unit mass shares the slope of counts per class, and a value quoted without that label cannot be compared with anything.

Fit by likelihood, not by regression on binned counts. The closed form is one line, and so is its standard error when the lower limit is fixed in advance; when the limit comes from a scan, take the standard error from a bootstrap that re-runs the scan, as in the section above. The four hundred replicate comparison above shows the closed form beating every binned alternative on bias and on error at both sample sizes.

Give the lower mass limit, the number of individuals above it, and what the exponent does as the limit moves. The limit is where the real disagreement between analyses lives, and a sensitivity curve costs one plot. A test of the exponent at a scanned limit needs the bootstrap, not the chi-square(1) table.

If a binned figure is wanted, and it usually is, plot the binned spectrum and fit the likelihood. The two are not in competition once the line on the figure is drawn from the likelihood estimate rather than from a regression through the points.

Honest limits

The data here are generated from a power law, so every method is being scored on its ability to recover a truth that exists. Real body mass distributions are not exactly power laws, and when the shape is wrong the likelihood estimator is not saved by being efficient: it returns the exponent of the best fitting power law, which may describe nothing. Testing whether a power law is the right family at all is a separate exercise, and the site has covered how weak that kind of comparison can be when two families are close.

The Kolmogorov-Smirnov criterion for the lower limit is a heuristic with known failure modes. It tends to pick a limit that discards more data than necessary when the departure from a power law is mild, and it has no defence at all against a distribution that is curved throughout, where it will find a local stretch that looks straight. The flat section of the curve above is doing more of the work than the minimum is.

The selectivity example uses a smooth retention curve that is still short of one at the limits the scan chooses: at the median limits in the test above the gear keeps only 0.67 and 0.82 of the animals it meets, which is why part of the rolled-over body enters the fit. Gear whose retention is still climbing across the whole fitted range is worse again: it biases the exponent in a way no choice of lower limit can remove, because there is no mass above which the sample is a fair draw.

The upper end is left alone throughout. A bounded power law was simulated and an unbounded one was fitted, which costs little here because the largest observed mass sits far below the bound. Samples whose largest individuals press against a real physiological or gear-imposed ceiling need the bounded likelihood, and the difference is then not small.

Finally, 400 replicates give each bias a Monte Carlo standard error of at most 0.0063 on the scale of the exponent. That is below every gap discussed above, including the 0.023 between six and ten bins at the larger sample size, so the ordering is solid. The fourth decimal place of any one of them is not.

References

Edwards AM, Robinson JPW, Plank MJ, Baum JK, Blanchard JL 2017 Methods in Ecology and Evolution 8(1):57-67 (10.1111/2041-210X.12641)

White EP, Enquist BJ, Green JL 2008 Ecology 89(4):905-912 (10.1890/07-1288.1)

Clauset A, Shalizi CR, Newman MEJ 2009 SIAM Review 51(4):661-703 (10.1137/070710111)

Sheldon RW, Prakash A, Sutcliffe WH 1972 Limnology and Oceanography 17(3):327-340 (10.4319/lo.1972.17.3.0327)

--- title: "Fitting a body size spectrum" description: "Log bins shift the fitted size spectrum exponent by exactly one, and the likelihood never sees a bin. Fitting individual body mass data in R, no histogram." date: "2026-08-08 14:00" date-modified: "2026-10-01" categories: [R, body size, power laws, statistics, ecology tutorial] image: thumbnail.png image-alt: "A log-log plot on a warm off-white background with two series of points from the same binned body mass sample, counts per class in red and counts per unit mass in green, and two fitted lines of visibly different steepness." --- *Corrected 1 October 2026: the post said the four binned quantities give four different slopes; counts per class and biomass per unit mass share one, so there are three, and What to report now says so.* *Corrected 27 September 2026: Honest limits said the example retention curve reaches one well below the chosen limit; it now says the curve is still short of one at the limits the scan chooses. The generalised Pareto remark now says that `alpha - 1`, not `alpha`, is the reciprocal of the shape parameter.* *Updated 25 September 2026: the post fitted the exponent but did not say how to test it. A new section, Testing the exponent when the lower limit is estimated, gives the exact test when the limit is the sample minimum and measures how often the usual test rejects a true exponent once the limit comes from a scan, with and without a bootstrap that re-runs the scan. What to report now takes the standard error from that bootstrap when the limit is scanned.* Weigh every individual in a benthic sample and the masses span four or five orders of magnitude, with thousands of tiny animals and a handful of large ones. The standard summary of that spread is the size spectrum: a claim that the number of individuals at mass `m` falls off as a power of `m`, and a single exponent that says how fast. The exponent gets compared across sites, across years and against fishing pressure, so it matters that two people analysing the same sample arrive at the same number. They frequently do not, and the reason is rarely the ecology. It is that the exponent depends on whether the counts were divided by the width of the bin they fell in, on how many bins there were, and on where the smallest mass in the fit was set. None of those is a property of the animals. This post separates the three. The bin width turns out to shift the fitted slope by exactly one, the number of bins turns out to bias it in a direction that can be measured, and the maximum likelihood estimator turns out not to need bins at all. ## A sample of individual masses Individual dry masses in milligrams, drawn from a bounded power law so that the answer is known. The density is proportional to `m^b` between a lower and an upper mass, and `b` is the quantity every method below is trying to recover. ```{r setup} #| message: false #| warning: false library(ggplot2) library(patchwork) 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)) } rbounded_pareto <- function(n, b, lo, hi) { u <- runif(n) ((hi^(b + 1) - lo^(b + 1)) * u + lo^(b + 1))^(1 / (b + 1)) } b_true <- -2.05 m_lo <- 0.5 m_hi <- 8000 n_ind <- 1200 set.seed(11) mass <- rbounded_pareto(n_ind, b_true, m_lo, m_hi) mass_max <- max(mass) ``` The sample holds `r sprintf("%d", n_ind)` individuals between `r sprintf("%.2f", min(mass))` and `r sprintf("%.0f", mass_max)` milligrams. Nothing in the analysis below uses `b_true` again until the end. ## The bin width is not a plotting choice Bins that double in width are the convention, because equal bins on the arithmetic scale put almost every animal in the first one. The bookkeeping question is what goes on the vertical axis: the raw count in each doubling class, or that count divided by the width of the class. ```{r bins} brk <- 2^(floor(log2(m_lo)):ceiling(log2(mass_max))) cls <- cut(mass, brk, right = FALSE) bm_in <- as.numeric(tapply(mass, cls, sum)) bm_in[is.na(bm_in)] <- 0 sp <- data.frame( lower = brk[-length(brk)], upper = brk[-1], n_in = as.numeric(table(cls)), biomass = bm_in ) sp$width <- sp$upper - sp$lower sp$mid <- sqrt(sp$lower * sp$upper) sp$dens <- sp$n_in / sp$width sp <- sp[sp$n_in > 0, ] fit_count <- lm(log(n_in) ~ log(mid), data = sp) fit_dens <- lm(log(dens) ~ log(mid), data = sp) slope_count <- unname(coef(fit_count)[2]) slope_dens <- unname(coef(fit_dens)[2]) slope_shift <- slope_count - slope_dens n_bin_used <- nrow(sp) ``` The two lines are fitted to the same `r sprintf("%d", n_bin_used)` bins holding the same `r sprintf("%d", n_ind)` animals. Counting gives a slope of `r sprintf("%.3f", slope_count)`. Dividing by bin width gives `r sprintf("%.3f", slope_dens)`. The gap between them is `r sprintf("%.4f", slope_shift)`, and it is not approximately one, it is one. The reason is that a doubling class starting at mass `a` has width `a`, so its width is proportional to its own position on the axis. Dividing the count by `a` subtracts `log(a)` from the vertical coordinate while the horizontal coordinate stays at `log(a)`, which tilts the whole line by exactly one unit of slope. Any binning scheme with widths proportional to mass does this; the factor is a property of the bins, not of the sample. ```{r fig-normalisation} #| fig-cap: "One sample of individual masses, binned once, plotted two ways." #| fig-alt: "Two log-log scatter panels side by side. The left panel plots raw counts per doubling class against class midpoint with a fitted line of slope near minus one. The right panel plots counts divided by class width against the same midpoints, with a visibly steeper fitted line of slope near minus two." pl <- ggplot(sp, aes(mid, n_in)) + geom_point(colour = te_forest, size = 2.4) + geom_smooth(method = "lm", formula = y ~ x, se = FALSE, colour = te_rust, linewidth = 0.7) + scale_x_log10() + scale_y_log10() + labs(x = "body mass (mg)", y = "individuals in class", title = "Counts per class", subtitle = sprintf("fitted slope %.2f", slope_count)) + theme_datasheet() pr <- ggplot(sp, aes(mid, dens)) + geom_point(colour = te_forest, size = 2.4) + geom_smooth(method = "lm", formula = y ~ x, se = FALSE, colour = te_rust, linewidth = 0.7) + scale_x_log10() + scale_y_log10(labels = scales::label_comma(drop0trailing = TRUE)) + labs(x = "body mass (mg)", y = "individuals per mg", title = "Counts per unit mass", subtitle = sprintf("fitted slope %.2f", slope_dens)) + theme_datasheet() pl + pr + plot_annotation(theme = theme_datasheet()) ``` This is the same bookkeeping that trips up dispersal kernels, where the density of distances carries a ring area that the density of positions does not, and the same failure to divide by the width of the measure produces a scale parameter that is wrong by a fixed factor. The bookkeeping continues into biomass. Total mass in a doubling class rises one further power of `m` above the count, and dividing that by class width takes one back off, so four quantities from one sample give only three different slopes in expectation: biomass per unit mass shares the slope of counts per class. ```{r spectra} slope_bm <- unname(coef(lm(log(biomass) ~ log(mid), data = sp))[2]) slope_bm_dens <- unname(coef(lm(log(biomass / width) ~ log(mid), data = sp))[2]) ``` Counts give `r sprintf("%.2f", slope_count)`, counts per unit mass give `r sprintf("%.2f", slope_dens)`, biomass gives `r sprintf("%.2f", slope_bm)` and biomass per unit mass gives `r sprintf("%.2f", slope_bm_dens)`. A published slope near `r sprintf("%.0f", -1)` and a published slope near `r sprintf("%.0f", -2)` can describe an identical sample. The number is uninterpretable without the sentence that says which of the four it is, and that sentence is missing from a lot of papers. ## The likelihood never sees a bin If the masses really do follow a power law above some lower limit, the exponent has a closed form. Writing the density as proportional to `m^b` with `b` below minus one, and setting `alpha = -b`, the maximum likelihood estimate from `n` observed masses above a lower limit `m0` is a single expression, and so is its standard error. ```{r mle} pareto_alpha <- function(x, m0) 1 + length(x) / sum(log(x / m0)) alpha_hat <- pareto_alpha(mass, m_lo) b_hat <- -alpha_hat alpha_se <- (alpha_hat - 1) / sqrt(length(mass)) nll_bounded <- function(b, x, lo, hi) { -sum(log((b + 1) * x^b / (hi^(b + 1) - lo^(b + 1)))) } b_bounded <- optimize(nll_bounded, c(-4, -1.01), x = mass, lo = m_lo, hi = m_hi)$minimum ``` The estimate is `r sprintf("%.3f", b_hat)` with a standard error of `r sprintf("%.3f", alpha_se)`, against a truth of `r sprintf("%.2f", b_true)`. Allowing for the upper bound leaves it at `r sprintf("%.3f", b_bounded)`: the largest animal in the tray sits far below the bound, so that refinement has nothing to correct here. Two things are worth noticing about that expression. It contains no bins, no midpoints and no regression, so the choice of doubling class or decade class or any other class cannot influence it. And it is the same estimator that sits behind the generalised Pareto fit used for extremes: a Pareto tail is a generalised Pareto with the scale tied to the threshold, and `alpha - 1` is the reciprocal of the generalised Pareto shape parameter. ## What the bins cost Comparing estimators needs repetition rather than one sample, so here are four hundred fresh samples at two sample sizes, each fitted by the closed form and by a log-binned regression at three bin counts. ```{r montecarlo} n_rep <- 400 binned_b <- function(x, m0, n_bin) { bb <- exp(seq(log(m0), log(max(x) * 1.0001), length.out = n_bin + 1)) ct <- as.numeric(table(cut(x, bb, right = FALSE))) wd <- diff(bb) md <- sqrt(bb[-length(bb)] * bb[-1]) keep <- ct > 0 unname(coef(lm(log(ct[keep] / wd[keep]) ~ log(md[keep])))[2]) } n_small <- 200 bin_grid <- c(6, 10, 20) set.seed(202) sweep_one <- function(n_draw) { out <- t(replicate(n_rep, { x <- rbounded_pareto(n_draw, b_true, m_lo, m_hi) c(-pareto_alpha(x, m_lo), vapply(bin_grid, function(k) binned_b(x, m_lo, k), 0)) })) colnames(out) <- c("likelihood", paste0("bins", bin_grid)) data.frame( method = colnames(out), n_draw = n_draw, bias = colMeans(out) - b_true, rmse = sqrt(colMeans((out - b_true)^2)), mc_se = apply(out, 2, sd) / sqrt(n_rep), row.names = NULL ) } mc <- rbind(sweep_one(n_small), sweep_one(n_ind)) mle_rmse_small <- mc$rmse[mc$method == "likelihood" & mc$n_draw == n_small] b20_rmse_small <- mc$rmse[mc$method == "bins20" & mc$n_draw == n_small] b20_bias_small <- mc$bias[mc$method == "bins20" & mc$n_draw == n_small] b06_bias_small <- mc$bias[mc$method == "bins6" & mc$n_draw == n_small] mle_bias_big <- mc$bias[mc$method == "likelihood" & mc$n_draw == n_ind] b20_bias_big <- mc$bias[mc$method == "bins20" & mc$n_draw == n_ind] rmse_ratio <- b20_rmse_small / mle_rmse_small mc_se_max <- max(mc$mc_se) b06_bias_big <- mc$bias[mc$method == "bins6" & mc$n_draw == n_ind] b10_bias_big <- mc$bias[mc$method == "bins10" & mc$n_draw == n_ind] bin_gap_big <- b10_bias_big - b06_bias_big ``` The closed form is essentially unbiased at both sample sizes, with a bias of `r sprintf("%+.3f", mle_bias_big)` at `r sprintf("%d", n_ind)` individuals. Every binned fit is biased towards a shallower spectrum, the bias grows with the number of bins rather than shrinking, and at `r sprintf("%d", n_small)` individuals the twenty bin version is off by `r sprintf("%+.3f", b20_bias_small)` on average with a root mean squared error `r sprintf("%.1f", rmse_ratio)` times the likelihood's. ```{r fig-cost} #| fig-cap: "Bias and error of four estimators of the spectrum exponent, over 400 samples." #| fig-alt: "Two panels of points joined by lines. The upper panel shows bias against method with the likelihood estimator sitting on zero and the three binned regressions rising above it, more steeply at the smaller sample size. The lower panel shows root mean squared error with the same ordering and a wider gap at the smaller sample size." mc$method <- factor(mc$method, levels = c("likelihood", paste0("bins", bin_grid)), labels = c("likelihood", sprintf("%d bins", bin_grid))) mc$label <- factor(ifelse(mc$n_draw == n_small, sprintf("%d individuals", n_small), sprintf("%d individuals", n_ind))) pa <- ggplot(mc, aes(method, bias, colour = label, group = label)) + geom_hline(yintercept = 0, colour = te_line, linewidth = 0.6) + geom_line(linewidth = 0.6) + geom_point(size = 2.6) + scale_colour_manual(values = c(te_forest, te_rust), name = NULL) + labs(x = NULL, y = "bias in b", title = "Binning biases the exponent upwards") + theme_datasheet() + theme(legend.position = "none") pb <- ggplot(mc, aes(method, rmse, colour = label, group = label)) + geom_line(linewidth = 0.6) + geom_point(size = 2.6) + scale_colour_manual(values = c(te_forest, te_rust), name = NULL) + labs(x = NULL, y = "root mean squared error", title = "And it costs precision as well") + theme_datasheet() + theme(legend.position = "bottom") pa / pb + plot_annotation(theme = theme_datasheet()) ``` The mechanism sits in the top bins. Under the truth, the expected count in the largest classes of a small sample is below one, and a bin cannot hold a fraction of an animal. It either holds none, in which case the logarithm is undefined and the bin is dropped, or it holds at least one, in which case it stands above where the line should pass. ```{r countfloor} expected_in <- function(bb, n_draw) { n_draw * (bb[-1]^(b_true + 1) - bb[-length(bb)]^(b_true + 1)) / (m_hi^(b_true + 1) - m_lo^(b_true + 1)) } set.seed(21) floor_stat <- replicate(n_rep, { x <- rbounded_pareto(n_small, b_true, m_lo, m_hi) bb <- exp(seq(log(m_lo), log(max(x) * 1.0001), length.out = 11)) ct <- as.numeric(table(cut(x, bb, right = FALSE))) ex <- expected_in(bb, n_small) thin <- ex < 1 c(sum(thin), sum(ct[thin] > 0), if (any(thin & ct > 0)) mean(ct[thin & ct > 0] / ex[thin & ct > 0]) else NA_real_) }) thin_bins <- mean(floor_stat[1, ]) thin_kept <- mean(floor_stat[2, ]) thin_ratio <- mean(floor_stat[3, ], na.rm = TRUE) ``` Across the four hundred small samples, an average of `r sprintf("%.1f", thin_bins)` of the ten bins expect fewer than one animal, `r sprintf("%.1f", thin_kept)` of them survive into the regression, and those survivors carry on average `r sprintf("%.1f", thin_ratio)` times the count the truth calls for. They sit at the far right of the horizontal range, where a point pulls hardest on a slope, and they can only ever be too high. Adding bins adds more of them, which is why the bias grew rather than shrank. ## Where the spectrum starts The lower limit is the remaining knob, and it is the one that moves the answer furthest. Gear misses small animals: a sieve, a net or a sorter retains large individuals reliably and small ones erratically, so the observed distribution rolls over at the bottom even when the underlying one does not. ```{r selectivity} n_pool <- 300000 m_pool_lo <- 0.05 retain_50 <- 1.5 retain_sd <- 0.4 set.seed(1) pool <- rbounded_pareto(n_pool, b_true, m_pool_lo, m_hi) p_retain <- 1 / (1 + exp(-(log(pool) - log(retain_50)) / retain_sd)) caught <- pool[runif(n_pool) < p_retain] b_naive <- -pareto_alpha(caught, min(caught)) ``` The retained sample holds `r sprintf("%d", length(caught))` individuals. Fitting the closed form from the smallest one in the tray gives `r sprintf("%.3f", b_naive)` against a truth of `r sprintf("%.2f", b_true)`, an error far larger than anything the binning produced. The estimator is doing exactly what it was asked; the request was wrong. The usual repair is to fit above a lower limit chosen so that the part of the sample above it looks Pareto, measured by the Kolmogorov-Smirnov distance between the empirical distribution above the limit and the fitted one. ```{r xminscan} ks_pareto <- function(x, m0) { z <- sort(x[x >= m0]) k <- length(z) a <- pareto_alpha(z, m0) theo <- 1 - (z / m0)^(-(a - 1)) list(alpha = a, n_above = k, ks = max(abs(c(theo - (seq_len(k) - 1) / k, seq_len(k) / k - theo)))) } n_grid <- 60 min_fit <- 100 lo_q <- 0.002 hi_q <- 0.92 xmin_grid <- exp(seq(log(unname(quantile(caught, lo_q))), log(unname(quantile(caught, hi_q))), length.out = n_grid)) scan_out <- do.call(rbind, lapply(xmin_grid, function(m0) { if (sum(caught >= m0) < min_fit) return(NULL) s <- ks_pareto(caught, m0) data.frame(xmin = m0, b = -s$alpha, ks = s$ks, n_above = s$n_above) })) pick <- which.min(scan_out$ks) xmin_ks <- scan_out$xmin[pick] b_ks <- scan_out$b[pick] n_above <- scan_out$n_above[pick] b_span <- range(scan_out$b) ``` Minimising that distance picks a lower limit of `r sprintf("%.1f", xmin_ks)` milligrams, keeps `r sprintf("%d", n_above)` individuals and returns `r sprintf("%.3f", b_ks)`, which lands on the truth. The more useful output is the rest of the curve: across the whole grid of candidate limits the exponent runs from `r sprintf("%.2f", b_span[2])` to `r sprintf("%.2f", b_span[1])`. That is most of a unit of exponent, from one tray of animals, produced by moving one knob that the methods section usually does not mention. ```{r fig-xmin} #| fig-cap: "Fitted exponent against the lower mass limit, with the limit chosen by minimum KS distance." #| fig-alt: "A line on a logarithmic horizontal axis showing the fitted exponent rising in magnitude from about minus one point three at the smallest lower limits to about minus two at larger ones, then flattening. A vertical dashed line marks the limit chosen by the Kolmogorov-Smirnov criterion and a horizontal dashed line marks the true exponent, and the two meet on the curve." ggplot(scan_out, aes(xmin, b)) + geom_hline(yintercept = b_true, colour = te_gold, linetype = "dashed", linewidth = 0.7) + geom_vline(xintercept = xmin_ks, colour = te_rust, linetype = "dashed", linewidth = 0.7) + geom_line(colour = te_forest, linewidth = 0.9) + scale_x_log10() + labs(x = "lower mass limit of the fit (mg)", y = "fitted exponent b", title = "One tray of animals, most of a unit of exponent", subtitle = sprintf("dashed gold: truth %.2f; dashed red: KS limit %.1f mg", b_true, xmin_ks)) + theme_datasheet() ``` Choosing a lower limit this way is the same operation as choosing a threshold for a generalised Pareto fit to extremes, with the same character: the criterion is a rule for picking, not evidence that the pick is right. The curve above is the honest report, and a single number pulled off it is not. ## Testing the exponent when the lower limit is estimated A common next question is whether the exponent differs from a particular value, for instance `b = -2`, the value at which each doubling class of mass holds the same biomass. When the lower limit is taken as the smallest mass in the sample, which is also the maximum likelihood estimate of the Pareto scale, the likelihood ratio test has an exact null distribution and no chi-square approximation is needed. Write `r0 = (alpha0 - 1) / (alpha_hat - 1)` for the hypothesised value `alpha0` (`alpha - 1` is the Pareto shape parameter and the lower limit is its scale, so `r0` is the ratio of hypothesised to estimated shape). Then `-2 log(LR) = 2 n (r0 - 1 - log(r0))`, and it depends on the data only through `W = 2 n r0 = 2 (alpha0 - 1) sum(log(x / min(x)))`. Under the null the values of `log(x / m)`, with `m` the true lower limit, are exponential with rate `alpha0 - 1`, and measuring them from the sample minimum instead leaves `W` distributed exactly as chi-square on `2(n - 1)` degrees of freedom (the minimum's own term is zero, and by the memoryless property of the exponential the other `n - 1` values, measured from the minimum, are again independent exponentials with rate `alpha0 - 1`). The ratio statistic is U-shaped in `W`, with its minimum at `W = 2n`, so an exact p-value adds the two tails of that chi-square beyond the two values of `W` that give the observed statistic. ```{r lr-exact} lr_pareto <- function(alpha_hat, k, alpha0) { r0 <- (alpha0 - 1) / (alpha_hat - 1) 2 * k * (r0 - 1 - log(r0)) } # P(LR >= lr) when W is chi-square on df_w: LR = W - 2n - 2n log(W / 2n) is # U-shaped in W with its minimum at W = 2n, so add the two outer tails of W lr_tail <- function(lr, n, df_w) { h <- function(v) v - 2 * n - 2 * n * log(v / (2 * n)) - lr w_lo <- uniroot(h, c(n * exp(-1 - lr / (2 * n)), 2 * n), tol = 1e-10)$root w_hi <- uniroot(h, c(2 * n, 2 * n * (2 + lr / n)), tol = 1e-10)$root pchisq(w_lo, df_w) + pchisq(w_hi, df_w, lower.tail = FALSE) } pareto_exact_test <- function(x, alpha0) { n <- length(x) w <- 2 * (alpha0 - 1) * sum(log(x / min(x))) lr <- lr_pareto(pareto_alpha(x, min(x)), n, alpha0) if (lr <= 0) return(c(lr = 0, w = w, p_exact = 1, p_chisq1 = 1)) c(lr = lr, w = w, p_exact = lr_tail(lr, n, 2 * (n - 1)), p_chisq1 = pchisq(lr, 1, lower.tail = FALSE)) } tray_test <- pareto_exact_test(mass, 2) n_exact <- c(10, 25, 100) r_exact <- 4000 alpha0 <- -b_true set.seed(3341) exact_check <- do.call(rbind, lapply(n_exact, function(n_draw) { sims <- replicate(r_exact, pareto_exact_test(rbounded_pareto(n_draw, b_true, m_lo, Inf), alpha0)) data.frame(n_draw = n_draw, size_chisq1 = mean(sims["p_chisq1", ] < 0.05), size_exact = mean(sims["p_exact", ] < 0.05), w_mean_ratio = mean(sims["w", ]) / (2 * (n_draw - 1))) })) exact_mcse <- sqrt(0.05 * 0.95 / r_exact) # the chi-square(1) test's own size is the same two-tail sum, taken at its critical value exact_chisq1 <- vapply(n_exact, function(n_draw) lr_tail(qchisq(0.95, 1), n_draw, 2 * (n_draw - 1)), 0) ``` On the tray of `r sprintf("%d", n_ind)` animals from the start of the post, the test of `b = -2` gives an exact p-value of `r sprintf("%.3f", tray_test["p_exact"])` and a chi-square(1) p-value of `r sprintf("%.3f", tray_test["p_chisq1"])`, so the tray does not separate `-2` from its true `r sprintf("%.2f", b_true)` at the 0.05 level. At that sample size the two p-values agree. They part company in small samples. At a nominal 0.05 the chi-square(1) reference has an exact size of `r sprintf("%.3f", exact_chisq1[1])` at n = `r sprintf("%d", n_exact[1])`, `r sprintf("%.3f", exact_chisq1[2])` at n = `r sprintf("%d", n_exact[2])` and `r sprintf("%.3f", exact_chisq1[3])` at n = `r sprintf("%d", n_exact[3])` (the same two-tail sum, taken at the chi-square(1) critical value). Over `r sprintf("%d", r_exact)` samples drawn with the null true it rejects in `r sprintf("%.3f", exact_check$size_chisq1[1])`, `r sprintf("%.3f", exact_check$size_chisq1[2])` and `r sprintf("%.3f", exact_check$size_chisq1[3])` of them, while the exact version gives `r sprintf("%.3f", exact_check$size_exact[1])`, `r sprintf("%.3f", exact_check$size_exact[2])` and `r sprintf("%.3f", exact_check$size_exact[3])`, within Monte Carlo error of 0.05 at every size (the standard error is `r sprintf("%.4f", exact_mcse)`). The simulated mean of `W` divided by `2(n - 1)` is `r sprintf("%.3f", min(exact_check$w_mean_ratio))` to `r sprintf("%.3f", max(exact_check$w_mean_ratio))`, as the chi-square null says it should be. That settles the case the textbook covers. It is not the case in this post, where the lower limit came from the Kolmogorov-Smirnov scan, and the scan looked at the same animals the test then uses. Clauset, Shalizi and Newman (2009) point out that the standard error `(alpha_hat - 1) / sqrt(n)` treats the lower limit as known, and they carry the uncertainty in the limit through a nonparametric bootstrap that re-runs the choice of limit on every resample. The same objection applies to the likelihood ratio: nothing in the chi-square(1) reference knows that the limit was picked. The size of that effect is measurable. The chunk below re-runs the post's own scan (the same grid of `r sprintf("%d", n_grid)` candidate limits, the same minimum of `r sprintf("%d", min_fit)` animals above the limit) in a faster form, checks that it picks the same limit on the gear-filtered tray, and then applies the chi-square(1) test of the true exponent after the scan to two kinds of sample. One is a clean power law with no lower tail to cut away, where a limit fixed at the true `m_lo` is available as a reference. The other is fresh samples drawn through the same gear as in the selectivity section, where the retention curve makes the bottom of the distribution not a power law at all. ```{r lr-scan} scan_ks <- function(x) { xs <- sort(x) n <- length(xs) lx <- log(xs) ends <- log(quantile(xs, c(lo_q, hi_q), names = FALSE)) g <- exp(seq(ends[1], ends[2], length.out = n_grid)) first <- findInterval(g, xs, left.open = TRUE) + 1 # first animal at or above each limit keep <- n - first + 1 >= min_fit g <- g[keep]; first <- first[keep]; lg <- log(g) best <- Inf; i_best <- NA for (i in seq_along(g)) { lz <- lx[first[i]:n] - lg[i] k <- length(lz) a <- 1 + k / sum(lz) # fitted minus empirical CDF just after each point; the KS distance is its largest # excursion on either side of each step gap <- 1 - exp(-(a - 1) * lz) - seq_len(k) / k d <- max(max(gap) + 1 / k, -min(gap)) if (d < best) { best <- d; i_best <- i } } z <- xs[first[i_best]:n] c(xmin = g[i_best], alpha = pareto_alpha(z, g[i_best]), n_above = length(z)) } xmin_fast <- unname(scan_ks(caught)["xmin"]) scan_matches_post <- xmin_fast == xmin_ks r_gear <- function(n_draw) { kept <- numeric(0) while (length(kept) < n_draw) { pool_r <- rbounded_pareto(40 * n_draw, b_true, m_pool_lo, m_hi) p_r <- 1 / (1 + exp(-(log(pool_r) - log(retain_50)) / retain_sd)) kept <- c(kept, pool_r[runif(length(pool_r)) < p_r]) } kept[seq_len(n_draw)] } n_scan <- c(300, n_ind) r_scan <- 1500 crit1 <- qchisq(0.95, 1) # with the limit known, W is chi-square on 2n, so the size of the chi-square(1) test can be computed exactly known_size <- vapply(n_scan, function(n_draw) lr_tail(crit1, n_draw, 2 * n_draw), 0) set.seed(3345) scan_size <- do.call(rbind, lapply(n_scan, function(n_draw) { clean <- replicate(r_scan, { x <- rbounded_pareto(n_draw, b_true, m_lo, Inf) s <- scan_ks(x) c(fixed = lr_pareto(pareto_alpha(x, m_lo), n_draw, alpha0) > crit1, scan = lr_pareto(s["alpha"], s["n_above"], alpha0) > crit1, alpha = s[["alpha"]], se = (s[["alpha"]] - 1) / sqrt(s[["n_above"]])) }) gear <- replicate(r_scan, { s <- scan_ks(r_gear(n_draw)) c(scan = lr_pareto(s["alpha"], s["n_above"], alpha0) > crit1, xmin = s["xmin"], b = -s["alpha"], se = (s[["alpha"]] - 1) / sqrt(s[["n_above"]])) }) data.frame(n_draw = n_draw, clean_fixed = mean(clean[1, ]), clean_scan = mean(clean[2, ]), gear_scan = mean(gear[1, ]), gear_xmin = median(gear[2, ]), gear_b = mean(gear[3, ]), clean_sd_ratio = sd(clean[3, ]) / mean(clean[4, ]), gear_sd_ratio = sd(-gear[3, ]) / mean(gear[4, ])) })) scan_size$gear_retain <- 1 / (1 + exp(-(log(scan_size$gear_xmin) - log(retain_50)) / retain_sd)) mcse <- function(p, reps) sqrt(p * (1 - p) / reps) scan_mcse_max <- max(mcse(unlist(scan_size[, c("clean_fixed", "clean_scan", "gear_scan")]), r_scan)) # distance of the fixed-limit rates from their exact size, in Monte Carlo standard errors fixed_z <- (scan_size$clean_fixed - known_size) / mcse(known_size, r_scan) ``` Run on the gear-filtered tray, the faster scan returns `r if (scan_matches_post) "exactly the same" else "a different"` limit of `r sprintf("%.1f", xmin_fast)` milligrams. With the limit fixed at the truth, the chi-square(1) test rejects the true exponent in `r sprintf("%.3f", scan_size$clean_fixed[1])` of `r sprintf("%d", r_scan)` clean samples at n = `r sprintf("%d", n_scan[1])` and `r sprintf("%.3f", scan_size$clean_fixed[2])` at n = `r sprintf("%d", n_scan[2])`, against an exact size of `r sprintf("%.3f", known_size[1])` and `r sprintf("%.3f", known_size[2])` (computed from the chi-square on `2n` degrees of freedom that `W` follows when the limit is known): `r sprintf("%.1f", abs(fixed_z[1]))` and `r sprintf("%.1f", abs(fixed_z[2]))` Monte Carlo standard errors away, so the simulation agrees with the algebra. After the scan, on the same samples, the rate is `r sprintf("%.3f", scan_size$clean_scan[1])` and `r sprintf("%.3f", scan_size$clean_scan[2])`. That is a modest excess on a clean power law, and quadrupling the sample barely moves it. On the gear-filtered samples the same test rejects the true exponent in `r sprintf("%.3f", scan_size$gear_scan[1])` and `r sprintf("%.3f", scan_size$gear_scan[2])` of samples, `r sprintf("%.1f", scan_size$gear_scan[1] / 0.05)` and `r sprintf("%.1f", scan_size$gear_scan[2] / 0.05)` times the nominal rate. The scan's median limit is `r sprintf("%.1f", scan_size$gear_xmin[1])` and `r sprintf("%.1f", scan_size$gear_xmin[2])` milligrams, and at the limit itself the gear keeps only `r sprintf("%.2f", scan_size$gear_retain[1])` and `r sprintf("%.2f", scan_size$gear_retain[2])` of the animals it meets, so part of the rolled-over body is inside the fit and the mean exponent is `r sprintf("%.3f", scan_size$gear_b[1])` and `r sprintf("%.3f", scan_size$gear_b[2])` against `r sprintf("%.2f", b_true)`. The Monte Carlo standard error of every rate in this paragraph and the last is at most `r sprintf("%.4f", scan_mcse_max)`. The bootstrap Clauset and colleagues recommend is the obvious repair: resample the animals with replacement, re-run the scan on each resample, and use the spread of the resampled exponents in place of the formula. Whether it restores the size can be measured the same way, at a cost, so the chunk below does it at the smaller sample size only. Two tests come from the same resamples: a Wald test with the bootstrap standard deviation as the standard error, and a percentile interval. ```{r lr-boot} n_boot_draw <- 300 r_boot_clean <- 200 r_boot_gear <- 500 b_boot <- 100 boot_one <- function(x) { s <- scan_ks(x) ab <- replicate(b_boot, scan_ks(x[sample.int(length(x), replace = TRUE)])["alpha"]) se_boot <- sd(ab) q_ab <- quantile(ab, c(0.025, 0.975), names = FALSE) c(wald = abs(s[["alpha"]] - alpha0) / se_boot > qnorm(0.975), perc = alpha0 < q_ab[1] | alpha0 > q_ab[2], ratio = se_boot / ((s[["alpha"]] - 1) / sqrt(s[["n_above"]])), alpha = s[["alpha"]], xmin = s[["xmin"]], se = se_boot, centre = mean(ab) - s[["alpha"]]) } set.seed(3346) boot_clean <- replicate(r_boot_clean, boot_one(rbounded_pareto(n_boot_draw, b_true, m_lo, Inf))) boot_gear <- replicate(r_boot_gear, boot_one(r_gear(n_boot_draw))) boot_size <- rbind(clean = rowMeans(boot_clean), gear = rowMeans(boot_gear)) boot_reps <- c(clean = r_boot_clean, gear = r_boot_gear) boot_se <- mcse(boot_size[, c("wald", "perc")], boot_reps) boot_mcse_max <- max(boot_se) # share of the gear-filtered excess over 0.05 that the bootstrap leaves in place boot_left <- (boot_size["gear", c("wald", "perc")] - 0.05) / (scan_size$gear_scan[1] - 0.05) boot_left_se <- sqrt(boot_se["gear", ]^2 + boot_left^2 * mcse(scan_size$gear_scan[1], r_scan)^2) / (scan_size$gear_scan[1] - 0.05) # clean samples: bootstrap rates below the chi-square(1) rate after the scan, in standard errors clean_gap_z <- (scan_size$clean_scan[1] - boot_size["clean", c("wald", "perc")]) / sqrt(boot_se["clean", ]^2 + mcse(scan_size$clean_scan[1], r_scan)^2) # where the gear-filtered Wald rejections come from g_alpha <- boot_gear["alpha", ] g_se <- boot_gear["se", ] g_rej <- boot_gear["wald", ] == 1 gear_diag <- c( bias_sd = (mean(g_alpha) - alpha0) / sd(g_alpha), wald_one_sd = mean(abs(g_alpha - alpha0) / sd(g_alpha) > qnorm(0.975)), cor_se = cor(g_alpha, g_se), cor_xmin = cor(g_alpha, boot_gear["xmin", ]), cor_xmin_se = cor(boot_gear["xmin", ], g_se), se_rejected = mean(g_se[g_rej]) / mean(g_se), below_truth = mean(g_alpha[g_rej] < alpha0), centre_sd = mean(boot_gear["centre", ]) / sd(g_alpha)) ``` With `r sprintf("%d", r_boot_clean)` clean and `r sprintf("%d", r_boot_gear)` gear-filtered samples and `r sprintf("%d", b_boot)` resamples per sample, the answer differs between the two kinds of data. On the clean power law both bootstrap rates are consistent with nominal: the Wald version rejects the true exponent in `r sprintf("%.3f", boot_size["clean", "wald"])` of the samples and the percentile version in `r sprintf("%.3f", boot_size["clean", "perc"])`, against `r sprintf("%.3f", scan_size$clean_scan[1])` for the chi-square(1) test after the scan at the same sample size, and neither is more than `r sprintf("%.1f", max(abs(boot_size["clean", c("wald", "perc")] - 0.05) / boot_se["clean", ]))` Monte Carlo standard errors from 0.05. The excess on clean data was small to begin with, though, and at this replication only the Wald rate is clearly below the chi-square(1) one (`r sprintf("%.1f", clean_gap_z["wald"])` standard errors of the difference, against `r sprintf("%.1f", clean_gap_z["perc"])` for the percentile rate). On the gear-filtered samples the bootstrap helps but does not finish the job: `r sprintf("%.3f", boot_size["gear", "wald"])` and `r sprintf("%.3f", boot_size["gear", "perc"])` against `r sprintf("%.3f", scan_size$gear_scan[1])`, both still above 0.05 by at least `r sprintf("%.1f", min((boot_size["gear", c("wald", "perc")] - 0.05) / boot_se["gear", ]))` Monte Carlo standard errors, so `r sprintf("%.0f", 100 * boot_left["wald"])` per cent of the excess over 0.05 is still there for the Wald version and `r sprintf("%.0f", 100 * boot_left["perc"])` per cent for the percentile version (each share with a Monte Carlo standard error of at most `r sprintf("%.2f", max(boot_left_se))`). The Monte Carlo standard error of the four bootstrap rates is at most `r sprintf("%.3f", boot_mcse_max)`. On average the bootstrap interval is wide enough. Over the `r sprintf("%d", r_scan)` samples of the previous chunk at n = `r sprintf("%d", n_scan[1])`, the real spread of the exponent after the scan is `r sprintf("%.2f", scan_size$clean_sd_ratio[1])` times the average formula standard error on clean samples and `r sprintf("%.2f", scan_size$gear_sd_ratio[1])` times on gear-filtered ones, and the bootstrap standard error is on average `r sprintf("%.2f", boot_size["clean", "ratio"])` and `r sprintf("%.2f", boot_size["gear", "ratio"])` times the formula's, so it covers that gap: the extra width is the uncertainty in the limit that the formula leaves out. On the gear-filtered samples the trouble is which samples get the narrow intervals. Across the `r sprintf("%d", r_boot_gear)` of them the fitted exponent is on average `r sprintf("%.2f", abs(gear_diag["bias_sd"]))` of its own standard deviation shallower than the truth, and a Wald test that used that one standard deviation for every sample would reject in `r sprintf("%.3f", gear_diag["wald_one_sd"])`. With each sample's own bootstrap standard error it rejects in `r sprintf("%.3f", boot_size["gear", "wald"])`, because the error and the interval move together: when the scan settles on a lower limit, more of the rolled-over body enters the fit, the exponent comes out shallower, and the bootstrap spread comes out narrower. Across samples the steepness of the fitted exponent correlates at `r sprintf("%.2f", gear_diag["cor_xmin"])` with the chosen limit and at `r sprintf("%.2f", gear_diag["cor_se"])` with its bootstrap standard error, and the chosen limit at `r sprintf("%.2f", gear_diag["cor_xmin_se"])` with the bootstrap standard error; `r if (gear_diag["below_truth"] == 1) sprintf("all %d", sum(g_rej)) else sprintf("%.0f per cent of the %d", 100 * gear_diag["below_truth"], sum(g_rej))` rejections are of an exponent shallower than the truth, and the rejected samples have a bootstrap standard error `r sprintf("%.2f", gear_diag["se_rejected"])` times the average. Resampling reproduces the shift rather than correcting it: the net offset of the resampled exponents from the sample's own estimate, averaged over the samples, is `r sprintf("%.2f", abs(gear_diag["centre_sd"]))` standard deviations, against the `r sprintf("%.2f", abs(gear_diag["bias_sd"]))` by which the estimates miss the truth, so a bootstrap estimate of the bias would find almost none of it. ```{r fig-lr-size} #| fig-cap: "Rejection rate of the true exponent at nominal 0.05: lower limit fixed, picked by the KS scan, and with the scan re-run inside a bootstrap." #| fig-alt: "Two stacked panels of points with horizontal error bars of two Monte Carlo standard errors on a rejection-rate axis, with a dashed gold line at the nominal 0.05; green points are samples of 300 animals and red points samples of 1200. In the upper panel, for a clean power law, the test with the limit fixed at the truth sits on or just right of the line, the chi-square(1) test after the KS scan sits near 0.09 at both sample sizes, and the bootstrap standard-error and percentile tests at 300 animals sit near 0.04 and 0.065 with bars that cross the line. In the lower panel, for gear-filtered samples, the chi-square(1) test after the scan sits near 0.2 at both sample sizes, the bootstrap standard-error test near 0.12 with its bar well clear of the line, and the bootstrap percentile test near 0.08 with its bar just clear of it." size_tab <- rbind( data.frame(sample = "clean power law", n_draw = scan_size$n_draw, test = "limit fixed at truth", rate = scan_size$clean_fixed, reps = r_scan), data.frame(sample = "clean power law", n_draw = scan_size$n_draw, test = "KS scan, chi-square(1)", rate = scan_size$clean_scan, reps = r_scan), data.frame(sample = "gear-filtered", n_draw = scan_size$n_draw, test = "KS scan, chi-square(1)", rate = scan_size$gear_scan, reps = r_scan), data.frame(sample = c("clean power law", "gear-filtered"), n_draw = n_boot_draw, test = "KS scan, bootstrap SE", rate = boot_size[, "wald"], reps = boot_reps), data.frame(sample = c("clean power law", "gear-filtered"), n_draw = n_boot_draw, test = "KS scan, bootstrap percentile", rate = boot_size[, "perc"], reps = boot_reps) ) size_tab$test <- factor(size_tab$test, levels = rev(c("limit fixed at truth", "KS scan, chi-square(1)", "KS scan, bootstrap SE", "KS scan, bootstrap percentile"))) size_tab$n_lab <- factor(sprintf("n = %d", size_tab$n_draw), levels = sprintf("n = %d", n_scan)) size_tab$se <- mcse(size_tab$rate, size_tab$reps) ggplot(size_tab, aes(rate, test, colour = n_lab)) + geom_vline(xintercept = 0.05, colour = te_gold, linetype = "dashed", linewidth = 0.7) + geom_errorbar(aes(xmin = rate - 2 * se, xmax = rate + 2 * se), width = 0.2, orientation = "y", linewidth = 0.6, position = position_dodge(width = 0.5)) + geom_point(size = 2.6, position = position_dodge(width = 0.5)) + facet_wrap(~ sample, ncol = 1, scales = "free_y") + scale_colour_manual(values = c(te_forest, te_rust), name = NULL) + labs(x = "rejection rate of the true exponent at nominal 0.05", y = NULL, title = "Picking the limit from the data inflates the test", subtitle = "bars: two Monte Carlo standard errors; dashed gold: nominal 0.05") + theme_datasheet() + theme(legend.position = "bottom") ``` In practice, then: if the lower limit is the sample minimum, the exact test above is the answer and needs nothing else. If the limit came from a scan, a chi-square(1) p-value at the chosen limit rejects a true exponent too often, a little on a clean power law and a lot when the gear has left a rolled-over body behind. The resample-and-rescan bootstrap is the better default: on a clean power law both of its rates are within Monte Carlo error of 0.05, but on gear-filtered data it removes only part of the excess, because the scan shifts the exponent, each resample reproduces the shift, and the most shifted samples get the narrowest intervals. The sensitivity curve still belongs in the report, because a test at one chosen limit says nothing about the others. These rates belong to one exponent, one scan grid and one retention curve, and a coarser grid, a larger minimum tail or a steeper retention curve will move them. The bootstrap was measured at one sample size and with `r sprintf("%d", b_boot)` resamples, against the 1000 or so Clauset and colleagues suggest; a noisier bootstrap standard error tends to make the Wald version reject more often, not less. ## What to report Say which spectrum the exponent belongs to. Counts per class, counts per unit mass and biomass per class differ by whole units of slope, biomass per unit mass shares the slope of counts per class, and a value quoted without that label cannot be compared with anything. Fit by likelihood, not by regression on binned counts. The closed form is one line, and so is its standard error when the lower limit is fixed in advance; when the limit comes from a scan, take the standard error from a bootstrap that re-runs the scan, as in the section above. The four hundred replicate comparison above shows the closed form beating every binned alternative on bias and on error at both sample sizes. Give the lower mass limit, the number of individuals above it, and what the exponent does as the limit moves. The limit is where the real disagreement between analyses lives, and a sensitivity curve costs one plot. A test of the exponent at a scanned limit needs the bootstrap, not the chi-square(1) table. If a binned figure is wanted, and it usually is, plot the binned spectrum and fit the likelihood. The two are not in competition once the line on the figure is drawn from the likelihood estimate rather than from a regression through the points. ## Honest limits The data here are generated from a power law, so every method is being scored on its ability to recover a truth that exists. Real body mass distributions are not exactly power laws, and when the shape is wrong the likelihood estimator is not saved by being efficient: it returns the exponent of the best fitting power law, which may describe nothing. Testing whether a power law is the right family at all is a separate exercise, and the site has covered how weak that kind of comparison can be when two families are close. The Kolmogorov-Smirnov criterion for the lower limit is a heuristic with known failure modes. It tends to pick a limit that discards more data than necessary when the departure from a power law is mild, and it has no defence at all against a distribution that is curved throughout, where it will find a local stretch that looks straight. The flat section of the curve above is doing more of the work than the minimum is. The selectivity example uses a smooth retention curve that is still short of one at the limits the scan chooses: at the median limits in the test above the gear keeps only `r sprintf("%.2f", scan_size$gear_retain[1])` and `r sprintf("%.2f", scan_size$gear_retain[2])` of the animals it meets, which is why part of the rolled-over body enters the fit. Gear whose retention is still climbing across the whole fitted range is worse again: it biases the exponent in a way no choice of lower limit can remove, because there is no mass above which the sample is a fair draw. The upper end is left alone throughout. A bounded power law was simulated and an unbounded one was fitted, which costs little here because the largest observed mass sits far below the bound. Samples whose largest individuals press against a real physiological or gear-imposed ceiling need the bounded likelihood, and the difference is then not small. Finally, `r sprintf("%d", n_rep)` replicates give each bias a Monte Carlo standard error of at most `r sprintf("%.4f", mc_se_max)` on the scale of the exponent. That is below every gap discussed above, including the `r sprintf("%.3f", bin_gap_big)` between six and ten bins at the larger sample size, so the ordering is solid. The fourth decimal place of any one of them is not. ## References Edwards AM, Robinson JPW, Plank MJ, Baum JK, Blanchard JL 2017 Methods in Ecology and Evolution 8(1):57-67 (10.1111/2041-210X.12641) White EP, Enquist BJ, Green JL 2008 Ecology 89(4):905-912 (10.1890/07-1288.1) Clauset A, Shalizi CR, Newman MEJ 2009 SIAM Review 51(4):661-703 (10.1137/070710111) Sheldon RW, Prakash A, Sutcliffe WH 1972 Limnology and Oceanography 17(3):327-340 (10.4319/lo.1972.17.3.0327) ## Related tutorials - [Peaks over threshold and the GPD](../peaks-over-threshold-and-the-gpd/) - [Allometry and log-log regression in R](../allometry-and-log-log-regression/) - [Gill-net selectivity and the fish it cannot see](../gill-net-selectivity-and-size-structure/) - [Levy walk or patchy habitat: which test to trust](../levy-walk-or-patchy-habitat/)