Estimating, Testing, and Simulating Abundance in a Mark-Recapture Experiment with recapr

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Estimating, Testing, and Simulating Abundance in a Mark-Recapture Experiment with recapr

Matt Tyers

2021-09-07

Abundance Estimates and Standard Errors

In a simple two-event mark-recapture experiment, point estimates of abundance can be calculated according to the Chapman, Petersen, or Bailey estimators using NChapman(), NPetersen(), or NBailey(), given the sample sizes in the first and second events, and the number of individuals recaptured in the second event. The Petersen estimator is the Maximum Likelihood Estimator (MLE), but the Chapman estimator generally performs better and is often recommended. Both assume sampling without replacement in the second sampling event. By contrast, the Bailey estimator assumes sampling with replacement in the second event.

The argument n1 denotes the sample size in the first event, with n2 denoting the sample size in the second event and m2 denoting the number of marked (recaptured) individuals in the second event.

library(recapr)
NChapman(n1=100, n2=150, m2=20)    # Abundance estimate

## [1] 725.2381

vChapman(n1=100, n2=150, m2=20)    # estimated variance

## [1] 16348.22

seChapman(n1=100, n2=150, m2=20)   # standard error

## [1] 127.8601

Confidence Intervals

Confidence intervals can be generated for the Chapman, Petersen, or Bailey estimators using ciChapman(), ciPetersen(), or ciBailey(). A point estimate is returned as well as two confidence intervals, one using a Wald-type Normal-distribution assumption, and the other calculated by means of bootstrapping capture histories. The bootstrapping interval is likely to have more appropriate bounds than the Normal interval, and has demonstrated better coverage in testing via simulation.

ciChapman(n1=100, n2=150, m2=20)

## $Nhat
## [1] 725.2381
## 
## $ciNorm
## [1] 474.6368 975.8394
## 
## $ciBoot
## [1]  524.8966 1172.1538

Simulation via Random Draws and Plotting Discrete Distributions

Vectors of random draws can be generated for the Chapman, Petersen, or Bailey estimators using rChapman(), rPetersen(), or rBailey() from an assumed value of total abundance (N). This may be useful for simulation. Sample sizes n1 and n2 may be specified, but capture probabilities p1 and/or p2 may be used instead. If so, sample size(s) will first be drawn from a binomial distribution for each random draw before drawing from the abundance estimator. This will result in a greater degree of dispersion, but may be appropriate in some cases.

A plotting function is also provided for vectors of discrete (non-continuous) data. Abundance estimates are calculated from count data, with the result having a non-integer but also non-continuous support. It is possible that plotdiscdensity() may be more appropriate for plotting a discrete (non-continuous) density than a kernel density plot or histogram, as the discontinuity is made explicit.

draws <- rChapman(length=10000, N=1500, n1=100, n2=120)
plotdiscdensity(draws)

Hypothesis Testing and Power Calculation

Approximate p-values can be returned using pChapman(), pPetersen(), or pBailey(), which use many random draws to simulate a null distribution. The null hypothesis abundance is specified in the nullN argument, along with sample sizes n1 and n2. The observed abundance estimate can be specified using estN, or else the number of recaptures can be used directly, as m2. The alternative hypothesis can be specified using the alternative argument, as alternative="less", "greater", or "2-sided".

In the example given below, the null-hypothesis abundance is 500, with 100 individuals observed in the first and second events, with 28 recaptures in the second event.

output <- pChapman(nullN=500, n1=100, n2=100, m2=28, alternative="less")
output

## $estN
## [1] 350.7586
## 
## $pval
## [1] 0.02068

plotdiscdensity(rChapman(length=100000, N=500, n1=100, n2=100))   # null distribution
abline(v=500, lwd=2, lty=2)              # Null hypothesis abundance plotted as a dashed line
abline(v=output$estN, lwd=2, col=2)      # Observed (estimated) abundance plotted as a red line

Power calculation can be done with powChapman(), powPetersen(), or powBailey(), which use random draws from an assumed true (alternative) distribution, given the sample sizes of both events. The nullN argument specifies the abundance used by the null hypothesis, and the trueN argument specifies the assumed true abundance used for the power calculation. The n1 and n2 arguments give the sample sizes, and alternative gives the direction of the alternative hypothesis (defaults to "less"), with alpha specifying the level of the test to use.

