Checking an N-mixture model

abundance

imperfect detection

model checking

ecology tutorial

Practical diagnostics for a fitted binomial N-mixture model: a parametric-bootstrap goodness-of-fit test with an overdispersion ratio, site-level quantile residuals, and sensitivity to the abundance ceiling, all in base R.

Author

Tidy Ecology

Published

2026-08-19

The companion tutorial on N-mixture reliability showed how much can go wrong beneath a clean-looking abundance estimate: the fit can be unbounded, the negative-binomial dispersion can be unidentified, and unmodelled detection heterogeneity can inflate abundance without warning. The lesson is that a fitted N-mixture model has to be checked rather than trusted. This tutorial assembles the routine checks, coded from scratch: a parametric-bootstrap goodness-of-fit test that also returns an overdispersion ratio, quantile residuals that expose lack of fit, and a sensitivity check on the numerical ceiling. It closes on an honest limit, that these tests can pass while the abundance is still wrong.

The estimator is the marginal likelihood with the latent abundance summed out, in the same vectorised form used throughout this series:

library(ggplot2)

nll_pois <- function(par, y, K) {
  lam <- exp(par[1]); p <- plogis(par[2])
  Nv <- 0:K; lp <- dpois(Nv, lam, log = TRUE); ll <- 0
  for (i in seq_len(nrow(y))) {
    yi <- y[i, ]; idx <- (max(yi) + 1):(K + 1); Ni <- Nv[idx]
    lb <- rowSums(vapply(yi, function(v) dbinom(v, Ni, p, log = TRUE),
                         numeric(length(Ni))))
    ll <- ll + log(sum(exp(lp[idx] + lb)))
  }
  if (!is.finite(ll)) return(1e10)
  -ll
}
fit_pois <- function(y, K, start = c(log(2), 0))
  optim(start, nll_pois, y = y, K = K, method = "BFGS")

Goodness of fit by parametric bootstrap

The natural fit statistic compares each observed count with what the fitted model expects. A single count has marginal mean equal to the product of abundance and detection, so a Freeman-Tukey discrepancy, which stabilises the variance of counts, sums the squared differences of square roots between observed and expected. Its value has no standard reference distribution here, because the mixture makes the counts more variable than a plain Poisson. The parametric bootstrap supplies the reference: simulate many data sets from the fitted model, refit each, and recompute the statistic, so the observed value is judged against the spread the model itself produces. The ratio of the observed statistic to the bootstrap mean is an overdispersion ratio, written c-hat: a value near one signals an adequate fit, and a value well above one flags variance the model has not captured (MacKenzie et al. 2002; Kery and Royle 2016).

# beta-binomial detection, used to build an overdispersed data set
rbb <- function(N, p, rho) {
  a <- p * (1 - rho) / rho; b <- (1 - p) * (1 - rho) / rho
  rbinom(length(N), N, rbeta(length(N), a, b))
}
ft <- function(y, lam, p) { E <- lam * p; sum((sqrt(y) - sqrt(E))^2) }

gof_boot <- function(y, K, B = 80, boot_seed = 1) {
  f <- fit_pois(y, K); lam <- exp(f$par[1]); p <- plogis(f$par[2])
  Tobs <- ft(y, lam, p); set.seed(boot_seed); Tb <- numeric(B)
  for (b in 1:B) {
    N  <- rpois(nrow(y), lam)
    ys <- matrix(rbinom(nrow(y) * ncol(y), rep(N, ncol(y)), p), nrow(y), ncol(y))
    fb <- fit_pois(ys, K); Tb[b] <- ft(ys, exp(fb$par[1]), plogis(fb$par[2]))
  }
  list(lam = lam, p = p, Tobs = Tobs, Tb = Tb,
       pval = mean(Tb >= Tobs), chat = Tobs / mean(Tb))
}

Tv <- 3L; K <- 40L; B <- 80L
set.seed(4280)                                   # a well-specified data set
R <- 110L; Ng <- rpois(R, 2.5)
yg <- matrix(rbinom(R * Tv, rep(Ng, Tv), 0.40), R, Tv)
set.seed(4261)                                   # detection heterogeneity added
Nb <- rpois(R, 2.5); yb <- sapply(1:Tv, function(t) rbb(Nb, 0.40, 0.30))

gg <- gof_boot(yg, K, B); gb <- gof_boot(yb, K, B)

