vignettes/multispecies_model.Rmd
multispecies_model.RmdThe previous sections have used wrapper functions to set up MizerParams objects that are appropriate for community and trait-based simulations. We now turn our attention to multispecies, or species-specific, models. These are potentially more complicated than the community and trait-based models and use the full power of the mizer package.
In multispecies type models multiple species are resolved. However, unlike in the trait-based model which also resolves multiple species, explicit species are represented. There are several advantages to this approach. As well as investigating the community as a whole (as was done for the community and trait-based models), we are able to investigate the dynamics of individual species. This means that species specific management rules can be tested and species specific metrics, such as yield, can be compared to reference levels.
A multispecies model can take more effort to set up. For example, each species will have different life-history parameters; there may be multiple gear types with different selectivities targeting different groups of species; the fishing effort of each gear may change with time instead of just being constant (which has been the case in the simulations we have looked at so far); the interactions between the species needs to be considered.
In later sections we build up a multispecies model for the North Sea. To effectively use mizer for a multispecies model we are going to have to take a closer look at the MizerParams class and the project() method. This will all be done in the context of examples so hopefully everything will be clear.
We also take a closer look at some of the summary plots and analyses that can be performed, for example, calculating a range of size-based indicators.
The MizerParams class is used for storing model parameters. We have already met the MizerParams class when we looked at community and trait-based models. However, direct handling of the class was largely hidden by the use of the wrapper functions. To set up a multispecies model we need to directly create and use the MizerParams class. This is probably the most complicated part of using the mizer package so we will take it slowly. For additional help you can look at the help page for the class by entering class ? MizerParams.
The MizerParams class stores the:
Note that the MizerParams class does not store any parameters that can vary through time, such as fishing effort or population abundance. These are stored in the MizerSim class which we will come to later in the section on running a simulation.
Although the MizerParams class seems complicated, it is relatively straightforward to set up and use. Objects of class MizerParams are created using the constructor method MizerParams() (this constructor method was called MizerParams() in previous version of mizer). This constructor method can take many arguments. However, creation is simplified because many of the arguments have default values.
In the rest of this section we look at the main arguments to the MizerParams() function. To help understand how the constructor is used and how the MizerParams class relates to the equations given in the model description section, there is an example section where we create example parameter objects using data that comes with the mizer package.
Although many of the arguments used when creating a MizerParams object are optional, there is one argument that must be supplied by the user: the species specific parameters. These are stored in a single data.frame object. The data.frame is arranged species by parameter, so each column is a parameter and each row has the parameters for one of the species in the model. Although it is possible to create the data.frame by hand in R, it is probably easier to create the data externally as a .csv file (perhaps using a suitable open source spreadsheet such as LibreOffice) and then read the data into R.
For each species in the model community there are certain parameters that are essential and that do not have default values. The user must provide values for these parameters. There are also some essential parameters that have default values, such as the selectivity function parameters, and some that are calculated internally using default relationships if not explicitly provided. These defaults are used if the parameters are not found in the data.frame. The species parameter information is stored as a data.frame which hold the information given in the MizerParams help page.
The essential columns of the species parameters data.frame that have no default values are: species, the names of the species in the community; w_inf, the asymptotic mass of the species; w_mat, the mass at maturation; beta and sigma, the predator-prey mass ratio parameters \(\beta\) and \(\sigma\) (see the section on predator/prey mass ratio) and the stock-recruitment parameters (by default, the stock-recruitment function is a Beverton-Holt type and so, unless a different SRR function is used, a r_max column must be provided - see the section on the stock-recruitment relationship for more details).
Essential columns that have default values are: k, the activity coefficient (default = 0); alpha, the assimilation efficiency (default = 0.6); erepro, the reproductive efficiency (default value = 1); w_min, the size of recruits (default value is the smallest size of the community size spectrum); sel_func, the name of the fishing selectivity function (default value = “knife_edge”); gear, the name of the fishing gear that catches the species (default value is the name of the species in the species column); catchability, the catchability of the fishing gear on that species (default value = 1). Additionally, columns that contain the selectivity function parameters are needed. As mentioned, the default selectivity function is a knife_edge function. This has an argument knife_edge_size that determines the knife-edge position and has a default value of w_min. If any of these columns are not included in the species parameter data.frame, the default values are used.
As mentioned above, there are some columns that are essential but if they are not provided, values for them are estimated using the values in the other columns. These columns are: h, the maximum food intake; gamma, the volumetric search rate; ks, the coefficient for standard metabolism; z0, the mortality from other sources, \(\mu_{b,i}.\) These parameters can be included as columns in the species parameters data.frame if they are available. If they are not provided then the MizerParams() construction method will try to calculate them.
The h column is calculated as:
\[\begin{equation} h = \frac{3 k_{vb}}{\alpha f_0} w_\infty^{1/3} \end{equation}\]
where \(k_{vb}\) is the von Bertalanffy \(K\) parameter and \(f_0\) is the feeding level of small individuals feeding mainly on the plankton. This mean that if an h column is not included in the species parameter data.frame, a column for k_vb is necessary. If it is not included then the MizerParams() method will fail with an error message. The calculation also requires a value of the feeding level of small individuals. This can be passed as an additional argument, f0, to the MizerParams() constructor and it has a default value of 0.6.
The gamma column is calculated using the equation for volumetric search rate. This calculation requires that the h column is available, either included in the species parameter data.frame, or calculated internally using the k_vb column as described above. This means that if you include a k_vb column in the data.frame, h and gamma will be calculated from it.
The z0 column (mortality from other sources) is calculated using the equation \(\mu_{b.i}(w) = \mu_b w_{\infty.i}^{1-n}.\) The values z0pre (equivalent to \(\mu_0\)) and z0exp (the power that \(w_\infty\) is raised to) can be passed as arguments to the MizerParams() constructor and have default values of 0.6 and n-1 (where n has a default value of 2/3) respectively.
The ks column is calculated as 20% of h, i.e. the standard metabolism coefficient is 20% of the maximum consumption. Standard metabolism is used in calculation of the growth of individuals as \(k_{s.i} w_i^p\) (see the section on allometric rates). The species parameter information is stored as a data.frame which hold the information described in the MizerParams help page.
In this section we take a look at how fishing is implemented and how fishing gears are set up within mizer.
In mizer, fishing mortality is imposed on species by fishing gears. The total fishing mortality is obtained by summing over the mortality from all gears, \[\begin{equation}
% {#eq:muf}
\mu_{f.i}(w) = \sum_g F_{g,i}(w),
\end{equation}\] where the fishing mortality \(F_{g,i}(w)\) imposed by gear \(g\) on species \(i\) at size \(w\) is calculated as: \[\begin{equation}
% {#eq:sel}
F_{g,i}(w) = S_{g,i}(w) Q_{g,i} E_{g}
\end{equation}\] where \(S\) is the selectivity by species, gear and size, \(Q\) is the catchability by species and gear and \(E\) is the fishing effort by gear. The selectivity at size has a range between 0 (not selected at that size) to 1 (fully selected at that size). Catchability is used as an additional scalar to make the link between gear selectivity, fishing effort and fishing mortality. For example, it can be set so that an effort of 1 gives a desired fishing mortality. In this way effort can then be specified relative to a ‘base effort’, e.g. the effort in a particular year. Fishing effort is not stored in the MizerParams object. Instead, effort is set when the simulation is run and can vary through time (see the section on running a simulation).
At the moment a species can only be selected by one fishing gear, although each gear can select more than one species (this is a limitation with the current package that will be developed in future releases).
The selectivity at size of each gear is given by a selectivity function. Some selectivity functions are included in the package. New functions can be defined by the user. Each gear has the same selectivity function for all the species it selects, but the parameter values for each species may be different, e.g. the lengths of species that a gear selects may be different.
The name of the selectivity function is given by the sel_func column in the species parameters data.frame. Each selectivity function has a range of arguments. Values for these arguments must be included as columns in the species parameters data.frame. The names of the columns must exactly match the names of the arguments. For example, the default selectivity function is knife_edge which has sudden change of selectivity from 0 to 1 at a certain size. The arguments for this selectivity function can be seen in the help page for this function. To see them enter:
It can be seen that the knife_edge() function has arguments w and knife_edge_size The first argument, w, is size (the function calculates selectivity at size). All selectivity functions must have w as the first argument. The values for the other arguments must be found in the species parameters data.frame. So for the knife_edge() function there should be a knife_edge_size column (but note that because knife_edge() is the default selectivity function, the knife_edge_size argument actually has a default value = w_mat). This can be seen in the example in the section with examples of making a MizerParams object. If the columns of the selectivity function arguments are not in the species parameter data.frame, an error is thrown when the MizerParams object is created.
Users are able to write their own size based selectivity function. The first argument to the function must be w and the function must return a vector of the selectivity (between 0 and 1) at size.
The name of the fishing gear is given in the gear column of the species parameter data.frame. If the gear column is not specified, the default gear name is simply the name of the species. This implies that each species is fished by a different gear. This approach can be used to explore the impacts of changing fishing mortality on individual species.
In size spectrum modelling recruitment refers to the flux of individuals that enter the size-spectrum at the smallest size group of that species (given by the parameter w_min in the species parameter data.frame). As can be seen in the section on reproduction, calculating the recruitment flux involves calculating the density independent recruitment, \(R_{p.i}\). The \(R_{p.i}\) is then modified by a stock-recruitment relationship (SRR) to impose some form of density-dependence. This then results in the density-dependent recruitment, \(R_i\) (see the section on the stock-recruitment relationship). Without this density dependence, the realised recruitment flux to the smallest size class is determined only by \(R_{p.i}\). The default SRR is a Beverton-Holt type function.
Similar to the fishing selectivity functions, any parameter used in the stock-recruitment function, other than \(R_{p.i}\), must be in the species parameter data.frame and the column must have the same name as the function argument. For example, the default stock-recruitment function has a second argument to it, r_max. Therefore the species parameter data.frame must have an r_max column.
Users are able to write their own stock-recruitment function. The first argument to the function must be rdi, which is the density independent recruitment \(R_{p.i}\), the second argument must be species_params, the dataframe with the species parameters. This function should then be assigned to the srr slot of the MizerParams object.
The interaction matrix describes the interaction of each pair of species in the model. This can be viewed as a proxy for spatial interaction e.g. to model predator-prey interaction that is not size based. The values in the interaction matrix are used to scale the encountered food and predation mortality (see the section on predator-prey preference). The matrix is square with every element being the interaction between a pair of species. The dimensions, nrows and ncolumns, therefore equal the number of species. The values are between 0 (species do not overlap and therefore do not interact with each other) to 1 (species overlap perfectly). If all the values in the interaction matrix are set to 1 then predator-prey interactions are determined entirely by size-preference.
The interaction matrix must be of type array or matrix. One way of creating your own is to enter the data using a spreadsheet (such as LibreOffice) and saving it as a .csv file. The data can be read into R using the command read.csv(). This reads in the data as a data.frame. We then need to convert this to a matrix using the as() function. An example of how to do this is given in the section with examples of making MizerParams objects.
It should be noted that the order of species in the interaction matrix has to be the same as the order in the species parameters data.frame. Although you can specify the dimnames of the interaction matrix, these names are overwritten by the species names from the species parameters data.frame inside the MizerParams constructor.
If an interaction matrix is not specified to the MizerParams() constructor the default interaction matrix is used. This has all values set to 1.
MizerParams() argumentsAs well as the essential species parameters data.frame and the interaction matrix, there are several other arguments to the MizerParams constructor. These have default values. The arguments can be seen in the MizerParams() help page.
Some of these parameters may be used to calculate the species specific parameters if they are not provided. For example, if there is no gamma column in the species parameter data.frame, then it is calculated using kappa and f0. This means that depending on which columns have been provided in the species parameters data.frame, some of the parameters from the the MizerParams() help page may not be used. For example, if the column z0 (the mortality from other sources) has been included, then the arguments z0pre and z0exp are not used.
Determining a value for the kappa argument can be difficult and may need to be estimated through some kind calibration process. The default value kappa is for the North Sea model.
MizerParams objectsAs mentioned in the preceding sections, an object of MizerParams is created by using the MizerParams() constructor method. You can see the help page for the constructor:
This shows that the constructor takes the following arguments:
object The species parameter data.frame (see the section on species parameters). This is compulsory with no default value.
inter The interaction matrix (see the section on predator-prey preference). The default is a matrix of 1s.
... Other model parameters.
In the rest of this section we demonstrate how to pull these elements together to make MizerParams objects.
The first step is to prepare the species specific parameter data.frame. As mentioned above, one way of doing this is to use a spreadsheet and save it as a .csv file. We will use this approach here. An example .csv file has been included in the package. This contains the species parameters for a multispecies North Sea model. This file is placed in the doc folder of the package installation. The location of the file can be found by running:
This file can be opened with most spreadsheets or a text editor for you to inspect. This can be loaded into R using the following code (after you have told R to look in the right directory):
This reads the .csv file into R in the form of a data.frame. You can check this with the class:
## [1] "data.frame"The example data.frame can be inspected by entering the name of the object.
## species w_inf w_mat beta sigma r_max k_vb
## 1 Sprat 33 13 51076 0.8 7.38e+11 0.681
## 2 Sandeel 36 4 398849 1.9 4.10e+11 1.000
## 3 N.pout 100 23 22 1.5 1.05e+13 0.849
## 4 Herring 334 99 280540 3.2 1.11e+12 0.606
## 5 Dab 324 21 191 1.9 1.12e+10 0.536
## 6 Whiting 1192 75 22 1.5 5.48e+11 0.323
## 7 Sole 866 78 381 1.9 3.87e+10 0.284
## 8 Gurnard 668 39 283 1.8 1.65e+12 0.266
## 9 Plaice 2976 105 113 1.6 4.08e+14 0.122
## 10 Haddock 3485 165 558 2.1 1.84e+12 0.271
## 11 Cod 40044 1606 66 1.3 8.26e+09 0.216
## 12 Saithe 16856 1076 40 1.1 1.12e+11 0.175You can see that there are \(12\) species and \(7\) columns of parameters: species, w_inf,w_mat,beta,sigma,r_max and k_vb.
Of these parameters, species, w_inf, w_mat, beta and sigma are essential and have no default values (as described in the section on species parameters). r_max is a SRR parameter. We are going to use the default Beverton-Holt type SRR which has r_max as an argument (see the section on the stock-recruitment relationship), making this column also essential. The final column, k_vb, will be used to calculate values for h and then gamma. This column is only essential here because the h and gamma are not included in the data.frame. It would also have been possible to include h and gamma columns in the data.frame and not include the k_vb column.
The values of the non-essential species specific parameters alpha, k, ks, z0, w_min and erepro are not included in the data.frame. This means that the default values will be automatically used when we create the MizerParams object.
Note that there are no columns describing the fishing selectivity. There is no sel_func column to determine the selectivity function. This means that the default selectivity function, knife_edge, will be used. As mentioned in the section on fishing gears, this function also needs another argument, knife_edge_size. This is not present in the data.frame and so it will be set to the default value of w_mat. Also, there is no catchability column so a default value for catchability of 1 will be used for all gears and species.
This species parameter data.frame is the minimum we need to create a MizerParams object as it contains only essential columns. To create the MizerParams object we pass the data.frame into the MizerParams() constructor method:
## Note: No sel_func column in species data frame. Setting selectivity to be 'knife_edge' for all species.## Note: No knife_edge_size column in species data frame. Setting knife edge selectivity equal to w_mat.## Note: No h column in species data frame so using f0 and k_vb to calculate it.## Note: No gamma column in species data frame so using f0, h, beta, sigma, lambda and kappa to calculate it.## Note: No z0 column in species data frame so using z0 = z0pre * w_inf ^ z0exp.## Note: No ks column in species data frame so using ks = h * 0.2.We have just created a MizerParams object:
## [1] "MizerParams"
## attr(,"package")
## [1] "mizer"As has been mentioned in the section on the community model and the section on the trait based model, a MizerParams object is made up of slots that store a wide range of model parameters. Each of these slots contains information on the parameters in the model. A description of these slots can be found by calling help() on the class: help("MizerParams-class"). The different slots can be accessed using the @ operator.
The slot species_params contains the species parameters data.frame that was passed in to the constructor. We can inspect this slot with:
## species w_inf w_mat beta sigma r_max k_vb gear k alpha
## Sprat Sprat 33 13 51076 0.8 7.38e+11 0.681 Sprat 0 0.6
## Sandeel Sandeel 36 4 398849 1.9 4.10e+11 1.000 Sandeel 0 0.6
## N.pout N.pout 100 23 22 1.5 1.05e+13 0.849 N.pout 0 0.6
## Herring Herring 334 99 280540 3.2 1.11e+12 0.606 Herring 0 0.6
## Dab Dab 324 21 191 1.9 1.12e+10 0.536 Dab 0 0.6
## Whiting Whiting 1192 75 22 1.5 5.48e+11 0.323 Whiting 0 0.6
## Sole Sole 866 78 381 1.9 3.87e+10 0.284 Sole 0 0.6
## Gurnard Gurnard 668 39 283 1.8 1.65e+12 0.266 Gurnard 0 0.6
## Plaice Plaice 2976 105 113 1.6 4.08e+14 0.122 Plaice 0 0.6
## Haddock Haddock 3485 165 558 2.1 1.84e+12 0.271 Haddock 0 0.6
## Cod Cod 40044 1606 66 1.3 8.26e+09 0.216 Cod 0 0.6
## Saithe Saithe 16856 1076 40 1.1 1.12e+11 0.175 Saithe 0 0.6
## erepro sel_func knife_edge_size h gamma z0
## Sprat 1 knife_edge 13 18.20276 3.190183e-11 0.18705957
## Sandeel 1 knife_edge 4 27.51606 1.503569e-11 0.18171206
## N.pout 1 knife_edge 23 32.83924 8.504085e-11 0.12926608
## Herring 1 knife_edge 99 35.03807 1.123212e-11 0.08647736
## Dab 1 knife_edge 21 30.67834 4.645215e-11 0.08735805
## Whiting 1 knife_edge 75 28.53952 7.390625e-11 0.05658819
## Sole 1 knife_edge 78 22.55847 3.115287e-11 0.06294752
## Gurnard 1 knife_edge 39 19.37727 2.948552e-11 0.06863713
## Plaice 1 knife_edge 105 14.62366 2.846500e-11 0.04171321
## Haddock 1 knife_edge 165 34.23910 4.037031e-11 0.03957464
## Cod 1 knife_edge 1606 61.58170 1.597275e-10 0.01753768
## Saithe 1 knife_edge 1076 37.39168 1.230565e-10 0.02340093
## ks w_min w_min_idx
## Sprat 3.640551 0.001 1
## Sandeel 5.503212 0.001 1
## N.pout 6.567848 0.001 1
## Herring 7.007614 0.001 1
## Dab 6.135668 0.001 1
## Whiting 5.707904 0.001 1
## Sole 4.511695 0.001 1
## Gurnard 3.875454 0.001 1
## Plaice 2.924733 0.001 1
## Haddock 6.847819 0.001 1
## Cod 12.316339 0.001 1
## Saithe 7.478336 0.001 1We can see that this contains the original species data.frame (with w_inf and so on), plus any default values that may not have been included in the original data.frame. For example, we can see that there are no columns for alpha and h and gamma etc.
Also note how the default fishing gears have been set up. Because we did not specify gear names in the original species parameter data.frame, each species is fished by a unique gear named after the species. This can be seen in the new gear column which holds the names of the fishing gears. Also, the selectivity function for each fishing gear has been set in the sel_func column to the default function, knife_edge(). A catchability column has been added with a default value of 1 for each of the species that the gear catches. An example of setting the catchability by hand can be seen in the section on the North Sea.
There is a summary() method for MizerParams objects which prints a useful summary of the model parameters:
As well as giving a summary of the species in the model and what gear is fishing what species, it gives a summary of the size structure of the community. For example there are \(100\) size classes in the community, ranging from \(0.001\) to \(4.4\times 10^{4}\) . These values are controlled by the arguments no_w, min_w and max_w respectively. For example, if we wanted 200 size classes in the model we would use:
So far we have created a MizerParams object by passing in only the species parameter data.frame argument. This means the interaction matrix will be set to the default value (see the section on predator/prey preference.) This is a matrix of 1s, implying that all species fully interact with each other, i.e. the species are spread homogeneously across the model area. For the North Sea this is not the case and so the model would be improved by also including an interaction matrix which describes the spatial overlap between species.
An example interaction matrix for the North Sea has been included in mizer as a .csv file. The location of the file can be found by running:
Take a look at it in a spreadsheet if you want. As mentioned above, to read this file into R we can make use of the read.csv() function. However, this time we want the first column of the .csv file to be the row names. We therefore use an additional argument to the read.csv() function: row.names.
The read.csv() function reads the data into a data.frame. We want the interaction matrix to be of class matrix so we need to make use of the as() function.
We now have an interaction matrix that can be passed to the MizerParams constructor along with the species parameter data.frame. To make the MizerClass object you just call the constructor method and pass in the arguments. We will use default values for the remainder of the arguments to the MizerParams() method.
This means that we only need to pass in two arguments to the constructor:
## Note: No sel_func column in species data frame. Setting selectivity to be 'knife_edge' for all species.## Note: No knife_edge_size column in species data frame. Setting knife edge selectivity equal to w_mat.## Note: No h column in species data frame so using f0 and k_vb to calculate it.## Note: No gamma column in species data frame so using f0, h, beta, sigma, lambda and kappa to calculate it.## Note: No z0 column in species data frame so using z0 = z0pre * w_inf ^ z0exp.## Note: No ks column in species data frame so using ks = h * 0.2.Note that the first argument must be the species parameters data.frame. The remaining arguments can be in any order but should be named. If we didn’t want to use default values for the other arguments we would pass them in to the constructor by name.
We now have all we need to start running projections. Before we get to that though, we’ll take a quick look at how different fishing gears can be set up.
In the above example, each species is caught by a different gear (named after the species it catches). This is the default when there is no gear column in the species parameter data.frame.
Here, we look at an example where we do specify the fishing gears. We take the original params_data species parameter data.frame that was read in above and bind an additional column, gear, to it. This gear column contains the name of the gear that catches the species in that row. Here we set up four different gears: Industrial, Pelagic, Beam and Otter trawl, that catch different combinations of species.
params_data_gears <- params_data
params_data_gears$gear <- c("Industrial","Industrial","Industrial",
"Pelagic","Beam","Otter",
"Beam","Otter","Beam",
"Otter","Otter","Otter")If you inspect the params_data object you will see a new column, gear, has been added to it.
We then make a new MizerParams object as before:
## Note: No sel_func column in species data frame. Setting selectivity to be 'knife_edge' for all species.## Note: No knife_edge_size column in species data frame. Setting knife edge selectivity equal to w_mat.## Note: No h column in species data frame so using f0 and k_vb to calculate it.## Note: No gamma column in species data frame so using f0, h, beta, sigma, lambda and kappa to calculate it.## Note: No z0 column in species data frame so using z0 = z0pre * w_inf ^ z0exp.## Note: No ks column in species data frame so using ks = h * 0.2.You can see the result by calling summary() on the params_gears object.
In this example the same gear now catches multiple stocks. For example, the Industrial gear catches Sprat, Sandeel and Norway Pout. Why would we want to set up the gears like this? In the next section we will see that to project the model through time you can specify the fishing effort for each gear through time. By setting the gears up in this way you can run different management scenarios of changing the efforts of the fishing gears rather than on individual species. It also means that after a simulation has been run you can examine the catches by gear.
Once the MizerParams object has been properly setup, it may be the case that one wishes put the system in steady state. Sometimes this can be done simply by running the model using project until it reaches steady state. However, this method is not sure to work, and there is a function called steady() that is more reliable. The function steady() must be supplied with a MizerParams object. It takes that MizerParams object, looks at the initial system state, computes the levels of reproduction of the different species, hold them fixed, and evolves the system until a steady state is reached (or more precisely, until the amount that the population abundances change during a time-step is below some small tolerance level). After this, the reproductive efficiency of each species is altered so that when the reproduction dynamics are turned back on (i.e., when we stop holding recruitment levels fixed), the values of the reproduction levels which we held the system fixed at will be realized. The steady function is not sure to converge, and the way it re-tunes the reproductive efficiency values may not be realistic, but the idea is to alter the other parameters in the system until steady() does arrive at a steady state with sensible reproductive efficiency values.
Now that we know how to create a multispecies model we shall discuss how to run a multispecies model.