1 Introduction

This document summarizes the new developments of the DLMtool R package to extend its capabilities to model fisheries with multiple stocks and multiple fishing fleets, and describes the application of the newly developed extension to three case-study fisheries from the California Department of Fish and Wildlife (CDFW).

DLMtool was developed as an R package designed to make management strategy evaluation (MSE) accessible to data-limited fisheries. The DLMtool operating model (OM) assumes a single stock (e.g., usually female component of the population) and a single fishing fleet (i.e., all fishing fleets are modeled together in aggregate). The MSEtool R package was developed as an extension to DLMtool for the purpose of allowing more complex features typical of more data-rich fisheries to be included in the modeling framework, such as quantitative stock assessments and complex spatial arrangements.

MultiMSE was developed to extend the modeling framework to include more complex biology (e.g., separate male and female populations, and hermaphroditism), interactions between multiple stocks (e.g., MICE models for evaluating ecosystem interactions), and options for fleet-specific management (e.g., separate recreational and commercial fleets). MultiMSE was developed as a proof-of-concept, to evaluate the viability and practicality of including more complex life-history patterns, and flexible management options within the DLMtool framework.

The aim of this project was to further develop the MultiMSE extension so that it can used to evaluate the performance of management options for CDFW fisheries with sexually dimorphic growth (male and female stocks), stock complexes (multiple species managed as a single unit), and fleet-specific management for fisheries that are exploited by fleets with different characteristics.

The purpose of this document is to describe the newly developed features of the model and the application of the multi-stock & multi-fleet model to three CDFW case study fisheries: 1) Rock Crab, 2) Lobster, and 3) Halibut, and to identify any remaining issues with applying the model to these fisheries.

1.1 New Package Structure

As mentioned above, the multi-stock & multi-fleet model was initially developed as a proof-of-concept and included in the MSEtool R package. Early on in this project it became clear that the development of the multi-stock & multi-fleet model required a substantial re-structure of the code in the DLMtool and MSEtool packages.

The new package structure now includes 3 R packages:

MSEtool (v3.0.0): The MSEtool package has been completely revised and now includes all functions for generating, simulating, and projecting single stock OMs (previously in DLMtool) and multi-stock OMs (previously in earlier version of MSEtool). This new version of MSEtool includes all existing functions for single stock OMs and replaces the earlier versions of the DLMtool R package. This version of MSEtool is only compatible with DLMtool v6+.
DLMtool (v6.0.0): This new version of DLMtool is now only a collection of data-limited management procedures. This version requires MSEtool v3+ as a dependency.
SAMtool (v1.0.0): The Stock Assessment Methods toolkit is a new package that contains the data-rich and data-moderate assessment methods and management procedures that were previously available in earlier versions of MSEtool. SAMtool also contains the operating model scoping function that has been updated and re-named RCM (Rapid Conditioning Model; previously SRA_scope). SAMtool requires MSEtool v3+ as a dependency.

1.2 Installing and Loading the New Packages

Development versions of the three packages are available to be installed from GitHub.

The Rtools software is required to install the new package versions from GitHub. We also recommend using the latest version of R (at least R 4.0.0).

While only MSEtool v6.0.0 is required to run the multi-stock, multi-fleet MSE, we recommend also installing both the new version of DLMtool and SAMtool so that the data-limited and data-rich management procedures are available to use.

The three packages can be installed directly from GitHub. This only needs to be done once (and repeated only if the development versions of the packages are updated):

# install.packages('devtools') # run this if devtools isn't installed
library(devtools)
install_github('Blue-Matter/MSEtool') 
install_github('Blue-Matter/DLMtool') 
install_github('Blue-Matter/SAMtool')

At the beginning of each new R session, the 3 packages need to be loaded:

# load package into R session (needs to be done for each new session)
library('MSEtool')

## Loading required package: snowfall

## Loading required package: snow

library('DLMtool')
library('SAMtool')

## 
## Attaching package: 'SAMtool'

## The following object is masked from 'package:MSEtool':
## 
##     userguide

We intend to publish the three new and updated packages on CRAN in the next few weeks, together with a new global package called openMSE. openMSE will act as an umbrella package that requires MSEtool, DLMtool, and SAMtool as dependencies; i.e., installing openMSE from CRAN will automatically install MSEtool, DLMtool, and SAMtool and make all MSE functions and MPs available.

2 Rock Crab Case Study

The rock crab case study involves three species, Brown, Red, and Yellow, which are managed as a single stock complex. Male and female rock crab have different growth patterns and mature at different sizes, and modeling all stocks together will involve 6 separate stocks.

Given that there are no interactions that are being modeled between the species, and the management options being explored by CDFW are focused on size regulations and effort controls, it may not be necessary to model the three stocks simultaneously. Wherever possible, we recommend limiting the number of stocks that are modeled together, as each additional stock significantly increases the run time of the model and may lead to increased convergence issues.

This case study includes the female and male stocks for the Red rock crab. This analysis can be repeated for the Brown and Yellow rock crab case studies.

The case study can be developed further if CDFW wishes to explore management options that will affect the dynamics of the three species simultaneously.

2.1 Set up `MOM` Object

2.1.1 Import Stocks and create Stock Objects

Note that the some redundant slots have now been removed from the Stock and Fleet objects in OMtool (e.g., M2, Mgrad, Amplitude, Period, SelYears, etc). The Excel files received from CDFW have been slightly modified to account for these changes, and will be distributed with this document.

 # directory where the OM Excel files are located
OMdir <- 'G:/My Drive/1_PROJECTS/CDFW_Multi_Fleet/RockCrabDLM'

# Import the Female and Male Red Crab Stocks
RedF_Stock <- XL2Stock(file.path(OMdir, 'OMRedF.xlsx'), msg=FALSE)
RedM_Stock <- XL2Stock(file.path(OMdir, 'OMRedM.xlsx'), msg=FALSE)

# a list of stock objects (female stock in first position)
Stocks <- list(RedF_Stock, RedM_Stock)

2.1.2 Import Fleet and create Fleet Object

A single fleet is imported:

Red_Fleet <- XL2Fleet(file.path(OMdir, 'OMRedF.xlsx'), msg=FALSE)

And a Fleets list is created, where the Red_Fleet is applied to both the female and male stocks:

Fleets <- list(
  list(Red_Fleet), # female stock
  list(Red_Fleet) # male stock
)

2.1.3 Set up 2-sex model

For a 2-sex model it is necessary to set up the model so that recruitment is determined by the spawning stock biomass of the female population. This is done using an argument to the multi-OM object (MOM) called SexPars:

SexPars <- list()
# Stock (Row) spawn from SSB contributed by Stock (Column)
SexPars$SSBfrom <- matrix(c(1,0,
                            1,0),
                          byrow=TRUE,nrow=2)

