A new integrated likelihood for estimating population size in dependent dual-record system

TL;DR

This paper introduces a novel integrated likelihood approach for estimating population size in dependent dual-record systems, addressing nonidentifiability issues and outperforming existing Bayesian methods in simulations and real data applications.

Contribution

Develops a new integrated likelihood method for dependent DRS that is computationally efficient and less sensitive to prior assumptions, improving estimation accuracy.

Findings

Proposed method outperforms existing Bayesian approaches in simulations.

The estimator is invariant and less sensitive to prior choices.

Application to real data demonstrates practical utility.

Abstract

Efficient estimation of population size from dependent dual-record system (DRS) remains a statistical challenge in capture-recapture type experiment. Owing to the nonidentifiability of the suitable Time-Behavioral Response Variation model (denoted as $M_{tb}$) under DRS, few methods are developed in Bayesian paradigm based on informative priors. Our contribution in this article is in developing integrated likelihood function from model $M_{tb}$ based on a novel approach developed by Severini (2007, Biometrika). Suitable weight function on nuisance parameter is derived under the assumption of availability of knowledge on the direction of behavioral dependency. Such pseudo-likelihood function is constructed so that the resulting estimator possess some desirable properties including invariance and negligible prior (or weight) sensitiveness. Extensive simulations explore the better performance of our proposed method in most of the situations than the existing Bayesian methods. Moreover, being a non-Bayesian estimator, it simply avoids heavy computational effort and time. Finally, illustration based on two real life data sets on epidemiology and economic census are presented.

Tables

Table 1. Table 1: 2 × 2 2 2 2\times 2 data structure from Dual-record-System (DRS) with corresponding cell probabilities mentioned in [ ] and p ⋅ ⋅ subscript 𝑝 ⋅ absent ⋅ p_{\cdot\cdot} =1.

	List 2
List 1	In	out	Total
In	$x_{11}\left[p_{11}\right]$	$x_{10}\left[p_{10}\right]$	$x_{1\cdot }\left[p_{1\cdot }\right]$
Out	$x_{01}\left[p_{01}\right]$	$x_{00}\left[p_{00}\right]$	$x_{0\cdot }\left[p_{0\cdot }\right]$
Total	$x_{\cdot 1}\left[p_{\cdot 1}\right]$	$x_{\cdot 0}\left[p_{\cdot 0}\right]$	$x_{\cdot \cdot }=N\left[p_{\cdot \cdot }\right]$

Table 2. Table 2: Summary results of the three methods (Lee, SEMWiG & N ^ t ​ b I superscript subscript ^ 𝑁 𝑡 𝑏 𝐼 \hat{N}_{tb}^{I} ) applied to the populations P1-P6 with N = 200 𝑁 200 N=200 and when ϕ > 1 italic-ϕ 1 \phi>1 is known.

Method		P1	P2	P3	P4	P5	P6
$\varphi =1.25$
$Lee$222Prior on $\varphi$ is chosen as U($1,2$) since $\varphi >1$ is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	193(15.38)	191(12.33)	195(6.90)	199(9.71)	188(14.05)	204(14.56)
	C.I.	$(171,226)$	$(176,212)$	$(187,208)$	$(179,227)$	$(175,208)$	$(179,237)$
$SEMWiG$333Prior on $\varphi$ is chosen as U($1,p^{-1}$) since $\varphi >1$ is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	197(12.79)	194(10.59)	198(5.89)	204(10.50)	190(10.38)	225(22.89)
	C.I.	$(194,201)$	$(192,197)$	$(197,200)$	$(200,208)$	$(189,191)$	$(220,230)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	191(13.16)	195(9.22)	199(5.07)	200(6.69)	190(12.98)	198(8.99)
	C.I.	$(173,208)$	$(181,209)$	$(188,211)$	$(186,215)$	$(174,202)$	$(183,215)$
$\varphi =1.50$
$Le{e}^a$	$\hat{N}$(RMSE)	173(28.33)	178(23.11)	189(11.48)	187(15.64)	174(26.09)	190(14.91)
	C.I.	$(159,196)$	$(167,194)$	$(183,201)$	$(172,209)$	$(166,189)$	$(172,216)$
$SEMWi{G}^b$	$\hat{N}$(RMSE)	176(25.80)	179(21.94)	191(10.53)	192(12.40)	173(26.98)	194(12.41)
	C.I.	$(174,178)$	$(178,181)$	$(189,193)$	$(190,194)$	$(172,175)$	$(192,197)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	175(26.99)	180(21.66)	191(10.56)	191(11.85)	175(26.54)	192(11.00)
	C.I.	$(157,192)$	$(166,195)$	$(181,202)$	$(177,208)$	$(160,188)$	$(177,210)$

Table 3. Table 3: Summary results of the three methods (Lee, SEMWiG & N ^ t ​ b I superscript subscript ^ 𝑁 𝑡 𝑏 𝐼 \hat{N}_{tb}^{I} ) applied to the populations P1-P6 with N = 500 𝑁 500 N=500 and when ϕ > 1 italic-ϕ 1 \phi>1 is known.

Method		P1	P2	P3	P4	P5	P6
$\varphi =1.25$
$Lee$444Prior on $\varphi$ is chosen as U($1,2$) since $\varphi >1$ is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	472(34.11)	478(24.33)	490(12.02)	488(17.03)	504(13.71)	527(32.80)
	C.I.	$(438,518)$	$(451,519)$	$(473,515)$	$(457,531)$	$(445,576)$	$(471,599)$
$SEMWiG$555Prior on $\varphi$ is chosen as U($1,p^{-1}$) since $\varphi >1$ is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	484(23.92)	485(21.25)	495(9.87)	502(16.32)	480(24.64)	514(22.89)
	C.I.	$(451,522)$	$(456,511)$	$(477,512)$	$(474,532)$	$(473,493)$	$(505,521)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	478(26.75)	486(19.60)	498(9.06)	497(13.85)	476(27.19)	496(13.35)
	C.I.	$(451,508)$	$(460,511)$	$(480,514)$	$(471,521)$	$(451,500)$	$(471,519)$
$\varphi =1.50$
$Le{e}^a$	$\hat{N}$(RMSE)	437(64.17)	455(47.37)	481(21.58)	483(28.67)	44556.40	492(20.01)
	C.I.	$(405,487)$	$(424,500)$	$(452,522)$	$(437,538)$	$(419,488)$	$(449,550)$
$SEMWi{G}^b$	$\hat{N}$(RMSE)	442(59.56)	450(51.53)	477(22.66)	476(27.53)	436(64.69)	480(25.67)
	C.I.	$(436,452)$	$(447,456)$	$(475,582)$	$(470,482)$	$(433,439)$	$(475,487)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	442(60.52)	455(47.24)	480(21.69)	479(23.62)	439(62.40)	479(24.25)
	C.I.	$(415,469)$	$(432,478)$	$(464,496)$	$(455,501)$	$(415,463)$	$(452,503)$

Table 4. Table 4: Summary results of the three methods (Lee, SEMWiG & N ^ t ​ b I superscript subscript ^ 𝑁 𝑡 𝑏 𝐼 \hat{N}_{tb}^{I} ) applied to the populations A1-A6 with N = 200 𝑁 200 N=200 and when ϕ < 1 italic-ϕ 1 \phi<1 is known.

