The Royle-Nichols model: abundance from detections

occupancy

imperfect detection

abundance

ecology tutorial

How repeat visits let an occupancy model estimate abundance, why the trick works, and why unmodelled detection heterogeneity is read as abundance.

Author

Tidy Ecology

Published

2026-08-24

You visited 300 sites six times each and wrote down, for every visit, whether you saw the species. You never counted anything. The Royle-Nichols model claims it can tell you the mean number of individuals per site from those ticks alone.

That claim deserves suspicion, and this post spends most of its length earning it and then taking it back. The machinery works, the arithmetic is exact, and the estimate is nonetheless only as good as an assumption the data cannot check.

library(ggplot2)
library(grid)

te_ink <- "#16241d"; te_forest <- "#275139"; te_gold <- "#c9b458"
te_rust <- "#b5534e"; te_sage <- "#93a87f"

theme_te <- function() {
  theme_minimal(base_size = 11) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "grey90", linewidth = 0.3),
          plot.title = element_text(colour = te_ink, face = "bold", size = 11),
          plot.subtitle = element_text(colour = "grey35", size = 9),
          axis.title = element_text(colour = te_ink, size = 9),
          axis.text = element_text(colour = "grey30", size = 8),
          legend.position = "bottom",
          legend.title = element_blank(),
          legend.text = element_text(size = 8))
}
two_panel <- function(p1, p2) {
  grid.newpage()
  pushViewport(viewport(layout = grid.layout(1, 2)))
  print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
  print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
}

The idea in one line

A site with more individuals is easier to detect. If each of the $N_i$ individuals at site $i$ is seen independently on a visit with probability $r$, then the site-level detection probability is not a constant but

\[p_i = 1 - (1 - r)^{N_i}.\]

Sites with $N_i = 0$ have $p_i = 0$: they are unoccupied. Sites with $N_i = 5$ are far easier than sites with $N_i = 1$. Put a Poisson prior on abundance, $N_i \sim \text{Poisson}(\lambda)$, and marginalise it away:

\[\Pr(d_i \mid \lambda, r) = \sum_{N = 0}^{\infty} \frac{e^{-\lambda} \lambda^N}{N!} \binom{K}{d_i} \left[1 - (1-r)^N\right]^{d_i} \left[(1-r)^N\right]^{K - d_i},\]

where $d_i$ is the number of visits out of $K$ on which the species was recorded. Occupancy falls out as a derived quantity: a site is occupied when $N_i > 0$, so $\psi = 1 - e^{-\lambda}$.

That is the whole model (Royle and Nichols 2003). Two parameters, one infinite sum truncated in practice, and no marked individuals anywhere.

set.seed(4275)
R <- 300L; K <- 6L
lambda_true <- 2.0; r_true <- 0.12
N <- rpois(R, lambda_true)
d <- rbinom(R, K, 1 - (1 - r_true)^N)

Nmax <- 80L; Nv <- 0:Nmax
rn_nll <- function(par, d, K) {
  lam <- exp(par[1]); r <- plogis(par[2])
  pd <- colSums(dpois(Nv, lam) *
    t(sapply(Nv, function(n) dbinom(0:K, K, 1 - (1 - r)^n))))
  -sum(log(pmax(pd[d + 1], 1e-300)))
}
occ_nll <- function(par, d, K) {
  psi <- plogis(par[1]); p <- plogis(par[2])
  pd <- psi * dbinom(0:K, K, p) + c(1 - psi, rep(0, K))
  -sum(log(pmax(pd[d + 1], 1e-300)))
}
ct <- list(maxit = 8000, reltol = 1e-12)
fit_nm <- function(nll, start, d, K) {
  o <- optim(start, nll, d = d, K = K, method = "Nelder-Mead", control = ct)
  optim(o$par, nll, d = d, K = K, method = "Nelder-Mead", control = ct)
}
f_rn <- fit_nm(rn_nll, c(0, 0), d, K)
f_oc <- fit_nm(occ_nll, c(0, 0), d, K)

se_rn <- sqrt(diag(solve(optimHess(f_rn$par, rn_nll, d = d, K = K))))
lam_hat <- exp(f_rn$par[1]); r_hat <- plogis(f_rn$par[2])
lam_ci <- exp(f_rn$par[1] + c(-1.96, 1.96) * se_rn[1])
psi_rn <- 1 - exp(-lam_hat)
psi_oc <- plogis(f_oc$par[1])
psi_true <- 1 - exp(-lambda_true)
naive <- mean(d > 0)

The truth here is $\lambda = 2.0$ individuals per site and $r = 0.12$ per individual per visit, which puts true occupancy at $\psi = 0.8647$. The species was recorded at least once at 66.3% of sites, so the naive estimate is already off by 0.20.

The fitted model returns $\hat\lambda = 1.906$ with a 95% interval of [1.371, 2.648] and $\hat{r} = 0.1333$, giving $\hat\psi = 0.8513$ against a truth of 0.8647. A standard occupancy model on the same ticks reports $\hat\psi = 0.7527$, which is 0.112 too low: with a single detection probability shared by all sites, the model has no way to know that the silent sites are silent because they are sparse rather than because they are empty.

