Beginner level tutorial - how to parameterise a Mizer model

In this tutorial you will learn

You are encouraged to run the commands in this tutorial in your own RStudio session. Before you start, you should make sure that you have the latest version of the mizerHowTo package by running the command

devtools::install_github("sizespectrum/mizerHowTo")

You should run the above command even if you already installed the package before, because it is currently changing all the time and reinstalling it will give you the latest version. If you happen to already have the latest version then the above command will know and not bother to reinstall.

Now you can load the package with

library(mizerHowTo)

That will print some messages that you can ignore.

Step 1 - What data are typically needed to parameterise a mizer model?

Figure 1: Illustration of the parameters needed (green boxes) and calibrated (blue boxes) by Mizer. The need for data is hierarchical: a model can be setup and calibrated with the information in red: knowledge of the asymptotic size and observations of biomass and fishing. The calibration can be refined by adding further information in life history parameters and by using knowledge of Fmsy to calibrate the reproductive efficiency (blue). Additional refinement can be done by specifying the interaction matrix, theta (blue). Other parameters can be adjusted but they are rarely known accurately on a species-by-species basis (black).

parameters’ description
parameter description
w_inf asymptotic weigth
w_mat maturation weight (determines when 50% of the population has matured using a sigmoid function)
beta preferred predator/prey mass ratio
sigma width of the feeding kernel
R_max Beverton-Holt density dependence parameter
k_vb von Bertalanffy growth parameter
l25 length at 25% catch
l50 length at 50% catch
a coefficient for age to size conversion
b constant for age to size conversion
catchability fisheries efficiency
h maximum intake rate
k metabolism constant
ks metabolism coefficient
z0 background mortality coefficient
gamma search volume (obtained from beta and sigma)
w_mat25 weight at which 25% of individuals are mature
erepro reproductive output scalar

The multispecies model in mizer allows one to resolve species-specific differences in life history and feeding parameters that are important for modelling particular ecosystems.

Let’s start with the minimal amount of information needed to run a mizer model.

First you are going to need a data frame of species specific parameters to input into the newMultispeciesParams() function. This function requires at least three columns of parameters:

  • species name (species)
  • asymptotic weight (w_inf)
  • maximum intake rate (h) or von Bertalanffy growth parameter (k_vb)

Von Bertalanffy Growth Function

The growth curve of a fish is commonly reported as a von Bertalanffy growth curve. The von Bertalanffy curve describes the length of a fish as a function of two parameters; the asymptotic length \(L_\infty\) and the von Bertalanffy growth coefficient \(K\) (k_vb in Mizer): \[ L_{(t)} = L_\infty (1-e^{-Kt}) \] where \(t\) is the age. \(L_\infty\) represents the asymptotic length and \(1/K\) is approximately the age at maturation. Mizer, however, uses two other parameters, \(W_\infty\) (w_inf in Mizer) and \(h\), where \(W_\infty\) represents the asymptotic weight and \(h\) the size-scale maximum consumption rate. Fortunately, the two sets of parameters are related, and one can derive \(h\) and \(W_\infty\) from von Bertalanffy \(K\) and \(L_\infty\).

The parameter \(h\) represents the maximum consumption rate of a fish, scaled by size. If we assume that stomach fullness (the feeding level) is that same for fish of all sizes, then \(h\) is related to the coefficient of the initial growth rate \(A\) (blue slanted dashed line) as \(h \approx 4A\). The value of \(h\) varies between 4 – 400, with a geometric mean around 22 g\(^{-0.25}\)/year across all fish species (Andersen 2019, chap. 2 and 11). A more precise estimate can be obtained from von Bertalanffy \(K\) and \(L_\infty\) as: \(h \approx 4.75 K\cdot L_\infty^{0.75}\) (with \(K\) given in 1/years and \(L_\infty\) in cm).

The data frame is arranged as species by parameter, so each column of the parameter data frame is a parameter and each row has the values of the parameters for one of the species in the model.

You will also need an interaction matrix that defines the overlapping interaction between each species (default is set to 1).

Interaction matrix