Method		A1	A2	A3	A4	A5	A6
$\varphi =0.60$
$Lee$666Prior on $\varphi$ is chosen as U($0.2,1.4$) since is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	236(43.96)	222(28.77)	214(15.51)	222(26.56)	234(38.12)	224(29.91)
	C.I.	$(186,316)$	$(195,271)$	$(200,244)$	$(188,291)$	$(198,295)$	$(186,297)$
$SEMWiG$777Prior on $\varphi$ is chosen as U($\hat{c},1$) since is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	213(13.44)	211(11.28)	208(7.90)	200(0.76)	216(16.71)	199(2.05)
	C.I.	$(200,230)$	$(203,219)$	$(203,211)$	$(193,209)$	$(210,222)$	$(193,204)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	214(15.72)	215(16.23)	210(11.26)	203(5.42)	220(21.52)	200(4.35)
	C.I.	$(197,230)$	$(205,226)$	$(206,215)$	$(194,215)$	$(211,231)$	$(190,214)$
$\varphi =0.80$
$Le{e}^a$	$\hat{N}$(RMSE)	196(7.78)	194(8.64)	199(3.89)	191(10.92)	199(9.26)	203(11.26)
	C.I.	$(169,242)$	$(184,220)$	$(193,213)$	$(178,247)$	$(186,229)$	$(175,258)$
$SEMWi{G}^b$	$\hat{N}$(RMSE)	192(8.43)	195(4.70)	195(4.48)	192(9.34)	199(3.45)	185(14.62)
	C.I.	$(188,198)$	$(189,200)$	$(193,197)$	$(186,196)$	$(194,204)$	$(182,190)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	196(8.25)	201(5.32)	201(3.25)	193(9.32)	201(4.39)	190(11.80)
	C.I.	$(181,208)$	$(190,212)$	$(195,207)$	$(183,206)$	$(193,215)$	$(178,202)$

Table 5. Table 5: Summary results of the three methods (Lee, SEMWiG & N ^ t ​ b I superscript subscript ^ 𝑁 𝑡 𝑏 𝐼 \hat{N}_{tb}^{I} ) applied to the populations A1-A6 with N = 500 𝑁 500 N=500 and when ϕ < 1 italic-ϕ 1 \phi<1 is known.

Method		A1	A2	A3	A4	A5	A6
$\varphi =0.60$
$Lee$888Prior on $\varphi$ is chosen as U($0.2,1.4$) since is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	625(136.80)	608(115.19)	553(56.30)	530(34.31)	619(124.43)	529(40.26)
	C.I.	$(485,827)$	$(505,744)$	$(511,619)$	$(472,671)$	$(510,775)$	$(464,669)$
$SEMWiG$999Prior on $\varphi$ is chosen as U($\hat{c},1$) since is known (see Chatterjee and Mukherjee, 2016b)	$\hat{N}$(RMSE)	521(25.48)	517(18.23)	517(17.16)	517(19.69)	527(27.73)	513(15.87)
	C.I.	$(490,549)$	$(505,529)$	$(514,521)$	$(499,539)$	$(512,542)$	$(491,538)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	541(43.00)	534(35.24)	528(28.34)	514(16.47)	548(49.80)	505(10.61)
	C.I.	$(519,570)$	$(520,550)$	$(522,535)$	$(500,528)$	$(535,566)$	$(486,520)$
$\varphi =0.80$
$Le{e}^a$	$\hat{N}$(RMSE)	472(35.42)	490(13.65)	509(10.34)	481(20.08)	531(39.79)	529(38.53)
	C.I.	$(431,565)$	$(460,560)$	$(485,545)$	$(448,545)$	$(475,633)$	$(459,608)$
$SEMWi{G}^b$	$\hat{N}$(RMSE)	478(25.27)	492(11.74)	498(5.28)	484(19.34)	490(12.15)	480(22.96)
	C.I.	$(455,507)$	$(477,510)$	$(489,508)$	$(465,506)$	$(485,496)$	$(462,493)$
${\hat{N}}_{tb}^I$	$\hat{N}$(RMSE)	487(17.81)	499(8.30)	502(5.29)	484(17.55)	506(10.19)	488(13.58)
	C.I.	$(465,512)$	$(484,515)$	$(492,511)$	$(470,498)$	$(490,521)$	$(469,505)$

Table 6. Table 6: Summary results of the proposed integrated likelihood estimate along with independence estimate, Lee’s full Bayes estimate and SEMWiG empirical Bayes estimates for Children Injury Data .

Class
(Estimate $\hat{c}$)		${\hat{N}}_{ind}$	Lee	SEMWiG	${\hat{N}}_{tb}^I$
Cyclists101010Methods developed in Chatterjee and Mukherjee (2016c, pp. 5-6) suggest recapture aversion	$\hat{N}$(s.e.)	495 (69.68)	292 (68.83)	244 (19.99)	303 (20.04)
(0.254)	C.I.	$(359,632)$	$(212,448)$	$(226,273)$	$(272,350)$
Passengersa	$\hat{N}$(s.e.)	249 (24.05)	193 (21.85)	181 (17.76)	192 (9.05)
(0.40)	C.I.	$(202,296)$	(163, 248)	$(170,208)$	$(175,212)$
Pedestriansa	$\hat{N}$(s.e.)	1323 (31.90)	1213 (99.69)	1110 (31.45)	1159 (13.53)
(0.592)	C.I.	$(1260,1385)$	$(1050,1424)$	$(1071,1148)$	$(1135,1186)$

Table 7. Table 7: Summary results of the proposed integrated likelihood estimate along with independence estimate, Lee’s full Bayes estimate and empirical Bayes estimate - SEMWiG for Handloom Data .

Population
(Estimate $\hat{c}$)		${\hat{N}}_{ind}$	Lee	SEMWiG	${\hat{N}}_{tb}^I$
Ward 2111111Chatterjee and Mukherjee (2016c, pp.10) suggests recapture proneness for this population	$\hat{N}$(s.e.)	159 (4.50)	187 (15.23)	168 (1.52)	164 (6.46)
(0.675)	C.I.	$(150,167)$	$(160,220)$	$(166,169)$	$(154,178)$
Ward 16121212Chatterjee and Mukherjee (2016c, pp.10) suggests recapture aversion for this population	$\hat{N}$(s.e.)	270 (21.53)	230 (37.56)	211 (9.67)	213 (7.95)
(0.382)	C.I.	$(228,312)$	$(185,318)$	$(202,225)$	$(200,230)$

Equations

$${\hat{N}}_{Nour}=x_0+\frac{2x_{11}{x}_{10}{x}_{01}}{(x_{11}^2+x_{10}{x}_{01})}.$$

$$L(N,p_{1\cdot },p,c)$$

$$L_{tb}(N,p_{1\cdot },p,\varphi )$$

$$L^I(\theta )$$

$$L_U^I(N)={\int }_{\Psi }L(\theta ,\psi |\text{text@underlinex})\propto (N+1{)}^{-1}(N-x_{1\cdot }+1{)}^{-1},$$

