Rejoinder to the discussion on "The central role of the identifying assumption in population size estimation"

We organize the discussants' major comments into the following categories: sensitivity analyses, zero counts, model selection, the marginal no-highest-order interaction (NHOI) assumption, and the usefulness of our proposed framework.

Discussion on "The central role of the identifying assumption in population size estimation" by Serge Aleshin-Guendel, Mauricio Sadinle, and Jon Wakefield

In this discussion response, we consider some practical implications of the authors' consideration of the no-highest-order interaction (NHOI) model for multiple systems estimation, which permits the authors to derive the explicit (albeit untestable) identifying assumption related to the unobserved (or missing) individuals. In particular, we discuss several aspects, from the standard process of model selection to potential poor predictive performance due to over-fitting and the implications of data reduction. We discuss these aspects in relation to the case study presented by the authors relating to the number of civilian casualties within the Kosovo war, and conduct further preliminary simulations to investigate these issues further. The results suggest that the NHOI models considered, despite having a potentially useful theoretical result in relation to the underlying identifying assumption, may perform poorly in practice.

shinyrecap: A Shiny Application for Population Size Estimation from Capture-Recapture Data (Preprint)

Background

Population size estimates (PSE) provide critical information in determining resource allocation for HIV services geared toward those at high risk of HIV, including female sex workers, men who have sex with men, and people who inject drugs. Capture-recapture (CRC) is often used to estimate the size of these often-hidden populations. Compared with the commonly used 2-source CRC, CRC relying on 3 (or more) samples (3S-CRC) can provide more robust PSE but involve far more complex statistical analysis.

Objective

This study aims to design and describe the Shiny application (shinyrecap), a user-friendly interface that can be used by field epidemiologists to produce PSE.

Methods

shinyrecap is built on the Shiny web application framework for R. This allows it to seamlessly integrate with the sophisticated CRC statistical packages (eg, Rcapture, dga, LCMCR). Additionally, the application may be accessed online or run locally on the user’s machine.

Results

The application enables users to engage in sample size calculation based on a simulation framework. It assists in the proper formatting of collected data by providing a tool to convert commonly used formats to that used by the analysis software. A wide variety of methodologies are supported by the analysis tool, including log-linear, Bayesian model averaging, and Bayesian latent class models. For each methodology, diagnostics and model checking interfaces are provided.

Conclusions

Through a use case, we demonstrated the broad utility of this powerful tool with 3S-CRC data to produce PSE for female sex workers in a subnational unit of a country in sub-Saharan Africa.

A need for speed in Bayesian population models: a practical guide to marginalizing and recovering discrete latent states

Bayesian population models can be exceedingly slow due, in part, to the choice to simulate discrete latent states. Here, we discuss an alternative approach to discrete latent states, marginalization, that forms the basis of maximum likelihood population models and is much faster. Our manuscript has two goals: (1) to introduce readers unfamiliar with marginalization to the concept and provide worked examples and (2) to address topics associated with marginalization that have not been previously synthesized and are relevant to both Bayesian and maximum likelihood models. We begin by explaining marginalization using a Cormack-Jolly-Seber model. Next, we apply marginalization to multistate capture-recapture, community occupancy, and integrated population models and briefly discuss random effects, priors, and pseudo-R² . Then, we focus on recovery of discrete latent states, defining different types of conditional probabilities and showing how quantities such as population abundance or species richness can be estimated in marginalized code. Last, we show that occupancy and site-abundance models with auto-covariates can be fit with marginalized code with minimal impact on parameter estimates. Marginalized code was anywhere from five to >1,000 times faster than discrete code and differences in inferences were minimal. Discrete latent states and fully conditional approaches provide the best estimates of conditional probabilities for a given site or individual. However, estimates for parameters and derived quantities such as species richness and abundance are minimally affected by marginalization. In the case of abundance, marginalized code is both quicker and has lower bias than an N-augmentation approach. Understanding how marginalization works shrinks the divide between Bayesian and maximum likelihood approaches to population models. Some models that have only been presented in a Bayesian framework can easily be fit in maximum likelihood. On the other hand, factors such as informative priors, random effects, or pseudo-R² values may motivate a Bayesian approach in some applications. An understanding of marginalization allows users to minimize the speed that is sacrificed when switching from a maximum likelihood approach. Widespread application of marginalization in Bayesian population models will facilitate more thorough simulation studies, comparisons of alternative model structures, and faster learning.

