In tcarruth/MSEtool: Management Strategy Evaluation Toolkit
library(MSEtool)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
knitr::opts_chunk$set(dpi=85)
options(width = 650)
Introduction
Management Strategy Evaluation can be a complex business even for single-stock and single-fleet operating models. Population dynamics, fleet dynamics, observation processes (e.g. catch over reporting), implementation error (e.g. catch-limit overages) and appropriate management procedures (MPs) all need to be specified.
In most fisheries settings, single-stock, single-fleet operating models are sufficient. Unless management advice is explicitly given on a fleet-by-fleet basis (e.g. two TACs each specific to recreational and commercial fisheries) then the fleet dynamics can be pooled into an aggregate fishery. Unless (A) the dynamics of one stock interact with another stock (e.g. a predator and prey) or (B) stocks with varying dynamics are managed as a complex (e.g. one TAC for all sculpins as in Alaska) then MPs can be tested for one individual stock at a time.
There are situations when managers wish to evaluate the robustness of a proposed management plan among fleets and/or stocks, acknowledging that there are interactions between these. To this end, MSEtool includes the function multiMSE(), a version of DLMtool::runMSE() that conducts MSE for multiple fleets/stocks and allows users to prescribe interactions between stocks using R models (including sex-specific population dyanmics and hermaphroditism).
To accomplish this, you will need to become familiar with a new class of operating model object MOM (multi operating model), a new class of MSE object MMSE (multi MSE), and optionally, a new class of management procedure MMP (multi MP).
Beyond learning about these new objects and functions, multiMSE() is extremely simple to use and hence misuse. I can't stress enough how important it is to have a clear justification for why a multi-stock and/or multi-fleet MSE is necessary for your particular management setting (it may not be). To code multiMSE() it was necessary to identify all of the various use-cases and it became clear that there are a very large number each with their own specific set of requirements and assumptions.
In this guide, I will attempt to clarify these various use cases and be explicit about the assumptions entailed.
Before going crazy with many fleet, many stock models, it might be good to temper your expectations. You can do an awful lot with multiMSE() but numerically solving the multi-stock, multi-fleet equations while maintaining MICE relationships is computationally costly. This is only reflected in the initialization phase of the MSE - which can take some time. A rough rule of thumb is that the minutes to initialize = (nsimulations * nstocks * nfleets) divided by (10 * ncores). So for a 100-simulation MSE with 4 stocks, 5 fleets and parallel computing set up over 4 cores you are talking about 2000 / 40 = 50 minute initialization. Faster solutions are on the way for non-MICE models that don't require R optimizaiton due to R MICE functions.
Caution: you are now wandering into a relatively unestablished and unchartered territory of MSE - good luck!
Remember you can always fire off an email to me if you have questions about multiMSE: t.carruthers@oceans.ubc.ca
For general help, Adrian Hordyk has produced a comprehensive guide for DLMtool including cheat sheets, and Quang Huynh has produced additional user guide vignettes for MSEtool that can be obtained using the standard R help functions:
DLMtool::userguide()
DLMtool::cheatsheets()
MSEtool::userguide()
The multi operating model object (class MOM)
The MOM object is not terribly different from a standard DLMtool/MSEtool OM object.
library(DLMtool) # load DLMtool library
DLMextra() # download lots of cool DLMtool/MSEtool objects from GitHub
library(DLMextra) # load all those objects
library(MSEtool) # load MSEtool (data rich MPs, MOM object class etc)
setup() # set up a cluster for parallel processing
slotNames('MOM') # what slots are there in an MOM object
class?MOM # check out the help file on MOM objects
The big difference between objects of class OM in DLMtool and those of class MOM is that instead of including a slot for each parameter of the stock (Stock), fleet (Fleet), observation (Obs) and implementation (Imp) models for OM objects, MOM objects include these in four slots that are lists of entire Stock, Fleet, Obs and Imp objects - after all it's multi-fleet and/or multi stock, right?!
Stocks, Fleets, Obs and Imps
There are four new slots in an MOM object (compared with an OM object): MOM@Stocks, MOM@Fleets, MOM@Obs, and MOM@Imps.
MOM@Stocks is a list of stock objects (nstocks long)MOM@Fleets is a hierarchical list of fleet objects (nstocks with fleets nested in stocks)MOM@Obs is a hierarchical list of observation error objects (nstocks with fleets nested in stocks)MOM@Imps is a hierarchical list of implementation error objects (nstocks with fleets nested in stocks)
Ideally you'd get your various dynamics from conditioned operating models either using the DLMtool:: StochasticSRA() approach or using the various MSEtool helper functions that read in stock assessment outputs from, for example, Stock Synthesis (SS2OM) or iSCAM (iSCAM2OM) assessments. However, for simplicity I'll demonstrate how to do multiMSE with some pre-existing Stock, Fleet, Obs, and Imp objects that were loaded with library(DLMextra):
avail('Stock') # what Stock objects are loaded?
avail('Fleet') # what Fleet objects are loaded?
avail('Obs') # what Obs objects are loaded?
avail('Imp') # what Imp objects are loaded?
Let's imagine that we lose our minds and decide to create a model that simulates interactions between bluefin tuna and herring (we think bluefin eat herring). We also want to simulate three fleets, a bluefin fleet, a herring fleet and a mixed fleet that somehow, against all the odds catches both bluefin and herring! About as hypothetical as it gets I know, but that is intentional - this is just a demonstration!
First, we start off by creating the list of stock objects, this is pretty easy:
# Stock 1 Stock 2
Stocks <- list(Bluefin_tuna, Herring)
A key constraint here is that the Stock objects have the same number of areas (e.g. Bluefin_tuna@nareas).
We now want to create our fleets for these stocks. Our bluefin-only fleet has stable effort (premade fleet type Generic_FlatE) our herring-only fleet has increasing effort (premade fleet type Generic_IncE) and our mythical mixed fishery has decreasing effort (premade fleet type Generic_DecE).
This is the tricky bit.
Nominally, the same fleet occupies the same position for each stock. For example, 'longline' is fleet 1 for both bluefin and herring. However, the fleet interactions with the stocks can (and often will) be different - e.g. the trend in exploitation rate of the longline fleet may be stable for bluefin but decreasing for herring (in this admittedly bizarre case). That is why multiMSE() requires the same number of fleets for each stock but you have to be able to specify unique fleet dynamics for each stock x fleet. For example longline_stable_exploitation could be Fleet 1 for bluefin and longline_decreasing_exploitation could be Fleet 1 for Herring - it is the same fleet but interacts differently with the two stocks.
A rare exception is if you have one fleet catching each stock. Then you can just have one fleet position for each stock and put the stock-specific fleet in position Fleet 1 for all stocks. Confused? I don't blame you! But this should get clearer as you work through this guide.
To specify exploitation dynamics across stocks and fleets you create a hierarchical list structure (Fleets nested in Stocks):
# Fleet 1 Fleet 2 Fleet 3
Fleets_both_Stocks <- list(Generic_FlatE, Generic_IncE, Generic_DecE) # for the sake of this exercise we assume identical fleet dynamics for both stocks
Fleets <- list( Fleets_both_Stocks, # Bluefin (position 1 in Stocks)
Fleets_both_Stocks) # Herring (position 2 in Stocks)
Since each of these Fleets is going to generate fishery data for each stock, they each need an observation model by stock and fleet. Just like the Fleets list, this is hierarchical - Fleets nested in Stocks. We are going to assume that we get good data (Precise and Unbiased) for bluefin and bad data (Imprecise and Biased) for herring.
# Fleet 1 Fleet 2 Fleet 3
Bluefin_obs <- list(Precise_Unbiased, Precise_Unbiased, Precise_Unbiased)
Herring_obs <- list(Imprecise_Biased, Imprecise_Biased, Imprecise_Biased)
Obs <- list( Bluefin_obs, # Bluefin (position 1 in Stocks)
Herring_obs) # Herring (position 2 in Stocks)
Lastly, we need to specify implementation error models for these, we are going to assume this is the same for both stocks and is fleet specific: perfect for the first 2 fleets and includes overages for the mixed fleet (fleet 3).
# Fleet 1 Fleet 2 Fleet 3
Imp_both_Stocks <- list(Perfect_Imp, Perfect_Imp, Overages)
Imps <- list( Imp_both_Stocks, # Bluefin (position 1 in Stocks)
Imp_both_Stocks) # Herring (position 2 in Stocks)
Each of the stocks already has a prespecified range of depletion:
Bluefin_tuna@D
Herring@D
But how do we attribute the magnitude of historical exploitation among the various fleets?
Catch fractions by fleet for each stock
To allow multiMSE() to run we need an extra Thing 1, and should have an additional Thing 2 after that.
Thing 1 is the current (most recent year) catch fractions among fleets for each stock. This information is required to calculate the relative catchabilities among the fleets.
Catch fractions are a list (nstock long) of matrices. Suppose we want 80% of current bluefin catches to come from fleet 1, 1% from fleet 2 and 19% from fleet 3, and also 60% of herring catches to come from Fleet 2, 5% from fleet 1 and the remaining 35% of herring catches from Fleet 3, this is how you would specify that:
nsim = 4
# Fleet1 Fleet2 Fleet3
bluefin_catchfrac<-c(0.8, 0.01, 0.19)
herring_catchfrac<-c(0.05, 0.6, 0.35)
CatchFrac <- list(
matrix( rep(bluefin_catchfrac, each=nsim), nrow=nsim),
matrix( rep(herring_catchfrac, each=nsim), nrow=nsim)
)
CatchFrac[[1]]
The observant types among you will have noticed that we now need to specify the number of MSE simulations we are going to be doing (nsim). For now this is just 12 for demonstration purposes, but ultimately this would be >148 for a real MSE.
The reason we have done this is to allow for uncertainty in the catch fractions - the user could make catch fractions vary among fleets by simulation. For now, these are the same among simulations, which is the simplest thing, albeit a bit boring.
Specifying inter-stock relationships
You'll remember that I said there was a Thing 2 that we should add to our multi operating model. Ask yourself why we are modelling bluefin and herring together?
Because we want to look at an interaction between these stocks. If we didn't, we would just model these stocks individually, and there would be only one stock in the Stocks slot and only one list of fleets in the Fleets slot (the remaining reason for modelling many stock without interactions is if they are managed by a single MP as a stock complex - we'll come to this in a bit).
If we didn't expect to set fleet-specific advice AND we weren't modelling bluefin and herring interactions, we should just be using the regular runMSE() from DLMtool and MSEtool and calculate the aggregate fleet dynamics (aggregate size vulnerability and trend in apical fishing mortality rate, for example).
But we're demonstrating multiMSE() here, and are therefore going to specify at least one relationship between the stocks (noting that you may wish to evaluate MPs that provide advice for a stock complex of multiple stock objects with no interactions - more below). The slot MOM@Rel is where you put these relationships as a list of R models with quite specific formatting.
Since we think bluefin are predators of herring, we are going to propose a relationship between bluefin abundance and the natural mortality rate of herring. Normally we would have derived this relationship empirically from data (in this case perhaps herring tagging survival estimates regressed on bluefin abundance estimates from a stock assessment).
However for the purposes of this demonstration we're just going to cook something up out of thin air: nominally we are going to assume that herring have a natural mortality rate of 0.2 when there are no bluefin, and 0.4 when bluefin tuna are at unfished levels. To invent this relationship we are going to need to calculate bluefin unfished biomass from the operating model, make some fake data and fit an R model.
Here goes:
ages <- 1:60
bf_M <- mean(Bluefin_tuna@M) # natural mortality rate
bf_Linf <- mean(Bluefin_tuna@Linf) # asymptotic length
bf_K <- mean(Bluefin_tuna@K) # von Bertalanffy growth parameter K
bf_t0 <- mean(Bluefin_tuna@t0) # von B. theoretical age at length zero
bf_R0 <- Bluefin_tuna@R0 # unfished recruitment
bf_surv <- bf_R0*exp(-bf_M*(ages-1)) # survival
bf_wt_age <- bf_Linf*(1-exp(-bf_K*(ages-bf_t0))) # weight at age
bf_unfished <- sum(bf_R0*bf_surv*bf_wt_age) # approxiate estimate of unfished biomass
M_err <- rlnorm(100,0,0.05)
M_2 <- seq(0.2,0.4,length.out=100) * M_err # made up herring (stock 2) M values
B_1 <- seq(0,bf_unfished,length.out=100) # made up bluefin tuna abundance levels
dat <- data.frame(M_2,B_1)
bfB_hM <- lm(M_2~B_1,dat=dat) # a linear model predicting M for stock 2 from biomass B for stock 1
summary(bfB_hM)
Rel <- list(bfB_hM)
plot(dat$B_1, dat$M_2,pch=19, xlab = "Bluefin biomass", ylab = "Herring natural mortality")
lines(dat$B_1, predict(bfB_hM,newdat=dat), col = 'red', lwd = 2)
points(dat$B_1, simulate(bfB_hM,nsim=1,newdat=dat)[, 1], col = 'green')
legend('topleft',legend = c("made up data", "best fit", "simulated"), text.col = c('black','red','green'))
That is a lot to take in.
We derived a rough level of unfished bluefin biomass, made up some data (normally we would hope to have collected these) with log-normal error in herring M. We then fitted an R model (red line) and also demonstrated that we could use the simulate() R function to simulate (green points) some new data based on that fit. Lastly, we placed our fitted model in a mysterious slot Rel. More on that in a minute.
Here are the important things to know about inter-stock relationships listed in the 'Rel' slot:
Any R model can be used that:
- is compatible with the function
simulate()
- has specific coding for independent (e.g. the bluefin biomass) and dependent variables (e.g. the herring natural mortality rate)
The coding of independent variables goes like this:
- B = total stock biomass
- SSB = total spawning stock biomass
- N = total stock numbers
The coding of dependent variables goes like this:
- M = Natural mortality rate of mature individuals
- K = von Bertalanffy growth parameter K
- Linf = asymtotic size Linf
- t0 = von Bertalanffy theoretical age at length-0 t0
- a = weight length parameter a (W=aL^b)
- b = weight length parameter b (W=aL^b)
- hs = steepness of the Beverton-Holt or Ricker stock recruitment function
The underscore_then_number denotes the stock. So SSB_5 is the spawning stock biomass of the fifth stock in the Stocks slot, Linf_2 is the asymptotic length of stock 2 in the stocks slot.
Currently there can only be one dependent variable but many independent variables so this model is possible:
M_3 ~ B_1 + B_2 - B_4
You can't have transformed independent variables but you can transform the independent variables so this model is possible:
hs_3 ~ log(B_2) * log(B_1)
but not this:
log(hs_3) ~ B_2 + B_1
And another thing, the order you place these in the Rel slot determines the order in which they operate. This may not be consequential yet, but plans are in the works to let an dependent variable in one relationship be the independent in the next.
The idea behind the Rel slot of the MOM object is to open up the option of including ecosystem driven relationships in a terse 'models of intermediate complexity' (MICE) format.
A note of caution before you go hog-wild with the Rel slot: it is probably fairly easy to set up a set of relationships and stock depletions that are impossible to solve in the intialization of the operating models in multiMSE. I haven't had this issue yet, but do give some thought about the proposed relationships before you, for example, make herring M = 5 when bluefin depletion is 0.1 and then specify bluefin depletion as 0.1 and herring at 0.8 ...
Constructing the MOM object
Now we have six lists:
- Stocks
- Fleets
- Obs
- Imps
- CatchFrac
- Rel
We can construct an object of class MOM and see what it looks like:
MOM_BH <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Rel=Rel, nsim=nsim)
MOM_BH <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Rel = Rel, nsim = nsim)
plot(MOM_BH)
Single Stock MOM
If you wish to specify only a single stock but multiple fleets (e.g. just bluefin above) that is pretty easy:
Stocks <- list( Bluefin_tuna )
Fleets <- list( list(Generic_FlatE, Generic_IncE, Generic_DecE) )
Obs <- list( list(Precise_Unbiased, Precise_Unbiased, Precise_Unbiased) )
Imps <- list( list(Perfect_Imp, Perfect_Imp, Overages) )
CatchFrac <- list( matrix( rep(c(0.8, 0.01, 0.19), each=nsim), nrow=nsim) )
MOM_1S <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, nsim=nsim)
Stocks is just a list 1 stock object long. CatchFrac is a list 1 stock long of catch fractions by simulation and fleet. For any hierarchical list (Fleets, Obs, Imps) there is just one position (one stock) in which you have to place a list of Fleet, Obs and Imp objects.
Single Fleet MOM
If you want to model multiple stocks but are happy with aggregate fleet dynamics for each stock, that is fairly straight forward:
# Stock 1 Stock 2
Stocks <- list( Bluefin_tuna, Herring )
# Bluefin Herring
Fleets <- list( list(Generic_FlatE), list(Generic_FlatE) )
Obs <- list( list(Precise_Unbiased), list(Imprecise_Biased) )
Imps <- list( list(Perfect_Imp), list(Overages) )
MOM_1F <- new('MOM', Stocks, Fleets, Obs, Imps, Rel = Rel, nsim=nsim)
Note we no longer have to specify CatchFrac as we do not have to calculate fleet-specific catchabilties using observed recent catch fractions.
A unique fleet for each stock
If you wish to model many stocks but in each case you are able to aggregate fleet dynamics (a single fleet), then you do exactly as above for the single fleet case but put a stock specific-fleet in position 1:
# Stock 1 Stock 2 Stock 3
Stocks <- list( Bluefin_tuna, Herring, Mackerel )
# Bluefin fleet Herring fleet Mackerel fleet
Fleets <- list( list(Generic_FlatE), list(Generic_IncE), list(Generic_DecE) )
Obs <- list( list(Precise_Unbiased), list(Imprecise_Biased), list(Precise_Biased) )
Imps <- list( list(Perfect_Imp), list(Overages), list(Overages) )
CatchFrac <- list( matrix(1, nrow=nsim), matrix(1, nrow=nsim), matrix(1, nrow=nsim) )
MOM_FS <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Rel = Rel, nsim=nsim)
One way or another we have our multi operating model, now we need some MPs to test...
Specifying MPs and various use cases
Let us assume you have specified a multi operating model with more than one stock and more than one fleet. There is a surprisingly large number of options for managing this simulated resource:
- 'complex': an single MP is specified that provides a single recommendation that is applied to all stocks (i.e. a stock complex) and all fleets (e.g. a single TAC is shared across stocks and fleets)
- 'by stock': an MP is specified by stock (e.g. a TAC is divvied up among fleets for a particular stock)
- 'by fleet' an MP is specified for each stock and each fleet based on stock-fleet specific data (e.g. based on fleet-specific data a TAC is provided for each fleet for a particular stock)
- 'multi MP' a single MP is specified that takes all the individual data by stock and fleet and provides an individual recomemndation for each fleet and stock.
Complex: MPs for stock complexes
To configure management of this type you need only provide a vector of MP names:
avail('MP')
MPs_c <- c("DCAC", "DBSRA", "DD")
When you do this, the following assumptions are applied:
- Data from all stocks and fleets are combined. Catches are aggregated. Biological quantities such as depletion are averages weighted by vulnerable biomass. Fleet specific quantities such as mean length in the catch are weighted by vulnerable biomass (among stocks) and by catches (among fleets).
- Among stocks, overall TACs are distributed in proportion to vulnerable biomass.
- Among fleets within stock, the stock TAC is distributed according to
CatchFrac.
- Input controls such as Effort and Size limits apply to all fleets and stocks simultaneously (effort controls are phrased in terms of percentage of today's effort and therefore affect all fleets equally).
A future priority is the development of Linear Programming tools to provide plausible predictions of the distribution of catches among stocks and fleets given boundary constraints on fishing mortality rates and profitability.
By stock: Stock-specific MPs
To specify an MP for each stock, you just list a vector by stock:
MPs_bs<-list( c("DCAC", "DBSRA", "DD"), # Stock 1
c("DCAC", "DCAC", "SPMSY") ) # Stock 2
Note that each position in the vector is a management scenario (we are describing three management systems in the code above). MP1 for Stock 1 and Stock 2 is DCAC for both stocks and hence we are testing a forward projection where DCAC is used on the data of each stock to provide advice for each stock. MP2 is DBSRA for Stock 1 and DCAC for Stock 2, so the second management system evaluates this combination simultaneously, and so on.
If a user wishes to evaluate a full-cross of DCAC and DBSRA MPs for both stocks then they would need four management systems:
MPs<-list( c("DCAC", "DBSRA","DCAC", "DBSRA"), # Stock 1
c("DCAC", "DCAC", "DBSRA","DBSRA") ) # Stock 2
Or, less laboriously, for five stocks and a total of 32 management systems:
myMPs <- c("DCAC", "DBSRA")
nstocks <- 5
MPs <- split(t(expand.grid(list(myMPs)[rep(1,nstocks)])), 1:nstocks)
The 'by stock' mode includes the following assumptions:
- Data are aggregated over fleets. Biological information such as depletion is taken from a single fleet. Fleet specific informaiton such as mean length in the catch is weighted by fleet catches.
- TAC advice is divided between the fleets according to
CatchFrac.
By fleet: fleet and stock specific MPs
Although not very likely, it is possible that users wish to take the individual data of a particular stock and fleet and provide management advice for that fleet and stock.
Believe it or not, the following describes just two management systems that will be tested by closed-loop projection:
# Fleet 1 Fleet 2
# MP1 MP2 MP1 MP2
MPs_bf<-list( list(c("DCAC", "DBSRA"), c("MCD", "AvC")), # Stock 1
list(c("CurE", "DD"), c("AvC","SPMSY")) ) # Stock 2
Management system 2 involved the simultaneous tesing of DBSRA (Stock 1, Fleet 1), AvC (Stock 1, Fleet2), DD (Stock 2, Fleet 1) and SPMSY (Stock 2, Fleet 2).
Multi-fleet and/or Multi-stock management plans
It is conceivable that managers would like to control fishery exploitation holistically, for example accounting for the depletion levels of various stocks on the management of any individual stock. In this case, stock and fleet specific data must be submitted to an MP that then provides advice individually for all stock and fleets.
To do this users need a special class of MP called MMP (multi management procedure) and the user then just specifies a vector of MMPs similarly to the case of 'complex' above:
MPs_mmp <- c("MMP_1", "MMP_2", "MMP_3")
Here is now to code an MMP. The vanilla class of method MP in DLMtool and MSEtool, takes in data in an object of class Data and provides management recommendations in an object of class Rec. The MMP works exactly the same way but it accepts Data objects in a hierarchical list (Fleets nested in Stocks) and provides MP recommendations in a hierarchical list of the same dimensions (Fleets nested in Stocks).
The reason you wish to do this is to allow the data from all fleets and stocks to impact advice of any fleet and stock. However for the sake of demonstrating the format of MMP, I'm going to just make an MMP object that sets TACs to the average historical catches of each fleet and stock.
mydaft_MMP <- function(DataList, reps=1){
nStocks <- length(DataList) # First level is stocks
nFleets <- length(DataList[[1]]) # Second level is fleets (same dimensions among stocks)
RecList <- new('list') # The hierarchical list we are going to put recommendations in
for(ss in 1:nStocks){
RecList[[ss]] <- new('list') # List of recommendations by fleet within stock
for(ff in 1:nFleets){
Rec <- new("Rec") # New blank recommendations object
Rec@TAC <- apply(DataList[[ss]][[ff]]@Cat,1,mean, na.rm = T) # TAC is average historical catch
RecList[[ss]][[ff]]<-Rec # Store recommendations object in RecList
}
}
RecList # Return the RecList
}
class(mydaft_MMP) <- 'MMP' # Assign our new function the correct class
The important thing to note is that a very complex set of rules could have been included in MMP that, for example account for gradients in abundance indices for all species on the change in TAC for any particular species.
MMPs can be considered as a multi-stock, multi-fleet management plans.
Running an multi MSE
Let's take the fantastical bluefin-herring MOM operating model object (MOM_BH) that we created above and run a multi MSE for some stock-specific MPs and plot the results:
MPs_BH <- list( c("DCAC", "DBSRA", "AvC", "MCD"), # Bluefin
c("DCAC", "DBSRA", "AvC", "MCD")) # Herring
MMSE_BH <- multiMSE(MOM_BH, MPs=MPs)
plot(MMSE_BH)
Understanding the MMSE object
The MMSE object produced from multiMSE has a number of slots. The easiest way to understand these is to interrogate the help documentation:
class?MMSE
Other options
Sex - specific multistock operating models
Many fish population exhibit sex-specific dynamics that may have impacts on the performance of candidate management procedures, such as sexual dimorphism and single-sex spawing biomass limitations. Other sex-specific phenomenon that may be important include sequential hermaphroditism and sex-specific vulnerability to exploitation.
In this example we place a female 'stock' in position 1, a male 'stock' in position 2 and use the the SexPars slot of the MOM object to calculate stock-recruitment according to a user-specified fraction of Spawning Stock Biomass originating from these stocks.
MaF<-MaM <- Mackerel # Copy Mackerel dynamics to female and male objects
MaF@Name <- "Mackerel_Female"
MaM@Name <- "Mackerel_Male"
MaF@Linf <- MaM@Linf*1.2 # females grow to 20% higher asymptotic lengths
MaF@M <- MaM@M*1.2 # female natural mortality rate is 20% higher
MaF@K <- MaM@K*1.2 # female fish grow 20% faster
Stocks <- list(MaF,MaM)
Fleets <- list(list(Generic_DecE),list(Generic_DecE)) # A single fleet model
Obs <- list(list(Perfect_Info),list(Perfect_Info))
Imps <- list(list(Perfect_Imp),list(Perfect_Imp))
SexPars <- list()
SexPars$SSBfrom <- matrix(c(1,1,0,0),nrow=length(Stocks)) # Stock (Row) spawn from SSB contributed by Stock (Column)
SexPars$SSBfrom
MOMsex <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, nsim = nsim, SexPars = SexPars)
The list item SSBfrom is a matrix with nstock to rows and nstock from columns. In this case, both rows (i.e. both stocks) get all of their Spawning Biomass from stock 1. There are 1s in column 1 and zeros in column 2.
If you wanted to model two sex-specific stocks where both recruit according to the female stock, the Stock and SSBfrom objects would look like this:
SnF<-SnM <- Snapper # Copy Mackerel dynamics to female and male objects
SnF@Name <- "Snapper_Female"
SnM@Name <- "Snapper_Male"
SnF@Linf <- SnM@Linf*1.1 # females grow to 10% higher asymptotic lengths
Stocks <- list(MaF, MaM, SnF, SnM)
SexPars$SSBfrom <- matrix(c(1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0),nrow=length(Stocks))
SexPars$SSBfrom
Stocks 1 and 2 (female and male Mackerel) recruit according to spawning biomass from stock 1 (female Mackerel SSB) while Stocks 3 and 4 (female and male Snapper) recruit according to spawning biomass from stock 3 (female snapper SSB).
If, for some reason, you want SSB to be calculated by more than one stock, then:
Stocks <- list(MaF, MaM, SnF, SnM)
SexPars$SSBfrom <- matrix(c(1,1,0,0,1,1,0,0,0,0,1,1,0,0,1,1),nrow=length(Stocks))
SexPars$SSBfrom
In this case, the recruitment of male and female mackerel is determined by the aggregate SSB of the two sexes (the same applied for the male and female snapper).
If you want to mix a sex-specific stock with non-sex specific stocks, you just put a 1 in the diagnonal for non-sex specific stocks. E.g:
Stocks <- list(Bluefin_tuna, MaF, MaM, Herring)
SexPars$SSBfrom <- matrix(c(1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1),nrow=length(Stocks))
SexPars$SSBfrom
When specified this way, Bluefin tuna(row 1) and Herring (row 4) recruitment is only determined by their stock-specific spawning biomass (rows 1 and 4, respectively), but Male and Female Mackerel recruitment (rows 2 and 3) depends only on female mackerel spawning stock biomass (column 2).
You can put any fraction you like in these:
SexPars$SSBfrom <- matrix(c(1,0,0,0,0,1,1,0,0,0.1,0.1,0,0,0,0,1),nrow=length(Stocks))
SexPars$SSBfrom
Total SSB for male and female mackerel is now calculated from all of female SSB and 10% of male SSB.
Sequential Hermaphroditism
Assuming you have specified sex-specific stock-recruitment using SexPars, you have the option of adding sequential hermaphroditism with the SexPars$Herm list.
To specify protogyny (Female - Male) where stock 1 is male and stock 2 is female, you include a vector H_1_2 that is the fraction male (Stock 1) at age:
H_1_2<-array(rep(c(0,0,0,0,0,0,0.05,0.1,0.2,0.35,0.65,0.8,0.9,1),each=nsim),c(nsim,14))
MOMsex@SexPars$Herm<-list(H_1_2=H_1_2)
These list objects have to be named H_stockfirst_stocksecond (stocks are sexes in this model) where first and second controls the order of the sequential hermaphroditism. These are always vectors that start with a zero (none of stockfirst) and end in 1 (all of stocksecond).
To specify protandry (Male - Female) you do the same as above but place the Male 'stock' object in stock position 1 and the female in stock position 2.
The H_1_2 type notation works for any position of stocks so:
H_3_11<-array(rep(c(0,0,0,0,0,0,0.05,0.1,0.2,0.35,0.65,0.8,0.9,1),each=nsim),c(nsim,14))
MOMsex@SexPars$Herm<-list(H_1_2 = H_1_2, H_3_11 = H_3_11)
Has two sequential hermaphroditic relationships, one as before exchanges between 'stocks' 1 and 2 and the vector is the fraction at age of 'stock' 1. The second is between 'stocks' 3 and 11 and the vector is the fraction at age of 'stock' 3.
Custom TAC allocations
Recall, that when more than one fleet is specified, fleet-specific catchabilities are calculated from the fraction of current stock catches for each fleet. This is the CatchFrac slot of an MOM object, e.g:
nsim = 12
# Fleet1 Fleet2 Fleet3
bluefin_catchfrac<-c(0.8, 0.01, 0.19)
herring_catchfrac<-c(0.05, 0.6, 0.35)
CatchFrac <- list(
matrix( rep(bluefin_catchfrac, each=nsim), nrow=nsim),
matrix( rep(herring_catchfrac, each=nsim), nrow=nsim)
)
CatchFrac[[1]]
By default, in multi-fleet operating models, future TACs in projected years are allocated to fleet according to CatchFrac. You can however overide this using the MOM slotAllocation. Lets say you want the same OM initialization using CatchFrac but you wish to halve the allocation to fleet 1 and redistribute to equally to fleets 2 and 3 in future years. You can do this with something like:
# Fleet1 Fleet2 Fleet3
bluefin_catchfrac<-c(0.8, 0.01, 0.19)
herring_catchfrac<-c(0.05, 0.6, 0.35)
CatchFrac <- list(
matrix( rep(bluefin_catchfrac, each=nsim), nrow=nsim),
matrix( rep(herring_catchfrac, each=nsim), nrow=nsim)
)
# Fleet1 Fleet2 Fleet3
bluefin_allocation<-c(0.4, 0.21, 0.39)
herring_allocation<-c(0.025, 0.6125, 0.3625)
Allocation <- list(
matrix( rep(bluefin_allocation, each=nsim), nrow=nsim),
matrix( rep(herring_allocation, each=nsim), nrow=nsim)
)
MOM_BH_A <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Allocation=Allocation, nsim=nsim)
Custom effort allocations
When specifying effort control MPs, by default future effort recomemdnations are scaled as a multiple of effort in the current year (e.g. this was 1 in 2018, so an MP halving effort would make this 0.5 for all fleets).
If users wish to control the relative effort of the various fleets in the future they can do this similarly to the Allocation slot above.
# Fleet1 Fleet2 Fleet3
bluefin_fleet_efforts <- c(1, 0.5, 0.0001)
herring_fleet_efforts <- c(1, 1, 1)
Efactor <- list(
matrix( rep(bluefin_fleet_efforts, each=nsim), nrow=nsim),
matrix( rep(herring_fleet_efforts, each=nsim), nrow=nsim)
)
MOM_BH_A <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Efactor=Efactor, nsim=nsim)
In this case all simulations are the same. Effort of the first bluefin fleet and all herring fleets will be as the default: phrased in terms of 100% of current year (e.g. 2018) fishing effort. However, MPs providing effort advice will now provide half the recommended amount to bluefin tuna fleet 2 and nearly zero effort to bluefin fleet 3.
This relatively basic effort distribution code will be improved in later package versions, but for now allows users to investigate alternative management regimes in which fleets can be reduced, increased or removed completely.
Acknowledgements
MSEtool development is funded by the Canadian Department of Fisheries and Oceans and benefits from ongoing collaborations with the Natural Resources Defense Council and number of Canadian government scientists including Robyn Forrest, Sean Anderson and Daniel Duplisea.
The multiMSE function relies heavily on very efficient and robust OMx operating model code developed by Adrian Hordyk in the DLMtool package.
The DLMtool package was developed at the University of British Columbia in collaboration with the Natural Resources Defense Council. DLMtool development has been funded by the Gordon and Betty Moore Foundation, the Packard Foundation, U.S. National Oceanic and Atmospheric Administration, Fisheries and Oceans Canada, the Walton Foundation, Resources Legacy Fund, the Natural Resources Defense Council, the United Nations Food and Agricultural Organization and the California Department of Fish and Wildlife.
tcarruth/MSEtool documentation built on Oct. 19, 2020, 6:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(MSEtool)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>") knitr::opts_chunk$set(dpi=85) options(width = 650)
Introduction
Management Strategy Evaluation can be a complex business even for single-stock and single-fleet operating models. Population dynamics, fleet dynamics, observation processes (e.g. catch over reporting), implementation error (e.g. catch-limit overages) and appropriate management procedures (MPs) all need to be specified.
In most fisheries settings, single-stock, single-fleet operating models are sufficient. Unless management advice is explicitly given on a fleet-by-fleet basis (e.g. two TACs each specific to recreational and commercial fisheries) then the fleet dynamics can be pooled into an aggregate fishery. Unless (A) the dynamics of one stock interact with another stock (e.g. a predator and prey) or (B) stocks with varying dynamics are managed as a complex (e.g. one TAC for all sculpins as in Alaska) then MPs can be tested for one individual stock at a time.
There are situations when managers wish to evaluate the robustness of a proposed management plan among fleets and/or stocks, acknowledging that there are interactions between these. To this end, MSEtool includes the function multiMSE(), a version of DLMtool::runMSE() that conducts MSE for multiple fleets/stocks and allows users to prescribe interactions between stocks using R models (including sex-specific population dyanmics and hermaphroditism).
To accomplish this, you will need to become familiar with a new class of operating model object MOM (multi operating model), a new class of MSE object MMSE (multi MSE), and optionally, a new class of management procedure MMP (multi MP).
Beyond learning about these new objects and functions, multiMSE() is extremely simple to use and hence misuse. I can't stress enough how important it is to have a clear justification for why a multi-stock and/or multi-fleet MSE is necessary for your particular management setting (it may not be). To code multiMSE() it was necessary to identify all of the various use-cases and it became clear that there are a very large number each with their own specific set of requirements and assumptions.
In this guide, I will attempt to clarify these various use cases and be explicit about the assumptions entailed.
Before going crazy with many fleet, many stock models, it might be good to temper your expectations. You can do an awful lot with multiMSE() but numerically solving the multi-stock, multi-fleet equations while maintaining MICE relationships is computationally costly. This is only reflected in the initialization phase of the MSE - which can take some time. A rough rule of thumb is that the minutes to initialize = (nsimulations * nstocks * nfleets) divided by (10 * ncores). So for a 100-simulation MSE with 4 stocks, 5 fleets and parallel computing set up over 4 cores you are talking about 2000 / 40 = 50 minute initialization. Faster solutions are on the way for non-MICE models that don't require R optimizaiton due to R MICE functions.
Caution: you are now wandering into a relatively unestablished and unchartered territory of MSE - good luck!
Remember you can always fire off an email to me if you have questions about multiMSE: t.carruthers@oceans.ubc.ca
For general help, Adrian Hordyk has produced a comprehensive guide for DLMtool including cheat sheets, and Quang Huynh has produced additional user guide vignettes for MSEtool that can be obtained using the standard R help functions:
DLMtool::userguide() DLMtool::cheatsheets() MSEtool::userguide()
The multi operating model object (class MOM)
The MOM object is not terribly different from a standard DLMtool/MSEtool OM object.
library(DLMtool) # load DLMtool library DLMextra() # download lots of cool DLMtool/MSEtool objects from GitHub library(DLMextra) # load all those objects library(MSEtool) # load MSEtool (data rich MPs, MOM object class etc) setup() # set up a cluster for parallel processing
slotNames('MOM') # what slots are there in an MOM object
class?MOM # check out the help file on MOM objects
The big difference between objects of class OM in DLMtool and those of class MOM is that instead of including a slot for each parameter of the stock (Stock), fleet (Fleet), observation (Obs) and implementation (Imp) models for OM objects, MOM objects include these in four slots that are lists of entire Stock, Fleet, Obs and Imp objects - after all it's multi-fleet and/or multi stock, right?!
Stocks, Fleets, Obs and Imps
There are four new slots in an MOM object (compared with an OM object): MOM@Stocks, MOM@Fleets, MOM@Obs, and MOM@Imps.
MOM@Stocksis a list of stock objects (nstocks long)MOM@Fleetsis a hierarchical list of fleet objects (nstocks with fleets nested in stocks)MOM@Obsis a hierarchical list of observation error objects (nstocks with fleets nested in stocks)MOM@Impsis a hierarchical list of implementation error objects (nstocks with fleets nested in stocks)
Ideally you'd get your various dynamics from conditioned operating models either using the DLMtool:: StochasticSRA() approach or using the various MSEtool helper functions that read in stock assessment outputs from, for example, Stock Synthesis (SS2OM) or iSCAM (iSCAM2OM) assessments. However, for simplicity I'll demonstrate how to do multiMSE with some pre-existing Stock, Fleet, Obs, and Imp objects that were loaded with library(DLMextra):
avail('Stock') # what Stock objects are loaded? avail('Fleet') # what Fleet objects are loaded? avail('Obs') # what Obs objects are loaded? avail('Imp') # what Imp objects are loaded?
Let's imagine that we lose our minds and decide to create a model that simulates interactions between bluefin tuna and herring (we think bluefin eat herring). We also want to simulate three fleets, a bluefin fleet, a herring fleet and a mixed fleet that somehow, against all the odds catches both bluefin and herring! About as hypothetical as it gets I know, but that is intentional - this is just a demonstration!
First, we start off by creating the list of stock objects, this is pretty easy:
# Stock 1 Stock 2 Stocks <- list(Bluefin_tuna, Herring)
A key constraint here is that the Stock objects have the same number of areas (e.g. Bluefin_tuna@nareas).
We now want to create our fleets for these stocks. Our bluefin-only fleet has stable effort (premade fleet type Generic_FlatE) our herring-only fleet has increasing effort (premade fleet type Generic_IncE) and our mythical mixed fishery has decreasing effort (premade fleet type Generic_DecE).
This is the tricky bit.
Nominally, the same fleet occupies the same position for each stock. For example, 'longline' is fleet 1 for both bluefin and herring. However, the fleet interactions with the stocks can (and often will) be different - e.g. the trend in exploitation rate of the longline fleet may be stable for bluefin but decreasing for herring (in this admittedly bizarre case). That is why multiMSE() requires the same number of fleets for each stock but you have to be able to specify unique fleet dynamics for each stock x fleet. For example longline_stable_exploitation could be Fleet 1 for bluefin and longline_decreasing_exploitation could be Fleet 1 for Herring - it is the same fleet but interacts differently with the two stocks.
A rare exception is if you have one fleet catching each stock. Then you can just have one fleet position for each stock and put the stock-specific fleet in position Fleet 1 for all stocks. Confused? I don't blame you! But this should get clearer as you work through this guide.
To specify exploitation dynamics across stocks and fleets you create a hierarchical list structure (Fleets nested in Stocks):
# Fleet 1 Fleet 2 Fleet 3 Fleets_both_Stocks <- list(Generic_FlatE, Generic_IncE, Generic_DecE) # for the sake of this exercise we assume identical fleet dynamics for both stocks Fleets <- list( Fleets_both_Stocks, # Bluefin (position 1 in Stocks) Fleets_both_Stocks) # Herring (position 2 in Stocks)
Since each of these Fleets is going to generate fishery data for each stock, they each need an observation model by stock and fleet. Just like the Fleets list, this is hierarchical - Fleets nested in Stocks. We are going to assume that we get good data (Precise and Unbiased) for bluefin and bad data (Imprecise and Biased) for herring.
# Fleet 1 Fleet 2 Fleet 3 Bluefin_obs <- list(Precise_Unbiased, Precise_Unbiased, Precise_Unbiased) Herring_obs <- list(Imprecise_Biased, Imprecise_Biased, Imprecise_Biased) Obs <- list( Bluefin_obs, # Bluefin (position 1 in Stocks) Herring_obs) # Herring (position 2 in Stocks)
Lastly, we need to specify implementation error models for these, we are going to assume this is the same for both stocks and is fleet specific: perfect for the first 2 fleets and includes overages for the mixed fleet (fleet 3).
# Fleet 1 Fleet 2 Fleet 3 Imp_both_Stocks <- list(Perfect_Imp, Perfect_Imp, Overages) Imps <- list( Imp_both_Stocks, # Bluefin (position 1 in Stocks) Imp_both_Stocks) # Herring (position 2 in Stocks)
Each of the stocks already has a prespecified range of depletion:
Bluefin_tuna@D Herring@D
But how do we attribute the magnitude of historical exploitation among the various fleets?
Catch fractions by fleet for each stock
To allow multiMSE() to run we need an extra Thing 1, and should have an additional Thing 2 after that.
Thing 1 is the current (most recent year) catch fractions among fleets for each stock. This information is required to calculate the relative catchabilities among the fleets.
Catch fractions are a list (nstock long) of matrices. Suppose we want 80% of current bluefin catches to come from fleet 1, 1% from fleet 2 and 19% from fleet 3, and also 60% of herring catches to come from Fleet 2, 5% from fleet 1 and the remaining 35% of herring catches from Fleet 3, this is how you would specify that:
nsim = 4 # Fleet1 Fleet2 Fleet3 bluefin_catchfrac<-c(0.8, 0.01, 0.19) herring_catchfrac<-c(0.05, 0.6, 0.35) CatchFrac <- list( matrix( rep(bluefin_catchfrac, each=nsim), nrow=nsim), matrix( rep(herring_catchfrac, each=nsim), nrow=nsim) ) CatchFrac[[1]]
The observant types among you will have noticed that we now need to specify the number of MSE simulations we are going to be doing (nsim). For now this is just 12 for demonstration purposes, but ultimately this would be >148 for a real MSE.
The reason we have done this is to allow for uncertainty in the catch fractions - the user could make catch fractions vary among fleets by simulation. For now, these are the same among simulations, which is the simplest thing, albeit a bit boring.
Specifying inter-stock relationships
You'll remember that I said there was a Thing 2 that we should add to our multi operating model. Ask yourself why we are modelling bluefin and herring together?
Because we want to look at an interaction between these stocks. If we didn't, we would just model these stocks individually, and there would be only one stock in the Stocks slot and only one list of fleets in the Fleets slot (the remaining reason for modelling many stock without interactions is if they are managed by a single MP as a stock complex - we'll come to this in a bit).
If we didn't expect to set fleet-specific advice AND we weren't modelling bluefin and herring interactions, we should just be using the regular runMSE() from DLMtool and MSEtool and calculate the aggregate fleet dynamics (aggregate size vulnerability and trend in apical fishing mortality rate, for example).
But we're demonstrating multiMSE() here, and are therefore going to specify at least one relationship between the stocks (noting that you may wish to evaluate MPs that provide advice for a stock complex of multiple stock objects with no interactions - more below). The slot MOM@Rel is where you put these relationships as a list of R models with quite specific formatting.
Since we think bluefin are predators of herring, we are going to propose a relationship between bluefin abundance and the natural mortality rate of herring. Normally we would have derived this relationship empirically from data (in this case perhaps herring tagging survival estimates regressed on bluefin abundance estimates from a stock assessment).
However for the purposes of this demonstration we're just going to cook something up out of thin air: nominally we are going to assume that herring have a natural mortality rate of 0.2 when there are no bluefin, and 0.4 when bluefin tuna are at unfished levels. To invent this relationship we are going to need to calculate bluefin unfished biomass from the operating model, make some fake data and fit an R model.
Here goes:
ages <- 1:60 bf_M <- mean(Bluefin_tuna@M) # natural mortality rate bf_Linf <- mean(Bluefin_tuna@Linf) # asymptotic length bf_K <- mean(Bluefin_tuna@K) # von Bertalanffy growth parameter K bf_t0 <- mean(Bluefin_tuna@t0) # von B. theoretical age at length zero bf_R0 <- Bluefin_tuna@R0 # unfished recruitment bf_surv <- bf_R0*exp(-bf_M*(ages-1)) # survival bf_wt_age <- bf_Linf*(1-exp(-bf_K*(ages-bf_t0))) # weight at age bf_unfished <- sum(bf_R0*bf_surv*bf_wt_age) # approxiate estimate of unfished biomass M_err <- rlnorm(100,0,0.05) M_2 <- seq(0.2,0.4,length.out=100) * M_err # made up herring (stock 2) M values B_1 <- seq(0,bf_unfished,length.out=100) # made up bluefin tuna abundance levels dat <- data.frame(M_2,B_1) bfB_hM <- lm(M_2~B_1,dat=dat) # a linear model predicting M for stock 2 from biomass B for stock 1 summary(bfB_hM) Rel <- list(bfB_hM)
plot(dat$B_1, dat$M_2,pch=19, xlab = "Bluefin biomass", ylab = "Herring natural mortality") lines(dat$B_1, predict(bfB_hM,newdat=dat), col = 'red', lwd = 2) points(dat$B_1, simulate(bfB_hM,nsim=1,newdat=dat)[, 1], col = 'green') legend('topleft',legend = c("made up data", "best fit", "simulated"), text.col = c('black','red','green'))
That is a lot to take in.
We derived a rough level of unfished bluefin biomass, made up some data (normally we would hope to have collected these) with log-normal error in herring M. We then fitted an R model (red line) and also demonstrated that we could use the simulate() R function to simulate (green points) some new data based on that fit. Lastly, we placed our fitted model in a mysterious slot Rel. More on that in a minute.
Here are the important things to know about inter-stock relationships listed in the 'Rel' slot:
Any R model can be used that:
- is compatible with the function
simulate() - has specific coding for independent (e.g. the bluefin biomass) and dependent variables (e.g. the herring natural mortality rate)
The coding of independent variables goes like this:
- B = total stock biomass
- SSB = total spawning stock biomass
- N = total stock numbers
The coding of dependent variables goes like this:
- M = Natural mortality rate of mature individuals
- K = von Bertalanffy growth parameter K
- Linf = asymtotic size Linf
- t0 = von Bertalanffy theoretical age at length-0 t0
- a = weight length parameter a (W=aL^b)
- b = weight length parameter b (W=aL^b)
- hs = steepness of the Beverton-Holt or Ricker stock recruitment function
The underscore_then_number denotes the stock. So SSB_5 is the spawning stock biomass of the fifth stock in the Stocks slot, Linf_2 is the asymptotic length of stock 2 in the stocks slot.
Currently there can only be one dependent variable but many independent variables so this model is possible:
M_3 ~ B_1 + B_2 - B_4
You can't have transformed independent variables but you can transform the independent variables so this model is possible:
hs_3 ~ log(B_2) * log(B_1)
but not this:
log(hs_3) ~ B_2 + B_1
And another thing, the order you place these in the Rel slot determines the order in which they operate. This may not be consequential yet, but plans are in the works to let an dependent variable in one relationship be the independent in the next.
The idea behind the Rel slot of the MOM object is to open up the option of including ecosystem driven relationships in a terse 'models of intermediate complexity' (MICE) format.
A note of caution before you go hog-wild with the Rel slot: it is probably fairly easy to set up a set of relationships and stock depletions that are impossible to solve in the intialization of the operating models in multiMSE. I haven't had this issue yet, but do give some thought about the proposed relationships before you, for example, make herring M = 5 when bluefin depletion is 0.1 and then specify bluefin depletion as 0.1 and herring at 0.8 ...
Constructing the MOM object
Now we have six lists:
- Stocks
- Fleets
- Obs
- Imps
- CatchFrac
- Rel
We can construct an object of class MOM and see what it looks like:
MOM_BH <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Rel=Rel, nsim=nsim)
MOM_BH <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Rel = Rel, nsim = nsim) plot(MOM_BH)
Single Stock MOM
If you wish to specify only a single stock but multiple fleets (e.g. just bluefin above) that is pretty easy:
Stocks <- list( Bluefin_tuna ) Fleets <- list( list(Generic_FlatE, Generic_IncE, Generic_DecE) ) Obs <- list( list(Precise_Unbiased, Precise_Unbiased, Precise_Unbiased) ) Imps <- list( list(Perfect_Imp, Perfect_Imp, Overages) ) CatchFrac <- list( matrix( rep(c(0.8, 0.01, 0.19), each=nsim), nrow=nsim) ) MOM_1S <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, nsim=nsim)
Stocks is just a list 1 stock object long. CatchFrac is a list 1 stock long of catch fractions by simulation and fleet. For any hierarchical list (Fleets, Obs, Imps) there is just one position (one stock) in which you have to place a list of Fleet, Obs and Imp objects.
Single Fleet MOM
If you want to model multiple stocks but are happy with aggregate fleet dynamics for each stock, that is fairly straight forward:
# Stock 1 Stock 2 Stocks <- list( Bluefin_tuna, Herring ) # Bluefin Herring Fleets <- list( list(Generic_FlatE), list(Generic_FlatE) ) Obs <- list( list(Precise_Unbiased), list(Imprecise_Biased) ) Imps <- list( list(Perfect_Imp), list(Overages) ) MOM_1F <- new('MOM', Stocks, Fleets, Obs, Imps, Rel = Rel, nsim=nsim)
Note we no longer have to specify CatchFrac as we do not have to calculate fleet-specific catchabilties using observed recent catch fractions.
A unique fleet for each stock
If you wish to model many stocks but in each case you are able to aggregate fleet dynamics (a single fleet), then you do exactly as above for the single fleet case but put a stock specific-fleet in position 1:
# Stock 1 Stock 2 Stock 3 Stocks <- list( Bluefin_tuna, Herring, Mackerel ) # Bluefin fleet Herring fleet Mackerel fleet Fleets <- list( list(Generic_FlatE), list(Generic_IncE), list(Generic_DecE) ) Obs <- list( list(Precise_Unbiased), list(Imprecise_Biased), list(Precise_Biased) ) Imps <- list( list(Perfect_Imp), list(Overages), list(Overages) ) CatchFrac <- list( matrix(1, nrow=nsim), matrix(1, nrow=nsim), matrix(1, nrow=nsim) ) MOM_FS <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Rel = Rel, nsim=nsim)
One way or another we have our multi operating model, now we need some MPs to test...
Specifying MPs and various use cases
Let us assume you have specified a multi operating model with more than one stock and more than one fleet. There is a surprisingly large number of options for managing this simulated resource:
- 'complex': an single MP is specified that provides a single recommendation that is applied to all stocks (i.e. a stock complex) and all fleets (e.g. a single TAC is shared across stocks and fleets)
- 'by stock': an MP is specified by stock (e.g. a TAC is divvied up among fleets for a particular stock)
- 'by fleet' an MP is specified for each stock and each fleet based on stock-fleet specific data (e.g. based on fleet-specific data a TAC is provided for each fleet for a particular stock)
- 'multi MP' a single MP is specified that takes all the individual data by stock and fleet and provides an individual recomemndation for each fleet and stock.
Complex: MPs for stock complexes
To configure management of this type you need only provide a vector of MP names:
avail('MP') MPs_c <- c("DCAC", "DBSRA", "DD")
When you do this, the following assumptions are applied:
- Data from all stocks and fleets are combined. Catches are aggregated. Biological quantities such as depletion are averages weighted by vulnerable biomass. Fleet specific quantities such as mean length in the catch are weighted by vulnerable biomass (among stocks) and by catches (among fleets).
- Among stocks, overall TACs are distributed in proportion to vulnerable biomass.
- Among fleets within stock, the stock TAC is distributed according to
CatchFrac. - Input controls such as Effort and Size limits apply to all fleets and stocks simultaneously (effort controls are phrased in terms of percentage of today's effort and therefore affect all fleets equally).
A future priority is the development of Linear Programming tools to provide plausible predictions of the distribution of catches among stocks and fleets given boundary constraints on fishing mortality rates and profitability.
By stock: Stock-specific MPs
To specify an MP for each stock, you just list a vector by stock:
MPs_bs<-list( c("DCAC", "DBSRA", "DD"), # Stock 1 c("DCAC", "DCAC", "SPMSY") ) # Stock 2
Note that each position in the vector is a management scenario (we are describing three management systems in the code above). MP1 for Stock 1 and Stock 2 is DCAC for both stocks and hence we are testing a forward projection where DCAC is used on the data of each stock to provide advice for each stock. MP2 is DBSRA for Stock 1 and DCAC for Stock 2, so the second management system evaluates this combination simultaneously, and so on.
If a user wishes to evaluate a full-cross of DCAC and DBSRA MPs for both stocks then they would need four management systems:
MPs<-list( c("DCAC", "DBSRA","DCAC", "DBSRA"), # Stock 1 c("DCAC", "DCAC", "DBSRA","DBSRA") ) # Stock 2
Or, less laboriously, for five stocks and a total of 32 management systems:
myMPs <- c("DCAC", "DBSRA") nstocks <- 5 MPs <- split(t(expand.grid(list(myMPs)[rep(1,nstocks)])), 1:nstocks)
The 'by stock' mode includes the following assumptions:
- Data are aggregated over fleets. Biological information such as depletion is taken from a single fleet. Fleet specific informaiton such as mean length in the catch is weighted by fleet catches.
- TAC advice is divided between the fleets according to
CatchFrac.
By fleet: fleet and stock specific MPs
Although not very likely, it is possible that users wish to take the individual data of a particular stock and fleet and provide management advice for that fleet and stock.
Believe it or not, the following describes just two management systems that will be tested by closed-loop projection:
# Fleet 1 Fleet 2 # MP1 MP2 MP1 MP2 MPs_bf<-list( list(c("DCAC", "DBSRA"), c("MCD", "AvC")), # Stock 1 list(c("CurE", "DD"), c("AvC","SPMSY")) ) # Stock 2
Management system 2 involved the simultaneous tesing of DBSRA (Stock 1, Fleet 1), AvC (Stock 1, Fleet2), DD (Stock 2, Fleet 1) and SPMSY (Stock 2, Fleet 2).
Multi-fleet and/or Multi-stock management plans
It is conceivable that managers would like to control fishery exploitation holistically, for example accounting for the depletion levels of various stocks on the management of any individual stock. In this case, stock and fleet specific data must be submitted to an MP that then provides advice individually for all stock and fleets.
To do this users need a special class of MP called MMP (multi management procedure) and the user then just specifies a vector of MMPs similarly to the case of 'complex' above:
MPs_mmp <- c("MMP_1", "MMP_2", "MMP_3")
Here is now to code an MMP. The vanilla class of method MP in DLMtool and MSEtool, takes in data in an object of class Data and provides management recommendations in an object of class Rec. The MMP works exactly the same way but it accepts Data objects in a hierarchical list (Fleets nested in Stocks) and provides MP recommendations in a hierarchical list of the same dimensions (Fleets nested in Stocks).
The reason you wish to do this is to allow the data from all fleets and stocks to impact advice of any fleet and stock. However for the sake of demonstrating the format of MMP, I'm going to just make an MMP object that sets TACs to the average historical catches of each fleet and stock.
mydaft_MMP <- function(DataList, reps=1){ nStocks <- length(DataList) # First level is stocks nFleets <- length(DataList[[1]]) # Second level is fleets (same dimensions among stocks) RecList <- new('list') # The hierarchical list we are going to put recommendations in for(ss in 1:nStocks){ RecList[[ss]] <- new('list') # List of recommendations by fleet within stock for(ff in 1:nFleets){ Rec <- new("Rec") # New blank recommendations object Rec@TAC <- apply(DataList[[ss]][[ff]]@Cat,1,mean, na.rm = T) # TAC is average historical catch RecList[[ss]][[ff]]<-Rec # Store recommendations object in RecList } } RecList # Return the RecList } class(mydaft_MMP) <- 'MMP' # Assign our new function the correct class
The important thing to note is that a very complex set of rules could have been included in MMP that, for example account for gradients in abundance indices for all species on the change in TAC for any particular species.
MMPs can be considered as a multi-stock, multi-fleet management plans.
Running an multi MSE
Let's take the fantastical bluefin-herring MOM operating model object (MOM_BH) that we created above and run a multi MSE for some stock-specific MPs and plot the results:
MPs_BH <- list( c("DCAC", "DBSRA", "AvC", "MCD"), # Bluefin c("DCAC", "DBSRA", "AvC", "MCD")) # Herring MMSE_BH <- multiMSE(MOM_BH, MPs=MPs) plot(MMSE_BH)
Understanding the MMSE object
The MMSE object produced from multiMSE has a number of slots. The easiest way to understand these is to interrogate the help documentation:
class?MMSE
Other options
Sex - specific multistock operating models
Many fish population exhibit sex-specific dynamics that may have impacts on the performance of candidate management procedures, such as sexual dimorphism and single-sex spawing biomass limitations. Other sex-specific phenomenon that may be important include sequential hermaphroditism and sex-specific vulnerability to exploitation.
In this example we place a female 'stock' in position 1, a male 'stock' in position 2 and use the the SexPars slot of the MOM object to calculate stock-recruitment according to a user-specified fraction of Spawning Stock Biomass originating from these stocks.
MaF<-MaM <- Mackerel # Copy Mackerel dynamics to female and male objects MaF@Name <- "Mackerel_Female" MaM@Name <- "Mackerel_Male" MaF@Linf <- MaM@Linf*1.2 # females grow to 20% higher asymptotic lengths MaF@M <- MaM@M*1.2 # female natural mortality rate is 20% higher MaF@K <- MaM@K*1.2 # female fish grow 20% faster Stocks <- list(MaF,MaM) Fleets <- list(list(Generic_DecE),list(Generic_DecE)) # A single fleet model Obs <- list(list(Perfect_Info),list(Perfect_Info)) Imps <- list(list(Perfect_Imp),list(Perfect_Imp)) SexPars <- list() SexPars$SSBfrom <- matrix(c(1,1,0,0),nrow=length(Stocks)) # Stock (Row) spawn from SSB contributed by Stock (Column) SexPars$SSBfrom MOMsex <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, nsim = nsim, SexPars = SexPars)
The list item SSBfrom is a matrix with nstock to rows and nstock from columns. In this case, both rows (i.e. both stocks) get all of their Spawning Biomass from stock 1. There are 1s in column 1 and zeros in column 2.
If you wanted to model two sex-specific stocks where both recruit according to the female stock, the Stock and SSBfrom objects would look like this:
SnF<-SnM <- Snapper # Copy Mackerel dynamics to female and male objects SnF@Name <- "Snapper_Female" SnM@Name <- "Snapper_Male" SnF@Linf <- SnM@Linf*1.1 # females grow to 10% higher asymptotic lengths Stocks <- list(MaF, MaM, SnF, SnM) SexPars$SSBfrom <- matrix(c(1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0),nrow=length(Stocks)) SexPars$SSBfrom
Stocks 1 and 2 (female and male Mackerel) recruit according to spawning biomass from stock 1 (female Mackerel SSB) while Stocks 3 and 4 (female and male Snapper) recruit according to spawning biomass from stock 3 (female snapper SSB).
If, for some reason, you want SSB to be calculated by more than one stock, then:
Stocks <- list(MaF, MaM, SnF, SnM) SexPars$SSBfrom <- matrix(c(1,1,0,0,1,1,0,0,0,0,1,1,0,0,1,1),nrow=length(Stocks)) SexPars$SSBfrom
In this case, the recruitment of male and female mackerel is determined by the aggregate SSB of the two sexes (the same applied for the male and female snapper).
If you want to mix a sex-specific stock with non-sex specific stocks, you just put a 1 in the diagnonal for non-sex specific stocks. E.g:
Stocks <- list(Bluefin_tuna, MaF, MaM, Herring) SexPars$SSBfrom <- matrix(c(1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1),nrow=length(Stocks)) SexPars$SSBfrom
When specified this way, Bluefin tuna(row 1) and Herring (row 4) recruitment is only determined by their stock-specific spawning biomass (rows 1 and 4, respectively), but Male and Female Mackerel recruitment (rows 2 and 3) depends only on female mackerel spawning stock biomass (column 2).
You can put any fraction you like in these:
SexPars$SSBfrom <- matrix(c(1,0,0,0,0,1,1,0,0,0.1,0.1,0,0,0,0,1),nrow=length(Stocks)) SexPars$SSBfrom
Total SSB for male and female mackerel is now calculated from all of female SSB and 10% of male SSB.
Sequential Hermaphroditism
Assuming you have specified sex-specific stock-recruitment using SexPars, you have the option of adding sequential hermaphroditism with the SexPars$Herm list.
To specify protogyny (Female - Male) where stock 1 is male and stock 2 is female, you include a vector H_1_2 that is the fraction male (Stock 1) at age:
H_1_2<-array(rep(c(0,0,0,0,0,0,0.05,0.1,0.2,0.35,0.65,0.8,0.9,1),each=nsim),c(nsim,14)) MOMsex@SexPars$Herm<-list(H_1_2=H_1_2)
These list objects have to be named H_stockfirst_stocksecond (stocks are sexes in this model) where first and second controls the order of the sequential hermaphroditism. These are always vectors that start with a zero (none of stockfirst) and end in 1 (all of stocksecond).
To specify protandry (Male - Female) you do the same as above but place the Male 'stock' object in stock position 1 and the female in stock position 2.
The H_1_2 type notation works for any position of stocks so:
H_3_11<-array(rep(c(0,0,0,0,0,0,0.05,0.1,0.2,0.35,0.65,0.8,0.9,1),each=nsim),c(nsim,14)) MOMsex@SexPars$Herm<-list(H_1_2 = H_1_2, H_3_11 = H_3_11)
Has two sequential hermaphroditic relationships, one as before exchanges between 'stocks' 1 and 2 and the vector is the fraction at age of 'stock' 1. The second is between 'stocks' 3 and 11 and the vector is the fraction at age of 'stock' 3.
Custom TAC allocations
Recall, that when more than one fleet is specified, fleet-specific catchabilities are calculated from the fraction of current stock catches for each fleet. This is the CatchFrac slot of an MOM object, e.g:
nsim = 12 # Fleet1 Fleet2 Fleet3 bluefin_catchfrac<-c(0.8, 0.01, 0.19) herring_catchfrac<-c(0.05, 0.6, 0.35) CatchFrac <- list( matrix( rep(bluefin_catchfrac, each=nsim), nrow=nsim), matrix( rep(herring_catchfrac, each=nsim), nrow=nsim) ) CatchFrac[[1]]
By default, in multi-fleet operating models, future TACs in projected years are allocated to fleet according to CatchFrac. You can however overide this using the MOM slotAllocation. Lets say you want the same OM initialization using CatchFrac but you wish to halve the allocation to fleet 1 and redistribute to equally to fleets 2 and 3 in future years. You can do this with something like:
# Fleet1 Fleet2 Fleet3 bluefin_catchfrac<-c(0.8, 0.01, 0.19) herring_catchfrac<-c(0.05, 0.6, 0.35) CatchFrac <- list( matrix( rep(bluefin_catchfrac, each=nsim), nrow=nsim), matrix( rep(herring_catchfrac, each=nsim), nrow=nsim) ) # Fleet1 Fleet2 Fleet3 bluefin_allocation<-c(0.4, 0.21, 0.39) herring_allocation<-c(0.025, 0.6125, 0.3625) Allocation <- list( matrix( rep(bluefin_allocation, each=nsim), nrow=nsim), matrix( rep(herring_allocation, each=nsim), nrow=nsim) ) MOM_BH_A <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Allocation=Allocation, nsim=nsim)
Custom effort allocations
When specifying effort control MPs, by default future effort recomemdnations are scaled as a multiple of effort in the current year (e.g. this was 1 in 2018, so an MP halving effort would make this 0.5 for all fleets).
If users wish to control the relative effort of the various fleets in the future they can do this similarly to the Allocation slot above.
# Fleet1 Fleet2 Fleet3 bluefin_fleet_efforts <- c(1, 0.5, 0.0001) herring_fleet_efforts <- c(1, 1, 1) Efactor <- list( matrix( rep(bluefin_fleet_efforts, each=nsim), nrow=nsim), matrix( rep(herring_fleet_efforts, each=nsim), nrow=nsim) ) MOM_BH_A <- new('MOM', Stocks, Fleets, Obs, Imps, CatchFrac, Efactor=Efactor, nsim=nsim)
In this case all simulations are the same. Effort of the first bluefin fleet and all herring fleets will be as the default: phrased in terms of 100% of current year (e.g. 2018) fishing effort. However, MPs providing effort advice will now provide half the recommended amount to bluefin tuna fleet 2 and nearly zero effort to bluefin fleet 3.
This relatively basic effort distribution code will be improved in later package versions, but for now allows users to investigate alternative management regimes in which fleets can be reduced, increased or removed completely.
Acknowledgements
MSEtool development is funded by the Canadian Department of Fisheries and Oceans and benefits from ongoing collaborations with the Natural Resources Defense Council and number of Canadian government scientists including Robyn Forrest, Sean Anderson and Daniel Duplisea.
The multiMSE function relies heavily on very efficient and robust OMx operating model code developed by Adrian Hordyk in the DLMtool package.
The DLMtool package was developed at the University of British Columbia in collaboration with the Natural Resources Defense Council. DLMtool development has been funded by the Gordon and Betty Moore Foundation, the Packard Foundation, U.S. National Oceanic and Atmospheric Administration, Fisheries and Oceans Canada, the Walton Foundation, Resources Legacy Fund, the Natural Resources Defense Council, the United Nations Food and Agricultural Organization and the California Department of Fish and Wildlife.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.