Capture heterogeneity: Mt, Mb and Mh in R

capture-recapture

abundance

model diagnostics

ecology tutorial

Closed capture-recapture in R: the failures of constant capture probability cost different amounts. Trap shyness inflates abundance, heterogeneity cuts it.

Author

Tidy Ecology

Published

2026-08-08

Modified

2026-10-15

Corrected 15 October 2026: the post said that a single violation can be diagnosed in a study of ordinary size; the one design measured is a well sampled study, and the claim is now limited to data of that quality. It also said that all four generators share a mean capture probability, which the trap shy one does not, and that heterogeneity did not bias the estimate in the mark resight post, which reports a small upward bias.

Corrected 7 October 2026: the Mt fit now uses a gradient search, since the simplex search stopped short of the maximum in some studies, which moves the trap shy selection slightly towards Mt. The post also said the number of animals seen is similar under all four generators, that the models AIC picks by mistake give similar estimates, and that the selection result would not survive a halved capture probability; it now says fewer are seen under heterogeneity, that trap shy studies sent to Mt keep most of the upward bias, and that the halved case, which was not run, would weaken the result.

Corrected 25 September 2026: the post said the data can show which assumption failed in a study of ordinary size. That holds while one assumption fails at a time; a new subsection, Two violations at once, measures trap shyness and heterogeneity together, where the comparison throws out the model that happens to be right, and the direction of the trap-response bias is now given for shyness and happiness separately.

Closed population capture-recapture starts from an assumption nobody believes: that every animal in the population has the same probability of being caught on every occasion. That is model M0, and it is the model whose arithmetic is short enough to write out. The alternatives form a standard family, and the three that matter are named after the way the assumption fails. Mt lets the probability differ between occasions. Mb lets it change after an animal has been caught once. Mh lets it differ between individuals.

The site already has M0 and the standard closed population estimators. What it does not have is a measurement of what happens when M0 is fitted to data generated by each of the others, which is the situation every real study is in, because none of the three violations can be ruled out from the field.

The three do not cost the same, and two of them push in opposite directions. This post generates data under each, fits all four models by writing the likelihoods out, and measures both the bias and whether model selection can find the right one.

Four generators

Six occasions, two hundred animals, and a capture history matrix.

library(ggplot2)

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          plot.title       = element_text(colour = te_ink, face = "bold"),
          axis.text        = element_text(colour = te_body))
}

n_occ  <- 6
n_true <- 200
p_flat  <- 0.30
p_by_occ <- c(0.15, 0.45, 0.20, 0.40, 0.25, 0.35)
p_first <- 0.35
p_again <- 0.15

gen_m0 <- function(N, p) matrix(rbinom(N * n_occ, 1, p), N, n_occ)
gen_mt <- function(N, pv) sapply(pv, function(p) rbinom(N, 1, p))
gen_mb <- function(N, p, recap) {
  X <- matrix(0, N, n_occ); caught <- rep(FALSE, N)
  for (j in seq_len(n_occ)) {
    X[, j] <- rbinom(N, 1, ifelse(caught, recap, p))
    caught <- caught | X[, j] == 1
  }
  X
}
gen_mh <- function(N, a, b) {
  pi_ind <- rbeta(N, a, b)
  matrix(rbinom(N * n_occ, 1, rep(pi_ind, n_occ)), N, n_occ)
}

The constant, time varying and heterogeneous generators share a mean capture probability of 0.30; the trap shy one has a first capture probability a little above it and a recapture probability well below it. The number of animals ever seen is similar under the constant, time varying and trap shy generators; under heterogeneity it is lower, because animals with a low capture probability are often never caught at all. What differs is where the variation sits.

Four likelihoods

Each model has the same structure: a combinatorial term for which animals were seen, the probability of every observed history, and the probability of the all zero history raised to the number never seen. The unknown abundance enters as the number of animals that produced the all zero history.

nll_m0 <- function(par, x, n_seen) {
  N <- n_seen + exp(par[1]); p <- plogis(par[2])
  -(lgamma(N + 1) - lgamma(N - n_seen + 1) +
      sum(x) * log(p) + (N * n_occ - sum(x)) * log(1 - p))
}
nll_mt <- function(par, X, n_seen) {
  N <- n_seen + exp(par[1]); pv <- plogis(par[-1]); cj <- colSums(X)
  -(lgamma(N + 1) - lgamma(N - n_seen + 1) +
      sum(cj * log(pv) + (n_seen - cj) * log(1 - pv)) +
      (N - n_seen) * sum(log(1 - pv)))
}
nll_mb <- function(par, X, n_seen) {
  N <- n_seen + exp(par[1]); p <- plogis(par[2]); recap <- plogis(par[3])
  first <- apply(X, 1, function(h) which(h == 1)[1])
  n_first  <- n_seen
  risk_p   <- sum(first) + (N - n_seen) * n_occ
  n_recap  <- sum(X) - n_seen
  risk_c   <- sum(n_occ - first)
  -(lgamma(N + 1) - lgamma(N - n_seen + 1) +
      n_first * log(p) + (risk_p - n_first) * log(1 - p) +
      n_recap * log(recap) + (risk_c - n_recap) * log(1 - recap))
}
nll_mh <- function(par, x, n_seen) {
  N <- n_seen + exp(par[1]); a <- exp(par[2]); b <- exp(par[3])
  hist_p <- function(k) lbeta(a + k, b + n_occ - k) - lbeta(a, b)
  -(lgamma(N + 1) - lgamma(N - n_seen + 1) +
      sum(hist_p(x)) + (N - n_seen) * hist_p(0))
}

fit_all <- function(X) {
  seen <- X[rowSums(X) > 0, , drop = FALSE]
  n_seen <- nrow(seen); x <- rowSums(seen)
  f0 <- optim(c(log(5), 0), nll_m0, x = x, n_seen = n_seen)
  ft <- optim(c(log(5), rep(0, n_occ)), nll_mt, X = seen, n_seen = n_seen,
              method = "BFGS", control = list(maxit = 2000))
  fb <- optim(c(log(5), 0, 0), nll_mb, X = seen, n_seen = n_seen,
              control = list(maxit = 2000))
  fh <- optim(c(log(5), 0, 0), nll_mh, x = x, n_seen = n_seen,
              control = list(maxit = 3000))
  singles <- sum(x == 1)
  c(M0 = n_seen + exp(f0$par[1]), Mt = n_seen + exp(ft$par[1]),
    Mb = n_seen + exp(fb$par[1]), Mh = n_seen + exp(fh$par[1]),
    jackknife = n_seen + singles * (n_occ - 1) / n_occ, seen = n_seen,
    aic_M0 = 2 * f0$value + 4, aic_Mt = 2 * ft$value + 2 * (1 + n_occ),
    aic_Mb = 2 * fb$value + 6, aic_Mh = 2 * fh$value + 6)
}