For the well-specified data the observed statistic of 157.9 sits inside the bootstrap spread, giving a p-value of 0.237 and a c-hat of 1.05: the model fits, and the abundance estimate of 2.51 matches the truth of 2.5. The overdispersed data tell a different story. The test rejects with a p-value of 0.000 and a c-hat of 1.25, and the same fit has quietly pushed abundance to 8.09, more than three times the truth, while reporting a detection probability of only 0.11. The test earns its keep: it catches the inflation the point estimate hides.

Left panel histogram of bootstrap goodness-of-fit statistics with a dashed vertical line for the observed value falling within the histogram. Right panel histogram with the dashed observed-value line lying to the right of the entire histogram. — Figure 1: Bootstrap distribution of the Freeman-Tukey statistic under the fitted model, with the observed value marked (dashed). For the well-specified data (left) the observed value falls inside the bulk, so the fit is adequate. For the overdispersed data (right) it lies far into the upper tail, so the test rejects and the overdispersion ratio exceeds one.

Quantile residuals

A single test statistic reduces the whole fit to one number. Residuals show where it fails. For discrete data the right tool is the randomised-quantile residual (Dunn and Smyth 1996): evaluate the fitted cumulative distribution just below and at each observation, draw a uniform value between them, and map it through the normal quantile function, so a correct model yields residuals that are standard normal. Knape et al. (2018) show that residuals formed on the site totals are more informative for N-mixture models than per-count residuals, because the extra correlation from detection heterogeneity accumulates across visits. The site total is convenient here: the sum of independent binomial counts that share an abundance is itself binomial in the total number of trials, so the total at a site with abundance N is Binomial with N times the number of visits as its size, and its marginal distribution follows by averaging over the Poisson abundance.

# marginal CDF of the site total S = sum of visits; S given N is Binomial(T*N, p)
Fss <- function(s, lam, p, Tv, K) { Nv <- 0:K; sum(dpois(Nv, lam) * pbinom(s, Tv * Nv, p)) }
rq_site <- function(y, lam, p, K, rseed = 99) {
  Tv <- ncol(y); S <- rowSums(y); set.seed(rseed)
  lo <- vapply(S, function(s) Fss(s - 1, lam, p, Tv, K), numeric(1))
  hi <- vapply(S, function(s) Fss(s,     lam, p, Tv, K), numeric(1))
  qnorm(pmin(pmax(runif(length(S), lo, hi), 1e-6), 1 - 1e-6))
}
rg <- rq_site(yg, gg$lam, gg$p, K); rb <- rq_site(yb, gb$lam, gb$p, K)

The residuals for the well-specified fit have a standard deviation of 1.03, close to one, with 6.4 per cent beyond the two-sigma band, near the expected five. The overdispersed fit spreads wider, a standard deviation of 1.22 and 10.9 per cent beyond two sigma, and the quantile plot shows the tails pulling away from the reference line.

Left panel quantile plot with points lying along the dashed one-to-one line. Right panel quantile plot with points departing from the line at both tails, indicating heavier tails than a normal distribution. — Figure 2: Normal quantile plots of site-level randomised-quantile residuals. Under the well-specified model (left) the points follow the reference line. Under the overdispersed model (right) the extreme residuals fan away from the line, the signature of variance the model has not absorbed.

Sensitivity to the ceiling

The reliability tutorial showed that when the fit is near unbounded, the abundance estimate climbs with the numerical ceiling K instead of settling. That makes ceiling sensitivity a fast, routine check: refit over a few values of K and confirm the estimate is stable.

set.seed(4254)                                   # low detection, unreliable
Rl <- 120L; Tl <- 2L; Nl <- rpois(Rl, 3.0)
yl <- matrix(rbinom(Rl * Tl, rep(Nl, Tl), 0.12), Rl, Tl)
Kseq <- c(30, 60, 120, 300)
lam_reliable   <- sapply(Kseq, function(K) exp(fit_pois(yg, K)$par[1]))
lam_unreliable <- sapply(Kseq, function(K) exp(fit_pois(yl, K, c(log(3), qlogis(0.2)))$par[1]))

The well-specified fit returns 2.51 at every ceiling from 30 to 300, so the answer is real. The low-detection fit gives 14.3, 34, 78 and 226 as the ceiling rises, an estimate governed by an arbitrary bound rather than by the data. Any abundance that moves with K should be discarded.