SexPars$SSBfrom

##      [,1] [,2]
## [1,]    1    0
## [2,]    1    0

SexPars$SSBfrom is a 2x2 matrix specifying recruitment to females (row 1) and males (row 2) comes from female spawning biomass (column 1).

Note that it is important that the SexPars$SSBfrom matches the order of the stocks in Stocks (ie., female in first position in list Stocks, and first column in SexPars$SSBfrom.

2.1.4 Create `MOM` Object

The multi-stock operating model is an object of class MOM. Similar to the OM object, it requires a list of Stock, Fleet, Observation (Obs), and Implementation (Imp) objects, as well as the number of simulations to run.

For simplicity, we will use the generic observation object:

ObsList <- list(
  list(Generic_Obs), # Obs for female
  list(Generic_Obs) # Obs for male (the same)
)

And no implementation error:

ImpList <- list(
  list(Perfect_Imp),
  list(Perfect_Imp)
)

As this is for demonstration and testing purposes, the number of simulations is set quite low so that it runs quickly:

nsim <- 10

By default, the operating model calculates depletion in terms of spawning biomass relative to the equilibrium unfished spawning biomass $(\text{SB}_0)$. To calculate depletion in terms of vulnerable biomass, use the custom parameters argument:

cpars <- list()
# skip this line to calculate depletion in terms of SSB
cpars$control$D <-'VB'

The MOM object is created using the new function:

MOM <- new("MOM",
           Stocks = Stocks,
           Fleets = Fleets,
           Obs = ObsList,
           Imps = ImpList,
           CatchFrac = NULL,
           SexPars=SexPars,
           cpars=cpars,
           nsim = nsim,
           interval=1)

2.2 Run the Historical Simulations

The simulations for the historical and projection period can be run separately in OMtool. This is important as running the historical simulations can be quite time-consuming in some cases where the model has to optimize for depletion for a large number of stocks and fleets.

2.2.1 Simulate Historical Fishery

The simulation of the historical fishery is done with the SimulateMOM function (here we are using parallel processing; set parallel to FALSE if there are any problems with using parallel processing on your machine):

RCrabHist <- SimulateMOM(MOM, parallel = TRUE)

To avoid re-running the historical simulations each time, they can be saved to disk:

saveRDS(RCrabHist, file.path(OMdir, 'RCrabHist.rda'))

2.3 Define Management Procedures

The management options that have been proposed by CDFW for the rock crab fishery are static size limits, male only fishery (release of all female crabs), and effort reductions.

Permits for the crab fishery are limited, but there is currently latent effort in the fishery, so effort could increase. Modeling the changes in fishing effort in the future projections requires either a bio-economic model that describes how effort changes in response to catch rates, or simplifying assumptions for the fishing effort in the future years.

The latter approach is used in this case study. Two scenarios are included for future fishing effort:

Effort remains at the current level
Effort increases over the next 5 years to twice the current level

These scenarios are included as examples only, and without further discussion with CDFW and fishery stakeholders, they should not be considered credible scenarios for future effort.

2.3.1 Current Size Limit

The curE MP assumes that the current size limit of 4.25 inch (108 mm) remains the same and fishing effort remains at the current level.

A MP is defined that increases effort over the first 5 years of the projection period to twice the current level (current size limit remains unchanged):

incE <- function(x, Data, ...) {
  rec <- new('Rec') # create recommendation object
  curYr <- max(Data@Year) - Data@LHYear # current year
  if (curYr <=5) {
    eff_inc <- 1 + 0.2 * curYr # effort increases 20% per year for first 5 years
  } else {
    eff_inc <- 2 # 100% increase after 5 years
  }
  rec@Effort <- eff_inc # populate Effort slot
  rec
}
class(incE) <- 'MP' # assign to class `MP`

2.3.2 Male Only Fishery & Current Size Limit

In a male only fishery, there is no retention of female catch. The MP does this by setting the retention curve to 0 for the female stock. Note that under this assumption, the female catch still suffers the discard mortality (if any) from crabs that are selected but not retained. To modify this assumption, the MP should change the selectivity curve instead of the retention curve.

Two MPs are defined, one with effort remaining at current levels, and the other with an increase in effort:

MOnly_curE <- function(x, Data, ...) {
  rec <- new('Rec') # create recommendation object
  # set length-at-5% retention to maximum size class (ie nothing is retained)
  rec@LR5 <- max(Data@CAL_bins) 
  rec@LFR <- rec@LR5 + 1
  rec@Effort <- 1
  rec
}
class(MOnly_curE) <- 'MP'

MOnly_incE <- function(x, Data, ...) {
  rec <- new('Rec') # create recommendation object
  # set length-at-5% retention to maximum size class (ie nothing is retained)
  rec@LR5 <- max(Data@CAL_bins) 
  rec@LFR <- rec@LR5 + 1
  
  curYr <- max(Data@Year) - Data@LHYear # current year
  if (curYr <=5) {
    eff_inc <- 1 + 0.2 * curYr # effort increases 20% per year for first 5 years
  } else {
    eff_inc <- 2 # 100% increase after 5 years
  }
  rec@Effort <- eff_inc # populate Effort slot
  rec
}
class(MOnly_incE) <- 'MP'

2.3.3 Increased Size Limit

CDFW proposed evaluating the impact of an increased size limit for red crab - increasing from 108 mm to 140 mm.

As the stocks are independent, this MP can be applied separately to the Brown and Yellow crab stocks to evaluate the impact a 140 mm size limit would have on these stocks.

Again, two MPs are defined:

newSL_curE <- function(x, Data, ...) {
  rec <- new('Rec') # create recommendation object
  # set retention to ~138
  rec@LR5 <- 138 
  rec@LFR <- 140
  
  rec@Effort <- 1 # populate Effort slot
  rec
}
class(newSL_curE) <- 'MP'

newSL_incE <- function(x, Data, ...) {
  rec <- new('Rec') # create recommendation object
  # set retention to ~138
  rec@LR5 <- 138 
  rec@LFR <- 140
  
  curYr <- max(Data@Year) - Data@LHYear # current year
  if (curYr <=5) {
    eff_inc <- 1 + 0.2 * curYr # effort increases 20% per year for first 5 years
  } else {
    eff_inc <- 2 # 100% increase after 5 years
  }
  rec@Effort <- eff_inc # populate Effort slot
  rec
}
class(newSL_incE) <- 'MP'

2.3.4 Effort Limit

The final MP examines the impact of a reduction in fishing effort. The MP has been defined to reduced fishing effort in the future years by an arbitrary 25% from the current level.

Similar to the previous MPs, as the stocks are independent this MP can be applied separately to the Brown and Yellow crab stocks to evaluate the impact of reduced fishing effort on these stocks.

reduceE <- function(x, Data, ...) {
  rec <- new('Rec') # create recommendation object

  rec@Effort <- 0.75 # populate Effort slot
  rec
}
class(reduceE) <- 'MP'

2.4 Run Forward Projections

The forward projections are done with the ProjectMOM function. This function requires the historical simulations returned by SimulateMOM (class multiHist) and a list of management procedures.

Some of the MPs are stock-specific (ie. MOnly_curE and MOnly_incE) so the MPs must be specified as a stock-specific list:

MPlist <- list(
  # MPs for female stock
  c('curE', 'incE', 'newSL_curE', 'newSL_incE', 'MOnly_curE', 'MOnly_incE', 'reduceE'), 
   # MPs for male stock
  c('curE','incE', 'newSL_curE', 'newSL_incE', 'curE', 'incE', 'reduceE') 
)

Notice that the curE and incE MPs are repeated for the male stock in the same position that the female stock has the no-retention MPs.

# load the saved historical simulations
RCrabHist <- readRDS(file.path(OMdir, 'RCrabHist.rda')) 
RedCrab_MSE <- ProjectMOM(RCrabHist, MPs=MPlist)
saveRDS(RedCrab_MSE, file.path(OMdir, 'RedCrab_MSE.rda'))

2.5 Examine Results

The ProjectMOM function returns an object of class MMSE, which contains all of the information from the MSE run. This object is similar to the MSE object (which has been slightly modified in the new version of MSEtool), with the key difference that slots B_BMSY, F_FMSY, etc now have dimensions for the multiple stocks and/or fleets.

The MMSE object may continue to be revised over time. The current slots in the object are:

slotNames('MMSE')

##  [1] "Name"     "nyears"   "proyears" "nMPs"     "MPs"      "MPcond"  
##  [7] "MPrefs"   "nsim"     "nstocks"  "nfleets"  "Snames"   "Fnames"  
## [13] "Stocks"   "Fleets"   "Obss"     "Imps"     "OM"       "Obs"     
## [19] "B_BMSY"   "F_FMSY"   "B"        "SSB"      "VB"       "FM"      
## [25] "C"        "TAC"      "SSB_hist" "CB_hist"  "FM_hist"  "Effort"  
## [31] "PAA"      "CAA"      "CAL"      "CALbins"  "MSY_P"    "FMSY_P"  
## [37] "SSBMSY_P" "Misc"

The MMSE object contains a lot of information and the results can be visualized in many different ways.

Currently, the main way to plot the results for the MMSE object is to use the generic plot function. This creates a series of projection plots, showing a) Spawning biomass $(\text{SB})$ relative to $\text{SB}_\text{MSY}$ for each stock and each MP, b) fishing mortality $(F)$ relative to $F_\text{MSY}$ for each stock and MP, c) total yield relative to the first projection year for each stock and MP, and d) average yield (absolute) for each stock and fleet (and MP).

