Outline and notation

We consider the following situation (see Figure 1 for illustration): The target population (consisting of N individuals) is followed over time during a time interval of total duration T. At each time-point a given individual may or may not be registered. We use ‘registered’ in a general sense here, it can for example refer to an instance of care at a medical facility; an arrest by the police; the purchase of a medically prescribed drug, etc. We assume that each individual in the target population, if not previously registered, has a fixed probability π of being registered at each time-point. We will refer to π as registration probability in the following. A registration (independent of the nature of the registry), must contain an identifier of the person (PIN) and the time-point of registration. We note that there might be a delay between the occurrence of the event we would like to track, and the time point of registration of this event. Consider for example an attempt to estimate the prevalence of a certain disease in the population. If health care records are used to estimate this, there might be a delay between the time of acquiring the disease and the first contact with the health care system. We will assume that such delays, if they occur, do not introduce systematic errors in the analysis and hence can be ignored.

Next we will detail how such register data can be used to estimate the total size N of the population. First we will describe the novel ML approach and then an approach that has been previously used with data of thistype.

Maximum likelihood estimation of N

Here we describe how the data in the register can be used to estimate N using ML (see Figure 1 for illustration). First the observation interval is divided into M non-overlapping epochs with equal duration d. That is, the j-th epoch, I _j is given by I _j =[t ₀+d(j−1),t ₀+d j) where t ₀ denotes the time point when the observations start and j=1,2,…,M. We will assume that time is discrete and that there are d opportunities to be registered in an interval of duration d. For each of the M epochs we count the number of new persons being registered in that epoch. That is, we count the persons registered during epoch I _j that were not registered in any of the previous epochs. These numbers will be called n ₁,n ₂,…,n _M and are used in the following to represent both random variables and samples of these random variables. It is n ₁,n ₂,…,n _M that constitute the data we will use to estimate the total number of persons in the population (i.e. N). The number of intervals M>1 is a parameter in the suggested method and we will study how the estimate of N depend on M below.

Given the log-likelihood function (Eq. 2) we take as our estimate of population size the value of N that maximizes that expression (henceforth denoted N _{M L} ). Note that we simultaneously estimate the registration probability p, but since this is not the parameter of main interest, we will focus on the estimator of N.

For the case M=2 the maximum of the likelihood function is given by $N_{\mathit{\text{ML}}2}=n_1^2/(n_1-n_2)$. However, in the general case, the maximum likelihood estimates are easiest found by maximizing Eq. 2 numerically, for example using Newton’s method. It should be noted that the ML estimator is not applicable for all samples. For the case M=2, for example, n ₁ must be greater than n ₂. General conditions for when the ML estimator is applicable have been derived before [18]. As we will show later, in cases where the ML estimator give good estimates of N, it is very unlikely to get a sample for which the ML estimator does not exist.

Truncated poisson estimate of N

Given a single source of data with possible multiple registrations per person there are alternative estimates of populations size that can be (and have been) used. The simplest of these estimates can be derived if one assumes that the sampled data follow a zero-truncated Poisson distribution (this is just a re-weighted standard Poisson distribution in which the zero bin is not observed). We note that if the probability π of being registered at each time point is independent of previous registrations, and there are T time points in total, then the number of registrations for each subject is distributed as a binomial variable with parameters T and π. Given reasonable values of these parameters, the distribution of the number of registrations will closely follow a Poisson distribution with parameter λ=defT π (e.g. [19]). Consider for example that observations are made daily for a year, then T=365, and if the probability of being registered at least once in a year is 0.3 then π≃0.001 which gives λ≃0.36. These are values for which the Poisson approximation of the Binomial distribution is very good and, consequently, a zero truncated Poisson distribution is a natural choice to model the distribution of the registrations.

By the law of large numbers we know that N _C will be close to N when N is large. We will investigate the convergence further below. This estimator (N _C ) was first derived in [17], using a more general formalism and we will refer to N _C also as Chao’s estimator.