A note on the optimiser, which matters more than it should. Fitting this likelihood with BFGS on the logit and log scales converges happily, reports success, and lands on a worse answer with a singular Hessian. Nelder-Mead restarted from its own solution, with optimHess for the standard errors, is slower and correct. The #| label on that chunk is doing nothing exciting; the restart is.

Where the abundance actually comes from

Nothing in the data says “abundance”. The only thing the model sees is how often each site produced a detection. So which feature of those numbers carries the information?

Start with an exact identity. For $N \sim \text{Poisson}(\lambda)$,

\[\mathbb{E}\left[(1-r)^N\right] = \sum_{N=0}^{\infty} \frac{e^{-\lambda}\lambda^N}{N!}(1-r)^N = e^{-\lambda r},\]

which is the probability of a non-detection on a single visit. With $K = 1$ the likelihood depends on $\lambda$ and $r$ only through their product. Every pair with $\lambda r = 0.24$ fits identically. The parameters are not merely imprecise, they are unidentified: the likelihood surface is an exact ridge.

prod_target <- lambda_true * r_true
lam_grid <- exp(seq(log(0.4), log(9), length.out = 60))
ridge <- data.frame(lambda = lam_grid, r = prod_target / lam_grid)
ridge <- ridge[ridge$r < 0.95, ]
nll_at <- function(lam, r, K, dd) rn_nll(c(log(lam), qlogis(r)), dd, K)

d1 <- rbinom(R, 1L, 1 - (1 - r_true)^N)
ridge$nll_k1 <- mapply(nll_at, ridge$lambda, ridge$r, MoreArgs = list(K = 1L, dd = d1))
ridge$nll_k6 <- mapply(nll_at, ridge$lambda, ridge$r, MoreArgs = list(K = K, dd = d))
spread_k1 <- diff(range(ridge$nll_k1))
spread_k6 <- diff(range(ridge$nll_k6))

Along that ridge the one-visit negative log-likelihood varies by 1.705e-13, which is zero to machine precision. The six-visit likelihood varies by 288.86 over the same curve and has a clear minimum.

pair <- data.frame(lambda = c(1, 4), r = c(0.24, 0.06))
pair$empty <- dpois(0, pair$lambda)
pair$p_typ <- 1 - (1 - pair$r)^round(pair$lambda)
pair$p_one <- 1 - exp(-pair$lambda * pair$r)

The repeat visits are what break the tie, and it is worth seeing why. Take two points on the ridge. With $\lambda = 1$ and $r = 0.24$, 36.8% of sites hold nothing and are silent forever, while a typical occupied site is detected about 24% of the time. With $\lambda = 4$ and $r = 0.06$, only 1.8% of sites are empty, and a typical site is detected about 22% of the time. Both worlds produce a detection on 21.3% of single visits, which is why one visit cannot separate them. Six visits can, because the first world leaves a large pile of all-zero sites and the second leaves almost none. Abundance is inferred from the shape of the detection-frequency table, not from any count.

obs <- as.numeric(table(factor(d, levels = 0:K)))
exp_occ <- R * (plogis(f_oc$par[1]) * dbinom(0:K, K, plogis(f_oc$par[2])) +
                c(1 - plogis(f_oc$par[1]), rep(0, K)))
exp_rn <- R * colSums(dpois(Nv, lam_hat) *
  t(sapply(Nv, function(n) dbinom(0:K, K, 1 - (1 - r_hat)^n))))

df1 <- data.frame(d = rep(0:K, 2), count = c(exp_occ, exp_rn),
                  model = rep(c("standard occupancy", "Royle-Nichols"), each = K + 1))
p1 <- ggplot() +
  geom_col(data = data.frame(d = 0:K, count = obs), aes(d, count),
           fill = te_sage, alpha = 0.55, width = 0.75) +
  geom_line(data = df1, aes(d, count, colour = model), linewidth = 0.7) +
  geom_point(data = df1, aes(d, count, colour = model), size = 1.6) +
  scale_colour_manual(values = c("standard occupancy" = te_rust,
                                 "Royle-Nichols" = te_forest)) +
  scale_x_continuous(breaks = 0:K) +
  labs(title = "Both models fit the frequencies", x = "detections out of six",
       y = "sites", subtitle = "bars: observed") +
  theme_te()

df2 <- data.frame(lambda = rep(ridge$lambda, 2),
                  nll = c(ridge$nll_k1 - min(ridge$nll_k1),
                          ridge$nll_k6 - min(ridge$nll_k6)),
                  visits = rep(c("one visit", "six visits"), each = nrow(ridge)))
p2 <- ggplot(df2, aes(lambda, nll, colour = visits)) +
  geom_line(linewidth = 0.8) +
  geom_vline(xintercept = lambda_true, linetype = "dashed",
             colour = te_ink, linewidth = 0.3) +
  scale_colour_manual(values = c("one visit" = te_gold, "six visits" = te_forest)) +
  scale_x_log10() +
  labs(title = "One visit cannot separate them", x = "abundance (log scale)",
       y = "negative log-likelihood, relative",
       subtitle = "along the curve where abundance times detection is fixed") +
  theme_te()
two_panel(p1, p2)