In the example given below, the power is calculated for a one-tailed test of a null abundance of 500, assuming a true abundance of 400. The test powers are then calculated and plotted for assumed abundances from 250 to 450. If the true abundance is 325, a one-tailed test of \(H_0: N \geq 500\) will have a power of roughly 90%.

powChapman(nullN=500, trueN=400, n1=100, n2=100, nsim=1000)

## [1] 0.333

Ntotry <- seq(from=250, to=450, by=25)
power <- sapply(Ntotry, function(x)
  powChapman(nullN=500, trueN=x, n1=100, n2=100, nsim=1000))
plot(Ntotry, power)

Sample Size Recommendation

Given a best-guess at the true abundance and possibly the sample size in one sampling event, a recommendation for the sample size(s) can be calculated from the desired confidence and relative accuracy with the methods of Robson & Regier (1964) using n2RR(). Desired estimate confidence and accuracy (elsewhere termed “precision”) of 95% and 10%, respectively, is analogous to estimation “such that the estimated abundance will be within 10% of the true value 95% of the time”.

Recommendations for sample size are provided for two scenarios: if the size of the other sample is as specified, and if the two sample sizes are equal, which is the most efficient for sampling in terms of total sample size (\(n_1+n_2\)). The example below gives the full output for all allowed values of confidence and relative accuracy.

n2RR(N=1000, n1=100)

## $conf_0.99
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.50                  305        181
## 2    0.25                  550        270
## 3    0.20                  639        308
## 4    0.15                  746        364
## 5    0.10                  863        455
## 6    0.05                  961        622
## 7    0.01                  999        891
## 
## $conf_0.95
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.50                  180        136
## 2    0.25                  387        210
## 3    0.20                  484        244
## 4    0.15                  617        298
## 5    0.10                  780        385
## 6    0.05                  933        555
## 7    0.01                  998        862
## 
## $conf_0.9
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.50                  118        109
## 2    0.25                  293        177
## 3    0.20                  387        210
## 4    0.15                  525        260
## 5    0.10                  711        344
## 6    0.05                  908        511
## 7    0.01                  996        839
## 
## $conf_0.85
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.50                   81         90
## 2    0.25                  233        155
## 3    0.20                  320        186
## 4    0.15                  455        234
## 5    0.10                  652        314
## 6    0.05                  882        477
## 7    0.01                  995        820
## 
## $conf_0.8
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.50                   57         76
## 2    0.25                  189        139
## 3    0.20                  268        168
## 4    0.15                  395        213
## 5    0.10                  596        289
## 6    0.05                  856        448
## 7    0.01                  994        803
## 
## $conf_0.75
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.50                   41         65
## 2    0.25                  155        125
## 3    0.20                  225        153
## 4    0.15                  343        195
## 5    0.10                  543        267
## 6    0.05                  827        422
## 7    0.01                  992        785

Output can be simplified by providing values of confidence and relative accuracy in the conf and acc arguments.

n2RR(N=1000, n1=100, conf=c(0.9,0.95), acc=c(0.15,0.1,0.05))

## $conf_0.9
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.15                  525        260
## 2    0.10                  711        344
## 3    0.05                  908        511
## 
## $conf_0.95
##   Rel acc n2 from specified n1 n if n1=n2
## 1    0.15                  617        298
## 2    0.10                  780        385
## 3    0.05                  933        555

An alternative approach using simulation is provided with plotn2sim(), in which the relative accuracy is simulated for a range of sample size values for the second event, at different levels of confidence. The values to plot for n2 can be set by n2range and n2step, giving the range endpoints and step size, respectively.

plotn2sim(N=1000, n1=100)

Simulation is also used in plotn1n2simmatrix(), in which the relative accuracy is calculated for a range of sample size values for both sampling events for a given level of confidence, in this case 95%.

plotn1n2simmatrix(N=1000)

Tests of Consistency for the Estimator

Testing the Need for Partial Stratification