The honest limit: passing is necessary, not sufficient

A clean bill of health from these checks is reassuring but not conclusive. The bootstrap test catches gross overdispersion, as it did above with a c-hat of 1.25. It is far less sensitive to modest detection heterogeneity, which the reliability tutorial showed can still bias abundance. A mild version makes the point.

With mild heterogeneity the goodness-of-fit test passes comfortably, a p-value of 0.725 and a c-hat of 0.98, yet abundance is estimated at 2.91, 17 per cent above the truth. The abundance estimate is more sensitive to the assumption than the fit statistics are, exactly the warning of Barker et al. (2018) and Link et al. (2018). Worse, overdispersion in detection and overdispersion in abundance are partly confounded (Knape et al. 2018), so a model that passes only after a negative-binomial term is added can still carry a biased abundance. The practical stance is the one the reliability tutorial reached: run these checks, report c-hat and residuals alongside the estimate, confirm stability in K, and treat a passing test as a floor rather than a guarantee. When an absolute abundance genuinely matters, identify detection directly with marked individuals rather than leaning on counts alone.

References

Knape, J. and others. 2018. Methods in Ecology and Evolution 9(10):2102-2114 (10.1111/2041-210X.13062).

Duarte, A., Adams, M.J. and Peterson, J.T. 2018. Ecological Modelling 374:51-59 (10.1016/j.ecolmodel.2018.02.007).

Dunn, P.K. and Smyth, G.K. 1996. Journal of Computational and Graphical Statistics 5(3):236-244 (10.1080/10618600.1996.10474708).

Barker, R.J., Schofield, M.R., Link, W.A. and Sauer, J.R. 2018. Biometrics 74(1):369-377 (10.1111/biom.12734).

Link, W.A., Schofield, M.R., Barker, R.J. and Sauer, J.R. 2018. Ecology 99(7):1547-1551 (10.1002/ecy.2362).

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

Kery, M. and Royle, J.A. 2016. Applied Hierarchical Modeling in Ecology, Volume 1. Academic Press. ISBN 978-0-12-801378-6.

