Chapter 26 Developing Custom Management Procedures
DLMtool was designed to be extensible in order to promote the development of new Management Procedures. In this chapter we design a series of new Management Procedures that include spatial controls and input controls in the form of size limit restrictions.
If you wish, you can also add your newly developed MPs to the DLMtool package so they are accessible to other uses. Of course you will be credited as the author. Please contact us for details how to do this.
As we saw before, real data are stored in a class of objects Data.
The DLMtool MSE function generates simulated data and puts it in exactly the same format as real data. This is highly desirable because it means that the same MP code that is tested in the MSE can then be used to make management recommendations.
If an MP is coded incorrectly it may catastrophically fail MSE testing and will therefore be excluded from use in management.
26.1 The Anatomy of an MP
Let’s examine an existing output MP to identify the MP data requirements.
## [1] "AvC" "BK" "BK_CC" "BK_ML" "CC1" "CC2" "CC3" "CC4"
## [9] "CC5" "CompSRA" "CompSRA4010" "CurC" "DAAC" "DBSRA" "DBSRA_40" "DBSRA4010"
## [17] "DCAC" "DCAC_40" "DCAC_ML" "DCAC4010" "DCACs" "DD" "DD4010" "DepF"
## [25] "DynF" "Fadapt" "Fdem" "Fdem_CC" "Fdem_ML" "Fratio" "Fratio_CC" "Fratio_ML"
## [33] "Fratio4010" "GB_CC" "GB_slope" "GB_target" "Gcontrol" "HDAAC" "ICI" "ICI2"
## [41] "Iratio" "Islope1" "Islope2" "Islope3" "Islope4" "IT10" "IT5" "Itarget1"
## [49] "Itarget2" "Itarget3" "Itarget4" "ITM" "L95target" "Lratio_BHI" "Lratio_BHI2" "Lratio_BHI3"
## [57] "LstepCC1" "LstepCC2" "LstepCC3" "LstepCC4" "Ltarget1" "Ltarget2" "Ltarget3" "Ltarget4"
## [65] "MCD" "MCD4010" "Rcontrol" "Rcontrol2" "SBT1" "SBT2" "SPmod" "SPMSY"
## [73] "SPslope" "SPSRA" "SPSRA_ML" "YPR" "YPR_CC" "YPR_ML" "avgMP" "TCPUE"
## [81] "THC"Since we’ve seen it used as a default MP in lots of the examples above, lets learn more about DCAC
We can even see all the code for this MP by simply typing the name of the MP into the console (this is a fantastic advantage of using R - there is complete transparency about package functions):
## function(x, Data, reps = 100, plot=FALSE) {
## rundcac <- DCAC_(x, Data, reps, updateD=TRUE)
## TAC <- TACfilter(rundcac$dcac)
##
## if (plot) DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))
##
## Rec <- new("Rec")
## Rec@TAC <- TAC
## Rec
## }
## <bytecode: 0x000001d484647608>
## <environment: namespace:DLMtool>
## attr(,"class")
## [1] "MP"“Crikey that looks complicated!” might be your first reaction. However this output MP function is easily demystified.
Like all MPs it has four arguments: x, Data, reps and plot (the last argument was added recently and is optional).
The argument x is the position in the Data object. When real data are stored in a Data object, there is only one position - there is only one real data set.
However, in MSE we conduct many simulations and x refers to simulated data from simulation number x. Any single parameters such as natural mortality rate (Mort) are a vector (nsim long). See Data@Mort[x] in the DCAC code. Any time series such as annual catches or relative abundance indices, are a matrix of nsim rows and nyears columns.
A range of objects of class Data are available:
## [1] "Atlantic_mackerel" "China_rockfish" "Cobia" "Example_datafile" "Gulf_blue_tilefish"
## [6] "ourReefFish" "Red_snapper" "SimulatedData" "Simulation_1" "China_rockfish2"
## [11] "Data" "Madeup" "Recs"For simplicity lets use a Data object with just one simulation, Simulation_1 and rename it Data
Since there is only one simulation in this data set (1 position) we can now see a single value of natural mortality rate:
## [1] 0.2244735And a matrix of catches with only 1 row:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 4.275057 12.43761 14.63192 35.31725 28.69802 30.84651 24.14059 36.78335 29.27517 38.18088 59.30242 56.08995 37.96849
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 51.84985 60.76729 41.53713 39.31114 57.9673 57.2248 72.37596 80.69301 76.63558 59.37687 48.82643 50.38803 80.71158
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39]
## [1,] 53.35875 74.31955 48.83082 43.348 65.17864 49.94281 47.38492 45.23197 80.92438 51.91477 29.36491 43.44834 49.46923
## [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 50.33217 49.53639 39.28779 27.31767 38.97092 51.1054 37.34677 37.33128 24.24094 23.47756 21.08158We could generate a single TAC recommendation from these data using DCAC by specifying position 1 (there is only 1 simulation) and by setting reps=1 (we want a single DCAC TAC recommendation)
## TAC (median)
## 37.35232If we wanted a stochastic estimate of the TAC we could increase the number of reps:
26.2 A Constant Catch MP
We’ve now got a better idea of the anatomy of an MP. It is a function that must accept three arguments (we will ignore plot for now):
- x: a simulation number
- Data: an object of class
Data - reps: the MP can provide a sample of TACs
repslong.
Let’s have a go at designing our own custom MP that can work with DLMtool. We’re going to develop an MP that sets the TAC as the ‘3rd highest catch’.
We decide to call our function THC
THC<-function(x, Data, reps){
# Find the position of third highest catch
THCpos<-order(Data@Cat[x,],decreasing=T)[3]
# Make this the mean TAC recommendation
THCmu<-Data@Cat[x,THCpos]
# A sample of the THC is taken according to a fixed CV of 10%
TACs <- THCmu * exp(rnorm(reps, -0.1^2/2, 0.1)) # this is a lognormal distribution
Rec <- new("Rec") # create a 'Rec'object
Rec@TAC <- TACs # assign the TACs to the TAC slot
Rec # return the Rec object
}To recap that’s just seven lines of code:
THC<-function(x, Data, reps){
THCpos<-order(Data@Cat[x,],decreasing=T)[3]
THCmu<-Data@Cat[x,THCpos]
Rec <- new("Rec")
Rec@TAC <- THCmu * exp(rnorm(reps, -0.1^2/2, 0.1))
Rec
}We can quickly test our new MP for the example Data object
## [1] 79.77229 84.19343 83.35881 77.35208 87.12422 83.23256 74.00345 95.07166 75.44146 80.00517Now that we know it works, to make the function compatible with the DLMtool package we have to assign it the class ‘MP’ so that DLMtool recognizes the function as a management procedure
If we want to run the MSE in parallel we need to export the newly created function to the cluster:
26.3 A More Complex MP
The THC MP is simple and frankly not a great performer (depending on depletion, life-history, adherence to TAC recommendations).
Let’s innovate and create a brand new MP that could suit a catch-data-only stock like Indian Ocean Longtail tuna!
It may be possible to choose a single fleet and establish a catch rate that is ‘reasonable’ or ‘fairly productive’ relative to current catch rates. This could be for example, 40% of the highest catch rate observed for this fleet or, for example, 150% of current cpue levels.
It is straightforward to design an MP that will aim for this target index level by making adjustments to the TAC.
We will call this MP TCPUE, short for target catch per unit effort:
TCPUE<-function(x,Data,reps){
mc<-0.05 # max change in TAC
frac<-0.3 # target index is 30% of max
nyears<-length(Data@Ind[x,]) # number of years of data
smoothI<-smooth.spline(Data@Ind[x,]) # smoothed index
targetI<-max(smoothI$y)*frac # target
currentI<-mean(Data@Ind[x,(nyears-2):nyears]) # current index
ratio<-currentI/targetI # ratio currentI/targetI
if(ratio < (1 - mc)) ratio <- 1 - mc # if currentI < targetI
if(ratio > (1 + mc)) ratio <- 1 + mc # if currentI > targetI
Rec <- new("Rec")
Rec@TAC <- Data@MPrec[x] * ratio * exp(rnorm(reps, -Data@CV_Ind[x]^2/2, Data@CV_Ind[x]))
Rec
}The TCPUE function simply decreases the past TAC (stored in Data@MPrec) if the index is lower than the target and increases the TAC if the index is higher than the target.
All that is left is to make it compatible with DLMtool:
26.4 Beyond the Catch Limit
All management procedures return an object of class ‘Rec’ that contains 13 slots:
## [1] "TAC" "Effort" "Spatial" "Allocate" "LR5" "LFR" "HS" "Rmaxlen" "L5" "LFS"
## [11] "Vmaxlen" "Fdisc" "Misc"We’ve already seen the TAC slot in the previous exercise. The remaining slots relate to various forms of input control:
- Effort (total allowable effort (TAE) relative to last historical year)
- Spatial - Fraction of each area that is open
- Allocate - Allocation of effort from closed areas to open areas
- LR5 - Length at 5% retention
- LFR - Length at 100% retention
- HS - Upper slot limit
- Rmaxlen - Retention of the maximum length class
- L5 - Length at 5% selection (e.g a change in gear type)
- LFS - Length at 100% selection (e.g a change in gear type)
- Vmaxlen - Selectivity of the maximum length class
- Fdisc - Update the discard mortality if required
- Misc - An optional slot for storing additional information
The curE MP just keeps effort constant at current levels:
## function(x, Data, reps, plot=FALSE) {
## # current effort
## rec <- new("Rec") # create recommendation object
## rec@Effort <- 1 #* Data@MPeff[x]
## if (plot) curE_plot(x, rec, Data)
## rec
## }
## <bytecode: 0x000001d4843e9ca0>
## <environment: namespace:DLMtool>
## attr(,"class")
## [1] "MP"Note that only the Effort slot in the Rec object is populated in this case.
To highlight the differences among Input control MPs examine spatial control MP MRreal that closes area 1 to fishing and reallocates fishing to the open area 2:
## function(x, Data, reps, plot=FALSE) {
## # A Marine reserve in area 1 with spatial reallocation of effort
##
## rec <- new("Rec") # create recommendation object
## rec@Allocate <- 1
## rec@Spatial <- c(0, rep(1, Data@nareas-1))
##
## if (plot) barplot(rec@Spatial, xlab="Area", ylab="Fraction Open", ylim=c(0,1),
## names=1:Data@nareas)
## return(rec)
## }
## <bytecode: 0x000001d4fbd54448>
## <environment: namespace:DLMtool>
## attr(,"class")
## [1] "MP"In contrast MRnoreal does not reallocate fishing effort:
## function(x, Data, reps, plot=FALSE) {
## # A Marine reserve in area 1 with no spatial reallocation of effort
##
## rec <- new("Rec") # create recommendation object
## rec@Allocate <- 0
## rec@Spatial <- c(0, rep(1, Data@nareas-1))
##
## if (plot) barplot(rec@Spatial, xlab="Area", ylab="Fraction Open", ylim=c(0,1),
## names=1:Data@nareas)
## return(rec)
## }
## <bytecode: 0x000001d4fbd5c170>
## <environment: namespace:DLMtool>
## attr(,"class")
## [1] "MP"The MP matlenlim only specifies the parameters of length retention using an estimate of length at 50% maturity (Stock@L50):
## function(x, Data, reps, plot=FALSE) {
## # Knife-edge vulnerability at estimated length-at-maturity
## rec <- new("Rec") # create recommendation object
## rec@LFR <- Data@L50[x] # new length at full retention
## rec@LR5 <- rec@LFR * 0.95 # new length at 5% retention
## if(plot) size_lim_plot(x, Data, rec)
## # other slots aren't specified so remain unchanged
## rec
## }
## <bytecode: 0x000001d4fc43cfe0>
## <environment: namespace:DLMtool>
## attr(,"class")
## [1] "MP"26.4.1 An Example Effort Control
Here we will copy and modify the MP we developed earlier to specify a new version of the target catch per unit effort MP (TCPUE) that provides effort recommendations:
TCPUE_e<-function(x,Data,reps){
mc<-0.05 # max change in TAC
frac<-0.3 # target index is 30% of max
nyears<-length(Data@Ind[x,]) # number of years of data
smoothI<-smooth.spline(Data@Ind[x,]) # smoothed index
targetI<-max(smoothI$y)*frac # target
currentI<-mean(Data@Ind[x,(nyears-2):nyears]) # current index
ratio<-currentI/targetI # ratio currentI/targetI
if(ratio < (1 - mc)) ratio <- 1 - mc # if currentI < targetI
if(ratio > (1 + mc)) ratio <- 1 + mc # if currentI > targetI
rec <- new("Rec")
rec@Effort <- Data@MPeff[x] * ratio
rec
}There have been surprisingly few changes to make TCPUE an input control MP that sets total allowable effort.
- We have had to use stored recommendations of effort in the
Data@MPeffslot, and - The final line of the MP is our input control recommendatation that only modified the Effort.
That is all.
Again, we need to assign our new function to class MP and export it to the cluster:
Let’s test the two MPs and see how they peform:
## Exporting custom MPs in global environment## Running MSE in parallel on 10 processors## MSE completed
## PNOF B50 LTY VY
## TCPUE 50.3 60.4 52.9 89.6
## TCPUE_e 58.4 82.0 90.0 2.1