> library(NCStats) # bin.ci, rcumsum
> library(NCStats) # bin.ci, rcumsum
The data for this box is very simple — number of fish in seven age-classes. Because of this simplicity the data can simply be entered into R using c() (for concatenation). The data could have been entered into an external file and read into R (I demonstrate this in the Box6.3 vignette). The data for this box are entered as follows:
> age <- 1:7
> catch <- c(257,407,147,32,17,5,4)
> rbind(age,catch) # for viewing purposes only.
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
age 1 2 3 4 5 6 7
catch 257 407 147 32 17 5 4
The author’s of the box limited the first calculation to only those fish that were age-2 or older. With this constraint the annual mortality calculation requires the total number of fish that were age-2 and older. In this example, the catch of fish age-2 and older can be found in positions two through seven of the catch vector. Certain positions of a vector can be obtained by appending the desired positions within square brackets immediately after the name of the vector. This is illustrated with the following examples,
> catch[2] # 2nd position [1] 407 > catch[c(2,4)] # 2nd & 4th positions [1] 407 32 > catch[2:4] # 2nd through 4th positions [1] 407 147 32 > catch[-2] # all but the 2nd position [1] 257 147 32 17 5 4
The total catch of age-2 and older fish is thus
> N <- sum(catch[2:7]) # sum in positions 2 thru 7 > N [1] 612
The annual mortality rate is then the number of age-2 fish divided by this total (i.e., the proportion of all fish age-2 and older that were age-2). Thus,
> A <- 100*catch[2]/N # annual mortality estimate > A [1] 66.50327 > se.A <- sqrt(A*(100-A)/N) # SE of A > se.A [1] 1.907862
The Heincke method essentially computes the proportion of a certain age out of a total number in an age range. This is nothing more then computing a binomial proportion — in fact the SE computed above is based on binomial theory. A confidence interval for a binary proportion can be computed with bin.ci() from the NCStats package. This function requires the number of "successes" (i.e., the age in question) as the first argument and the total number of trials (i.e., all fish in the age categories) as the second argument. Thus, a (default) 95% CI for the annual mortality rate computed above is computed with,
> 100 * bin.ci(catch[2], N) 95% LCI 95% UCI 62.67124 70.12941
The authors of Box 6.2 mention that the estimated annual mortality rate, using the method above, for age-3 and older fish is 72%. This begs the question of whether there is an efficient method of computing the estimated annual mortality rate for each age (and older). Fortunately, this can be accomlished fairly efficently with rcumsum() from the NCStats package. This function computes the reverse cumulative summation for a vector. In other words, it will compute the total number of fish age-1 and older, then the total number of fish age-2 and older, then age-3 and older, etc. For example,
> N <- rcumsum(catch) > N [1] 869 612 205 58 26 9 4
The estimated annual mortality rates can then be calculated by dividing the catch vector by these reverse cumulative sums,
> 100 * catch/N [1] 29.57422 66.50327 71.70732 55.17241 65.38462 55.55556 100.00000
Appropriate CIs are computed with,
> 100 * bin.ci(catch, N) 95% LCI 95% UCI 26.63515 32.69308 62.67124 70.12941 65.18616 77.42990 42.45208 67.25015 46.22057 80.58777 26.66513 81.12215 51.01092 100.00000
Finally, if you want to put all of this together into a nice display then try,
> data.frame(age, catch, N, A = 100 * catch/N, 100 * bin.ci(catch, N)) age catch N A X95..LCI X95..UCI 1 1 257 869 29.57422 26.63515 32.69308 2 2 407 612 66.50327 62.67124 70.12941 3 3 147 205 71.70732 65.18616 77.42990 4 4 32 58 55.17241 42.45208 67.25015 5 5 17 26 65.38462 46.22057 80.58777 6 6 5 9 55.55556 26.66513 81.12215 7 7 4 4 100.00000 51.01092 100.00000