Two panels. The left shows bars of observed site counts by number of detections with two fitted curves tracking them. The right shows a flat line for one visit and a U-shaped curve with a minimum for six visits. — Figure 1: Left: observed detection frequencies with the expected frequencies from the standard occupancy model and the Royle-Nichols model. Right: the negative log-likelihood along the curve where the product of abundance and per-individual detection is held constant, for one visit and for six.

The part that should worry you

Here is the problem. Abundance is not the only thing that makes some sites easier than others. Wind, canopy, observer skill, time of day, and the animal’s own behaviour all produce sites where detection is high and sites where it is low. The Royle-Nichols model has one slot for site-level variation in detection, and that slot is labelled abundance. Anything that lands in it comes out wearing that label.

So simulate a world with no abundance at all. Every occupied site holds the species and nothing else is known about it; detection probability varies across occupied sites according to a Beta distribution. There is no $N$, no $\lambda$, no counting to be done.

set.seed(4273)
z_het <- rbinom(R, 1L, 0.60)
p_het <- rbeta(R, 1.2, 4.0)
d_het <- rbinom(R, K, z_het * p_het)

mix_nll <- function(par, d, K) {
  psi <- plogis(par[1]); p1 <- plogis(par[2]); p2 <- plogis(par[3]); w <- plogis(par[4])
  pd <- psi * (w * dbinom(0:K, K, p1) + (1 - w) * dbinom(0:K, K, p2)) +
        c(1 - psi, rep(0, K))
  -sum(log(pmax(pd[d + 1], 1e-300)))
}
aic <- function(f, k) 2 * f$value + 2 * k
h_oc <- fit_nm(occ_nll, c(0, 0), d_het, K)
h_rn <- fit_nm(rn_nll, c(0, 0), d_het, K)
h_mx <- fit_nm(mix_nll, c(0, qlogis(0.4), qlogis(0.1), 0), d_het, K)
lam_het <- exp(h_rn$par[1])
se_het <- sqrt(diag(solve(optimHess(h_rn$par, rn_nll, d = d_het, K = K))))
lam_het_ci <- exp(h_rn$par[1] + c(-1.96, 1.96) * se_het[1])

Fitted to data containing no individuals to count, the Royle-Nichols model reports $\hat\lambda = 0.705$ individuals per site, with a 95% interval of [0.563, 0.882] that excludes nothing of interest and certainly does not warn you. Worse, it wins: AIC 768.5 against 779.4 for the standard occupancy model. Model selection rewards the model that invented a quantity out of nothing, because that quantity is a perfectly good description of heterogeneity, which is really there. Only its name is wrong.

The mirror image is worse still. Return to the genuine Royle-Nichols data, where abundance really does drive detection, and fit a two-point finite mixture instead: occupied sites are either easy or hard, in some proportion, with no abundance anywhere in the specification.

g_mx <- fit_nm(mix_nll, c(0, qlogis(0.4), qlogis(0.1), 0), d, K)
psi_mx <- plogis(g_mx$par[1])
aic_rn <- aic(f_rn, 2); aic_mx <- aic(g_mx, 4)

prof_psi <- function(psi_fix, d, K) {
  f <- function(q) {
    p1 <- plogis(q[1]); p2 <- plogis(q[2]); w <- plogis(q[3])
    pd <- psi_fix * (w * dbinom(0:K, K, p1) + (1 - w) * dbinom(0:K, K, p2)) +
          c(1 - psi_fix, rep(0, K))
    -sum(log(pmax(pd[d + 1], 1e-300)))
  }
  o <- optim(c(qlogis(0.4), qlogis(0.1), 0), f, method = "Nelder-Mead", control = ct)
  optim(o$par, f, method = "Nelder-Mead", control = ct)$value
}
prof_rn <- function(psi_fix, d, K) {
  lam_fix <- -log(1 - psi_fix)
  f <- function(q) rn_nll(c(log(lam_fix), q), d, K)
  optimise(f, c(-8, 4))$objective
}
psi_grid <- seq(0.70, 0.999, length.out = 26)
dev_mx <- 2 * (sapply(psi_grid, prof_psi, d = d, K = K) - g_mx$value)
dev_rn <- 2 * (sapply(psi_grid, prof_rn, d = d, K = K) - f_rn$value)
dev_mx_at_truth <- 2 * (prof_psi(psi_true, d, K) - g_mx$value)

On the same detection histories the mixture reports $\hat\psi = 1.000$, meaning every single site is occupied, against the Royle-Nichols answer of 0.851 and a truth of 0.865. The AIC gap between them is 0.1, which is nothing. And the mixture’s answer is not an optimiser failure: its profile deviance is monotone up to the boundary, sitting at 0.01 when $\psi$ is pinned at the true value, well inside the 3.84 cutoff. The likelihood genuinely prefers “everything is occupied, some sites are just very hard” over the truth.

