// library(ascii)
// setwd("c://aaaWork//web//fishR//bookex//AIFFD//Box6_8")
// Sweave("Box6_8a.Rnw",driver=RweaveAsciidoc)








AIFFD Box 6-8  Vignette
=======================
:Author:    Derek H. Ogle
:Email:     dogle@northland.edu
:Date:      30-June-2009
:Revision:  1

The raw data for the example in Box 6.8 does not exist.  Indeed enough information is given in the box to complete the calculations but most of that information in in summarized form.  To illustrate how to achieve some of the results from raw data I simulated the length data from the results shown for the length frequency graph (the simulation process is illustrated link:Box6_8Aa.html[here]).  The age data that would be required to estimate L~infinity~ and K cannot be simulated from the given information.  Thus, the process and the results shown in Box 6.8 will not be perfectly reproduced here.  The results should be a close facsimile.

== Required Packages and Setting Options
----
> library(NCStats)  # Subset

----

== Preparing Data
The working directory is set, the link:box6_8.txt[] data file is read, and the structure of the data frame is observed with,
----
> setwd("c://aaaWork//web//fishR//bookex//AIFFD//Box6_8")
> d <- read.table("box6_8.txt", header = TRUE)
> str(d)
'data.frame':	1559 obs. of  1 variable:
 $ len: num  90.9 90.9 90.4 90.5 90.7 ... 
----

Note that the authors restricted their analysis to fish 150-mm and larger because those were full recruited to the gear.  A new data frame with just these fish is constructed with,
----
> d1 <- Subset(d, len >= 150)
> str(d1)
'data.frame':	1267 obs. of  1 variable:
 $ len: num  150 150 150 150 150 ... 
----
This data frame should now be attached for use,
----
> attach(d1)
	The following object(s) are masked from d1 ( position 3 ) :

	 len 


	The following object(s) are masked from d1 ( position 4 ) :

	 len  
----

== Beverton-Holt Method for Estimating Z from Mean Length
The authors illustrate the use of equation 6.8 to estimate Z.  This method requires knowledge of L~infinity~ and K from the von Bertalanffy growth model (see link:../Box5_5/Box5_4a.html[Box 5.4 Vignette] for how to fit a von Bertalanffy model to length-at-age data).  In this box we are given this summary information,
----
> Linf <- 636
> K <- 0.226

----

In addition, the mean length and the minimum length (that is fully vulnerable to the gear) must be found from the raw data,
----
> Lmean <- mean(len)
> Lmean
[1] 255.6757 
> Lmin <- min(len)
> Lmin
[1] 150 
----

Thus, using equation 6.8 the instantaneous mortality (Z) is estimated to be,
----
> Z1 <- K * (Linf - Lmean)/(Lmean - Lmin)
> Z1
[1] 0.8133687 
----

== Hoenig Method for Estimating Z from Median Length
Hoenig's (1983) modification of the Beverton and Holt's method is shown in equation 6.9 in the book.  This method still requires L~infinity~, K, and Lminimum but it also requires the median length.  The median length is found with,
----
> Lmedian <- median(len)
> Lmedian
[1] 250 
----
and the estimate of Z using equation 6.9 is,
----
> Ymedian <- -log(1 - Lmedian/Linf)
> Ymin <- -log(1 - Lmin/Linf)
> Z2 <- 0.693 * K/(Ymedian - Ymin)
> Z2
[1] 0.6798504 
----

== Estimating Z from Length-Converted Catch Curve
Estimating Z from a length-converted catch curve as described by Pauly (1984) requires finding the midpoint of length bins and the number of fish in each of those bins.  In essence, one needs to know the basic calcualtions for constructing a length-frequency histogram.  In R, histograms are constructed with +*[red]#hist()#*+.  This function only __requires__ a vector of values to construct the histogram as the first argument.  However, you can control the bins used by the histogram algorithm by sending break points for the bins to the +*[red]#breaks=#*+ argument.  In addition, +*[red]#hist()#*+ defaults to the bins being right-inclusive -- i.e., a value of 10 would be included in the 9-10 bin rather than the 10-11 bin -- whereas most fisheries biologists are used to using left-inclusive bins.  Left-inclusive bins can be used by including the +*[red]#right=FALSE#*+ argument to +*[red]#hist()#*+.  In this example, I want to use left-inclusive bins that start at 150-mm, end at 530-mm, and are 10-mm wide.  These bins are constructed with,
----
> bins <- seq(150, 530, by = 10)

----

The results of +*[red]#hist()#*+ can be saved to an object so that the midpoints and the frequencies in each interval can be accessed for Pauly's method.  This histogram is constructed with,
----
> h <- hist(len, breaks = bins, right = FALSE, xlab = "Total Length (mm)")

----
image::Box6_8a-012.png[]

The midpoint values and frequencies are then extracted with,
----
> h$mids
 [1] 155 165 175 185 195 205 215 225 235 245 255 265 275 285 295 305 315 325 335 345 355 365 375 385 395 405 415 425 435 445 455 465 475 485 495 505