$$\pi (\psi |N)$$

$$L_J^I(N)$$

$${\hat{\gamma }}_{\theta }=\hat{\gamma }+O(N^{-1/2})O(|\theta -\hat{\theta }|)$$

$${\overline{L}}^I(\theta )$$

$$E\{{\ell }_{\psi }(\theta ,\psi );\hat{\theta },\gamma \}$$

$${\hat{\psi }}_{\theta }$$

$${\overline{L}}_{tb}^I(N)$$

$$N\le \frac{x_0(r_1+s_1-1)}{r_2+r_1}-1$$

$${\overline{L}}_{M_{tb}}^I(N,c,p)$$

$${\ell }_{p_{1\cdot }}(N,p_{1\cdot },c,p)$$

$${\ell }_c(N,{\psi }_1,{\psi }_2)$$

$${\ell }_p(N,p_{1\cdot },c,p)$$

$$E({\ell }_{p_{1\cdot }}(N,p_{1\cdot },c,p):N_0,p_{1\cdot ;0},c_0,p_0)$$

$$E({\ell }_c(N,p_{1\cdot },c,p):N_0,p_{1\cdot ;0},c_0,p_0)$$

$$E({\ell }_p(N,p_{1\cdot },c,p):N_0,p_{1\cdot ;0},c_0,p_0)$$

$$p_{1\cdot }$$

$$c$$

$$p$$

$$L_{tb}(N,\gamma )$$

$$L_{tb}(N,\gamma |\text{text@underlinex})\pi (\gamma )$$

$${\overline{L}}_{tb}^I(N)$$

$$\frac{{\overline{L}}_{tb}^I(N+1)}{{\overline{L}}_{tb}^I(N)}$$

$$\iff s_2+s_{2}N-x_0$$

$$\iff (N+1)(r_1+s_1-s_2)$$

$$\iff N$$

$$N$$

$$\frac{N+1}{x_0}$$

$$N^2{r}_2-N(x_{1\cdot }{r}_2+x_{01}{s}_1-r_2-x_0)-x_{1\cdot }(x_0+r_2-s_1)$$

Full text

Keywords : Capture-recapture, Direction of behavioral dependence, Human population, Time-Behavioral Response Variation model, Unrelated nuisance parameter.

1 INTRODUCTION

Dual-record System (DRS) is a special type of data-structure obtained from a capture-recapture type experiment, especially designed for estimating a specified population size, say $N$, from two sampling occasions. Federal agencies are generally interested to know the actual size of a specified population. Often census operations fail to extract the true knowledge on the size of the population for various reasons. Therefore, either any other contemporary census counting data is used or another survey is conducted independently after the census operation in order to estimate the true size $N$. Application of this technique is common in various interdisciplinary platforms, such as epidemiological study, socio-economic census, study of episodic events, etc. To estimate $N$ from these two lists of information (each of them is supposed to be incomplete), idea from a popular technique, called capture-recapture, is borrowed. Wolter (1986) sketched different capture-recapture models for N from DRS based on the pioneering work of Otis et al. (1978) who introduced several likelihood models for different plausible situations. In practice for homogeneous group, model $M_{t}$ has received much attention from both the frequentist and Bayesian statisticians due to its simplicity. $M_{t}$ accounts for time(t) variation effect and assumes causal independence between the sources of information. In the capture-recapture literature on model $M_{t}$, several Likelihood (e.g. Bishop et al., 1975; Huggins, 1989), pseudo-likelihood (e.g. Bolfarine et al., 1992; Chatterjee and Mukherjee, 2016a) and Bayesian methods (e.g. George and Robert, 1992; Xu et al., 2014) are available. In this context, the popular Lincoln-Petersen estimator is identical with the estimator derived from conditional likelihood (Chatterjee and Mukherjee, 2016a). But the specific assumption of causal independence in $M_{t}$ may seriously mislead in most of the situations for human population, especially when capture probabilities also vary with behavioral response (Chandrasekar and Deming, 1949). When both the time (t) variation effect and behavior response (b) effect acts together, model $M_{tb}$ is appropriate. Moreover, this model can be treated as most general statistical form of capture-recapture model for homogeneous population. Gosky and Ghosh (2011) found the model $M_{tb}$ as the most robust model in estimating $N$ based on comparative simulation study in Bayesian paradigm over all the models proposed in Otis et al. (1978). The underlying behavior response effect (say, $\phi$) classifies a given population as recapture prone (when $\phi>1$) or recapture averse (when $\phi<1$). Usually, demographic study exhibits recapture prone type dependence. On the contrary, study on drug abused population usually reveals that underlying list-dependence is negative, i.e., drug abused population is recapture averse. However, Otis et al. (1978) addressed the non-identifiability problem related to this model. Lloyd*′s (1994) martingale approach, Chao et al.′s (2000) quasi-likelihood approach and Yang and Chao′*s (2005) univariate Markovian method successfully solve the nonidentifiability for number of sampling occasions ($T$) strictly more than two. Lee et. al. (2003) successfully developed a full Bayesian technique with little informative prior but their demonstration is in the spirit of large number of sampling occasions, which is seldom used for human population. Later, Wang et al. (2015) proposed a hierarchical Bayesian $M_{tb}$ model for multiple lists with the assumption that the odds of recapture bears a constant relationship to the odds of initial capture. Capture-recapture type experiment with $T=2$ is commonly exercized in most of the applications for human population and estimation of $N$ is found to be an challenging task in the presence of causal dependence between sampling occasions. Chatterjee and Mukherjee (2016b) discusses some issues related to full Bayes method with non-informative prior in this context. They also developed some empirical Bayes strategies in DRS considering the present problem in a missing data framework. Generally, in Bayesian paradigm, difficulty may arise as the resulting estimator for $N$ may be very sensitive to the choice of prior(s). However, Nour (1982) proposed an estimator in DRS with the assumption equivalent to recapture proneness but avoiding Bayesian technique. Literature (Nour, 1982; Chatterjee and Mukherjee, 2016b) reveals that if correct available directional knowledge on $\phi$ is applied, inference on $N$ could be improved certainly in both of Bayesian and non-Bayesian paradigms for model $M_{tb}$. However, an efficient classification strategy (recapture proneness or aversion) of the given population in terms of the behavioral nature is proposed in Chatterjee and Mukherjee (2016c). Appropriateness and the challenging identifiability problem of the present model under DRS motivate us to consider the problem of $N$ estimation in the present paper. Here, we proposed a novel integrated likelihood method as a suitable non-Bayesian strategy to meet our goal particularly when the underlying population is correctly known as recapture prone or averse.

