# Population Size Estimation with Many Lists and Heterogeneity: A Conditional
Log-Linear Model Among the Unobserved

Mateo Dulce Rubio  Edward Kennedy

(Carnegie Mellon University)

###### Abstract

We contribute a general and flexible framework to estimate the size of a
closed population in the presence of $K$ capture-recapture lists and
heterogeneous capture probabilities. Our novel identifying strategy leverages
the fact that it is sufficient for identification that a subset of the $K$
lists are not arbitrarily dependent within the subset of the population
unobserved by the remaining lists, conditional on covariates. This
identification approach is interpretable and actionable, interpolating between
the two predominant approaches in the literature as special cases:
(conditional) independence across lists and log-linear models with no highest-
order interaction. We derive nonparametric doubly-robust estimators for the
resulting identification expression that are nearly optimal and approximately
normal for any finite sample size, even when the heterogeneous capture
probabilities are estimated nonparametrically using machine learning methods.
Additionally, we devise a sensitivity analysis to show how deviations from the
identification assumptions affect the resulting population size estimates,
allowing for the integration of domain-specific knowledge into the
identification and estimation processes more transparently. We empirically
demonstrate the advantages of our method using both synthetic data and real
data from the Peruvian internal armed conflict to estimate the number of
casualties. The proposed methodology addresses recent critiques of capture-
recapture models by allowing for a weaker and more interpretable identifying
assumption and accommodating complex heterogeneous capture probabilities
depending on high-dimensional or continuous covariates.

##  1 Introduction

Estimating the size of a population is a critical challenge across various
domains, such as biology for estimating animal populations and conservation
efforts (Schwarz and Seber,, 1999; Burgar et al.,, 2018), public health for
estimating disease prevalence (Hook and Regal,, 1995; Böhning et al.,, 2020),
and human rights for estimating the number of victims of specific crimes
(Manrique-Vallier et al.,, 2013; Silverman,, 2020). In many high-stakes
applications, we observe capture-recapture data, where each capture episode
involves a subsample from the population and records it in a list. In such
data, the same unit can be observed multiple times, while some units are never
captured. Intuitively, the intersection between capture lists can be exploited
to estimate the total population size, as smaller populations tend to have
units appearing in multiple lists more frequently than larger populations.

Various approaches have been used for estimating the size of a population from
multiple-list capture-recapture data. Among these methodologies, one popular
approach assumes that capture episodes are independent and that all units in
the target population have the same probability of being observed. Under these
assumptions, the observed data can be used to estimate the number of
unobserved units by extrapolating the estimated probabilities of each observed
capture profile (Bishop et al.,, 2007; Otis et al.,, 1978). However, the
assumptions of independence and homogeneity can be implausible in real-world
applications, which limit the applicability and robustness of the resulting
estimators. For instance, in settings where demographic or biological
attributes make some units more or less easy to capture, assuming homogeneous
capture probabilities may bias the estimated population size in arbitrary
ways. Therefore, there is a pressing need for innovative statistical
methodologies in the face of recent and challenging applications, such as
estimating the population size of HIV patients (Wesson et al.,, 2024),
casualties during armed conflicts (Manrique-Vallier et al.,, 2021), or victims
of modern slavery (Binette and Steorts,, 2022), particularly in the presence
of heterogeneous capture probabilities and complex dependencies between
capture episodes.

To move away from the assumptions of homogeneity and independence, more
general approaches have been proposed, such as Bayesian latent class models
(Manrique-Vallier,, 2016), and methods leveraging semiparametric efficiency
theory (You et al.,, 2021; Das et al.,, 2023). However, these methods rely on
identification assumptions that can be difficult to understand, making it
challenging to argue for their validity in practical applications.
Furthermore, they sometimes conflate issues of identification and estimation
(Aleshin-Guendel et al.,, 2024), or use only part of the available information
collected during the capture episodes. Despite these advances, there remains a
need for a general framework that allows for an interpretable, useful, and
plausible identifying assumption and accommodates heterogeneous capture
probabilities in a flexible manner, even when this heterogeneity depends on
high-dimensional or continuous attributes. Developing such a general and
flexible framework for population size estimation is the ultimate goal of our
work.

###  1.1 Contributions

Concretely, we contribute to the population size estimation literature in four
key ways:

  * •

We introduce a novel conditional log-linear model that incorporates
heterogeneous capture probabilities according to a set of available
covariates, allowing for high-dimensional and continuous attributes.

  * •

We propose a new identification strategy by noting that it is sufficient for
identification that a subset of the capture-recapture lists are not
arbitrarily dependent within the subset of the population unobserved by the
remaining lists, conditional on covariates. This strategy is both
interpretable and actionable, generalizing the two predominant approaches in
the literature: conditional independence between capture episodes and
traditional log-linear models with no highest-order interaction.

  * •

We derive nonparametric doubly-robust estimators for the capture probability
and the total population size. The proposed one-step estimators are shown to
be nearly minimax optimal and approximately normal for any finite sample size,
providing robust and reliable estimates even under flexible estimation of the
heterogeneous capture probabilities using machine learning models.

  * •

Finally, we introduce a sensitivity analysis framework that reveals how
violations of the identification assumptions impact the population size
estimates. This allows researchers to incorporate domain-specific knowledge
into the estimation process more transparently and rigorously.

To validate our methodology, we conduct empirical evaluations using both
synthetic data and real data from the Peruvian internal armed conflict. Our
results highlight the advantages of our approach in accurately estimating
population size in real-world applications, allowing for a weaker and more
interpretable identifying assumption and accommodating complex heterogeneous
capture probabilities.

##### Paper outline

Section 2 describes the setup of the capture-recapture problem and the
predominant identification strategies used in previous work. Section 3
presents our conditional log-linear model for heterogeneous capture
probabilities and introduces our novel identification strategy to express the
capture probability using the observed data. In Section 4, we derive a one-
step estimator using semiparametric efficiency theory and show that is nearly
optimal and approximately normal for any finite sample size, even under
flexible estimation of the heterogeneous capture probabilities. In Section 5,
we devise a sensitivity analysis and construct analogous one-step estimators
that bound the capture probability, while Section 6 extends our results to the
total population size. Section 7 illustrates the advantages of the proposed
methods, over a natural plug-in estimator, in simulated and real data from the
Peruvian internal armed conflict during 1980-2000. Finally, Section 8
concludes and discusses the broad applicability of the proposed framework for
population size estimation across various domains.

##  2 Problem Setup

Consider a closed population of unknown size $n$ that we aim to estimate.
Suppose that we have access to $K$ partial lists, each capturing different
subsets of this population. Then, each uniquely identified subject
$i\in\\{1,\dots,n\\}$ has an associated capture profile vector:

$Y_{i}=(Y_{i1},\dots,Y_{iK}),$

where $Y_{ik}\in\\{0,1\\}$ indicates whether unit $i$ appears on list
$k\in\\{1,\dots,K\\}$. In addition, let us assume that we observe a set of
covariates $X_{i}\in\mathbb{R}^{d}$ for each subject $i$. Crucially, note that
we only observe the $N\leq n$ units for whom $Y_{i}\neq 0$, where $0$ denotes
the $K$-dimensional zero vector. That is, we only observe those that appear in
at least one of the $K$ lists. Our goal is then to use these $N$ capture-
recapture observations to estimate the total population size $n$.

In this setup, the total number of unique observed units,
$N=\sum_{i=1}^{n}\mathbbm{1}(Y_{i}\neq 0)$, is a random variable that follows
a binomial distribution: $N\sim\text{Bin}(n,\psi)$, where
$\psi=\mathbb{P}(Y\neq 0)$ is the capture probability. Treating the population
size $n$ as a fixed parameter satisfying $\mathbb{E}[N]=\psi n$, we can derive
a straightforward estimator for $n$ based on the observed sample size $N$:

$\widehat{n}=\frac{N}{\widehat{\psi}}.$ (1)

Therefore, estimating the population size $n$ reduces to estimating the
inverse capture probability $\psi^{-1}$ from the observed data.

However, identifying and estimating $\psi=1-\mathbb{P}(Y=0)$ is notably
challenging precisely because we do not observe any of those units with
$Y_{i}=0$. Intuitively, we aim to estimate $\mathbb{P}(Y=0)$ by extrapolating
the information from those with $Y\neq 0$ (Manrique-Vallier et al.,, 2021),
which may be arbitrarily different populations. This is a hopeless task unless
we introduce additional identifying assumptions that allow us to use the
information from the observed distribution to identify and estimate the
population parameter $\psi$ using the observed data.111We refer the reader to
Aleshin-Guendel et al., (2024) for a detailed discussion on the role of the
identification assumption in population size estimation.

Formally, we assume that the covariate and conditional list membership
distributions are the same for all units, and that these probabilities do not
depend on other units’ variables. This implies that the random vectors
$(X_{i},Y_{i})$ are iid with respect to some (unknown) underlying distribution
$\mathbb{P}$. Crucially, the capture-recapture setting suffers from biased
sampling in the sense that the observed data $\\{(X_{i},Y_{i})\\}_{i=1}^{N}$
are actually iid draws from the conditional distribution

$\mathbb{Q}(Y=y,X=x)\equiv\mathbb{P}(Y=y,X=x\mid Y\neq 0).$

That is, the observed distribution $\mathbb{Q}$ captures the joint
distribution of $(X,Y)$ conditional on them being observed, i.e., given $Y\neq
0$.

To deal with such biased sampling, let $q_{y}(x)=\mathbb{Q}(Y=y\mid X=x)$
denote the conditional probability that the list membership vector
$Y=(Y_{1},\dots,Y_{K})$ equals the specific capture profile
$y=(y_{1},\dots,y_{k})\in\\{0,1\\}^{K}$, given covariates $X=x$. These measure
the capture heterogeneity in the target population, and we will refer to the
set $\\{q_{y}(X)=\mathbb{Q}(Y=y\mid X)\\}_{y\neq 0}$ as $q$-probabilities (Das
et al.,, 2023). Then, note that for any $y\neq 0$

$\mathbb{P}(Y=y\mid X=x)=\mathbb{P}(Y=y\mid X=x,Y\neq 0)\mathbb{P}(Y\neq 0\mid
X=x)=q_{y}(x)\gamma(x),$ (2)

where $\gamma(x)=\mathbb{P}(Y\neq 0\mid X=x)$ is the conditional capture
probability for a subject with covariates $X=x$. Finally, note that we can
express $\psi^{-1}=\mathbb{P}(Y\neq 0)^{-1}$ as the harmonic mean of
$\gamma(X)$ under distribution $\mathbb{Q}$ (Das et al.,, 2023; Johndrow et
al.,, 2019):

$\frac{1}{\psi}=\mathbb{E}_{\mathbb{Q}}\left[\frac{1}{\gamma(X)}\right]=\int\frac{1}{\gamma(x)}\mathbb{Q}(X=x)dx,$
(3)

where $\mathbb{Q}(X=x)=\mathbb{P}(X=x\mid Y\neq 0)$. Therefore, if we can
identify the conditional capture probability $\gamma(X)$ from the available
data, we can use the previous identity (3) to identify the inverse capture
probability $\psi^{-1}$ using the observed data from $\mathbb{Q}$.

###  2.1 Identification strategies in capture-recapture settings

In this section we summarize two popular identification strategies used in the
capture-recapture literature to identify the capture probability and the total
population size: conditional independence given covariates, and no highest-
order interaction in a log-linear model with many lists. Intuitively, since
the capture-recapture setting suffers from biased sampling, it is necessary to
introduce additional assumptions that enable us to extrapolate the information
from distribution $\mathbb{Q}$, observed for those that appear at least on one
list, to the distribution $\mathbb{P}$, in particular for those who are
unobserved. Crucially, if lists are irreparably dependent we will not be able
to identify the target population size (Aleshin-Guendel et al.,, 2024).

For instance, in the $K=2$ list case without covariates, assuming that the
probability of appearing in one list is independent of appearing on the other
list allows for the total population size to be estimated directly using the
Lincoln-Petersen estimator (Lincoln,, 1930; Petersen,, 1896). This assumption
is untestable and can be implausible in practical settings, but is relaxed by
measuring covariates that potentially affect the probability of being
captured. Then, using the identification assumption of conditional
independence between lists, the population size $n$ can be correctly
identified and estimated from observed data (Rivest and Baillargeon,, 2014;
Das et al.,, 2023). Nonetheless, this is still considered a strict
identification assumption, implicitly assuming that all relevant variables can
be measured, and typically only allowing for the use of two lists.

An alternative common approach when $K$ capture-recapture lists are available
is to assume an underlying log-linear model. These models assume that the
probability of observing each capture profile vector $y$ satisfies (Fienberg,,
1972)

$\log(\mathbb{P}(Y=y))=\alpha_{0}+\sum_{y^{\prime}\neq
0}\alpha_{y^{\prime}}\prod_{k:\leavevmode\nobreak\ y^{\prime}_{k}=1}y_{k},$
(4)

for some parameters $\alpha_{y}$, where the sum is over all possible non-zero
capture profiles $y^{\prime}\in\\{0,1\\}^{K}$. In this setting, log-linear
models are unidentifiable given that there are $2^{K}$ parameters $\alpha_{y}$
to be estimated from $2^{K}-1$ observed capture profiles; all but $y=0$.
Therefore, the standard approach is to assume no highest-order interaction,
$\alpha_{1}=0$, which assumes that the probability $\mathbb{P}(Y=1)$ can be
expressed using only lower-level interaction terms captured by the
$\alpha_{y}$ parameters, for $y\neq 1$ (Fienberg,, 1972; Bishop et al.,,
2007). We provide a more detailed interpretation of the no highest-order
interaction assumption in terms of conditional odds ratios in the Appendix.
However, these models require overlap between all combinations of lists and
typically assume homogeneous capture probabilities, meaning that all units
have the same probability of being observed, which can be considered strong
assumptions. For instance, in settings where demographic or biological
attributes make some units more or less easy to capture, assuming homogeneous
capture probabilities may bias the estimated population size in complicated
ways (Manrique-Vallier et al.,, 2021).

In summary, there are two predominant approaches for identifying the total
population size from capture-recapture data: (i) measuring all possible
features and invoking a conditional independence assumption (Das et al.,,
2023), or (ii) collecting multiple lists and using a log-linear model with no
highest-order interaction term (You et al.,, 2021). In general, for any number
of $K$ lists, assuming that two lists are conditionally independent given the
other lists and covariates is a stricter identification assumption than the no
highest-order interaction in a conditional log-linear model. The latter
identifying assumption accounts for more complex capture-recapture
dependencies across lists, but is harder to interpret for a number of lists
larger than $K>3$. In addition, standard log-linear models requires positivity
for all capture profiles and cannot account for heterogeneous capture
probabilities.

In the next section, we demonstrate how the two predominant strategies can be
integrated by using conditional log-linear models combined with a weaker and
more interpretable identification assumption. The proposed identification
strategy generalizes the two predominant approaches noting that it is
sufficient for identification that a subset of the $K$ lists has no highest-
order interaction, conditional on not being observed by the remaining lists
and covariates. For instance, the assumption

$Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid X,Y_{3}=0,\dots,Y_{K}=0,$

is sufficient to identify the capture probability $\psi$. Note that the no
highest-order interaction assumption and conditioning only on not being
observed by the remaining lists do not impose any restrictions on the observed
distribution $\mathbb{Q}$, allowing for flexible nonparametric data-driven
modeling of $\psi$.

##  3 Conditional Log-Linear Model Among the Unobserved

