datalowSA: A package to facilitate application of data poor stock assessments (datalowSA).

(This package is well described in https://rdrr.io/github/haddonm/datalowSA/f/README.md)

Surplus production model (spm https://rdrr.io/github/haddonm/datalowSA/man/spm.html) is a component of the package datalowSA. “It is one of the simplest analytic methods available that provides for a full fish stock assessment of the population dynamics of the stock being examined.”. “SPM represents the dynamics using a Schaefer or a Fox model. The outputs include predicted Biomass, year, catch, cpue, predicted cpue, contributions to q, ssq, and depletion levels. Generally it would be more sensible to use simpspm when fitting a Schaefer model and simpfox when fitting a Fox model as those functions are designed to generate only the predicted cpue required by the functions ssq and negLL, but the example shows how it could be used. The function spm is used inside displayModel and could be used alone, to generate a full-list of model outputs after the model has been fitted.”.

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 at https://www.rstudio.com/products/rstudio/download/.

The supported version of SPM is part of a set of stock assessment tools within Dr Haddon’s datalowSA package. Therefore one needs to first install datalowSA from within R.

Note for R4.0.0.0+:

R4.0.0.0+ is available recently. Under this version, Rtools40 and RStudio1.2.5042+ (if you use RStudio) are required. If you get the error messages of missing packages such as “backports” or “fs”, you need to install the packages manually.

Set up checking/installing packages functions

# 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(github_pkg){
  pkg = tail(unlist(strsplit(github_pkg,"/")),1)
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)>0) 
    devtools::install_github(github_pkg,build_vignettes=T)
  sapply(pkg, require, character.only = TRUE)
}

Installing packages from GitHub directly.

Check and install “devtools”

Check_Install.packages("devtools")
## [1] "'devtools' has been installed already!"
## Loading required package: devtools
## Loading required package: usethis
## Error: package or namespace load failed for 'usethis' in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
##  namespace 'rlang' 1.0.5 is already loaded, but >= 1.0.6 is required
## devtools 
##    FALSE
# IsExst_Devtools = check.packages("devtools")
# 
# if(!IsExst_Devtools) install.packages("devtools") else print("'devtools' has been installed already!")

Check and install datalowSA from github

# Check if existed and if not, install the required package TMBhelper:
github_pkg = "haddonm/datalowSA"

Github_CheckInstall.packages(github_pkg=github_pkg)
## Loading required package: datalowSA
## datalowSA 
##      TRUE
#install_github("jabbamodel/JABBA")

NB Once you have installed datalowSA, you don’t need to run the install.packages code chunk again. Luckily it will just give an error stating it won’t re-install if your version is up to date. 

library(datalowSA)

It should be noted that there are also quite good vignettes for this package by typing browseVignettes(package = “datalowSA”). In that file for SPM, Dr Haddon describes several steps of which we only show a few below with our example dataset.

Reading in the data

An 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.

Catches <- read.csv("Catches.csv")
print(head(Catches))
##   Year Fleet Value
## 1 1960     1     0
## 2 1961     1     0
## 3 1962     1     0
## 4 1963     1     0
## 5 1964     1     0
## 6 1965     1     0
CPUEs <- read.csv("CPUE.csv")
print(head(CPUEs))
##   Year Fleet    Value  CV
## 1 1971     1 2.641477 0.3
## 2 1972     1 2.224941 0.3
## 3 1973     1 2.526852 0.3
## 4 1974     1 1.630703 0.3
## 5 1975     1 5.053905 0.3
## 6 1976     1 2.844623 0.3

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-1. For this package, we provide the catch data in tons and cpue data in tons per effort unit.

Making input data into the model

Input data expected as a column of years, catches and cpue (as we have set up in Example). The model also needs some idea of the population growth rate (r), K and the initial biomass (Binit) of the population if the data are not collected from start of the fishery.

In this example,

inputDat = merge(Catches[,c("Year","Value")],
                 CPUEs[,c("Year","Value")],by="Year")
year <- inputDat$Year
catch <- inputDat$Value.x
cpue <- inputDat$Value.y
dat <- makespmdata(cbind(year,catch,cpue))
#pars <- c(0.241,51049)  # from cMSY simulation's mean_r and mean_K
pars <- c(0.281,45236)  # from CMSY actual medium r and medium K. The data collected from the beginning of the fishery, therefore no Binit given.

Running the model

The syntax for using spm function as:

spm(inp=pars,indat=dat,schaefer = TRUE)

However as mentioned above, the function spm is used inside ‘displayModel’ and displayModel takes a set of parameters and the spmdat matrix and plots the predicted depletion, catch, the surplus production, and the CPUE and the model fit to CPUE.