All the model parameters in $M_{tb}$ except the interest parameter $N$ are regarded as nuisance parameters (say, $\psi\in\Psi$). In these contexts, some useful likelihood-based inference procedures through the construction of pseudo-likelihood functions by eliminating the nuisance parameters are discussed in Severini (2000). This elimination of nuisance parameter may be done by maximization (profile likelihood) or conditioning (conditional likelihood) or integrating the likelihood function over $\Psi$ with respect to some weight function (integrated likelihood). Integrated likelihood has an advantage that it always exists unlike other pseudo-likelihood methods. Salasar et al. (2014) analysed integrated likelihood approach with uniform and Jeffrey’s prior for eliminating nuisance parameters in $M_{t}$. In integrated likelihood method the main challenge is to choose a suitable prior weight function on the nuisance parameters. Severini (2007) presents a novel approach for selecting a weight function so that the resulting integrated likelihood is useful for non-Bayesian inference and also posses some nice statistical properties. Recently, Chatterjee and Mukherjee (2016a) has developed an integrated likelihood for $M_{t}$, with the help from Severini (2007). In this article, we extend the work of Chatterjee and Mukherjee (2016a) for our present complex model $M_{tb}$-DRS. But the main challenge here is to successfully overcome the current model identifiability problem by suitably choosing informative priors so that some desirable properties hold. In summary, this article is framed to provide an alternative or supplement to the few existing methods in the literature of traditional homogeneous two-sample capture-recapture data (i.e. DRS) when the two lists are thought to be behaviorally dependent.

In the next section, we discuss data structure for DRS and relevant model $M_{tb}$. In section 3, at first we discuss the integrated likelihood method using weight function as uniform and Jeffrey’s densities. Later we propose a novel integrated likelihood method through the construction of unrelated nuisance parameter and informative priors. Evaluation of our proposed method by comparing with some other available Bayesian methods is carried out by an extensive simulation study in section 4. Thereafter, we illustrate our method by applying them to two real life data sets. Finally in section 5, we summarize our findings and provide some comments about the usefulness of the proposed integrated likelihood.

2 ANALYSES ON DUAL-RECORD SYSTEM: PRELIMINARIES

2.1 Dual-record data structure

The idea of Dual collection came from very popular capture-recapture analysis in wildlife management to estimate population size. Let us consider a human population U of size N is to be estimated. Any attempt to enlist all the individuals in $U$ is believed to be incomplete as it fails to capture all individuals in that population. In this paper we consider two common basic assumptions that (1) population is closed within the time of two sources gathering information, (2) individuals are homogeneous with respect to capture probabilities. To estimate the true N, minimum two sources of information covering that population is needed. When information is collected by two sources and classify all the captured individuals in U according to a multinomial fashion (see Table 1), then it is known as Dual-record System or Dual Collection. This type of classification is obtained by matching the individuals captured by the first (list 1) and second sources (list 2). The total number of distinct captured individuals by the two lists is $x_{0}$ (say), then $x_{0}=x_{1\cdot}+x_{\cdot 1}-x_{11}$. Clearly, the number of missed individuals $x_{00}$ by both systems is unknown and that makes the total population size N($=x_{\cdot\cdot}$) unknown. Expected Proportions or probabilities for each cell are also given in Table 1 and these notations will be followed throughout this paper.

Combining all the information, estimate of N could be obtained assuming different conditions on the individual’s capture probabilities leading to different models. A very common practice, across all fields of applications, is to assume casual independence between two lists’ probabilities. Hence, the conditional likelihood estimate from the resulting multinomial model (denoted as $M_{t}$) is $\hat{N}_{ind}=[x_{1\cdot}x_{\cdot 1}/x_{11}]$, which is popularly known as Lincoln-Petersen estimator (Otis et al, 1978) or dual system estimator (DSE); see Wolter (1986) and Chatterjee and Mukherjee (2016a) for details. But this model is highly criticized by several statisticians and practitioners due to the failure of its independence assumption in real life applications. In demographic study, violation of causal independence often is commonly observed due to positive correlation between two sources (or lists) of counts (see Chandrasekar and Deming, 1949; Nour, 1982;). Assuming such positive dependency, Nour(1982) proposed an estimate of $N$ as

[TABLE]

Negative dependence is observed in case of epidemiological surveillance of rare or critical disease, like HIV, drug abusing, etc.

2.2 Model $M_{tb}$

Causal independence assumption is criticised in surveys and censuses of human populations. The concern is that an individual’s probability of capture in List 1 may be change in response to capture in the second list. An individual who is captured in first attempt may have more (or less) chance to be included in the List 2 than the individual who has not been captured in first attempt. The change may occur due to different causes (see Wolter, 1986). This change is grossly known as behavioral response variation. When this chance is more, the corresponding individuals are treated as recapture prone, when chance is less, individuals are treated to be recapture averse. When this feature is combined with the time variation effect, one would get a complex model denoted as $M_{tb}$. To model this behavioral response variation generally, let us consider the following notations:

Prob(An individual present in List 1) = $p_{1\cdot}$,

Prob(An individual present in List 2 $|$ not present in List 1) = $p_{01}/(1-p_{1\cdot})$ = $p$ and

Prob(An individual is captured in List 2 $|$ captured in List 1) = $p_{11}/p_{1\cdot}$ = $c$.

Therefore, the likelihood function for model $M_{tb}$ in DRS is

[TABLE]

for $N>x_{0}$, $0<p_{1\cdot},p,c<1$, consists lesser number of sufficient statistics ($x_{11},x_{01},x_{10}$) than the parameters ($N,p_{1\cdot},p,c$) (see Otis et al., 1978). One can consider a popular assumption that recapture probability at second sample, c, is equal to a constant multiple of the conditional probability $p$, hence, $c=\phi p$. Chao et al. (2000) adopted this relation from Lloyd (1994) to get rid of the problem. Then reparameterized version of likelihood (1) becomes

[TABLE]

and $\phi$ is termed as the behavioral response effect characterizing the behavioral dependency of an individual belongs to the population at the time of second capture. In Equation (2), $\phi$ is orthogonal to $N$. It can be noticed that, in DRS, the dimension cannot be reduced by any reparametrization and therefore, identifiability problem persists in Equation (2) also. $\phi$ and $p$ are not estimable separately but their product $c$ is rather estimable. However, this second form of parameterization in Equation (2) may be of interest in lieu of Equation (1) (see Chao et al., 2000) as Equation (2) is characterized by the parameter $\phi$, which has a clear implication to define the nature of underlying behavioral dependence among two sources, i.e. whether the given population is recapture prone or averse. Replacement of $p$ with $c/\phi$ in Equation (1) is another version of parametrization and in this form, $\phi$ is not orthogonal to $N$.

3 INTEGRATED LIKELIHOOD METHOD

Let us consider a statistical model with likelihood function $L(\lambda|\underline{\textbf{x}})$ with $\lambda=(\theta,\psi)$, where $\theta(\in\Theta)$ is parameter of interest and $\psi(\in\Psi)$ represents nuisance parameter. Both the $\theta$ and $\psi$ may be vector valued. Presence of more nuisance parameters in the model affects the comparative inferential study based on the likelihood (see Severini, 2000). Now our aim is to find a function that can summarize the set of likelihoods $\mathcal{L}^{*}=\{L(\theta,\psi|\underline{\textbf{x}}):\psi\in\Psi\}$ over $\Psi$. That summerized function $L^{*}(\theta)$ of $\theta$, is some extent used as if the inference frame has $\theta$ as the full parameter and therefore has likelihood function, $L^{*}(\theta)$. We refer such functions $L^{*}(\theta)$ here as pseudo likelihood function of $\theta$. Construction of such pseudo-likelihood is performed by elimination nuisance parameter which can be handled by integrated likelihood method. In this approach nuisance parameter is eliminated through integration or rather it can be said that the set of likelihoods $\mathcal{L}^{*}=\{L(\theta,\psi|\underline{\textbf{x}}):\psi\in\Psi\}$ is summarized over $\Psi$ by a weighted average with respect to a chosen function on $\Psi$, say, $\pi(\psi|\theta)$ defined on $\Psi$. Hence, integrated likelihood function with respect to the weight $\pi(\psi|\theta)$ is