Such matrix determines the availability of a prey to a predator. Values ranges from 0 to 1 and are assigned to each combination of predator/prey in the ecosystem (canibalism included). 1 is translated by full availability of a prey to a predator whereas 0 is none. The interaction value may vary due to different time overlap between prey and predator (e.g. diurne and nocturne animals) or different spatial overlap (e.g. depths and biomes).

All other parameters either have default values or will be calculated from the supplied parameters.

Let’s start with a small made up example assuming we have already found some parameters from FishBase or the literature. Here we consider two interacting species and a background resource spectrum (we will come back to this). Using the North Sea as an example, let’s start with the iconic Atlantic Cod and a key local prey species Sandeel.

To set up this model we can look on Fishbase (www.fishbase.se or use the rfishbase package) to find the asymptotic weights of these two species. A search on fishbase reveals some parameters from empirical von Bertalanffy growth curves which can be used to set up the life history parameters - w_inf and k_vb. Often von Bertalanffy relationships are based on length and not weight. The length parameters can be converted to weight using length-weight regressions

Length-weight regression

The asymptotic length and weight are related as: \[W_\infty = a \cdot L_\infty^b\] Species-specific conversion coefficients \(a\) and \(b\) are found in Froese (2006). Reasonable general values are: \(a = 0.01\ \mathrm{g/cm}^{-3}\) and \(b = 3\).

We find an asymptotic length of 132 cm and k_vb of \(0.2\)/year for cod in the North Sea here:https://www.fishbase.se/popdyn/PopGrowthList.php?ID=69&GenusName=Gadus&SpeciesName=morhua&fc=183

And for sandeel, 18.5 cm and k_vb = 0.4/year https://www.fishbase.se/popdyn/PopGrowthList.php?ID=37&GenusName=Ammodytes&SpeciesName=marinus&fc=402

We then need to convert these asymptotic lengths to weights using length-weight regression parameters, which also can be found on fishbase or in the literature.

For cod the estimation is \(w_{inf} = 24600\) g https://www.fishbase.se/popdyn/LWRelationshipList.php?ID=69&GenusName=Gadus&SpeciesName=morhua&fc=183

For sandeel the estimation is \(w_{inf} = 22.83\) g https://www.fishbase.se/popdyn/LWRelationshipList.php?ID=37&GenusName=Ammodytes&SpeciesName=marinus&fc=402

If w_inf or l_inf are not available you could use maximum observed sizes. These values are often systematically larger than estimates of asymptotic weight, and we recommend you check the literature or size-at-age data for your system to check reliability of fishbase estimates.

We also know from the literature that the preferred predator-prey mass ratio for Cod is approximately 100 and we might guess that it is approximately 10000 for Sandeel as they feed on prey much smaller than themselves.

Predator:prey mass ratio PPMR

The PPMR gives an indication of the prey size predators preferably feed upon. A small PPMR means that the predator feeds on prey close to their size (e.g PPMR of 10 equal a prey 10 times smaller than the predator) whereas a large PPMR means that the predator feeds on prey considerably smaller than their own size (e.g baleen whales). The preferred PPMR in Mizer is noted beta.

For simplicity, we will assume defaults for all other parameters and that both species occur in the same environment throughout their lives (interaction matrix = 1).

smallExample <- data.frame("species" = c("sandeel", "cod"),
                           "w_inf" = c(23,24600), 
                           "k_vb" = c(0.1,0.2),
                           "beta" = c(10000,100))
knitr::kable(smallExample, caption = "species' parameters")
species’ parameters
species w_inf k_vb beta
sandeel 23 0.1 10000
cod 24600 0.2 100 
smallInter <- matrix(c(1, 1, 1, 1), ncol = 2, 
                     dimnames = list(smallExample$species, smallExample$species))
knitr::kable(smallInter, caption = "species' interactions")
species’ interactions
sandeel cod
sandeel 1 1
cod 1 1

The smallExample shows the format of the data frame for the species parameters and the smallInter shows the format of the interaction matrix. Both of these are required as inputs for creating a mizerParams object. Many mizer users will collect this information in a spreadsheet prior to reading in the data into mizer.