Schaefer Model

To fit catch and cpue data using Shaefer model and the model set default LRP as 0.2 and TRP as 0.48. The model results are plotted.

ans <- displayModel(pars,dat,schaefer = TRUE)

str(ans)
## List of 11
##  $ Dynamics :List of 5
##   ..$ outmat    : num [1:31, 1:8] 45236 45137 44966 44050 43385 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:31] "1971" "1972" "1973" "1974" ...
##   .. .. ..$ : chr [1:8] "ModelB" "Catch" "CPUE" "PredCE" ...
##   ..$ q         : num 4.96e-05
##   ..$ msy       : num 3178
##   ..$ parameters: num [1:2] 2.81e-01 4.52e+04
##   ..$ sumout    : Named num [1:9] 2.81e-01 4.52e+04 4.52e+04 3.18e+03 1.00 ...
##   .. ..- attr(*, "names")= chr [1:9] "r" "K" "B0" "msy" ...
##  $ BiomProd : num [1:100, 1:2] 100 556 1012 1468 1924 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "x" "y"
##  $ rmseresid: num 1.17
##  $ MSY      : num 3178
##  $ Bmsy     : num 22440
##  $ Dmsy     : num 0.496
##  $ Blim     : num 9218
##  $ Btarg    : num 21528
##  $ Ctarg    : num 3170
##  $ Dcurr    : Named num 0.289
##   ..- attr(*, "names")= chr "2001"
##  $ rmse     : logi NA

To optimise the parameters and apply Schaefer spm modeling again with optimised parameters, and to output model results in text.

bestSP <- optim(par=pars,fn=ssq,callfun=simpspm,indat=dat)
ans <- displayModel(bestSP$par,dat,schaefer=TRUE)

str(ans)
## List of 11
##  $ Dynamics :List of 5
##   ..$ outmat    : num [1:31, 1:8] 22670 22571 22446 21620 21380 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:31] "1971" "1972" "1973" "1974" ...
##   .. .. ..$ : chr [1:8] "ModelB" "Catch" "CPUE" "PredCE" ...
##   ..$ q         : num 0.000108
##   ..$ msy       : num 4238
##   ..$ parameters: num [1:2] 7.48e-01 2.27e+04
##   ..$ sumout    : Named num [1:9] 7.48e-01 2.27e+04 2.27e+04 4.24e+03 1.00 ...
##   .. ..- attr(*, "names")= chr [1:9] "r" "K" "B0" "msy" ...
##  $ BiomProd : num [1:100, 1:2] 100 328 556 784 1012 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "x" "y"
##  $ rmseresid: num 1.14
##  $ MSY      : num 4238
##  $ Bmsy     : num 11271
##  $ Dmsy     : num 0.497
##  $ Blim     : num 4432
##  $ Btarg    : num 10815
##  $ Ctarg    : num 4229
##  $ Dcurr    : Named num 0.422
##   ..- attr(*, "names")= chr "2001"
##  $ rmse     : logi NA

Fox Model

To fit catch and cpue data using the Fox model and set the default LRP to 0.2 and TRP to 0.48, and plot the resulst:

ans <- displayModel(pars,dat,schaefer = FALSE)

str(ans)
## List of 11
##  $ Dynamics :List of 5
##   ..$ outmat    : num [1:31, 1:8] 45236 45137 44966 44050 43390 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:31] "1971" "1972" "1973" "1974" ...
##   .. .. ..$ : chr [1:8] "ModelB" "Catch" "CPUE" "PredCE" ...
##   ..$ q         : num 3.4e-05
##   ..$ msy       : num 4676
##   ..$ parameters: num [1:2] 2.81e-01 4.52e+04
##   ..$ sumout    : Named num [1:9] 2.81e-01 4.52e+04 4.52e+04 4.68e+03 1.00e-08 ...
##   .. ..- attr(*, "names")= chr [1:9] "r" "K" "B0" "msy" ...
##  $ BiomProd : num [1:100, 1:2] 100 556 1012 1468 1924 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "x" "y"
##  $ rmseresid: num 1.48
##  $ MSY      : num 4676
##  $ Bmsy     : num 16513
##  $ Dmsy     : num 0.365
##  $ Blim     : num 9218
##  $ Btarg    : num 21528
##  $ Ctarg    : num 4492
##  $ Dcurr    : Named num 0.795
##   ..- attr(*, "names")= chr "2001"
##  $ rmse     : logi NA

To optimise the parameters and do Fox spm modeling again with optimised parameters, ando to output model results in text.