The jackknife on the last line is the classical estimator for individual heterogeneity from Burnham and Overton, and it is one line: add to the number of animals seen a correction driven by how many were caught exactly once.

What each violation costs

n_rep <- 60
scenarios <- list(
  "constant (M0)"          = function() gen_m0(n_true, p_flat),
  "time varying (Mt)"      = function() gen_mt(n_true, p_by_occ),
  "trap shy (Mb)"          = function() gen_mb(n_true, p_first, p_again),
  "heterogeneous (Mh)"     = function() gen_mh(n_true, 1.2, 2.8))

set.seed(3)
results <- lapply(scenarios, function(g) as.data.frame(t(replicate(n_rep, fit_all(g())))))

summary_row <- function(r) {
  aic <- as.matrix(r[, c("aic_M0", "aic_Mt", "aic_Mb", "aic_Mh")])
  picked <- c("M0", "Mt", "Mb", "Mh")[apply(aic, 1, which.min)]
  data.frame(M0 = median(r$M0), Mt = median(r$Mt), Mb = median(r$Mb),
             Mh = median(r$Mh), jackknife = median(r$jackknife),
             seen = median(r$seen),
             chosen = round(100 * mean(picked == "M0"), 0))
}
tbl <- do.call(rbind, lapply(results, summary_row))
tbl$truth_recovered <- round(100 * (tbl$M0 / n_true - 1), 1)
print(round(tbl[, c("seen", "M0", "Mt", "Mb", "Mh", "jackknife",
                    "truth_recovered")], 1))

                    seen    M0    Mt    Mb    Mh jackknife truth_recovered
constant (M0)      177.0 201.5 201.4 200.1 203.4     228.9             0.8
time varying (Mt)  177.5 200.3 197.4 215.7 201.3     228.0             0.1
trap shy (Mb)      184.5 266.3 262.7 198.4 271.1     272.6            33.2
heterogeneous (Mh) 148.0 156.1 156.0 157.6 192.6     187.8           -22.0

Read the M0 column against a true abundance of 200.

Time variation is nearly free. Capture probability swinging between 0.15 and 0.45 across the six occasions leaves the constant probability estimator at 200, an error of 0.1 per cent. The reason is that the number of animals never seen depends on the product of the per occasion failure probabilities, and that product is not much changed by moving probability between occasions.

A trap response is not free, and when it is shyness it goes up. With a recapture probability of 0.15 against a first capture probability of 0.35, M0 reports 266, an overestimate of 33.2 per cent. Animals that avoid the traps after their first capture look like animals that were never there to catch.

Heterogeneity is not free and it goes down. The Bayesian closed capture-recapture post on this site states that this is the failure worth worrying about; here is the number: M0 reports 156, an underestimate of 22.0 per cent. The mechanism is that the animals that are easy to catch are caught repeatedly, which makes the population look thoroughly sampled, while the shy animals contribute almost nothing and are quietly assumed not to exist.

long <- do.call(rbind, lapply(names(results), function(nm) {
  r <- results[[nm]]
  do.call(rbind, lapply(c("M0", "Mt", "Mb", "Mh", "jackknife"), function(k)
    data.frame(scenario = nm, fitted = k, value = r[[k]])))
}))
long$scenario <- factor(long$scenario, levels = names(scenarios))
long$fitted   <- factor(long$fitted,
                        levels = c("M0", "Mt", "Mb", "Mh", "jackknife"))

ggplot(long, aes(fitted, value)) +
  geom_hline(yintercept = n_true, colour = te_ink, linetype = "dashed",
             linewidth = 0.6) +
  geom_boxplot(fill = te_gold, colour = te_ink, alpha = 0.7,
               outlier.size = 0.5, linewidth = 0.35) +
  facet_wrap(~ scenario) +
  coord_cartesian(ylim = c(100, 330)) +
  labs(x = "model fitted", y = "estimated abundance",
       title = "Which violation, and which way it pushes") +
  theme_datasheet() +
  theme(strip.text = element_text(colour = te_ink, face = "bold"))

Four panels of boxplots. Under constant capture every model sits on the dashed line except the jackknife, which sits above it. Under time variation Mb is displaced upward and the jackknife more so, while M0, Mt and Mh stay on the line. Under a trap response every model except Mb sits well above the line. Under heterogeneity M0, Mt and Mb sit well below it, Mh straddles it and the jackknife sits a little below it. — Figure 1: Estimated abundance from the four models and the jackknife, for data generated under each of the four models, over 60 replicates. The dashed line is the true abundance of 200.

Can model selection find the right one

If the bias depends this much on which violation is present, everything rests on being able to tell them apart from the capture histories.

pick_rate <- function(r) {
  aic <- as.matrix(r[, c("aic_M0", "aic_Mt", "aic_Mb", "aic_Mh")])
  picked <- c("M0", "Mt", "Mb", "Mh")[apply(aic, 1, which.min)]
  round(100 * table(factor(picked, levels = c("M0", "Mt", "Mb", "Mh"))) / nrow(r))
}
selection <- do.call(rbind, lapply(results, pick_rate))
selection

                   M0  Mt Mb Mh
constant (M0)      75   7 15  3
time varying (Mt)   0 100  0  0
trap shy (Mb)       0   5 95  0
heterogeneous (Mh)  0   2  0 98

It can, on data this size, while only one assumption fails at a time. The generating model has the lowest AIC in the large majority of replicates in every scenario, though a wrong pick can be costly: the few trap shy studies that AIC sends to Mt get an estimate about as inflated as M0’s, because Mt has no term for the trap response. Two hundred animals over six occasions with a mean capture probability of 0.30 is a generous study, and for a single violation the message is a positive one: the diagnosis is available in the data, and skipping the comparison is a choice rather than a limitation. The next subsection shows where that stops.

Two violations at once

Every scenario above breaks one assumption. Real animals can be unequal and trap shy at the same time, and Mb has a property that decides what happens then: its abundance estimate uses only the first capture of each animal. The recaptures inform the recapture probability and nothing else, so for abundance Mb is a removal estimator, and removal under individual heterogeneity runs low for the reason worked out in Declaring eradication after empty traps: the animals still uncaught after several occasions are the hard ones.

The generator below draws each animal’s first capture probability from the same Beta(1.2, 2.8) as the heterogeneous scenario and multiplies it by a fixed factor once the animal has been caught: 0.43, the ratio of the trap shy scenario, for shy animals and 1.5 for trap happy ones. The last lines fit the removal likelihood to the expected first capture counts, N B(a + 1, b + j - 1) / B(a, b) on occasion j; the result is the answer Mb gives when the first captures of a 200 animal study fall exactly on their expected values.