To illustrate this next step we will use pre-existing set of species’ specific parameters for a previously published North Sea model, that used more detailed fisheries dependent and independent data to calculate species parameters.

These files are pre-loaded in this package. More information on how to set this up is here: https://sizespectrum.org/mizer/articles/a_multispecies_model_of_the_north_sea.html

North Sea species’ parameters
species w_inf w_mat beta sigma k_vb sel_func l25 l50 a b catchability gear
Sprat 33.0 13 51076 0.8 0.681 sigmoid_length 7.6 8.1 0.007 3.014 1.2953333 sigmoid_gear
Sandeel 36.0 4 398849 1.9 1.000 sigmoid_length 9.8 11.8 0.001 3.320 0.0651055 sigmoid_gear
N.pout 100.0 23 22 1.5 0.849 sigmoid_length 8.7 12.2 0.009 2.941 0.3138000 sigmoid_gear
Dab 324.0 21 191 1.9 0.536 sigmoid_length 11.5 17.0 0.010 2.986 0.9780000 sigmoid_gear
Herring 334.0 99 280540 3.2 0.606 sigmoid_length 10.1 20.8 0.002 3.429 0.1815000 sigmoid_gear
Gurnard 668.0 39 283 1.8 0.266 sigmoid_length 19.8 29.0 0.004 3.198 0.4625057 sigmoid_gear
Sole 866.0 78 381 1.9 0.284 sigmoid_length 16.4 25.8 0.008 3.019 0.3738333 sigmoid_gear
Whiting 1192.0 75 22 1.5 0.323 sigmoid_length 19.8 29.0 0.006 3.080 0.2426667 sigmoid_gear
Plaice 2976.0 105 113 1.6 0.122 sigmoid_length 11.5 17.0 0.007 3.101 0.1848333 sigmoid_gear
Haddock 4316.5 165 558 2.1 0.271 sigmoid_length 19.1 24.3 0.005 3.160 0.3015000 sigmoid_gear
Saithe 39658.6 1076 40 1.1 0.175 sigmoid_length 35.3 43.6 0.007 3.075 0.3930000 sigmoid_gear
Cod 39851.3 1606 66 1.3 0.216 sigmoid_length 13.2 22.9 0.005 3.173 0.2666675 sigmoid_gear
North Sea interaction matrix
Sprat Sandeel N.pout Dab Herring Gurnard Sole Whiting Plaice Haddock Saithe Cod
Sprat 0.7291292 0.0340844 0.0635452 0.3624155 0.2741698 0.1751558 0.2979556 0.2652592 0.3706598 0.0813555 0.0168132 0.3375797
Sandeel 0.0340844 0.6811988 0.0489243 0.0973666 0.0588821 0.0599265 0.0602086 0.0751001 0.0780185 0.0939573 0.0160902 0.0994345
N.pout 0.0635452 0.0489243 0.7966043 0.0908880 0.2975507 0.3062414 0.0167902 0.2998989 0.0785582 0.5491755 0.2949894 0.3250226
Dab 0.3624155 0.0973666 0.0908880 0.8081777 0.2896396 0.2204120 0.3804746 0.3338973 0.5649221 0.1316806 0.0313820 0.4164780
Herring 0.2741698 0.0588821 0.2975507 0.2896396 0.6589010 0.2751063 0.2001414 0.3737393 0.2779187 0.3483547 0.1262059 0.4047793
Gurnard 0.1751558 0.0599265 0.3062414 0.2204120 0.2751063 0.8801050 0.1067789 0.3710990 0.1649212 0.3573544 0.1235199 0.3518328
Sole 0.2979556 0.0602086 0.0167902 0.3804746 0.2001414 0.1067789 0.7155805 0.1922746 0.3913732 0.0344780 0.0124205 0.2576123
Whiting 0.2652592 0.0751001 0.2998989 0.3338973 0.3737393 0.3710990 0.1922746 0.7092823 0.2950381 0.3916479 0.1022817 0.4406088
Plaice 0.3706598 0.0780185 0.0785582 0.5649221 0.2779187 0.1649212 0.3913732 0.2950381 0.7192239 0.1124851 0.0329494 0.3504367
Haddock 0.0813555 0.0939573 0.5491755 0.1316806 0.3483547 0.3573544 0.0344780 0.3916479 0.1124851 0.8583072 0.2616747 0.3957734
Saithe 0.0168132 0.0160902 0.2949894 0.0313820 0.1262059 0.1235199 0.0124205 0.1022817 0.0329494 0.2616747 0.6638355 0.2089450
Cod 0.3375797 0.0994345 0.3250226 0.4164780 0.4047793 0.3518328 0.2576123 0.4406088 0.3504367 0.3957734 0.2089450 0.7865471

