---
title: "Catch Curves (with selectivity) Installation"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    theme: cerulean
  pdf_document: default
  word_document: default
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Catch curves (with selectivity)

Catch curve analysis is a method for estimating the total mortality of a stock (<em>Z</em>). The rate at which individuals die can be estimated from the slope ofa regression of the logarithms of relative numbers present in each age class on age. It can be used whenever there is one or more years of catch-at-age data (or at-length data, if they can be converted to age). The data can be fishery-dependent or -independent so long as they are representative of the population's relative age/length structure. Fishing mortality (<em>F</em>) can be estimated as <em>Z</em>-<em>M</em> given an estimate of natural mortality (<em>M</em>) from another source (e.g., the literature, from a marine protected area, or from a tagging study).

The method is based on the equilibrium assumption (i.e., the total population receives constant recruitment and there is a constant mortality rate each year), which implies that all cohorts would be identical and the numbers in each age class would decline exponentially. Thus, if one assumes that fishing mortality and natural mortality have been constant, and one has the age-composition of the catch for a single year, if one plots the log-transformed numbers-at-age  against age, then the slope of a linear regression through the data will provide an estimate of total mortality (<em>Z</em>).

Conventional catch-curve methods rely on the strong assumptions of constant fishing and natural mortality rates above some fully selected age that is usually estimated by visually inspecting a plot of log-catch at age vs age. Rather than just ignoring the earlier weakly selected ages below the (assumed) fully selected age, as in the classical catch-curve, it is possible to estimate selectivity parameters as well as the total mortality rate. However, this requires the use of an equilibrium based age-structured model rather than the log-transformed linear regression (see equations in Wayte and Klaer, 2010).

Catch curves with selectivity use a simple age-structured model to estimate selectivity from the age data. This provides an estimate of fully-selected fishing mortality rather than an average fishing mortality applied to all included age classes (as in the classic catch curve). This also circumvents the need to select a minimum age from which to fit the regression (to avoid the effects of selectivity), as is required for the classical catch curve.

As opposed to classic catch curves, catch curves with selectivity provide an estimate of fully-selected fishing mortality, rather than an average fishing mortality applied to all included age classes. 

Although there are several packages that can be used to run catch curves with selectivity, the one by Dr Malcolm Haddon is supported here, because his version has been applied in several cases in Australia (refs needed).

See the Appendix for further description and equations, as supplied by Dr. Haddon and referencing Wayte and Klaer (2010) <a href="https://www.sciencedirect.com/science/article/pii/S0165783610002055?via%3Dihub"><em>(Fisheries Research 106: 310-320)</em></a>.

# Installation
This package is run from R as found on CRAN. R can be installed from https://cran.r-project.org/. 

The recommended way to run R is using RStudio. After you have installed R, then install RStudio. It can be found https://www.rstudio.com/products/rstudio/download/.

The supported version of Catch Curve Analysis is part of a set of stock assessment tools within Dr Haddon's <em>datalowSA</em> package. See the link https://github.com/haddonm/datalowSA. 

Therefore one needs to first install <em>  <a href="https://rdrr.io/github/haddonm/datalowSA/">datalowSA</em></a> from within R. 


Set up checking/installing packages functions

```{r checking and installing functions, eval=TRUE}
# Ordinary checking/installing function
Check_Install.packages <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg,dependencies = TRUE) else print(paste0("'",pkg,"' has been installed already!"))
  sapply(pkg, require, character.only = TRUE)
}

# github checking/installing function
Github_CheckInstall.packages <- function(pkg,github_pkg,build_vignettes){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    devtools::install_github(pkg,github_pkg,build_vignettes)
  sapply(pkg, require, character.only = TRUE)
}

```

Installing packages from GitHub directly. 

Check and install "devtools"

```{r install_from_github, eval=TRUE}

Check_Install.packages("devtools")

```

Check and install <em>datalowSA</em> from github

``` {r packages install_from_github, eval=TRUE}
# Check if existed and if not, install the required package TMBhelper:
github_pkg = "haddonm/datalowSA"
Github_CheckInstall.packages("datalowSA",github_pkg=github_pkg,build_vignettes=TRUE)

#install_github("haddonm/datalowSA")

```

<p style="color:red"> NB Once you have installed <em>datalowSA</em>, you don't need to run the install.packages code chunk again. Luckily it will just give an error stating it won't reinstall if your version is up to date.

Note that you may get an error, because the current version you have of Rcpp may preclude the installation of <em>datalowSA</em>. If this occurs, try reinstalling Rcpp from CRAN (Tools -> Install Packages). </p>


```{r load package, eval=TRUE}
library(datalowSA)
```

It should be noted that there are also good vignettes for this package by typing <em>browseVignettes(package = "datalowSA")</em>. 