For a simple Petersen-type estimator of abundance to be consistent, one of three conditions must be true: either there must be complete mixing in the time between events, or the probability of capture must be equal for all individuals in the first event, or the probability of capture must be equal for all individuals in the second event. This is typically investigated by means of a set of \(\chi^2\) tests for each condition, with failure to reject the null hypothesis of any of the tests indicating that a simple Petersen-type estimator is reasonable to use. If the null hypotheses are rejected in all tests, a completely stratified or partially stratified (Darroch-type) estimator should be used.

The consistencytest() function is provided for the instance that sampling in both events is stratified in such a way that movement between strata may occur, such as spatial or temporal stratification. The stratification used in events 1 and 2 do not need to be the same. Arguments n1 and n2 accept vectors of sample sizes in each event by strata. Recaptures of marked individuals can be specified in one of two ways: either as vectors of strata membership in m2strata1 and m2strata2, or as a matrix in stratamat. Arguments m2strata1 and m2strata2 accept vectors of the first- and second-event stratum membership for each recaptured individual, with only entries of 1, 2, 3, ... accepted. Alternatively, stratamat accepts a matrix in which each cell represents the number of recaptures for each combination of event 1 and event 2 strata, with rows corresponding to event 1 strata and columns corresponding to event 2 strata.

In the example below, there were two strata in event 1 (with sample sizes of 284 and 199) and three strata in event 2, and 30 individuals were marked in stratum 1 of event 1 and were recaptured in stratum 1 of event 2. These tests yield strong evidence of incomplete mixing in the time between events and of unequal marked to unmarked ratios among recapture strata, but also no evidence against the probability of re-sighting a released animal being independent of its stratum of origin. Therefore, use of a simple Petersen-type estimator may be considered justified if closure is assumed.

mat <- matrix(c(30,15,1,0,22,15), nrow=2, ncol=3, byrow=TRUE)
consistencytest(n1=c(284,199), n2=c(347,3616,1489), stratamat=mat)

## MIXING TEST 
## H0: Movement probabilities from stratum i to stratum j are the same among sections (all theta_ij = theta_j) 
##  
##    1  2  3 not recaptured
## 1 30 15  1            238
## 2  0 22 15            162
## 
##  X-squared:  44.43179   df:  3   p-val:  1.221827e-09 
##  
## EQUAL PROPORTIONS TEST (SPAS terminology) 
## H0: Marked to unmarked ratio among recapture strata is constant (Sum_i a_i*theta_ij/U_j = k) 
##  
##              1    2    3
## recaptured  30   37   16
## unmarked   317 3579 1473
## 
##  X-squared:  125.4408   df:  2   p-val:  5.766008e-28 
##  
## COMPLETE MIXING TEST (SPAS terminology) 
## H0: Probability of re-sighting a released animal is independent of its stratum of origin (Sum_i theta_ij*p_j = d) 
##  
##                  1   2
## recaptured      46  37
## not recaptured 238 162
## 
##  X-squared:  0.3185941   df:  1   p-val:  0.5724538 
##

Since the decision to use a simple Petersen-type estimator is based on a failure to reject, the power of the tests used must be considered. The powconsistencytest() function uses simulation as well as Cohen’s non-central \(\chi^2\) method to calculate the power of the consistency tests, given anticipated values of sample sizes by strata and assumed movement probabilities.

The pmat argument is a matrix of assumed movement probabilities between strata, with rows corresponding to first-event strata and columns corresponding to second-event strata, with an additional column corresponding to the probability of NOT being recaptured in the second event. Values are conditioned on row, that is, by first-event strata. For example, a pmat with a first row equal to (0.05, 0.1, 0.15, 0.7) would imply a 5% chance that individuals captured in the first-event strata 1 will be recaptured in second-event strata 1, and a 70% chance that individuals captured in the first-event strata 1 will not be recaptured in event 2.

Because of the row-wise scaling, specifying a row equal to (0.05,0.1, 0.15, 0.7) would be equivalent to that row having values (1, 2, 3, 14).

mat <- matrix(c(1,2,3,10,3,2,1,10), nrow=2, ncol=4, byrow=TRUE)
powconsistencytest(n1=c(100,200), n2=c(100,100,100), pmat=mat)