This estimate was derived in [10] using a different approach.

Monte Carlo simulations

To evaluate the performance of the estimators we performed Monte Carlo simulations. In these simulations we generated data for 512 ‘time steps’ (e.g. days) where each of N individuals has a probability p of being registered. Simulations where run for different values of N and different values of p (described in the Results). All simulation results were based on 5000 realizations. For the ML estimator the maximum of the likelihood (Eq. 2) was found numerically (for M>2) using Newtons method. In some cases the method did not converge, something that often was due to a failure of automatically specifying initial conditions within the region of convergence. In such cases a small change to the initial conditions is often sufficient for convergence. Since there is an explicit expression for the case M=2, the corresponding estimate is a good starting point for the algorithm. Although it is possible (in particular for small samples) to get data for which the ML estimator does not exists [18], for our purposes problems of convergence were very rare in all cases where the estimator will be useful in practice (see Results). An indication that a sample is not suited for the ML estimator is if n ₁≤n ₂ (for the case of M=2), and if so other methods of estimation must be used.

The ML estimator was derived assuming that all the individuals in the population have the same registration probability, something that translates to p in the simulations is the same for individuals. To model the effect of heterogeneity in the target population we run simulations where we allowed p to vary from individual to individual. In particular, the ps were randomly sampled from a normal distribution with a fixed mean value. The parameter used to control population heterogeneity was the standard deviation of the normal distribution. This way the effect of increasing population heterogeneity can be studied systematically from a constant p (zero standard deviation) to ps that vary a lot between individuals (large standard deviation). We also run simulations where we sampled ps from uniform distributions but the results were similar to those obtained using a normal distribution and are therefore not shown.

Empirical data

Apart from the evaluation done on simulated data, we also wanted to use the suggested method on ‘real’ data. Here it was deemed important to use a data set where N could be assumed known, since then we can evaluate how close to the true value the different estimates are. We approached this by using as our target ‘population’ the accidentally deceased opiate misusers (fatal overdoses) in Sweden during 2006–2011. Most such cases occur outside of hospitals and the causes of death are therefore investigated in detailed forensic autopsies (including a toxicological analysis). We can therefore use the central Swedish Causes of Death registry held by The National Board of Health and Welfare (Socialstyrelsen) to count these cases and get a reasonably precise value of N.

To identify the deceased persons making up the target population we proceeded as follows. We selected all persons between 23 and 60 years of age at the time of death who were listed in the causes of death registry with heroin or methadone poisoning (ICD-10 codes T40.1 and T40.3) as contributing cause of death. This choice of substances should be obvious. However, there are also a substantial number of cases with morphine poisoning (T40.2) but these were not included as they were judged to constitute a mixed group of individuals, where only some are opiate misusers in the sense relevant here. The upper age limit was to minimize the number of cases where methadone had been medically prescribed and used in pain management therapy. The lower age limit was to enable a search of five year medical history for all individuals (we are tacitly assuming that most younger persons in this population had started their opiate misuse by the time they were 18 years old). We furthermore excluded all cases that were judged to be suicides. After applying these criteria there were 486 persons remaining who can be assumed to have misused opiates on a regular basis and these 486 will constitute our target population. To evaluate the population size estimators we will now see how well we can estimate the size of this population from another set of data. For this we will use health care data as described next.

For each of the N=486 individuals in the population we extracted information on previous contacts with the Swedish health care system by using the National Patient Registry (held by The National Board of Health and Welfare). This registry contains all instances of inpatient care in Sweden with good coverage and contains the main diagnosis (classified according to the WHO ICD-10 classification system) and time points for care. For the target population we extracted for all occurrences of inpatient care under diagnoses indicative of substance misuse. In particular, we used the following ICD-10 diagnoses: F11-F16, F18, and F19 and restricted the search to five years immediately preceding the time of death. During these time intervals 262 of the 485 persons (54%) had been receiving care under these diagnoses. We note that this might be taken as support for the assumption that the 485 persons constitute a relatively homogeneous population with respect to substance misuse. There were a total of 1165 hospital records from these 262 persons and most individuals had more than one care occasion. These health-care records will be used below to estimate N.