## Reading in the data

To run the Catch Curve Analysis with Selectivity requires the following inputs:

* **Numbers-at-age** * - for a single year, as a vector. This must be absolute numbers (as opposed to proportions-at-age) and a representative (across time) random sample.
* **A single natural mortality estimate** * 
* **The maximum age (within the available data)** * 
* **Age-at-50%-selectivity, <em>a</em><sub>50</sub>** * 
* **Delta (age-at-95% selectivity (<em>a</em><sub>95</sub>) less the age-at-50% selectivity (<em>a</em><sub>50</sub>), ie (<em>a</em><sub>95</sub>-<em>a</em><sub>50</sub>)** *  
* **An estimate of current fishing mortality (used as an initial value)** * 

(Note that delta feeds into the equation for selectivity-at-age: 

Selectivity at age is modelled using (1+exp(-ln(19)(a - <em>a</em><sub>50</sub>)/delta))<sup>-1</sup> instead of (1+exp(-ln(19)(a -<em>a</em><sub>50</sub>)/(<em>a</em><sub>95</sub>-<em>a</em><sub>50</sub>)))<sup>-1</sup>
to reduce the correlation between two estimated parameters (<em>a</em><sub>50</sub> and <em>a</em><sub>95</sub>), and to avoid <em>a</em><sub>50</sub> values larger than <em>a</em><sub>95</sub> values. It is relatively easy to get a <em>a</em><sub>95</sub> value less than the <em>a</em><sub>50</sub> value, and this leads to a selectivity pattern that declines instead of increasing with age.


In general, it goes without saying that one requires a minimum number of ages in the data set to ensure that there are more data points than parameter values. The selectivity-included method in Dr. Haddon's package estimates three parameters and it is a reasonable rule of thumb to have at least three data points per parameter estimated - so, at least nine ages here.

A simulated example dataset has been created for the whole toolbox. Let us now read in these data. We assume you have saved the data to the same directory as your working directory.  Otherwise you can change our code to point to the directory where the data file is and call it your working directory using the setwd command.

The simulated catch-at-age data, <em>TotalCatchAge</em>, has years (1960-2001) in rows (2-43), with 25 age bins (from ages 1 to 25) in columns (B to Z). The first 11 years of this file have no data, so we will focus on year 1971 onwards.

Let us read in the data:

```{r load data }

TotCatAge <- read.csv("TotalCatchAge.csv")
print(TotCatAge[12:42,])
```

On inspection, there is no catch within many older ages from 1988. On the other hand, the earlier years also appear to have very few counts relative to those in later years.

We can visualise the age composition for specified years using the function <em>compyear</em>. In our case, we will focus on 1980-1995 (rows 21 to 36) inclusive. (Note that, when plotting, we will have to remove the year column from the data).

``` {r visualise age data}
naa <- TotCatAge
yrs <- 21:36
naa <- naa[yrs,2:26] # We must remove the first column containing the years
ages <- 1:25 
#plotprep() 
compyear(naa,1980:1995,ages,plots=c(4,4),freq=TRUE,border=3) 
```

As per the advice in Dr. Haddon's vignettes: when our simulated example age data are considered using these frequency plots, it is clear that the notion of a exponential decline is only approximated by these data. Such deviations from the theoretical expectation derive from the fact that real populations are not in equilibrium and recruitment is not constant. So the dynamics appear to be realistic although the assumptions about equilibrium do not appear to be valid. Despite this we can use these data to illustrate the use of catch curves.

The sample sizes from each year are printed on the above plot, but can also be found by using the following.

``` {r find sample sizes}
sample_sizes <- apply(TotCatAge[,2:26],1,sum)
sample_sizes
```

We can see a strong year cohort progressing through from 1990 onwards, which violates the assumption of constant recruitment. Despite the large sample sizes in the 1990s, let us instead focus our analysis on years 1986 and 1988 in the first instance. These are years with high sample sizes (2840 and 2974, respectively), but without visual evidence of strong recruitment pulses. 

Our dataset was created using Stock Synthesis based on an age-at-maturity of roughly 3-4 years and a rate of natural mortality of 0.2 yr <sup>-1</sup> for females, and 0.25 yr <sup>-1</sup> for males. 

Regarding the input parameters:

- Natural mortality estimate = 0.225 yr <sup>-1</sup> (average of males and females)
- The maximum age in the sample data = 25.
- Age-at-50%-selectivity = 4.45 years (average of males and females: 3.97 for females, 4.93 for males)
- Age-at-95%-selectivity = 19.07 years (average of males and females: 13.14 for females, 19.07 for males)
- Delta = 19.07-4.45 = 14.62 years
- An estimate of current fishing mortality (used as an initial value) = <p style="color:red">???STILL NEEDED </p>

### Fitting the model to our data  

<em>selectCC</em> is the function that sets up and runs the catch curve (with selectivity) model and simulations. 

The main input data expected is a vector of numbers-at-age for a single year. As stated above, we will try years 1986 and 1988 (rows 27 and 29 of <em>TotCatAge</em>), separately.

We will assign the same names to the numbers-at-age vector and the other inputs as per datalowSA.

In creating our catch-at-age vector, there are a few potential trip-points. The first is that the data (numbers-at-age, <strong>numaa</strong>) MUST be a vector, not a list. The second is that the rownames (or names) of the vector should equate to a character list of the actual ages.

```{r Create catch-at-age vector for 1986}
numaa <- TotCatAge[27,2:26]
names(numaa) <- as.character(1:25)
numaa <- unlist(numaa)
print(numaa)
```

The model needs the natural mortality and maximum sampled age as a globals list which we can call <em>glb</em>, copying the example in <em>datalowSA</em>. Let's create the list for <em>glb</em> (which we could of course have instead read in).

```{r CreateGlb}
NatM <- 0.225
maxage <- 25
glb <- list(M = NatM,maxage = maxage)
print(glb)
```

Finally, the parameters relating to selectivity:

```{r Create selectivity pars}
pars <- c(A50=4.45,delta=14.62,fcur=0.6)
print(pars)
```
 
## Running the model

Now that we have the data and inputs set up, the next is to run the model itself, using <em>selectCC</em>. 

```{r runModel}
library(datalowSA)
out <- selectCC(glb$M,glb$maxage,numaa,pars,plot=TRUE)
```

This shows the fitted catch curve.

We can examine the model output to check the fitted selectivity parameters (<em>age-at-50%-maturity</em>, <em>A50</em>, <em>delta</em>, and current fishing mortality, <em>f</em><sub>cur</sub>) and predicted catch-at-age:

``` {print output}
print(out)
```

We can see that, for 1986, this results in a fully selected fishing mortality estimate of 0.1342896 yr<sup>-1</sup> (and a logistically-increasing selectivity-at-age, as indicated in the resultant plot). 

Using instead the catch-at-age data from 1988:

```{r Create catch-at-age vector for 1992, run model and print output}
numaa <- TotCatAge[29,2:26]
names(numaa) <- as.character(1:25)
numaa <- unlist(numaa)
print(numaa)
out <- selectCC(glb$M,glb$maxage,numaa,pars,plot=TRUE)
print(out)
```

For 1988, fully selected fishing mortality is estimated to be 0.2494271 yr<sup>-1</sup> (again with a logistically-increasing selectivity-at-age, indicated in the resultant plot). 

(Note: using 1990's data (row 31) results in a fishing mortality estimate of 0.1949766 yr<sup>-1</sup> for all ages. The fitted proportion-at-age plot shows a negative exponential decline, which equates to all age classes being fully selected. This is the result of the strong cohort in the earliest age class. Using 1992's data (row 33), we see an increase, a peak and then a decline in proportion-at-age (reflecting the strong cohort having since aged by 2 years), which yielded a more logistic increase in selectivity with age. The fully-selected fishing mortality was estimated to be very high, at 0.952 yr<sup>-1</sup>, again as a result of this strong year cohort).

Other years can be tested, as desired. It is important to reiterate that the data do not result in a constant fishing mortality over time.

## Interpreting results in a management context 

Catch curve analysis results in an estimate of total mortality, and, assuming a fixed natural mortality, fishing mortality. It does not provide an estimate of stock status per se, except where consistently low fishing mortality is used as a proxy for sustainability (although this would require many years of age-composition data). The estimates of fishing mortality could be compared to a reference point such as the fishing mortality at maximum sustainable yield, if this were available.

As only a single year of numbers-at-age is required, catch curve analysis suits fisheries which have occasional age-structured data, and whose catch-at-age distribution is not expected to vary strongly through time. As such, catch-curve analysis is frequently applied to data-poor fisheries, where a variety of disparate information sources (each of low quality for stock assessment) may be available. Future work could formalize a process by which catch curves are combined with expert opinion (i.e., interviews with fishers) and other data-poor analyses (e.g., changes in mean length) into a quantitative estimate of stock status (Thorson and Prager 2011).

## Key references

### Catch Curve Methodology

Chapman, D.G. and D.S. Robson. 1960. The analysis of a catch curve. <em>Biometrics</em> 16:354-368.

Dunn, A., Francis, R.I.C.C. and I.J. Doonan. 2002. Comparison of the Chapman-Robson and regression estimators of Z from catch-curve data when non-sampling stochastic error is present. <em>Fisheries Research</em> 59:149-159. https://doi.org/10.1016/S0165-7836(01)00407-6

Gulland, J.A. 1971. The fish resources of the ocean. West Byfleet, UK: Fishing News Books. 

Smith, M.W., Then, A.Y., Wor, C., Ralph, G., Pollock, K.H. and J.M. Hoenig. 2012. Recommendations for catch-curve analysis. <em>North American Journal of Fisheries Management</em> 32:956-967. http://dx.doi.org/10.1080/02755947.2012.711270
 
Thorson, J.T. and M.H. Prager. 2011. Better catch curves: Incorporating age-specific natural mortality and logistic selectivity. <em>Transactions of the American Fisheries Society</em> 140:356-366. http://dx.doi.org/10.1080/00028487.2011.557016

Wayte, S.E. and N.L. Klaer. 2010. An effective harvest strategy using improved catch-curves. <em>Fisheries Research</em> 106:310-320. https://doi.org/10.1016/j.fishres.2010.08.012

### Catch Curve Evaluation

Allen, M.S. 1997. Effects of variable recruitment on catch-curve analysis for crappie populations. <em>North American Journal of Fisheries Management</em> 17(1):202-205. https://doi.org/10.1577/1548-8675(1997)017<0202:EOVROC>2.3.CO;2

Griffiths, S.P. 2010. Stock assessment and efficacy of size limits on longtail tuna (<em>Thunnus tonggol</em>) caught in Australian waters. <em>Fisheries Research</em> 102(3):248-257. https://doi.org/10.1016/j.fishres.2009.12.004 

Oyarzún, C., Cortés, N. and E. Leal. 2013. Age, growth and mortality of southern rays bream <em>Brama australis</em> (Bramidae) off the southeastern Pacific coast. <em>Revista de biología marina y oceanografía</em> 48(3). http://dx.doi.org/10.4067/S0718-19572013000300014 

See references in http://derekogle.com/fishR/examples/oldFishRVignettes/CatchCurve.pdf.

## Appendix: Catch Curve Equations

### Catch-Curves with Selectivity

(reproduced with permission from Dr. Haddon's vignettes)

A simple model of the expected explonential decline of numbers-at-age is needed:

$$\begin{matrix}
   {{N}_{0}}= & 1 & a=0  \\
   {{N}_{a}}= & {{N}_{a-1}}\exp \left( -{{s}_{a-1}}F-M \right) & 0<a<{{a}_{\max }}  \\
   {{N}_{a}}= & \frac{{{N}_{a-1}}\exp \left( -{{s}_{a-1}}F-M \right)}{\left( 1-\exp \left( -{{s}_{a}}F-M \right) \right)} & a={{a}_{\max }}  \\
\end{matrix} $$

where $N_a$ is the numbers at age $a$, $s_a$ is the selectivity at age $a$, $F$ is the fully- selected average fishing mortality over the years for which agedtaa are available, $M$ is the assumed rate of natural mortality, and $a_{max}$ is the maximum age. Selectivity is assumed to be asymptotic (so this method may not be suitable for data collected using gillnets or hook gear):

$$s_a = \frac{1}{1+exp(-\ln(19)(a-a_{50})/\delta)}$$
where $a$ is age, $a_{50}$ is the age-at-50%-selectivity and $\delta$ is equal to $a_{95} - a_{50}$.

The use of 19 is something of an enigma: if you do not use 19, the model does not predict the 95% of the selectivity, but it is not a simple relationship.

It should also be noted that the selectivity curve used does not need to be the particular version of the logistic that is used herein. Instead one can use other logistic equations, and as they do not use the rather arcane -Ln(19) they are sometimes easier to understand

Three parameters are estimated, both the selectivity paramters, $a_{50}$ and $\delta$, and the average fully-selected $F$ value. $\delta$ is used instead of $a_{95}$ as this reduces the correlation between the two selectivity paarameters and helps avoid $a_{50}$ values larger than the  $a_{95}$ values. These parameters can be estimated through minimizing the multinomial negative log-likelihood:

$$-veLL=-\sum\limits_{y}\sum\limits_{a}{{O}_{y,a}}log\left(\frac{{\hat{N}_{a}}}{\sum{{\hat{N}_{a}}}} \right)$$

where $O_{y,a}$ is the observed proportion of numbers-at-age $a$ in the catch during year $y$, and ${\hat{N}_{a}}$ is the expected or predicted proportion of numbers-at-age $a$ in the catch during year $y$. The predicted numbers-at-age are calculated using the Baranov catch equation:

$${\hat{N}_a}=\frac{{N_a}{s_a}F}{(s_aF+M)}\left(1-e^{-(s_aF+M)} \right )$$