[37] 515 525 
> h$counts
 [1] 106  99  82  81  56  45  27  36  43  52  65  73  63  42  44  40  37  31  36  15  19  18  13   8  20  19   9  12   9  11  13  13   4  10   6   3
[37]   3   4 
----

The midpoint values are converted to t-prime values with,
----
> tprime <- -log(1-h$mids/Linf)
> round(tprime,3)                # compare to values in Box 6.8
 [1] 0.279 0.300 0.322 0.344 0.366 0.389 0.413 0.437 0.461 0.486 0.512 0.539 0.566 0.594 0.623 0.653 0.684 0.715 0.748 0.782 0.817 0.853 0.891 0.930
[25] 0.970 1.013 1.057 1.103 1.152 1.203 1.257 1.314 1.374 1.438 1.506 1.580 1.659 1.746 
----
and the counts are log-transformed with,
----
> logN <- log(h$counts)

----

The linear regression denoted by equation 6.10 in the book is then fit with,
----
> lm1 <- lm(logN ~ tprime)

----
and the slope coefficient (i.e., second coefficient) is extracted with,
----
> coef(lm1)
(Intercept)      tprime 
   4.973924   -2.228366  
> b <- coef(lm1)[2]
> b
   tprime 
-2.228366  
----

Which results in the following estimate of Z,
----
> Z3 <- K * (1 - b)
> Z3
   tprime 
0.7296106  
----

== Empirical Estimates of L~infinity~ and K
The authors (on page 249) provide three methods for estimating L~infinity~ and K if length-at-age data do not exist.  These methods are illustrated here.

=== Three Largest Fish Method Estimate of L~infinity~
Finding the mean lengths of the three largest fish in the length vector is a bit of work which is best shown in several steps -- (1) order the vector positions from smallest to largest fish, (2) find positions of three largest fish, (3) find lengths of those three fish, and (4) find mean of those lengths.  This is illustrated below,
----
> ord.len <- order(len)
> ord.top3 <- ord.len[(length(len) - 2):length(len)]
> len.top3 <- len[ord.top3]
> len.top3
[1] 520.1 520.4 520.7 
> Lmean.top3 <- mean(len.top3)
> Lmean.top3
[1] 520.4 
----

Pauly (1984) then suggests dividing this mean by 0.95 to produce an estimate of L~infinity~,
----
> Linf1 <- Lmean.top3/0.95
> Linf1
[1] 547.7895 
----
This is dramatically lower then the L~infinity~ provided by the authors.

=== Froese and Binohlan (2000) Emprical Equation Method Estimate of L~infinity~
Froese and Binohlan (2000) developed an empirical equation for estimating L~infinity~ from the maximum length observed for the population.  The equation is log(Linf) = 0.044 + 0.984log(Lmax).  The maximum value in these data and the estimate of L~infinity~ using this equation are,
----
> Lmax <- max(len)
> Lmax
[1] 520.7 
> Linf2 <- exp(0.044 + 0.984 * log(Lmax))
> Linf2
[1] 492.3015 
----
This is even more dramatically lower then the L~infinity~ provided by the authors.

=== Wetherall (1987) Method Estimate of L~infinity~
The method of Wetherall (1987) requires knowing the lower limit of each length interval,
----
> Lx <- h$breaks
> Lx
 [1] 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500
[37] 510 520 530 
----
In addition, the method requires knowing the mean lengths of all fish in each interval and all larger intervals.  For example, we must find the mean length of fish in the 150-mm and all larger intervals, in the 160-mm and all larger intervals, the 170-mm and all larger intervals, and so on.  This type of calculation is best handled with a loop as follows,
----
> Lmeans <- NULL                                 # initiate the Lmeans vector
> for (i in 1:length(Lx)) {                      # loop through length intervals
+   Lmeans <- c(Lmeans,mean(len[len>=Lx[i]]))
+ }
> Lmeans
 [1] 255.6757 265.2842 275.0573 283.8116 293.1231 299.9441 305.5528 308.8864 313.2161 318.3559 324.6812 333.0746 343.6345 354.1260 361.9154 370.8235
[17] 379.8128 389.1072 397.7955 409.3823 414.7160 421.6949 428.7274 433.9826 437.1434 445.2034 454.0010 458.4500 464.4684 469.0507 474.6804 482.0465
[33] 491.4133 494.6346 503.4625 511.3400 515.9714 520.3000      NaN 
----

The difference between these means and the lower limit of the length interval must then be found,
----
> Ldiff <- Lmeans - Lx

----

and the linear regression of these differences on the lower limits of the length intervals is then fit with,
----
> lm2 <- lm(Ldiff ~ Lx)

----
The coefficients of this model fit are extracted with,
----
> coef(lm2)
(Intercept)          Lx 
152.4397441  -0.2763679  
> a <- coef(lm2)[1]
> b <- coef(lm2)[2]

----
and L~infinity~ is estimated with,
----
> Linf3 <- -a/b
> Linf3
(Intercept) 
   551.5827  
----
Again, this is dramatically lower then the L~infinity~ provided by the authors.