The second type of data you are going to need is a data set of catch and/or spawning stock biomass (\(SSB\)) of the selected species, so you can compare the model output’s to real data. The fisheries time-series of the North Sea are also available in this repository. We will use these data later in our second tutorial on calibration.

Because the North Sea is heavily fished we also need information on fishing intensity and other parameters have been entered into the species parameter file that relate to the type of fishing gear selectivity (which is assumed to be species-specific).

In mizer, fishing mortality rates at size for each gear are calculated as

\[F = catchability\cdot selectivity \cdot effort\]

Selectivity is determined using the l_25, l_50, and sel_func parameters but simpler approaches can also be used. See https://sizespectrum.org/mizer/reference/setFishing.html One effort is set during the simulation for each gear in the ecosystem, which is determined by gear. Finally, catchability is catchability.

For the North Sea we assumed \(catchability.effort\) could be estimated from the fishing mortality rates of fully selected sizes/ages of fish from single-species stock assessments.

data-raw/DATASET.R contains the code extracting the data from the ICES stock assessment database in a usable format. Fishing mortality data is averaged over 2014-2019 as it is a relatively stable period in catches and has the maximum amount of data across all species concerned. The catchAvg object (loaded with the package) contains the average catch data over 2014-2019 and is used later to compare model output versus empirical data.

Step 2 - How to convert the data into a valid mizerParams object

In this section you will:

  • Learn to format raw data into a Mizer compatible format

Inputing the previous data frame and interaction matrix into the newMultispeciesParams() function output a fully fleshed mizerParams object

param <- newMultispeciesParams(smallExample, smallInter)
## Note: No h provided for some species, so using f0 and k_vb to calculate it.
## Note: Because you have n != p, the default value is not very good.
## Note: No ks column so calculating from critical feeding level.
## Note: Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
## Note: Using f0, h, lambda, kappa and the predation kernel to calculate gamma.
species_params(param)
##         species w_inf k_vb  beta   w_mat w_min alpha interaction_resource
## sandeel sandeel    23  0.1 10000    5.75 0.001   0.6                    1
## cod         cod 24600  0.2   100 6150.00 0.001   0.6                    1
##                 n   p         q pred_kernel_type sigma         h k        ks
## sandeel 0.6666667 0.7 0.7166667        lognormal     2  2.126856 0 0.2407671
## cod     0.6666667 0.7 0.7166667        lognormal     2 45.821729 0 4.1110890
##                 z0        gamma     w_mat25 m erepro R_max
## sandeel 0.21098033 3.661644e-12    5.151761 1      1   Inf
## cod     0.02063033 1.002286e-10 5510.144528 1      1   Inf

You can see that mizer has estimated several missing parameters. [Go through each of these and take a look at figure 1]

Let’s do the same for our North Sea data set and look at what parameters we can find in the mizerParams object:

params_uncalibrated <- newMultispeciesParams(nsParams, interNS, kappa = 1e11,
                                             initial_effort = 1)
## Note: No h provided for some species, so using f0 and k_vb to calculate it.
## Note: Because you have n != p, the default value is not very good.
## Note: No ks column so calculating from critical feeding level.
## Note: Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
## Note: Using f0, h, lambda, kappa and the predation kernel to calculate gamma.

We set initial_effort = 1 as we have already set the catchability column of the species data frame to the estimated fishing mortality, and in mizer the fishing mortality of fully selected fish is the product of catchability and effort.