In this section, we give a new approach that incorporates heterogeneous
capture probabilities into a conditional log-linear model where list
membership varies across covariates $X$. We introduce a new identification
approach by assuming there is a subset of $J\leq K$ lists that are not
arbitrarily dependent given covariates, in the subset of the population
unobserved by the remaining $K-J$ lists. We show that this identifying
assumption is sufficient to express the target capture probability using the
observed data, and that it generalizes the two most common identification
strategies used in previous literature. In short, we derive an identifying
expression for $\gamma(X)$ from the proposed conditional log-linear model
among the unobserved, which we use to identify the capture probability using
(3) and propose efficient estimators for $\psi^{-1}$. Finally, we discuss how
to interpret and use the identifying assumption introduced in our work.

In more detail, we first fix a subset $Y_{1},\dots,Y_{J}$ of $J$ lists ($J\leq
K$).222For simplicity, we will index the subset of $J$ lists as
$Y_{1},\dots,Y_{J}$. But we do not assume these are captured before lists
$Y_{J+1},\dots,Y_{K}$ or in any particular order. Then, we use a conditional
log-linear model to express the (log) conditional probability that
$(Y_{1},\dots,Y_{J})=(y_{1},\dots,y_{J})$, conditional on covariates $X$ and
on not being present in the other $K-J$ lists. That is, for those unobserved
by lists $Y_{J+1},\dots,Y_{K}$. Mathematically,

$\log\left[\mathbb{P}\left\\{(Y_{1},\dots,Y_{J})=(y_{1},\dots,y_{J})\mid
X,(Y_{J+1},\dots,Y_{K})=0\right\\}\right]=\alpha_{0}(X)+\sum_{(y^{\prime}_{1},\dots,y^{\prime}_{J})\neq
0}\alpha_{y^{\prime}}(X)\prod_{j:\leavevmode\nobreak\ y^{\prime}_{j}=1}y_{j},$
(5)

where the $\alpha_{y}(X)$ are arbitrary. Notably, treating $\alpha_{y}(X)$ as
parameters depending on the covariates $X$ allows us to account for
heterogeneous capture probabilities of the subjects in the target population,
including high-dimensional and continuous attributes.

######  Remark 1.

To ensure that our log-linear model is well-defined, formally it is required
that $\mathbb{P}(X,(Y_{J+1},\dots,Y_{K})=0)>0$. This holds by design in the
capture-recapture setting because we assume that there is still a portion of
the population not captured by the available lists, and in particular by the
$J+1,\dots,K$ lists. That is, $\mathbb{P}((Y_{J+1},\dots,Y_{K})=0\mid X)>0$.
Otherwise, we could use the observed records directly as the total population
size.

###  3.1 Identification in a conditional log-linear model

To identify the conditional capture probability we follow the standard
identification strategy in log-linear models and assume that the highest-order
interaction parameter in (5) is zero, but only conditionally on covariates $X$
(Fienberg,, 1972). Recall that the highest-order coefficient $\alpha_{1}(X)$
in a log-linear model captures higher order interactions between the capture
episodes, and it is used to encode certain lack of dependence in the
population distribution $\mathbb{P}$. This assumption of no arbitrary
dependence is necessary because, as explained above, otherwise the observed
data distribution $\mathbb{Q}$ would be uninformative about the target
estimand $\psi=\mathbb{P}(Y\neq 0)$.

Building on the saturated log-linear model (5), we can express the highest-
order interaction coefficient, $\alpha_{1}(X)$, in terms of the
$q$-probabilities and the conditional capture probability $\gamma(X)$, using
the relationship between the true distribution $\mathbb{P}$ and the observed
distribution $\mathbb{Q}$ in (2):

$\alpha_{1}(X)=\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{J+|y|}\log(q_{y}(X)\gamma(X))+(-1)^{J}\log(1-\gamma(X)),$ (6)

where $|y|=\sum_{j=1}^{J}y_{j}$ is the number of positive lists in the capture
profile $y=(y_{1},\dots,y_{J})$ ($L^{1}-$norm) (Bishop et al.,, 2007; You et
al.,, 2021; Binette and Steorts,, 2022). Then, assuming $\alpha_{1}(X)=0$, we
obtain the following point-identification result for the capture probability
$\psi$. We defer all the proofs in our work to the Appendix.

######  Proposition 1.

(Identification Result) Assume no highest-order interaction in a conditional
log-linear model for $J$ lists conditional on covariates and on not being
captured by the remaining $K-J$ lists. Assume positivity $q_{y}(X)>0$, with
probability 1. Then, we can identify the inverse capture probability
$\psi^{-1}$ in terms of the observed data by

$\frac{1}{\psi}=\mathbb{E}_{\mathbb{Q}}\left[\frac{1}{\gamma(X)}\right]=\mathbb{E}_{\mathbb{Q}}\left[1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\right)\right],$ (7)

where the $q$-probabilities
$q_{y}(X)=\mathbb{Q}((Y_{1},\dots,Y_{J},Y_{J+1},\dots,Y_{K})=(y_{1},\dots,y_{J},0,\dots,0)\mid
X)$ are computed for the complete $K$-dimensional capture profiles $y$ of the
form $y=(y_{1},\dots,y_{J},0,\dots,0)$.

Our identification strategy is weaker than assuming a log-linear model for all
$K$ lists with no highest-order interaction in two ways. First, when $J=K$ we
recover the standard log-linear model (4), conditional on covariates $X$,
which assumes a lack of dependence between all available $K$ lists in the
target population. Notably, our proposed approach only requires that $J\leq K$
lists are not arbitrarily dependent in the subset of the population unobserved
by the rest of the $K-J$ lists, conditional on covariates. Moreover, our
identifying assumption captures the right intuition that identification is
easier when (i) there are more relevant covariates $X$ available, or (ii)
there is a larger number $K$ of lists that already capture a portion of the
target population.

Secondly, note that we only require positivity of the $q$-probabilities for
the $2^{J}-1$ capture profiles of the form
$y=(y_{1},\dots,y_{J},0,\dots,0)\neq 0$. On the other hand, using a
conditional log-linear model on all of the $K$ lists needs the $2^{K}-1$
$q$-probabilities for all possible capture profiles $y=(y_{1},\dots,y_{K})\neq
0$ to be bounded away from zero with probability one. In the presence of many
lists, this is a strong assumption due to potential lack of overlap between
all possible combinations of lists. For instance, if $J=3$ and $K=5$, our
approach assumes positivity for seven $q$-probabilities, while traditional
log-linear models need positivity to hold for $31$ of them. This has been
studied in the capture-recapture literature using data models with structural
zeros (Bishop et al.,, 2007). However, such methods typically introduce
additional parametric assumptions mixing the identification and the estimation
of the total population size. Our method overcomes this problem
nonparametrically and does not require overlap between all $K$ lists, but only
between the $J$ lists present in the conditional log-linear model.

Another special case arises when $J=2$. In this scenario, even in the presence
of a large $K$, our identification strategy reduces to assuming conditional
independence between two lists in the population unobserved by
$Y_{3},\dots,Y_{K}$, given covariates. The obtained identification expression
resembles the result in Das et al., (2023), which can be seen as a
nonparametric conditional Lincoln-Peterson estimator for heterogeneous
populations. Crucially, unlike the Das et al., (2023) estimator, our method
uses the information from all $K$ lists in both the identification and
estimation procedures. We summarize this special case in the following
Corollary as we believe it can be useful in many settings where there are two
larger lists that are believed to be (conditionally) independent.

######  Corollary 1.

When $J=2$, our identification assumption is equivalent to assuming that two
lists are conditionally independent given covariates for those not captured by
$Y_{3},\dots,Y_{K}$:

$Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid X,Y_{3}=0,\dots,Y_{K}=0.$ (8)

In this case, identification of the capture probability $\psi$ simplifies to

$\frac{1}{\psi}=\mathbb{E}_{\mathbb{Q}}\left[\frac{1}{\gamma(X)}\right]=\mathbb{E}_{\mathbb{Q}}\left[\frac{q_{10}(X)q_{01}(X)}{q_{11}(X)}+1\right],$
(9)

where
$q_{y_{1}y_{2}}(X)=\mathbb{Q}(Y_{1}=y_{1},Y_{2}=y_{2},Y_{3}=0,\dots,Y_{K}=0\mid
X)$.

######  Remark 2.

The conditional log-linear model and the no highest-order interaction
assumption are used to identify the target parameter $\psi^{-1}$ by assuming a
lack of dependence of some of the capture episodes in the population
distribution $\mathbb{P}$. Nonetheless, this does not impose any restrictions
on the observed distribution $\mathbb{Q}$, such that our model is
nonparametric in that we do not constrain the nuisance functions
$q_{y}(x)=\mathbb{Q}(Y=y\mid X=x)$ needed to identify $\psi$. Therefore, we
can leverage semiparametric efficiency theory to propose efficient estimators
without introducing additional potentially restrictive assumptions, allowing
for flexible estimation of the nuisance heterogeneous $q$-probabilities.

To summarize, we have introduced a novel identification strategy assuming that
a subset of $J\leq K$ lists are not arbitrarily dependent within the
population not captured by the rest of the lists, conditional on covariates.
This identification assumption generalizes the two predominant approaches in
the literature: conditional independence between capture episodes and
traditional log-linear models with no highest-order interaction between $K$
lists. Moreover, it is interpretable and actionable in the sense that to make
identification of the capture probability easier, an analyst can either
measure more relevant covariates or collect additional lists. In the next
section, we use the obtained identifying expression to derive the efficiency
bound to estimate $\psi^{-1}$ and propose a “doubly-robust” estimator that
closely attains such bound and is approximately normal in finite samples.

##  4 Efficient Nonparametric Estimation

In this section, we build on the identifying expression from Proposition 1
arising from the identification assumption of no highest-order interaction in
a conditional log-linear model among the unobserved. We explicitly derive the
Efficient Influence Function for (7) and the corresponding efficiency bound
for estimating $\psi^{-1}=\mathbb{E}_{\mathbb{Q}}[\gamma(X)^{-1}]$. We then
introduce a novel one-step estimator for the capture probability and show that
it is nearly optimal in a “doubly-robust” fashion and approximately Gaussian
for any sample size $N$.

###  4.1 Efficient influence function & efficiency bound

Before introducing our proposed estimator for the inverse capture probability,
we first benchmark the estimation problem using the semiparametric efficiency
bound. Importantly, the efficiency bound gives the best possible (scaled) Mean
Squared Error that can be achieved for estimating $\psi^{-1}$ without
introducing additional, potentially restrictive, parametric assumptions on the
nuisance functions $q_{y}(X)=\mathbb{Q}(Y=y\mid X)$. Therefore, in the absence
of a correct parametric model for the $q$-probabilities, having an estimator
that attains such efficiency bound is critical to allow for the flexible
estimation of the nuisance functions capturing heterogeneous capture
probabilities in a data-driven way.

Formally, the efficiency bound for estimating the inverse capture probability
can be computed as the variance of the Efficient Influence Function (EIF) for
$\psi^{-1}$. Therefore, we first use standard machinery from semiparametric
efficiency theory to derive the corresponding EIF for the inverse capture
probability from the identifying expression in equation (7) (cf. Kennedy,
(2022); van der Vaart, (2002); You et al., (2021)).

######  Lemma 1.

Assume (i) no highest-order interaction in a conditional log-linear model for
$J$ lists, conditional on covariates and on being unobserved by the remaining
$K-J$ lists, and (ii) positivity for the relevant $q$-probabilities. Then, the
Efficient Influence Function for
$\psi^{-1}=\mathbb{E}_{\mathbb{Q}}[\gamma(X)^{-1}]$ is given by

$\phi(X,Y)=\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}+1-\frac{1}{\psi},$ (10)

where $\gamma(X)^{-1}=1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\right)$ is the identifying expression for the
conditional capture probability. As in Proposition 1, the $q$-probabilities
are computed for the $K$-dimensional capture profiles $y$ of the form
$y=(y_{1},\dots,y_{J},0,\dots,0)$.

In the following Theorem, we use this result to compute the semiparametric
efficiency bound to estimate the target estimand of the inverse capture
probability, $\psi^{-1}=\mathbb{E}_{\mathbb{Q}}[\gamma(X)^{-1}]$. Crucially,
this benchmarks our estimation problem because it gives the best possible
scaled MSE that can be achieved for estimating $\psi^{-1}$ without introducing
additional parametric assumptions.

######  Theorem 1.

Under the same assumptions as in Lemma 1, the semiparametric efficiency bound
for estimating $\psi^{-1}=\mathbb{E}_{\mathbb{Q}}[\gamma(X)^{-1}]$ takes the
form

$\text{Var}(\phi(X,Y))=\text{Var}\left(\frac{1}{\gamma(X)}\right)+\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\sum_{(y_{1},\dots,y_{J})\neq
0}\frac{1}{q_{y}(X)}-1\right],$ (11)

where $\gamma(X)^{-1}=1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\right)$ and the $q$-probabilities are computed
for the $K$-dimensional capture profiles $y=(y_{1},\dots,y_{J},0,\dots,0)$.

Notably, the efficiency bound benchmarking our estimation problem is composed
by two terms. The first term captures the heterogeneity in the conditional
capture probabilities $\gamma(X)$, such that if these are more homogeneous we
can estimate the marginal capture probability better. The second term involves
the magnitude of these conditional probabilities and the magnitude of the
individual $q$-probabilities for the capture profiles
$y=(y_{1},\dots,y_{J},0,\dots,0)$. Intuitively, as those unobserved by
$Y_{J+1},\dots,Y_{K}$ are more easily observed under any combination of the
lists $Y_{1},\dots,Y_{J}$, we can estimate $\psi^{-1}$ more efficiently.

Note that for a standard log-linear model on $K$ lists, the resulting
efficiency bound would include the $q$-probabilities for all possible capture
profiles $y=(y_{1},\dots,y_{K})$. Therefore, when there are rare capture
profiles (i.e., $q_{y}(X)$ is small for some $y$’s), the efficiency bound is
larger and estimation is more difficult. This shows how our approach not only
allows for weaker identification assumptions, but also makes estimation
easier. This is illustrated empirically in our data analysis section using
real data from the Peruvian internal conflict.

###  4.2 Nonparametric estimation

We now build over the previous results to propose a one-step double machine-
learning estimator that allows for flexible estimation of the
$q$-probabilities, for example using machine learning models. We provide
finite-sample guarantees for our one-step estimator showing that it is close
to an empirical average of the EIF, which attains the semiparametric
efficiency bound derived in (11), with high probability (nearly optimal).
Moreover, we show that it is approximately normal allowing for standard
inference procedures. We complement our analysis by comparing the proposed
estimator to a natural, but typically sub-optimal plug-in estimator for
$\psi^{-1}$.

####  4.2.1 Plug-in estimator

To motivate the proposed one-step estimator we first introduce a simple plug-
in estimator for $\psi^{-1}$. Concretely, from the identification expression
in equation (7), a straightforward plug-in estimator to estimate
$\psi^{-1}=\mathbb{E}_{\mathbb{Q}}[\gamma(X)^{-1}]$ from observed data
$\\{(X_{i},Y_{i})\\}_{i=1}^{N}$ is computed by replacing the conditional
capture probability $\gamma(X)$ by an estimate $\widehat{\gamma}(X)$, and by
estimating the expectation $\mathbb{E}_{\mathbb{Q}}$ with the empirical
average $\mathbb{Q}_{N}$. Mathematically:

$\widehat{\frac{1}{\psi_{pi}}}=\mathbb{Q}_{N}\left(\frac{1}{\widehat{\gamma}(X_{i})}\right)=\frac{1}{N}\sum_{i=1}^{N}\left(1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(\widehat{q}_{y}(X_{i}))\right)\right),$ (12)