--- title: "Checking an N-mixture model" description: "Practical diagnostics for a fitted binomial N-mixture model: a parametric-bootstrap goodness-of-fit test with an overdispersion ratio, site-level quantile residuals, and sensitivity to the abundance ceiling, all in base R." date: "2026-08-19 10:00" categories: [R, abundance, imperfect detection, model checking, ecology tutorial] image: thumbnail.png image-alt: "Overlaid bootstrap distributions of a goodness-of-fit statistic for a passing and a failing N-mixture fit, with the observed value marked inside the bulk for one and out in the tail for the other." --- The [companion tutorial on N-mixture reliability](../n-mixture-identifiability-and-reliability/) showed how much can go wrong beneath a clean-looking abundance estimate: the fit can be unbounded, the negative-binomial dispersion can be unidentified, and unmodelled detection heterogeneity can inflate abundance without warning. The lesson is that a fitted N-mixture model has to be checked rather than trusted. This tutorial assembles the routine checks, coded from scratch: a parametric-bootstrap goodness-of-fit test that also returns an overdispersion ratio, quantile residuals that expose lack of fit, and a sensitivity check on the numerical ceiling. It closes on an honest limit, that these tests can pass while the abundance is still wrong. The estimator is the marginal likelihood with the latent abundance summed out, in the same vectorised form used throughout this series: ```{r} #| label: machinery library(ggplot2) nll_pois <- function(par, y, K) { lam <- exp(par[1]); p <- plogis(par[2]) Nv <- 0:K; lp <- dpois(Nv, lam, log = TRUE); ll <- 0 for (i in seq_len(nrow(y))) { yi <- y[i, ]; idx <- (max(yi) + 1):(K + 1); Ni <- Nv[idx] lb <- rowSums(vapply(yi, function(v) dbinom(v, Ni, p, log = TRUE), numeric(length(Ni)))) ll <- ll + log(sum(exp(lp[idx] + lb))) } if (!is.finite(ll)) return(1e10) -ll } fit_pois <- function(y, K, start = c(log(2), 0)) optim(start, nll_pois, y = y, K = K, method = "BFGS") ``` ```{r} #| label: helpers #| echo: false theme_te <- function(base_size = 12) { theme_minimal(base_size = base_size) + theme( panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e3e0d5", linewidth = 0.3), plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), axis.text = element_text(colour = "#40493c"), axis.title = element_text(colour = "#275139"), plot.title = element_text(colour = "#275139", face = "bold", size = rel(0.98)), legend.position = "top", legend.title = element_blank(), legend.text = element_text(colour = "#40493c") ) } te_forest <- "#275139"; te_moss <- "#7cc28a"; te_bark <- "#9a7b4f"; te_slate <- "#4a6b8a" two_panel <- function(pa, pb) { grid::grid.newpage() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(pa, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pb, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) } ``` ## Goodness of fit by parametric bootstrap The natural fit statistic compares each observed count with what the fitted model expects. A single count has marginal mean equal to the product of abundance and detection, so a Freeman-Tukey discrepancy, which stabilises the variance of counts, sums the squared differences of square roots between observed and expected. Its value has no standard reference distribution here, because the mixture makes the counts more variable than a plain Poisson. The parametric bootstrap supplies the reference: simulate many data sets from the fitted model, refit each, and recompute the statistic, so the observed value is judged against the spread the model itself produces. The ratio of the observed statistic to the bootstrap mean is an overdispersion ratio, written c-hat: a value near one signals an adequate fit, and a value well above one flags variance the model has not captured (MacKenzie et al. 2002; Kery and Royle 2016). ```{r} #| label: gof #| echo: true # beta-binomial detection, used to build an overdispersed data set rbb <- function(N, p, rho) { a <- p * (1 - rho) / rho; b <- (1 - p) * (1 - rho) / rho rbinom(length(N), N, rbeta(length(N), a, b)) } ft <- function(y, lam, p) { E <- lam * p; sum((sqrt(y) - sqrt(E))^2) } gof_boot <- function(y, K, B = 80, boot_seed = 1) { f <- fit_pois(y, K); lam <- exp(f$par[1]); p <- plogis(f$par[2]) Tobs <- ft(y, lam, p); set.seed(boot_seed); Tb <- numeric(B) for (b in 1:B) { N <- rpois(nrow(y), lam) ys <- matrix(rbinom(nrow(y) * ncol(y), rep(N, ncol(y)), p), nrow(y), ncol(y)) fb <- fit_pois(ys, K); Tb[b] <- ft(ys, exp(fb$par[1]), plogis(fb$par[2])) } list(lam = lam, p = p, Tobs = Tobs, Tb = Tb, pval = mean(Tb >= Tobs), chat = Tobs / mean(Tb)) } Tv <- 3L; K <- 40L; B <- 80L set.seed(4280) # a well-specified data set R <- 110L; Ng <- rpois(R, 2.5) yg <- matrix(rbinom(R * Tv, rep(Ng, Tv), 0.40), R, Tv) set.seed(4261) # detection heterogeneity added Nb <- rpois(R, 2.5); yb <- sapply(1:Tv, function(t) rbb(Nb, 0.40, 0.30)) gg <- gof_boot(yg, K, B); gb <- gof_boot(yb, K, B) ``` For the well-specified data the observed statistic of `r sprintf("%.1f", gg$Tobs)` sits inside the bootstrap spread, giving a p-value of `r sprintf("%.3f", gg$pval)` and a c-hat of `r sprintf("%.2f", gg$chat)`: the model fits, and the abundance estimate of `r sprintf("%.2f", gg$lam)` matches the truth of 2.5. The overdispersed data tell a different story. The test rejects with a p-value of `r sprintf("%.3f", gb$pval)` and a c-hat of `r sprintf("%.2f", gb$chat)`, and the same fit has quietly pushed abundance to `r sprintf("%.2f", gb$lam)`, more than three times the truth, while reporting a detection probability of only `r sprintf("%.2f", gb$p)`. The test earns its keep: it catches the inflation the point estimate hides. ```{r} #| label: fig-gof #| echo: false #| fig-width: 8.6 #| fig-height: 3.9 #| fig-cap: "Bootstrap distribution of the Freeman-Tukey statistic under the fitted model, with the observed value marked (dashed). For the well-specified data (left) the observed value falls inside the bulk, so the fit is adequate. For the overdispersed data (right) it lies far into the upper tail, so the test rejects and the overdispersion ratio exceeds one." #| fig-alt: "Left panel histogram of bootstrap goodness-of-fit statistics with a dashed vertical line for the observed value falling within the histogram. Right panel histogram with the dashed observed-value line lying to the right of the entire histogram." gdf <- rbind(data.frame(T = gg$Tb, case = "well-specified"), data.frame(T = gb$Tb, case = "overdispersed")) obs <- data.frame(case = c("well-specified", "overdispersed"), Tobs = c(gg$Tobs, gb$Tobs)) mk_hist <- function(cs, fill, obsv) { ggplot(subset(gdf, case == cs), aes(T)) + geom_histogram(bins = 22, fill = fill, colour = "#f5f4ee", linewidth = 0.2) + geom_vline(xintercept = obsv, linetype = "dashed", colour = te_bark, linewidth = 0.9) + labs(title = if (cs == "well-specified") "Well-specified model" else "Overdispersed data", x = "Freeman-Tukey statistic", y = "Bootstrap frequency") + theme_te() } two_panel(mk_hist("well-specified", te_forest, gg$Tobs), mk_hist("overdispersed", te_slate, gb$Tobs)) ``` ## Quantile residuals A single test statistic reduces the whole fit to one number. Residuals show where it fails. For discrete data the right tool is the randomised-quantile residual (Dunn and Smyth 1996): evaluate the fitted cumulative distribution just below and at each observation, draw a uniform value between them, and map it through the normal quantile function, so a correct model yields residuals that are standard normal. Knape et al. (2018) show that residuals formed on the site totals are more informative for N-mixture models than per-count residuals, because the extra correlation from detection heterogeneity accumulates across visits. The site total is convenient here: the sum of independent binomial counts that share an abundance is itself binomial in the total number of trials, so the total at a site with abundance N is Binomial with N times the number of visits as its size, and its marginal distribution follows by averaging over the Poisson abundance. ```{r} #| label: resid #| echo: true # marginal CDF of the site total S = sum of visits; S given N is Binomial(T*N, p) Fss <- function(s, lam, p, Tv, K) { Nv <- 0:K; sum(dpois(Nv, lam) * pbinom(s, Tv * Nv, p)) } rq_site <- function(y, lam, p, K, rseed = 99) { Tv <- ncol(y); S <- rowSums(y); set.seed(rseed) lo <- vapply(S, function(s) Fss(s - 1, lam, p, Tv, K), numeric(1)) hi <- vapply(S, function(s) Fss(s, lam, p, Tv, K), numeric(1)) qnorm(pmin(pmax(runif(length(S), lo, hi), 1e-6), 1 - 1e-6)) } rg <- rq_site(yg, gg$lam, gg$p, K); rb <- rq_site(yb, gb$lam, gb$p, K) ``` The residuals for the well-specified fit have a standard deviation of `r sprintf("%.2f", sd(rg))`, close to one, with `r sprintf("%.1f", 100 * mean(abs(rg) > 1.96))` per cent beyond the two-sigma band, near the expected five. The overdispersed fit spreads wider, a standard deviation of `r sprintf("%.2f", sd(rb))` and `r sprintf("%.1f", 100 * mean(abs(rb) > 1.96))` per cent beyond two sigma, and the quantile plot shows the tails pulling away from the reference line. ```{r} #| label: fig-resid #| echo: false #| fig-width: 8.6 #| fig-height: 3.9 #| fig-cap: "Normal quantile plots of site-level randomised-quantile residuals. Under the well-specified model (left) the points follow the reference line. Under the overdispersed model (right) the extreme residuals fan away from the line, the signature of variance the model has not absorbed." #| fig-alt: "Left panel quantile plot with points lying along the dashed one-to-one line. Right panel quantile plot with points departing from the line at both tails, indicating heavier tails than a normal distribution." qq_df <- function(r) { n <- length(r); data.frame(th = qnorm(ppoints(n)), sr = sort(r)) } mk_qq <- function(r, ttl, col) { ggplot(qq_df(r), aes(th, sr)) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = te_bark) + geom_point(colour = col, size = 1.3, alpha = 0.75) + coord_cartesian(ylim = c(-4, 4)) + labs(title = ttl, x = "Theoretical quantile", y = "Residual quantile") + theme_te() } two_panel(mk_qq(rg, "Well-specified model", te_forest), mk_qq(rb, "Overdispersed data", te_slate)) ``` ## Sensitivity to the ceiling The reliability tutorial showed that when the fit is near unbounded, the abundance estimate climbs with the numerical ceiling `K` instead of settling. That makes ceiling sensitivity a fast, routine check: refit over a few values of `K` and confirm the estimate is stable. ```{r} #| label: ksens #| echo: true set.seed(4254) # low detection, unreliable Rl <- 120L; Tl <- 2L; Nl <- rpois(Rl, 3.0) yl <- matrix(rbinom(Rl * Tl, rep(Nl, Tl), 0.12), Rl, Tl) Kseq <- c(30, 60, 120, 300) lam_reliable <- sapply(Kseq, function(K) exp(fit_pois(yg, K)$par[1])) lam_unreliable <- sapply(Kseq, function(K) exp(fit_pois(yl, K, c(log(3), qlogis(0.2)))$par[1])) ``` The well-specified fit returns `r sprintf("%.2f", lam_reliable[1])` at every ceiling from `r Kseq[1]` to `r Kseq[length(Kseq)]`, so the answer is real. The low-detection fit gives `r sprintf("%.1f", lam_unreliable[1])`, `r sprintf("%.0f", lam_unreliable[2])`, `r sprintf("%.0f", lam_unreliable[3])` and `r sprintf("%.0f", lam_unreliable[length(lam_unreliable)])` as the ceiling rises, an estimate governed by an arbitrary bound rather than by the data. Any abundance that moves with `K` should be discarded. ## The honest limit: passing is necessary, not sufficient A clean bill of health from these checks is reassuring but not conclusive. The bootstrap test catches gross overdispersion, as it did above with a c-hat of `r sprintf("%.2f", gb$chat)`. It is far less sensitive to modest detection heterogeneity, which the reliability tutorial showed can still bias abundance. A mild version makes the point. ```{r} #| label: mild #| echo: false set.seed(7001) Nm <- rpois(R, 2.5); ym <- sapply(1:Tv, function(t) rbb(Nm, 0.40, 0.08)) gm <- gof_boot(ym, K, B) ``` With mild heterogeneity the goodness-of-fit test passes comfortably, a p-value of `r sprintf("%.3f", gm$pval)` and a c-hat of `r sprintf("%.2f", gm$chat)`, yet abundance is estimated at `r sprintf("%.2f", gm$lam)`, `r sprintf("%.0f", 100 * (gm$lam / 2.5 - 1))` per cent above the truth. The abundance estimate is more sensitive to the assumption than the fit statistics are, exactly the warning of Barker et al. (2018) and Link et al. (2018). Worse, overdispersion in detection and overdispersion in abundance are partly confounded (Knape et al. 2018), so a model that passes only after a negative-binomial term is added can still carry a biased abundance. The practical stance is the one the reliability tutorial reached: run these checks, report c-hat and residuals alongside the estimate, confirm stability in `K`, and treat a passing test as a floor rather than a guarantee. When an absolute abundance genuinely matters, identify detection directly with marked individuals rather than leaning on counts alone. ## Related tutorials - [The abundance-detection trade-off and N-mixture reliability](../n-mixture-identifiability-and-reliability/) - [N-mixture models: abundance from repeated counts](../n-mixture-models-abundance/) - [Bayesian N-mixture models](../bayesian-n-mixture-abundance/) - [Checking a bootstrap](../checking-a-bootstrap/) ## References Knape, J. and others. 2018. Methods in Ecology and Evolution 9(10):2102-2114 (10.1111/2041-210X.13062). Duarte, A., Adams, M.J. and Peterson, J.T. 2018. Ecological Modelling 374:51-59 (10.1016/j.ecolmodel.2018.02.007). Dunn, P.K. and Smyth, G.K. 1996. Journal of Computational and Graphical Statistics 5(3):236-244 (10.1080/10618600.1996.10474708). Barker, R.J., Schofield, M.R., Link, W.A. and Sauer, J.R. 2018. Biometrics 74(1):369-377 (10.1111/biom.12734). Link, W.A., Schofield, M.R., Barker, R.J. and Sauer, J.R. 2018. Ecology 99(7):1547-1551 (10.1002/ecy.2362). MacKenzie, D.I., Nichols, J.D., Lachman, G.B., Droege, S., Royle, J.A. and Langtimm, C.A. 2002. Ecology 83(8):2248-2255. Kery, M. and Royle, J.A. 2016. Applied Hierarchical Modeling in Ecology, Volume 1. Academic Press. ISBN 978-0-12-801378-6.