plot(RedCrab_MSE, maxcol=7, curyr=2018)

2.5.1 Performance Metrics and Trade-Off Plots

The example performance metrics functions included in MSEtool (e.g., functions of class PM and plotting functions such as TradePlot) are designed for single stock OMs. Because there are so many different ways to calculate performance statistics for multiple stock/fleet MSEs, custom performance metric methods must be written for objects of class MMSE.

The following are some example custom performance metric functions based on the performance metrics used by CDFW in single stock/fleet MSEs:

2.5.1.1 Performance Limit: Probability Female SB > 0.5SBMSY after 1 Mean Generation Time

This performance metric calculates the probability that the spawning biomass of the female stock $\left(\text{SB}\right)$ (stock 1) is greater than 0.5 $\text{SB}_{\text{MSY}}$ after the mean generation time (MGT; or 10 years in MGT < 10):

# PM 1: Prob B>0.5BMSY after 1 MGT (or 10 years) throught to last projection year
P50mgt <- function (MMSE = NULL, Ref = 0.5, Yrs = NULL) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  
  mgt <- ceiling(MMSE@OM[[1]][[1]]$MGT)
  mgt[mgt<10] <- 10 #if MGT is less than 10 years, default to 10 years
  capMGT <- round(mean(mgt))
  mgt <- mean(mgt)
  
  PMobj <- new("PMobj")
  PMobj@Name <- paste0("Overfished Limit 1: Probabilty Female SB > ",Ref*100, "% SB_MSY (after ", capMGT, " years)")
  PMobj@Caption <-  paste0("Overfished Limit 1: \nProb. Female SB > ", Ref*100, "% SB_MSY (after ", capMGT, " years)")
  PMobj@Ref <- Ref

  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
    
  PMobj@Stat <- matrix(NA, nrow=MMSE@nsim, ncol=MMSE@nMPs)
  PMobj@Prob <- apply(MMSE@B_BMSY[,1, ,mgt:Yrs[2]] > PMobj@Ref, 1:2, mean) 
  PMobj@Mean <- calcMean(PMobj@Prob)
  PMobj@MPs <- MPnames
  PMobj
}
class(P50mgt) <- 'PM'

2.5.1.2 Performance Limit: Probability Female SB > 0.5SBMSY in Last 10 Years

This performance metric calculates the probability that the female $\text{SB}$ (stock 1) is greater than 0.5 $\text{SB}_{\text{MSY}}$ in the last 10 years of the projection period:

# PM 2: B>0.5BMSY in last 10 years
P50Last10 <- function(MMSE=NULL, Ref=0.5, Yrs=-10) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  
  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
  
  PMobj <- new("PMobj")
  PMobj@Name <- "Overfished Limit 2: Probability Female SB > 50% SB_MSY (last 10 years)"
  PMobj@Caption <- "Overfished Limit 2: \nProb. Female SB > 50% SB_MSY (last 10 years)"
  PMobj@Ref <- Ref
  PMobj@Stat <- MMSE@B_BMSY[,1, ,Yrs[1]:Yrs[2]]
  PMobj@Prob <- calcProb(PMobj@Stat > PMobj@Ref, MMSE)
  PMobj@Mean <- calcMean(PMobj@Prob)
  PMobj@MPs <- MPnames
  PMobj
}
class(P50Last10) <- 'PM'

2.5.1.3 Performance Object: Average Long-Term Yield (Last 10 Years)

This performance metric function calculates the overall average long-term yield relative to the reference yield (ie. summing both stocks):