Population size estimation with interval censored counts and external information: Prevalence of multiple sclerosis in Rome

We discuss Bayesian log-linear models for incomplete contingency tables with both missing and interval censored cells, with the aim of obtaining reliable population size estimates. We also discuss use of external information on the censoring probability, which may substantially reduce uncertainty. We show in simulation that information on lower bounds and external information can each improve the mean squared error of population size estimates, even when the external information is not completely accurate. We conclude with an original example on estimation of prevalence of multiple sclerosis in the metropolitan area of Rome, where five out of six lists have interval censored counts. External information comes from mortality rates of multiple sclerosis patients.

Multiple Systems Estimation (or Capture-Recapture Estimation) to Inform Public Policy

Estimating population sizes has long been of interest, from the estimation of the human or ecological population size within regions or countries to the hidden number of civilian casualties in a war. Total enumeration of the population, for example, via a census, is often infeasible or simply impractical. However, a series of partial enumerations or observations of the population is often possible. This has led to the ideas of capture-recapture methods, which have been extensively used within ecology to estimate the size of wildlife populations, with an associated measure of uncertainty, and are most effectively applied when there are multiple capture occasions. Capture-recapture ideology can be more widely applied to multiple data-sources, by the linkage of individuals across the multiple lists. This is often referred to as Multiple Systems Estimation (MSE). The MSE approach has been preferred when estimating "capture-shy" or hard-to-reach populations, including those caught up in the criminal justice system; or homeless; or trafficked; or civilian casualties of war. Motivated by a range of public policy applications of MSE, each briefly introduced, we discuss practical problems with potentially substantial methodological implications. They include: "period" definition; "case" definition; when an observed count is not a true count of the population of interest but an upper bound due to mismatched definitions; exact or probabilistic matching of "cases" across different lists; demographic or other information about the "case" which may influence capture-propensities; required permissions to access extant-lists; list-creation by research-teams or interested parties; referrals (if presence on list A results - almost surely - in presence on list B); different mathematical models leading to widely different estimated population sizes; uncertainty in estimation; computational efficiency; external validation; hypothesis-generation; and additional independent external information. Returning to our motivational applications, we focus on whether the uncertainty which qualified their estimates was sufficiently narrow to orient public policy; and, if not, what options were available and/or taken to reduce the uncertainty or to seek external validation. We also consider whether MSE was hypothesis-generating: in the sense of having spawned new lines of inquiry.

Estimation of a closed population size of tadpoles in temporary pond

The practice of capture-recapture to estimate the diversity is well known to many animal groups, however this practice in the larval phase of anuran amphibians is incipient. We aimed at evaluating the Lincoln estimator, Venn diagram and Bayes theorem in the inference of population size of a larval phase anurocenose from lotic environment. The adherence of results was evaluated using the Kolmogorov-Smirnov test. The marking of tadpoles for later recapture and methods measurement was made with eosin methylene blue. When comparing the results of Lincoln-Petersen estimator corresponding to the Venn diagram and Bayes theorem, we detected percentage differences per sampling, i.e., the proportion of sampled anuran genera is kept among the three methods, although the values are numerically different. By submitting these results to the Kolmogorov-Smirnov test we have found no significant differences. Therefore, no matter the estimator, the measured value is adherent and estimates the total population. Together with the marking methodology, which did not change the behavior of tadpoles, the present study helps to fill the need of more studies on larval phase of amphibians in Brazil, especially in semi-arid northeast.

MCMC model determination for discrete graphical models