gen_mbh <- function(N, a, b, mult) {
  pi_ind <- rbeta(N, a, b); recap <- pmin(1, mult * pi_ind)
  X <- matrix(0, N, n_occ); caught <- rep(FALSE, N)
  for (j in seq_len(n_occ)) {
    X[, j] <- rbinom(N, 1, ifelse(caught, recap, pi_ind))
    caught <- caught | X[, j] == 1
  }
  X
}
fit_cover <- function(X) {
  seen <- X[rowSums(X) > 0, , drop = FALSE]; n_seen <- nrow(seen); x <- rowSums(seen)
  f0 <- optim(c(log(5), 0), nll_m0, x = x, n_seen = n_seen)
  se <- sqrt(solve(optimHess(f0$par, nll_m0, x = x, n_seen = n_seen))[1, 1])
  c(fit_all(X), m0_covers = abs(log(n_true - n_seen) - f0$par[1]) < 1.96 * se)
}
picked_value <- function(r) {
  aic <- as.matrix(r[, c("aic_M0", "aic_Mt", "aic_Mb", "aic_Mh")])
  picked <- c("M0", "Mt", "Mb", "Mh")[apply(aic, 1, which.min)]
  as.matrix(r)[cbind(seq_len(nrow(r)), match(picked, names(r)))]
}

set.seed(58)
both <- list(
  "shy, heterogeneous"   = as.data.frame(t(replicate(n_rep,
    fit_cover(gen_mbh(n_true, 1.2, 2.8, p_again / p_first))))),
  "happy, heterogeneous" = as.data.frame(t(replicate(n_rep,
    fit_cover(gen_mbh(n_true, 1.2, 2.8, 1.5))))))
selection_both <- do.call(rbind, lapply(both, pick_rate))
both_est <- sapply(both, function(r) c(picked = median(picked_value(r)),
  M0 = median(r$M0), Mb = median(r$Mb), Mh = median(r$Mh),
  M0_cover = mean(r$m0_covers)))

first_cap <- n_true * exp(lbeta(1.2 + 1, 2.8 + seq_len(n_occ) - 1) - lbeta(1.2, 2.8))
nll_removal <- function(par, u) {
  n_u <- sum(u); N <- n_u + exp(par[1]); p <- plogis(par[2])
  -(lgamma(N + 1) - lgamma(N - n_u + 1) + n_u * log(p) +
      (sum((seq_len(n_occ) - 1) * u) + (N - n_u) * n_occ) * log(1 - p))
}
fit_limit <- optim(c(log(20), 0), nll_removal, u = first_cap)
mb_limit <- sum(first_cap) + exp(fit_limit$par[1])
selection_both

                     M0 Mt Mb  Mh
shy, heterogeneous    0 12 87   2
happy, heterogeneous  0  0  0 100

round(both_est, 2)

         shy, heterogeneous happy, heterogeneous
picked               158.54               202.00
M0                   202.70               150.12
Mb                   158.44               158.76
Mh                   231.79               202.00
M0_cover               0.97                 0.02

round(mb_limit, 1)

[1] 158.8

With shy animals that also differ, AIC picks Mb in 87 per cent of 60 studies, and the estimate of whichever model AIC picks has a median of 159 against a true 200. That shortfall owes nothing to the trap response. The removal fit to expected first captures gives 159, and no recapture factor appears anywhere in it; Mb on the trap happy data has a median of 159, the same number. The comparison picks the model that explains the recaptures and inherits the heterogeneity bias through the first captures.

AIC picks M0 in none of these studies. Its median is 203, and its Wald interval covers the true abundance in 97 per cent of them. At these two strengths the upward push of shyness and the downward pull of heterogeneity cancel. Nothing in the histories tells the analyst that they have, and the model the comparison throws out is the one that happens to be right.

With trap happy animals the extra recaptures look like heterogeneity: AIC picks Mh in 100 per cent of studies, and the picked estimate has a median of 202. A model with both effects, Mbh, belongs to the mixture family of Pledger (2000); it is not in the candidate set here, and this post does not measure how it does.

How much heterogeneity is too much

The Mh scenario above used one particular amount of individual variation. Holding the mean capture probability fixed and varying its spread shows how quickly the constant probability estimator gives way.

nll_two_pt <- function(par, x, n_seen) {
  N <- n_seen + exp(par[1]); w <- plogis(par[2])
  p1 <- plogis(par[3]); p2 <- plogis(par[4])
  hist_p <- function(k) log(w * p1^k * (1 - p1)^(n_occ - k) +
                              (1 - w) * p2^k * (1 - p2)^(n_occ - k))
  -(lgamma(N + 1) - lgamma(N - n_seen + 1) +
      sum(hist_p(x)) + (N - n_seen) * hist_p(0))
}
fit_hetero <- function(X) {
  s <- X[rowSums(X) > 0, , drop = FALSE]; n_seen <- nrow(s); x <- rowSums(s)
  f0 <- optim(c(log(5), 0), nll_m0, x = x, n_seen = n_seen)
  fh <- optim(c(log(5), 0, 0), nll_mh, x = x, n_seen = n_seen,
              control = list(maxit = 3000))
  f2 <- optim(c(log(5), 0, -1, 1), nll_two_pt, x = x, n_seen = n_seen,
              control = list(maxit = 5000))
  c(M0 = n_seen + exp(f0$par[1]), beta = n_seen + exp(fh$par[1]),
    two_pt = n_seen + exp(f2$par[1]),
    jackknife = n_seen + sum(x == 1) * (n_occ - 1) / n_occ,
    nll_beta = fh$value, nll_two = f2$value)
}

set.seed(21)
mean_p <- p_flat
shapes <- c(50, 6, 2.4, 1.2, 0.6)
sweep_h <- do.call(rbind, lapply(shapes, function(a) {
  b  <- a * (1 - mean_p) / mean_p
  cv <- sqrt((a * b) / ((a + b)^2 * (a + b + 1))) / mean_p
  r  <- as.data.frame(t(replicate(50, fit_hetero(gen_mh(n_true, a, b)))))
  data.frame(cv = cv, M0 = median(r$M0), beta = median(r$beta),
             two_pt = median(r$two_pt), jackknife = median(r$jackknife))
}))
print(round(sweep_h, 1))

   cv    M0  beta two_pt jackknife
1 0.1 197.7 200.7  200.6     225.0
2 0.3 189.4 201.6  201.2     220.8
3 0.5 172.1 196.0  183.9     203.3
4 0.7 156.8 188.7  172.3     187.9
5 0.9 134.8 198.0  151.6     163.8

At a coefficient of variation of 0.33 in capture probability, which is modest for a real population, the constant probability estimate is already 5 per cent low. At 0.88 it is 33 per cent low.

The jackknife is the standard insurance against this, and the first row is its premium. On data with essentially no heterogeneity it returns 225 against a true 200, an overestimate of 12 per cent. Insurance that costs a 12 per cent bias when the risk is absent is worth buying only when the risk is real, which is an argument for looking rather than for defaulting either way.