[TABLE]

see Severini (2000) for detailed discussion. One advantage of integrated likelihood over pseudo-likelihoods (conditional, marginal) is that it is always possible to construct unlike conditional or marginal likelihood. However, it is not necessary to choose $\pi(\psi|\theta)$ as a proper density function in this context. But one drawback includes the plausible subjectiveness due to $\pi(\psi|\theta)$. The basic aim always remains to choose a suitable $\pi(\psi|\theta)$ such that $L^{I}(\theta)$ could be efficiently useful for non-Bayesian likelihood inference. The mle of $L^{I}(\theta)$ is to be treated as the resulting estimate of $N$ from this method.

In capture-recapture context, for fixed $\theta$, Jeffrey’s and uniform priors are the two most popular non-informative prior densities on $\psi$ (Salasar, 2014). When $\theta=N$, $\psi=(p_{1\cdot},c,p)$ and uniform prior $\pi(\psi|N)\propto 1$ is chosen for $\psi$, then from Equation (3), integrate likelihood from likelihood (1) becomes

[TABLE]

for $N>x_{0}$. Clearly, $L_{U}^{I}(N)$ is strictly decreasing over its domain. Hence, $L_{U}^{I}(N)$ fails. Again, if we consider $\pi(\psi|N)$ as Jeffrey’s prior, then $\pi(\psi|\theta)\propto\sqrt{|\mathcal{I}_{N}(\psi)|}$, where $\mathcal{I}_{N}(\psi)$ is $3\times 3$ Fisher’s information matrix. Therefore,

[TABLE]

Hence, using the weight in Equation (4) in Equation (3), one would have the integrate likelihood corresponds to likelihood (1) as

[TABLE]

which is of $O(N^{-1})$ and hence, strictly decreasing over its domain $N>x_{0}$. Note that both of the above integrated likelihoods based on non-informative priors fails, because these priors could not add subjective information for $\psi$ in the likelihood (1) so that it became well-behaved and produce reasonable estimate for $N$.

3.1 Proposed Integrated Likelihood approach

In order to construct an integrated likelihood function to be useful for the present likelihood model, suitable informative prior should be used. There are two alternative parameterizations Equations (1) and (2) of the model $M_{tb}$ in section 2.2 and we develope our theoretical findings for both parameterizations on different $\psi$. To begin, let us start by fetching the idea of strongly unrelated parameters defined in Severini (2007).

Definition 1**.**

Two parameters $\gamma$ and $\theta$ are said to be strongly unrelated if

[TABLE]

holds.

Henceforth, we consider $\gamma$ as a nuisance parameter with the same dimension as $\psi$ and $\gamma$ is strongly unrelated to $\theta$, the interest parameter here. Therefore, we follow the prior choosing mechanism developed in Severini (2007). This mechanism helps us to find a strongly unrelated nuisance parameter $\gamma$, in terms of $\psi$ and $\theta$ for model $M_{tb}$, in such a way so that $\gamma$ and $\theta$ would become independent under $\pi(\psi|\theta)$, i.e., $\pi(\psi|\theta)=\pi(h(\gamma,\theta)|\theta)=\pi(\gamma|\theta)=\pi(\gamma)$. Therefore, the task is to find such parameter $\gamma$ and then choose a prior density $\pi(\gamma)$ for $\gamma$ that does not depend on $\theta$. Hence, the integrated likelihood function for $\theta$ with respect to $\pi(\gamma)$ is

[TABLE]

Construction of such nuisance parameter $\gamma$ is discussed below.

We consider the following equation from Severini(2007):

[TABLE]

from which one can solve $\gamma$ as $\gamma(\theta,\psi;\hat{\theta})$. Severini (2007) proved that $\hat{\gamma}=\hat{\psi}$ and $\gamma$ is strongly unrelated to $\theta$ i.e. $\hat{\gamma}_{\theta}=\hat{\gamma}+O(n^{-1/2})O(|\theta-\hat{\theta}|)$. Then solution $\gamma(\theta,\psi;\hat{\theta})$ is called zero-score-expectation parameter. Now, one can choose any suitable prior $\pi(\gamma)$ for $\gamma$ as $\overline{L}^{I}(\theta)$ in Equation (5) does not heavily depend on the chosen prior whereas for orthogonal parameters, proposed integrated likelihood (5) may be dependent to the choice of prior.

One can find $\gamma$ in different way. The aim is to find a function $\gamma(\theta,\psi)$ such that $\hat{\gamma}_{\theta}=\hat{\gamma}+O(n^{-1/2})O(|\theta-\hat{\theta}|)$. Hence, one can find such a parameter $\gamma=g(\theta,\psi)$, which implies $\psi=h(\theta,\gamma)$ for some h, if exists so that $\hat{\psi}_{\theta}=h(\theta,\hat{\gamma}_{\theta})=h(\theta,\hat{\gamma})+O(n^{-1/2})$. So, for any value of $\theta$, $\hat{\psi}_{\theta}$ depends on the data only through $\hat{\gamma}$. In many situations, this $\hat{\gamma}$ does not exist. So, we consider $\gamma$ as a function of $\hat{\theta}$ in addition to ($\theta,\psi$). It can be written as $\gamma=g(\theta,\psi;\hat{\theta})$ which implies $\psi=h(\theta,\gamma;\hat{\theta})$. As $\hat{\psi}_{\theta}=h(\theta,\hat{\gamma}_{\theta};\hat{\theta})$, then we must have

[TABLE]

when $\hat{\gamma}_{\theta}=\hat{\gamma}+O(n^{-1/2})O(|\theta-\hat{\theta}|)$. Hence, one wish to find such a function $h(\cdot)$ so that Equation (7) holds.

Let us address some desirable properties of $\overline{L}^{I}(\theta)$ in the context of model $M_{tb}$. By construction, unrelated parameter $\gamma$ is less related to $\theta$ than any nuisance parameter. Hence, prior sensitiveness to $\theta$ can be reduced by the use of $\pi(\gamma)$ in lieu of $\pi(\psi|\theta)$. Moreover, score unbiasedness and information unbiasedness criteria are incorporated in the present construction and the resulting likelihood estimate is invariant with respect to the different parameterizations. For detailed discussion of these properties, we refer Severini (2007).

Now we discuss the application of the proposed idea of integrated likelihood method to model $M_{tb}$ in the context of DRS. We construct the relevant unrelated nuisance parameter $\gamma$, then choose $\pi(\gamma)$ satisfying posterior unbiasedness condition. Firstly, we present the consequent theorems, results and properties for the parametrization in likelihood (1). Further, we also analyse the parametrization in likelihood (2) following similar course of theorems and results as for parametrization in likelihood (1).

Implementation for the likelihood (1)