In this paper we compare two alternative MCMC samplers for the Bayesian analysis of discrete graphical models; we present both a hierarchical and a nonhierarchical version of them. We first consider the MC3 algorithm by Madigan and York (1995) for which we propose an extension that allows for a hierarchical prior on the cell counts. We then describe a novel methodology based on a reversible jump sampler. As a prior distribution we assign, for each given graph, a hyper-Dirichlet distribution on the matrix of cell probabilities. Two applications to real data are presented.

A default prior distribution for contingency tables with dependent factor levels

A default prior distribution is proposed for the Bayesian analysis of contingency tables. The prior is specified to allow for dependence between levels of the factors. Different dependence structures are considered, including conditional autoregressive and distance correlation structures. To demonstrate the prior distribution, a dataset is considered which involves estimating the number of injecting drug users in the eleven National Health Service board regions of Scotland using an incomplete contingency table where the dependence structure relates to geographical regions.

Incomplete contingency tables with censored cells with application to estimating the number of people who inject drugs in Scotland

Estimating the size of hidden or difficult to reach populations is often of interest for economic, sociological or public health reasons. In order to estimate such populations, administrative data lists are often collated to form multi-list cross-counts and displayed in the form of an incomplete contingency table. Log-linear models are typically fitted to such data to obtain an estimate of the total population size by estimating the number of individuals not observed by any of the data-sources. This approach has been taken to estimate the current number of people who inject drugs (PWID) in Scotland, with the Hepatitis C virus diagnosis database used as one of the data-sources to identify PWID. However, the Hepatitis C virus diagnosis data-source does not distinguish between current and former PWID, which, if ignored, will lead to overestimation of the total population size of current PWID. We extend the standard model-fitting approach to allow for a data-source, which contains a mixture of target and non-target individuals (i.e. in this case, current and former PWID). We apply the proposed approach to data for PWID in Scotland in 2003, 2006 and 2009 and compare with the results from standard log-linear models. © 2013 The Authors. Statistics in Medicine published by John Wiley & Sons, Ltd.

Injecting drug users in Scotland, 2006: Listing, number, demography, and opiate-related death-rates

Using Bayesian capture-recapture analysis, we estimated the number of current injecting drug users (IDUs) in Scotland in 2006 from the cross-counts of 5670 IDUs listed on four data-sources: social enquiry reports (901 IDUs listed), hospital records (953), drug treatment agencies (3504), and recent Hepatitis C virus (HCV) diagnoses (827 listed as IDU-risk). Further, we accessed exact numbers of opiate-related drugs-related deaths (DRDs) in 2006 and 2007 to improve estimation of Scotland's DRD rates per 100 current IDUs. Using all four data-sources, and model-averaging of standard hierarchical log-linear models to allow for pairwise interactions between data-sources and/or demographic classifications, Scotland had an estimated 31700 IDUs in 2006 (95% credible interval: 24900-38700); but 25000 IDUs (95% CI: 20700-35000) by excluding recent HCV diagnoses whose IDU-risk can refer to past injecting. Only in the younger age-group (15-34 years) were Scotland's opiate-related DRD rates significantly lower for females than males. Older males' opiate-related DRD rate was 1.9 (1.24-2.40) per 100 current IDUs without or 1.3 (0.94-1.64) with inclusion of recent HCV diagnoses. If, indeed, Scotland had only 25000 current IDUs in 2006, with only 8200 of them aged 35+ years, the opiate-related DRD rate is higher among this older age group than has been appreciated hitherto. There is counter-balancing good news for the public health: the hitherto sharp increase in older current IDUs had stalled by 2006.

Disease burden of prosthetic joint infections after hip and knee joint replacement in Finland during 1999-2004: capture-recapture estimation