# Calculate the Average LTY Yield - Last 10 Years
LTAvY <- function (MMSE = NULL, Ref = 1, Yrs = -10) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  
  ns <- MMSE@nstocks
  nf <- MMSE@nfleets
  nsim <- MMSE@nsim
  refY <- array(NA, dim=c(nsim, ns, nf))
  for (s in 1:ns){
    for (f in 1:nf) {
      refY[,s,f] <- MMSE@OM[[s]][[f]]$RefY
    }
  }
  RefY <- apply(refY, 1, sum)
  
  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
  
  Catch <- apply(MMSE@C[,,,,Yrs[1]:Yrs[2], drop=FALSE], c(1,4,5), sum)
  RefY <- array(RefY, dim=dim(Catch))
  RelCatch <- Catch/RefY
  
  PMobj <- new("PMobj")
  PMobj@Name <- "Long-term Average Yield"
  PMobj@Caption <- "Long-term Yield:\nAvg. Yield rel. to Maximum Yield (last 10 yrs)"
  PMobj@Ref <- Ref
  PMobj@Stat <- RelCatch
  PMobj@Prob <- apply(PMobj@Stat,c(1,2),mean,na.rm=T)  #avg for last 10 years for each MP/each sim
  PMobj@Mean <-apply(PMobj@Prob,2,mean,na.rm=T)
  PMobj@MPs <- MPnames
  PMobj
}
class(LTAvY) <- 'PM'

2.5.1.4 Performance Objective: Average Short-Term Yield (First 10 Years)

This performance metric function calculates the overall average short-term yield relative to the reference yield (ie. summing both stocks):

# Calculate the Average STY Yield - First 10 Years
STAvY <- function (MMSE = NULL, Ref = 1, Yrs = 10) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  
  ns <- MMSE@nstocks
  nf <- MMSE@nfleets
  nsim <- MMSE@nsim
  refY <- array(NA, dim=c(nsim, ns, nf))
  for (s in 1:ns){
    for (f in 1:nf) {
      refY[,s,f] <- MMSE@OM[[s]][[f]]$RefY
    }
  }
  RefY <- apply(refY, 1, sum)
  
  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
  
  Catch <- apply(MMSE@C[,,,,Yrs[1]:Yrs[2], drop=FALSE], c(1,4,5), sum)
  RefY <- array(RefY, dim=dim(Catch))
  RelCatch <- Catch/RefY
  
  PMobj <- new("PMobj")
  PMobj@Name <- "Short-term Average Yield"
  PMobj@Caption <- "Short-term Yield:\nAvg. Yield rel.to Maximum Yield (first 10 yrs)"
  PMobj@Ref <- Ref
  PMobj@Stat <- RelCatch
  PMobj@Prob <- apply(PMobj@Stat,c(1,2),mean,na.rm=T)  #avg for last 10 years for each MP/each sim
  PMobj@Mean <-apply(PMobj@Prob,2,mean,na.rm=T)
  PMobj@MPs <- MPnames
  PMobj
}
class(STAvY) <- 'PM'

2.5.1.5 Performance Objective: Average Yield (all but first 10 years)

This performance metric function calculates the overall average yield relative to the reference yield (ie. summing both stocks) over all but the first 10 years (assuming the MOM has a 50 year projection period):

# Calculate the Average Yield - All but first 10 Years
AvY <- function (MMSE = NULL, Ref = 1, Yrs = -40) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  
  ns <- MMSE@nstocks
  nf <- MMSE@nfleets
  nsim <- MMSE@nsim
  refY <- array(NA, dim=c(nsim, ns, nf))
  for (s in 1:ns){
    for (f in 1:nf) {
      refY[,s,f] <- MMSE@OM[[s]][[f]]$RefY
    }
  }
  RefY <- apply(refY, 1, sum)
  
  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
  
  Catch <- apply(MMSE@C[,,,,Yrs[1]:Yrs[2], drop=FALSE], c(1,4,5), sum)
  RefY <- array(RefY, dim=dim(Catch))
  RelCatch <- Catch/RefY
  
  PMobj <- new("PMobj")
  PMobj@Name <- "Average Yield"
  PMobj@Caption <- "Avg. Yield rel. to Max Yield (Yrs 11-50)"
  PMobj@Ref <- Ref
  PMobj@Stat <- RelCatch
  PMobj@Prob <- apply(PMobj@Stat,c(1,2),mean,na.rm=T)  #avg for last 10 years for each MP/each sim
  PMobj@Mean <-apply(PMobj@Prob,2,mean,na.rm=T)
  PMobj@MPs <- MPnames
  PMobj
}
class(AvY) <- 'PM'

2.5.1.6 Performance Objective: Probability Female SB > 0.8SBMSY in Last 10 Years

# B>0.8BMSY in last 10 years
P80Last10 <- function (MMSE = NULL, Ref = 0.8, Yrs = -10) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  PMobj <- new("PMobj")
  PMobj@Name <- "Female Spawning Biomass relative to SBMSY (Last 10 yrs)"
  if (Ref != 1) {
    PMobj@Caption <- paste0("Long-term Healthy Biomass: \nProb. SB > ", Ref, " SB_MSY (Years ", 
                            Yrs[1], " - ", Yrs[2], ")")
  } else {
    PMobj@Caption <- paste0("Long-term Healthy Biomass: Prob. SB > SBMSY (Years ", Yrs[1], 
                            " - ", Yrs[2], ")")
  }
  PMobj@Ref <- Ref

  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
    
  PMobj@Stat <- matrix(NA, nrow=MMSE@nsim, ncol=MMSE@nMPs)
  PMobj@Prob <- apply(MMSE@B_BMSY[,1, ,Yrs[1]:Yrs[2]] > PMobj@Ref, 1:2, mean) 
  PMobj@Mean <- calcMean(PMobj@Prob)
  PMobj@MPs <- MPnames
  PMobj
}
class(P80Last10) <- 'PM'

2.5.1.7 Performance Objective: Probability Average Annual Variability in Yield < 20%

# AAVY with new Caption
VarY <- function (MMSE = NULL, Ref = 0.2, Yrs = NULL) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  
  PMobj <- new("PMobj")
  PMobj@Name <- paste0("Average Annual Variability in Yield (Years ", 
                       Yrs[1], "-", Yrs[2], ")")
  PMobj@Caption <- paste0("Variation in Yield:\nProb. AAVY < ", Ref * 100, "% (Years ", Yrs[1], "-", Yrs[2], ")")
  
  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
  
  Catch <- apply(MMSE@C[,,,,Yrs[1]:Yrs[2], drop=FALSE], c(1,4,5), sum) # total catch
  y1 <- Yrs[1]:(Yrs[2] - 1)
  y2 <- (Yrs[1] + 1):Yrs[2]
  AAVY <- apply(((((Catch[, , y1] - Catch[, , y2])/Catch[, , y2])^2)^0.5), 
                c(1, 2), mean)
  PMobj@Stat <- AAVY
  PMobj@Ref <- Ref
  PMobj@Prob <- calcProb(PMobj@Stat < Ref, MMSE)
  PMobj@Mean <- calcMean(PMobj@Prob)
  PMobj@MPs <- MPnames
  PMobj
}
class(VarY) <- 'PM'

