Nested Dirichlet Process for population size estimation from multi-list recapture data

The estimation of the size of a closed population from multi-list recapture data has been studied in many settings, for example estimation of census undercount [Chao and Tsay, 1998][Darroch et al., 1993], estimation of deaths in armed conflict [Ball et al., 2003][Manrique-Vallier et al., 2013][Manrique-Vallier et al., 2019], estimation of drug injectors [Overstall et al., 2014] and estimation of human trafficking victims [Heijden, 2016]. In general, each record in the multi-list recapture data has descriptive features, like time, location, gender, age, etc. To reduce uncertainty of list capture probabilities imposed by hierarchical structure such as location difference, it’s necessary to account for heterogeneity of response patterns in estimating the population size. One way to account for this heterogeneity is stratification. Stratification by every available location or time category may result in too many strata [Ball et al., 2003] [Manrique-Vallier et al., 2019]. Combining categories based on expert opinion is highly subjective. In this paper, we put heterogeneity caused by location and time into the model by building a non-parametric multi-layer latent class model based on the non-parametric one layer latent class model for multiple recapture data (LCMCR) proposed by Manrique-Vallier [2016]. In Manrique-Vallier (2016)’s paper, the latent layer models individual heterogeneity and individuals in the same latent class are independently captured by each data source. To reflect the hierarchical structure of the data, we add an additional layer on top of the individual layer to capture the top group (location-time) differences and to allow dependence among individuals in the same top latent layer.
Our work is motivated by the problem of estimating the number of casualties in the Syrian conflict. This conflict spanned many years and many regions of Syria, with times and periods of relative calm and heated conflict. We consider data from four reporting sources, with different reporting patterns across time, region, and type of death.
Many techniques estimate the population size by modeling list dependency. A class of generalized linear models, known as log-linear models [Bishop et al., 1975] assume the expected log of cell count is linearly related to a set of list interactions. Averaging over Bayesian decomposable graphical models, which represent graphical models of list dependency, is also a classical method to estimate the population size from multi-list recapture data [Madigan et al., 1995][Madigan and York, 1997]. Those methods treat all individuals the same which may not be proper in some cases. For example civilian and military deaths in the Syrian conflict data are captured differently by some data lists. Rasch models and extensions on them [Rasch, 1993] [Darroch et al., 1993] [Agresti, 1994] [Fienberg et al., 1999] incorporate individual heterogeneity into the log-linear model. A more flexible method, mixture models has also been used to capture individual heterogeneity [Manrique-Vallier and Fienberg, 2008] [Manrique-Vallier, 2016]. One strong assumption in the one layer latent class model is that individuals are independent, which might not be proper for data with hierarchical structure.
A popular alternative to the one layer latent class model for solving the individual dependence problem in nested data is multi-level latent class models [Vermunt, 2003] [Teh et al., 2006] [Rodrlguez et al., 2008]. However, multi-level latent class models haven’t been applied in population size estimation for multi-list recapture data. In over-time and across-location multi-list recapture data, we want the top layer to capture location-times that having similar recording patterns and the bottom layer to capture hidden classes of individuals within the top layer latent class. To realize this goal, both the hierarchical Dirichlet process (HDP) [Teh et al., 2006] and the nested Dirichlet process (NDP) [Rodrlguez et al., 2008] are great candidates. NDP allows both mixture components and weights to change within different top layer classes, but HDP only differs in weights. Due to the complicated and potentially highly distinct list dependencies among top layer latent classes, recording patterns might differ very much between one class containing locations with intense conflicts and one class containing locations with much less conflict. Therefore, we use the nested Dirichlet Process in this paper. NDP is usually applied in clustering nested data, like data with topic hierarchies [Blei et al., 2003][Blei et al., 2010][Fox et al., 2011]. In this paper, we extend it for multiple recapture data to identify more accurate hidden homogeneous classes among top level groups and among individuals within top level groups to better estimate population size from multi-list recapture data.
The article is organized as follows. In Section 2, we talk about the data and problem that motivates us for this paper. Then we review the one layer non-parametric latent class model for multiple recapture data (LCMCR), introduce our proposed approach nested latent class model for multiple recapture data (NLCMCR), and we apply MCMC inference for parameter estimation in Section 3. In Section 4, we do simulations to compare results from the one layer latent class model (LCMCR) and our nested latent class model (NLCMCR). In Section 5, we apply the NLCMCR in a sample of Syrian conflict data to estimate population size. In Section 6, we discuss this paper and make conclusions.