bestSP <- optim(par=pars,fn=ssq,callfun=simpfox,indat=dat)
ans <- displayModel(bestSP$par,dat,schaefer=FALSE)

str(ans)
## List of 11
##  $ Dynamics :List of 5
##   ..$ outmat    : num [1:31, 1:8] 30574 30475 30309 29405 28798 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:31] "1971" "1972" "1973" "1974" ...
##   .. .. ..$ : chr [1:8] "ModelB" "Catch" "CPUE" "PredCE" ...
##   ..$ q         : num 9.06e-05
##   ..$ msy       : num 3741
##   ..$ parameters: num [1:2] 3.33e-01 3.06e+04
##   ..$ sumout    : Named num [1:9] 3.33e-01 3.06e+04 3.06e+04 3.74e+03 1.00e-08 ...
##   .. ..- attr(*, "names")= chr [1:9] "r" "K" "B0" "msy" ...
##  $ BiomProd : num [1:100, 1:2] 100 408 716 1023 1331 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "x" "y"
##  $ rmseresid: num 1.12
##  $ MSY      : num 3741
##  $ Bmsy     : num 11182
##  $ Dmsy     : num 0.366
##  $ Blim     : num 6256
##  $ Btarg    : num 14568
##  $ Ctarg    : num 3592
##  $ Dcurr    : Named num 0.389
##   ..- attr(*, "names")= chr "2001"
##  $ rmse     : logi NA

Explanation of the results

spm generally produces six plots for each of the model run. They include “ExploitB Depletion”,“Scaled CPUE”,“LN Residuals”,“Catch”,“Surplus Production - Biomass” and “Surplus Production - Depletion”. From these graphs one can already get quite a lot of the required management information, including MSY and the depletion rate. We have two models examples and both use 0.2B0 and 0.482B0 as limit and target reference point proxies.

There are two options to set the parameters for the modeling for each of the models (Schaefer and Fox). Option 1 is to use given parameters (r, K); Option 2 is to optimise the parameters given the input values for r and K.

  • ExploitB Depletion displays the trajectory of the exploitable biomass depletion rate. The red line is the LRP and the green line is the TRP. In Schaefer model results, two options are similar and

  • Catch displays the catch series with MSY as the red line.

  • Scaled CPUE displays the model fit to the CPUE data. The red line is the model predictions and the black dots are the CPUE data points.

  • Surplus Production - Biomass displays surplus production vs biomass. Vertical red lines represent Blim, Btarg and BMSY and the horizontal red dash line is MSY.

  • LN Residuals displays the residuals (scaled observe CPUE / scaled predicted CPUE) with a value of RMSE (root mean square error)

  • Surplus Production - Depletion displays the surplus production vs deletion function. The vertical lines represent the LRF (red), the TRP (green) and MSY point (dark green). The horizontal lines represent MSY (red) and Ctarg (green).

Key references

SPM Methodology

Schaefer, M.B. 1954. Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Bulletin, Inter-American Tropical Tuna Commission 1:25-56. https://www.iattc.org/PDFFiles/Bulletins/_English/Vol-1-No-2-1954-SCHAEFER,%20MILNER%20B._Some%20aspects%20of%20the%20dynamics%20of%20populations%20important%20to%20the%20management%20of%20the%20commercial%20marine%20fisheries.pdf

Schaefer, M.B. 1957. A study of the dynamics of the fishery for yellowfin tuna in the Eastern Tropical Pacific Ocean. Bulletin, Inter-American Tropical Tuna Commission 2:247-285. https://www.iattc.org/PDFFiles/Bulletins/_English/Vol-2-No-6-1957-SCHAEFER,%20MILNER%20B_A%20study%20of%20the%20dynamics%20of%20the%20fishery%20for%20yellowfin%20tuna%20in%20the%20eastern%20tropical%20Pacific%20Ocean.pdf

Walters, C.J., Martell, S.J.D. and J. Korman 2006. A stochastic approach to stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 63:212-223. https://doi.org/10.1139/f05-213

SPM Evaluation

Froese, R., Demirel, N., Coro, G., Kleisner, K. and H. Winker. 2017. Estimating fisheries reference points from catch and resilience. Fish and Fisheries. 18:506-526. https://doi.org/10.1111/faf.12190.

datalowSA package

Haddon, M., Burch, P., Dowling, N. and R. Little, R. 2019. Reducing the Number of Un-defined Species in Future Status of Australian Fish Stocks Reports: Phase Two - training in the assessment of data-poor stocks. FRDC Final Report 2017/102. CSIRO Oceans and Atmosphere and Fisheries Research Development Corporation. Hobart. 125 p.