2.5.1.8 Performance Objective: Probability Average Annual Variability in Effort< 20%

VarE <- function (MMSE = NULL, Ref = 0.2, Yrs = NULL) {
  if (class(MMSE)!='MMSE') stop('Require object of class `MMSE`')
  Yrs <- ChkYrs(Yrs, MMSE)
  PMobj <- new("PMobj")
  PMobj@Name <- paste0("Average Annual Variability in Effort (Years ", 
                       Yrs[1], "-", Yrs[2], ")")
  PMobj@Caption <- paste0("Variation in Effort:\nProb. AAVE < ", Ref * 100, "% (Years ", 
                          Yrs[1], "-", Yrs[2], ")")
  y1<- Yrs[1]:(Yrs[2]-1) # year index
  y2<-(Yrs[1]+1):Yrs[2] 
  
  Effort <- apply(MMSE@Effort[,,,,Yrs[1]:Yrs[2], drop=FALSE], c(1,4,5), mean)
  AAVE <- apply(((((Effort[, , y1] - Effort[, , y2])/Effort[, , y2])^2)^0.5), c(1, 2), mean)
  
  # MP names by stock/fleet (female MP first, then male)
  MPs <- MMSE@MPs
  if (class(MPs) == 'list') {
    # MPs by stock
    if (class(MPs[[1]][[1]]) == 'list') {
      # by fleet as well
      stop("Function needs updating")
    } else {
      stocks <- 1:MMSE@nstocks
      df <- data.frame(MPs)
      colnames(df) <- stocks
      MPnames <- rep('', nrow(df))
      for (i in 1:nrow(df)) {
        mps <- as.character(df[i,] )
        MPnames[i] <- paste(unique(mps), collapse="-")
      }
    }
  }
  
  PMobj@Stat <- AAVE
  PMobj@Ref <- Ref
  PMobj@Prob <- calcProb(PMobj@Stat < Ref, MMSE)  # probability AAVE < 0.2 
  
  PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
  PMobj@MPs <- MPnames
  PMobj
}
class(VarE) <- 'PM'

2.5.1.9 Trade-Off Plot

The TradePlot function can then be used to plot the custom performance metrics (note the legend currently does not work for results from MMSE objects due to different MPs that are applied to stocks/fleets):

Performance Limits:

TradePlot(RedCrab_MSE, 'P50mgt', 'P50Last10', Lims=c(0.8,0.8))

Performance Objectives:

TradePlot(RedCrab_MSE, 
          'LTAvY', 'STAvY',
          'P80Last10', 'LTAvY',
          'VarY', 'VarE',
          Lims=0)

## Recycling limits

Custom plotting functions specific to the application will need to be developed if stock or fleet specific performance metrics are required.

3 Lobster Case Study

The lobster case study has two stocks (female and male) and two fleets (commercial and recreational). Some of the explanatory material in this section is repeated from the previous section so that the description of the case study can be self contained.

3.1 Set up `MOM` Object

3.1.1 Import Stocks and create Stock Objects

 # directory where the OM Excel files are located
OMdir <- 'G:/My Drive/1_PROJECTS/CDFW_Multi_Fleet/LobsterDLM'

# Import the Female and Male LobsterStocks
# (ignore any warnings about missing slots)
LobsterF_Stock <- XL2Stock(file.path(OMdir, 'Lobster_Stock_Female.xlsx'), msg=FALSE)
LobsterM_Stock <- XL2Stock(file.path(OMdir, 'Lobster_Stock_Male.xlsx'), msg=FALSE)

# a list of stock objects (female stock in first position)
Stocks <- list(LobsterF_Stock, LobsterM_Stock)

3.1.2 Import Fleets and create Fleet Objects

Two fleets, commercial and recreational, are imported:

Com_Fleet <- XL2Fleet(file.path(OMdir, 'Lobster_Fleet_Commercial.xlsx'), msg=FALSE)

## Warning in XL2Fleet(file.path(OMdir, "Lobster_Fleet_Commercial.xlsx"), msg
## = FALSE): SelYears AbsSelYears L5Lower L5Upper LFSLower LFSUpper VmaxLower
## VmaxUpper are not valid slots in object class Fleet

Rec_Fleet <- XL2Fleet(file.path(OMdir, 'Lobster_Fleet_Recreational.xlsx'), msg=FALSE)

## Warning in XL2Fleet(file.path(OMdir, "Lobster_Fleet_Recreational.xlsx"), :
## SelYears AbsSelYears L5Lower L5Upper LFSLower LFSUpper VmaxLower VmaxUpper are
## not valid slots in object class Fleet

Note the warning message about invalid slots SelYears, AbsSelYears, etc. These slots have no been removed from the Fleet object. Time-varying selectivity is passed in to the model with cpars feature.

The time-varying selectivity parameters are the same for both the recreational and commercial fleets, so one selectivity matrix is created here and applied to both stocks. Note that these values are interpreted here as time-varying selectivity rather than retention (i.e., after 1955, the length at 5% capture was ~ 75 - 83 mm):

nsim <- 10 # using low number of simulations for testing 
proyears <- 30 # 30-year projection period
nyears <- Com_Fleet@nyears 

L5_y <- matrix(NA, nrow=nsim, ncol=proyears+nyears) # empty matrices
LFS_y <- matrix(NA, nrow=nsim, ncol=proyears+nyears)
Vmaxlen_y <- matrix(NA, nrow=nsim, ncol=proyears+nyears)
  
# from SelYears in Fleet tab
# Yr index where selectivity changes
SLdata <- readxl::read_xlsx(file.path(OMdir, 'Lobster_Fleet_Commercial.xlsx'))

## New names:
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * `` -> ...6
## * ...

sel_years <- as.numeric(SLdata[which(SLdata$Slot == "SelYears"), 2:5])

L5Lower <- as.numeric(SLdata[which(SLdata$Slot == "L5Lower"), 2:5])
L5Upper<- as.numeric(SLdata[which(SLdata$Slot == "L5Upper"), 2:5])

LFSLower <- as.numeric(SLdata[which(SLdata$Slot == "LFSLower"), 2:5])
LFSUpper<- as.numeric(SLdata[which(SLdata$Slot == "LFSUpper"), 2:5])

VmaxLower <- as.numeric(SLdata[which(SLdata$Slot == "VmaxLower"), 2:5])
VmaxUpper <- as.numeric(SLdata[which(SLdata$Slot == "VmaxUpper"), 2:5])

