REVIEW 4 major objections 4 minor 29 references
This paper claims that separating loss calibration from detection estimation in two stages gives unbiased abundance estimates and valid intervals when surveys both miss and destroy individuals.
Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →
T0 review · deepseek-v4-flash
2026-08-01 11:52 UTC pith:PGSLHVGL
load-bearing objection The two-stage point estimator is a real contribution, but the robust variance derivation drops a needed cross term and the Bayesian simulation results are hard to trust; still worth sending to referees. the 4 major comments →
Two-Stage Estimation of Population Abundance with Robust Inference under Interacting Survey Protocols
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
Core claim
The central claim is that the interaction between two survey protocols can be handled by estimating loss and detection separately rather than jointly. The model assumes the high-accuracy method is perfect, and each session of the low-accuracy method leaves a fraction p_intact of individuals alive to be seen later; with Z sessions the reference count shrinks multiplicatively by p_intact^Z. Because of the marginal Poisson representation of the N-mixture model, stage 1 is a standard Poisson regression of control and high-accuracy counts on covariates plus Z log p_intact; it produces a counterfactual predicted abundance λ̂_pred = exp(X^T β̂) representing the population had no loss occurred. Stag
What carries the argument
The load-bearing object is the two-stage estimator built on a marginal Poisson representation. Stage 1 fits the Poisson GLM log λ = X^T β + Z log p_intact to the control and post-combing shaving counts; Z log p_intact acts as an offset-like term whose coefficient yields the log per-session survival probability. From it comes the counterfactual abundance λ̂_pred = exp(X^T β̂), the 'no loss' prediction for each host. Stage 2 fits y_low ~ Poisson(λ̂_pred(1 − (1 − p_low)^Z)) to estimate the per-session detection probability. The 'hamburger' robust variance (Eq. 25), A^{-1}(B + C Vβ C^T)A^{-⊤}, is a sandwich estimator with an extra layer that adds the first-stage variance Vβ; it is what converts
Load-bearing premise
The reference (shaving) method must be perfect—it detects every parasite that survives the combing; if shaving misses any, the loss probability p_intact absorbs those misses and both stages of the estimator are biased.
What would settle it
A controlled experiment in which hosts carry a known number of parasites (established by exhaustive dissection or artificial infestation): if the shaving method recovers fewer than the known surviving total, the perfect-reference assumption fails and the two-stage estimator will show bias in p_intact and p_low. A second, simpler falsifier: simulate data where shaving detection probability is, say, 0.9 instead of 1; coverage of the proposed 95% intervals should drop below nominal.
If this is right
- If the framework is correct, ecologists can pair a cheap, loss-inducing survey with a high-accuracy reference survey and obtain unbiased abundance estimates and trustworthy confidence intervals.
- The point estimate of detection probability no longer conflates missed detections with destroyed individuals, so naive calibration (dividing combed counts by post-combing shaved counts) is shown to overestimate p_low and underestimate abundance.
- The robust intervals remain at or above nominal coverage under Poisson over-dispersion and negative-binomial misspecification, whereas the joint Bayesian hierarchical model can undercover by more than 40 percentage points (e.g., 53% coverage when the nominal is 95%).
- Because the robust variance applies directly to a logistic-regression formulation, the detection model can be extended with covariates without a new variance derivation.
- In the motivating data, the method estimates p_intact ≈ 0.7 and p_low ≈ 0.04, implying a single combing session removes roughly 30% of ectoparasites; this motivates correcting combing-based abundance estimates by about 1/(1−(1−p_low)^Z) adjusted for loss.
Where Pith is reading between the lines
- This reader's extension: the perfect-reference assumption is the natural pressure point; if shaving itself has detection error, p_intact absorbs it and the whole correction is biased. A validation experiment with known parasite loads would settle this.
- The two-stage loss-then-detection template should transfer to any paired cheap/destructive method and expensive/reference method—soil cores, insect traps, or collection methods that kill specimens—where the reference count is reduced by prior sampling.
- The paper leaves the joint-model undercoverage theoretically unexplained; a testable extension is to compare alternative joint parameterizations (e.g., weakly informative priors or reparameterized latent abundances) to see whether the feedback gap is prior-driven or structural.
- For spatio-temporal or autocorrelated surveys, the authors suggest GMM/correlated estimating equations; a concrete next step is to simulate spatially correlated counts and check whether the hamburger variance still covers, since the current derivation assumes i.i.d. units.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a two-stage estimator for population abundance when a low-accuracy survey (combing) both misses individuals and causes sample loss, and a high-accuracy reference method (shaving) is assumed perfect. Stage 1 fits a Poisson GLM to control counts and post-treatment high-accuracy counts to estimate abundance parameters and the per-session survival probability pintact, then forms counterfactual predictions. Stage 2 fits a Poisson model for low-accuracy counts conditional on these predicted counterfactual abundances to estimate the detection probability plow. The paper derives a sandwich-type 'hamburger' variance estimator intended to propagate first-stage uncertainty, and compares the method to a Bayesian joint hierarchical model in simulations and an ectoparasite dataset.
Significance. If the derivation is correct, the method offers ecologists a computationally simple two-stage procedure for combining a cheap, loss-inducing survey with a high-accuracy reference, with frequentist intervals that account for generated regressors. The point-estimation logic is sound under the stated model: Stage 1 is a standard Poisson GLM with offset, Stage 2 is a Poisson MLE conditional on a consistent first-stage estimate. The simulations cover a wide grid of parameter values and two misspecification scenarios, and the real-data application is practically relevant. However, the paper's central comparison to the Bayesian joint model is not credible as presented, and the robust variance derivation requires an independence assumption that is not stated. The perfect-detection assumption for the reference method is load-bearing and is not validated or subjected to sensitivity analysis.
major comments (4)
- [§1.1, Eqs. (1)–(2)] The unbiasedness claim rests on the reference (shaving) method having detection probability exactly 1. If shaving detects parasites with probability d<1 (constant across groups), the first-stage Poisson GLM identifies λ·d, not λ, and the counterfactual prediction is biased downward by d. The Stage 2 MLE for plow is then inflated by roughly 1/d, and test-unit abundance estimates are biased low by about d. This is not covered by the misspecification simulations (overdispersion, negative binomial), and the exact d=1 assumption is not validated in the real-data example. The text only says the method is 'assumed to detect nearly all remaining parasites' (§1.1). The paper should either extend the model to include d (if identifiable), or provide a sensitivity analysis showing how estimates change as d varies below 1.
- [§3, Tables 1–2] The Bayesian joint model is correctly specified (proper priors, true data-generating process), yet Table 1 shows systematic bias that increases with sample size (e.g., N=300, pintact=0.7, plow=0.2: Bayesian means 0.747 and 0.216 vs. true 0.7 and 0.2). A correctly specified Bayesian model should be consistent and approximately unbiased for large N; the reported pattern suggests MCMC non-convergence, poor mixing, or an implementation error, rather than a property of the joint model itself. The claim that the proposed two-stage method outperforms the joint model is a central contribution, but the simulation evidence is not credible without convergence diagnostics (e.g., trace plots, effective sample size, multiple chains) or a resolution of the discrepancy. The Discussion's statement that 'the theoretical reasons are not fully understood' is insufficient.
- [Appendix, Eq. (25)] The asymptotic variance formula Vη = A^{-1}(B + C Vβ C^T) A^{-T} omits the cross-covariance terms between the second-stage scores ui and the first-stage influence functions ψi. Expanding Eq. (23) gives Var(u_i + C ψ_i) = B + C Vβ C^T + C Cov(ψ_i,u_i) + Cov(u_i,ψ_i) C^T. The paper does not state or prove that Cov(ψ_i,u_i)=0. If yhigh and ylow are conditionally independent binomial thinnings of the same latent abundance, they are not independent unconditionally, and the cross term is generally nonzero. If the intended model is multinomial thinning (detected and surviving categories are disjoint), that must be stated explicitly, and the resulting independence between yhigh and ylow should be derived. As written, the derivation is incomplete.
- [§3, Tables 3–4] The misspecification simulations are framed as support for 'robust' inference, but the proposed variance estimator is a sandwich-type correction for first-stage parameter uncertainty; it does not correct bias from misspecified first-stage mean structure. The simulations only alter the variance (overdispersion, negative binomial) while keeping the mean model correct. They do not evaluate the more serious misspecifications: omitted covariates, nonlinearity in the abundance model, or the reference method detection error discussed above. The claim that the method 'remains valid under certain forms of model misspecification' should be narrowed, and the word 'robust' should not be taken to imply robustness to mean-model misspecification.
minor comments (4)
- [Table 5] The reported SEs for pintact and plow (0.131 and 0.196) appear to be on the logit scale, while the 95% intervals are back-transformed to the probability scale. This should be stated explicitly in the table caption or text; as printed, the SEs are inconsistent with the intervals for a probability-scale parameter.
- [Section 4] Table numbering is inconsistent: the real-data point estimates and intervals are referred to as 'Table 3', but Tables 3 and 4 are already used for the misspecification coverage results. The real-data table should be renumbered (e.g., Table 5).
- [Throughout] Minor typographical issues: 'Humberger-type' should be 'hamburger-type'; 'Bayesian joint hierarchical model' is used with inconsistent word order; Table 3 in Section 4 says 'Table 3 summarizes' but refers to the wrong table; reference to Murphy and Topel (2002) lacks volume/page details.
- [Eq. (13)] The log-likelihood contribution in Eq. (13) uses y_i^low as an exponent, which is unconventional notation; should be written as y_i^low · log μ_i or with a clear definition. Also, the dependence on Z_j in q_X(p) is implicit; the definition of q_X(p) in Eq. (11) should state that it depends on Z_i.
Circularity Check
No circularity: the two-stage estimator uses separate data and estimating equations, and the robust variance explicitly propagates first-stage uncertainty.
full rationale
The paper's derivation chain is not circular. Stage 1 identifies (β, p_intact) from control counts and post-loss reference counts under varying effort Z (Eqs. 5–6, 9), while Stage 2 identifies p_low from the low-accuracy counts y_low conditional on the first-stage counterfactual intensity λ_pred (Eq. 8). These stages use distinct outcome variables: Stage 1 uses y_ctrl and y_high; Stage 2 uses y_low. The second-stage estimating equation is therefore not a restatement of the first-stage fit, and p_low is identified by separate variation in y_low. The robust variance in Eq. (25) explicitly adds the term C V_β C^T, i.e., it treats the first-stage estimate as random rather than assuming it known; in the Poisson-marginalized model the cross-covariance between the stage-1 and stage-2 score functions is zero, so the usual two-step M-estimator variance is appropriate. The perfect-detection assumption in Eqs. (1)–(2) is a substantive identification assumption: if shaving detection is imperfect, the estimates are biased, but that is a misspecification risk, not a circular reduction of the result to its inputs. The only self-citation (Katahira et al. 2022) is used as the data source and as an empirical comparison, not as a load-bearing theoretical premise; the core methodology rests on standard GLM/M-estimation theory and the paper's own closed-form derivations. Hence no circular step is present, and the appropriate score is 0.
Axiom & Free-Parameter Ledger
free parameters (3)
- β (abundance regression coefficients) =
Real-data estimate: Intercept -1.974, Sex 1.509, species dummies 3.784/0.747/1.609/4.489/-1.504, TL 0.079, HBL -0.559 (T
- p_intact (per-session survival probability under the low-accuracy method) =
Real-data estimate: 0.700 (Table 5)
- p_low (per-session detection probability of the low-accuracy method) =
Real-data estimate: 0.041 (Table 5)
axioms (4)
- domain assumption The reference (shaving) method has perfect detection (Eqs. 1–2).
- domain assumption Latent abundance is Poisson with log-linear mean, and detection/survival follow independent binomial thinning per session (Eqs. 2–4).
- standard math Regularity conditions A1–A5 in the appendix (M-estimator smoothness, invertibility, CLT) hold.
- ad hoc to paper Cov(u_i(η0,β0), ψ_i)=0 between second-stage scores and first-stage influence functions.
read the original abstract
Estimating population abundance from field surveys is often complicated by interference between multiple survey protocols. In this paper, we propose a two-stage estimation framework for abundance models in which detection processes interact, leading to both missed detections and sample loss caused by survey procedures. Our approach separates the calibration of sample loss from the estimation of detection probability, thereby avoiding the feedback and weak identifiability that can arise in fully joint hierarchical models. We further derive a sandwich-type robust variance estimator that propagates first-stage uncertainty into the second stage and remains valid under certain forms of model misspecification. Simulation studies demonstrated that the proposed method provides more reliable uncertainty quantification than a Bayesian hierarchical joint model, which tends to underestimate uncertainty even under correct specification. We illustrate the practical utility of the method using ectoparasite abundance data from the invasive Pallas's squirrel \textit{Callosciurus erythraeus}.
Reference graph
Works this paper leans on
-
[1]
Barker, R. J., Schofield, M. R., Link, W. A., & Sauer, J. R. (2018). On the reliability of N-mixture models for count data. Biometrics, 74(1), 369--377. https://doi.org/10.1111/biom.12734
-
[2]
Boulanger, J., McLELLAN, B. N., Woods, J. G., Proctor, M. F., & Strobeck, C. (2004). Sampling design and bias in DNA‐based capture‐mark–recapture population and density estimates of grizzly bears. The Journal of Wildlife Management, 68(3), 457--469. https://doi.org/10.2193/0022-541X(2004)068[0457:SDABID]2.0.CO;2
-
[3]
Chinchio, E., Crotta, M., Romeo, C., Drewe, J. A., Guitian, J., & Ferrari, N. (2020). Invasive alien species and disease risk: An open challenge in public and animal health. PLoS Pathogens, 16(10), e1008922. https://doi.org/10.1371/journal.ppat.1008922
-
[4]
Eicker, F. (1963). Asymptotic normality and consistency of the least squares estimators. The Annals of Mathematical Statistics, 34(2), 447--456. https://doi.org/10.1214/aoms/1177704156
arXiv 1963
-
[5]
Goedknegt, M. A., Feis, M. E., Wegner, K. M., Luttikhuizen, P. C., Buschbaum, C., Camphuysen, K. C., et al. (2016). Parasites and marine invasions: Ecological and evolutionary perspectives. Journal of Sea Research, 113, 11--27. https://doi.org/10.1016/j.seares.2015.12.003
-
[6]
Hone, J. (2008). On bias, precision and accuracy in wildlife aerial surveys. Wildlife Research, 35(4), 253--257. https://doi.org/10.1071/WR07144
-
[7]
Hopkins, G. H. E. (1949). The host-associations of the lice of mammals. Proceedings of the Zoological Society of London, 119, 387--604. https://doi.org/10.1111/j.1096-3642.1949.tb00888.x
arXiv 1949
-
[8]
Huber, P. J. (1967). The behavior of maximum likelihood estimates under nonstandard conditions. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, pp. 221--233. University of California Press
1967
-
[9]
Huber, P. J. (1973). Robust regression: Asymptotics, conjectures and Monte Carlo. The Annals of Statistics, 1(5), 799--821. https://doi.org/10.1214/aos/1176342503
arXiv 1973
-
[10]
Katahira, H., Eguchi, Y., Hirose, S., Ohtani, Y., Banzai, A., Ohkubo, Y., & Shimamoto, T. (2022). Spillover and spillback risks of ectoparasites by an invasive squirrel Callosciurus erythraeus in Kanto region of Japan. International Journal for Parasitology: Parasites and Wildlife, 19, 1--8. https://doi.org/10.1016/j.ijppaw.2022.07.006
-
[11]
F., & Swihart, R
Kellner, K. F., & Swihart, R. K. (2014). Accounting for imperfect detection in ecology: a quantitative review. PLOS ONE, 9(10), e111436
2014
-
[12]
Kéry, M., & Schaub, M. (2011). Bayesian Population Analysis Using WinBUGS: A Hierarchical Perspective. Academic Press
2011
-
[13]
Knape, J., Arlt, D., Barraquand, F., Berg, ., Chevalier, M., P \"a rt, T., Ruete, A., & \.Z mihorski, M. (2018). Sensitivity of binomial N-mixture models to overdispersion: The importance of assessing model fit. Methods in Ecology and Evolution, 9(10), 2102--2114. https://doi.org/10.1111/2041-210X.13062
-
[14]
Link, W. A., Schofield, M. R., Barker, R. J., & Sauer, J. R. (2018). On the robustness of N-mixture models. Ecology, 99(7), 1547--1551. https://doi.org/10.1002/ecy.2362
-
[15]
Mulero-Pázmány, M., Jenni-Eiermann, S., Strebel, N., Sattler, T., Negro, J. J., & Tablado, Z. (2017). Unmanned aircraft systems as a new source of disturbance for wildlife: A systematic review. PLOS ONE, 12(6), e0178448. https://doi.org/10.1371/journal.pone.0178448
-
[16]
Murphy, K. M., & Topel, R. H. (2002). Estimation and inference in two-step econometric models. Journal of Business & Economic Statistics Journal of Business and Economic Statistics https://doi.org/10.1198/073500102753410417
-
[17]
Newey, W. K., & McFadden, D. (1994). Large sample estimation and hypothesis testing. In R. F. Engle & D. L. McFadden (Eds.), Handbook of Econometrics, Vol. 4, pp. 2111--2245. Elsevier. https://doi.org/10.1016/S1573-4412(05)80005-4
-
[18]
Poulin, R. (2017). Invasion ecology meets parasitology: Advances and challenges. International Journal for Parasitology: Parasites and Wildlife, 6(3), 361--363. https://doi.org/10.1016/j.ijppaw.2017.03.006
-
[19]
Royle, J. A. (2004). N-mixture models for estimating population size from spatially replicated counts. Biometrics, 60(1), 108--115. https://doi.org/10.1111/j.0006-341X.2004.00142.x
Pith/arXiv arXiv 2004
-
[20]
Ryer, C. H., & Olla, B. L. (1999). Light-induced changes in the prey consumption and behavior of two juvenile planktivorous fish. Marine Ecology Progress Series, 181, 41--51. https://doi.org/10.3354/meps181041
-
[21]
Seki, Y., & Sato, T. (2022). Habitat selection of invasive alien Pallas's squirrels (Callosciurus erythraeus) in an urban habitat with small fragmented green spaces. Mammalia, 86(1), 37--43. https://doi.org/10.1515/mammalia-2021-0071
-
[22]
Thaweepworadej, P., & Evans, K. L. (2023). Squirrel and tree-shrew responses along an urbanisation gradient in a tropical mega-city---reduced biodiversity, increased hybridisation of Callosciurus squirrels, and effects of habitat quality. Animal Conservation, 26(1), 46--60. https://doi.org/10.1111/acv.12797
-
[23]
Veech, J. A., Ott, J. R., & Troy, J. R. (2016). Intrinsic heterogeneity in detection probability and its effect on N-mixture models. Methods in Ecology and Evolution, 7(9), 1019–1028. https://doi.org/10.1111/2041-210X.12566
-
[24]
White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817--838. https://doi.org/10.2307/1912934
doi:10.2307/1912934 1980
-
[25]
Zhang, L., Rohr, J., Cui, R., Xin, Y., Han, L., Yang, X., et al. (2022). Biological invasions facilitate zoonotic disease emergences. Nature Communications, 13(1), 1762. https://doi.org/10.1038/s41467-022-29378-2
-
[26]
Zigler, C. M., Watts, K., Yeh, R. W., Wang, Y., Coull, B. A., & Dominici, F. (2013). Model feedback in Bayesian propensity score estimation. Biometrics, 69(1), 263--273. https://doi.org/10.1111/j.1541-0420.2012.01830.x
arXiv 2013
-
[27]
Frontiers in Handwriting Recognition (ICFHR), 2014 14th International Conference on , pages=
Real-time segmentation of on-line handwritten arabic script , author=. Frontiers in Handwriting Recognition (ICFHR), 2014 14th International Conference on , pages=. 2014 , organization=
2014
-
[28]
Soft Computing and Pattern Recognition (SoCPaR), 2014 6th International Conference of , pages=
Fast classification of handwritten on-line Arabic characters , author=. Soft Computing and Pattern Recognition (SoCPaR), 2014 6th International Conference of , pages=. 2014 , organization=
2014
-
[29]
arXiv preprint arXiv:1804.09028 , year=
Estimate and Replace: A Novel Approach to Integrating Deep Neural Networks with Existing Applications , author=. arXiv preprint arXiv:1804.09028 , year=
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.