Let us consider the parameter of interest $\theta=N$. Following theorem finds strongly unrelated parameter corresponding to nuisance parameter $\psi=(p_{1\cdot},c,p)$. Proof is given in the Appendix.

Theorem 1**.**

Theorem 3.1. For parametrization in (1), consider $\theta=N$ and $\psi=(p_{1\cdot},c,p)$. Then using Equation (6), strongly unrelated parameter is $\gamma$=($\gamma_{1},\gamma_{2},\gamma_{3}$), where $\gamma_{1}=(N/\hat{N}_{ind})p_{1\cdot}$, $\gamma_{2}=c$ and $\gamma_{3}=\frac{p(1-p_{1\cdot})}{(\hat{N}_{ind}/N)-p_{1\cdot}}$.

Since current model suffers from non-identifiability, so non-informative priors for resultant unrelated parameters $\gamma$ (in Theorem 1) would not work satisfactorily, as in the case of ordinary integrated likelihoods under uniform and Jeffrey’s priors. Thus, we consider some informative prior for $\gamma$ subject to the condition that hyper-parameters satisfy some relations that lead to a integrated likelihood.

Result 1**.**

In connection with Theorem 1, if the prior $\pi(\gamma)$ is of the form $\pi(\gamma)=\pi(\gamma_{1})\pi(\gamma_{2})\pi(\gamma_{3})$ and $\pi(\gamma_{1})=GB1(b_{1}=\frac{N}{\hat{N}_{ind}},r_{1},s_{1})$, $\pi(\gamma_{2})=Unif(0,1)$ and $\pi(\gamma_{3})=GB1(b_{2}=\frac{1-p_{1\cdot}}{(\hat{N}_{ind}/N)-p_{1\cdot}},r_{2},s_{2})$, for any positive real numbers $r_{2}$, $s_{2}$, $r_{1}$, $s_{1}$ satisfying $r_{2}+s_{2}=s_{1}$, where GB1() stands for Generalized Beta distribution of Type 1, then using Equation (5) integrated likelihood of $N$ for model $M_{tb}$ becomes

[TABLE]

for $N\geq x_{0}$, where $\Gamma(a)$ denotes the Gammaa function of $a$ equivalent to $\int_{0}^{\infty}e^{-x}x^{a-1}$. Thus, for the given values of $r_{1}$, $s_{1}$ and $s_{2}$, $\overline{L}_{tb}^{I}(N)$, in Equation (8), is non-decreasing in $N$ for

[TABLE]

and $\overline{L}_{tb}^{I}(N)$ converges to 0 as $N\rightarrow\infty$.

Proof. Given in the Appendix.

The rationale behind the assumption $r_{2}+s_{2}=s_{1}$ primarily is to make the resulting integrate likelihood well-behaved, and secondarily, to reduce dimension of the hyper-parameters in $\pi(\gamma)$. Now, on the basis of the relationship between $\gamma$ and $\psi$, mentioned in Theorem 1, we suggest the values of hyper-parameters $r_{2}$, $s_{2}$ and $r_{1}$ using posterior unbiasedness of $\gamma$ (after considering $\gamma$ as function of $N$). We suggest $r_{2}=bx_{12}$ and $s_{2}=b(N^{*}-x_{0})$ for some real positive $b$, where $N^{*}$ is some other working estimate of $N$. So, $s_{1}=b(N^{*}-x_{1\cdot})$. Now, for fixed $N$, $\gamma_{1}=(N/\hat{N}_{ind})p_{1\cdot}=(N/\hat{N}_{ind})(x_{1\cdot}/N)=(x_{1\cdot}/\hat{N}_{ind})$. Hence, from the posterior unbiasedness condition regarding unrelated parameter $\gamma_{1}$, we have $(x_{1\cdot}/\hat{N}_{ind})\simeq E_{\pi}(\gamma_{1})=(N/\hat{N}_{ind})\{(r_{1}+x_{1\cdot})/(r_{1}+s_{1}+N)\}$ and this implies $r_{1}=x_{1\cdot}s_{1}/(N-x_{1\cdot})=r_{1}(N,s_{1})$, say, after some algebraic manipulation. Thus, $r_{1}$ depends on $N,s_{1}$ and therefore, the right hand side of the condition in Equation (9) also becomes dependent on $N$. So, $N\leq\frac{x_{0}\{r_{1}(N)+s_{1}-1\}}{\{r_{2}+r_{1}(N)\}}-1\Leftrightarrow\frac{x_{0}\{r_{1}(N)+s_{1}-1\}}{\{r_{2}+r_{1}(N)\}}-N\geq 1$ and $\frac{x_{0}\{r_{1}(N)+s_{1}-1\}}{\{r_{2}+r_{1}(N)\}}-N$ is non-increasing in $N$ for fixed $s_{1}$ and $r_{2}$. Hence, the result stated in the following theorem, obtained from a very easy algebraic simplification, discusses the possibility of existence of the corresponding maximum likelihood estimate. Proof is sketched in the Appendix.

Theorem 2**.**

For $r_{1}=x_{1\cdot}s_{1}/(N-x_{1\cdot})=r_{1}(N,s_{1})$, $\exists$ a real number, say $N_{0}=N_{0}(\underline{\textbf{x}},s_{1},r_{2})$, for which inequality in Equation (9) is equivalent to $N\leq N_{0}<\infty$.

Hence from Theorem 2, it can be said that $\overline{L}_{tb}^{I}(N)$ is increasing in $N$ for $N\leq N_{0}$ and hence, estimate based on our proposed integrated lkelihood method is $\hat{N}_{tb}^{I}=[N_{0}]+1$, if $N_{0}$ is not an integer and $\hat{N}_{tb}^{I}=[N_{0}]$ and $[N_{0}]+1$ if $N_{0}$ is an integer. The expression for $N_{0}$ turns out to be mathematically intractable and thus obtaining an explicit solution is not possible. We explore several properties of the resulting estimator, as pointed out in Severini(2007), through computation in section 4.1. To implement the above prior specification, we suggest $N^{*}=(\hat{N}_{Nour}+\hat{N}_{ind})/2$ and $b=(1+(\hat{N}_{ind}-((x_{0}+\hat{N}_{ind})/2)-1)^{-1})/2$, when we don’t know anything about the plausible direction of $\phi$. On the other hand, if we know that $\phi>1$ (recapture proneness), then $N^{*}=\hat{N}_{Nour}$ and $b=1$. If $\phi<1$ (recapture aversion), suggested $N^{*}=\hat{N}_{ind}$ and $b=(\hat{N}_{ind}-((x_{0}+\hat{N}_{ind})/2)-1)^{-1}$.

Remark 1**.**

If anyone is interested to consider $\theta=(N,c,p)$, where $c=\phi p$, relevant for likelihood (1) then strongly unrelated parameters can be obtained from the relevant log-likelihood functions as before. Now, for $\psi=p_{1\cdot}$, if the prior is taken on the associated unrelated parameter $\gamma$ as $\pi(\gamma)=GB1(b=\frac{N}{\hat{N}_{ind}},r,s)$, for any positive real numbers $r$ and $s$; then integrated likelihood for $N(\geq x_{0})$ reduces to

[TABLE]

