Estimates of abundance (\(N\)) are commonly of interest in both basic and applied ecological studies, but estimating \(N\) is much more difficult than you might initially expect. There are a variety of methods that can be used.
To estimate N with closed-population, mark-recapture surveys, one typically conducts several capture events in a short window of time during which the population is assumed to be closed to changes in abundance, i.e., birth, death, emigration, or immigration events don’t take place during the sampling window. On any given capture event, one can expect to capture some proportion of the actual population and the number of unique individuals caught (the count, \(C\)) is expected to be a function of population size, \(N\), and the probability of capture, \(p\), for the given species, study location, and methods being used. Captures can be physical captures or might be sightings (e.g., from camera traps) or DNA samples (e.g., from tissue samples such as hair).
\[E(C) = N \cdot p\]
The key to successful estimation of abundance relies on estimating \(p\) well. And the canonical estimator is:
\[ \hat N = \frac{C}{\hat p}\]
Ball and urn studies are commonly used as a motivating example and work as follows. Reach in to an urn with 100 white balls and remove a sample (\(n_1 = 30\)). Mark the balls in the sample and replace them. Take another random sample (\(n_2 = 36\)) and count the number of marked (\(m_2 = 10\)) and unmarked balls (\(u_2 = 26\)).
We expect the proportion of marked balls to total balls in the population to be the same as the proportion of marked and unmarked balls observed in our sample.
\[\frac{n_1}{\hat N} = \frac{m_2}{n_2}, e.g., \frac{30}{\hat N} = \frac{10}{36}\] \[\hat N = \frac{n_1 \cdot n_2}{m_2}\]
This relationship is the basis of the Lincoln-Petersen Estimator for 2 occasions, which for our example above is as follows.
\[\hat N = \frac{30 \cdot 36}{10} = 108\]
You can re-arrange this equation to see that \(\hat N = \frac{C}{\hat p}\). If you consider \(n_1\) to be \(C\) and you estimate \(\hat p_1 = \frac{n_1}{\hat N} = \frac{m_2}{n_2}\)
\[P(n_1, n_2, m_2 | N, p_1, p_2) = \frac{N!}{m_2!(n_1-m_2)!(n_2-m_2)!(N-r)!} (p_1 p_2)^{m_2} (p_1 q_2)^{n_1-m_2} (q_1 p_2)^{n_2-m_2}(q_1 q_2)^{N-r}\]
In the equation above, \(q_i = 1- p_i\), \(r = n_1 + n_2 - m_2\), the number of unique individuals captured over the 2 occasions, and \(N-r\) is the number of animals that were not caught. Notice that \(N-r\) is the key quantity you need to estimate in an abundance study (more on that below). Also, note that we will refer to \(r\) more commonly in our work and in MARK as \(M_{t+1}\), i.e., the number of unique animals marked when the study is completed.
Closed-form maximum likelihood estimators (MLEs) for a 2-occasion study are as follows.
\(\hat p_1 = \frac{m_2}{n_2} = \frac{n_1}{\hat N}\)
\(\hat p_2 = \frac{m_2}{n_1} = \frac{n_2}{\hat N}\)
Biased-adjusted estimator for \(N\): \(\hat N = \frac{(n_1 + 1)(n_2+1)}{m_2+1}-1\)
\(\widehat{var}(\hat N) = \frac{(n_1+1)(n_2+1)(n_1-m_2)(n_2-m_2)}{(m_2+1)^2(m_2+2)^2}\)
Example
Using the probability equation above for a situation with the following known parameter values: \(N=250, p_1 = 0.6, p_2 = 0.3\), the expected frequencies of each encounter history are as follows.
Encounter History | E(frequency) | Probability |
---|---|---|
11 | 45 | \(N p_1 p_2 = 250 \cdot 0.6 \cdot 0.3\) |
10 | 105 | \(N p_1 (1-p_2)= 250 \cdot 0.6 \cdot (1-0.3)\) |
01 | 30 | \(N (1-p_1) p_2)= 250 \cdot (1-0.6) \cdot 0.3\) |
00 | 70 | \(N (1-p_1) (1-p_2)= 250 \cdot (1-0.6) \cdot (1-0.3)\) |
Notice that for this situation,
\(n_1 = 45 + 105 = 150\)
\(n_2 = 45 + 35 = 80\)
\(m_2 = 45\)
\(\frac{m2}{n_2} = \frac{45}{75} = 0.6 = p_1\)
\(\frac{m2}{n_1} = \frac{45}{150} = 0.3 = p_2\)
Population is closed to additions (birth & immigration) and to losses (death & emigration)
Marks are not lost or overlooked or misread by researchers
All animals are equally likely to be captured in each sample (this is the assumption that we will almost always need to relax). We’ll relax it to be: capture probabilities are appropriately modeled.
For most real biological applications, we need to:
Consider possible heterogeneity in capture probability,
Carefully consider the concept of closure, and
Use more than 2 occasions.
Essentially, there are 3 broad sources of variation that might influence \(p\) as nicely laid out in the Wildlife Monograph by Otis et al. (1978).
temporal variation in capture probability (‘t’),
behavioral responses to trapping (‘b’), where animals might become trap happy or trap shy after an initial capture/handling event, and
heterogeneity in capture probabilities for different individuals (‘h’), where different individuals are repeatedly harder or easier to catch than others.
As you’ve seen when reading chapter 14 of CW, which was written by Paul Lukacs, University of Montana, there is a great variety of approaches to estimating abundance with mark-recapture methods. We are going to focus on the full likelihood approaches to evaluating possible sources of variation in \(p\) and \(c\) while estimating \(N\). In all the approaches we’ll use, we’ll work with individually marked animals so that individual encounter histories are known, which is important to evaluating evidence for (1) behavioral responses to capture and (2) individual heterogeneity in capture probability.
As explained on page 14-4 of CW, we can express the likelihood as follows.
\[\mathcal{L}(N,\boldsymbol{p}|data) = \frac{N!}{(N-M_{t+1})!} \prod_{h} P[h]^{n_h} \cdot P[{\text {not encountered}}]^{N-M_{t+1}}\],
where \(\boldsymbol{p}\) is the vector of values for capture probability, \(p\), \(M_{t+1}\) is the number of unique individuals captured (or the minimum value that you know for \(N\) by the end of the study), \(h\) refers to an encounter history, e.g., 101, and \(n_h\) is the number of individuals with encounter history \(h\).
For reasons explained on page 14-4 of CW, the likelihoods we’ll use in Program MARK will use \(f_0\) in place of \(N-M_{t+1}\), where \(f_0\) is the number of individuals in the population that were never captured, i.e., the number (or frequency) of individuals that have histories that are all 0. Our goal in modeling mark-recapture data for a closed population is focused on estimating \(f_0\), which we add to \(M_{t+1}\) to obtain \(\hat N\), i.e., \(\hat N = M_{t+1} + \hat f_{0}\).
In MARK, we’ll use the log link function for \(f_0\), which constrains the values of \(\hat f_{0}\) to be \(\ge 0\). \(\hat N\) is presented as a derived parameter in model output as it’s derived as \(\hat N = M_{t+1} + \hat f_{0}\). Because \(M_{t+1}\) is known and we estimated \(f_0\), \(\widehat{var}(\hat N) = \widehat{var}(\hat {f}_0)\). With this construction, your values for \(\hat N\) will never be \(<M_{t+1}\), which is sensible as \(M_{t+1}\) is, as noted above and repeated here for emphasis, the minimum value that you know for \(N\) by the end of the study. The likelihood with \(f_0\) substituted for \(N-M_{t+1}\) is as follows.
\[\mathcal{L}(N,\boldsymbol{p}|data) = \frac{N!}{(N-M_{t+1})!} \prod_{h} P[h]^{n_h} \cdot P[{\text {not encountered}}]^{f_0}\]
With this general form of a likelihood in hand, you can then use maximum-likelihood estimation for a given model structure and data set to obtain estimates of capture probabilities and of \(f_0\). You can also try other model structures and compare the \(AIC_c\) values for competing models of interest. For example, you might be concerned that individuals become trap shy after being captured and handled. If so, you could run a model in which the probability of initial capture, \(p\), and the probability of being recaptured, \(c\), are allowed to differ from one another. You could then compare the \(AIC_c\) score for that model to the score for a model that constrains \(p=c\) but that allows \(p\) to vary by capture occasion.
\[P(h = 101) = p \cdot (1-p) \cdot p\]
\[P(h = 101) = p_1 \cdot (1-p_2) \cdot p_3\]
\[P(h = 101) = p \cdot (1-c) \cdot c\]
Note that for \(M_b\) recapture information is used only in the estimation of \(c\) (a nuisance parameter). That is, recaptures don’t provide information with respect to estimation of \(p\) or \(f_0\) (or therefore, \(N\). Thus, once \(M_b\) has been identified as appropriate, the estimation is done as if the study were a removal study where one attempts to deplete the population’s number of unmarked individuals over multiple capture occasions. It is critical in this type of study that the number of new animals captured on each occasion decreases over successive occasion. If this does not occur sufficiently, the model can fail to produce results. In particular, failure can readily occur if you have temporal changes in \(p\) such that few individuals are caught early in the study and many new individuals are caught late in the study.
4a. \(M_{bh}\) - generalized removal estimator. Recall that for \(M_b\), estimation of \(N\) is based on information on initial captures, and thus, \(M_b\) is a form of removal estimator. If we have heterogeneity, we expect (1) the average probability of first capture to decrease with each occasion and (2) the most rapid decrease to occur over the first few occasions. Fundamentally, you are depleting the pool of unmarked individuals that are most catchable as the study progresses. On the first occasion, no animal has been caught before so \(\bar p_1\) is high relative to, say, \(\bar p_6\) because by occasion 6 most of the animals that are easy to catch have already been caught. With this approach, you try to fit the data using a version of \(M_b\) with various numbers of \(\bar p's\). The idea is to fit a sequence of increasingly general models to the average value for initial capture probability and to evaluate the results.
Step 1. Evaluate the fit of \(M_b\), which uses a single average p for all occasions: \(\bar p_1 = \bar p_2 = \bar p_3 = \ldots \bar p_t\).
Step 2. Allow \(\bar p_1\) to be different from the others, i.e., you are catching most of the highly catchable animals on the 1st occasions, and the remaining pool of unmarked animals is quite homogeneous with respect to \(p\), e.g., \(\bar p_1 \ne \bar p_2 = \bar p_3 = \ldots \bar p_t\).
Step 3, 4, … Generalize further, e.g., \(\bar p_1 \ne \bar p_2 \ne \bar p_3 = \ldots \bar p_t\), and use \(AIC_c\) scores to compare versions.
4b. and 4c. - Finite mixture models: Full likelihood models have been developed that model the heterogeneity in \(p\) as coming from a finite mixture of values. For example, imagine a situation without time variation in \(p\) and a population that consists of 2 types of individuals: 1 type with \(p_A\) and another with \(p_B\). To estimate \(N\) for this situation, you estimate what proportion of the population is in mixture A, designated by a parameter \(\pi\), and what proportion is of type B, which will be \(1-\pi\). You can also try models with a mixture of >2 types of individuals, e.g., \(\pi_A, \pi_B\), and \(\pi_C\), where \(\pi_C = 1 - (\pi_A + \pi_B)\). I designated this set as “4b. and 4c.” because in Program MARK, you can choose to work with heterogeneity models that either include \(p\) and \(c\) or just \(p\) (no behavioral response to initial capture). This class of models can be quite useful as a way of incorporating heterogeneity in \(p\) and are well worth learning about if you venture into abundance estimation work in your career.
Quoting from page 14-9 of CW, “A subtlety of the closed population models is that the last \(p\) parameter is not identifiable unless a constraint is imposed. When no constraint is imposed on the last \(p_i\), the likelihood is maximized with the last \(p=1\), giving the estimate \(\hat{N} = M_{t+1}\).”
That quote is followed a bit further down with this important statement on page 14-10. “Thus, it is diagnostic to check to see whether \(\hat{N} = M_{t+1}\), and if so, to see if the last \(p_i\) estimate equals 1. If they are, then you’ve forgotten to constrain (the last) p.”
In some of the model structures reviewed above, the constraint exists and no problems exist. For \(M_0\) and \(M_b\), all the \(p_i\) are set equal to one another. For \(M_t\), the constraint is achieved because \(p_i = c_i\). If you run \(M_{t,b}\), you need to model the \(p_i\) and \(c_i\) using an additive structure in which there is a constant offset between the 2.
The following quote from Paul Lukacs on page 14-34 of CW is important to consider:
“the single best way to minimize the bias caused by individual heterogeneity is to get \(p\) as high as possible – the ‘big law’ of capture-recapture design. When \(p\) is high there is little room for variation and little chance that an individual is not detected. Bill Link (2003, 2004) demonstrated that different models of the form of individual heterogeneity can lead to very different estimates of abundance and fit the data equally well. The magnitude of the differences in abundance estimates is related to \(p\); when \(p\) is small the differences can be large. Therefore,to have much hope of estimating abundance with little bias, capture probability must be relatively high.”
It’s also important to realize that you can’t investigate sources of variation in \(p\) without having >2 occasions and some of the methods will require more. For example, quoting again from page 14-34 of CW, “A lower level of >5 occasions is likely necessary to achieve reasonable performance for either the continuous- or discrete-mixture approaches(White & Cooch 2017).”
Link, W.A. 2003. Non-identifiability of population size from capture-recapture data with heterogeneous detection probabilities. Biometrics 59:1125-1132.
Link, W.A. Individual heterogeneity and identifiability in capture–recapture models. 2004. Animal Biodiversity and Conservation 27:87–91.
Otis, D.L., K.P. Burnham, G.C. White, and D.R. Anderson. 1978. Statistical Inference from Capture Data on Closed Animal Populations. Wildlife Monographs 62:3-135.
White, G.C., and E.G. Cooch. 2017. Population abundance estimation with heterogeneous encounter probabilities using numerical integration. Journal of Wildlife Management 81:322-336.