Human Rights Data Analysis Group (HRDAG) is a non-profit organization that applies rigorous science to the analysis of human rights violations around the world. One of its project is to estimate the total number of killings during the Syrian conflict based on multi-list recapture data. The Syrian conflict data we are using contains documented, identifiable victims who were killed during Syrian conflict from March 2011 to March 2016. Each death record has variables describing this person, which include the person’s name, death date, governornate (region in Syria), gender, age. Deaths were recorded by four data sources $(S=4)$ investigating deaths in the Syrian conflict, namely Syrian Center for Statistics and Research (SCSR), Damascus Center for Human Rights Studies (DCHRS), Syrian Network for Human Rights (SNHR) and Violations Documentation Center (VDC). Each record might be captured by more than one data source, thus the number of capture patterns is $2^{S}-1=15$ excluding the undocumented killings, with $S$ as number of data sources. Due to data confidential, in this paper we randomly generate a sample of $n=36226$ from all the documented killings. More details about the full documented victims have been discussed in Price et al. [2013a], Price et al. [2013b], and Price et al. [2014]. The number of killings recorded under each pattern in this sampled Syrian conflict data is summarized in Table 1. We can see that $n_{1000}=6039$ deaths are captured by VDC only, $n_{1010}=652$ are captured by VDC and DCHRS, not by SNHR and SCSR, and $n_{1111}=4252$ are captured by all four data sources. Estimating the undocumented killings is equivalent to estimating $n_{0000}$, and is the goal of our number inference.

In the Syrian conflict data, documented killings are from 14 governorates across the country. Since our sampled data are generated randomly from the full data set, recording patterns in our sampled data are similar as the full documented Syrian conflict data. From Figure 1, we can see that recording patterns within governorate change overtime. For example from 04/2011 to 12/2012 deaths captured by all four sources overtake records in other patterns in Rural Damascus. From 01/2013 to 08/2014, more deaths are captured by VDC, SNHR and SCSR together, but not by DCHRS. From 03/2015 to 12/2015, $n_{1101}$ is larger than others or more deaths were captured by VDC, DCHRS and SCSR, but not by SNHR. Some sources capture more killings than others in some governorates, for example, most killings were recorded by VDC in Tartus. Meanwhile, the documented number of killings recorded in different governorates differs much too. All those differences are not hard to explain if we think about the location of each governorate, when and where a small or a big event happened. With those findings, we believe that it’s not a good idea to combine all the death records simply over all time and governorates like what we did in Table 1 to estimate the total number of killings. Due to the long period and many governorates in this data set, it’s also a challenge to do proper stratification subjectively. Therefore, our nested model is important for detecting higher level (e.g. governorate-time) strata in this problem.

In this section, we generate multiple systems recapture data from a two layer latent class model, then we estimate the population size in three different ways:

1. 1.

   LCMCR: Latent class model for multiple recapture data.

2. 2.

   Multi-LCMCR: Fit LCMCR on each top layer latent class which is known in the simulated data, and then sum up population size estimations for these two sub-groups to get the overall population size.

3. 3.

   NLCMCR: nested latent class model for multiple recapture mode.

For the simulated data, we use $S=4$ data sources, $J=100$ (e.g. 100 location-times) top layer groups, $N=10000$, and number of individuals under each top layer group ($N_{j}$) ranges from 2 to 602 with a standard deviation of 116. Other parameters for the simulated data are listed in Table 2. We simulate data by assuming groups within each top layer have similar recording patterns. About $40\%$ of top groups belongs to the first latent class and $60\%$ in the second latent class. From Table 2, we can see that all four data sources have strong capture probability and they have many overlapping recordings when the top layer latent class is 1 ($z^{(2)}=1$) and individuals are more likely to be captured by the first data source only or captured by both the first and second data sources when they are in the top layer latent class of 2 ($z^{(2)}=2$). Therefore, we can see an obvious nested structure in the simulated data.

Figure 5 is boxplot of posterior estimation of population size under three different models: LCMCR, Multi-LCMCR and NLCMCR. The yellow line is the true population size. We can see that the estimation in LCMCR (red boxplot) is more biased. However, if we use LCMCR to estimate the population size separately by the top layer latent classes, the estimation is much better (green boxplot). This result makes sense because the data within each top layer latent class is from a one layer latent class model. If we stratify the data based on the true top layer latent class, we’ll get very good population size estimation using LCMCR for each strata, thus ideal estimation for the overall population size. The estimation from NLCMCR (blue boxplot) is very close to the result from the Multi-LCMCR method. This means that NLCMCR successfully detected the capture pattern differences among top level groups and dependency among individuals in the same top layer latent class. We can also see that NLCMCR gives much smaller uncertainty for the population size estimation than LCMCR.
Table 3 summarizes estimates and $95\%$ credible intervals for latent class proportions and list capture probabilities by using the one layer latent class model LCMCR and the nested model NLCMCR. We can see that both LCMCR and NLCMCR learn the number of latent classes as 4 which is the truth. Comparing parameters got from NLCMCR with the true values from which we simulation the data in Table 2, the estimates are almost unbiased with small uncertainty. To compare parameters got from LCMCR with the true values in Table 2, we first multiply the top layer proportion by the bottom layer proportion in Table 2 to get the true proportion for the corresponding latent class. We get proportions for the four latent classes $(0.4*0.8=0.32,0.4*(0.2)=0.08,0.6*0.6=0.36,0.6*0.4=0.24)$, which differ with the latent class proportions got from LCMCR $(0.53,0.26,0.17,0.03)$ clearly. Since the model LCMCR doesn’t group individuals into their true latent classes, the list capture probabilities within each latent class don’t reflect the true capture probabilities in Table 2 as well. Overall, both population size estimation and parameter estimation strongly suggest the importance of accounting for hierarchical structure in the multiple recapture data when hierarchy is present in the data.