## MIXING TEST 
## H0: Movement probabilities from stratum i to stratum j are the same among sections (all theta_ij = theta_j) 
## Sample size (first event):  300    Significance level:  0.05 
## 
##  Null hypothesis cross-classification probabilities:  
##        1      2      3 Not recap
## 1 0.0486 0.0417 0.0347    0.2083
## 2 0.0972 0.0833 0.0694    0.4167
## 
##  Alternative hypothesis cross-classification probabilities:  
##        1      2      3 Not recap
## 1 0.0208 0.0417 0.0625    0.2083
## 2 0.1250 0.0833 0.0417    0.4167
## 
##  Null hypothesis expected counts:  
##       1    2     3 Not recap
## 1 14.58 12.5 10.42      62.5
## 2 29.17 25.0 20.83     125.0
## 
##  Alternative hypothesis expected counts:  
##       1    2     3 Not recap
## 1  6.25 12.5 18.75      62.5
## 2 37.50 25.0 12.50     125.0
## 
## Power:  0.9496757    Power (from simulation):  0.9541 
##  
##  
## EQUAL PROPORTIONS TEST 
## H0: Equal probability of capture among n1 strata (all p_i equal) 
## Sample size (second event):  300    Significance level:  0.05 
## 
##  Null hypothesis capture probabilities:  
##                 1      2      3
## Recaptured 0.1250 0.1250 0.1250
## Unmarked   0.2083 0.2083 0.2083
## 
##  Alternative hypothesis capture probabilities:  
##                 1      2      3
## Recaptured 0.1458 0.1250 0.1042
## Unmarked   0.1875 0.2083 0.2292
## 
##  Null hypothesis expected counts:  
##               1    2    3
## Recaptured 37.5 37.5 37.5
## Unmarked   62.5 62.5 62.5
## 
##  Alternative hypothesis expected counts:  
##                1    2     3
## Recaptured 43.75 37.5 31.25
## Unmarked   56.25 62.5 68.75
## 
## Power:  0.3532643    Power (from simulation):  0.3515 
##  
##  
## COMPLETE MIXING TEST 
## H0: Equal probability of recapture among n2 strata (all p_j equal) 
## Sample size (first event):  300    Significance level:  0.05 
## 
##  Null hypothesis capture probabilities:  
##               1      2
## Marked   0.1250 0.2500
## Unmarked 0.2083 0.4167
## 
##  Alternative hypothesis recapture probabilities:  
##               1      2
## Marked   0.1250 0.2500
## Unmarked 0.2083 0.4167
## 
##  Null hypothesis expected counts:  
##             1   2
## Marked   37.5  75
## Unmarked 62.5 125
## 
##  Alternative hypothesis expected counts:  
##             1   2
## Marked   37.5  75
## Unmarked 62.5 125
## 
## Power:  0.05    Power (from simulation):  0.0526 
##  
##

Testing the Need for Complete Stratification

In the event that there is no movement observed between strata, or the population is stratified in such a way that there can be no movement between strata (such as by age or sex), a stratified estimator may be needed. The strattest() function provides a means of conducting the typical \(\chi^2\) tests to investigate the necessity of a stratified estimator.

Since strata membership remains the same in both events, usage is simpler, with n1, n2, and m2 accepting vectors of first- and second-event sample sizes and recapture numbers by strata. In the example below, there is strong evidence of unequal capture probabilities in the first event but not in the second, and very strong evidence of differences in respective stratum capture probabilities between the first and second events.

strattest(n1=c(100,100), n2=c(50,200), m2=c(20,15))

## 
##  Equal proportions test 
## H0: Equal probability of capture among n1 strata (all p_i equal) 
##  
##             1   2
## recaptured 20  15
## unmarked   30 185
## 
##  X-squared:  32.44394   df:  1   p-val:  1.226811e-08 
##  
## 
##  Complete mixing test 
## H0: Equal probability of recapture among n2 strata (all p_j equal) 
##  
##                 1  2
## recaptured     20 15
## not recaptured 80 85
## 
##  X-squared:  0.5541126   df:  1   p-val:  0.4566422 
##  
## 
##  First vs. second event test 
## H0: Probabilities of capture are the same between first and second events (all p_i equal to respective p_j) 
##  
##                1   2
## first event  100 100
## second event  50 200
## 
##  X-squared:  43.66013   df:  1   p-val:  3.906508e-11 
##