where $\widehat{q}_{y}(X_{i})$ are estimates of the $q$-probabilities,
$q_{y}(X_{i})=\mathbb{Q}(Y=y\mid X=X_{i})$, for each capture profile
$y=(y_{1},\dots,y_{J},0,\dots,0)\neq 0$.

Note that in constructing the plug-in estimator for the (inverse) capture
probability, one estimates conditional capture probabilities for each unit,
$\widehat{\gamma}(X_{i})$, and their empirical harmonic mean (i.e., with
respect to $\mathbb{Q}_{N}$) is computed to obtain $\widehat{\psi}_{pi}^{-1}$.
Crucially, a correct parametric model for the $q$-probabilities is often not
available and a natural alternative is to estimate them using flexible
methods, such as machine learning models, especially when the heterogeneous
capture probabilities depend on high-dimensional or continuous covariates.
However, such data-adaptive methods can cause $\widehat{\psi}_{pi}^{-1}$ to
exhibit slower than $\sqrt{N}$-convergence rates and to have an unknown
limiting distribution. The first limitation makes the plug-in estimator
inefficient, while the second prohibits us from making any inference about
$\widehat{\psi}_{pi}^{-1}$.

####  4.2.2 One-step estimator

To overcome these limitations of the plug-in estimator in semi- and non-
parametric settings, we leverage semiparametric efficiency theory and propose
the following one-step estimator using the Efficient Influence Function for
$\psi^{-1}$ derived in Lemma 1. Intuitively, the one-step estimator is a
debiased version of the plug-in estimator obtained from adding an estimate of
the mean of the EIF to the plug-in estimator (cf. van der Vaart, (2002);
Kennedy, (2022)). Concretely, the proposed one-step estimator
$\widehat{\psi}^{-1}_{os}$ takes the form:

$\widehat{\frac{1}{\psi_{os}}}=\mathbb{Q}_{N}\left(\left(\frac{1}{\widehat{\gamma}(X_{i})}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y_{i}=y)}{\widehat{q}_{y}(X_{i})}+1\right).$
(13)

As before, to construct $\widehat{\psi}_{os}^{-1}$ we first estimate the
nuisance $q$-probabilities from the available data, and then combine those
into capture probabilities. Using these estimates, we then compute the
empirical mean of the (uncentered) EIF to obtain $\widehat{\psi}_{os}^{-1}$.

######  Remark 3.

We assume the nuisance functions $q_{y}(X)$ are estimated using separate
samples from the one used to compute the empirical mean $\mathbb{Q}_{N}$ in
the proposed estimators for $\psi^{-1}$. These cross-fitting/sample-splitting
approaches control empirical process terms in our estimators, without
introducing potentially restrictive assumptions (Chernozhukov et al.,, 2018;
Kennedy,, 2022).

We now show that the proposed one-step estimator is nearly optimal,
approximately normal, and exhibits a doubly-robust-type second order
remainder. Specifically, we show that these properties hold for any number $N$
of observed units even when the $q$-probabilities are estimated flexibly using
machine learning models. We restrict ourselves to providing finite sample
guarantees because this is often the relevant setting in capture-recapture
problems where it is not possible to let the sample size $N\to\infty$, as we
aim to estimate a finite population of size $n\geq N$.

######  Theorem 2.

For any sample size $N$ the proposed one-step estimator
$\widehat{\psi}^{-1}_{os}$ satisfies

$\Big{|}(\widehat{\psi}^{-1}_{os}-\psi^{-1})-\mathbb{Q}_{N}(\phi)\Big{|}\leq\eta,$
(14)

for any $\eta>0$, with probability at least

$1-\left(\frac{1}{\eta^{2}}\right)\mathbb{E}_{\mathbb{Q}}\left[\widehat{R}_{2}^{2}+\frac{1}{N}\left\lVert\widehat{\phi}-\phi\right\rVert^{2}\right],$

where $\phi$ is the Efficient Influence Function (10) with second-order
reminder given by

$\begin{split}\widehat{R}_{2}=\mathbb{E}_{\mathbb{Q}}&\left[\left(\frac{1}{\widetilde{\gamma}(X)}-1\right)\left\\{\sum_{y_{i}\neq
y_{j}\neq
0}(-1)^{|y_{i}|+|y_{j}|}\frac{(q_{y_{i}}(X)-\widehat{q}_{y_{i}}(X))}{\widetilde{q}_{y_{i}}(X)}\frac{(q_{y_{j}}(X)-\widehat{q}_{y_{j}}(X))}{\widetilde{q}_{y_{j}}(X)}\right.\right.\\\
&+\left.\left.\sum_{|y|\neq 0,\leavevmode\nobreak\
even}\left(\frac{q_{y}(X)-\widehat{q}_{y}(X)}{\widetilde{q}_{y}(X)}\right)^{2}\right\\}\right],\end{split}$
(15)

where the sums are taken with respect to capture profiles of the form
$y=(y_{1},\dots,y_{J},0,\dots,0)$ and for some value $\widetilde{q}(X)$ of the
$q$-probabilities between $q(X)$ and $\widehat{q}(X)$.

The preceding theorem shows that, for any given sample size $N$, the proposed
one-step estimator $\widehat{\psi}^{-1}_{os}$ is nearly optimal, meaning that
it approximates the empirical average of the EIF, $\mathbb{Q}_{N}(\phi)$,
centered at the true inverse capture probability. This approximation holds for
any predetermined error tolerance $\eta>0$, with high probability depending on
the second-order remainder $\widehat{R}_{2}$ and the $L_{2}$-error in
estimating the influence function itself. Crucially, both components are
driven by the error in estimating the nuisance $q$-probabilities, but the
latter to a lesser extent by the normalization by $N$. Note that Theorem 2
also provides a conventional convergence in probability result, since for any
fixed $\eta$ the probability goes to $1$ whenever $\widehat{R}_{2}$ goes to
$0$ as the sample size $N$ increases.

Moreover, the previous result implies that the proposed estimator enjoys a
finite-sample “doubly-robust” property in the sense that $\widehat{R}_{2}^{2}$
can be small even if some of the nuisance $q$-probabilities are not well
estimated. In particular, our one-step estimator is nearly optimal with
probability depending on the square of the errors in estimating the $q_{y}$’s.
For example, $R_{2}$ will be small (i.e., of order $1/\sqrt{N}$) if each of
those estimation errors is of order $N^{-1/4}$, which can be achieved for many
flexible nonparametric estimators (Kennedy,, 2022). Thus, when a correct
parametric model for the $q$-probabilities is not available, the proposed one-
step estimator allows for the estimation of the nuisance functions using a
data-driven approach while retaining desirable statistical guarantees. This
result is detailed in the next Corollary.

######  Corollary 2.

Assume that $\widehat{q}_{y}(X)\geq\epsilon$ and $\left\lVert
q_{y}(X)-\widehat{q}_{y}(X)\right\rVert\leq\xi_{N}$, for all capture profiles
$y=(y_{1},\dots,y_{J},0,\dots,0)\neq 0$. Then,

$\mathbb{P}\left(\Big{|}(\widehat{\psi}^{-1}_{os}-\psi^{-1})-\mathbb{Q}_{N}(\phi)\Big{|}\leq\eta\right)\geq
1-\left(\frac{C}{\eta^{2}}\right)\left(\xi_{N}^{2}+\frac{1}{N}\right)$ (16)

for any sample size N and any $\eta>0$, where $C$ is a constant independent of
N.

In summary, we have that the proposed one-step estimator
$\widehat{\psi}^{-1}_{os}$ is close to a sample average of the EIF with high
probability when $\widehat{R}_{2}^{2}$ is small (nearly optimal), which can be
achieved even in large nonparametric models. Notably, this approximation to a
sample average of EIF values further implies that our estimator is
approximately Gaussian, and therefore we can derive approximately valid
confidence intervals using standard tools. This key result is formalized in
the following Theorem.

######  Theorem 3.

$\widehat{\psi}^{-1}_{os}-\psi^{-1}$ follows an approximately Gaussian
distribution, with the difference in cumulative distributions

$\Bigg{|}\mathbb{P}\left(\frac{\widehat{\psi}^{-1}_{os}-\psi^{-1}}{\widehat{\sigma}/\sqrt{N}}\leq
t\right)-\Phi(t)\Bigg{|}$

uniformly bounded above by

$\frac{C}{\sqrt{N}}\mathbb{E}_{\mathbb{Q}}\left[\frac{\rho}{\widetilde{\sigma}^{3}}\right]+\frac{1}{\sqrt{2\pi}}\left(\sqrt{N}\mathbb{E}_{\mathbb{Q}}\left[\frac{|\widehat{R}_{2}|}{\widetilde{\sigma}}\right]+|t|\mathbb{E}_{\mathbb{Q}}\left[\Big{|}\frac{\widehat{\sigma}}{\widetilde{\sigma}}-1\Big{|}\right]\right),$
(17)

where $\widehat{\sigma}^{2}=\widehat{var}(\widehat{\phi})$,
$\widetilde{\sigma}^{2}=var(\widehat{\phi}|Z^{n})$,
$\rho=\mathbb{E}[|\widehat{\phi}-\mathbb{Q}(\widehat{\phi})|^{3}|Z^{n}]$, and
$C<1/2$ is the Berry-Esseen constant.

The previous Berry-Esseen-type result implies that usual ($1-\alpha$)-level
confidence intervals of the form

$\widehat{CI}(\psi^{-1})=\left[\widehat{\psi}^{-1}_{os}\pm
z_{1-\alpha/2}\widehat{\sigma}/\sqrt{N}\right],$ (18)

can be used for inference on the inverse capture probability with nearly-valid
finite sample coverage guarantees, where
$\widehat{\sigma}^{2}=\widehat{var}(\widehat{\phi})$ is the unbiased empirical
variance of the estimated efficient influence function. Therefore, using the
proposed identification strategy, the derived one-step estimator for the
capture probability is not only nearly optimally efficient, but is also
approximately Normal in finite samples. The first property means that it
approximates a sample average of EIF values with high probability if the
nuisance errors are small in a doubly-robust fashion, while the latter
provides a general approach for inference on $\psi^{-1}$ and to construct
confidence intervals with non-asymptotic guarantees in a Berry-Esseen sense.

Finally, given that our identification assumption generalizes the two
predominant approaches in the capture-recapture literature, the derived one-
step estimator can be adapted to these settings while retaining its robust
statistical properties. For instance, when $J=K$, the proposed one-step
estimator can be used to estimate the inverse capture probability by
considering all capture profiles $y=(y_{1},\dots,y_{K})$ and corresponding
$q$-probabilities $\mathbb{Q}(Y=y\mid X)$. In the $J=2$ setting, our doubly-
robust estimator simplifies to

$\widehat{\frac{1}{\psi_{os}}}=\mathbb{Q}_{N}\left(\left(\frac{1}{\widehat{\gamma}(X)}-1\right)\left\\{\frac{Y_{1}(1-Y_{2})}{\widehat{q}_{10}(X)}+\frac{(1-Y_{1})Y_{2}}{\widehat{q}_{01}(X)}-\frac{Y_{1}Y_{2}}{\widehat{q}_{11}(X)}\right\\}+1\right),$
(19)

where $\widehat{q}_{y_{1}y_{2}}(X)$ are estimates of
$q_{y_{1}y_{2}}(X)=\mathbb{Q}(Y_{1}=y_{1},Y_{2}=y_{2},Y_{3}=0,\dots,Y_{K}=0\mid
X)$.

##  5 Sensitivity Analysis: Bounded Highest-Order Interaction

In this section, we relax the no highest-order interaction assumption and
propose a sensitivity analysis that allows for arbitrary dependence between
the $J$ capture episodes in our conditional log-linear model among the
unobserved (5). Specifically, relaxing the identifying assumption leads to a
partial identification of $\psi^{-1}$, resulting in an interval that contains
the true target parameter as a function of a parameter bounding the highest-
order interaction coefficient. We then leverage the one-step estimator derived
for the no highest-order interaction case to propose similar one-step
estimators for the upper and lower bounds of interest. Our estimators are
shown to be nearly optimal and approximately normal, even when flexibly
estimating the $q$-probabilities, under an additional margin condition (26).
Finally, we derive approximately valid confidence intervals that can be
combined to bound the true capture probability target parameter with high
confidence. The proposed sensitivity analysis reveals how deviations from the
no highest-order interaction and positivity assumptions impact the resulting
population size estimates, allowing analysts to incorporate domain-specific
knowledge into the identification and estimation steps more transparently and
rigorously.

###  5.1 Partial identification

We relax the identifying assumption of no highest-order interaction in our
saturated conditional log linear model among the unobserved (5) and replace it
by the assumption that the highest-order interaction coefficient,
$\alpha_{1}(X)$, is uniformly bounded by a parameter $\delta\geq 0$.
Mathematically, in a log-linear model for $J$ lists conditional on covariates
and on being unobserved by the remaining $K-J$ lists, the relaxed
identification approach assumes that

$\Big{|}\alpha_{1}(X)\Big{|}=\Bigg{|}\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{J+|y|}\log(q_{y}(X)\gamma(X))+(-1)^{J}\log(1-\gamma(X))\Bigg{|}\leq\delta,$
(20)

where $|y|=\sum_{j=1}^{J}y_{j}$ is the number of positive lists in the capture
profile $y=(y_{1},\dots,y_{J},0,\dots,0)$. Note that $\delta=0$ recovers the
no highest-order interaction assumption and the previously proposed estimator
would still be valid. On the other hand, setting $\delta>0$ is indeed a
relaxation of our assumption allowing for a non-negligible highest-order
interaction between the $J$ capture episodes. The Appendix provides an
interpretation of the bounded highest-order interaction assumption in terms of
conditional odds ratios.

By setting a bound on the highest-order interaction parameter, this approach
yields bounds on the inverse capture probability $\psi^{-1}$. Essentially,
this means that while we can no longer identify the exact value of the capture
probability, we achieve a partial identification of $\psi$ using an interval
depending on the bound $\delta$. From the relaxed identifying assumption
$|\alpha_{1}(X)|\leq\delta$, we obtain a partial identification expression
detailed in the following Proposition 2.

######  Proposition 2.

(Partial-Identification Result) Assume positivity and bounded highest-order
interaction coefficient, $|\alpha_{1}(X)|\leq\delta$, in a log-linear model
for $J$ lists, conditional on covariates and on not being captured by the
other $K-J$ lists. Then, we can partially identify the inverse capture
probability $\psi^{-1}$ using the observed data with the lower and upper
bounds

$\frac{1}{\psi_{\ell}}=\mathbb{E}_{\mathbb{Q}}\left[\min\left\\{1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))-\delta\right),\frac{1}{\epsilon}\right\\}\right],$
(21)
$\frac{1}{\psi_{u}}=\mathbb{E}_{\mathbb{Q}}\left[\min\left\\{1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))+\delta\right),\frac{1}{\epsilon}\right\\}\right],$
(22)

such that $\psi^{-1}\in\left[\psi_{\ell}^{-1},\psi_{u}^{-1}\right]$. As
before, $y$ denotes the complete $K$-dimensional capture profile
$y=(y_{1},\dots,y_{J},0,\dots,0)$.

Intuitively, the proposed bounds are derived from a partial identification
expression for the conditional capture probability

$\frac{1}{\gamma(X)}\in\left[1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\pm\delta\right)\right],$