d_sw <- do.call(rbind, lapply(c("M0", "beta", "two_pt", "jackknife"), function(k)
  data.frame(cv = sweep_h$cv, value = sweep_h[[k]], estimator = k)))
d_sw$estimator <- factor(d_sw$estimator, levels = c("M0", "beta", "two_pt", "jackknife"),
                         labels = c("M0, constant p", "Mh, beta mixture",
                                    "Mh, two point mixture", "jackknife"))

ggplot(d_sw, aes(cv, value, colour = estimator)) +
  geom_hline(yintercept = n_true, colour = te_ink, linetype = "dashed",
             linewidth = 0.6) +
  geom_line(linewidth = 1) + geom_point(size = 2.2) +
  scale_colour_manual(values = c(te_rust, te_forest, te_gold, "#7b6a4f"),
                      name = NULL) +
  labs(x = "coefficient of variation of individual capture probability",
       y = "median estimated abundance",
       title = "Where the constant probability estimate goes") +
  theme_datasheet() +
  theme(legend.position = "bottom", legend.text = element_text(size = 9))

Four lines against increasing variation in capture probability. The constant probability line falls steadily from 200 to about 135. The jackknife starts at 225 and falls to about 164. The beta mixture stays near the dashed line at 200. — Figure 2: Median estimated abundance against the coefficient of variation of individual capture probability, holding the mean at 0.30, over 50 replicates per point.

Two heterogeneity models, one dataset, two answers

The two_pt column of that sweep is a second Mh model: instead of a beta distribution of capture probabilities it assumes two kinds of animal, easy and hard, in unknown proportions. Both are respectable choices, and on the same data they do not agree.

set.seed(404)
pair <- as.data.frame(t(replicate(40, fit_hetero(gen_mh(n_true, 1.2, 2.8)))))
gap_n   <- median(abs(pair$beta - pair$two_pt))
gap_nll <- median(abs(pair$nll_beta - pair$nll_two))

Over 40 datasets the two mixing distributions differ by a median of 27 animals, which is 13 per cent of the true abundance, while their log likelihoods differ by a median of 0.73. The data cannot separate them and the answers are not the same. Link proved in 2003 that this is not a numerical accident: different mixing distributions can give identical probabilities for every observable capture frequency and still imply different numbers of animals never caught. The unseen part of the population is where the disagreement lives, and no amount of fieldwork on the seen part settles it.

ggplot(pair, aes(beta, two_pt)) +
  geom_abline(slope = 1, intercept = 0, colour = te_ink, linetype = "dashed",
              linewidth = 0.6) +
  geom_hline(yintercept = n_true, colour = te_gold, linewidth = 0.7) +
  geom_vline(xintercept = n_true, colour = te_gold, linewidth = 0.7) +
  geom_point(colour = te_forest, size = 2.2, alpha = 0.85) +
  labs(x = "abundance from the beta mixture",
       y = "abundance from the two point mixture",
       title = "Same histories, two defensible answers") +
  theme_datasheet()

Scatter of forty points, most of them well below the dashed equality line, so the two point mixture returns a lower abundance than the beta mixture on nearly every dataset. Solid lines mark the true abundance of 200 on both axes: the cloud straddles the vertical line and sits mostly below the horizontal one. — Figure 3: Abundance from a beta mixture against abundance from a two point mixture, fitted to the same 40 simulated datasets. The dashed line marks equality and the solid lines mark the true abundance.

What to report

Fit more than M0 and report the comparison. The measurement above says a single violation can be diagnosed in a study as well sampled as this one, so a paper with data of that quality that reports only M0 has declined to look rather than been unable to. The comparison is evidence about which assumption fails, not proof: when two fail together it can settle firmly on a model that is still wrong.

Say which direction the residual risk runs. Trap shyness and individual heterogeneity bias abundance in opposite directions (trap happiness pushes the same way as heterogeneity), so a study that has ruled out one of them has not narrowed the interval on the other.

Report the number of animals caught exactly once alongside the estimate. It is the statistic every heterogeneity correction is built from, and a reader can see from it how much of the answer is extrapolation.

If a heterogeneity model is used, say which mixing distribution and acknowledge that the choice moves the answer. Reporting a beta mixture estimate to three significant figures without that sentence overstates what the data contain.

Honest limits

The population is genuinely closed here: no births, no deaths, no movement in or out. Real six occasion studies are closed by assumption and the assumption is usually a compromise about the length of the trapping session. Losses to the population and a trap response are also hard to tell apart from capture histories alone.

Outside the one subsection on two violations at once, only one violation acts at a time, and that subsection tries a single strength of heterogeneity with two trap responses. The named combinations Mth, Mbh and Mtbh exist because real studies have several at once; none of them was fitted here.

The behavioural response in the single violation scenarios is trap shyness, where recapture is less likely than first capture. Trap happiness is at least as common with baited traps and reverses the sign, so the direction reported for the trap shy scenario is the direction of that scenario, not a property of behavioural response as such.

The two hundred animals and mean capture probability of 0.30 make this a well sampled study. Halve the capture probability and the number of animals never seen grows, every estimator’s variance grows with it, and the model selection result reported here would weaken, though that case was not run for this post.

There is a reconciliation to make with the mark resight post on this site, where individual heterogeneity left the point estimate with only a small upward bias but destroyed the coverage of the interval. The difference is structural. In mark resight the number of marked animals is known, so heterogeneity affects how well the marked fraction is estimated. Here abundance is inferred from the shape of the capture frequency distribution, and heterogeneity changes that shape, so it moves the estimate itself.

References

Burnham KP, Overton WS 1978 Biometrika 65(3):625-633 (10.1093/biomet/65.3.625)

Pledger S 2000 Biometrics 56(2):434-442 (10.1111/j.0006-341X.2000.00434.x)

Link WA 2003 Biometrics 59(4):1123-1130 (10.1111/j.0006-341X.2003.00129.x)