obs_h <- as.numeric(table(factor(d_het, levels = 0:K)))
ex <- function(par, kind) {
  if (kind == "occ") R * (plogis(par[1]) * dbinom(0:K, K, plogis(par[2])) +
                          c(1 - plogis(par[1]), rep(0, K)))
  else if (kind == "rn") R * colSums(dpois(Nv, exp(par[1])) *
      t(sapply(Nv, function(n) dbinom(0:K, K, 1 - (1 - plogis(par[2]))^n))))
  else R * (plogis(par[1]) * (plogis(par[4]) * dbinom(0:K, K, plogis(par[2])) +
      (1 - plogis(par[4])) * dbinom(0:K, K, plogis(par[3]))) +
      c(1 - plogis(par[1]), rep(0, K)))
}
df3 <- data.frame(d = rep(0:K, 3),
                  count = c(ex(h_oc$par, "occ"), ex(h_rn$par, "rn"), ex(h_mx$par, "mix")),
                  model = rep(c("standard occupancy", "Royle-Nichols", "two-point mixture"),
                              each = K + 1))
p3 <- ggplot() +
  geom_col(data = data.frame(d = 0:K, count = obs_h), aes(d, count),
           fill = te_sage, alpha = 0.55, width = 0.75) +
  geom_line(data = df3, aes(d, count, colour = model), linewidth = 0.7) +
  scale_colour_manual(values = c("standard occupancy" = te_rust,
                                 "Royle-Nichols" = te_forest,
                                 "two-point mixture" = te_gold)) +
  scale_x_continuous(breaks = 0:K) +
  labs(title = "No abundance in these data", x = "detections out of six", y = "sites",
       subtitle = "Royle-Nichols still reports individuals per site") +
  theme_te()

df4 <- data.frame(psi = rep(psi_grid, 2), dev = c(dev_rn, dev_mx),
                  model = rep(c("Royle-Nichols", "two-point mixture"),
                              each = length(psi_grid)))
p4 <- ggplot(df4, aes(psi, dev, colour = model)) +
  geom_line(linewidth = 0.8) +
  geom_hline(yintercept = qchisq(0.95, 1), linetype = "dotted", colour = "grey40") +
  geom_vline(xintercept = psi_true, linetype = "dashed", colour = te_ink,
             linewidth = 0.3) +
  scale_colour_manual(values = c("Royle-Nichols" = te_forest,
                                 "two-point mixture" = te_gold)) +
  coord_cartesian(ylim = c(0, 12)) +
  labs(title = "Same data, two answers", x = "occupancy", y = "profile deviance",
       subtitle = "dashed: truth. dotted: the 3.84 cutoff") +
  theme_te()
two_panel(p3, p4)

Two panels. The left shows observed bars with three near-identical fitted curves. The right shows a U-shaped deviance curve with a clear minimum and a monotone curve running to the right-hand boundary. — Figure 2: Left: expected detection frequencies for data with heterogeneous detection but no abundance, under three models. Right: profile deviance for occupancy on genuine Royle-Nichols data, under the Royle-Nichols model and under a two-point mixture.

What this means when you have real ticks

The two failures are the same failure seen from opposite sides. Site-level variation in detection is real and estimable; its cause is not. Royle (2006) makes the point directly for occupancy models, and Link (2003) proved the general version for closed-population abundance: different mixing distributions can produce near-identical data and quite different inferences, and no quantity of data resolves the choice, because the choice is not in the data.

So the honest reading of a Royle-Nichols fit is narrow. The $\hat\lambda$ is an estimate of mean abundance if you are willing to assert that abundance is the only source of site-level heterogeneity in detection, that individuals are detected independently, and that the abundance distribution is Poisson. The first of those three is an ecological claim about your field protocol, not a statistical one, and it is exactly the claim that fails when you survey along a gradient of vegetation density.

Model selection does not rescue you. On the very data set generated above, the standard occupancy model has an AIC of 908.6 against the Royle-Nichols model’s 912.3, so the AIC prefers the model that is wrong about the process and wrong about occupancy by 0.112. What the fit statistics can and cannot decide here is the subject of checking an occupancy variant, and it is not a comfortable read.

What does help is design. Count something, even crudely, at a subset of sites, and abundance becomes identified by data rather than by assumption. Measure the covariates you believe drive detection and put them in the detection model, so that the heterogeneity slot is not the only place left for them to go. Neither of these is an analysis decision.

References

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

MacKenzie D, Nichols J, Lachman G, Droege S, Royle J, Langtimm C 2002 Ecology 83(8):2248-2255

Royle J 2006 Biometrics 62(1):97-102 (doi:10.1111/j.1541-0420.2005.00439.x)

Royle J, Nichols J 2003 Ecology 84(3):777-790

Kery M, Royle J 2016 Applied Hierarchical Modeling in Ecology. Academic Press (ISBN 978-0-12-801378-6)