for (i in seq_along(sel_years)) {
  yr.ind <- sel_years[i]
  yr.ind2 <- sel_years[i+1]
  if (is.na(yr.ind2)) yr.ind2 <- nyears+proyears
  
  l5s  <- runif(nsim*10, L5Lower[i], L5Upper[i]) # sample L5
  lfs <- runif(nsim*10, LFSLower[i], LFSUpper[i]) # sample LFS
    chk <- lfs>l5s # make sure LFS > L5
  l5s <- l5s[chk][1:nsim]
  lfs <- lfs[chk][1:nsim]
  
  L5_y[,yr.ind:yr.ind2] <- l5s
  LFS_y[,yr.ind:yr.ind2] <- lfs
  Vmaxlen_y[,yr.ind:yr.ind2] <- runif(nsim, VmaxLower[i], VmaxUpper[i]) 
}

# make list
temp <- list()
temp$L5_y <- L5_y
temp$LFS_y <- LFS_y
temp$Vmaxlen_y <- Vmaxlen_y

# add to cpars
# female
cpars <- list()
cpars[[1]] <- list() 
cpars[[1]][[1]] <- temp # commercial
cpars[[1]][[2]] <- temp # recreational

 # male
cpars[[2]] <- list()
cpars[[2]][[1]] <- temp # commercial
cpars[[2]][[2]] <- temp # recreational

Note also that historical MPAs are not currently implemented in the multi-stock model (this will be added at a later date):

Com_Fleet@MPA <- FALSE 
Rec_Fleet@MPA <- FALSE

A Fleets list is created, where the two fleets are applied to both the female and male stocks:

Fleets <- list(
  list(Com_Fleet, Rec_Fleet), # female stock
  list(Com_Fleet, Rec_Fleet) # male stock
)

3.1.3 Set up 2-sex model

SexPars <- list()
# Stock (Row) spawn from SSB contributed by Stock (Column)
SexPars$SSBfrom <- matrix(c(1,0,
                            1,0),
                          byrow=TRUE,nrow=2)

SexPars$SSBfrom

##      [,1] [,2]
## [1,]    1    0
## [2,]    1    0

SexPars$SSBfrom is a 2x2 matrix specifying recruitment to females (row 1) and males (row 2) comes from female spawning biomass (column 1).

Note that it is important that the SexPars$SSBfrom matches the order of the stocks in Stocks (ie., female in first position in list Stocks, and first column in SexPars$SSBfrom.

3.1.4 Allocate Catches between Fleets

The CatchFrac argument is used to allocated the catch fraction among the fleets in the MOM object.

The data files from CDFW suggest that recent recreational catches are about 55% of the commercial catch, or 36% of the total catch. Assuming that the recent catch is made up of equal number of male and female lobster, the CatchFrac argument is populated as:

FemaleCatch <- c(0.64, 0.36) # commercial and recreational fleets
MaleCatch <- c(0.64, 0.36)

CatchFrac <- list(
   matrix(rep(FemaleCatch, each=nsim), nrow=nsim),
   matrix(rep(MaleCatch, each=nsim), nrow=nsim)
)

3.1.5 Create `MOM` Object

For simplicity, we will use the generic observation object for both fleets and both stocks:

ObsList <- list(
  list(Generic_Obs, Generic_Obs), # Obs for female (2 fleets)
  list(Generic_Obs, Generic_Obs) # Obs for male (the same)
)

And no implementation error on either fleets or stocks:

ImpList <- list(
  list(Perfect_Imp, Perfect_Imp),
  list(Perfect_Imp, Perfect_Imp)
)

The MOM object is created using the new function:

MOM <- new("MOM",
           Stocks = Stocks,
           Fleets = Fleets,
           Obs = ObsList,
           Imps = ImpList,
           CatchFrac = CatchFrac,
           SexPars=SexPars,
           cpars=cpars,
           nsim = nsim,
           interval=1,
           proyears=proyears)

3.2 Run the Historical Simulations

3.2.1 Simulate Historical Fishery

LobsterHist <- SimulateMOM(MOM, parallel = TRUE)

To avoid re-running the historical simulations each time, they can be saved to disk:

saveRDS(LobsterHist, file.path(OMdir, 'LobsterHist.rda'))

3.3 Define Management Procedures

The management options that have been proposed by CDFW for the lobster fishery are static size limits, male only fishery (release of all female lobster), and effort reductions.

The recreational lobster fishery is managed with daily bag limits (7) and a minimum size limit. Direct modeling of the bag limits is not currently supported in the MSE framework and will require further development (and possibly additional data) to implement.

For these scenarios, fishing effort for both the recreational and commercial fleets is assumed to remain at the current (2017) level, unless changed by a management procedure.

3.3.1 Current Management

The curE management procedure assumes that the current size limit remains the same for the two fleets, and effort remains at the 2017 level.

3.3.2 Fleet-Specific Size Limits

The ComIncSL MP increases the size of retention of the commercial fleet to an arbitrary 110 mm (for example only) while keeping the recreational size limit at the current value of 82 mm. Effort remains at the current (2017) level for both fleets.

ComIncSL <- function(x, Data, ...) {
  rec <- new('Rec')
  rec@LR5 <- 109
  rec@LFR <- 110
  rec
}
class(ComIncSL) <- 'MP'

3.3.3 Commercial Effort Limit

The ComEL MP reduces fishing effort of the commercial fleet to an arbitrary 75% of the most recent level. Size limits and effort for the recreational fleet remains unchanged.

ComEL <- function(x, Data, ...) {
  rec <- new("Rec")
  rec@Effort <- 0.75
  rec
}
class(ComEL) <- 'MP'

3.3.4 Commercial TAC

The ComTAC MP sets an annual TAC for the commercial fleet. The current effort limit is also applied, so that fishing effort for the commercial fleet cannot increase higher than the 2017 level. The effort limit also remains for the recreational fleet, and size limits are unchanged for both fleets.

ComTAC <- function(x, Data, ...) {
  rec <- new('Rec')
  rec@Effort <- 1 # effort capped at current level
    yrs <- min(Data@Year):(Data@Year[Data@Year==Data@LHYear[1]])
  yr.ind <- match(yrs, Data@Year)
  histCatch <- Data@Cat[x, yr.ind]
  meanC <- mean(histCatch, na.rm = T)
  rec@TAC <- meanC # TAC set to historical average catch
  rec
}
class(ComTAC) <- 'MP'

3.3.5 Male Only Fishery

Finally, the MOnly MP sets retention to zero for female lobster for both the commercial and recreational fleets, and keeps fishing effort fixed to the current level:

MOnly <- function(x, Data, ...) {
  rec <- new('Rec') 
  # set length-at-5% retention to maximum size class (ie nothing is retained)
  rec@LR5 <- max(Data@CAL_bins) 
  rec@LFR <- rec@LR5 + 1
  rec@Effort <- 1
  rec
}
class(MOnly) <- 'MP'

3.4 Run Forward Projections

The forward projections are done with the ProjectMOM function. This function requires the historical simulations returned by SimulateMOM (class multiHist) and a list of management procedures.

The management procedures are specified in byFleet mode; MPs applied to each stock and fleet separately. Notice that the recreational fleet assumes current size limit and effort fixed at current level (curE) in all cases except the last MP, where there is no retention of female lobster in the recreational fleet.

MP_by_fleet <- list(
    # MPs for female stock
    list(
      c('curE', 'ComIncSL', 'ComEL', 'ComTAC','MOnly'), # Comm Fleet
      c('curE', 'curE', 'curE', 'curE', 'MOnly') # Rec fleet
      ),
     # MPs for male stock
    list(
      c('curE', 'ComIncSL', 'ComEL', 'ComTAC', 'curE'), # Comm Fleet
      c('curE', 'curE', 'curE', 'curE', 'curE') # Rec fleet
    )
)

# load the saved historical simulations
LobsterHist <- readRDS(file.path(OMdir, 'LobsterHist.rda')) 

# run the forward projections
Lobster_MSE <- ProjectMOM(LobsterHist, MPs=MP_by_fleet)

# save the MSE run
saveRDS(Lobster_MSE, file.path(OMdir, 'Lobster_MSE.rda'))

3.5 Examine Results

Similar to the rock crab case study, the results are plotted with the generic plot function, which creates a series of projection plots, showing a) Spawning biomass $(\text{SB})$ relative to $\text{SB}_\text{MSY}$ for each stock and each MP, b) fishing mortality $(F)$ relative to $F_\text{MSY}$ for each stock and MP, c) total yield relative to the first projection year for each stock and MP, and d) average yield (absolute) for each stock and fleet (and MP):