with an additional constraint on the parameter $\delta$ to ensure that
$\gamma(X)\in[\epsilon,1]$, under the positivity assumption on the conditional
capture probability. This parametrization of the identification assumptions
allows to examine how variations in the $\delta$ and $\epsilon$ parameters can
affect estimates of the total population size. This provides a straightforward
and flexible approach to integrating domain-specific insights into the
identification and estimation procedures.

###  5.2 One-step estimator

We now turn to the problem of estimating the proposed bounds that partially
identify the capture probability under the bounded highest-order interaction
condition. With a slight abuse of notation, we use $\psi_{\delta}^{-1}$ to
refer to either the lower or upper bound derived in Proposition 2. By their
symmetry, the results presented in this section apply to both
$\psi_{\ell}^{-1}$ and $\psi_{u}^{-1}$, so we present them succinctly as
applying to the generic bound $\psi_{\delta}^{-1}$ defined by

$\frac{1}{\psi_{\delta}}=\mathbb{E}_{\mathbb{Q}}\left[\min\left\\{\frac{1}{\gamma_{\delta}(X)},\frac{1}{\epsilon}\right\\}\right]=\mathbb{E}_{\mathbb{Q}}\left[\min\left\\{1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\pm\delta\right),\frac{1}{\epsilon}\right\\}\right],$

where $\gamma_{\delta}(X)^{-1}=1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\pm\delta\right)$ denotes the conditional capture
probability associated with the corresponding bound of interest and
$y=(y_{1},\dots,y_{J},0,\dots,0)$.

Our previous strategy was to use the Efficient Influence Function for
$\psi^{-1}$ to construct a one-step estimator for the capture probability
under the assumption of no highest-order interaction. However, since the min
function appearing in the identifying expressions in Proposition 2 is not
smooth, we cannot directly follow the same approach and derive the EIF to
propose estimators for the target bounds $\psi^{-1}_{\ell}$ and
$\psi^{-1}_{u}$. To overcome this, note that we can rewrite our generic bound
as

$\begin{split}\frac{1}{\psi_{\delta}}&=\mathbb{E}_{\mathbb{Q}}\left[\frac{1}{\gamma_{\delta}(X)}\mathbbm{1}\left(\frac{1}{\gamma_{\delta}(X)}\leq\frac{1}{\epsilon}\right)\right]+\frac{1}{\epsilon}\mathbb{E}_{\mathbb{Q}}\left[\mathbbm{1}\left(\frac{1}{\gamma_{\delta}(X)}>\frac{1}{\epsilon}\right)\right]\\\
&=\mathbb{E}_{\mathbb{Q}}\left[\left(\frac{1}{\gamma_{\delta}(X)}-\frac{1}{\epsilon}\right)\mathbbm{1}\left(\frac{1}{\gamma_{\delta}(X)}-\frac{1}{\epsilon}\leq
0\right)\right]+\frac{1}{\epsilon},\end{split}$ (23)

which suggest to use a flexible estimator for the first term involving the
conditional capture probability, and a plug-in estimator for the indicator
function to deal with its non-smoothness.

Concretely, we propose the following one-step estimator for
$\psi_{\delta}^{-1}$:

$\widehat{\frac{1}{\psi_{\delta}}}=\mathbb{Q}_{N}\left(\left(\widehat{\varphi}_{\delta}(X,Y)-\frac{1}{\epsilon}\right)\mathbbm{1}\left(\frac{1}{\widehat{\gamma}_{\delta}(X)}-\frac{1}{\epsilon}\leq
0\right)\right)+\frac{1}{\epsilon},$ (24)

where