Similarly, the powstrattest() function provides power estimates for the tests used in strattest(), given assumed values of total abundance by strata and anticipated values of sample sizes by strata, using simulation as well as Cohen’s non-central \(\chi^2\) method.

powstrattest(N=c(1000,2000), n1=c(100,200), n2=c(200,200))

## $Test1
## First-event capture probabilities 
##  
## Detecting capture probabilities for each strata:  0.1 0.1 
## As compared to null values:  0.1 0.1 
## sample size:  400  and significance level:  0.05 
##  
## power:  0.05 
## power (from simulation):  0.04948 
## 
## $Test2
## Second-event capture probabilities 
##  
## Detecting capture probabilities for each strata:  0.2 0.1 
## As compared to null values:  0.1333333 0.1333333 
## sample size:  300  and significance level:  0.05 
##  
## power:  0.6707468 
## power (from simulation):  0.65945 
## 
## $Test3
## First-event vs. second-event capture probabilities 
##  
## Detecting capture probabilities for each strata:  0.6666667 0.5 
## As compared to null values:  0.5714286 0.5714286 
## sample size:  700  and significance level:  0.05 
##  
## power:  0.9928497 
## power (from simulation):  0.99563

Partially and Completely Stratified Estimators

If sampling is stratified in both sampling events and some movement exists between strata but mixing is incomplete in the time between events and tests indicate unequal capture probabilities in both events, a partially stratified (Darroch-type) estimator should be used. The NDarroch() function provides a means for doing this, and argument usage is the same as that in consistencytest().

The function returns abundance estimates and standard errors for each stratum in both sampling events, as well as an overall abundance estimate and standard error. Implementation is currently using Darroch’s method for calculations, and will only accept non-singular input matrices.

mat <- matrix(c(59,30,1,45,280,38,0,42,25), nrow=3, ncol=3, byrow=TRUE)
NDarroch(n1=c(484,1468,399), n2=c(847,6616,2489), stratamat=mat)

## $Nhat_strata1
## [1]  3584.554 16769.123 33699.792
## 
## $se_Nhat_strata1
## [1] 1914.125 6902.991 8979.480
## 
## $Nhat_strata2
## [1]  6360.409 18069.681 29623.380
## 
## $se_Nhat_strata2
## [1]  848.4842 4331.6190 7791.7721
## 
## $Nhat
## [1] 54053.47
## 
## $se_Nhat
## [1] 4890.999

If no movement can occur between strata, a completely stratified estimator can be used. Functions Nstrat(), vstrat(), sestrat(), and cistrat(), are provided to compute abundance estimates, estimated variance, standard error, and confidence intervals. In all cases, n1, n2, and m2 accept vectors of first- and second-event sample sizes and recapture numbers by strata, and the additional estimator accepts the type of Petersen-type estimator to use: either "Chapman", "Petersen", or "Bailey". Additionally, rstrat() generates random draws from given vectors of N, n1, and n2, or with vectors of capture probabilities p1 and/or p2.

Nstrat(n1=c(100,200), n2=c(100,500), m2=c(10,10))

## [1] 10080

vstrat(n1=c(100,200), n2=c(100,500), m2=c(10,10))

## [1] 6513699

sestrat(n1=c(100,200), n2=c(100,500), m2=c(10,10))

## [1] 2552.195

cistrat(n1=c(100,200), n2=c(100,500), m2=c(10,10))

## $Nstrat
## [1] 10080
## 
## $Nhat_by_strat
## [1]  926.3636 9153.6364
## 
## $ciNorm_strat
## [1]  5077.79 15082.21
## 
## $ciBoot_strat
## [1]  6649.363 20866.843

draws <- rstrat(length=10000, N=c(500,1000), n1=c(100,200), n2=c(100,500))
plotdiscdensity(draws)

draws <- rstrat(length=100000, N=c(5000,10000), n1=c(500,200), n2=c(500,200))
plotdiscdensity(draws)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.