You will note that in the case of the North Sea, we provide another parameter kappa. kappa sets the carrying capacity of the ecosystem and is set to reflect the entire volume of the North Sea - hence the very large kappa value.

Background resource spectrum

In Mizer, a background resource spectrum is simulated to represent all the food sources that are not explicilty modelled as fish. It is especially useful to feed small-sized fish that are not piscivores. The carrying capacity of the background spectrum is determined by kappa and its regeneration rate by r_pp. High kappa implies more food and therefore faster growth rate and/or higher sustainable fish biomass.

mizerParams object
species w_inf w_mat beta sigma k_vb l25 l50 a b catchability h ks z0 gamma erepro R_max
Sprat 33.0 13 51076 0.8 0.681 7.6 8.1 0.007 3.014 1.2953333 14.46675 1.593753 0.1870596 0 1 Inf
Sandeel 36.0 4 398849 1.9 1.000 9.8 11.8 0.001 3.320 0.0651055 25.62741 2.936414 0.1817121 0 1 Inf
N.pout 100.0 23 22 1.5 0.849 8.7 12.2 0.009 2.941 0.3138000 31.20422 3.372902 0.1292661 0 1 Inf
Dab 324.0 21 191 1.9 0.536 11.5 17.0 0.010 2.986 0.9780000 34.87720 3.781368 0.0873580 0 1 Inf
Herring 334.0 99 280540 3.2 0.606 10.1 20.8 0.002 3.429 0.1815000 28.36363 2.920263 0.0864774 0 1 Inf
Gurnard 668.0 39 283 1.8 0.266 19.8 29.0 0.004 3.198 0.4625057 20.64990 2.193126 0.0686371 0 1 Inf
Sole 866.0 78 381 1.9 0.284 16.4 25.8 0.008 3.019 0.3738333 24.73805 2.567302 0.0629475 0 1 Inf
Whiting 1192.0 75 22 1.5 0.323 19.8 29.0 0.006 3.080 0.2426667 31.77220 3.301616 0.0565882 0 1 Inf
Plaice 2976.0 105 113 1.6 0.122 11.5 17.0 0.007 3.101 0.1848333 16.94072 1.740765 0.0417132 0 1 Inf
Haddock 4316.5 165 558 2.1 0.271 19.1 24.3 0.005 3.160 0.3015000 41.46028 4.196598 0.0368503 0 1 Inf
Saithe 39658.6 1076 40 1.1 0.175 35.3 43.6 0.007 3.075 0.3930000 59.95343 5.700788 0.0175943 0 1 Inf
Cod 39851.3 1606 66 1.3 0.216 13.2 22.9 0.005 3.173 0.2666675 69.40226 6.511735 0.0175659 0 1 Inf

Step 3 - Exploring the model

In this section you will:

  • Develop a feel for how the various parameters affect the steady state of your model.

The easiest way to start exploring the effect of the model parameters on its steady state is to use an interactive web application that is currently being developed in the mizerExperimental package. So we start by installing that package:

devtools::install_github("sizespectrum/mizerExperimental")

You should run the above command even if you already installed the package before, because it is currently changing all the time and reinstalling it will give you the latest version. If you happen to already have the latest version then the above command will know and not bother to reinstall.

Next we load the package with

library(mizerExperimental)

We now start the web app with the command

params <- tuneParams(params_uncalibrated)

This will open a new tab in a webbrowser with quite a lot on it. Among many others there is a button labelled “Instructions” at the top left. If you click that you will be guided through the main elements of the page. However in this workshop we will demonstrate the use of the tuneParams() app in person once the current breakout room session ends.

When you press the “Done” button in the app, the MizerParams object with all the changes you have made to it will be returned. You can now continue to work with it in the console. For example you could try

plot(params)

Step 4 - Tuning stock-recruitment

In this section you will:

  • Try to get coexistence between your species manually

First, let’s change the reproductive efficiency in the model to the same value for all species and then use the project() function with our mizerParams object to project forward in time and see if species manage to coexist together.