$\varphi_{\delta}(X,Y)=\left(\frac{1}{\gamma_{\delta}(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}+1,$ (25)

is the (uncentered) EIF for $\mathbb{E}[\gamma_{\delta}(X)^{-1}]$. Note that
the EIF (25) parallels that in Lemma 1 for $\psi^{-1}$ in the no highest-order
interaction assumption setting. This is a consequence of the analogous
conditional capture probabilities $\gamma(X)$ and $\gamma_{\delta}(X)$
involved in the derivation of the corresponding EIFs. This similarity also
endows the proposed one-step estimator $\widehat{\psi}_{\delta}^{-1}$ with
statistical properties akin to those of the one-step estimator derived for the
scenario with no highest-order interaction.

Specifically, in the following we show that the proposed estimator is nearly
optimal and approximately Gaussian with finite sample guarantees. Intuitively,
the proposed one-step estimator behaves like an infeasible estimator that has
full knowledge of the conditional capture probability $\gamma_{\delta}(X)$ and
uses it to compute the indicator function in (23). For such “oracle” behavior
to be met, a sufficient condition is the so-called margin condition, which has
been used to deal with non-smooth estimands in causal inference (Luedtke and
Van Der Laan,, 2016; Kennedy et al.,, 2020; Levis et al.,, 2023) and
classification settings (Audibert and Tsybakov,, 2007). Formally, we assume
that there exists $\beta>0$ and a nonnegative constant $C$, such that for any
$t\geq 0$

$\mathbb{P}\left(\Bigg{|}\frac{1}{\gamma_{\delta}(X)}-\frac{1}{\epsilon}\Bigg{|}\leq
t\right)\leq Ct^{\beta}.$ (26)

The margin condition controls the probability that the conditional capture
probability is near the bound $\frac{1}{\epsilon}$ making the indicator
function non-differentiable, with $\beta$ controlling how quickly this
probability vanishes and becomes negligible for the correct estimation of
$\psi_{\delta}^{-1}$. This key result is formalized in the following Theorem.

######  Theorem 4.

Under the margin condition for $\beta>0$ and nonnegative $C$, the one-step
estimator $\widehat{\psi}^{-1}_{\delta}$ satisfies

$\mathbb{P}\left(\Big{|}(\widehat{\psi}^{-1}_{\delta}-\psi^{-1}_{\delta})-\mathbb{Q}_{N}(\phi_{\delta})\Big{|}\leq\eta\right)\geq
1-\frac{1}{\eta^{2}}\mathbb{E}\left[\widehat{R}_{2,\delta}^{2}+\frac{1}{N}\left\lVert\widehat{\phi}_{\delta}-\phi_{\delta}\right\rVert^{2}\right]$
(27)

for any sample size $N$ and $\eta>0$, where

$\phi_{\delta}(X,Y)=\left(\varphi_{\delta}(X,Y)-\frac{1}{\epsilon}\right)\mathbbm{1}\left(\frac{1}{\gamma_{\delta}(X)}-\frac{1}{\epsilon}\leq
0\right)+\frac{1}{\epsilon}-\frac{1}{\psi_{\delta}}$ (28)

and the reminder term is bounded by

$\widehat{R}_{2,\delta}\leq\sum_{y_{i}\neq y_{j}\neq 0}\left\lVert
q_{y_{i}}-\widehat{q}_{y_{i}}\right\rVert\left\lVert
q_{y_{j}}-\widehat{q}_{y_{j}}\right\rVert+\sum_{|y|\neq 0,\leavevmode\nobreak\
even}\left\lVert
q_{y}-\widehat{q}_{y}\right\rVert^{2}+C\left\lVert\widehat{\gamma}-\gamma\right\rVert_{\infty}^{1+\beta},$
(29)

where the sums are taken with respect to capture profiles
$y=(y_{1},\dots,y_{J},0,\dots,0)$.

Similarly to the setting assuming no highest-order interaction, this result
shows that the proposed one-step estimator for the generic bound
$\widehat{\psi}_{\delta}^{-1}$ is close to the empirical average of its
corresponding EIF (28). This holds for any given error tolerance $\eta>0$ and
sample size $N$ with high probability. Crucially, this probability now depends
on the margin condition, in addition to the second-order remainder
$\widehat{R}_{2}$ and the $L_{2}$-error in estimating the EIF. Therefore, the
proposed estimator enjoys a similar doubly-robust-type property if $\beta\geq
1$, depending on the square of the errors in estimating the $q_{y}$’s.
Analogous to Corollary 2, a sufficient condition for
$R_{2}=o_{\mathbb{P}}(1/\sqrt{N})$ is that each of the estimation errors
$\left\lVert\widehat{q}_{y}(X)-q_{y}(X)\right\rVert$ be of order $N^{-1/4}$,
which can be obtained for many off-the-shelf machine learning models
(Kennedy,, 2022). This allows to estimate the nuisance functions in a data-
driven way without sacrificing desired statistical properties. Note that the
previous theorem also gives standard conditions for
$\widehat{\psi}_{\delta}^{-1}$ to be $\sqrt{N}$-consistent and asymptotically
normal.

Moreover, Theorem 4 provides that the proposed one-step estimator
$\widehat{\psi}^{-1}_{\delta}$ is not only close to a sample mean of EIF
values, but is also approximately Gaussian in finite samples. We detail this
result in the following Theorem and use such property to derive approximately
valid confidence intervals for the inverse capture probability.

######  Theorem 5.

$\widehat{\psi}^{-1}_{\delta}-\psi^{-1}_{\delta}$ follows an approximately
Gaussian distribution, with the difference in cumulative distributions
uniformly bounded above by

$\Bigg{|}\mathbb{P}\left(\frac{\widehat{\psi}^{-1}_{\delta}-\psi^{-1}_{\delta}}{\widehat{\sigma}_{\delta}/\sqrt{N}}\leq
t\right)-\Phi(t)\Bigg{|}\leq\frac{1}{\sqrt{2\pi}}\left(\sqrt{N}\mathbb{E}_{\mathbb{Q}}\left[\frac{|\widehat{R}_{2,\delta}|}{\widetilde{\sigma}_{\delta}}\right]+|t|\mathbb{E}_{\mathbb{Q}}\left[\Big{|}\frac{\widehat{\sigma}_{\delta}}{\widetilde{\sigma}_{\delta}}-1\Big{|}\right]\right)+\frac{C}{\sqrt{N}}\mathbb{E}_{\mathbb{Q}}\left[\frac{\rho_{\delta}}{\widetilde{\sigma}^{3}_{\delta}}\right],$

where $\widehat{\sigma}^{2}_{\delta}=\widehat{var}(\widehat{\phi}_{\delta})$,
$\widetilde{\sigma}^{2}_{\delta}=var(\widehat{\phi}_{\delta}|Z^{n})$,
$\rho_{\delta}=\mathbb{E}[|\widehat{\phi}_{\delta}-\mathbb{Q}(\widehat{\phi}_{\delta})|^{3}|Z^{n}]$,
$C<1/2$ is the Berry-Esseen constant, and $\widehat{R}_{2,\delta}$ is bounded
by (29).

The previous theorem implies that the following are approximately valid
$(1-\alpha)$-level confidence intervals for the bounds on the capture
probability $\psi^{-1}\in[\psi^{-1}_{\ell},\psi^{-1}_{u}]$:

$\widehat{CI}(\psi^{-1}_{\ell})=\left[\widehat{\psi}^{-1}_{\ell}\pm
z_{1-\alpha/2}\widehat{\sigma}_{\ell}/\sqrt{N}\right],\quad\quad\widehat{CI}(\psi^{-1}_{u})=\left[\widehat{\psi}^{-1}_{u}\pm
z_{1-\alpha/2}\widehat{\sigma}_{u}/\sqrt{N}\right],$ (30)

where $\widehat{\psi}_{\ell}^{-1}$ and $\widehat{\psi}_{u}^{-1}$ are the
proposed one-step estimators (24), and
$\widehat{\sigma}^{2}_{\ell}=\widehat{var}(\widehat{\phi_{\ell}})$,
$\widehat{\sigma}^{2}_{u}=\widehat{var}(\widehat{\phi_{u}})$, are the unbiased
empirical variances of the corresponding EIFs estimated from expression (28).
Therefore, an asymptotically valid ($1-\alpha$)-level confidence interval for
$\psi^{-1}$ can be derived by combining the lower limit of the CI for
$\widehat{\psi}^{-1}_{\ell}$ and the upper limit of the CI for
$\widehat{\psi}^{-1}_{u}$:

$\psi^{-1}\in\left[\widehat{\psi}_{\ell}^{-1}-z_{1-\alpha/2}\widehat{\sigma}_{\ell}/\sqrt{N},\quad\widehat{\psi}_{u}^{-1}+z_{1-\alpha/2}\widehat{\sigma}_{u}/\sqrt{N}\right],$
(31)

where $\\{\widehat{\sigma}_{\ell},\widehat{\sigma}_{u}\\}$ are defined as
before. Alternatively, for more precise inference, the CIs in (30) can be
combined into a confidence interval for the partially identified parameter
$\psi^{-1}$ using the methods proposed in Imbens and Manski, (2004).

In summary, the proposed sensitivity analysis accommodates a more general
dependence between capture episodes in a log-linear model at the expense of
only achieving partial identification of the capture probability. We argue
that the proposed bounds are still relevant to population size estimation
because they provide an informative interval containing the true capture
probability, relying on a weaker identification assumption than the no
highest-order interaction. Moreover, the $\delta$ and $\epsilon$ parameters
allow analysts to incorporate domain knowledge about possible bounds on the
highest-order interaction and positivity assumptions, and facilitates the
examination of how variations in these bounds can affect estimates of the
total population size. As before, our approach can also be used for
sensitivity analysis in standard log linear models ($J=K$) and conditional
independence settings ($J=2$).

##  6 Inference for total population size

So far, we have focused on constructing flexible and robust estimators for the
the inverse capture probability $\psi^{-1}$ with finite sample guarantees. In
this section, we shift our attention to estimating the total population size
$n$, which is typically the ultimate goal in a capture-recapture setting. To
achieve this, we use $\widehat{n}=N/\widehat{\psi}$ as natural estimator for
$n$, given an estimator for the capture probability $\widehat{\psi}$, where
$N\sim\text{Bin}(n,\psi)$ is the observed sample size of uniquely captured
subjects. Therefore, in the following we leverage the one-step estimators
proposed above and provide estimators and approximately valid confidence
intervals for $n$.

###  6.1 Population size under no highest-order interaction

Concretely, in the no highest-order interaction scenario, we propose to
estimate the total population size $n$ via
$\widehat{n}_{os}=N/\widehat{\psi}_{os}$, where $\widehat{\psi}_{os}^{-1}$ is
the one-step estimator proposed in (13). We then derive a confidence interval
for the total population size $n$ using Theorem 4 in Das et al., (2023),
provided that $\widehat{\psi}_{os}^{-1}$ is approximated by a sample average.
We adapt this result to our setting in the following Theorem 6.

######  Theorem 6.

Consider a log-linear model for $J$ lists, conditional on covariates and on
not being captured by the remaining $K-J$ lists. Assuming no highest-order
interaction term, the following is an approximately $(1-\alpha)$ confidence
interval for the population size $n$ centered at
$\widehat{n}_{os}=N/\widehat{\psi}_{os}$

$\widehat{CI}(n)=\left[\widehat{n}_{os}\pm
z_{\alpha/2}\sqrt{\widehat{n}_{os}\left(\widehat{\psi}_{os}\widehat{\sigma}^{2}+\frac{1-\widehat{\psi}_{os}}{\widehat{\psi}_{os}}\right)}\right],$
(32)

where $\widehat{\sigma}^{2}=\widehat{var}(\widehat{\phi})$ is the unbiased
empirical variance of the estimated EIF (10) and $\widehat{\psi}_{os}^{-1}$ is
the proposed one-step estimator (13).

###  6.2 Bounds on the population size under partial identification

The relaxed assumption of a bounded highest-order interaction coefficient in
the saturated conditional log-linear model among the unobserved (5) implies
lower and upper bounds on the total population size $n$ given by

$n_{\ell}=\frac{N}{\psi_{\ell}},\quad\quad\quad n_{u}=\frac{N}{\psi_{u}},$
(33)

such that $n\in[n_{\ell},n_{u}]$. Analogously as before, we estimate such
bounds using the derived one-step estimators (24) for $\psi_{\ell}^{-1}$ and
$\psi_{u}^{-1}$, respectively. Furthermore, in the following Theorem, we
leverage that our one-step estimators can be approximated by sample averages
and derive approximately valid confidence intervals for the estimated bounds
on the population size $n$.

######  Theorem 7.

Consider a log-linear model for $J$ lists, conditional on covariates and on
not being captured by the remaining $K-J$ lists. The following are approximate
$(1-\alpha)$ confidence intervals for the bounds on the population size
$n\in[n_{\ell},n_{u}]$, assuming a bounded highest-order interaction term:

$\widehat{CI}(n_{\ell})=\left[\frac{N}{\widehat{\psi}_{\ell}}\pm
z_{\alpha/2}\sqrt{N\left(\widehat{\sigma}^{2}_{\ell}+\frac{1-\widehat{\psi}_{\ell}}{\widehat{\psi}_{\ell}^{2}}\right)}\right],$
(34) $\widehat{CI}(n_{u})=\left[\frac{N}{\widehat{\psi}_{u}}\pm
z_{\alpha/2}\sqrt{N\left(\widehat{\sigma}^{2}_{u}+\frac{1-\widehat{\psi}_{u}}{\widehat{\psi}_{u}^{2}}\right)}\right],$
(35)

where $\widehat{\psi}_{\ell}$ and $\widehat{\psi}_{u}$ are one-step estimators
(24) for the corresponding lower and upper bounds in Proposition 2, and
$\widehat{\sigma}^{2}_{\ell}=\widehat{var}(\widehat{\phi_{\ell}})$,
$\widehat{\sigma}^{2}_{u}=\widehat{var}(\widehat{\phi_{u}})$, are the unbiased
empirical variances of the estimated influence functions derived in (28).

Finally, the previous theorem implies that

$\widehat{CI}(n)=\left[\frac{N}{\widehat{\psi}_{\ell}}-z_{\alpha/2}\sqrt{N\left(\widehat{\sigma}^{2}_{\ell}+\frac{1-\widehat{\psi}_{\ell}}{\widehat{\psi}_{\ell}^{2}}\right)},\quad\frac{N}{\widehat{\psi}_{u}}+z_{\alpha/2}\sqrt{N\left(\widehat{\sigma}^{2}_{u}+\frac{1-\widehat{\psi}_{u}}{\widehat{\psi}_{u}^{2}}\right)}\right]$
(36)

is an approximate $(1-\alpha)$ confidence interval for the total population
size $n$. Similarly as before, these bound estimates can be combined for more
precise inference following Imbens and Manski, (2004) methodology.

Note that $\delta=0$ recovers the CI (32) for the corresponding point-
identified population size $n$. In both cases, the coverage guarantee of the
derived confidence intervals is mainly driven by the error term through
$\mathbb{E}[|\widehat{R}_{2}^{2}|]$, under standard boundedness assumptions.
This highlights the importance of efficient estimation of
$\\{\psi_{\ell},\psi_{u}\\}$ (or $\psi$ under no highest-order interaction)
and the approximation of $\\{\widehat{\psi}_{\ell},\widehat{\psi}_{u}\\}$
($\widehat{\psi}$) by sample averages. Crucially, our proposed one-step
estimators can attain an error $\widehat{R}_{2}^{2}$ of order $1/\sqrt{N}$ as
shown above, even when the $q$-probabilities are estimated flexibly using
machine learning models.

##  7 Data Analysis: Synthetic and Real Data from the Peruvian Conflict

In the following, we demonstrate the advantages of the proposed methods using
both simulated and real capture-recapture data from the Peruvian internal
armed conflict (1980-2000). We compare our one-step estimator with a simple,
but typically suboptimal, plug-in estimator and show that our method dominates
in several regimes. Moving forward, we use data from the Peruvian Truth
Commission to estimate the number of casualties and highlight the role of the
identification strategy adopted on the resulting estimated population size.
Finally, we apply the proposed sensitivity analysis to estimate bounds on the
population size for different values of the highest-order interaction term
bound $\delta$, providing intervals that partially identify the number of
casualties and illustrating the sensitivity of the estimates to the
identification assumptions.

###  7.1 Synthetic capture-recapture data with 3 lists and no highest-order
interaction assumption

To empirically demonstrate the properties of the proposed one-step estimator
(13), we begin by simulating a scenario with latent heterogeneous capture
probabilities $q_{y}$ drawn from a uniform distribution with an empirical
lower bound of approximately $0.05$. We then simulate the estimation step of
the $q$-probabilities by setting

$\widehat{q}_{y}=expit(logit(q_{y})+\epsilon_{y}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \epsilon_{y}\sim\mathcal{N}(bn^{-\alpha},n^{-2\alpha}),$

where the $\alpha$ term gives the rate of convergence of the estimated
$\widehat{q}_{y}$ to the true latent $q_{y}$ (because the RMSE scales as
$n^{-\alpha}$), and $b$ controls the bias of the estimates. We vary $\alpha$
and $b$ to analyze the behavior of the resulting $\widehat{\psi}^{-1}$ for
different biases and rates of convergence of the nuisance functions,
simulating different scenarios where an analyst is able to estimate the
$q$-probabilities with varying precision. All the simulations are done $1,000$
times for a true $n=10,000$ and $\psi^{-1}=1.41$ ($\psi=0.7$), yielding an
observed $N\approx 7,000$.

From our experiments we can observe that, when the $q$-probabilities are
estimated at better rates closer to $\alpha=0.5$ (the parametric rate), both
the plug-in and the one-step estimators perform well with decreasing bias and
MSE, and coverage at the expected nominal level (Figures 1 and 2). The
advantages of the proposed one-step estimator are larger in regimes of slow
rate of convergence, especially for $\widehat{q}_{y}$ estimates with large
noise $b=10$ (Fig. 2). For instance, for $\alpha$ smaller than the non-
parametric rate of $\alpha=0.25$, the plug-in estimator
$\widehat{\psi}_{pi}^{-1}$ has worst bias, MSE, and coverage than our one-step
estimator $\widehat{\psi}_{os}^{-1}$. Note that the plug-in estimator does not
have a well-defined variance formula, thus we use the estimated variance from
the one-step estimator to construct corresponding confidence intervals.
Finally, it should also be noticed that in the large noise regime the one-step
estimator suffers from large variance which dominates the behavior of its MSE.

(a) Bias $|\psi^{-1}-\widehat{\psi}^{-1}|$.

(b) MSE $(\psi^{-1}-\widehat{\psi}^{-1})^{2}$.

(c) Empirical coverage of $\psi^{-1}$.

Figure 1: Results for $\psi^{-1}$ using synthetic data from three lists and
small noise on $\widehat{q}_{y}$ ($b=1)$. True $n=10,000$ and
$\psi^{-1}=1.41$. Average over $1,000$ simulations.

(a) Bias $|\psi^{-1}-\widehat{\psi}^{-1}|$.

(b) MSE $(\psi^{-1}-\widehat{\psi}^{-1})^{2}$.

(c) Empirical coverage of $\psi^{-1}$.

Figure 2: Results for $\psi^{-1}$ on synthetic data from three lists and large
noise on $\widehat{q}_{y}$ ($b=10)$. True $n=10,000$ and $\psi^{-1}=1.41$.
Average over $1,000$ simulations.

Overall, the proposed one-step estimator for the (inverse) capture probability
$\psi^{-1}$ exhibits a better performance than the plug-in estimator,
especially when the nuisance $q$-probabilities cannot be well estimated. This
is often the case in the absence of a reliable parametric model for $q_{y}$
and in complex capture-recapture settings where the capture probabilities vary
with respect to continuous or high-dimensional variables. In these scenarios,
the proposed one-step estimator proves useful allowing for flexible estimation
of the $q$-probabilities and robust estimation of $\psi$.

###  7.2 Application to the Peruvian internal armed conflict

Building on the theoretical and empirical advantages of our one-step
estimator, we apply the proposed methodology to estimate the size of the
victim population in the Peruvian internal armed conflict. This capture-
recapture dataset was used by the Peruvian Truth and Reconciliation Commission
to estimate the number of victims of killings and disappearances during the
conflict between 1980 and 2000 (Ball et al.,, 2003). The dataset consists of
seven lists coming from three different sources of information documenting the
conflict: (i) the Public Defender’s Office (CP), (ii) the Truth and
Reconciliation Commission (CVR), and (iii) five human rights and non-
governmental organizations (ODH). The data includes information on $24,692$
uniquely identified victims of killings or disappearances, their capture
profiles, and complementary demographic variables (Ball et al.,, 2003;
Rendon,, 2019).

The original report by the Peruvian Truth and Reconciliation Commission
estimated approximately 69,000 victims during the conflict using log-linear
models and the no highest-order interaction assumption (Ball et al.,, 2003).
However, recent academic articles have estimated the victim population size in
48,000 using an alternative methodological approach (Rendon,, 2019), between
58,000 and 66,000 using Bayesian latent class models (Manrique-Vallier et
al.,, 2019), and in 69,000 victims using a two-list doubly robust estimator
under conditional independence (Das et al.,, 2023). Following Das et al.,
(2023), we estimate the $q$-probabilities from the available demographic
variables of sex, age and geographical location using the Random Forest
methodology (Breiman,, 2001).

First, we follow the approach in previous literature and consolidate the five
lists from different NGOs into a single list (ODH in (Ball et al.,, 2003)) and
use a standard conditional log-linear model on the resulting $K=3$ lists: CP,
CVR, and ODH. We analyze how the estimated victim population size varies for
different positivity lower bounds $\epsilon$ truncating the estimated
$\widehat{q}_{y}$’s in Figure 3(a). Without imposing a lower bound on the
estimated $q$-probabilities (i.e., no positivity constraints), we obtain
uninformative impossible estimates of $\widehat{n}=30^{\prime}715,584$ ($95\%$
CI: $[24,692,79^{\prime}973,134]$), larger than the Peruvian population in
1980. In this setting, the plug-in estimator returns $\widehat{n}=\infty$.
Notably, by imposing a $1\%$ lower bound on all the $q$-probabilities, the
proposed one-step estimator estimates $\widehat{n}=121,733$ ($95\%$ CI:
$[24,692,456,336]$),333The provided CI is not symmetric around $\widehat{n}$
since we set the lower bound to be the observed $N$ whenever it is estimated
to be a lower value. while using $10\%$ as the lower bound yields
$\widehat{n}=52,215$ ($95\%$ CI: $[51,073,53,356]$, not shown). Using a lower
bound of $4\%$, the results are in line with previous estimates of the victims
of the Peruvian internal conflict: $\widehat{n}=69,970$ ($95\%$ CI:
$[64,528,75,413]$) (Ball et al.,, 2003; Manrique-Vallier et al.,, 2019; Das et
al.,, 2023). This behavior highlights the crucial role of the positivity
assumption of the nuisance functions $q_{y}$’s on the estimated population
size, which is often overlooked in empirical applications of capture-recapture
methods.

We demonstrate the advantages of our methodology by employing a conditional
log-linear model on the DP and CVR lists, conditional on not being captured by
the remaining five ODH lists. In this setting, the no highest-order
interaction identifying assumption is equivalent to assuming conditional
independence between these two lists and therefore we use the estimator (19)
from Corollary 1. The identification assumption that the DP and CVR lists are
conditionally independent has been previously used in the literature to
estimate the number of casualties given their scope and capture strategies
(Das et al.,, 2023). Moreover, we relax this assumption by conditioning on
being unobserved by the other five lists from the ODH organizations. The
results, show in Figure 3(b), demonstrate that this approach is less
sensitivity to the positivity bound giving informative results for $\epsilon$
as small as $1\%$: $\widehat{n}=67,177$ ($95\%$ CI: $[50,287,84,067]$).
Notably, the point estimates quickly concentrates around approximately
$68,000$ victims in line with the previous literature on the Peruvian internal
armed conflict. In contrast, the plug-in estimator in this case returns
$\widehat{n}=85,350$ ($95\%$ CI: $[68,446,102,253]$) for $\epsilon=0.01$.

(a) Three lists (CP, CVR, ODH) with no highest-order interaction conditional
on covariates.

(b) Two independent lists (CVR, DP) conditional on being unobserved by ODH and
covariates.

Figure 3: Estimated number of casualties in the Peruvian internal conflict for
different identification strategies vs. positivity lower bound on the
estimated $q$-probabilities.

Next, we use the proposed sensitivity analysis framework to estimate bounds
that partially identify the population size and analyze how sensitive the
results are to deviations from the no highest-order interaction assumption. We
present the estimated lower and upper bounds, and the resulting interval for
$n$, when we vary the bound $\delta$ on the highest-order interaction term in
Figure 4. As expected, setting $\delta=0$ recovers the no highest-order
interaction case, while larger values of $\delta$ yield wider bounds, as they
allow for more complex dependence structures between the different capture
episodes making identification more difficult. Figure 4(a) depicts the
sensitivity analysis for the conditional log-linear model with three lists.
For instance, using $\delta=0.2$ the resulting estimated bounds are
$\widehat{n}_{\ell}=63,991$ ($95\%$ CI: $[58,593,69,389]$) and
$\widehat{n}_{u}=74,121$ ($95\%$ CI: $[68,324,79,918]$). That is, using a
conditional log-linear model with highest-order interaction bounded by
$\delta=0.2$, we estimate the population size to be between $58,593$ and
$79,918$. On the other hand, Figure 4(b) shows that the conditional log-linear
model among the unobserved using two lists (CVR and DP) is less sensitive to
the no highest-order interaction assumption, making the estimated population
size more robust to deviations from the conditional independence assumption.
With $\delta=0.2$, the number of casualties is estimated between $53,027$ and
$72,622$, with the lower and upper bounds $\widehat{n}_{\ell}=58,950$ ($95\%$
CI: $[53,027,64,872]$) and $\widehat{n}_{u}=66,348$ ($95\%$ CI:
$[60,074,72,622]$). Note that for this analysis we use fixed positivity bounds
$\epsilon=4\%$ and $\epsilon=2\%$ for the conditional log-linear models with
three and two lists, respectively, given our previous results in the no
highest-order interaction setting.

(a) Log-linear model with three lists (CVR, DP and ODH) conditional on
covariates ($\epsilon=4\%$).

(b) Two lists (CVR and DP) conditional on being unobserved by ODH and
covariates ($\epsilon=2\%$).

Figure 4: Estimated bounds for the number of casualties in the Peruvian
internal conflict for different values of the bound $\delta$ on the highest-
order interaction in a conditional log-linear model.

In summary, our empirical analysis of real data from the Peruvian internal
armed conflict demonstrates the key role of the identification assumptions in
causality estimation. In traditional log-linear models using smaller values
for the positivity assumption makes the estimation more difficult (as
theoretically predicted by the efficiency bound in Theorem 1), resulting in
very wide confidence intervals. However, imposing a large positivity bound may
bias the estimation in complicated ways. On the other hand, using larger
values for the highest-order interaction bound accounts for more complex
dependence structures between the different lists, but makes identification
more difficult and yields wider bounds. Crucially, our framework allows us to
study how these bounds affect the estimated population size and to incorporate
domain knowledge of plausible values. Using the proposed identification
strategy on the two largest lists in the Peruvian internal armed conflict DP
and CVR, conditional on not being captured by the remaining five ODH lists, we
estimate approximately $68,000$ victims during the Peruvian internal armed
conflict in line with previous literature.

##  8 Conclusions

In this paper, we introduced a general and flexible framework for estimating
the total size of a closed population from capture-recapture data, addressing
key challenges such as heterogeneous capture probabilities, low overlap
between lists, and complex dependencies across capture-recapture episodes. We
derived a one-step estimator for the capture probability within a novel
conditional log-linear model on the unobserved, accommodating heterogeneous
and dependent capture probabilities. This approach extends standard
assumptions of no highest-order interaction and conditional independence,
commonly used in previous literature. Our analysis demonstrates that our
estimators are nearly optimal, approximating a sample average of efficient
influence function values with high probability when nuisance errors are small
in a doubly-robust fashion. Additionally, our estimators are shown to be
approximately normal in a Berry-Esseen sense which yield confidence intervals
with non-asymptotic guarantees. These properties hold for any number of
observed units, even when the $q$-probabilities are flexibly estimated using
machine learning models in a data-driven way without compromising statistical
guarantees. Furthermore, we extended our results to estimate and perform
inference on the total population size $n$.

Notably, our identification strategy is weaker than assuming a log-linear
model for all $K$ lists with no highest-order interaction, requiring only that
$J\leq K$ lists are not arbitrarily dependent within the unobserved
population, conditional on covariates. This approach is both interpretable and
actionable: to make the identification of the capture probability easier, an
analyst can measure more relevant covariates or collect additional lists.
Moreover, it only requires positivity for $2^{J}-1$ of the $q$-probabilities,
which is substantially fewer than the $2^{K}-1$ possible capture profiles.
Finally, our sensitivity analysis framework allows easy incorporation of
domain knowledge about identification assumption parameters and analysis of
their effects on the estimated population size.

For instance, in our analysis of simulated and real data from the Peruvian
Truth Commission, we highlighted the critical role of the identifying approach
used in population size estimation, emphasizing the need for a more flexible
and rigorous methodology. Standard log-linear models were found to be highly
sensitive to the identification assumptions: smaller values for the positivity
parameter made the estimation more difficult, resulting in very wide
confidence intervals, while larger values for the highest-order interaction
bound made identification more difficult and yielded wider bounds. In
contrast, our proposed identification strategy, using a conditional log-linear
model on the two largest lists ($J=2$) conditioned on not being captured by
the remaining five ODH lists and covariates, provided a stable and transparent
estimation of the victim population, consistent with previous literature.

Our methodology represents a significant contribution to the population size
estimation literature, offering a flexible, transparent, and rigorous
approach. The proposed methodology addresses recent critiques of capture-
recapture models, offering a weaker and more interpretable identifying
assumption and accommodating complex heterogeneous capture probabilities
depending on high-dimensional or continuous covariates. We believe that our
framework has broad applicability in various domains where understanding the
size of a hidden population is critical, such as biology, public health, human
rights, and many others.

## References

  * Aleshin-Guendel et al., (2024) Aleshin-Guendel, S., Sadinle, M., and Wakefield, J. (2024).  The central role of the identifying assumption in population size estimation.  Biometrics, 80(1):ujad028.
  * Audibert and Tsybakov, (2007) Audibert, J.-Y. and Tsybakov, A. B. (2007).  Fast learning rates for plug-in classifiers.  The Annals of Statistics, 35(2):608 – 633.
  * Ball et al., (2003) Ball, P., Asher, J., Sulmont, D., and Manrique, D. (2003).  How many peruvians have died? an estimate of the total number of victims killed or disappeared in the armed internal conflict between 1980 and 2000.  Report, American Association for the Advancement of Science (AAAS).
  * Bickel et al., (1993) Bickel, P. J., Klaassen, C. A., Bickel, P. J., Ritov, Y., Klaassen, J., Wellner, J. A., and Ritov, Y. (1993).  Efficient and adaptive estimation for semiparametric models, volume 4.  Springer.
  * Binette and Steorts, (2022) Binette, O. and Steorts, R. C. (2022).  On the reliability of multiple systems estimation for the quantification of modern slavery.  Journal of the Royal Statistical Society Series A: Statistics in Society, 185(2):640–676.
  * Bishop et al., (2007) Bishop, Y. M., Fienberg, S. E., and Holland, P. W. (2007).  Discrete multivariate analysis: Theory and practice.  Springer Science & Business Media.
  * Böhning et al., (2020) Böhning, D., Rocchetti, I., Maruotti, A., and Holling, H. (2020).  Estimating the undetected infections in the covid-19 outbreak by harnessing capture–recapture methods.  International Journal of Infectious Diseases, 97:197–201.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine learning, 45:5–32.
  * Burgar et al., (2018) Burgar, J. M., Stewart, F. E., Volpe, J. P., Fisher, J. T., and Burton, A. C. (2018).  Estimating density for species conservation: Comparing camera trap spatial count models to genetic spatial capture-recapture models.  Global Ecology and Conservation, 15:e00411.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68.
  * Das et al., (2023) Das, M., Kennedy, E. H., and Jewell, N. P. (2023).  Doubly robust capture-recapture methods for estimating population size.  Journal of the American Statistical Association, pages 1–13.
  * Fienberg, (1972) Fienberg, S. E. (1972).  The multiple recapture census for closed populations and incomplete 2k contingency tables.  Biometrika, 59(3):591–603.
  * Hook and Regal, (1995) Hook, E. B. and Regal, R. R. (1995).  Capture-recapture methods in epidemiology: methods and limitations.  Epidemiologic reviews, 17(2):243–264.
  * Imbens and Manski, (2004) Imbens, G. W. and Manski, C. F. (2004).  Confidence intervals for partially identified parameters.  Econometrica, 72(6):1845–1857.
  * Johndrow et al., (2019) Johndrow, J., Lum, K., and Manrique-Vallier, D. (2019).  Low-risk population size estimates in the presence of capture heterogeneity.  Biometrika, 106(1):197–210.
  * Kennedy, (2022) Kennedy, E. H. (2022).  Semiparametric doubly robust targeted double machine learning: a review.  arXiv preprint arXiv:2203.06469.
  * Kennedy et al., (2020) Kennedy, E. H., Balakrishnan, S., and G’Sell, M. (2020).  Sharp instruments for classifying compliers and generalizing causal effects.
  * Kennedy et al., (2023) Kennedy, E. H., Balakrishnan, S., and Wasserman, L. (2023).  Semiparametric counterfactual density estimation.  Biometrika, 110(4):875–896.
  * Levis et al., (2023) Levis, A. W., Bonvini, M., Zeng, Z., Keele, L., and Kennedy, E. H. (2023).  Covariate-assisted bounds on causal effects with instrumental variables.  arXiv preprint arXiv:2301.12106.
  * Lincoln, (1930) Lincoln, F. C. (1930).  Calculating waterfowl abundance on the basis of banding returns.  Number 118. US Department of Agriculture.
  * Luedtke and Van Der Laan, (2016) Luedtke, A. R. and Van Der Laan, M. J. (2016).  Statistical inference for the mean outcome under a possibly non-unique optimal treatment strategy.  Annals of statistics, 44(2):713.
  * Manrique-Vallier, (2016) Manrique-Vallier, D. (2016).  Bayesian population size estimation using dirichlet process mixtures.  Biometrics, 72(4):1246–1254.
  * Manrique-Vallier et al., (2021) Manrique-Vallier, D., Ball, P., and Sadinle, M. (2021).  Capture-recapture for casualty estimation and beyond: recent advances and research directions.  Statistics in the Public Interest: In Memory of Stephen E. Fienberg, pages 15–31.
  * Manrique-Vallier et al., (2019) Manrique-Vallier, D., Ball, P., and Sulmont, D. (2019).  Estimating the number of fatal victims of the peruvian internal armed conflict, 1980-2000: an application of modern multi-list capture-recapture techniques.  arXiv preprint arXiv:1906.04763.
  * Manrique-Vallier et al., (2013) Manrique-Vallier, D., Price, M. E., and Gohdes, A. (2013).  Multiple systems estimation techniques for estimating casualties in armed conflicts.  Counting civilian casualties: An introduction to recording and estimating nonmilitary deaths in conflict, 165.
  * Otis et al., (1978) Otis, D. L., Burnham, K. P., White, G. C., and Anderson, D. R. (1978).  Statistical inference from capture data on closed animal populations.  Wildlife monographs, (62):3–135.
  * Petersen, (1896) Petersen, C. G. J. (1896).  The yearly immigration of young plaice in the limfjord from the german sea.  Rept. Danish Biol. Sta., 6:1–48.
  * Rendon, (2019) Rendon, S. (2019).  Capturing correctly: A reanalysis of the indirect capture–recapture methods in the peruvian truth and reconciliation commission.  Research & Politics, 6(1):2053168018820375.
  * Rivest and Baillargeon, (2014) Rivest, L.-P. and Baillargeon, S. (2014).  Capture-recapture methods for estimating the size of a population: dealing with variable capture probabilities.  Statistics in Action: A Canadian Outlook, 40:289–304.
  * Schwarz and Seber, (1999) Schwarz, C. J. and Seber, G. A. (1999).  Estimating animal abundance: review iii.  Statistical Science, 14(4):427–456.
  * Silverman, (2020) Silverman, B. W. (2020).  Multiple-systems analysis for the quantification of modern slavery: classical and bayesian approaches.  Journal of the Royal Statistical Society Series A: Statistics in Society, 183(3):691–736.
  * van der Vaart, (2002) van der Vaart, A. W. (2002).  Semiparametric statistics.  Lecture Notes in Math., (1781).
  * Wesson et al., (2024) Wesson, P., Das, M., Chen, M., Hsu, L., McFarland, W., Kennedy, E., and Jewell, N. P. (2024).  Evaluating a targeted minimum loss-based estimator for capture-recapture analysis: An application to hiv surveillance in san francisco, california.  American Journal of Epidemiology, 193(4):673–683.
  * You et al., (2021) You, Y., van der Laan, M., Collender, P., Cheng, Q., Hubbard, A., Jewell, N. P., Hu, Z. T., Mejia, R., and Remais, J. (2021).  Estimation of population size based on capture recapture designs and evaluation of the estimation reliability.  arXiv preprint arXiv:2105.05373.

## Appendix

### A. No highest-order interaction in conditional log-linear models

In the following, we build some intuition on the interpretation of the
highest-order interaction term in (conditional) log-linear models and the
corresponding identification assumption. We anchor its interpretation to
deviations from a conditional independence assumption and show how it can
account for more general capture-recapture settings. Without loss of
generality, we state the different identification assumptions in terms of the
first two lists $Y_{1},Y_{2}$.

First, note that in the $K=2$ list case the highest-order interaction term can
be expressed as

$\exp(\alpha_{1}(X))=\frac{p_{00}(X)p_{11}(X)}{p_{10}(X)p_{01}(X)},$ (37)

where $p_{y}(X)=\mathbb{P}(Y=y\mid X)$ is the true probability of having
capture profile $y=(y_{1},y_{2})$, given covariates $X$. Assuming no highest-
order interaction implies

$\displaystyle\exp(0)=1=\frac{p_{00}(X)p_{11}(X)}{p_{10}(X)p_{01}(X)}$
$\displaystyle\iff\frac{p_{11}(X)}{p_{01}(X)}=\frac{p_{10}(X)}{p_{00}(X)}$
$\displaystyle\iff odds(Y_{1}=1\mid Y_{2}=1,X)=odds(Y_{1}=1\mid Y_{2}=0,X).$

This condition states that the odds of appearing on list 1 conditional on
appearing on list 2 and covariates $X$ is equal to the odds of appearing on
list 1 conditional on not appearing on list 2 and covariates $X$. Thus, in the
$K=2$ list case, the no highest-order interaction assumption is equivalent to
conditional independence between lists: $Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid
X$.

Similarly, for the $K=3$ lists case we have that

$\displaystyle\exp(\alpha_{1}(X))$
$\displaystyle=\frac{p_{100}(X)p_{010}(X)p_{001}(X)p_{111}(X)}{p_{000}(X)p_{110}(X)p_{101}(X)p_{011}(X)}$
$\displaystyle=\frac{odds(Y_{1}\mid Y_{2}=0,Y_{3}=0,X)odds(Y_{1}\mid
Y_{2}=1,Y_{3}=1,X)}{odds(Y_{1}\mid Y_{2}=1,Y_{3}=0,X)odds(Y_{1}\mid
Y_{2}=0,Y_{3}=1,X)}$
$\displaystyle=\frac{OR_{Y_{1},Y_{2}}(Y_{3}=1,X)}{OR_{Y_{1},Y_{2}}(Y_{3}=0,X)}.$

To understand this expression, first note that under the conditional
independence assumption $Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid Y_{3},X$, we have
that $OR_{Y_{1},Y_{2}}(Y_{3}=1,X)=OR_{Y_{1},Y_{2}}(Y_{3}=0,X)=1$. Moving away
from this assumption, which can be restrictive in several capture-recapture
settings, the assumption of no highest-order interaction implies that these
odds ratios are equal

$OR_{Y_{1},Y_{2}}(Y_{3}=1,X)=OR_{Y_{1},Y_{2}}(Y_{3}=0,X),$

but not necessarily equal to 1.

In the $K=4$ case, the highest-order interaction term in a saturated log-
linear model takes the following form, where
$odds(Y_{1}\mid(y_{2}y_{3}y_{4}),X)=odds(Y_{1}\mid
Y_{2}=y_{2},Y_{3}=y_{3},Y_{4}=y_{4},X)$,

$\displaystyle\exp(\alpha_{1}(X))$
$\displaystyle=\frac{p_{0000}(X)p_{1100}(X)p_{1010}(X)p_{1001}(X)p_{0110}(X)p_{0101}(X)p_{0011}(X)p_{1111}(X)}{p_{1000}(X)p_{0100}(X)p_{0010}(X)p_{0001}(X)p_{1110}(X)p_{1101}(X)p_{1011}(X)p_{0111}(X)}$
$\displaystyle=\frac{odds(Y_{1}|(100),X)odds(Y_{1}|(010),X)odds(Y_{1}|(001),X)odds(Y_{1}|(111),X)}{odds(Y_{1}|(000),X)odds(Y_{1}|(110),X)odds(Y_{1}|(101),X)odds(Y_{1}|(011),X)}$
$\displaystyle=\frac{OR_{Y_{1},Y_{2}}(Y_{3}=0,Y_{4}=0,X)OR_{Y_{1},Y_{2}}(Y_{3}=1,Y_{4}=1,X)}{OR_{Y_{1},Y_{2}}(Y_{3}=1,Y_{4}=0,X)OR_{Y_{1},Y_{2}}(Y_{3}=0,Y_{4}=1,X)}.$

Importantly, a similar conditional independence assumption
$Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid Y_{3},Y_{4},X$ implies that
$OR_{Y_{1},Y_{2}}(Y_{3}=y_{3},Y_{4}=y_{4},X)=1$ for all combinations of
$(y_{3},y_{4})$. The no highest-order interaction assumption relaxes this
condition, but still imposes certain symmetry in the ORs such that

$\frac{OR_{Y_{1},Y_{2}}(Y_{3}=1,Y_{4}=1,X)}{OR_{Y_{1},Y_{2}}(Y_{3}=0,Y_{4}=1,X)}=\frac{OR_{Y_{1},Y_{2}}(Y_{3}=1,Y_{4}=0,X)}{OR_{Y_{1},Y_{2}}(Y_{3}=0,Y_{4}=0,X)}.$

Finally, for the general case of $K$ lists, we have that the highest-order
interaction satisfies the identity

$\exp((-1)^{K}\alpha_{1}(X))=\frac{\prod_{\\{y\leavevmode\nobreak\
even\\}}p_{y}(X)}{\prod_{\\{y\leavevmode\nobreak\
odd\\}}p_{y}(X)}=\frac{\prod_{\\{(y_{3},...,y_{K})\leavevmode\nobreak\
even\\}}OR_{Y_{1},Y_{2}}((y_{3},...,y_{K}),X)}{\prod_{\\{(y_{3},...,y_{K})\leavevmode\nobreak\
odd\\}}OR_{Y_{1},Y_{2}}((y_{3},...,y_{K}),X)},$ (38)

where “$y$ even/odd” denotes, with a slight abuse of notation, that the number
of non-zero entries of the capture profile vector $y$ ($L_{1}$ norm) is
even/odd. Analogous to the previous cases, imposing the conditional
independence assumption of $Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid
Y_{3},...,Y_{K},X$, implies that all the odds ratios appearing in (38) are
equal to $1$, which can be too restrictive in real-world applications for
complex capture-recapture settings with stronger dependencies across capture
episodes. On the other hand, the no highest-order interaction approach can be
seen as a relaxation of the conditional independence condition assuming
certain symmetry on the ORs such that

$\prod_{\\{(y_{3},...,y_{K})\leavevmode\nobreak\
even\\}}OR_{Y_{1},Y_{2}}((y_{3},...,y_{k}),X)=\prod_{\\{(y_{3},...,y_{k})\leavevmode\nobreak\
odd\\}}OR_{Y_{1},Y_{2}}((y_{3},...,y_{k}),X).$

This condition will be satisfied if all ORs are equal to one (i.e., if the
conditional independence assumption holds) or if deviations from conditional
independence “cancel out” across terms. This can be the case, for example, if
all but two ORs are equal to $1$ and
$OR_{Y_{1},Y_{2}}((1,y_{4},...,y_{k}),X)=OR_{Y_{1},Y_{2}}((0,y_{4},...,y_{k}),X)=2$.
The proposed sensitivity analysis using a bounded highest-order interaction
assumption further relaxes the conditional independence assumption and this
idealized symmetry in the odds ratios, and allows for these to vary up to the
bound $\exp(\delta)$. This accounts for more arbitrary dependencies between
capture episodes such that the deviation from this symmetry in ORs is bounded
by $\exp(\delta)$.

In general, for any number of $K$ lists in a conditional log-linear model,
assuming that two lists are conditionally independent given the other lists
and covariates is a stricter identification assumption than the no highest-
order interaction. The latter identifying assumption accounts for more complex
capture-recapture dependencies across lists, but is harder to interpret for a
number of lists larger than $K>3$. The proposed identification strategy
introduced in our work generalizes these two approaches noting that it is
sufficient for identification that a subset of the $K$ lists has no highest-
order interaction, conditional on not being observed by the remaining lists
and covariates. For instance,

$Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid X,Y_{3}=0,\dots,Y_{K}=0,$

is sufficient to identify the capture probability $\psi$.

######  Remark 4.

The conditional independence assumption $Y_{1}\perp\\!\\!\\!\perp Y_{2}\mid
X,Y_{3},\dots,Y_{K}$ further imposes some restrictions on the coefficients of
the saturated log linear model implied by the condition that
$OR_{Y_{1},Y_{2}}((y_{3},...,y_{k}),X)=1$ for all $(y_{3},...,y_{k})$. These
extra modeling conditions arising from the conditional independence
identifying assumption (which are testable from the available data) constrain
the models that can be used for estimating $\mathbb{Q}$. Therefore, it
restricts analysts to semiparametric modeling of the capture probability to
distributions that satisfy these additional conditions. On the other hand, the
no highest-order interaction assumptions and conditioning only on not being
observed by the remaining lists do not impose any restrictions on the observed
distribution $\mathbb{Q}$, allowing for flexible nonparametric data-driven
modeling of $\psi$.

### B. Proofs

###### Proof of Proposition 1.

The highest-order interaction term in a saturated conditional log-linear model
on $J$ lists, conditional on being unobserved by the other $K-J$ lists and
covariates, can be expressed as

$\alpha_{1}(X)=\sum_{(y_{1},\dots,y_{J})}(-1)^{J+|y|}\log\left[\mathbb{P}\\{(Y_{1},\dots,Y_{J})=(y_{1},\dots,y_{J})\mid
X,(Y_{J+1},\dots,Y_{K})=0\\}\right].$

Note that we can write the log conditional probability for each capture
profile $(y_{1},\dots,y_{J})$ as

$\log\left[\frac{\mathbb{P}\\{(Y_{1},\dots,Y_{J},Y_{J+1},\dots,Y_{K})=(y_{1},\dots,y_{J},0,\dots,0)\mid
X\\}}{\mathbb{P}\\{(Y_{J+1},\dots,Y_{K})=0\mid X\\}}\right],$

and the terms $\log[\mathbb{P}\\{(Y_{J+1},\dots,Y_{K})=0\mid X\\}]$ will
cancel out from the same number of even and odd capture profiles in the sum in
$\alpha_{1}(X)$. Then we can use the relationship between the population and
the observed distributions $\mathbb{P}(Y=y\mid X)=q_{y}(X)\gamma(X)$, for all
$y\neq 0$, to rewrite the highest-order interaction coefficient in terms of
the $q$-probabilities and the conditional capture probability $\gamma(X)$:

$\displaystyle\alpha_{1}(X)$
$\displaystyle=\sum_{(y_{1},\dots,y_{J})}(-1)^{J+|y|}\log[\mathbb{P}(Y=(y_{1},\dots,y_{J},0,\dots,0)\mid
X)]$ $\displaystyle=\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{J+|y|}\log[q_{(y_{1},\dots,y_{J},0,\dots,0)}(X)\gamma(X)]+(-1)^{J}\log(1-\gamma(X)).$

Assuming no highest-order interaction $\alpha_{1}(X)=0$, we solve for
$\gamma(X)$ and use that $\psi^{-1}$ can be expressed as the empirical
harmonic mean $\mathbb{E}_{\mathbb{Q}}[\gamma(X)^{-1}]$ to obtain the
identifying expression.

∎

###### Proof of Lemma 1.

We show that the proposed expression

$\phi(X,Y)=\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}+1-\frac{1}{\psi}$

is the Efficient Influence Function for $\psi^{-1}$ by computing the remainder
term of its von Mises expansion (Bickel et al.,, 1993). Let $\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu=\psi(\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern 1.5mu)$ be the plug-
in estimator of the capture probability for a generic distribution $\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern 1.5mu$. Then, the
remainder of the von Mises expansion is given by

$\displaystyle R_{2}(\mathbb{Q},\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern
1.5mu)=\frac{1}{\mkern 1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu}-\frac{1}{\psi}+\mathbb{E}_{\mathbb{Q}}\left(\phi(X,Y;\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern 1.5mu)\right)$
$\displaystyle=\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu}-\frac{1}{\psi}+\mathbb{E}_{\mathbb{Q}}\left[\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}(X)}+1-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}}\left[\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{q_{y}(X)}{\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}(X)}+1-\frac{1}{\gamma(X)}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}}\left[\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\left(\frac{q_{y}(X)}{\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}(X)}-1\right)+\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}+1-\frac{1}{\gamma(X)}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}}\left[\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\left(\frac{q_{y}(X)}{\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}(X)}-1\right)+\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu(X)}-1\right)-\left(\frac{1}{\gamma(X)}-1\right)\right].$