plot(Lobster_MSE, curyr=2017)

The performance limits and objectives can be plotted the same way as done for the rock crab case study, using the TradePlot function and the custom PM functions:

Performance Limits:

TradePlot(Lobster_MSE, 'P50mgt', 'P50Last10', Lims=c(0.8,0.8))

Performance Objectives:

TradePlot(Lobster_MSE, 
          'LTAvY', 'STAvY',
          'P80Last10', 'LTAvY',
          'VarY', 'VarE',
          Lims=0)

## Recycling limits

4 Halibut Case Study

The halibut case study is different from the previous two examples. In this case the MOM object is created by directly importing the output of the Stock Synthesis 3 stock assessment.

This case study is based on the North California stock. The analysis can be repeated for the South stock using the SS3 output for that stock.

Some of the explanatory material in this section is repeated from the previous section so that the description of the case study can be self contained.

4.1 Import the `MOM` Object

The MOM object is imported directly from the SS3 output usingthe SS2OM function:

 # directory where the SS files are located
SSdir <- 'G:/My Drive/1_PROJECTS/CDFW_Multi_Fleet/Halibut/SS3/North'

MOM <- SS2MOM(SSdir, 
              Obs=Generic_Obs,
              Imp=Perfect_Imp,
              nsim=10) # 10 simulations for quick testing

## Registered S3 method overwritten by 'gdata':
##   method         from  
##   reorder.factor gplots

## -- Using function SS_output of package r4ss version 1.40.0 to extract 
## data from SS file structure --

## Reading directory: G:/My 
## Drive/1_PROJECTS/CDFW_Multi_Fleet/Halibut/SS3/North

## Some stats skipped because the .par file not found.

## -- End of r4ss operations --

## 2 -sex and 5 -fleet model detected.

The MOM object has been populated with the information for 2 stocks (female and male) and 5 fishing fleets:

length(MOM@Stocks)

## [1] 2

sapply(1:length(MOM@Stocks), function(x) slot(MOM@Stocks[[x]], 'Name'))

## [1] "Female" "Male"

length(MOM@Fleets)

## [1] 2

sapply(1:length(MOM@Fleets[[1]]), function(x) slot(MOM@Fleets[[1]][[x]], 'Name'))

## [1] "Bottom.Trawl" "Gillnet"      "CPFV"         "Other.Rec"    "Comm.HL"

The SexPars and CatchFrac slots are also populated:

MOM@SexPars

## $SSBfrom
##      [,1] [,2]
## [1,]    1    0
## [2,]    1    0

# Catch fraction for females 
MOM@CatchFrac[[1]]

##            [,1]        [,2]      [,3]      [,4]      [,5]
##  [1,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [2,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [3,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [4,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [5,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [6,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [7,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [8,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
##  [9,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551
## [10,] 0.2004257 0.000161206 0.1670437 0.3960143 0.2363551

4.2 Run the Historical Simulations

4.2.1 Simulate Historical Fishery

HalibutHist <- SimulateMOM(MOM, parallel = TRUE)

To avoid re-running the historical simulations each time, they can be saved to disk:

saveRDS(HalibutHist, file.path(SSdir, 'HalibutHist.rda'))

4.3 Define Management Procedures

A range of management options have been proposed by CDFW for the halibut fishery, including fleet-specific size limits, total catch and total effort limits for the commercial fleets, and modifications to the selectivity of the commercial fishing gear.

The recreational halibut fishery is managed with daily bag limits (7) and a minimum size limit. Direct modeling of the bag limits is not currently supported in the MSE framework and will require further development (and possibly additional data) to implement.

For these scenarios, fishing effort for both the recreational and commercial fleets is assumed to remain at the current (2019) level, unless changed by a management procedure.

The following example MPs have been designed based on the list of management options provided by CDFW for the halibut fishery. Not all of the options suggested by CDFW were able to be modeled at this time (e.g., bag limits). The MPs will be developed further based on feedback from CDFW.

4.3.1 Current Management

The curE management procedure assumes that the current size limit remains the same for the two fleets, and effort remains at the 2019 level.

4.3.2 Fleet-Specific Size Limits & TAC

These MPs specify an increase to the size of retention for the recreational and commercial hook & line fleets, remove the size limits for the trawl and gillnet fleets, and set a TAC using the index of abundance from the trawl logbook CPUE data.

Not that the model currently uses simulated data. It can be updated in the future to include the real trawl CPUE logbook data.

# MP to increase size limit by 15%
IncSL <- function(x, Data, ...) {
  rec <- new("Rec")
  rec@LR5 <- Data@OM$LR5[x] * 1.15 # increase real LR5 by 15%
  rec@LFR <- Data@OM$LFR[x] * 1.15 # increase by 15%
  rec
}
class(IncSL) <- 'MP'