Since the individuals were all deceased at the time of project initiation, ethics committee approval was not necessary according to the regional research ethics board (Etikprövningsnä mnden i Stockholm).

Variability and bias of the ML estimator

Ideally, an estimator of population size should always give the right answer. However, since the estimator is computed from data (assumed to be random) it will both be variable (i.e. will give different results for different realizations) and it might have a bias (i.e. the expected value of the estimator might differ from the true value). In this section we study how the variance and bias of the estimator depend on total population size (N) and the fraction of N that is sampled (observed).

To study the variability and bias of the estimator we made Monte Carlo simulations for N=500,1000,5000, and 50000 and for each N we varied the fraction of observed units from 0.1 to 0.8 in steps of 0.05. To characterize the variability we used the estimated coefficient of variation (CV) (sample standard deviation divided by sample mean). Results are shown in Figure 2A. For all Ns the CV starts with a high value and then decreases as the fraction sampled increases.

To investigate the bias of the estimator we divided the sample means from the simulations with their corresponding theoretical value (i.e. with the corresponding Ns). The resulting measure will be referred to as relative bias. Note that a relative bias of one means that the estimator is unbiased. Figure 2B shows the relative bias for the parameter values investigated. For sufficiently large fraction of sampled units, the bias goes to zero independently of population size N. However, for smaller fractions sampled, the bias can be substantial and for the two smallest populations (N=500 and N=1000) the bias is a non-monotonous function of sample size; it starts of with a substantial negative bias that later becomes positive before going to zero (relative bias going to one).

If we consider these bias and variance results together they show that whether the ML estimator is useful or not depends strongly on both N and the sample size. For example, if the true N is in the order of 50000 then already with a sampled fraction of about 15% (i.e. a sample size of 7500) very good estimates are typically obtained. If, on the other hand, the true N is around 500 the sampled fraction should not be less than 55% (sample size of 275) for the estimator to be reasonably accurate.

Population heterogeneity

So far we have assumed that the population is homogeneous with respect to the probability of being registered in each time interval (at least before being registered the first time). This translates to a constant p for all ‘subjects’ in the simulations. In the real world this assumption is obviously unrealistic and in this section we study the effects of relaxing it (see Methods).

Figure 3 shows the effect of making the units different with respect to the probability of being registered. Over the range studied, the heterogeneity does not influence the performance of the estimator greatly. Indeed, as shown in the inset in Figure 3B, the variability in registration probability between individuals can be substantial without leading to detrimental performance. Of course, if the variability in p (between units) become sufficiently large the performance will eventually break down. In particular, when many units start having zero probability of being registered the estimator will be biased.

Choosing M

In applications of the capture-and-remove method in population biology the number of observation occasions (captures) might often be determined by the design of the study. In epidemiology, on the other hand, we typically have a certain freedom to choose the number of ‘observation occasions’ (epochs, denoted M herein). Indeed, we might have access to daily registrations for a year or more and it is therefore of interest to see how dividing the total time interval into epochs of different length will influence the performance of the estimator. We investigated this by running Monte Carlo simulations for different values of M. We used M=2,4,8,16,32, and since the total time interval consisted in 512 observation units (‘days’) we used exactly the same data in all these five cases. Note that the case M=2 is special in that a closed-form expression exists whereas for M>2 the ML estimate has to be found iteratively.

