Recently I taught a class on estimating fisheries exploitation rates with tagging studies. A common study design to estimate exploitation rates is to release tagged fish into the system and allow anglers to return the tags when they harvest a fish. The simplest estimate of exploitation rate is obtained by dividing the number of tags returned by the number of tagged fish released into the fishery. However, three key sources of uncertainty are mortality of fish due to the tagging process, tag loss, and non-reporting by anglers. Studies can be designed to estimate these quantities by combining both passive and active tags (e.g. Kerns et al. 2013) or performing additional side studies. Often tag loss, tag reporting, and non-reporting rates are adopted from other studies and systems and used to correct exploitation rates for these unknowns. The following method is a parametric bootstrap that incorporates external estimates of tag loss, tag mortality and non-reporting rates into the exploitation rate estimate. Furthermore, it incorporates the uncertainty in these quantities using the bootstrap procedure.

The process that generates the data (i.e. the number of tags returned by anglers) can be thought of as a series of nested stochastic events. After the tagged fish are released, there is a binomial process that determines how many tags fall out, followed by a binomial process that determines how many fish die, followed by a binomial process that determines how many fish get harvested, followed by a binomial process that determines how many tags get reported. If we assume that tag loss and tagging mortality occur simultaneously before harvest, the model can be represented as:

where *T _{mort}* is the tag-related mortality rate,

*T*is the tag loss rate,

_{loss}*U*is the exploitation rate, and

*R*is the tag reporting rate. However, because

*U*is unknown and is the value that we are interested in estimating, the model must be rearranged for the bootstrap procedure. First we must calculate the

*U*point estimate as:

Then the bootstrap simulation is represented as:

where *U _{boot}* is the values of

*U*for one bootstrap iteration. Further uncertainty in the

*T*,

_{mort}*T*, and

_{loss}*R*rates can be incorporated into the bootstrap by drawing a random number from a beta distribution for each rate for each iteration. The following R code demonstrates this procedure.

#------------------------------------------------------------------------------# # This part builds a function to calculate the shape parameters of the # beta distribution from a specified mean and standard deviation. #------------------------------------------------------------------------------# betapar = function(mu,sd){ shape1 = ((1-mu)/(sd^2) - 1/mu)*mu^2 shape2 = shape1*(1/mu - 1) c(shape1=shape1,shape2=shape2) } #------------------------------------------------------------------------------# # This part specifies the input parameters for the simulation and tests # the betapar function #------------------------------------------------------------------------------# # Input values tagloss.mu = 0.07 tagloss.sd = 0.05 tagmort.mu = 0.12 tagmort.sd = 0.07 reprate.mu = 0.66 reprate.sd = 0.18 # test function (the first values is shape1(alpha) and the second is shape2 (beta)) betapar(tagloss.mu,tagloss.sd) mean(rbeta(10000,betapar(tagloss.mu,tagloss.sd)[1],betapar(tagloss.mu,tagloss.sd)[2])) sd(rbeta(10000,betapar(tagloss.mu,tagloss.sd)[1],betapar(tagloss.mu,tagloss.sd)[2])) # Mean and standard deviation match input values pretty well # Data tagged = 200 returned = 50 #------------------------------------------------------------------------------# # This part calculates estimates #------------------------------------------------------------------------------# fun = function(tagged, returned){ Adj_taggs = tagged * (1-tagloss.mu) * (1-tagmort.mu) True_catch = returned/reprate.mu U_point_est = True_catch/Adj_taggs #----------------------------------- Tl = rbeta(1,betapar(tagloss.mu,tagloss.sd)[1],betapar(tagloss.mu,tagloss.sd)[2]) Tm = rbeta(1,betapar(tagmort.mu,tagmort.sd)[1],betapar(tagmort.mu,tagmort.sd)[2]) Rr = rbeta(1,betapar(reprate.mu,reprate.sd)[1],betapar(reprate.mu,reprate.sd)[2]) Tadj = rbinom(1,tagged, (1-Tl)*(1-Tm)) Ctrue = rbinom(1,Tadj,U_point_est) Ret = rbinom(1,Ctrue,Rr) Uboot = Ret/Rr/Tadj Uboot} #------------------------------------------------------------------------------# # This part runs bootstrap simulation for a range of sample sizes #------------------------------------------------------------------------------# tvec = c(20,200,2000,20000) rvec = c(5,50,500,5000) results = matrix(NA,10000,4) for(i in 1:4){ results[,i] = replicate(10000,fun(tvec[i],rvec[i]),simplify=T) } boxplot(results,ylim=c(0,1))