# MP to remove size limit and set TAC 
noSL_TAC <- function(x, Data, reps) {
  rec <- new("Rec")
  rec@LR5 <- Data@OM$L5[x] # set retention to selectivity
  rec@LFR <- Data@OM$LFR[x]
  rec@Rmaxlen <- Data@OM$Rmaxlen[x]
  
  # Set TAC using `ITarget1`
  runItarget <- Itarget_(x, Data, reps, plot=FALSE, 
                         yrsmth=5, xx=0, Imulti=1.5)
  
  rec@TAC <-  runItarget$TAC
  rec
}
class(noSL_TAC) <- 'MP'

4.3.3 Total Allowable Effort Limits

The ItargetE1 MP maintains all current selectivity & retention patterns and sets a TAE using the index of total abundance.

4.4 Run Forward Projections

The forward projections are done with the ProjectMOM function. This function requires the historical simulations returned by SimulateMOM (class multiHist) and a list of management procedures.

The management procedures are specified in byfleet mode; MPs applied to each stock and fleet separately. Because there are no MPs with sex-specific management, the MPs are repeated for the male stock.

The fleet names are:

Fnames <- sapply(1:5, function(x)
  slot(MOM@Fleets[[1]][[x]], 'Name'))
Fnames

## [1] "Bottom.Trawl" "Gillnet"      "CPFV"         "Other.Rec"    "Comm.HL"

The MPs are specified as follows:

Current effort and size limits for all fleets
Higher size limit for recreational and commercial hook and line, and remove size limit and apply TAC for the trawl and gillnet fleets. Effort for the hook and line and recreational fleets remains fixed at current level.
A TAE for trawl and gillnet fleets. Other fleets remain at current effort.

MP_by_fleet <- list(
    # MPs for female stock
    list(
      c('curE', 'noSL_TAC', 'ItargetE1'), # Bottom trawl
      c('curE', 'noSL_TAC', 'ItargetE1'), # Gillnet
      c('curE', 'IncSL', 'curE'), # CPFV
      c('curE', 'IncSL', 'curE'), # Other.Rec
      c('curE', 'IncSL', 'curE') # Comm.HL
      ),
     # MPs repeated for male stock
    list(
      c('curE', 'noSL_TAC', 'ItargetE1'), # Bottom trawl
      c('curE', 'noSL_TAC', 'ItargetE1'), # Gillnet
      c('curE', 'IncSL', 'curE'), # CPFV
      c('curE', 'IncSL', 'curE'), # Other.Rec
      c('curE', 'IncSL', 'curE') # Comm.HL
    )
)

# load the saved historical simulations
HalibutHist <- readRDS(file.path(SSdir, 'HalibutHist.rda'))

# Run forward projections
Halibut_MSE <- ProjectMOM(HalibutHist, MPs=MP_by_fleet)

# Save MSE
saveRDS(Halibut_MSE, file.path(SSdir, 'Halibut_MSE.rda'))

4.5 Examine Results

Similar to the previous two case studies, the results are plotted with the generic plot function, which creates a series of projection plots, showing a) Spawning biomass $(\text{SB})$ relative to $\text{SB}_\text{MSY}$ for each stock and each MP, b) fishing mortality $(F)$ relative to $F_\text{MSY}$ for each stock and MP, c) total yield relative to the first projection year for each stock and MP, and d) average yield (absolute) for each stock and fleet (and MP):

plot(Halibut_MSE, curyr=2019)

The performance limits and objectives can be plotted the same way as done for the rock crab case study, using the TradePlot function and the custom PM functions:

Performance Limits:

TradePlot(Halibut_MSE, 'P50mgt', 'P50Last10', Lims=c(0.8,0.8))

Performance Objectives:

TradePlot(Halibut_MSE, 
          'LTAvY', 'STAvY',
          'P80Last10', 'LTAvY',
          'VarY', 'VarE',
          Lims=0)

## Recycling limits

Developing the MSEtool package for Management Strategy Evaluation with Multiple Stocks and Fleets

Final Report for the California Department of Fish and Wildlife Multi-Stock Multi-Fleet Case-studies

Adrian Hordyk

adrian@bluematterscience.com

30 October, 2020

1 Introduction

1.1 New Package Structure

1.2 Installing and Loading the New Packages

2 Rock Crab Case Study

2.1 Set up MOM Object

2.1.1 Import Stocks and create Stock Objects

2.1.2 Import Fleet and create Fleet Object

2.1.3 Set up 2-sex model

2.1.4 Create MOM Object

2.2 Run the Historical Simulations

2.2.1 Simulate Historical Fishery

2.3 Define Management Procedures

2.3.1 Current Size Limit

2.3.2 Male Only Fishery & Current Size Limit

2.3.3 Increased Size Limit

2.3.4 Effort Limit

2.4 Run Forward Projections

2.5 Examine Results

2.5.1 Performance Metrics and Trade-Off Plots

2.5.1.1 Performance Limit: Probability Female SB > 0.5SBMSY after 1 Mean Generation Time

2.5.1.2 Performance Limit: Probability Female SB > 0.5SBMSY in Last 10 Years

2.5.1.3 Performance Object: Average Long-Term Yield (Last 10 Years)

2.5.1.4 Performance Objective: Average Short-Term Yield (First 10 Years)

2.5.1.5 Performance Objective: Average Yield (all but first 10 years)

2.5.1.6 Performance Objective: Probability Female SB > 0.8SBMSY in Last 10 Years

2.5.1.7 Performance Objective: Probability Average Annual Variability in Yield < 20%

2.5.1.8 Performance Objective: Probability Average Annual Variability in Effort< 20%

2.5.1.9 Trade-Off Plot

3 Lobster Case Study

3.1 Set up MOM Object

3.1.1 Import Stocks and create Stock Objects

3.1.2 Import Fleets and create Fleet Objects

3.1.3 Set up 2-sex model

3.1.4 Allocate Catches between Fleets

3.1.5 Create MOM Object

3.2 Run the Historical Simulations

3.2.1 Simulate Historical Fishery

3.3 Define Management Procedures

3.3.1 Current Management

3.3.2 Fleet-Specific Size Limits

3.3.3 Commercial Effort Limit

3.3.4 Commercial TAC

3.3.5 Male Only Fishery

3.4 Run Forward Projections

3.5 Examine Results

4 Halibut Case Study

4.1 Import the MOM Object

4.2 Run the Historical Simulations

4.2.1 Simulate Historical Fishery

4.3 Define Management Procedures

4.3.1 Current Management

4.3.2 Fleet-Specific Size Limits & TAC

4.3.3 Total Allowable Effort Limits

4.4 Run Forward Projections

4.5 Examine Results

2.1 Set up `MOM` Object

2.1.4 Create `MOM` Object

3.1 Set up `MOM` Object

3.1.5 Create `MOM` Object

4.1 Import the `MOM` Object