In drfinlayscott/mizer: Dynamic Multi-Species Size Spectrum Modelling

knitr::opts_chunk$set(echo = TRUE)
library("mizer")
library("tidyverse")

We now consider the case where we know some life-history parameters of some important species in our ecosystem, and want to study these embedded in a full set of background species. We use a trait-based model for the background species and then add our foreground species with their specified parameters.

Let's look at hake and mullet for example for our foreground species. Let's assume we know the following parameters for them:

species_params <- tribble(
    ~species, ~w_inf, ~w_mat, ~beta, ~sigma, ~a,      ~b,    ~k_vb, ~l50, ~l25,
    "Hake",   10470,  250,    11,    1.1,    0.00667, 3.035, 0.178, 16.6, 16,
    "Mullet", 250,    16.5,   283,   1.8,    0.00624, 3.15,  0.6,   15.5, 13.2
)
species_params$sel_func <- "sigmoid_length"

Here w_inf is the asymptotic size, w_mat the maturity size, beta and sigma the parameters of the log-normal predation kernel, a and b the parameter in the length to weight conversion $w = a\, l^b$, k_vb the von Bertalanffy growth curve parameter and l50 and l25 the parameters of the sigmoid_length() selectivity function.

In practice you would probably not type them into R as above but create a spreadsheet with the parameters and then read that in with read_csv().

We now start by creating a scaling model for the background species.

params <- set_scaling_model(
    min_w_pp = 10^-12,
    no_sp = 16, 
    min_w_inf = 1, 
    max_w_inf = 10^6,
    min_egg = 10^-5,
    min_w_mat = 10^-0.6,
    w_pp_cutoff = 2 * 10^-0.6,
    no_w = 200,
    bmort_prop = 0.1) %>% 
    steady(t_max = 100) %>% 
    markBackground()

Note how we used the pipe operator %>% to pipe the result of set_scaling_model() into steady() and then the result from that into markBackground(). The mizer functions are pipe-friendly in that their first parameter is always the MizerParams object or the MizerSim object.

The call to markBackground() just specifies that these are background species. This has the effect for example that they are grey when plotted:

plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

and we will see other consequences below.

We now add our foreground species with the parameters specified in the data frame species_params.

params <- addSpecies(params, species_params)

plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

We see that addSpecies() adds the new species at a very low abundance compared to the background species. This is on purpose, because at that low abundance there is very little backreaction from the background species, and therefore we are in a regime where we understand the steady-state size distribution of the species analytically. Therefore mizer was able to find initial values for the size distributions that are close to steady state values. We see that little changes as we run the system to steady state:

params <- steady(params)
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

We make a copy of this for later purposes.

params2 <- params

Now we can increase the abundance of the background species.

params <- rescaleAbundance(params, factor = 16)
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

Note how rescaleAbundance() knew to only rescale the foreground species and leave the background species alone. Note also, that the introduction of the foreground species has led to a bulge in the total community spectrum. However the abundance of the background species was originally chosen so that the community abundance would be close to the Sheldon power-law spectrum. We therefore retune their abundance so that again the community spectrum is as close to the Sheldon power-law as possible.

params <- retuneBackground(params)
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

We see that as we introduce the foreground species, the background species become less important.

The above plot does not yet show a steady state. We still need to call steady() to find the steady state.

params <- steady(params)
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

Oh dear! The steady state looks quite different. The background species have increased a lot in their abundance. Perhaps the answer is to just retune their abundances and run to steady state again.

params <- params %>% retuneBackground() %>% steady()
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

It is getting worse! Clearly this iteration is unstable. This is because we changed the system by too big a step in one go when we rescaled the abundances of the foreground species by a factor of 16. We need to do this in smaller steps.

params <- params2 %>% 
    rescaleAbundance(factor = 2) %>% 
    retuneBackground() %>% 
    steady() %>% 
    rescaleAbundance(factor = 2) %>% 
    retuneBackground() %>% 
    steady() %>% 
    rescaleAbundance(factor = 2) %>% 
    retuneBackground() %>% 
    steady() %>% 
    rescaleAbundance(factor = 2) %>% 
    retuneBackground() %>% 
    steady()
plotSpectra(params, power = 2, total = TRUE, 
            ylim = c(10^-6, NA), wlim = c(10^-4, NA))

We can rescale the foreground species to whatever abundance we like. So we need some observation to choose the right abundance. In our system we have for both of our foreground species observations of the total biomass of all individuals above a certain cutoff size. We put that information into the species_params data frame.

species_params(params)$biomass_observed <- NA
species_params(params)$cutoff_size <- NA
species_params(params)["Hake", "biomass_observed"] <- 0.013
species_params(params)["Hake", "cutoff_size"] <- 0.3
species_params(params)["Mullet", "biomass_observed"] <- 0.0018
species_params(params)["Mullet", "cutoff_size"] <- 0.25

Now I want to introduce a very convenient gadget that allows you to tune the parameters of the system interactively. This is not yet very well documented, so I will demonstrate.