The result follows from computing the second-order Taylor expansion for

$\frac{1}{\gamma(X)}-1=\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{1+|y|}\log(q_{y}(X))\right)$

around $\\{\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}(X)\\}_{y\neq 0}$, such that

$\displaystyle\frac{1}{\gamma}-1$ $\displaystyle=\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu}-1\right)+\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\left(\frac{q_{y}}{\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern 1.5mu_{y}}-1\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\left(\frac{1}{\widetilde{\gamma}}-1\right)\left[\sum_{y_{i}\neq y_{j}\neq
0}(-1)^{|y_{i}|+|y_{j}|}\frac{(q_{y_{i}}-\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y_{i}})}{\widetilde{q}_{y_{i}}}\frac{(q_{y_{j}}-\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y_{j}})}{\widetilde{q}_{y_{j}}}+\sum_{|y|\neq 0,\leavevmode\nobreak\
even}\left(\frac{q_{y}-\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}}{\widetilde{q}_{y}}\right)^{2}\right]$

for some value $\widetilde{q}$ of the $q$-probabilities between $q_{y}$ and
$\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern 1.5mu_{y}$, where we
drop the $X$ argument from $q_{y}(X)$ and $\gamma(X)$ for ease of exposition.
In the last expression, the sums are taken with respect to all possible
capture profiles of the form $y=(y_{1},\dots,y_{J},0,\dots,0)$.