We evaluated the Finnish Hospital Infection Program (SIRO) conducting incidence surveillance for prosthetic joint infection (PJI) from 1999 to 2004. We estimated its sensitivity using capture-recapture methods and assessed the disease burden of PJIs after hip (THA) and knee (TKA) arthroplasties (N = 13 482). The following three data sources were used: SIRO, the Finnish Arthroplasty Register (FAR), and the Finnish Patient Insurance Center (FPIC), which were cross-matched, and 129 individual PJIs were identified. After adjusting for the positive predictive value of SIRO (91%) a log-linear model including an interaction term between FAR and FPIC provided an estimated PJI rate of 1.6% [95% confidence interval (CI): 1.2-2.4] for THA and 1.3% (1.1-1.6) for TKA. Sensitivity for SIRO varied from 36% to 57%. The annual disease burden was 2.1 PJIs per 100 000 population after THA and 1.5 after TKA. The true disease burden of PJIs may be heavier than the rates from national sentinel surveillance systems usually suggest.

Assessing threats to species at risk using stage-structured state-space models: mortality trends in skate populations

Population models are needed to assess the threats to species at risk and to evaluate alternative management actions. Data to support modeling is limited for many species at risk, and commonly used approaches generally assume stationary vital rates, a questionable assumption given widespread ecosystem change. We describe a modeling approach that can be applied to time series of length composition data to estimate vital rates and test for changes in these rates. Our approach uses stage-structured population models fit within a Bayesian state-space model. This approach simultaneously allows for both process and observation uncertainty, and it facilitates incorporating prior information on population dynamics and on the monitoring process. We apply these models to populations of winter skate (Leucoraja ocellata) that have been designated as "endangered" or "threatened." These models indicate that natural mortality has decreased for juveniles and increased for adults in these populations. The declines observed in these populations had been attributed to unsustainable rates of bycatch in fisheries for other groundfishes; our analyses indicate that increased natural mortality of adults is also an important factor contributing to these declines. Adult natural mortality was positively related to grey seal (Halichoerus grypus) abundance, suggesting the hypothesis that increased adult mortality reflected increased predation by expanding grey seal herds. Population projections indicated that the threatened population would be expected to stabilize at a low level of abundance if all fishery removals were eliminated, but that the endangered population would likely continue to decline even in the absence of fishery removals. We note that time series of size distributions are available for most marine fish populations monitored by research surveys, and we suggest that a similar approach could be used to extract information from these time series in order to estimate mortality rates and changes in these rates.

Capture-recapture estimation using finite mixtures of arbitrary dimension

Reversible jump Markov chain Monte Carlo (RJMCMC) methods are used to fit Bayesian capture-recapture models incorporating heterogeneity in individuals and samples. Heterogeneity in capture probabilities comes from finite mixtures and/or fixed sample effects allowing for interactions. Estimation by RJMCMC allows automatic model selection and/or model averaging. Priors on the parameters stabilize the estimates and produce realistic credible intervals for population size for overparameterized models, in contrast to likelihood-based methods. To demonstrate the approach we analyze the standard Snowshoe hare and Cottontail rabbit data sets from ecology, a reliability testing data set.

Estimating current injectors in Scotland and their drug-related death rate by sex, region and age-group via Bayesian capture--recapture methods

Using Bayesian capture-recapture methods, we estimate current injectors in Scotland in 2003, and, thereby, injectors' drug-related death rates for the period 2003-2005. Four different data sources are considered [Hepatitis C Virus (HCV) database, hospital admissions, social enquiry reports, and drug misuse database reports by General Practices or Drug Treatment Agencies] which provide covariate information on sex, region (Greater Glasgow versus elsewhere in Scotland) and age group (15-34 years and 35+ years).We quantified Scotland's current injectors in 2003 at 27,400 (95% highest probability density interval: 20,700-32,100) by incorporating underlying model uncertainty in terms of the possible interactions present between data sources and/or covariates. The posterior probability was 72% that Scotland had more current injectors in 2003 than in 2000. Detailed comparison with 2000 gave evidence of importantly changed numbers of current injectors for different covariate classes.In addition, and of particular social interest, is the estimation of injectors' drug-related death rates. Expert information was used to construct upper and lower bounds on the number of drug-related deaths pertaining to injectors, which were then used to provide bounds on injectors' drug-related death rates. Failure to incorporate expert information could result in over-estimation of drug-related death rates for subclasses of injectors.