We found that when the fraction of the total population sampled was large enough to give a relative variability (coefficient of variation) of 0.3 or smaller then the estimates obtained for M≥8 were virtually identical (correlation coefficients between estimates obtained for different values of M were ≥0.98). However, compared to the case of M=2, estimates differed and even if the correlation coefficient was relatively high (around 0.8) the estimates with M=2 had slightly larger coefficients of variation and a larger bias. Indeed, Figure 4 shows the simulation results for N=1000 and compares the performance of the estimator with M=2 to that with M=16. It is clear that the size of the effects of M are rather small but that in general an estimator with M≥8 should be used.

For smaller fractions of the population size the different estimators have less regular behavior, and correlation between different estimates might be as low as 0.35. Varying M and comparing the resulting estimates might therefore constitute a check that a sufficient fraction of the total population has been sampled.

Truncated poisson estimation

A simple alternative to the estimator we are suggesting when having data from one source is to assume that the samples come from a zero-truncated distribution and use all, or parts of the data to estimate the unobserved zero bin (e.g. [10]). For the type of data considered here a powerful and often used estimator goes under the name of Chao’s estimator. For the truncated Poisson distribution this estimator is easy to derive (Methods) but it was originally derived in a more general setting [17]. In this section we investigate the performance of this estimator (N _C ) and will see how it fails when we introduce simple dependencies between the probabilities of beingregistered.

We note that there are alternative ‘truncated Poisson’ estimators and one that has been applied was suggested in [10] and is also derived in the Methods. We made extensive tests using also this estimator and on our simulated data it behaved very similar to N _C but had a slightly larger variance. We will therefore focus on N _C in thefollowing.

Variability and bias of Chao’s estimator

The variability and bias of Chao’s estimator was determined analogously to that of the ML estimator. The results are shown in Figure 5. It is clear that when the Poisson assumption is fulfilled, Chao’s estimator outperforms the ML estimator. This is in particular the case when only a small fraction of the population is sampled. For the case N=500, for example, Chao’s estimator gives reasonably reliable estimates already when the sample size is around 125. However, when the assumptions are violated, Chao’s estimator can have a substantial bias, even when the fraction sampled islarge.

Systematic error in Chao’s estimator

One crucial assumption for any estimator relying on a truncated Poisson distribution is that the probability of being registered at any time is independent of previous registrations. In epidemiological data in general, and in health care records in particular, it is unlikely that this assumption holds. We will investigate this in an example data set below.

Here we will study the behavior of Chao’s estimator when the probability of being registered two or more times is different from the probability of being registered once. In the Monte Carlo simulations we started with all first-time registration probabilities (p ₁) equal and then once a unit was registered we changed (or not) the future registration probability p ₂. The effect of this will be that the distribution of the registrations no longer follow a Poisson distribution. As shown in Figure 6 the effect on Chao’s estimate of violating the Poisson assumption can be substantial. When a registered unit is less likely to become registered again, Chao’s estimator has a large positive bias, and this even for large sample sizes. When the probability of being registered more than once is larger than the probability of first registration, the estimator has a strong negative bias.

The reason for this large bias can be seen from the following considerations. Assume that that p ₂=p ₁/k, where k is an integer parameter. Then if N is large and p ₁ is small

$\frac{E{\left(h_1\right)}^2}{E\left(h_2\right)}\simeq \mathit{\text{kN}},$

which implies that Chao’s estimator will have an arbitrarily large positive bias.

Application to real data

We next apply N _{M L} and N _C to a real data set. We wanted to use the estimators on a population where we already knew (or at least had a very good approximation of) the size from other sources. This can be viewed as an attempt to validate the methods using real data; only by knowing the true size can we determine which of the methods that give the best estimate. As explained in Methods, the cases of fatal opiate overdoses in Sweden from 2006 to 2011 is a candidate for such a population. We took the total size of this population to be number of deceased as reported by The National Board of Health and Welfare. It totaled 486 persons. Now we will now see how well we can ‘estimate’ this number from the available health care records from the deceased persons. There were 262 of these individuals that had received inpatient care under the diagnoses of interest (see Methods) during this time.