Note that the first two terms in $R_{2}(\mathbb{Q},\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern 1.5mu)$ cancel out
with the linear terms in the Taylor expansion for $\gamma^{-1}-1$. Then, the
remainder is given by

$R_{2}(\mathbb{Q},\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern
1.5mu)=\mathbb{E}_{\mathbb{Q}}\left[\left(\frac{1}{\widetilde{\gamma}}-1\right)\left\\{\sum_{y_{i}\neq
y_{j}\neq 0}(-1)^{|y_{i}|+|y_{j}|}\frac{(q_{y_{i}}-\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y_{i}})}{\widetilde{q}_{y_{i}}}\frac{(q_{y_{j}}-\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y_{j}})}{\widetilde{q}_{y_{j}}}+\sum_{|y|\neq 0,\leavevmode\nobreak\
even}\left(\frac{q_{y}-\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}}{\widetilde{q}_{y}}\right)^{2}\right\\}\right]$

The last equality shows that the remainder term $R_{2}(\mathbb{Q},\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern 1.5mu)$ only depends
on products of errors of the nuisance functions, which implies in particular
that
$\frac{d}{d\epsilon}R_{2}(\mathbb{Q},\mathbb{Q}_{\epsilon})\big{|}_{\epsilon=0}=0$
for any parametric submodel $\mathbb{Q}_{\epsilon}$. Therefore, by Lemma 2 in
Kennedy et al., (2023), this implies that the proposed $\phi(X,Y)$ is the
Influence Function for $\psi^{-1}$. Moreover, the influence function is unique
and efficient provided that our nonparametric model does not restrict the
tangent space (van der Vaart,, 2002; Bickel et al.,, 1993), concluding our
proof.

∎

###### Proof of Theorem 1.

The efficiency bound for $\psi^{-1}=\mathbb{E}[\gamma(X)^{-1}]$ can be
computed as the variance of its efficient influence function $\phi(X,Y)$.
Using the law of total variance, we have that

$\text{Var}(\phi(X,Y))=\text{Var}\left(\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}+1-\frac{1}{\psi}\right)=T_{1}+T_{2},$

where

$\displaystyle T_{1}$
$\displaystyle=\text{Var}\left(\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}\leavevmode\nobreak\
\Bigg{|}\leavevmode\nobreak\ X\right]\right),$ $\displaystyle T_{2}$
$\displaystyle=\mathbb{E}\left[\text{Var}\left(\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}\leavevmode\nobreak\
\Bigg{|}\leavevmode\nobreak\ X\right)\right].$

The first term can be simplified to

$\displaystyle T_{1}$
$\displaystyle=\text{Var}\left(\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbb{E}[\mathbbm{1}(Y=y)\mid X]}{q_{y}(X)}\right)$
$\displaystyle=\text{Var}\left(\left(\frac{1}{\gamma(X)}-1\right)\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\right)=\text{Var}\left(\frac{1}{\gamma(X)}\right).$

We now analyze the second term $T_{2}$:

$\displaystyle
T_{2}=\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\text{Var}\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\frac{\mathbbm{1}(Y=y)}{q_{y}(X)}\leavevmode\nobreak\
\Bigg{|}\leavevmode\nobreak\ X\right)\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\left\\{\sum_{(y_{1},\dots,y_{J})\neq
0}\frac{\text{Var}(\mathbbm{1}(Y=y)|X)}{q_{y}(X)^{2}}+2\sum_{y_{i}\neq
y_{j}\neq
0}\frac{(-1)^{|y_{i}|+|y_{j}|}}{q_{y_{i}}(X)q_{y_{j}}(X)}Cov(\mathbbm{1}(Y=y_{i}),\mathbbm{1}(Y=y_{j})|X)\right\\}\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\left\\{\sum_{(y_{1},\dots,y_{J})\neq
0}\frac{q_{y}(X)(1-q_{y}(X))}{q_{y}(X)^{2}}+2\sum_{y_{i}\neq y_{j}\neq
0}\frac{(-1)^{|y_{i}|+|y_{j}|}}{q_{y_{i}}(X)q_{y_{j}}(X)}\left(\mathbb{E}[\mathbbm{1}_{y_{i}}\mathbbm{1}_{y_{j}}|X]-q_{y_{i}}(X)q_{y_{j}}(X)\right)\right\\}\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\left\\{\sum_{(y_{1},\dots,y_{J})\neq
0}\frac{1-q_{y}(X)}{q_{y}(X)}-2\sum_{y_{i}\neq y_{j}\neq
0}\frac{(-1)^{|y_{i}|+|y_{j}|}}{q_{y_{i}}(X)q_{y_{j}}(X)}q_{y_{i}}(X)q_{y_{j}}(X)\right\\}\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\left\\{\sum_{(y_{1},\dots,y_{J})\neq
0}\frac{1}{q_{y}(X)}-\sum_{(y_{1},\dots,y_{J})\neq 0}1-2\sum_{y_{i}\neq
y_{j}\neq 0}(-1)^{|y_{i}|+|y_{j}|}\right\\}\right]$
$\displaystyle=\mathbb{E}\left[\left(\frac{1}{\gamma(X)}-1\right)^{2}\left\\{\sum_{(y_{1},\dots,y_{J})\neq
0}\frac{1}{q_{y}(X)}-1\right\\}\right],$

where the last equality follows from

$\displaystyle\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{2|y|}+2\left(\sum_{y_{i}\neq y_{j}\neq
0}(-1)^{|y_{i}|+|y_{j}|}\right)=\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|}\right)^{2}=1.$

Putting the two terms together yields the result. ∎

###### Proof of Theorem 2.

First note that by the definition of our one-step estimator and the von Mises
expansion of $\psi$, we have that

$\displaystyle\frac{1}{\widehat{\psi}_{os}}-\frac{1}{\psi}-\mathbb{Q}_{N}(\phi)$
$\displaystyle=\frac{1}{\widehat{\psi}_{pi}}+\mathbb{Q}_{N}(\widehat{\phi})-\frac{1}{\psi}-\mathbb{Q}_{N}(\phi)$
$\displaystyle=-\mathbb{Q}(\widehat{\phi})+\widehat{R}_{2}+\mathbb{Q}_{N}(\widehat{\phi})-\mathbb{Q}_{N}(\phi)$
$\displaystyle=\widehat{R}_{2}+(\mathbb{Q}_{N}-\mathbb{Q})(\widehat{\phi}-\phi),$

where $\widehat{R}_{2}$ is the remainder term of the von Mises expansion (15)
derived in the proof of Lemma 1. The last equality uses the property that the
EIF has mean zero, $\mathbb{Q}(\phi)=0$

Following the proof of Theorem 2 in Das et al., (2023) we have that