and it fails to produce the mle of $\theta=(N,c,p)$. Perhaps the failure is due to the scarcity of enough information to make inference about the two parameters $N$ and $\phi$, which are actually orthogonal to each other. For details about parameter orthogonality, see Cox and Reid (1987).

4 NUMERICAL ILLUSTRATION

4.1 Simulation study

In this section we conduct a simulation study to evaluate the performance of our proposed approach and understand its efficiency in applying the method using the available directional knowledge on $\phi$. This study is designed as follows. Let us simulate hypothetical populations corresponding to six pairs of capture probabilities ($p_{1\cdot},p_{\cdot 1}$)={(0.50, 0.65), (0.60, 0.70), (0.80, 0.70), (0.70, 0.55), (0.55, 0.75), (0.70, 0.50)} for each case of recapture prone (represented by $\phi=1.25,1.50$) and recapture averse (represented by $\phi=0.60,0.80$) situations. We denote six populations corresponding to the six pairs of ($p_{1\cdot},p_{\cdot 1}$) as P1, P2, P3, P4, P5 and P6 respectively for recapture prone situations and associated results are presented in Tables 2 and 3 for true population size $N=200$ and $500$ respectively. Results of the other six populations comprising the same six pairs of ($p_{1\cdot},p_{\cdot 1}$) reflecting recapture averse situations, namely A1-A6 are shown in Tables 4 and 5 for the same two $N$ values respectively. 1000 data sets on ($x_{1\cdot},x_{\cdot 1},x_{11}$) are generated from each of the 12 populations. Our proposed integrated likelihood estimate have been obtained for each data set. Finally, estimate of $\hat{N}_{tb}^{I}$ is obtained by averaging over 1000 posterior means. Based on those 1000 estimates, the sample RMSE (Root Mean Square Error) and $95\%$ confidence interval (C.I.) are computed. In addition to that, we also compute similar statistics for Lee’s and SEMWiG method, proposed respectively in Lee et al. (2003) and Chatterjee and Mukherjee (2016b), to compare the performance of our newly proposed non-Bayesian method. However, Lee et al. (2003) illustrated their approach in the context of animal capture-recapture experiment with a large number of sampling occasions. To compute the estimates using Lee and SEMWiG methods in this article, we use same priors as considered in Chatterjee and Mukherjee (2016b) when information is available on the directional nature of $\phi$. Details of their computation strategy, particularly for DRS, can be found in Chatterjee and Mukherjee (2016b). All of these comparison results are summarized in Tables 2-5.

From the Tables 2-5 it can be clearly noticed that our proposed estimates are more efficient than Lee’s in all the situations, in terms of accuracy, RMSE and shorter length of confidence interval. When $\phi$ is far below $1$ (i.e. for $\phi=0.60$), SEMWiG produce slightly better results for small population. This discrepancy increases when $N$ is larger. In all other situations, performance of $\hat{N}_{tb}^{I}$ is better than SEMWiG in all the 12 populations except P1 for both the values of $\phi$. Precisely, except P1 and the populations with $\phi=0.60$, we found that $\hat{N}_{tb}^{I}$* $<$ SEMWiG $<$ Lee*, in terms of RMSE and SEMWiG $<$ $\hat{N}_{tb}^{I}$ $<$ Lee in terms of length of the associated interval estimates of $N$. Further, invariance property of the estimates obtained through our proposed integrated likelihood is verified as the estimates for true $N=500$ are 2.5 times higher than that for true $N=200$. Lastly, another important feature of $\hat{N}_{tb}^{I}$ is that being a non-Bayesian pseudo-likelihood based inferential strategy, it does not incur serious computational effort as for existing Bayesian strategies - Lee’s and SEMWiG.

4.2 Real data Example I: Children Injury Data

In Epidemiological study, use of capture-recapture experiment is very popular but more than two lists are hardly ever found. The simple estimate $\hat{N}_{ind}$ assuming list-independence is widely employed in this domain, even sometimes without judging its relevancy. Here we consider a work by Jarvis et al. (2000), in which authors illustrate the serious drawbacks in the use of this estimator specifically for injury related data. The problem was to get the count of children under 15 years of age from addresses in Northumbria who were seriously injured in local Motor Vehicles Accidents (MVA) between 1 April, 1990 and 31 March, 1995. One source was Stats19 data covering all road traffic accidents in Northumbria causing injuries to children that had been reported to the police and another was the Hospital Episode data (HES) covering admissions of children. The associated DRS data are presented in Jarvis et al. (2000, Table 4, pp. 48) for three different classes - Cyclists, Passengers and Pedestrians. Jarvis et al. argued that children injured in MVAs as pedestrians or cyclists rarely enter insurance claims for which they have to inform police for case diary. Sometimes the police, in establishing whether an injury is serious, are recommended to contact the hospital to find out whether the child is admitted or not. The common estimates under independence ($\hat{N}_{ind}$) for these three classes are shown in third column of Table 6. It is noted that $\hat{N}_{ind}$’s are more than twice the total number of cases actually observed ($x_{0}$). Also, value of the estimate $\hat{c}$ for these three classes are 0.25, 0.40 and 0.59 respectively, which are substantially small. All these direct to the possibility of list dependency (indicating recapture aversion, due to very small amount of recapture) and this motivate us to include this example in our illustration. These three classes have a common feature that $x_{1\cdot}<x_{\cdot 1}$ (the next example of Handloom Data has the opposite feature). We present summary of results in Table 6 for our proposed integrated likelihood method along with Lee et al’s Bayes and Chatterjee and Mukherjee’s empirical Bayes SEMWiG estimates for comparison.

For the data on all the classes in Table 6, all the competing estimators (Lee, SEMWiG & $\hat{N}_{tb}^{I}$) agree with the negative departure from independence, i.e. recapture aversion and the resulting dual system estimates $\hat{N}_{ind}$ produce large estimates which seriously overestimates $N$. Lee’s estimates posses larger variation than all other estimates and hence its credible intervals are too wide. Estimate $\hat{N}_{tb}^{I}$ has better efficiency than all other estimates. The proposed $\hat{N}_{tb}^{I}$ with the suggested directional knowledge (i.e. $\phi<1$) produces estimate between Lee and SEMWiG. Further, in most of the all cases, it has smaller variance and tighter confidence bounds, as expected.

4.3 Real data Example II: Handloom Data

Let us consider a new data from a survey aimed to estimate the undercount in the census of handloom workers residing at Gangarampur in South Dinajpur district of state West Bengal in India. The survey was post enumeration type (i.e. PES) and conducted in November 2013, three months after the census operation (see SOSU (2014)). This data sets is also used in Chatterjee and Mukherjee (2016c) for real data illustration of their proposed behavioral dependence classification methods. Their classification strategies reveal that Ward 2 and Ward 16 has the nature of recapture proneness and aversion respectively. For details on the associated data sets in DRS format and possible threat of list-dependence through behavioral response variation, we refer to Chatterjee and Mukherjee (2016c). Being quite assured about the homogeneity within wards from the experts of Textile Directorate, we apply the model $M_{tb}$ for these data and compute the estimates following our proposed integrated likelihood method (in section 3.1). Here also we also compute Lee’s Bayes estimate and the empirical Bayes estimate - SEMWiG for comparison in dependence situation. We also report the summary results if list-independence is assumed in order to measure the extent of deviation of other dependent estimates from independence.