--- title: "The Royle-Nichols model: abundance from detections" description: "How repeat visits let an occupancy model estimate abundance, why the trick works, and why unmodelled detection heterogeneity is read as abundance." date: "2026-08-24 09:00" categories: [R, occupancy, imperfect detection, abundance, ecology tutorial] image: thumbnail.png image-alt: "Two panels: fitted detection frequencies against observed counts, and a likelihood ridge that flattens with a single visit." --- You visited 300 sites six times each and wrote down, for every visit, whether you saw the species. You never counted anything. The Royle-Nichols model claims it can tell you the mean number of individuals per site from those ticks alone. That claim deserves suspicion, and this post spends most of its length earning it and then taking it back. The machinery works, the arithmetic is exact, and the estimate is nonetheless only as good as an assumption the data cannot check. ```{r} #| label: setup #| message: false library(ggplot2) library(grid) te_ink <- "#16241d"; te_forest <- "#275139"; te_gold <- "#c9b458" te_rust <- "#b5534e"; te_sage <- "#93a87f" theme_te <- function() { theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "grey90", linewidth = 0.3), plot.title = element_text(colour = te_ink, face = "bold", size = 11), plot.subtitle = element_text(colour = "grey35", size = 9), axis.title = element_text(colour = te_ink, size = 9), axis.text = element_text(colour = "grey30", size = 8), legend.position = "bottom", legend.title = element_blank(), legend.text = element_text(size = 8)) } two_panel <- function(p1, p2) { grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) } ``` ## The idea in one line A site with more individuals is easier to detect. If each of the $N_i$ individuals at site $i$ is seen independently on a visit with probability $r$, then the site-level detection probability is not a constant but $$p_i = 1 - (1 - r)^{N_i}.$$ Sites with $N_i = 0$ have $p_i = 0$: they are unoccupied. Sites with $N_i = 5$ are far easier than sites with $N_i = 1$. Put a Poisson prior on abundance, $N_i \sim \text{Poisson}(\lambda)$, and marginalise it away: $$\Pr(d_i \mid \lambda, r) = \sum_{N = 0}^{\infty} \frac{e^{-\lambda} \lambda^N}{N!} \binom{K}{d_i} \left[1 - (1-r)^N\right]^{d_i} \left[(1-r)^N\right]^{K - d_i},$$ where $d_i$ is the number of visits out of $K$ on which the species was recorded. Occupancy falls out as a derived quantity: a site is occupied when $N_i > 0$, so $\psi = 1 - e^{-\lambda}$. That is the whole model (Royle and Nichols 2003). Two parameters, one infinite sum truncated in practice, and no marked individuals anywhere. ```{r} #| label: model set.seed(4275) R <- 300L; K <- 6L lambda_true <- 2.0; r_true <- 0.12 N <- rpois(R, lambda_true) d <- rbinom(R, K, 1 - (1 - r_true)^N) Nmax <- 80L; Nv <- 0:Nmax rn_nll <- function(par, d, K) { lam <- exp(par[1]); r <- plogis(par[2]) pd <- colSums(dpois(Nv, lam) * t(sapply(Nv, function(n) dbinom(0:K, K, 1 - (1 - r)^n)))) -sum(log(pmax(pd[d + 1], 1e-300))) } occ_nll <- function(par, d, K) { psi <- plogis(par[1]); p <- plogis(par[2]) pd <- psi * dbinom(0:K, K, p) + c(1 - psi, rep(0, K)) -sum(log(pmax(pd[d + 1], 1e-300))) } ct <- list(maxit = 8000, reltol = 1e-12) fit_nm <- function(nll, start, d, K) { o <- optim(start, nll, d = d, K = K, method = "Nelder-Mead", control = ct) optim(o$par, nll, d = d, K = K, method = "Nelder-Mead", control = ct) } f_rn <- fit_nm(rn_nll, c(0, 0), d, K) f_oc <- fit_nm(occ_nll, c(0, 0), d, K) se_rn <- sqrt(diag(solve(optimHess(f_rn$par, rn_nll, d = d, K = K)))) lam_hat <- exp(f_rn$par[1]); r_hat <- plogis(f_rn$par[2]) lam_ci <- exp(f_rn$par[1] + c(-1.96, 1.96) * se_rn[1]) psi_rn <- 1 - exp(-lam_hat) psi_oc <- plogis(f_oc$par[1]) psi_true <- 1 - exp(-lambda_true) naive <- mean(d > 0) ``` The truth here is $\lambda = `r sprintf("%.1f", lambda_true)`$ individuals per site and $r = `r sprintf("%.2f", r_true)`$ per individual per visit, which puts true occupancy at $\psi = `r sprintf("%.4f", psi_true)`$. The species was recorded at least once at `r sprintf("%.1f%%", 100 * naive)` of sites, so the naive estimate is already off by `r sprintf("%.2f", psi_true - naive)`. The fitted model returns $\hat\lambda = `r sprintf("%.3f", lam_hat)`$ with a 95% interval of [`r sprintf("%.3f", lam_ci[1])`, `r sprintf("%.3f", lam_ci[2])`] and $\hat{r} = `r sprintf("%.4f", r_hat)`$, giving $\hat\psi = `r sprintf("%.4f", psi_rn)`$ against a truth of `r sprintf("%.4f", psi_true)`. A standard occupancy model on the same ticks reports $\hat\psi = `r sprintf("%.4f", psi_oc)`$, which is `r sprintf("%.3f", psi_true - psi_oc)` too low: with a single detection probability shared by all sites, the model has no way to know that the silent sites are silent because they are sparse rather than because they are empty. A note on the optimiser, which matters more than it should. Fitting this likelihood with BFGS on the logit and log scales converges happily, reports success, and lands on a worse answer with a singular Hessian. Nelder-Mead restarted from its own solution, with `optimHess` for the standard errors, is slower and correct. The `#| label` on that chunk is doing nothing exciting; the restart is. ## Where the abundance actually comes from Nothing in the data says "abundance". The only thing the model sees is how often each site produced a detection. So which feature of those numbers carries the information? Start with an exact identity. For $N \sim \text{Poisson}(\lambda)$, $$\mathbb{E}\left[(1-r)^N\right] = \sum_{N=0}^{\infty} \frac{e^{-\lambda}\lambda^N}{N!}(1-r)^N = e^{-\lambda r},$$ which is the probability of a non-detection on a single visit. With $K = 1$ the likelihood depends on $\lambda$ and $r$ **only through their product**. Every pair with $\lambda r = 0.24$ fits identically. The parameters are not merely imprecise, they are unidentified: the likelihood surface is an exact ridge. ```{r} #| label: ridge prod_target <- lambda_true * r_true lam_grid <- exp(seq(log(0.4), log(9), length.out = 60)) ridge <- data.frame(lambda = lam_grid, r = prod_target / lam_grid) ridge <- ridge[ridge$r < 0.95, ] nll_at <- function(lam, r, K, dd) rn_nll(c(log(lam), qlogis(r)), dd, K) d1 <- rbinom(R, 1L, 1 - (1 - r_true)^N) ridge$nll_k1 <- mapply(nll_at, ridge$lambda, ridge$r, MoreArgs = list(K = 1L, dd = d1)) ridge$nll_k6 <- mapply(nll_at, ridge$lambda, ridge$r, MoreArgs = list(K = K, dd = d)) spread_k1 <- diff(range(ridge$nll_k1)) spread_k6 <- diff(range(ridge$nll_k6)) ``` Along that ridge the one-visit negative log-likelihood varies by `r sprintf("%.3e", spread_k1)`, which is zero to machine precision. The six-visit likelihood varies by `r sprintf("%.2f", spread_k6)` over the same curve and has a clear minimum. ```{r} #| label: why pair <- data.frame(lambda = c(1, 4), r = c(0.24, 0.06)) pair$empty <- dpois(0, pair$lambda) pair$p_typ <- 1 - (1 - pair$r)^round(pair$lambda) pair$p_one <- 1 - exp(-pair$lambda * pair$r) ``` The repeat visits are what break the tie, and it is worth seeing why. Take two points on the ridge. With $\lambda = 1$ and $r = 0.24$, `r sprintf("%.1f%%", 100 * pair$empty[1])` of sites hold nothing and are silent forever, while a typical occupied site is detected about `r sprintf("%.0f%%", 100 * pair$p_typ[1])` of the time. With $\lambda = 4$ and $r = 0.06$, only `r sprintf("%.1f%%", 100 * pair$empty[2])` of sites are empty, and a typical site is detected about `r sprintf("%.0f%%", 100 * pair$p_typ[2])` of the time. Both worlds produce a detection on `r sprintf("%.1f%%", 100 * pair$p_one[1])` of single visits, which is why one visit cannot separate them. Six visits can, because the first world leaves a large pile of all-zero sites and the second leaves almost none. Abundance is inferred from the **shape of the detection-frequency table**, not from any count. ```{r} #| label: fig-ridge #| fig-width: 8 #| fig-height: 3.4 #| fig-cap: "Left: observed detection frequencies with the expected frequencies from the standard occupancy model and the Royle-Nichols model. Right: the negative log-likelihood along the curve where the product of abundance and per-individual detection is held constant, for one visit and for six." #| fig-alt: "Two panels. The left shows bars of observed site counts by number of detections with two fitted curves tracking them. The right shows a flat line for one visit and a U-shaped curve with a minimum for six visits." obs <- as.numeric(table(factor(d, levels = 0:K))) exp_occ <- R * (plogis(f_oc$par[1]) * dbinom(0:K, K, plogis(f_oc$par[2])) + c(1 - plogis(f_oc$par[1]), rep(0, K))) exp_rn <- R * colSums(dpois(Nv, lam_hat) * t(sapply(Nv, function(n) dbinom(0:K, K, 1 - (1 - r_hat)^n)))) df1 <- data.frame(d = rep(0:K, 2), count = c(exp_occ, exp_rn), model = rep(c("standard occupancy", "Royle-Nichols"), each = K + 1)) p1 <- ggplot() + geom_col(data = data.frame(d = 0:K, count = obs), aes(d, count), fill = te_sage, alpha = 0.55, width = 0.75) + geom_line(data = df1, aes(d, count, colour = model), linewidth = 0.7) + geom_point(data = df1, aes(d, count, colour = model), size = 1.6) + scale_colour_manual(values = c("standard occupancy" = te_rust, "Royle-Nichols" = te_forest)) + scale_x_continuous(breaks = 0:K) + labs(title = "Both models fit the frequencies", x = "detections out of six", y = "sites", subtitle = "bars: observed") + theme_te() df2 <- data.frame(lambda = rep(ridge$lambda, 2), nll = c(ridge$nll_k1 - min(ridge$nll_k1), ridge$nll_k6 - min(ridge$nll_k6)), visits = rep(c("one visit", "six visits"), each = nrow(ridge))) p2 <- ggplot(df2, aes(lambda, nll, colour = visits)) + geom_line(linewidth = 0.8) + geom_vline(xintercept = lambda_true, linetype = "dashed", colour = te_ink, linewidth = 0.3) + scale_colour_manual(values = c("one visit" = te_gold, "six visits" = te_forest)) + scale_x_log10() + labs(title = "One visit cannot separate them", x = "abundance (log scale)", y = "negative log-likelihood, relative", subtitle = "along the curve where abundance times detection is fixed") + theme_te() two_panel(p1, p2) ``` ## The part that should worry you Here is the problem. Abundance is not the only thing that makes some sites easier than others. Wind, canopy, observer skill, time of day, and the animal's own behaviour all produce sites where detection is high and sites where it is low. The Royle-Nichols model has one slot for site-level variation in detection, and that slot is labelled abundance. Anything that lands in it comes out wearing that label. So simulate a world with no abundance at all. Every occupied site holds the species and nothing else is known about it; detection probability varies across occupied sites according to a Beta distribution. There is no $N$, no $\lambda$, no counting to be done. ```{r} #| label: het-data set.seed(4273) z_het <- rbinom(R, 1L, 0.60) p_het <- rbeta(R, 1.2, 4.0) d_het <- rbinom(R, K, z_het * p_het) mix_nll <- function(par, d, K) { psi <- plogis(par[1]); p1 <- plogis(par[2]); p2 <- plogis(par[3]); w <- plogis(par[4]) pd <- psi * (w * dbinom(0:K, K, p1) + (1 - w) * dbinom(0:K, K, p2)) + c(1 - psi, rep(0, K)) -sum(log(pmax(pd[d + 1], 1e-300))) } aic <- function(f, k) 2 * f$value + 2 * k h_oc <- fit_nm(occ_nll, c(0, 0), d_het, K) h_rn <- fit_nm(rn_nll, c(0, 0), d_het, K) h_mx <- fit_nm(mix_nll, c(0, qlogis(0.4), qlogis(0.1), 0), d_het, K) lam_het <- exp(h_rn$par[1]) se_het <- sqrt(diag(solve(optimHess(h_rn$par, rn_nll, d = d_het, K = K)))) lam_het_ci <- exp(h_rn$par[1] + c(-1.96, 1.96) * se_het[1]) ``` Fitted to data containing no individuals to count, the Royle-Nichols model reports $\hat\lambda = `r sprintf("%.3f", lam_het)`$ individuals per site, with a 95% interval of [`r sprintf("%.3f", lam_het_ci[1])`, `r sprintf("%.3f", lam_het_ci[2])`] that excludes nothing of interest and certainly does not warn you. Worse, it wins: AIC `r sprintf("%.1f", aic(h_rn, 2))` against `r sprintf("%.1f", aic(h_oc, 2))` for the standard occupancy model. Model selection rewards the model that invented a quantity out of nothing, because that quantity is a perfectly good description of heterogeneity, which is really there. Only its name is wrong. The mirror image is worse still. Return to the genuine Royle-Nichols data, where abundance really does drive detection, and fit a two-point finite mixture instead: occupied sites are either easy or hard, in some proportion, with no abundance anywhere in the specification. ```{r} #| label: mirror g_mx <- fit_nm(mix_nll, c(0, qlogis(0.4), qlogis(0.1), 0), d, K) psi_mx <- plogis(g_mx$par[1]) aic_rn <- aic(f_rn, 2); aic_mx <- aic(g_mx, 4) prof_psi <- function(psi_fix, d, K) { f <- function(q) { p1 <- plogis(q[1]); p2 <- plogis(q[2]); w <- plogis(q[3]) pd <- psi_fix * (w * dbinom(0:K, K, p1) + (1 - w) * dbinom(0:K, K, p2)) + c(1 - psi_fix, rep(0, K)) -sum(log(pmax(pd[d + 1], 1e-300))) } o <- optim(c(qlogis(0.4), qlogis(0.1), 0), f, method = "Nelder-Mead", control = ct) optim(o$par, f, method = "Nelder-Mead", control = ct)$value } prof_rn <- function(psi_fix, d, K) { lam_fix <- -log(1 - psi_fix) f <- function(q) rn_nll(c(log(lam_fix), q), d, K) optimise(f, c(-8, 4))$objective } psi_grid <- seq(0.70, 0.999, length.out = 26) dev_mx <- 2 * (sapply(psi_grid, prof_psi, d = d, K = K) - g_mx$value) dev_rn <- 2 * (sapply(psi_grid, prof_rn, d = d, K = K) - f_rn$value) dev_mx_at_truth <- 2 * (prof_psi(psi_true, d, K) - g_mx$value) ``` On the same detection histories the mixture reports $\hat\psi = `r sprintf("%.3f", psi_mx)`$, meaning every single site is occupied, against the Royle-Nichols answer of `r sprintf("%.3f", psi_rn)` and a truth of `r sprintf("%.3f", psi_true)`. The AIC gap between them is `r sprintf("%.1f", abs(aic_mx - aic_rn))`, which is nothing. And the mixture's answer is not an optimiser failure: its profile deviance is monotone up to the boundary, sitting at `r sprintf("%.2f", dev_mx_at_truth)` when $\psi$ is pinned at the true value, well inside the `r sprintf("%.2f", qchisq(0.95, 1))` cutoff. The likelihood genuinely prefers "everything is occupied, some sites are just very hard" over the truth. ```{r} #| label: fig-limit #| fig-width: 8 #| fig-height: 3.4 #| fig-cap: "Left: expected detection frequencies for data with heterogeneous detection but no abundance, under three models. Right: profile deviance for occupancy on genuine Royle-Nichols data, under the Royle-Nichols model and under a two-point mixture." #| fig-alt: "Two panels. The left shows observed bars with three near-identical fitted curves. The right shows a U-shaped deviance curve with a clear minimum and a monotone curve running to the right-hand boundary." obs_h <- as.numeric(table(factor(d_het, levels = 0:K))) ex <- function(par, kind) { if (kind == "occ") R * (plogis(par[1]) * dbinom(0:K, K, plogis(par[2])) + c(1 - plogis(par[1]), rep(0, K))) else if (kind == "rn") R * colSums(dpois(Nv, exp(par[1])) * t(sapply(Nv, function(n) dbinom(0:K, K, 1 - (1 - plogis(par[2]))^n)))) else R * (plogis(par[1]) * (plogis(par[4]) * dbinom(0:K, K, plogis(par[2])) + (1 - plogis(par[4])) * dbinom(0:K, K, plogis(par[3]))) + c(1 - plogis(par[1]), rep(0, K))) } df3 <- data.frame(d = rep(0:K, 3), count = c(ex(h_oc$par, "occ"), ex(h_rn$par, "rn"), ex(h_mx$par, "mix")), model = rep(c("standard occupancy", "Royle-Nichols", "two-point mixture"), each = K + 1)) p3 <- ggplot() + geom_col(data = data.frame(d = 0:K, count = obs_h), aes(d, count), fill = te_sage, alpha = 0.55, width = 0.75) + geom_line(data = df3, aes(d, count, colour = model), linewidth = 0.7) + scale_colour_manual(values = c("standard occupancy" = te_rust, "Royle-Nichols" = te_forest, "two-point mixture" = te_gold)) + scale_x_continuous(breaks = 0:K) + labs(title = "No abundance in these data", x = "detections out of six", y = "sites", subtitle = "Royle-Nichols still reports individuals per site") + theme_te() df4 <- data.frame(psi = rep(psi_grid, 2), dev = c(dev_rn, dev_mx), model = rep(c("Royle-Nichols", "two-point mixture"), each = length(psi_grid))) p4 <- ggplot(df4, aes(psi, dev, colour = model)) + geom_line(linewidth = 0.8) + geom_hline(yintercept = qchisq(0.95, 1), linetype = "dotted", colour = "grey40") + geom_vline(xintercept = psi_true, linetype = "dashed", colour = te_ink, linewidth = 0.3) + scale_colour_manual(values = c("Royle-Nichols" = te_forest, "two-point mixture" = te_gold)) + coord_cartesian(ylim = c(0, 12)) + labs(title = "Same data, two answers", x = "occupancy", y = "profile deviance", subtitle = "dashed: truth. dotted: the 3.84 cutoff") + theme_te() two_panel(p3, p4) ``` ## What this means when you have real ticks The two failures are the same failure seen from opposite sides. Site-level variation in detection is real and estimable; its **cause** is not. Royle (2006) makes the point directly for occupancy models, and Link (2003) proved the general version for closed-population abundance: different mixing distributions can produce near-identical data and quite different inferences, and no quantity of data resolves the choice, because the choice is not in the data. So the honest reading of a Royle-Nichols fit is narrow. The $\hat\lambda$ is an estimate of mean abundance **if** you are willing to assert that abundance is the only source of site-level heterogeneity in detection, that individuals are detected independently, and that the abundance distribution is Poisson. The first of those three is an ecological claim about your field protocol, not a statistical one, and it is exactly the claim that fails when you survey along a gradient of vegetation density. Model selection does not rescue you. On the very data set generated above, the standard occupancy model has an AIC of `r sprintf("%.1f", aic(f_oc, 2))` against the Royle-Nichols model's `r sprintf("%.1f", aic_rn)`, so the AIC prefers the model that is wrong about the process and wrong about occupancy by `r sprintf("%.3f", psi_true - psi_oc)`. What the fit statistics can and cannot decide here is the subject of [checking an occupancy variant](../checking-an-occupancy-variant/), and it is not a comfortable read. What does help is design. Count something, even crudely, at a subset of sites, and abundance becomes identified by data rather than by assumption. Measure the covariates you believe drive detection and put them in the detection model, so that the heterogeneity slot is not the only place left for them to go. Neither of these is an analysis decision. ## References Link W 2003 Biometrics 59(4):1123-1130 (doi:10.1111/j.0006-341X.2003.00129.x) MacKenzie D, Nichols J, Lachman G, Droege S, Royle J, Langtimm C 2002 Ecology 83(8):2248-2255 Royle J 2006 Biometrics 62(1):97-102 (doi:10.1111/j.1541-0420.2005.00439.x) Royle J, Nichols J 2003 Ecology 84(3):777-790 Kery M, Royle J 2016 Applied Hierarchical Modeling in Ecology. Academic Press (ISBN 978-0-12-801378-6) ## Related tutorials - [Imperfect detection and occupancy](../imperfect-detection-occupancy/) - [Two-species occupancy and interactions](../two-species-occupancy-interactions/) - [False positives in occupancy models](../false-positive-occupancy-models/) - [Checking an occupancy variant](../checking-an-occupancy-variant/)