$\mathbb{E}\left[\Bigg{|}\frac{1}{\widehat{\psi}_{os}}-\frac{1}{\psi}-\mathbb{Q}_{N}(\phi)\Bigg{|}^{2}\right]=\mathbb{E}\left[\Big{|}\widehat{R}_{2}+(\mathbb{Q}_{N}-\mathbb{Q})(\widehat{\phi}-\phi))\Big{|}^{2}\right]\leq\mathbb{E}[\widehat{R}_{2}^{2}]+\frac{1}{N}\mathbb{E}\left[\left\lVert\widehat{\phi}-\phi\right\rVert^{2}\right].$

The result follows by applying Markov’s inequality for error tolerance
$\eta>0$

$\displaystyle\mathbb{P}\left[\Bigg{|}\frac{1}{\widehat{\psi}_{os}}-\frac{1}{\psi}-\mathbb{Q}_{N}(\phi)\Bigg{|}\leq\eta\right]$
$\displaystyle\geq
1-\frac{1}{\eta^{2}}\mathbb{E}\left[\Bigg{|}\frac{1}{\widehat{\psi}_{os}}-\frac{1}{\psi}-\mathbb{Q}_{N}(\phi)\Bigg{|}^{2}\right]$
$\displaystyle\geq
1-\frac{1}{\eta^{2}}\mathbb{E}\left[\widehat{R}_{2}^{2}+\frac{1}{N}\left\lVert\widehat{\phi}-\phi\right\rVert^{2}\right].$

∎

###### Proof of Theorem 3.

The proof follows the same logic as the proof of Theorem 3 in Das et al.,
(2023). In detail, note that by the definition of the proposed one-step
estimator we have that
$\widehat{\psi}^{-1}_{os}-\widehat{\psi}-\widehat{R}_{2}=(\mathbb{Q}_{N}-\mathbb{Q})\widehat{\phi}$
is a sample average of a fixed function conditional on the training sample
$Z^{N}=\\{(X_{i},Y_{i}\\}_{i=1}^{N}$. Therefore, by Berry-Essen Theorem we
have that

$\Phi(t^{\prime})-\frac{C\rho}{\widetilde{\sigma}^{3}\sqrt{N}}\leq\mathbb{P}\left(\frac{\widehat{\psi}^{-1}_{os}-\widehat{\psi}-\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}}\leq
t^{\prime}\leavevmode\nobreak\ \Bigg{|}\leavevmode\nobreak\
Z^{N}\right)\leq\Phi(t^{\prime})+\frac{C\rho}{\widetilde{\sigma}^{3}\sqrt{N}},$

for any $t^{\prime}$ and $N$, where
$\widetilde{\sigma}^{2}=\text{Var}(\widehat{\phi}\mid Z^{N})$ and
$\rho=\mathbb{E}[|\widehat{\phi}-\mathbb{Q}\widehat{\phi}|^{3}\mid Z^{N}]$.

Taking
$t^{\prime}=\frac{\widehat{\sigma}}{\widetilde{\sigma}}t-\frac{\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}}$,
where $\widehat{\sigma}^{2}=\widehat{\text{Var}}(\widehat{\phi})$ is the
estimated unconditional variance of the estimated EIF, we can rewrite the
previous bounds in terms of the estimation error
$\widehat{\psi}^{-1}_{os}-\widehat{\psi}$:

$\Phi\left(\frac{\widehat{\sigma}}{\widetilde{\sigma}}t-\frac{\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}}\right)-\frac{C\rho}{\widetilde{\sigma}^{3}\sqrt{N}}\leq\mathbb{P}\left(\frac{\widehat{\psi}^{-1}_{os}-\widehat{\psi}}{\widehat{\sigma}/\sqrt{N}}\leq
t\leavevmode\nobreak\ \Bigg{|}\leavevmode\nobreak\
Z^{N}\right)\leq\Phi\left(\frac{\widehat{\sigma}}{\widetilde{\sigma}}t-\frac{\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}}\right)+\frac{C\rho}{\widetilde{\sigma}^{3}\sqrt{N}}.$

Finally, note that using the mean value theorem we can bound the difference in
cumulative distributions by

$\displaystyle\Bigg{|}\Phi\left(\frac{\widehat{\sigma}}{\widetilde{\sigma}}t-\frac{\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}}\right)-\Phi(t)\Bigg{|}=\Bigg{|}\Phi^{\prime}(t_{N})\left[\frac{\widehat{\sigma}}{\widetilde{\sigma}}t-\frac{\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}}-t\right]\Bigg{|}\leq\frac{1}{\sqrt{2\pi}}\left(\Big{|}\frac{\widehat{\sigma}}{\widetilde{\sigma}}-1\Big{|}|t|+\frac{|\widehat{R}_{2}|}{\widetilde{\sigma}/\sqrt{N}}\right),$

for some
$t_{N}\in(t,\frac{\widehat{\sigma}}{\widetilde{\sigma}}t-\frac{\widehat{R}_{2}}{\widetilde{\sigma}/\sqrt{N}})$,
where the last inequality uses that $\text{sup}_{t}\Phi^{\prime}(t)\leq
1/\sqrt{2\pi}$ and the triangle inequality. The result follows by iterated
expectations:

$\Bigg{|}\mathbb{P}\left(\frac{\widehat{\psi}^{-1}_{os}-\psi^{-1}}{\widehat{\sigma}/\sqrt{N}}\leq
t\right)-\Phi(t)\Bigg{|}\leq\frac{1}{\sqrt{2\pi}}\left(|t|\mathbb{E}_{\mathbb{Q}}\left[\Big{|}\frac{\widehat{\sigma}}{\widetilde{\sigma}}-1\Big{|}\right]+\sqrt{N}\mathbb{E}_{\mathbb{Q}}\left[\frac{|\widehat{R}_{2}|}{\widetilde{\sigma}}\right]\right)+\frac{C}{\sqrt{N}}\mathbb{E}_{\mathbb{Q}}\left[\frac{\rho}{\widetilde{\sigma}^{3}}\right].$

∎

###### Proof of Proposition 2.

Following the proof of Proposition 1, we can express the highest-order
coefficient of the proposed log-linear model as

$\alpha_{1}(X)=\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{J+|y|}\log(q_{y}(X)\gamma(X))+(-1)^{J}\log(1-\gamma(X)),$

for $q$-probabilities $q_{y}(X)=\mathbb{Q}(Y=(y_{1},\dots,y_{J},0,\dots,0)\mid
X)$. Then, a bounded highest-order coefficient $|\alpha_{1}(X)|\leq\delta$
implies that the conditional capture probability is partially identified by

$\frac{1}{\gamma(X)}\in\left[1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\pm\delta\right)\right].$

Next, note that the parameter $\delta$ must be set such that
$\gamma(X)\in[\epsilon,1]$ to ensure the validity of the conditional capture
probabilities it yields. Consequently, we derive the following upper and lower
bounds on $\gamma(X)^{-1}$ to partially identify our target parameter while
incorporating such constraint on $\delta$:

$\frac{1}{\gamma(X)}\in\left[\frac{1}{\gamma_{\ell}(X)},\frac{1}{\gamma_{u}(X)}\right]=\left[\min\left\\{1+\exp\left(\sum_{(y_{1},\dots,y_{J})\neq
0}(-1)^{|y|+1}\log(q_{y}(X))\pm\delta\right),\frac{1}{\epsilon}\right\\}\right].$
(39)

Finally, we establish bounds on the inverse capture probability
$\psi^{-1}\in[\psi_{\ell}^{-1},\psi_{u}^{-1}]$, defined by
$\psi_{\ell}^{-1}=\mathbb{E}_{\mathbb{Q}}\left[\gamma_{\ell}^{-1}(X)\right]$
and $\psi_{u}^{-1}=\mathbb{E}_{\mathbb{Q}}\left[\gamma_{u}^{-1}(X)\right]$,
derived from the lower and upper bounds in (39), respectively.

∎

###### Proof of Theorem 4.

We follow the same approach as Lemma 1 to derive the remainder term in the von
Mises expansion for $\psi^{-1}_{\delta}$. We omit the arguments of the
functions for clarity. Then, We have that for a generic distribution $\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern 1.5mu$,

$\displaystyle R_{2,\delta}(\mathbb{Q},\mkern
1.5mu\overline{\mkern-1.5mu\mathbb{Q}\mkern-1.5mu}\mkern
1.5mu)=\frac{1}{\mkern 1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\psi_{\delta}}+\mathbb{E}_{\mathbb{Q}}\left(\mkern
1.5mu\overline{\mkern-1.5mu\phi\mkern-1.5mu}\mkern 1.5mu_{\delta}\right)$
$\displaystyle=\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\psi_{\delta}}+\mathbb{E}_{\mathbb{Q}}\left[\left(\mkern
1.5mu\overline{\mkern-1.5mu\varphi\mkern-1.5mu}\mkern
1.5mu_{\delta}-\frac{1}{\epsilon}\right)\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq
0\right)+\frac{1}{\epsilon}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu_{\delta}}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}}\left[\left(\mkern
1.5mu\overline{\mkern-1.5mu\varphi\mkern-1.5mu}\mkern
1.5mu_{\delta}-\frac{1}{\epsilon}\right)\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq
0\right)+\frac{1}{\epsilon}-\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\right)\mathbbm{1}\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\leq
0\right)-\frac{1}{\epsilon}\right]$
$\displaystyle=\mathbb{E}_{\mathbb{Q}}\left[\left(\mkern
1.5mu\overline{\mkern-1.5mu\varphi\mkern-1.5mu}\mkern
1.5mu_{\delta}-\frac{1}{\gamma_{\delta}}\right)\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq
0\right)+\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\right)\left[\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq
0\right)-\mathbbm{1}\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\leq
0\right)\right]\right].$

The first term can be bounded above using an analogous remainder term from
Lemma 1 given the symmetry of the EIFs for $\gamma$ and $\gamma_{\delta}$:

$\displaystyle\mathbb{E}_{\mathbb{Q}}\left[\left(\mkern
1.5mu\overline{\mkern-1.5mu\varphi\mkern-1.5mu}\mkern
1.5mu_{\delta}-\frac{1}{\gamma_{\delta}}\right)\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq 0\right)\right]$
$\displaystyle\leq\mathbb{E}_{\mathbb{Q}}\left[\Bigg{|}\mkern
1.5mu\overline{\mkern-1.5mu\varphi\mkern-1.5mu}\mkern
1.5mu_{\delta}-\frac{1}{\gamma_{\delta}}\Bigg{|}\right]$
$\displaystyle\leq\mathbb{E}_{\mathbb{Q}}\left[\sum_{y_{i}\neq y_{j}\neq
0}\left\lVert q_{y_{i}}-\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y_{i}}\right\rVert\left\lVert q_{y_{j}}-\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y_{j}}\right\rVert+\sum_{|y|\neq 0,\leavevmode\nobreak\
even}\left\lVert q_{y}-\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{y}\right\rVert^{2}\right],$

where the sums are taken over capture profiles of the form
$y=(y_{1},\dots,y_{J},0,\dots,0)$.

For the second term, we have that

$\displaystyle\mathbb{E}_{\mathbb{Q}}$
$\displaystyle\left[\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\right)\left[\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq
0\right)-\mathbbm{1}\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\leq
0\right)\right]\right]$
$\displaystyle\leq\mathbb{E}_{\mathbb{Q}}\left[\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\Bigg{|}\Bigg{|}\mathbbm{1}\left(\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}-\frac{1}{\epsilon}\leq
0\right)-\mathbbm{1}\left(\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\leq
0\right)\Bigg{|}\right]$
$\displaystyle\leq\mathbb{E}_{\mathbb{Q}}\left[\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\Bigg{|}\mathbbm{1}\left(\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\Bigg{|}\leq\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}\Bigg{|}\right)\right]$
$\displaystyle\leq\mathbb{E}_{\mathbb{Q}}\left[\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}\Bigg{|}\mathbbm{1}\left(\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\Bigg{|}\leq\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}\Bigg{|}\right)\right]$
$\displaystyle\leq\left\lVert\frac{1}{\gamma_{\delta}}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}\right\rVert_{\infty}\mathbb{P}\left(\Bigg{|}\frac{1}{\gamma_{\delta}}-\frac{1}{\epsilon}\Bigg{|}\leq\left\lVert\frac{1}{\gamma_{\delta}}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}\right\rVert_{\infty}\right)$ $\displaystyle\leq
C\left\lVert\frac{1}{\gamma_{\delta}}-\frac{1}{\mkern
1.5mu\overline{\mkern-1.5mu\gamma\mkern-1.5mu}\mkern
1.5mu_{\delta}}\right\rVert_{\infty}^{1+\beta}$

where the second inequality follows from Lemma 1 in Kennedy et al., (2020) and
the last inequality by the margin condition (26).

Finally, following similarly to Theorem 2, we have by Markov’s inequality that
for any error tolerance $\eta>0$

$\displaystyle\mathbb{P}\left[\Bigg{|}\frac{1}{\widehat{\psi}_{\delta}}-\frac{1}{\psi_{\delta}}-\mathbb{Q}_{N}(\phi_{\delta})\Bigg{|}\leq\eta\right]\geq
1-\frac{1}{\eta^{2}}\mathbb{E}\left[\widehat{R}_{2,\delta}^{2}+\frac{1}{N}\left\lVert\widehat{\phi}_{\delta}-\phi_{\delta}\right\rVert^{2}\right],$

concluding our proof.

∎

###### Proof of Theorem 5.

The proof follows the same logic that the proof of Theorem 3 noting that the
proposed one step estimator takes the form
$\widehat{\psi}^{-1}_{\delta}-\widehat{\psi}-\widehat{R}_{2,\delta}=(\mathbb{Q}_{N}-\mathbb{Q})\widehat{\phi}_{\delta}$,
which is a sample average of a fixed function conditional on the training
sample, where $\widehat{R}_{2,\delta}$ is the remainder term derived in
Theorem 4.

∎

###### Proof of Theorem 6.

By construction of our one-step estimator we have that

$\displaystyle\frac{1}{\widehat{\psi}_{os}}-\frac{1}{\psi}$
$\displaystyle=\frac{1}{\widehat{\psi}_{pi}}+\mathbb{Q}_{N}(\widehat{\phi})-\frac{1}{\psi}=-\mathbb{Q}(\widehat{\phi})+\widehat{R}_{2}+\mathbb{Q}_{N}(\widehat{\phi}),$

where $\widehat{R}_{2}$ is the remainder term of the von Mises expansion for
$\psi^{-1}$ (15) and
$\mathbb{Q}(\widehat{\phi})=\int\widehat{\phi}(z)d\mathbb{Q}(z)$. Therefore,
applying Theorem 4 in Das et al., (2023) we have that

$\widehat{CI}(n)=\left[\widehat{n}_{os}\pm
z_{\alpha/2}\sqrt{\widehat{n}_{os}\left(\widehat{\psi}_{os}\widehat{\sigma}^{2}+\frac{1-\widehat{\psi}_{os}}{\widehat{\psi}_{os}}\right)}\right],$

is an approximately $(1-\alpha)$ confidence interval for the population size
$n$ centered at $\widehat{n}_{os}=N/\widehat{\psi}_{os}$, where
$\widehat{\sigma}^{2}=\widehat{\text{Var}}(\widehat{\phi})$ is the unbiased
empirical variance of the estimated influence function.

∎

###### Proof of Theorem 7.

The results follow from Theorem 6 noting that our estimator for the generic
bound is of the required form

$\frac{1}{\widehat{\psi}_{\delta}}-\frac{1}{\psi}=-\mathbb{Q}(\widehat{\phi}_{\delta})+\widehat{R}_{2,\delta}+\mathbb{Q}_{N}(\widehat{\phi}_{\delta}),$

derived in Theorem 4.

∎