Table 7 says that the proposed $\hat{N}_{tb}^{I}$ finds around 164 handloom workers assuming $\phi>1$. For the other sampled ward, when we incorporate the recapture aversion suggestion in our proposed method, it implies that approximately 213 workers are residing, which is very close to the SEMWiG estimate. Under the consideration of recapture proneness, efficiency of proposed estimator is better than Lee’s but little smaller than SEMWiG. When recapture aversion is assumed, $\hat{N}_{tb}^{I}$ is most efficient.

5 CONCLUSIONS

List-independence assumption does not hold satisfactorily in many instances. Various data from epidemiological studies and undercount or overcoount in demographic data motivate us to use a suitable model by avoiding the assumption of list-independence. As far as homogeneous human population size estimation is concerned, two-sample capture-recapture experiment is appropriate along with $M_{tb}$ modelling. Here, we consider integrated likelihood method as a non-Bayesian strategy which has a potential to overcome the non-identifiability in $M_{tb}$-DRS model. We have shown that general integrated likelihoods using common non-informative priors fail to produce estimate for population size $N$. To overcome this shortcoming, here we proposed an integrated likelihood based on suitable prior for unrelated nuisance parameter with the help from a novel idea by Severini (2007). This pseudo-likelihood mechanism produces efficient estimates satisfying several properties including invariance, less prior sensitiveness, etc. To get rid of the aforesaid nonidentifiability problem, choice of suitable informative prior is always a statistical challenge. However in many instances (e.g. epidemiology or demography or economic), direction of the underlying behavioral dependency (i.e., whether the given population is likely to be recapture prone or averse) can be anticipated correctly. Therefore, choice of hyperparameters in priors are suggested depending upon the availability of such directional knowledge on $\phi$. Indeed, as per our knowledge, this article presents the first efficient non-Bayesian strategy for the complex $M_{tb}$-DRS model. Our simulation study supports the fact that this newly developed pseudo-likelihood method is either more efficient in some regular situations than the existing methods or it is comparable to them in other situations. Hence, this integrated likelihood method can be treated as an efficient and application worthy alternative estimation mechanism, which of course incurs very less computation burden than other existing methods for estimating $N$ in this context.

Appendix

Proof.

Proof of Theorem 1. The log-likelihood of the model $M_{tb}$ with parameterization in (1) is

$\ell(\theta,\psi)=\ell(N,p_{1\cdot},c,p)=\Sigma_{i=0}^{x_{0}-1}ln(N-i)+x_{1\cdot}ln(p_{1\cdot})+(N-x_{1\cdot})ln(1-x_{1\cdot})+x_{01}ln(p)+(N-x_{0})ln(1-p)+x_{11}ln(c)+x_{10}ln(1-c)$ and, therefore

[TABLE]

Now, if we take expectations on the partial derivatives of log-likelihood functions over the distribution fixing ($\theta,\psi$) at ($\theta_{0},\psi_{0}$), we would have

[TABLE]

Using Equation (6) we obtain the following result for strongly unrelated nuisance parameters ($\gamma_{1},\gamma_{2},\gamma_{3}$) from Equations (10)-(12) respectively as

$\gamma_{1}=(N/\hat{N}_{ind})p_{1\cdot},\hskip 14.45377pt\gamma_{2}=c,\hskip 14.45377pt\gamma_{3}=\frac{p(1-p_{1\cdot})}{(\hat{N}_{ind}/N)-p_{1\cdot}}.$ ∎

Proof.

Proof of Result 1. Nuisance parameters $\psi=(p_{1\cdot},c,p)$ of the model $M_{tb}$ are expressed by the unrelated parameters $\gamma$, stated in Theorem 1, as follows:

[TABLE]

Now, replacing $\psi=(p_{1\cdot},c,p)$ in likelihood (1) by $\gamma$ using equations (13)-(15), we have rewritten the likelihood (1) as

[TABLE]

provided $\gamma_{1}<(N/\hat{N})$. In order to proceed for integrated likelihood, we start by choosing Beta distributions for original nuisance parameters $p_{1\cdot}$ and $p$ with parameters ($r_{1},s_{1}$) and ($r_{2},s_{2}$) respectively, for positive real numbers $r_{1}$, $s_{1}$, $r_{2}$, $s_{2}$. Therefore, newly determined unrelated nuisance parameters $\gamma_{1}$ and $\gamma_{3}$ would follow the distributions $\pi(\gamma_{1})\equiv GB1(b_{1}=\frac{N}{\hat{N}_{ind}},r_{1},s_{1})$ and $\pi(\gamma_{3})\equiv GB1(b_{2}=\frac{1-p_{1\cdot}}{(\hat{N}_{ind}/N)-p_{1\cdot}},r_{2},s_{2})$ respectively, where GB1() stands for Generalized Beta distribution of Type 1. We also consider the prior on $c=\gamma_{2}$ as $\pi(\gamma_{2})\equiv Unif(0,1)$. So, $\pi(\gamma)$ is of the form $\pi(\gamma)=\pi(\gamma_{1})\pi(\gamma_{2})\pi(\gamma_{3}|\gamma_{1})$.

Hence, following some algebraic simplification $L(\theta,\gamma|\textbf{\text@underline{x}})\pi(\gamma)$ in Equation (5) becomes

[TABLE]

with $\theta=N$. Now, if we integrate Equation (16) w.r.t. $\gamma$, then we have

[TABLE]

assuming $r_{2}+s_{2}=s_{1}$. Thus, we have Equation (8).

Now,

[TABLE]

Therefore,

[TABLE]

Hence the result in Equation (9) since $r_{2}+s_{2}=s_{1}$. $\Box$ ∎

Proof.

Proof of Theorem 2. For $r_{1}=x_{1\cdot}s_{1}/(N-x_{1\cdot})=r_{1}(N,s_{1})$, $\exists$ a real number, say $N_{0}=N_{0}(\underline{\textbf{x}},s_{1},r_{2})$, for which inequality in Equation (9) is equivalent to $N\leq N_{0}<\infty$.

[TABLE]

Therefore, the above inequality simplifies to

[TABLE]

Hence there exists two roots of the quadratic equation (17), say $N^{\prime}_{0}$ and $N_{0}$, both functions of ($\underline{\textbf{x}},s_{1},r_{2}$) such that $N^{\prime}_{0}\leq N\leq N_{0}.$

Now, the constant part of the quadratic function in the left side of the inequality (17) remains positive iff $b<\frac{x_{0}}{N^{*}-x_{0}}$, since $r_{2}-s_{1}=-s_{2}=-b(N^{*}-x_{0})$ (see section 3.1). Hence, the smaller root $N^{\prime}_{0}$ becomes negative iff $b<\frac{x_{0}}{N^{*}-x_{0}}$. Unless $\frac{x_{0}}{N^{*}-x_{0}}\leq 1$ or equivalently, $N^{*}\geq 2x_{0}$, $N^{\prime}_{0}$ will be negative for both of the choices of $b$ values mentioned in section 3.1. $\Box$ ∎