species_params(params_uncalibrated)$erepro <- .01
sim_uncalibrated <- project(params_uncalibrated)
plotCalibration(sim_uncalibrated, stage = 1)

The top panel shows the different species size spectrum at the last time step of the simulation while the bottom panel shows the abundance per species through time.

These plots show that species do not coexist and several go extinct. This is because there was no external density dependence (R_max is set at Inf) and the largest species (Cod and Saithe) are out-competing the rest.

Stock recruitment relationship

In Mizer, density dependence is modelled as a stock-recruitment relationship (SRR). The default is such that the recruitment approaches a maximum as the independendent reproduction rate (RDI) increases (similar to a Holling type II function response). \[RDD = R_{max}\frac{RDI}{RDI + R_{max}}\] RDD is the dependent reproduction rate that will determine the recruitement per species in the model.

Figure 2: Illustration of the interplay between the reproductive output and the external density dependence. The model simulates a size spectrum of each species (right). From that spectrum it calculates the reproductive output (blue arrow). The reproductive output is reduced by imposing an external density dependence in the form of a stock recruitment relationship (left panel; orange line). The recruitment then determines the flow of new recruits into the size spectrum. The simulation continues until the size spectrum reaches an equilibrium (which may be a dynamic equilibrium). Changing Rmax up or down by a factor results in a change up or down of the spectrum by the same factor. Rmax therefore functions as the central parameter that adjusts the total biomass of a species.

To get coexistence one needs to guess reasonable \(R_{max}\) values which will stop out-competition from a few species. We find that the density dependence is positively related to body size, meaning that large individuals will need to have a stronger density-dependence applied to them (and therefore a lower \(R_{max}\) value). Let’s start again and set some guessed R_max values.

Question: What does Rmax need to be to achieve coexistence? Experiment with the code chunk below.

params_guessed <- params_uncalibrated
# penalise the large species with higher density dependence
species_params(params_guessed)$R_max <- resource_params(params_guessed)$kappa * species_params(params_guessed)$w_inf^-1

sim_guessed <- project(params_guessed)
# saveRDS(sim_guessed, "../../data-raw/HTM1_sim.rds")
plotCalibration(sim_guessed)

The ecosystem looks way better. Saithe’s largest individuals are having a hard time, but at least species coexist.

\(R_{max}\) affects the density-dependent reproduction rate (\(RDD\)) by limiting the maximum amount of recruits surviving from the density-independent reproduction rate (\(RDI\)). Let’s look at the \(RDI/RDD\) ratio to see how strong \(R_{max}\) acts on our different species.

Is the physiological recruitment, \(RDI\), much higher than the realised recruitment, \(RDD\)? High \(RDI/RDD\) ratio indicates strong density dependence, meaning that the carrying capacity is controlling the population rather than predation or competition. Larger species often require more of this density dependent control than smaller ones. If \(RDI/RDD\) is too low, the efficiency of reproduction (erepro) can be lowered to ensure species do not outcompete others or are over-resilient to fishing. The largest species that were the most limited by our new \(R_{max}\) do not show a strong density dependence. The medium-sized species are the most affected here.

Is the physiological recruitment, \(RDI\), much higher than the realised recruitment, \(RDD\)? High \(RDI/RDD\) ratio indicates strong density dependence, meaning that the carrying capacity is controlling the population rather than predation or competition. Larger species often require more of this density dependent control than smaller ones.

If \(RDI/RDD\) is too low, the efficiency of reproduction (erepro) can be lowered to ensure species do not out-compete others or are over-resilient to fishing.

The largest species that were the most limited by our new \(R_{max}\) do not show a strong density dependence. This is because they are anyway hard hit by fishing in this model. The medium-sized species are the most affected here.

Step 5 - Checking assumptions

In this section you will:

  • Lookout for tell-tale sign of something wrong happening in your ecosystem

Let’s go through an example we have made previously to provide some tips of what to look for once you run the model.

The figure below shows an example of a calibrated ecosystem. It is obtained from plotSummary().

Figure 3: summary diagnostic

Let’s compare it to our own little ecosystem

plotSummary(sim_guessed)