From the descriptions in Section 2, we know that the Syrian conflict multiple systems recapture data has a hierarchical structure over time and governornate. From the simulation results we know that hierarchical structure is important to consider when estimating the population size. Therefore, we take governornate-time as the top group (or second layer) and individuals as the first layer. Then we apply the nested Latent class model for recapture model (NLCMCR) to estimate the total number of killings from the sampled Syrian conflict data. We also apply LCMCR and Bayesian model averaging over decomposable graphical models (BMD)[Madigan and York, 1997][Madigan et al., 1995] to compare the results, which are summarized in Table 4. From Table 4 we see that the population size estimation results from NLCMCR and LCMCR are similar, but NLCMCR has much smaller confidence interval and higher lower credible interval bound. Bayesian model averaging over decomposable graphical model (BMD) gives a much smaller estimation.
Figure 6 summarizes the clustering result of the top layer (governornate-time). We can see that most of the records from Tartus, As-Suwayda and Latakia are clustered into the sixth class ($z^{(2)}=6$). From Figure 1 we know that it is because they have similar recording patterns and most deaths recorded in those governorates were from VDC. From Table 5, we can also see that individuals in this group are mainly come from the first individual latent class ($z^{(1)}=1$) with VDC having a much higher capture probability than other data sources (VDC: 0.38, SNHR: 0.038, DCHRS: 0.005, SCSR: 0.011). The second plot in Figure 9 gives us some idea about how individuals in each goverornate-time were captured, for example after 06/2014 in Tartus all four data sources have very small capture rate. Governorate-times having red inversed triangles indicates some individuals were captured because of stronger capture probabilities of VDC (0.502) and SNHRS (0.61).
We can see that triangles ($z^{(2)}=2$) in Figure 6 are mainly from all governorates, excluding Tartus, As-Suwayda and Latakia, from around 03/2013 to 08/2014. Individuals in this group are much less likely to be captured by DCHRS, which is also reflected in Table 5. Missing individuals in this group are either very likely to be captured by VDC, SNHRS and SCSR together not by DHCRS (VDC: 0.906, SNHRS: 0.841, DCHRS: 0.031, SCSR: 0.882), or very likely to be captured by none of them (VDC: 0.33, SNHRS: 0.254, DCHRS: 0.012, SCSR: 0.16). In the second plot of Figure 7, we can see that red triangles are larger than black triangles in Aleppo, Rural Damasucs, Idlib and Darra which indicates more individual are captured by VDC, SNHRS and SCSR in those governorates.
Similarly, we can summarize that undocumented individuals in Aleppo, Rural Damascus, Damascus, Idlib, Deir ez-Zor, Homs, Hama and Daraa from 03/2011 to 12/2012 (black circled gov-time ($z^{(2)}=1$) in Figure 6) are more likely to be missed by DCHRS (0.457) and SCSR (0.386) or by the four data sources together. In those governorates from 03/2015 to 12/2015 ($z^{(2)}=3$), most undocumented individuals are either not likely to be captured by all four data sources or only not likely to be captured by SNHRS (0.073). Red plus signs in Aleppo, Rural Damasucs, Idlib and Darra in the first plot of Figure 8 dominate, which is similar as the situation in the second plot of Figure 7, but with individuals captured by VDC, DCHRS and SCSR. In those governorates in 10/2012, 01, 02, 04, 05/2013 and 12/2014 ($z^{(2)}=4$), documented individuals mainly come from individuals who are captured by VDC and SNHRS. Undocumented individuals in gov-times clustered into $z^{(2)}=5$ are either not likely to be captured by all four data sources or only not likely to be captured by SCSR (0.395).

In this paper, we find similar capture patterns in some governorates and months in the sampled Syrian conflict recapture data, which means that heterogeneity is clear in this dataset. In order to combine heterogeneity with modeling and allow information sharing across and within strata, we extend the one layer latent class model with Dirichlet Process prior to a multi-layer latent class model with Nested Dirichlet Process prior to estimate population size from multi-list capture data (NLCMCR). In clustering problems, NDP is preferred when the data has a hierarchical structure as it allows dependence for objects within the same top layer latent class. In our multi-list recapture setting, NLCMCR retains its flexible property and uses heterogeneity from both top groups (e.g. location-time) and individuals to detect better latent classes in the data. NLCMCR clusters individuals into more accurate latent classes than LCMCR when the data has a nested structure. This is important in application because it’s useful to better summarize group properties based on individual characteristics. Meanwhile, the model NLCMCR gives similar population size estimation with the model LCMCR but with smaller uncertainty if the data has a nested structure. This narrows down uncertainty of our estimation of the true population size. It also reduces uncertainty of other parameter estimates, such as list capture probabilities in each latent class and latent class proportions. Overall, NLCMCR provides both impressed population size estimation and description of the capture patterns.