{'fraction_non_alphanumeric': 0.04171976551548363, 'fraction_numerical': 0.020405887600356824, 'mean_word_length': 4.525171624713959, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Doubly robust capture-recapture methods for estimating population size

Estimation of population size using incomplete lists (also called the capture-recapture problem) has a long history across many biological and social sciences. For example, human rights and other groups often construct partial and overlapping lists of victims of armed conflicts, with the hope of using this information to estimate the total number of victims. Earlier statistical methods for this setup either use potentially restrictive parametric assumptions, or else rely on typically suboptimal plug-in-type nonparametric estimators; however, both approaches can lead to substantial bias, the former via model misspecification and the latter via smoothing. Under an identifying assumption that two lists are conditionally independent given measured covariate information, we make several contributions. First, we derive the nonparametric efficiency bound for estimating the capture probability, which indicates the best possible performance of any estimator, and sheds light on the statistical limits of capture-recapture methods. Then we present a new estimator, and study its finite-sample properties, showing that it has a double robustness property new to capture-recapture, and that it is near-optimal in a non-asymptotic sense, under relatively mild nonparametric conditions. Next, we give a method for constructing confidence intervals for total population size from generic capture probability estimators, and prove non-asymptotic near-validity. Finally, we study our methods in simulations, and apply them to estimate the number of killings and disappearances attributable to different groups in Peru during its internal armed conflict between 1980 and 2000.


Introduction
Capture-recapture is a study design for estimating population size when only a fraction of the population is observed. This setup arises frequently, for example in studying ecological abundance, disease prevalence, and casualties in armed conflicts. Capture-recapture has a long history, dating back to at least Graunt in the 1600s (Hald, 2003), who used it to estimate plague prevalence in England. Similarly, in 1802 Laplace estimated the total population of France (Goudie and Goudie, 2007), and Petersen (1896) the abundance of plaice fish. More recently, it has been used in diverse settings ranging from estimating the number of pages on the web (Fienberg et al., 1999) to the total number of victims in a war (Ball et al., 2003), among many others.
The simplest capture-recapture setup, credited to Petersen (1896), consists of two independent lists with partial captures from the population of interest. There have been many generalizations over time. For our purposes, much of the previous work in capture-recapture can be viewed as falling within one of three streams. The first and oldest stream includes relatively simple data structures, e.g., involving no covariate information and relatively few lists (Darroch, 1958;Petersen, 1896;Schnabel, 1938). More recent advances in this stream include Burnham and Overton (1979) and Lee and Chao (1994). A second stream emerged to handle more intricate data structures, e.g., complex covariate information to help account for heterogeneity/dependence, largely using model-based approaches (Alho, 1990;Carothers, 1973;Fienberg, 1972;Huggins, 1989;Link, 2003;Pollock, 2002;Tilling and Sterne, 1999;Yip et al., 2001). However, the advantages of this second stream typically come at the expense of potentially restrictive parametric modeling assumptions, which when violated would induce bias. A third more recent stream addresses similar data structures as the second, but using more flexible nonparametric tools, e.g., local kernel or nonparametric Bayes or spline methods (Chen and Lloyd, 2000;Huggins andHwang, 2007, 2011;Kurtz, 2018;Manrique-Vallier, 2016;Stoklosa and Huggins, 2012;Yee et al., 2015;Zwane and van der Heijden, 2005). However, the work in this third stream has so far relied on interpretable but typically suboptimal plug-in estimators, which can suffer from nonparametric smoothing bias and slow rates of convergence (Robins et al., 2008;van der Laan and Robins, 2003;van der Vaart, 2014). We refer to Kurtz (2018) for a more detailed review of this third stream.
Our work takes the nonparametric perspective, but uses advances in efficiency theory to characterize optimality and improve simple plug-in estimators (Bickel et al., 1993;Kennedy, 2016;van der Laan and Robins, 2003). Under an identifying assumption that two lists are conditionally independent given measured covariate information (described in Section 2), we make several contributions.
• In Section 3 we derive the nonparametric efficiency bound for estimating the capture probability, which indicates the best possible performance of any estimator, and sheds light on the statistical limits of capture-recapture methods.
• In Section 4 we present a new doubly robust estimator, and study its finite-sample error, showing it is near-optimal in a non-asymptotic sense, under mild nonparametric conditions.
• In Section 5 we give a general method for constructing confidence intervals for total population size from generic capture probability estimators, and prove non-asymptotic near-validity.
• In Section 6 we study our methods in simulations, and apply them to estimate the number of killings and disappearances attributable to different groups in Peru during its internal armed conflict between 1980 and 2000.

Setup
Consider a finite population of n individuals, where the size n is unknown and to be estimated. We suppose there are K different lists of individuals from this population, yielding indicators Y ik ∈ {0, 1} of whether individual i ∈ {1, ..., n} appeared on list k ∈ {1, ..., K}. We let Y i = (Y i1 , ..., Y iK ) T denote the vector indicating list membership (i.e., capture profile) information for individual i. For example, in the K = 2 case, a profile Y i = (1, 0) T would mean that individual i appears on list 1 but not list 2. We consider the case where covariates X i ∈ R d are also available for each individual i = 1, ..., n. We assume an individual's chances of appearing on any given lists (and their covariates) do not depend on what happens with any other individuals, and also that the covariate and (conditional) list membership distributions are the same across individuals i = 1, ..., n. This implies that the random vectors Z i = (X i , Y i ) are independent and identically distributed according to some distribution P.
Remark 1. The setup above is commonly referred to as "heterogeneous" (Huggins, 1989;Pollock, 2002;Tilling and Sterne, 1999) since list membership Y can vary with covariates X. In other words, individuals with different covariates can have different chances of list membership.
Remark 2. In what follows we use the following standard notation. We let E Q denote an expectation under distribution Q, and let f 2 Q = f (z) 2 dQ(x) denote the corresponding squared L 2 (Q) norm; we let Q N denote the empirical measure under distribution Q. Finally we let a b mean a ≤ Cb for some universal constant C.
If every individual in the population appeared on at least one list (and could be uniquely identified), then the population size would of course be known without error; however in practice a possibly substantial fraction of individuals do not appear on any list. In other words, there are some individuals with Y = 0 that we do not observe. This means that, although the distribution P governs the capture profiles, we cannot sample from P directly. Instead we only see the N = n i=1 1(Y i = 0) individuals for whom Y ik = 1 for some k. This is illustrated in Figure 1.
List 3 (Y 3 = 1) Figure 1: Schematic of data structure for K = 3 lists. Observed data (i.e., those with Y = 0 in the union of the three lists) are represented with dark gray, while unobserved individuals with Y = 0 are represented with light gray. Individuals appearing in all three lists have Y = (1, 1, 1).
Hence the capture-recapture design is an example of biased sampling (Breslow et al., 2000;Qin, 2017;Vardi, 1985). In particular, the observed data Z i = (X i , Y i ), i = 1, ..., N are actually iid draws from a conditional distribution Q defined as is the (marginal) capture probability.
Remark 3. Some authors use N to denote the total population size and n for the observed number of captures; in contrast, we use N for the observed number and n for the total size, following the convention of saving upper case for random variables. Note the observed number of captures N is random since it depends on the random selection indicators, while the total population size n is fixed; specifically N ∼ Bin(n, ψ). Nonetheless much of our analysis will be conditional on the observed sample size N .
Recall our overall goal is to estimate the total population size n = N + i 1(Y i = 0). Since N ∼ Bin(n, ψ), the population size can be viewed as a fixed population parameter given by Intuitively, the lower the capture probability, the more the observed number N must be inflated to reflect the total population size. The quantity E(N ) in (3) can of course be unbiasedly estimated with the observed number of captures N ; therefore estimating population size essentially boils down to estimating the capture probability, which can then be used to inflate the observed N via the estimator Thus we turn towards the crucial question of how to efficiently estimate the capture probability (specifically, in the presence of high-dimensional and/or complex covariates X), before coming back to inference about n in Section 5. In the next section we discuss identification of ψ from the distribution Q from which we sample; this will require extra assumptions, since if lists are irreparably dependent we will have no information about those who are unobserved.

Identification
As mentioned in the previous section, without additional assumptions, the observed data distribution Q of list membership among those on at least one list is completely uninformative about the capture probability ψ = P(Y = 0). A variety of assumptions have been used to identify this and related quantities in previous literature; broadly, there must be some lack of dependence across lists in order to identify and estimate the odds and thus the overall population size. The popular Petersen estimator (Petersen, 1896) assumed independence between K = 2 lists and Darroch (1958) extended it to assume independence across K > 2 lists. Fienberg (1972) assumed a log-linear model for the expected number of observations with each capture profile, with one parameter necessarily set to zero (typically the highest-order interaction term across lists). You et al. (2021) have presented a generalizable framework adaptable to various identification assumptions including log-linear model assumption and independence between two lists conditional on the remaining list(s) for K > 2 list case.
When one has access to not only list membership but also covariate information, conditional versions of these assumptions can be used (Chao, 1987;Huggins and Hwang, 2007;Sekar and Deming, 1949;Tilling and Sterne, 1999). Importantly, this allows heterogeneous capture probabilities that vary across units, and thus relaxes identifying assumptions; this is conceptually similar to how measured confounders are exploited in observational studies for causal inference (Hernan and Robins, 2019). Sekar and Deming (1949) studied the conditional case in the discrete and low-dimensional setup; the continuous case has been studied using parametric models by Alho (1990); Chao (1987); Huggins (1989); Pledger (2000); Pollock et al. (1990);Tilling and Sterne (1999). Burnham and Overton (1979); Huggins and Hwang (2007) used non-parametric jackknife estimator and bandwidth selection respectively. Analogous to the assumption of Tilling and Sterne (1999) for the two list case, our main identifying assumption for the K list case is that there is a known pair among the K lists which are conditionally independent. Without loss of generality, we let the lists be ordered so that the conditional independent pair is list 1 and 2: where Y k denotes the capture indicator variable for list k for k = 1, . . . , K.
Assumption 1 says that the chance of appearing on list 1 is the same regardless of list 2 membership, among those with the same measured covariate values, i.e., that the list indicators Y 1 and Y 2 are conditionally independent given X. This assumption can be viewed as an important relaxation of a more standard assumption of marginal independence, particularly when lists cover different parts of a population.
For example, consider a toy setup where there are two regions, equally populous. Suppose people who live in region A have a 90% chance of appearing on list 1 and a 10% chance of appearing on list 2, while people who live in region B have the reverse: a 10% chance of appearing on list 1 and a 90% chance on list 2. Thus list 1 tends to capture region A people, and list 2 tends to capture region B people. In this case, even if conditional independence holds (i.e., the chance of appearing on list 1, within each region, is the same regardless of whether you appear on list 2), the lists will not be marginally independent. Intuitively, the reasoning behind this is straightforward. Once we know you are on list 2, we have some additional information about your chances of appearing on list 1: you are more likely to live in region B, and so less likely to be on list 1. More specifically, the chances of appearing on either list are both 50% marginally, but the chance of appearing on list 1 given that you are on list 2 is only 18%.
The conditional independence in Assumption 1 has been used relatively extensively in the capturerecapture literature; we refer to Alho et al. (1993);Brenner (1995); Huggins (1989); Pollock et al. (1990);Tilling (2001); Tilling and Sterne (1999) for more details and discussions of when this assumption may hold and when it may fail.
It is known (e.g., as in Tilling and Sterne, 1999) that under Assumption 1 the capture probability ψ = P(Y = 0) can be identified from the biased observed data distribution Q. Specifically, let denote the observational probability (under Q) of appearing on list 1, 2, and both, respectively. These probabilities will be referred to as the q-probabilities at various points throughout.
Remark 4. Note that when there are only K = 2 lists, it must be that q 1 (x) + q 2 (x) − q 12 (x) = 1 since each observed unit must appear on list 1, list 2, or both, according to the sampling distribution Q. In general when K > 2 one only has the requirement 0 ≤ q 1 (x) + q 2 (x) − q 12 (x) ≤ 1, since it is possible some individuals only appear on lists j ≥ 3 other than 1 and 2.
Remark 5. Note that when there are more than two lists and without loss of generality the conditionally independent list pair is 1 and 2, the remaining lists are aiding in the estimation by potentially increasing the number of observed individuals. This in turn leads to variance reduction as is discussed in appendix A.1.
For posterity we give the identification result for ψ in the following proposition (with a proof given in the Appendix).
Proposition 1. Under Assumption 1 and the positivity condition Q{q 12 (X) > 0} = 1, the conditional and marginal capture probabilities are identified from Q by In the following sections we give three main contributions. First we derive the efficiency bound for estimating the capture probability ψ under a nonparametric model that puts no parametric restrictions on the 'nuisance' functions (q 1 , q 2 , q 12 ); second, we construct novel estimators that attain the efficiency bound under weak nonparametric conditions (e.g., allowing the use of flexible machine learning tools); and third, we give a general method for building corresponding confidence intervals for the total population size n, given any asymptotically linear estimate ψ of the capture probability.
Remark 6. All subsequent results apply to the statistical parameter ψ, which we define from here on as the harmonic mean on the right-hand-side of (6). Under the identifying Assumption 1 (and positivity), ψ also represents the capture probability on the left-hand-side of (6), but our statistical results do not require this link, and apply to the harmonic mean in (6) regardless.
Remark 7. When there are more than two lists, it is possible that multiple list pairs satisfy assumption 1. This presents multiple different opportunities to identify and estimate the capture probability, and so yields a semiparametric model with testable implications. We leave exploration of this setup for potential future work.

Efficiency Bound
In this section, we derive the nonparametric efficiency bound for estimating the capture probability ψ using an iid sample from distribution Q. This gives a crucial benchmark against which one can compare candidate estimators: once an estimator is shown to attain this bound, no further improvements can be made (at least asymptotically) without adding extra assumptions. To the best of our knowledge, the only previous efficiency bounds in the capture-recapture literature are for low-dimensional parametric models, where standard results from maximum likelihood theory apply (Fienberg, 1972;Gimenez et al., 2005).
In order to derive the efficiency bound, we use tools from semiparametric theory (Bickel et al., 1993). A fundamental goal here is to characterize influence functions, and the efficient influence function in particular. The efficient influence function of a parameter acts as the derivative term in a distributional Taylor expansion of the parameter, viewed as a map on distributions; thus it can represent the change in the parameter after perturbing the distribution it takes as input. Practically, the efficient influence function for a parameter has several critical implications. First, as mentioned above, it leads to a minimax efficiency bound (van der Vaart, 2002a) and thus provides a benchmark for efficient estimation in flexible nonparametric models. Further, it can be used to construct efficient estimators that attain the bound under weak assumptions, and sheds light on the regularity conditions necessary for said efficiency, as will be shown in Section 4. More details on nonparametric efficiency theory can be found in Bickel et al. (1993), van der Vaart (2002a), and van der Laan and Robins (2003), among others; reviews can be found in Tsiatis (2006) and Kennedy (2016), for example.
Our first result gives the form of the efficient influence function for the capture probability, in an unrestricted nonparametric model. Lemma 1. Let g : R → R be any function differentiable at the true capture probability ψ defined in (6). Under a nonparametric model, the efficient influence function for the parameter g(ψ) is given by The proof is presented in the Appendix. Note that the first term in the efficient influence function is a product of the inverse conditional capture probability γ −1 with a term whose conditional expectation given X (under Q) equals 1. The efficient influence function will be bounded for example if q 12 (x) ≥ for some > 0 and all x (note that q 1 (x) ∧ q 2 (x) ≥ q 12 (x) so q 12 (x) ≥ implies all the q-probabilities are bounded below by ).
The variance of the efficient influence function acts as a nonparametric efficiency bound, in that no estimator can achieve a better mean squared error in a local minimax sense (van der Vaart, 2002b). The following theorem and corollary give the form of this bound and formalize the minimax result, respectively. All expectations and variances are under distribution Q unless noted otherwise.
Theorem 1. Let g : R → R be any function differentiable at ψ. The nonparametric efficiency bound for estimation of g(ψ) is given by var{f and q 0 (x) = 1 − q 1 (x) − q 2 (x) + q 12 (x) is the chance of appearing on neither list 1 nor 2.
The magnitude of the efficiency bound in Theorem 1 is driven by three main factors: (i) the magnitude of the conditional capture probabilities, (ii) the chance of appearing on both lists, and (iii) the heterogeneity in the conditional capture probabilities.
The dependence on (i) in the bound in Theorem 1 occurs through the term (1 − γ)/γ 2 , i.e., the odds of capture divided by the capture probability. The dependence on (ii) occurs through the odds (1 − q 12 )/q 12 as well as the probability ratio q 0 /q 12 . The dependence on the heterogeneity (iii) occurs through the var(1/γ) term. Note that the probabilities γ and q 12 in (i) and (ii) are related, but q 12 can be small even when the capture probability γ is not, depending on the size of q 1 and q 2 .
More specifically, all else equal, the variance bound increases with: (i) smaller capture probabilities γ, (ii) smaller chances of appearing on both lists q 12 , and (iii) greater heterogeneity in the capture probabilities γ. Therefore capture probabilities can be estimated most efficiently when capture is likely, when there is substantial overlap across lists, and when capture probabilities are more homogeneous.
Remark 9. For K = 2 lists the quantity q 0 (x) = 0 is exactly zero, but when K > 2 it can be positive.
In addition to informing what factors yield more or less efficient capture probability estimation, the variance in Theorem 1 also acts as a local minimax lower bound, as formalized in the following corollary.
Corollary 1. For any estimator g( ψ), it follows that where TV(Q, Q) is the total variation between Q and Q, ψ = ψ(Q) and ψ = ψ(Q) are the capture probabilities under Q and Q, respectively, and f g (ψ) and σ 2 = σ 2 (Q) are defined as in Theorem 1.
Corollary 1 shows that the worst-case mean squared error of any estimator of ψ, locally near the true Q, cannot be smaller than the efficiency bound, asymptotically and after scaling by N . This local minimax result gives an important benchmark for efficient estimation of the inverse capture probability: no estimator can have mean squared error uniformly better than the variance of the efficient influence function divided by N , without adding extra assumptions and/or structure to the nonparametric model we consider.
Remark 11. The local minimax result in Corollary 1 holds for any subconvex loss function : Importantly, in the next section we construct estimators that can achieve the nonparametric efficiency bound under weak conditions that allow for flexible estimation of the q-probabilities, e.g., using machine learning tools.

Setup
Recall we let Q N denote the empirical measure under Q, so that sample averages can be written with the short-hand The simplest estimator of the capture probability ψ is just a plug-in which replaces unknown quantities in the definition of ψ with estimates, i.e., by estimating the conditional capture probability γ(X) = q12(X) q1(X) q2(X) for every unit, and computing the harmonic mean of the values across the sample. This estimator has been used relatively extensively in previous work (Darroch, 1958;Fienberg, 1972;Huggins and Hwang, 2007;Tilling and Sterne, 1999). When the qprobabilities are estimated with correctly specified parametric models, the plug-in estimator ψ pi will be √ n-consistent and asymptotically normal under standard regularity conditions. However, when the covariates contain any continuous components and/or are high-dimensional, it is usually very unlikely an analyst would have enough a priori knowledge to be able to correctly specify a low-dimensional parametric model, let alone three (one for each q-probability nuisance function).
This difficulty of correct model specification suggests trying to flexibly estimate the q-probabilities, e.g., using logistic regression with model selection, or lasso, or nonparametric tools like random forests, neural nets, RKHS regression, etc. Unfortunately, when the plug-in estimator ψ pi is constructed from these kinds of data-adaptive methods, it in general loses the nice properties it has in the parametric setup. Specifically, without special tuning of particular methods, it in general would suffer from slower than √ n-convergence rates, and have an unknown limiting distribution, making it not only inefficient but also leaving no tractable way to do inference. This deficiency of plug-in estimators is by now relatively well-known in functional estimation problems (Chernozhukov et al., 2018;van der Laan and Robins, 2003;Wu et al., 2019); however we have not seen it highlighted in the capture-recapture setting. (We show some of these issues numerically with simulations in Section 6.1.) Luckily, the plug-in can be improved upon using tools from semiparametric efficiency theory (Bickel et al., 1993;Kennedy, 2016;Tsiatis, 2006;van der Laan and Robins, 2003;van der Vaart, 2002a). In what follows, we will present and study a novel doubly robust estimator, which can attain the efficiency bound from the previous section even when built from flexible data-adaptive regression tools.

Doubly Robust Estimator
As mentioned above, the plug-in estimator (8) has some important deficiencies in semi-and nonparametric settings. The plug-in estimator can be debiased by adding an estimate of the mean of the efficient influence function (Bickel et al., 1993;Kennedy, 2016;Tsiatis, 2006;van der Laan and Robins, 2003;van der Vaart, 2002a). This leads to our proposed doubly robust estimator where q j are estimates of the q-probabilities (e.g., via regression predictions).
Remark 12. In order to avoid potentially restrictive empirical process conditions, we estimate Q N and q j from separate independent samples. Specifically, we estimate the q-probability nuisance functions by fitting regressions in a training sample, independent of a test sample Q N . With iid data, one can always obtain such samples by splitting at random in half, or folds. This yields a loss in efficiency, but that can be fixed by swapping the samples/folds, computing the estimate on each, and averaging. This is referred to as cross-fitting, and has been used for example by Bickel and Ritov (1988); Chernozhukov et al. (2017); Robins et al. (2008); Zheng and van der Laan (2010). Here we analyze a single split procedure, merely to simplify notation; extending to averages across independent splits is straightforward.
In the following sub-section, we derive finite-sample error bounds and distributional approximations for our doubly robust estimator, which are valid for any sample size.

Non-asymptotic Error Bounds and Approximate Normality
In this section, we provide our three main theoretical results regarding error bounds for our proposed method. In particular we show that our estimator is nearly efficient, doubly robust, and approximately normal. Importantly, we show all these properties hold in finite samples, without resorting to asymptotics.
In the previous section, we derived the efficient influence function, which is the crucial component of the local minimax lower bound we gave in Corollary 1. This corollary shows the minimax optimal estimator has mean squared error that scales like the variance of the efficient influence function divided by N , so that an optimal estimator would be one that behaves like a sample average of the efficient influence function. Our first result shows that our proposed estimator does in fact behave like an average of the efficient influence function, depending on the size of nuisance error. In what follows, we use φ and φ to denote the efficient influence function φ(Z; Q) and its estimate φ(Z; Q) respectively.
Theorem 2. For any sample size N and error tolerance δ > 0, we have where R 2 is a second-order error term given by with the latter bound on R 2 holding as long as (q 12 ∧ q 12 ) ≥ .
Theorem 2 shows that our proposed estimator is within δ of a sample average of the efficient influence function, centered at the true (inverse) capture probability, with high probability, at every sample size. For a given observed number of captures N and error δ, this probability depends on two factors: (i) a second-order error term R 2 , which is driven by the error in estimating the nuisance q-probabilities; and (ii) the L 2 error in estimating the efficient influence function itself, divided by N , which also depends on estimation error of the q-probabilities, but in a weaker way due to the division by N . When R 2 goes to 0 as N increases, the probability above goes to 1 for any fixed δ. Hence, Theorem 2 implies usual asymptotic convergence in probability, but in addition it gives an error bound that is valid for any finite N .
For example, if the q-probabilities are estimated with errors upper bounded by cN −1/4 , then with at least 95% probability, our proposed estimator will be within 4c 2 5/N of the average efficient influence function. More generally, if E| R 2 | 1/ √ N , then our proposed estimator will be within 1/ √ N (up to constants) of this efficient average, with high probability. For example, if q 1 , q 2 , q 12 and γ belong to Holder classes H(β 1 ), H(β 2 ), H(β 12 ) and H(β 0 ) respectively, where H(s) is a Holder class with smoothness index s (Györfi et al., 2006;Tsybakov, 2008), then a sufficient condition for this kind of result, if the q-probabilities are estimated at minimax optimal rates, would be that the minimum smoothness is at least half the dimension of the covariates, i.e., min{β 0 , β 12 , β 1 , β 2 } ≥ d/2. Similarly, if the q-probabilities were s-sparse, then a sufficient condition would be that s √ n with lasso-style methods (Farrell, 2015). However, such N −1/4 nuisance errors are only sufficient conditions for 1/ √ N capture probability errors; one only needs the remainder error R 2 to be small enough, which could also be achieved if some combinations of q-probabilities are estimated well, even if others are not. We give more detail on this phenomenon in our next result.
Namely, beyond being close to an optimally efficient estimator, in the next result we show that our proposed estimator enjoys a finite-sample multiple robustness phenomenon, never before shown in capture-recapture problems. This phenomenon indicates that the overall estimation error can be small as long as some, but not all, nuisance probabilities are estimated with small error.
Corollary 2. Suppose (q 12 ∧ q 12 ) ≥ . Assume one of the following holds: As a consequence of Corollary 2, the proposed estimator is doubly robust, i.e., if either estimator of q 1 or q 2 has small error, and either estimator of q 12 and γ has small error, then the overall error of our proposed estimator (given by R 2 ) will be just as small, up to constants, even if the other estimators have large errors or are misspecified. (We note that, although this kind of robustness is sometimes called multiple robustness (Vansteelandt et al., 2008), we use the term doubly robust since the error structure is still second-order, i.e., involving products of errors, albeit with more terms). This property is very useful when one of the lists is difficult to estimate, for example, due to high-dimensional covariates, or q-probabilities that are complex functions of continuous covariates.
Remark 13. Note that when there are only K = 2 lists, then we have the relation q 1 (x)+q 2 (x)−q 12 (x) = 1 for all x. Hence, small bias for any two estimators of q 12 , q 1 and q 2 automatically implies small bias for the third. One might expect then that double robustness does not arise in the K = 2 list setting; however this is not quite true. To see why, note that it could be possible to estimate γ with small error, for example, even when some of the q-probability estimators are misspecified. These and other issues related to estimation of the conditional capture probability γ will be important to explore in future work.
When the remainder error R 2 is sufficiently small, Theorem 2 and Corollary 2 tell us we can approximate ψ −1 dr with a sample average of the efficient influence function. For the purposes of inference, this suggests a confidence interval of the form where σ is a variance term defined in Theorem 3. In the next Berry-Esseen-type result, we exploit this closeness with a sample average and further show that our proposed estimator, properly scaled, is approximately Gaussian. This will show that the above confidence interval gives nearly-valid finitesample coverage guarantees.
Theorem 3. Let σ 2 = var( φ) be the unbiased empirical variance of the estimated efficient influence function. Then ψ −1 dr − ψ −1 follows an approximately Gaussian distribution, with the difference in cumulative distribution functions uniformly bounded above by The above result shows that the estimation error scaled by σ/ √ N is approximately standard normal. The first term on the right hand side of (11) is the usual Berry-Esseen bound. The second term captures the effect of the nuisance estimation error R 2 . The third term is the estimation error in the variance. Since E| σ − σ| is bounded above by cN −1/2 (proof in the appendix), the overall error in the Gaussian approximation is driven by the second term, involving nuisance error R 2 . This will be the main driver of whether the interval has approximately correct coverage. We note that the above theorem implies convergence in distribution whenever E| R 2 | = o(1/ √ N ) (which can hold for a wide variety of flexible nonparametric estimators of the q-probabilities, as discussed after Theorem 2), but in addition gives a more precise error bound that holds for any finite sample size.
Note that Theorem 3 immediately implies that the error in coverage for the proposed confidence interval defined in (10) is no more than twice the error bound on the right-hand-side of (11), with t = z α/2 . Further, a Berry-Esseen-style bound similar to that of Theorem 3 (along with subsequent coverage guarantees and corollaries) can be obtained for any function g(·) of ψ −1 dr satisfying the conditions from Friedrich (1989). This implies the same kind of coverage guarantees for ψ, for example, using the confidence interval which can be motivated via the delta method. The error in the coverage of this estimated interval is twice the bound in Theorem 3, modulo some extra dependence on g. Importantly, the unbiased empirical variance σ 2 is a consistent estimator of the efficiency bound σ 2 = var(φ) in the sense that E|σ − σ| E φ − φ + N −1/2 . This shows the crucial result that our estimator is approximately minimax optimal in the sense of Corollary 1, if the nuisance error is small enough.
A natural consequence of the above theorem is the following corollary, which presents a simple bound on the error of the normal approximation, under some natural conditions on the nuisance error R 2 and variance.
Corollary 3. Assume σ 1, E| R 2 | N −2β and α > δ for some δ > 0. Then the coverage error for the proposed (1 − α) confidence interval defined in (10) is upper bounded by Therefore if β > 1/4 there exists some sample size N at which the coverage error is never more than , for any N > N .
Since this corollary is a special case of Theorem 3, mainly aimed at presenting the result in a simple form, we refer to the above discussion for more details. However we note that the condition that E| R 2 | N −2β would hold for example if the q-probabilities were estimated optimally when contained in Holder classes with smoothness index s, where β = s 2s+d (or under some conditions on sparsity, as discussed after Theorem 2). Then β > 1/4 would mean s > d/2, aligning with our earlier results.
In this section, we have given finite-sample error bounds and distributional approximations for our proposed estimator, which are valid for any sample size, allowing accurate estimation and approximately valid confidence guarantees, even in complex nonparametric models where the q-probabilities are estimated with flexible machine learning tools. In the next section, we consider a slightly modified version of the estimator which could further improve finite-sample properties.

Targeted Maximum Likelihood Estimator
We have seen that our proposed doubly robust estimator (9) is close in a finite-sample sense to an optimal sample average, and possesses crucial double robustness properties. However it is possible this estimator may not respect the bounds on the parameter space; for example ψ dr may fall outside [0, 1] if some of the estimates of the q-probabilities are small. A simple fix is to truncate the estimator ψ dr to always lie in [0, 1]. Here for completeness we discuss an alternative approach using targeted maximum likelihood estimation (TMLE) (van der Laan and Rose, 2011;van der Laan and Rubin, 2006), which is an iterative procedure that fluctuates nuisance estimates so that a plug-in estimator built from them also approximately solves an efficient influence function estimating equation. TMLE thus leads to estimators that are asymptotically equivalent to one-step bias-corrected estimators, but which could bring some finite-sample advantages.
Below we present an algorithm detailing the computation of a TMLE for ψ. At a high level, the procedure involves iterative updating of initial nuisance estimates, based on the quantities H 12,t , H 1,t , and H 2,t defined shortly, which are called clever covariates in the TMLE literature.
, where q * j are estimates obtained after convergence.
Interestingly, since the TMLE above is not a sample average like our main proposed estimator from the previous subsection, it is less clear how to derive finite-sample error bounds, despite this being the main motivation for TMLE. Since the estimates q * j (x) obtained after convergence satisfy Q N {φ(Z; Q * )} ≈ 0, the asymptotic behavior matches the doubly robust estimator in (9), but for describing finite-sample behavior we resort to simulations, detailed in Section 6.1.

Inference for Population Size
In the previous section we gave doubly robust estimators for the capture probability and studied finitesample properties. In this section we give a crucial result that shows how to obtain an approximate confidence interval for the population size, given a generic initial estimator of the (inverse) capture probability. Importantly, our results only require this initial estimator to be weakly approximated by a sample average, and otherwise are completely agnostic to how the capture probability is estimated. This appears in stark contrast to most of the literature on this topic, where the inferential procedures are very closely tied to specific model assumptions and estimator constructions.
This main inferential result is given in the following theorem.
Theorem 4. Suppose we are given an initial estimator ψ that satisfies for ϕ a generic influence function with mean zero and R 2 an error term. Let τ 2 = ψ ς 2 + 1− ψ ψ and τ 2 = ψ ς 2 + 1−ψ ψ (ψ R 2 + 1) 2 , where ς 2 = var( ϕ) is the unbiased empirical variance of the estimated influence function and ς = var( ϕ | Z n ) the true conditional variance. Then the (1 − α) confidence interval given by has coverage error upper bounded as where C is the Berry-Esseen constant and Theorem 4 gives a non-asymptotic upper bound on how much the coverage P( CI n n) of our proposed interval can deviate from its nominal (1 − α) level. Before describing the coverage guarantee, we first describe the proposed interval. The length of this interval is driven by three factors: (i) the estimated odds of not being captured (1 − ψ)/ ψ, (ii) the variance of the inverse capture probability estimator ς 2 , and (iii) the sample size N . As one would expect, higher odds of capture yield more precise inference about population size, all else equal, as does more efficient estimation of ψ. Specifically, the length of the interval shrinks to zero when the capture probability is very large, regardless of the sample size N . Also note that even if ψ were known, one would still have an interval of the form based on the fact that n = N/ √ n is approximately normal. Although sample size N appears in the numerator of the interval width (contrary to standard intervals), it only appears through its square root, showing that in an asymptotic regime where N → ∞, the width still grows at a slower rate than the sample size. Intuitively, one can think of this interval as taking n = N/ ψ and multiplying by 1 ± z α/2 / √ n, which does tend to zero as sample size N grows.
Remark 15. For K = 2 lists and in the absence of covariates, the confidence interval reduces to n ± z α/2 n(1− ψ) ψ q12 , which approximately resembles the Wald-type confidence interval for the Petersen estimator (Evans et al., 1996). Now we describe the coverage guarantee of Theorem 4. Importantly, the bound on the coverage error depends on a number of factors, as shown above appearing the sum of the three terms in (13). Under typical boundedness assumptions, the first and third terms would be of smaller order, and the second term would dominate. This second term is driven by the size of R 2 in terms of its mean absolute value, i.e., how well the initial estimator ψ is approximated by a sample average. If R 2 is not substantially smaller than 1/ √ n, then the confidence interval would not be guaranteed to cover the true population size n at its nominal level. This points to the importance of efficient estimation of ψ; for example, as shown in the previous section, our proposed estimator ψ dr can be approximated by a sample average up to smaller than 1/ √ n error, even in a nonparametric model when q-probabilities are estimated flexibly.
Remark 16. A unique feature of Theorem 4 is that it is valid for any estimator approximated by a sample average, regardless of what underlying identification or estimation assumptions were used in its construction. This means if another analyst did not believe the independent lists condition in Assumption 1, and instead constructed an estimate of the capture probability under a different identifying assumption, they could also use the above theorem to construct a confidence interval and assess its finite-sample coverage.
A natural consequence of Theorem 4 is the following corollary, which parallels Corollary 3 in giving a simple bound on the error of the normal approximation, under some natural conditions. Corollary 4. Assume τ 1, E| R 2 | N −2β and α > δ for some δ > 0. Then the coverage error for the proposed (1 − α) confidence interval defined in (10) is upper bounded by P CI n n − (1 − α) n (1−4β)/2 + 1 √ n .
Therefore if β > 1/4 there exists some population size n at which the coverage error is never more than , for any n > n .
Since the result in Corollary 4 is similar to that of Corollary 3, we refer there for related discussion. The main point is that, as long as our initial estimator is well-approximated by a sample average, no matter how it was constructed or what assumptions it relies on, our proposed confidence interval (12) will be approximately valid.

Simulation & Application
So far we have proposed doubly robust estimators for the capture probability, and a general approach for constructing confidence intervals for the total population size, all with non-asymptotic error guarantees. In this section we study the performance of our methods in simulated data, and apply them to estimate the total number of killings in the internal armed conflict in Peru during 1980-2000. The code used to generate the results is available at mqnjqrid/capture_recapture.

Simulation
Here we use a simulation set-up similar to that of Tilling and Sterne (1999). Specifically we take n = 5000 samples from where a takes values {−2.513, −1.758, −0.66} to ensure that the capture probability ψ takes values {0.3, 0.5, 0.8}, respectively. This gives sample sizes N approximately equal to {1500, 2500, 4000}.
Recall that under P, list membership Y 1 and Y 2 are conditionally independent, so the conditional capture probability is and q-probabilities are equal to q j (x) = P(Y j = 1 | X = x)/γ(x).
We construct estimates of the q-probabilities via q j (x) = expit[logit{q j (x)} + j ], where ∼ N (n −α , n −2α ). This allows us to carefully control the error of the q-probability estimators; since the root mean squared error scales like n −α , this can be viewed as the rate of convergence. We run 500 simulations for α ∈ {0.1, 0.2, 0.25, 0.3, 0.4, 0.5}. Note α values 0.5 and 0.25 correspond to the parametric (n −1/2 ) and nonparametric (n −1/4 ) rates, respectively. Figure 2 shows the estimated bias and the root mean square error (RMSE) of ψ, along with the coverage proportion for the confidence interval of the total population size.
Remark 17. For plug-in estimators, there is no well-defined variance formula (this is a main motivation for our doubly robust construction). Therefore to construct confidence intervals with the plug-in estimator, we used the estimated variance of the doubly robust estimator.
Overall, the simulations illustrate the phenomena expected from our theoretical results: when the q-probabilities are estimated with low error (i.e., α large), all the methods do well, whereas when the q-probabilities are difficult to estimate (i.e., α smaller) the proposed methods do substantially better in terms of bias, error, and coverage. For example, when the true capture probability is 50%, the simple plug-in estimator gives substantial bias as soon as α < 0.4 (i.e., when the q-probabilities are estimated at slower than n −2/5 rates). However, the bias of the proposed doubly robust estimator is relatively unaffected until α < 0.2, with the TMLE somewhere in between. The story is similar for the RMSE, which is largely driven by the bias in this problem. The coverage is approximately at the nominal 95% level as soon as α > 0.2 (i.e., when the q-probabilities are estimated at faster than n −1/5 rates), whereas the plug-in estimator substantially under-covers (e.g., nearly zero at α = 0.2) until α ≥ 0.4. We note that when population size or capture probability is small (e.g., capture probability substantially less than 50%), estimation becomes more challenging and the story is less clear about which method does better. For reference, results for population sizes varying from n = 5000 to n = 25000 (in the α = 0.25 case) are given in the Appendix in Figure 6.

Data Analysis
We apply our proposed methods to estimate the number of killings and disappearances attributable to different groups in Peru during its internal armed conflict between 1980 and 2000. We use data collected by the Truth and Reconciliation commission of Peru (Ball et al., 2003), as well as detailed geographic information, following Rendon (2019a).
There is an ongoing debate regarding the total number of killings and disappearances in the conflict, as well as about which groups are most responsible, e.g., the PCP-Shining Path versus the State or other groups. Ball et al. (2003) estimated approximately 69,000 total killings and disappearences, finding the Shining Path responsible for the majority. In contrast, Rendon (2019a) estimated approximately 48,000 killings and disappearances, with the State responsible for the majority, though many geographic strata were excluded. Most recently, Manrique-Vallier and Ball (2019) included a newly available list and estimated approximately 58,000-65,000 killings and disappearances, depending on choices of priors, with the Shining Path responsible for the majority.
Before describing our specific approach, we first describe the data and give some summary statistics. As explained in Ball et al. (2003), the data come from a few main sources: the Truth and Reconciliation Commission (CVR), the Public Defender Office (DP), and 4-5 other human rights groups and NGOs (ODH); we combined the DP and ODH lists since they have similar demographics. The data contains identifiers of people who have been killed or disappeared, as well as which of the source lists they appeared on, and covariates including age, gender, and geographic location of the killing or disappearance (measured via 58 geographic strata as in Ball et al. (2003), as well as bivariate latitude/longitude as in Rendon (2019a)). To avoid missing completely at random assumptions, we also included missingness indicators for victims with missing age (28% missing), gender (<1% missing), or location (11% missing). The lists of all the covariates is available in the appendix C. The total number of killings and disappearances across all lists was 24,692. Importantly, the lists capture different demographics, which points to the necessity of relaxing classical marginal independence via the conditional independence in Assumption 1. For example, the CVR list mostly includes victims who were killed, while the DP and ODH lists mostly include victims who disappeared, as shown in Figure 3. Similarly, geographic diversity varies across lists, as shown in Figure 4. For example, almost 60% of Shining Path victims in the DP and ODH lists come from two smaller districts (Chungui and Luis Carranzo) of Ayacucho, while in the CVR list the Shining Path victims are more uniformly spread across the country. More details on the data are available in Appendix C. Now we move to our analysis. Our goal was to estimate the number of killings and disappearances attributable to the State and Shining Path, as well as those that were not identified as either. We used our proposed doubly robust estimator (9) with five-fold cross-fitting, and we estimated the qprobabilities via random forests (using the ranger package in R). We truncated all q-probability estimates at 0.01. Figure 5 shows the estimated number of killings and disappearances along with 95% confidence intervals obtained using the interval (12). We estimate the total number of killings and disappearances across groups to be 68,874 (95% CI: 58,204), close to the estimates in Ball et al. (2003) and the diffuse prior-based estimate in Manrique-Vallier and Ball (2019) (which used an additional list). Overall we find the State responsible for more disappearances, and Shining Path responsible for more killings; however we estimate the number of killings and disappearances by unidentified perpetrators to be larger than that for either group. In terms of the overall killings and

Discussion
In this paper, we study estimation of population size and capture probability in the capture-recapture set-up where two lists are conditionally independent given measured covariates. We make four main contributions to the literature. First, we derive the nonparametric efficiency bound for estimating the capture probability, which indicates the best possible performance of any estimator, in a local asymptotic minimax sense. As far as we know this kind of lower bound result has not appeared in the literature, even in simple settings without covariates. Second, we present a new doubly robust estimator, and study its finite-sample properties; in addition to double robustness, we show that it is near-optimal in a non-asymptotic sense, under relatively mild nonparametric conditions. Third, we give a method for constructing confidence intervals for total population size from generic capture probability estimators, and prove non-asymptotic near-validity. And fourth, we study our methods in simulations, and apply them to estimate the number of killings and disappearances attributable to different groups in Peru during its internal armed conflict between 1980 and 2000.
There are many ways one could extend and build on the work in this paper. For example, rather than assuming a known pair of lists are conditionally independent given the covariates, one could instead take a sensitivity analysis and/or partial identification approach. For example, one could assume that a pair of lists is only nearly conditionally independent, up to some deviation δ, and estimate bounds on the capture probability and population size accordingly. This relies on weaker assumptions, with the trade-off of yielding less precise inferences. Another extension would be to flexibly estimate conditional capture probability or population size, given a continuous covariate such as age or time. For example, for the internal armed conflict in Peru it might be of interest to estimate the number of victims by age. This would require a non-trivial extension of the current methods, but would be important future work.

A Proofs of theorems and results
Proof of proposition 1. The goal is to express ψ in terms of the observed data distribution Q. By definition ψ = P(Y = 0) and γ(x) = P(Y = 0|X = x).
First we show that ψ is the harmonic mean of γ(x). It is easy to see that . Thus, We have shown a that the capture probability ψ is the harmonic mean of the conditional capture probabilities γ(x) under the observed data distribution Q.
Next, we express γ(x) in terms of the q-probabilities. Under assumption 1, for any fixed covariate value x, . For a capture history y = 0, Substituting the above relation in the assumption statement and using the notation of the q-probabilities, we get Proof of lemma 1. There are several ways one can derive an efficient influence function and corresponding efficiency bound. Here, we first find a putative influence function in the discrete case using a standard Gateaux derivative argument. Then we show that this influence function is actually the efficient one in a general nonparametric model (continuous or discrete or mixed), by checking that the corresponding remainder term in a von Mises expansion is second-order. To find a candidate influence function, we consider a special parametric submodel (i.e., deviation from Q) given by Q with density q = (1 − )q(z) + q(z) whereq =q(z) = 1(z =z) is a point mass at Z =z, and for which the pathwise derivative actually equals the influence function (in the discrete case) Hampel (1974); Mises (1947). We also let q s, (x) denote the analog of q s (x) under the submodel for s ∈ {1, 2, 12}, e.g., the marginal density for X under Q is q (x) = y q (z) = (1 − )q(x) + 1(x =x).

Now the above pathwise derivative equals
where the last line follows by the product rule. For the discrete case, we use the notation of the integral to denote summation over x.
For the derivatives appearing above, by the definition ofq, we have q (x) = ∂ ∂ q (x) = 1(x = x) − q(x). Similarly, by using derivative of product rule for q 1, (x), we get The last step follows from the product rule of derivatives. Finally, setting = 0, we get q 1 (x)| =0 =

1(x=x)
q(x) {Ỹ 1 − q 1 (x)}. The derivatives for q 2, and q 12, follow similarly. Thus, combining the above results and using the discrete nature of the distribution, we get Now that we have a candidate influence function that is valid in a discrete model, we evaluate the remainder term R 2 in a general submodel, to show that it is actually the efficient influence function in the general case as well.
Lettingψ = ψ(Q) for a generic distributionQ, the remainder of the so-called von Mises expansion (Bickel et al., 1993;van der Laan and Robins, 2003) is then given by: The above shows that the remainder of the von Mises expansion is second-order, i.e., involving products of differences between components of Q andQ. This fact implies the more general pathwise differentiability condition that for any smooth parametric submodel Q (Kennedy, 2020). Thus, the candidate influence function satisfies the above pathwise differentiability condition and hence, is an efficient influence function. Moreover, since the model is non-parametric, φ is the only efficient influence function (Bickel et al., 1993;Tsiatis, 2006;van der Vaart, 2002b).
Proof of Theorem 1. We will evaluate the variance of the efficient influence function φ ≡ φ(Z; Q). Recall, Z = (Y, X). Hence, we use the law of total variance by conditioning on X.
The last step follows using 1 ψ − 1 ψ = − φ(z)dQ(z) + R 2 from the proof of lemma 1. The formula of R 2 is presented in the proof of lemma 1.
For simplicity of notation in the proof we will use Q(·) = E(·|Y = 0, Z n ) = (·)dQ(z) to denote the observed population average function conditioned on the training sample Z n . Using that fact that Qφ = 0, we get We want to bound E N by presenting a bound on E(E 2 N ). We will show that the expected value of E is expected value of the error term R 2 .
Next, we will show that the second quantity is 0.
The second equality follows because Q N is evaluated on i.i.d. terms and Q( φ − φ) = E( φ − φ|Z n ) by definition. Next, we evaluate the variance of E N .
var(E N ) = E{var(E N |Z n )} + var{E(E N |Z n )}, by the law of total variance The last equality follows because R 2 and Q(·) are constants when conditioned on Z n . Now, Thus, combining the above results we get Hence, by using Markov's inequality combined with the above inequality gives the required result.
Proof of Theorem 3. The error in the estimation of ψ −1 is ψ −1 − ψ −1 . Following the proof of Theorem 2, it can be expressed as where Q φ denotes φ(z)dQ(z) as mentioned in the proof of Theorem 2. Let σ 2 = var( φ | Z n ), and note ψ −1 − ψ −1 − R 2 = (Q N − Q) φ is a sample average of a fixed function given the training sample Z n . Therefore by Berry-Esseen we have for any t and N that Now note by the mean value theorem, for some t n between t and t + σ where the second inequality used the fact that sup t Φ (t) ≤ 1/ √ 2π and the triangle inequality. Therefore This implies by iterated expectation that Therefore if σ 1 with probability one, E(ρ) < c for some c > 0, E| σ − σ| n −1/2 and E| R 2 | n −2β , then If β > 1/4, then there exists an N = N guaranteeing that the LHS is no more than .

Bound on E| σ − σ|
We will show that this bound is n −1/2 . For φ, we have defined the quantities σ 2 = var( φ | Z n ) and σ 2 = var( φ). We can further expand σ 2 as follows The second quantity σ 2 is the unbiased estimator of the variance of φ given the test sample i.e.
where the second equality follows from the unbiasedness property. We need the bound on | σ − σ| which can be expressed as a linear function of the absolute difference of the respective squares.

Thus, combining the results
Proof of Theorem 4. The estimate of population size, n = N/ ψ depends on two random quantities (i) the number of observations N and (ii) the estimate of the capture probability ψ.
Calculation of mean and variance of n First, we re-write N ψ as a sample average for ease of calculation. As mentioned in the proof of Theorem 2, we use Q ϕ to denote the population average conditioned on training sample Z n i.e., ϕ(z)dQ(z).
Thus, n − n − nψ R 2 is a sample average nP n ζ. Moreover, when we condition on the training sample Z n , the ζ's are i.i.d. We present the conditional mean and variance of ζ below.
By simple algebra, we can bound the quantity on the right hand side above as follows.
The last inequality follows using the definition of τ and triangle inequality. To evaluate the third term, we will use the relation Similar to the proof of Theorem 2, The last bound follows by the relation between τ and ς. Thus, Combining the results above, we get the following bound.

B Result with simulated data for varying total population size
The plug-in, the doubly robust and the TMLE are applied on the simulated data from section 6.1 total population size varying from 5000 to 25000. We focus on the case α = 0.25 since it is the nonparametric convergence rate. The plots are presented in figure 6. For each combination of (ψ, n), we simulated a dataset 500 times. The estimated bias and RMSE times √ n shows that the doubly robust and the targeted maximum likelihood estimators have convergence rate 1/ √ n. The plug-in estimator on the other hand, has a slower convergence rate. Similarly, the coverage of the total population size estimate for the plug-in estimator is much lower than the nominal coverage of 0.95 compared to the proposed methods' coverage.

C Peru Internal Conflict Data 1980-2000
The data is collected by the Truth and Reconciliation Commission of Peru (Ball et al., 2003). It was further expanded by Rendon (2019a) with the addition of geographical parameters (Rendon, 2019b). The original dataset contains the geographic location in the form of Peru's UBIGEO codes. Rendon (2019a) added the continuous geographical coordinates for each region. Ball et al. (2003) divided Peru into 58 stratas using the UBIGEO codes. We used the approach of Rendon (2019a) and calculated the latitude, longitude and area for a region (department or strata) using shape files from PERÚ (2014). For the latitude and longitude of a region, we averaged the latitude and longitude of the border of that region as in Rendon (2019a).
The victims that have no assigned strata are discarded by Ball et al. (2003). We present the statistics of these victims in List of covariates used to model the nuisance functions are as follows: • age: numeric variable and takes 0 for missing age.
• indicator for non-missing age: takes value 1 if age is present and 0 otherwise.
• gender: has levels male, female and others (including missing).
• situation: whether the victim was killed or disappeared.
• perpetrator: four levels indicating whether the individual is a victim of the State, Shining Path, others, or unidentified groups.
• indicator for non-missing department information: takes value 1 if department information is available for the individual and 0 otherwise.
• department latitude, longitude and area in hectares. For the individuals with missing department code, we use the average latitude, average longitude and median area of all the departments.  Figure 6: Estimated bias, RMSE (multiplied by √ n to study the rate), and population size empirical coverage, for simulated data with true population size n ∈ {5000, 10000, 15000, 20000, 25000}, across true capture probability ψ ∈ {0.8, 0.5, 0.3}, and q-probability error rate n −0.25 . • strata code: 59 possible levels. Details and construction of the 58 strata are available in Ball et al. (2003)). Those with missing strata take value 59.
• strata latitude, longitude and area in hectares. For the individuals with missing strata code, we use the average latitude, average longitude and median area of all the strata.
• indicator of non-missing strata code.

C.1 Results
We present the exact estimated number of victims using our proposed doubly robust estimation in Table 2. The difference in the estimated number of victims of the State and the Shining Path for the 25 departments and the seven geographic regions (Ball et al., 2003, see) are presented in Figure  7. The State has a significantly higher estimated number of victims in department Ayacucho and the Northern region. The Shining Path has higher estimated number of victims in departments Junin and Puno.