A Bayesian hierarchical model for the estimation of two incomplete surveillance data sets

A model-based approach to analyze two incomplete disease surveillance datasets is described. Such data typically consist of case counts, each originating from a specific geographical area. A Bayesian hierarchical model is proposed for estimating the total number of cases with disease while simultaneously adjusting for spatial variation. This approach explicitly accounts for model uncertainty and can make use of covariates.The method is applied to two surveillance datasets maintained by the Centers for Disease Control and Prevention on Rocky Mountain spotted fever (RMSF). An inference is drawn using Markov Chain Monte Carlo simulation techniques in a fully Bayesian framework. The central feature of the model is the ability to calculate and estimate the total number of cases and disease incidence for geographical regions where RMSF is endemic.The information generated by this model could significantly reduce the public health impact of RMSF and other vector-borne zoonoses, as well as other infectious or chronic diseases, by improving knowledge of the spatial distribution of disease risk of public health officials and medical practitioners. More accurate information on populations at high risk would focus attention and resources on specific areas, thereby reducing the morbidity and mortality caused by some of the preventable and treatable diseases.

A Bayesian Approach to the Multistate Jolly-Seber Capture-Recapture Model

This article considers a Bayesian approach to the multistate extension of the Jolly-Seber model commonly used to estimate population abundance in capture-recapture studies. It extends the work of George and Robert (1992, Biometrika79, 677-683), which dealt with the Bayesian estimation of a closed population with only a single state for all animals. A super-population is introduced to model new entrants in the population. Bayesian estimates of abundance are obtained by implementing a Gibbs sampling algorithm based on data augmentation of the missing data in the capture histories when the state of the animal is unknown. Moreover, a partitioning of the missing data is adopted to ensure the convergence of the Gibbs sampling algorithm even in the presence of impossible transitions between some states. Lastly, we apply our methodology to a population of fish to estimate abundance and movement.

Prior information in behavioral capture-recapture methods: demographic influences on drug injectors' propensity to be listed in data sources and their drug-related mortality

The authors present findings from a Bayesian analysis of Scotland's four primary capture-recapture data sources for 2000 that was carried out to estimate numbers of current injecting drug users by region (Greater Glasgow vs. elsewhere in Scotland), sex (male vs. female), and age group (15-34 years vs. > or =35 years). A secondary goal of the analysis was to obtain Bayesian estimates and credible intervals for the demographic influences on Scotland's drug-related death rate per 100 current injectors. Incorporation of informative priors altered the models with highest posterior probability. Expert opinion on how demography influenced Scottish drug injectors' propensity to be listed in different data sources was taken into account, along with external information about European injectors' drug-related death rates and male:female ratios. Higher drug-related mortality was confirmed in older drug injectors and those outside of Greater Glasgow. Female injectors' lower drug-related death rate was not sustained beyond 34 years of age. The authors recommend that demographic influences be accommodated in behavioral capture-recapture estimation, especially when it is a prelude to secondary analysis, such as the analysis of drug-related death rates presented here.

Bayesian Population Dynamics of Interacting Species: Great Gerbils and Fleas in Kazakhstan

We propose a discrete-time Bayesian hierarchical model for the population dynamics of the great gerbil-flea ecological system. The model accounts for the sampling variability arising from data originally collected for other purposes. The prior for the unknown population densities incorporates specific biological hypotheses regarding the interacting dynamics of the two species, as well as their life cycles, where density-dependent effects are included. Posterior estimates are obtained via Markov chain Monte Carlo. The variance of the observed density estimates is a quadratic function of the unknown density. Our study indicates the presence of a density-dependent growth rate for the gerbil population. For the flea population there is clear evidence of density-dependent over-summer net growth, which is dependent on the flea-to-gerbil ratio at the beginning of the reproductive summer. Over-winter net growth is favored by high density. We estimate that on average 35% of the gerbil population survives the winter. Our study shows that hierarchical Bayesian models can be useful in extracting ecobiological information from observational data.