--- title: "Capture heterogeneity: Mt, Mb and Mh in R" description: "Closed capture-recapture in R: the failures of constant capture probability cost different amounts. Trap shyness inflates abundance, heterogeneity cuts it." date: "2026-08-08 13:00" date-modified: "2026-10-15" categories: [R, capture-recapture, abundance, model diagnostics, ecology tutorial] image: thumbnail.png image-alt: "Estimated abundance against the amount of individual variation in capture probability on a warm off-white panel. The constant probability line falls steadily below the true value of 200, the beta mixture line stays on it throughout, and the jackknife line starts well above it before crossing down through it." --- *Corrected 15 October 2026: the post said that a single violation can be diagnosed in a study of ordinary size; the one design measured is a well sampled study, and the claim is now limited to data of that quality. It also said that all four generators share a mean capture probability, which the trap shy one does not, and that heterogeneity did not bias the estimate in the mark resight post, which reports a small upward bias.* *Corrected 7 October 2026: the Mt fit now uses a gradient search, since the simplex search stopped short of the maximum in some studies, which moves the trap shy selection slightly towards Mt. The post also said the number of animals seen is similar under all four generators, that the models AIC picks by mistake give similar estimates, and that the selection result would not survive a halved capture probability; it now says fewer are seen under heterogeneity, that trap shy studies sent to Mt keep most of the upward bias, and that the halved case, which was not run, would weaken the result.* *Corrected 25 September 2026: the post said the data can show which assumption failed in a study of ordinary size. That holds while one assumption fails at a time; a new subsection, Two violations at once, measures trap shyness and heterogeneity together, where the comparison throws out the model that happens to be right, and the direction of the trap-response bias is now given for shyness and happiness separately.* Closed population capture-recapture starts from an assumption nobody believes: that every animal in the population has the same probability of being caught on every occasion. That is model M0, and it is the model whose arithmetic is short enough to write out. The alternatives form a standard family, and the three that matter are named after the way the assumption fails. Mt lets the probability differ between occasions. Mb lets it change after an animal has been caught once. Mh lets it differ between individuals. The site already has M0 and the standard closed population estimators. What it does not have is a measurement of what happens when M0 is fitted to data generated by each of the others, which is the situation every real study is in, because none of the three violations can be ruled out from the field. The three do not cost the same, and two of them push in opposite directions. This post generates data under each, fits all four models by writing the likelihoods out, and measures both the bias and whether model selection can find the right one. ## Four generators Six occasions, two hundred animals, and a capture history matrix. ```{r setup} #| message: false #| warning: false library(ggplot2) te_paper <- "#f5f4ee" te_ink <- "#16241d" te_body <- "#2c3a31" te_forest <- "#275139" te_rust <- "#b5534e" te_gold <- "#c9b458" te_line <- "#dad9ca" theme_datasheet <- function() { theme_minimal(base_size = 12) + theme(plot.background = element_rect(fill = te_paper, colour = NA), panel.background = element_rect(fill = te_paper, colour = NA), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), text = element_text(colour = te_body), plot.title = element_text(colour = te_ink, face = "bold"), axis.text = element_text(colour = te_body)) } n_occ <- 6 n_true <- 200 p_flat <- 0.30 p_by_occ <- c(0.15, 0.45, 0.20, 0.40, 0.25, 0.35) p_first <- 0.35 p_again <- 0.15 gen_m0 <- function(N, p) matrix(rbinom(N * n_occ, 1, p), N, n_occ) gen_mt <- function(N, pv) sapply(pv, function(p) rbinom(N, 1, p)) gen_mb <- function(N, p, recap) { X <- matrix(0, N, n_occ); caught <- rep(FALSE, N) for (j in seq_len(n_occ)) { X[, j] <- rbinom(N, 1, ifelse(caught, recap, p)) caught <- caught | X[, j] == 1 } X } gen_mh <- function(N, a, b) { pi_ind <- rbeta(N, a, b) matrix(rbinom(N * n_occ, 1, rep(pi_ind, n_occ)), N, n_occ) } ``` The constant, time varying and heterogeneous generators share a mean capture probability of `r sprintf("%.2f", p_flat)`; the trap shy one has a first capture probability a little above it and a recapture probability well below it. The number of animals ever seen is similar under the constant, time varying and trap shy generators; under heterogeneity it is lower, because animals with a low capture probability are often never caught at all. What differs is where the variation sits. ## Four likelihoods Each model has the same structure: a combinatorial term for which animals were seen, the probability of every observed history, and the probability of the all zero history raised to the number never seen. The unknown abundance enters as the number of animals that produced the all zero history. ```{r likelihoods} nll_m0 <- function(par, x, n_seen) { N <- n_seen + exp(par[1]); p <- plogis(par[2]) -(lgamma(N + 1) - lgamma(N - n_seen + 1) + sum(x) * log(p) + (N * n_occ - sum(x)) * log(1 - p)) } nll_mt <- function(par, X, n_seen) { N <- n_seen + exp(par[1]); pv <- plogis(par[-1]); cj <- colSums(X) -(lgamma(N + 1) - lgamma(N - n_seen + 1) + sum(cj * log(pv) + (n_seen - cj) * log(1 - pv)) + (N - n_seen) * sum(log(1 - pv))) } nll_mb <- function(par, X, n_seen) { N <- n_seen + exp(par[1]); p <- plogis(par[2]); recap <- plogis(par[3]) first <- apply(X, 1, function(h) which(h == 1)[1]) n_first <- n_seen risk_p <- sum(first) + (N - n_seen) * n_occ n_recap <- sum(X) - n_seen risk_c <- sum(n_occ - first) -(lgamma(N + 1) - lgamma(N - n_seen + 1) + n_first * log(p) + (risk_p - n_first) * log(1 - p) + n_recap * log(recap) + (risk_c - n_recap) * log(1 - recap)) } nll_mh <- function(par, x, n_seen) { N <- n_seen + exp(par[1]); a <- exp(par[2]); b <- exp(par[3]) hist_p <- function(k) lbeta(a + k, b + n_occ - k) - lbeta(a, b) -(lgamma(N + 1) - lgamma(N - n_seen + 1) + sum(hist_p(x)) + (N - n_seen) * hist_p(0)) } fit_all <- function(X) { seen <- X[rowSums(X) > 0, , drop = FALSE] n_seen <- nrow(seen); x <- rowSums(seen) f0 <- optim(c(log(5), 0), nll_m0, x = x, n_seen = n_seen) ft <- optim(c(log(5), rep(0, n_occ)), nll_mt, X = seen, n_seen = n_seen, method = "BFGS", control = list(maxit = 2000)) fb <- optim(c(log(5), 0, 0), nll_mb, X = seen, n_seen = n_seen, control = list(maxit = 2000)) fh <- optim(c(log(5), 0, 0), nll_mh, x = x, n_seen = n_seen, control = list(maxit = 3000)) singles <- sum(x == 1) c(M0 = n_seen + exp(f0$par[1]), Mt = n_seen + exp(ft$par[1]), Mb = n_seen + exp(fb$par[1]), Mh = n_seen + exp(fh$par[1]), jackknife = n_seen + singles * (n_occ - 1) / n_occ, seen = n_seen, aic_M0 = 2 * f0$value + 4, aic_Mt = 2 * ft$value + 2 * (1 + n_occ), aic_Mb = 2 * fb$value + 6, aic_Mh = 2 * fh$value + 6) } ``` The jackknife on the last line is the classical estimator for individual heterogeneity from Burnham and Overton, and it is one line: add to the number of animals seen a correction driven by how many were caught exactly once. ## What each violation costs ```{r scenarios} #| warning: false n_rep <- 60 scenarios <- list( "constant (M0)" = function() gen_m0(n_true, p_flat), "time varying (Mt)" = function() gen_mt(n_true, p_by_occ), "trap shy (Mb)" = function() gen_mb(n_true, p_first, p_again), "heterogeneous (Mh)" = function() gen_mh(n_true, 1.2, 2.8)) set.seed(3) results <- lapply(scenarios, function(g) as.data.frame(t(replicate(n_rep, fit_all(g()))))) summary_row <- function(r) { aic <- as.matrix(r[, c("aic_M0", "aic_Mt", "aic_Mb", "aic_Mh")]) picked <- c("M0", "Mt", "Mb", "Mh")[apply(aic, 1, which.min)] data.frame(M0 = median(r$M0), Mt = median(r$Mt), Mb = median(r$Mb), Mh = median(r$Mh), jackknife = median(r$jackknife), seen = median(r$seen), chosen = round(100 * mean(picked == "M0"), 0)) } tbl <- do.call(rbind, lapply(results, summary_row)) tbl$truth_recovered <- round(100 * (tbl$M0 / n_true - 1), 1) print(round(tbl[, c("seen", "M0", "Mt", "Mb", "Mh", "jackknife", "truth_recovered")], 1)) ``` Read the M0 column against a true abundance of `r sprintf("%d", n_true)`. Time variation is nearly free. Capture probability swinging between `r sprintf("%.2f", min(p_by_occ))` and `r sprintf("%.2f", max(p_by_occ))` across the six occasions leaves the constant probability estimator at `r sprintf("%.0f", tbl["time varying (Mt)", "M0"])`, an error of `r sprintf("%.1f", tbl["time varying (Mt)", "truth_recovered"])` per cent. The reason is that the number of animals never seen depends on the product of the per occasion failure probabilities, and that product is not much changed by moving probability between occasions. A trap response is not free, and when it is shyness it goes up. With a recapture probability of `r sprintf("%.2f", p_again)` against a first capture probability of `r sprintf("%.2f", p_first)`, M0 reports `r sprintf("%.0f", tbl["trap shy (Mb)", "M0"])`, an overestimate of `r sprintf("%.1f", tbl["trap shy (Mb)", "truth_recovered"])` per cent. Animals that avoid the traps after their first capture look like animals that were never there to catch. Heterogeneity is not free and it goes down. The Bayesian closed capture-recapture post on this site states that this is the failure worth worrying about; here is the number: M0 reports `r sprintf("%.0f", tbl["heterogeneous (Mh)", "M0"])`, an underestimate of `r sprintf("%.1f", abs(tbl["heterogeneous (Mh)", "truth_recovered"]))` per cent. The mechanism is that the animals that are easy to catch are caught repeatedly, which makes the population look thoroughly sampled, while the shy animals contribute almost nothing and are quietly assumed not to exist. ```{r fig-bias} #| fig-cap: "Estimated abundance from the four models and the jackknife, for data generated under each of the four models, over 60 replicates. The dashed line is the true abundance of 200." #| fig-alt: "Four panels of boxplots. Under constant capture every model sits on the dashed line except the jackknife, which sits above it. Under time variation Mb is displaced upward and the jackknife more so, while M0, Mt and Mh stay on the line. Under a trap response every model except Mb sits well above the line. Under heterogeneity M0, Mt and Mb sit well below it, Mh straddles it and the jackknife sits a little below it." #| fig-width: 7.4 #| fig-height: 4.6 long <- do.call(rbind, lapply(names(results), function(nm) { r <- results[[nm]] do.call(rbind, lapply(c("M0", "Mt", "Mb", "Mh", "jackknife"), function(k) data.frame(scenario = nm, fitted = k, value = r[[k]]))) })) long$scenario <- factor(long$scenario, levels = names(scenarios)) long$fitted <- factor(long$fitted, levels = c("M0", "Mt", "Mb", "Mh", "jackknife")) ggplot(long, aes(fitted, value)) + geom_hline(yintercept = n_true, colour = te_ink, linetype = "dashed", linewidth = 0.6) + geom_boxplot(fill = te_gold, colour = te_ink, alpha = 0.7, outlier.size = 0.5, linewidth = 0.35) + facet_wrap(~ scenario) + coord_cartesian(ylim = c(100, 330)) + labs(x = "model fitted", y = "estimated abundance", title = "Which violation, and which way it pushes") + theme_datasheet() + theme(strip.text = element_text(colour = te_ink, face = "bold")) ``` ## Can model selection find the right one If the bias depends this much on which violation is present, everything rests on being able to tell them apart from the capture histories. ```{r selection} pick_rate <- function(r) { aic <- as.matrix(r[, c("aic_M0", "aic_Mt", "aic_Mb", "aic_Mh")]) picked <- c("M0", "Mt", "Mb", "Mh")[apply(aic, 1, which.min)] round(100 * table(factor(picked, levels = c("M0", "Mt", "Mb", "Mh"))) / nrow(r)) } selection <- do.call(rbind, lapply(results, pick_rate)) selection ``` It can, on data this size, while only one assumption fails at a time. The generating model has the lowest AIC in the large majority of replicates in every scenario, though a wrong pick can be costly: the few trap shy studies that AIC sends to Mt get an estimate about as inflated as M0's, because Mt has no term for the trap response. Two hundred animals over six occasions with a mean capture probability of `r sprintf("%.2f", p_flat)` is a generous study, and for a single violation the message is a positive one: the diagnosis is available in the data, and skipping the comparison is a choice rather than a limitation. The next subsection shows where that stops. ### Two violations at once Every scenario above breaks one assumption. Real animals can be unequal and trap shy at the same time, and Mb has a property that decides what happens then: its abundance estimate uses only the first capture of each animal. The recaptures inform the recapture probability and nothing else, so for abundance Mb is a removal estimator, and removal under individual heterogeneity runs low for the reason worked out in [Declaring eradication after empty traps](../declaring-eradication-after-negative-surveys/): the animals still uncaught after several occasions are the hard ones. The generator below draws each animal's first capture probability from the same Beta(1.2, 2.8) as the heterogeneous scenario and multiplies it by a fixed factor once the animal has been caught: `r sprintf("%.2f", p_again / p_first)`, the ratio of the trap shy scenario, for shy animals and 1.5 for trap happy ones. The last lines fit the removal likelihood to the expected first capture counts, `N B(a + 1, b + j - 1) / B(a, b)` on occasion `j`; the result is the answer Mb gives when the first captures of a `r sprintf("%d", n_true)` animal study fall exactly on their expected values. ```{r two-violations} #| warning: false gen_mbh <- function(N, a, b, mult) { pi_ind <- rbeta(N, a, b); recap <- pmin(1, mult * pi_ind) X <- matrix(0, N, n_occ); caught <- rep(FALSE, N) for (j in seq_len(n_occ)) { X[, j] <- rbinom(N, 1, ifelse(caught, recap, pi_ind)) caught <- caught | X[, j] == 1 } X } fit_cover <- function(X) { seen <- X[rowSums(X) > 0, , drop = FALSE]; n_seen <- nrow(seen); x <- rowSums(seen) f0 <- optim(c(log(5), 0), nll_m0, x = x, n_seen = n_seen) se <- sqrt(solve(optimHess(f0$par, nll_m0, x = x, n_seen = n_seen))[1, 1]) c(fit_all(X), m0_covers = abs(log(n_true - n_seen) - f0$par[1]) < 1.96 * se) } picked_value <- function(r) { aic <- as.matrix(r[, c("aic_M0", "aic_Mt", "aic_Mb", "aic_Mh")]) picked <- c("M0", "Mt", "Mb", "Mh")[apply(aic, 1, which.min)] as.matrix(r)[cbind(seq_len(nrow(r)), match(picked, names(r)))] } set.seed(58) both <- list( "shy, heterogeneous" = as.data.frame(t(replicate(n_rep, fit_cover(gen_mbh(n_true, 1.2, 2.8, p_again / p_first))))), "happy, heterogeneous" = as.data.frame(t(replicate(n_rep, fit_cover(gen_mbh(n_true, 1.2, 2.8, 1.5)))))) selection_both <- do.call(rbind, lapply(both, pick_rate)) both_est <- sapply(both, function(r) c(picked = median(picked_value(r)), M0 = median(r$M0), Mb = median(r$Mb), Mh = median(r$Mh), M0_cover = mean(r$m0_covers))) first_cap <- n_true * exp(lbeta(1.2 + 1, 2.8 + seq_len(n_occ) - 1) - lbeta(1.2, 2.8)) nll_removal <- function(par, u) { n_u <- sum(u); N <- n_u + exp(par[1]); p <- plogis(par[2]) -(lgamma(N + 1) - lgamma(N - n_u + 1) + n_u * log(p) + (sum((seq_len(n_occ) - 1) * u) + (N - n_u) * n_occ) * log(1 - p)) } fit_limit <- optim(c(log(20), 0), nll_removal, u = first_cap) mb_limit <- sum(first_cap) + exp(fit_limit$par[1]) selection_both round(both_est, 2) round(mb_limit, 1) ``` With shy animals that also differ, AIC picks Mb in `r sprintf("%.0f", selection_both["shy, heterogeneous", "Mb"])` per cent of `r sprintf("%d", n_rep)` studies, and the estimate of whichever model AIC picks has a median of `r sprintf("%.0f", both_est["picked", "shy, heterogeneous"])` against a true `r sprintf("%d", n_true)`. That shortfall owes nothing to the trap response. The removal fit to expected first captures gives `r sprintf("%.0f", mb_limit)`, and no recapture factor appears anywhere in it; Mb on the trap happy data has a median of `r sprintf("%.0f", both_est["Mb", "happy, heterogeneous"])`, the same number. The comparison picks the model that explains the recaptures and inherits the heterogeneity bias through the first captures. AIC picks M0 in `r if (selection_both["shy, heterogeneous", "M0"] == 0) "none" else sprintf("%.0f per cent", selection_both["shy, heterogeneous", "M0"])` of these studies. Its median is `r sprintf("%.0f", both_est["M0", "shy, heterogeneous"])`, and its Wald interval covers the true abundance in `r sprintf("%.0f", 100 * both_est["M0_cover", "shy, heterogeneous"])` per cent of them. At these two strengths the upward push of shyness and the downward pull of heterogeneity cancel. Nothing in the histories tells the analyst that they have, and the model the comparison throws out is the one that happens to be right. With trap happy animals the extra recaptures look like heterogeneity: AIC picks Mh in `r sprintf("%.0f", selection_both["happy, heterogeneous", "Mh"])` per cent of studies, and the picked estimate has a median of `r sprintf("%.0f", both_est["picked", "happy, heterogeneous"])`. A model with both effects, Mbh, belongs to the mixture family of Pledger (2000); it is not in the candidate set here, and this post does not measure how it does. ## How much heterogeneity is too much The Mh scenario above used one particular amount of individual variation. Holding the mean capture probability fixed and varying its spread shows how quickly the constant probability estimator gives way. ```{r sweep} #| warning: false nll_two_pt <- function(par, x, n_seen) { N <- n_seen + exp(par[1]); w <- plogis(par[2]) p1 <- plogis(par[3]); p2 <- plogis(par[4]) hist_p <- function(k) log(w * p1^k * (1 - p1)^(n_occ - k) + (1 - w) * p2^k * (1 - p2)^(n_occ - k)) -(lgamma(N + 1) - lgamma(N - n_seen + 1) + sum(hist_p(x)) + (N - n_seen) * hist_p(0)) } fit_hetero <- function(X) { s <- X[rowSums(X) > 0, , drop = FALSE]; n_seen <- nrow(s); x <- rowSums(s) f0 <- optim(c(log(5), 0), nll_m0, x = x, n_seen = n_seen) fh <- optim(c(log(5), 0, 0), nll_mh, x = x, n_seen = n_seen, control = list(maxit = 3000)) f2 <- optim(c(log(5), 0, -1, 1), nll_two_pt, x = x, n_seen = n_seen, control = list(maxit = 5000)) c(M0 = n_seen + exp(f0$par[1]), beta = n_seen + exp(fh$par[1]), two_pt = n_seen + exp(f2$par[1]), jackknife = n_seen + sum(x == 1) * (n_occ - 1) / n_occ, nll_beta = fh$value, nll_two = f2$value) } set.seed(21) mean_p <- p_flat shapes <- c(50, 6, 2.4, 1.2, 0.6) sweep_h <- do.call(rbind, lapply(shapes, function(a) { b <- a * (1 - mean_p) / mean_p cv <- sqrt((a * b) / ((a + b)^2 * (a + b + 1))) / mean_p r <- as.data.frame(t(replicate(50, fit_hetero(gen_mh(n_true, a, b))))) data.frame(cv = cv, M0 = median(r$M0), beta = median(r$beta), two_pt = median(r$two_pt), jackknife = median(r$jackknife)) })) print(round(sweep_h, 1)) ``` At a coefficient of variation of `r sprintf("%.2f", sweep_h$cv[2])` in capture probability, which is modest for a real population, the constant probability estimate is already `r sprintf("%.0f", 100 * (1 - sweep_h$M0[2] / n_true))` per cent low. At `r sprintf("%.2f", sweep_h$cv[5])` it is `r sprintf("%.0f", 100 * (1 - sweep_h$M0[5] / n_true))` per cent low. The jackknife is the standard insurance against this, and the first row is its premium. On data with essentially no heterogeneity it returns `r sprintf("%.0f", sweep_h$jackknife[1])` against a true `r sprintf("%d", n_true)`, an overestimate of `r sprintf("%.0f", 100 * (sweep_h$jackknife[1] / n_true - 1))` per cent. Insurance that costs a `r sprintf("%.0f", 100 * (sweep_h$jackknife[1] / n_true - 1))` per cent bias when the risk is absent is worth buying only when the risk is real, which is an argument for looking rather than for defaulting either way. ```{r fig-sweep} #| fig-cap: "Median estimated abundance against the coefficient of variation of individual capture probability, holding the mean at 0.30, over 50 replicates per point." #| fig-alt: "Four lines against increasing variation in capture probability. The constant probability line falls steadily from 200 to about 135. The jackknife starts at 225 and falls to about 164. The beta mixture stays near the dashed line at 200." #| fig-width: 6.8 #| fig-height: 4.2 d_sw <- do.call(rbind, lapply(c("M0", "beta", "two_pt", "jackknife"), function(k) data.frame(cv = sweep_h$cv, value = sweep_h[[k]], estimator = k))) d_sw$estimator <- factor(d_sw$estimator, levels = c("M0", "beta", "two_pt", "jackknife"), labels = c("M0, constant p", "Mh, beta mixture", "Mh, two point mixture", "jackknife")) ggplot(d_sw, aes(cv, value, colour = estimator)) + geom_hline(yintercept = n_true, colour = te_ink, linetype = "dashed", linewidth = 0.6) + geom_line(linewidth = 1) + geom_point(size = 2.2) + scale_colour_manual(values = c(te_rust, te_forest, te_gold, "#7b6a4f"), name = NULL) + labs(x = "coefficient of variation of individual capture probability", y = "median estimated abundance", title = "Where the constant probability estimate goes") + theme_datasheet() + theme(legend.position = "bottom", legend.text = element_text(size = 9)) ``` ## Two heterogeneity models, one dataset, two answers The `two_pt` column of that sweep is a second Mh model: instead of a beta distribution of capture probabilities it assumes two kinds of animal, easy and hard, in unknown proportions. Both are respectable choices, and on the same data they do not agree. ```{r identifiability} #| warning: false set.seed(404) pair <- as.data.frame(t(replicate(40, fit_hetero(gen_mh(n_true, 1.2, 2.8))))) gap_n <- median(abs(pair$beta - pair$two_pt)) gap_nll <- median(abs(pair$nll_beta - pair$nll_two)) ``` Over `r sprintf("%d", nrow(pair))` datasets the two mixing distributions differ by a median of `r sprintf("%.0f", gap_n)` animals, which is `r sprintf("%.0f", 100 * gap_n / n_true)` per cent of the true abundance, while their log likelihoods differ by a median of `r sprintf("%.2f", gap_nll)`. The data cannot separate them and the answers are not the same. Link proved in 2003 that this is not a numerical accident: different mixing distributions can give identical probabilities for every observable capture frequency and still imply different numbers of animals never caught. The unseen part of the population is where the disagreement lives, and no amount of fieldwork on the seen part settles it. ```{r fig-ident} #| fig-cap: "Abundance from a beta mixture against abundance from a two point mixture, fitted to the same 40 simulated datasets. The dashed line marks equality and the solid lines mark the true abundance." #| fig-alt: "Scatter of forty points, most of them well below the dashed equality line, so the two point mixture returns a lower abundance than the beta mixture on nearly every dataset. Solid lines mark the true abundance of 200 on both axes: the cloud straddles the vertical line and sits mostly below the horizontal one." #| fig-width: 6.2 #| fig-height: 4.2 ggplot(pair, aes(beta, two_pt)) + geom_abline(slope = 1, intercept = 0, colour = te_ink, linetype = "dashed", linewidth = 0.6) + geom_hline(yintercept = n_true, colour = te_gold, linewidth = 0.7) + geom_vline(xintercept = n_true, colour = te_gold, linewidth = 0.7) + geom_point(colour = te_forest, size = 2.2, alpha = 0.85) + labs(x = "abundance from the beta mixture", y = "abundance from the two point mixture", title = "Same histories, two defensible answers") + theme_datasheet() ``` ## What to report Fit more than M0 and report the comparison. The measurement above says a single violation can be diagnosed in a study as well sampled as this one, so a paper with data of that quality that reports only M0 has declined to look rather than been unable to. The comparison is evidence about which assumption fails, not proof: when two fail together it can settle firmly on a model that is still wrong. Say which direction the residual risk runs. Trap shyness and individual heterogeneity bias abundance in opposite directions (trap happiness pushes the same way as heterogeneity), so a study that has ruled out one of them has not narrowed the interval on the other. Report the number of animals caught exactly once alongside the estimate. It is the statistic every heterogeneity correction is built from, and a reader can see from it how much of the answer is extrapolation. If a heterogeneity model is used, say which mixing distribution and acknowledge that the choice moves the answer. Reporting a beta mixture estimate to three significant figures without that sentence overstates what the data contain. ## Honest limits The population is genuinely closed here: no births, no deaths, no movement in or out. Real six occasion studies are closed by assumption and the assumption is usually a compromise about the length of the trapping session. Losses to the population and a trap response are also hard to tell apart from capture histories alone. Outside the one subsection on two violations at once, only one violation acts at a time, and that subsection tries a single strength of heterogeneity with two trap responses. The named combinations Mth, Mbh and Mtbh exist because real studies have several at once; none of them was fitted here. The behavioural response in the single violation scenarios is trap shyness, where recapture is less likely than first capture. Trap happiness is at least as common with baited traps and reverses the sign, so the direction reported for the trap shy scenario is the direction of that scenario, not a property of behavioural response as such. The two hundred animals and mean capture probability of `r sprintf("%.2f", p_flat)` make this a well sampled study. Halve the capture probability and the number of animals never seen grows, every estimator's variance grows with it, and the model selection result reported here would weaken, though that case was not run for this post. There is a reconciliation to make with the mark resight post on this site, where individual heterogeneity left the point estimate with only a small upward bias but destroyed the coverage of the interval. The difference is structural. In mark resight the number of marked animals is known, so heterogeneity affects how well the marked fraction is estimated. Here abundance is inferred from the shape of the capture frequency distribution, and heterogeneity changes that shape, so it moves the estimate itself. ## References Burnham KP, Overton WS 1978 Biometrika 65(3):625-633 (10.1093/biomet/65.3.625) Pledger S 2000 Biometrics 56(2):434-442 (10.1111/j.0006-341X.2000.00434.x) Link WA 2003 Biometrics 59(4):1123-1130 (10.1111/j.0006-341X.2003.00129.x) ## Related tutorials - [Closed-population capture-recapture in R](../closed-capture-recapture/) - [Disturbed animals in repeated counts](../disturbed-animals-in-repeated-counts/) - [Removal and depletion sampling in R](../removal-sampling-depletion/) - [Trap-specific responses in spatial capture-recapture](../trap-specific-responses-in-scr/)