sizespectrum / mizer Goto Github PK

It is suggested that starvation mortality should be added, as indicated in the original mizer vingette, and in Hartvig and Andersen (2013, Theoretical Population Ecology).

This is now added as a pull request in the mizer development branch.

The starvation mortality is calculated at each time step and each size class and is defined as:
0, if e (assimilated intake - maintenance) >= 0

• e/ 0.1*m, if e <0
  The value 0.1 is just a fixed parameter, following Hartvig and Andersen (2013), used a proportionality constant or reserve proportion in total weight (assumed to be fixed for all size classes). Such formulation leads to instantaneous starvation mortality of 1 if energy deficit is 10% of the total body weight.

Currently in mizer negative e values are converted to 0 to ensure there is negative growth.
In the version with starvation mortality, this is still the case, but negative e values add an extra starvation mortality term

@Kenhasteandersen , thank you for adding the example parameter objects. I have added their help pages to the website at https://sizespectrum.org/mizer/dev/reference/index.html#section-example-parameter-sets. In the new version of mizer, parameter objects also carry initial state values. Would it be alright with you if I set those initial values to the steady state values?

Currently the handling of the colour and linetypes used in mizer plots is a mess. There are some good aspects to it, but they are undocumented and hidden in the mess. It would be really good if someone could clean this up.

Here is a description of how things are handled at the moment:

There are two slots in a MizerParams object:

• @slot linecolour A named vector of colour values, named by species or resource or other type (e.g. "fishing"). Used to give consistent colours in plots. Used not only for line colours but also for fill colours.
• @slot linetype A named vector of linetypes, named by species or resource or other type. Used to give consistent line types in plots.

When a new MizerParams object is created with emptyParams(), these slots get initialised with hard-coded values, except where the user included columns linetype and linecolour in the species_params data frame. This happens in lines 575 and following in the file MizerParams-class.R.

The mizer plots, which are created with ggplot2, then use the @linecolour slot via scale_colour_manual() and scale_linetype_manual().

The @linecolour and @linetype slots are also modified in the functions addSpecies(), removeSpecies() and renameSpecies(), for obvious reasons.

Linetypes and linecolours for unstructured resources are currently set in setResourceDynamics(), but because the code for unstructured resources is going to be taken out of mizer again, to be placed into a mizer extension, you can ignore that code. But the message is that extensions will need to be able to add extra entries into the linecolours and linetypes slots or whatever takes their place in your improved scheme.

There are two very desirable goals for any improvement:

• The setting of line colours and types should be documented

• The default colour palette should not be hardcoded but changeable by the user, possibly by setting an options variable.

In addition, it would be nice to make the code a bit more elegant and robust.

Equation (3.2) in the vignette gives an incorrect description of the boundary condition due to reproduction. It says that the density N(w_0) of individuals at w_0 is fixed in terms of R and g(w_0). In reality the amount R of reproduction only contributes to the rate of change in the density. This is correctly implemented in the code.

I do not have a good suggestion of how to change the vignette, because the way reproduction is actually implemented is difficult to briefly describe in the continuum formulation used in that section of the vignette. In reality R is spread out over the entire first weight class and contributes to the rate of change of the number of individuals in that weight class. This is difficult to describe in that part of the vignette because it does not yet introduce the weight classes. I would probably propose to just take out the equation and just describe in words how reproduction is implemented. Making a change is a good idea though, because I got really confused by the discrepancy between the description in the vignette and the actual implementation. It makes a big difference to the solution whether R fixes the boundary value or whether it only affects its rate of change.

Just a reminder about stochastic recruitment and primary production, as discussed

I think it would be very useful to have a page on our website that helps people who already know Ecopath with Ecosim to understand mizer.

Luckily, Peter West (@peterwest17), a Master student at the University of the Algarve, is doing his masters project on a comparison between EwE and Mizer by creating a mizer model for an ecosystem for which there is already an Ecopath model, described in Gamito, Sofia, and Karim Erzini. “Trophic Food Web and Ecosystem Attributes of a Water Reservoir of the Ria Formosa (South Portugal).” Ecological Modelling 181, no. 4 (February 2005): 509–20. https://doi.org/10.1016/j.ecolmodel.2004.02.024. The project is supervised by Karim Erzini and myself.

Peter is currently working on relating the quantities and rates in the Ecopath model to integrals over the size-resolved quantities and rates in mizer.

We need to change the naming of r_max -> R_max, to avoid confusion with the population growth rate (r_max).

This should happen before the release of version 1.1. This functionality will then later go into a mizer extension.

We have started to built a website with the mizer documentation. You can view its current state at http://gustavdelius.github.io/mizer/ . There is clearly much scope for improvement, but I think it is already useful. Do we want that kind of website for the official mizer?

We used the pkgdown package to create the website. This means it is created from the content in the github repository itself, making it easy for people to contribute. I have started to document how to develop the website at http://gustavdelius.github.io/mizer/articles/editing_website.html. There is also a "Edit this page" link at the top of most pages which allow anyone to suggest edits (this works via github pull requests).

Hi Finlay, you indicated that you wanted someone to adopt mizer. I have created an 'organisation' called 'sizespectrum' and have invited you as an owner. You could transfer the mizer project to that organisation and in that way you can share the maintainership of mizer with others, including me. What do you think? For how to effect a transfer, see https://help.github.com/articles/transferring-a-repository/

The colour palette we are currently using is taken from http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette because it works well for colour-blind people. It looks good for area plots, for example:

but the palette is quite horrible for line plots, for example:

Does anyone know of a good colour palette to use for line plots?

We are planning to have a large number of vignettes on the mizer website. Some of them will be quite large due to the inclusion of screenshots and plots. We do not want to include all of them in the package, because that would inflate the download size of the package. We expect most users to prefer to look at the online documentation anyway, rather than the vignettes downloaded with the package.

So the issue to be discussed here: which vignettes, if any, do we want to include in the package?

The next question is whether we want to include a vignette in pdf format. I think that would be quite nice for users who still like to print out things. However It will be a pain to maintain the pdf vignette and the website separately. So it would be nice if the pdf vignette could be created from selected vignettes on the website. In principle it is possible to combine vignettes and render to pdf using the method described at https://stackoverflow.com/questions/25824795/how-to-combine-two-rmarkdown-rmd-files-into-a-single-output . However the result does not look perfect and the cross-linking does not work. A better option is to wait for the resolution of r-lib/pkgdown#853

In order to be able to model detritus and carrion in shallow soft-bottom fishing grounds, we are adding unstructured resource components to mizer. The currently proposed design of the new feature is described at http://gustavdelius.github.io/mizer/reference/resource_dynamics.html. This is a call for comments.

I want to apply this model to an ocean ecosystem. I try to evaluate the response of some important economic species to fishing. Most species in the simulation study are pelagic predators. Data about the background resources (initial_n_pp) in the open ocean is difficult to obtain and I don't know how to estimate this parameter (carrying capacity). Many thanks!

Does anyone remember the purpose of the print_it = TRUE option in the mizer plot functions? It has the effect of triggering an print(p) where p is the plot. Because the function also visibly returns the plot, in interactive use the plot gets printed anyway, even with print_it = FALSE and with print_it = TRUE it gets printed twice.

Any objections to removing the print_it argument from all plot functions in order to not confuse new users and to simplify?

Currently the function that is used to project the plankton abundance forward by one time step is stored in the @plankton_dynamics slot of the MizerParams object as a function. This leads to problems when the params object is shared with others, see #91 . We should therefore change to saving the name of the function instead of the function itself. Presumably it would then also be a good idea to change the plankton_dynamics argument to the setPlankton() function to a string rather than a function.

Currently a misformed species parameter dataframe can lead to some unhelpful error messages.

Mizer is really mostly unit agnostic, i.e., the user gets to choose the units they like best when they calibrate the system. However the defaults are set up with grams and years in mind as weight and time units and in the newer versions of mizer some plots put these units on the axes.

The current situation is a bit confusing for users, as Chris Griffiths has just pointed out to me.

Should we make the situation less confusing by imposing the use of grams for sizes and for biomass and years for time? Or do we want to let the user make the choice? In the latter, do we want to store the choice of units in the params object so that the plot functions can label the axes appropriately? I am leaning towards the simpler solution of imposing units. Then the documentation of functions like getYield can be made clearer by saying that it returns the yield in grams per year.

A related issue is the question of whether all rates and abundances should be seen as per volume or for the entire study area. Here I think different users might have different preferences. I guess specifying rates and abundances per volume is a bit problematic because one needs to decide to which depth to go when determining the volume of the study area. Any comments?

One thing that I found really difficult to figure out when I was new to mizer was how to change the parameters in a model after I had set it up. Sometimes changing the species parameters in the MizerParams object would work, sometimes it would have no effect. In order to save future users from similar confusion, I would like to introduce a function changeParams() that users can call after changing a species parameter to change the model accordingly. So for example they could do:

params <- set_trait_model()
params@species_params$gamma[3] <- 1000
params <- changeParams(params)

and the result would be a model with an appropriately changed search volume.

Another thing that I found difficult to find out without looking very closely at the code documentation was whether some particular assumption about the form of some model function, like the allometric scaling assumptions for example, was hardwired into mizer or was something I could overrule. I therefore gave the changeParams() function additional arguments that allow the user to set up the most general mizer model possible.

Please take a look at http://gustavdelius.github.io/mizer/reference/changeParams.html to see the details of what I am proposing. I have split the overall changeParams() function into several smaller functions changing particular aspects of the model. I may not have made very good choices for the names of the functions, so any suggestions for changes will be appreciated.

Besides the benefit to new users, the introduction of these functions also makes the code much easier to work with, because they do not need to wade through the complicated old MizerParams() code to find where a particular slot was set up but can go directly to the appropriate function, because MizerParams() also uses these functions to fill the slots in the MizerParams object. Furthermore it is easier to know where to document the setting of the various parameters, because this can now go into the documentation of the appropriate function.

When calling the summary() function on a MizerParams object, I think the number of size bins in the plankton spectrum should refer to the number of non-empty bins (i.e. the number of bins containing plankton), not the length of the complete size spectrum w_full.

See patricksykes@7b8ebbc for how I think this might be improved.

Currently, if the coefficient ks of the allometric metabolic rate is not given in the species parameters, it is set to 0.2 * h. Given that the coefficient of the allometric growth rate is alpha * feeding_level * h - ks, species will not grow unless alpha * feeding_level > 0.2. Would it perhaps be more appropriate to set ks to 0.5 * alpha * feeding_level * h so that only half of the consumption is lost to metabolism? Any advice @Kenhasteandersen @juliablanchard @astaaudzi ?

It is known that if a stock falls below a certain threshold (5-10% of virgin biomass) its recovery often can be very slow or non-existent, because Allee effects or depensation start operating at very low densities. This would not be the case with the current assumptions in mizer, because even one fish left in a population could produce millions of eggs and they all would recruit bringing abundance up rapidly again. This means that we may overestimate population recovery rates and impacts of intensive fishing.

I am wondering whether it would make sense to allow users to setup their own SR functions, which would overwrite current SR in mizer? This way users could explore a wider range of scenarios. No need to modify the core mizer assumptions, just add flexibility

In the vignette, in section 5.4, there is an example of a trophic cascade. In the explanation it says: "...leading to an increase in their abundance. This in turn increases the predation mortality on their smaller prey, which reduces their abundance, and so on."

From this explanation one would expect the distance between the peak and the trough to be roughly equal to beta, the preferred predator-prey size ratio. However in Figure 6 we can see that that is not the case. The value of beta for that figure is beta=100, but the distance between peak and trough is only about 10.

I then changed beta to 10, expecting to see at least a decrease in the distance between peak and trough, but it stays roughly the same, as the attached image shows. The dotted line is that for beta=10. If anything, the distance between peak and trough has become wider by decreasing the predator/prey ratio, which appears to be in conflict with your explanation.

Do you know what is going on?

Here is the code used to produce that image:

params_knife <- set_community_model(z0 = 0.1, recruitment = 4e7, 
                                    alpha = 0.2, f0 = 0.7, knife_edge_size = 1000)
sim0 <- project(params_knife, effort = 0, t_max = 50)
sim1 <- project(params_knife, effort = 1, t_max = 50)
relative_abundance <- sim1@n[51,,] / sim0@n[51,,]
plot(x=sim0@params@w, y=relative_abundance, log="x", type="n", 
     xlab = "Size (g)", ylab="Relative abundance")
lines(x=sim0@params@w, y=relative_abundance)
lines(x=c(min(sim0@params@w),max(sim0@params@w)), y=c(1,1),lty=2)
abline(v=1000, lty=2)
abline(v=10, lty=2)

params_knife <- set_community_model(z0 = 0.1, recruitment = 4e7, beta = 10,
                                    alpha = 0.2, f0 = 0.7, knife_edge_size = 1000)
sim2 <- project(params_knife, effort = 0, t_max = 50)
sim3 <- project(params_knife, effort = 1, t_max = 50)
relative_abundance2 <- sim3@n[51,,] / sim2@n[51,,]
lines(x=sim0@params@w, y=relative_abundance2, lty=3, lwd=2)

While equation (3.13) in the mizer vignette says that when food supply is insufficient to cover the cost of metabolism and movement there is starvation mortality, no such starvation mortality is implemented in the mizer code. Shall we just state clearly in the vignette that it is not implemented?

Mizer parameters are usually calibrated to a particular system fished with a particular fishing effort. I think it would make sense if that information about the initial fishing effort that was used to calibrate the model was saved in the MizerParams object.

@astaaudzi has requested a function that returns the diet of a predator. I believe what she wants is a function like the following:

#' Get diet of predator resolved by prey species
#'
#' Calculates the rate at which a predator of a particular species and size consumes 
#' biomass of a prey species, where for this purpose we treat the plankton like another
#' prey species. Thus it performs the same integration as \code{getAvailEnergy()} but
#' does not aggregate over prey species, and multiplies by (1-feeding_level) to get the
#' consumed biomass rather than the available biomass.
#'
#' @param params A MizerParams object
#' @param n An array (species x size) with the abundance density of fish
#' @param n_pp A vector with the abundance of plankton
#' 
#' @return An array (predator species  x predator size x (prey species + plankton) )
#' @export
getDiet <- function(params, n, n_pp) {
    ...
}

Mizer should have a function for plotting the yield curve. @Kenhasteandersen mentioned that he found the yield curve to be useful while calibrating a model. Ken, do you have code to share for plotting the yield curve in mizer?

Mizer has a function getSSB(sim) that according to its documentation "Calculates the spawning stock biomass (SSB) through time of the species in the MizerSim class. SSB is calculated as the total mass of all mature individuals."

So it should multiply the total abundance by the maturity ogive (giving the proportion of individuals at size that are mature) and integrate over sizes. In practice however, it multiplies not by the proportion of individuals that are mature but by the proportion psi of the energy available for growth and reproduction that is put into reproduction. Thus it underestimates the spawning stock biomass quite substantially.

I propose that we split psi into its two factors, the maturity ogive giving the proportion of mature individuals (a sigmoidal function) and the reproduction allocation function that gives the proportion of the energy that a mature individual invests into reproduction (a power law). Then we can calculate the SSB properly which will make it easier to compare the mizer model output to observations of SSB.

I am posting this here in case I misunderstood what people mean by SSB.

Often one will want to combine the result of two different simulations on the same plot, for example to compare the consequences of different fishing regimes. This is not easy to achieve with mizer's existing plot functions. I have played with two different approaches to solving this problem:

1. Allowing an optional second MizerSim object to be passed to a plotting function, see plotYield().
2. Creating functions that return the dataframe needed to produce the plot, rather than the plot itself, see getBiomassFrame() and getSSBFrame(), and then providing a function displayFrames() that takes two of these frames and plots them.

Neither of these are perfect, so please make your own suggestions.

We should probably settle on one approach and then implement it consistently for all plotting functions where it is appropriate.

We should add a mizer_version slot to the MizerParams object which holds the version of the mizer code with which the params object was created. That way the results from simulations performed with that model will remain reproducible in the future by just checking out that version of the mizer code from github. It would also allow us to selectively enable new features of mizer only for params objects that support them, or to upgrade an old params object to a new version of mizer. Thanks to @astaudzi for bringing up this issue.

We might also want to add a repo and a ref slot, so that users can indicate that a params object is to be used with code in another repository or from a particular commit or tag. This can then be used together with install_github()

While we are at it, how about also adding a title and a description slot, to make it easier for people to remember what their params objects were for and a creation_time slot because why not and it can be populated automatically and may help people later to know which is their latest version.

Because at a later stage we might want to start a repository of mizer models contributed by different users, I also propose to add an author and an orcid slot, the former for a human-readable author name, the latter for their ORCID id.

Finally it would be good to know whether a params object is a community model, a trait-based model, a multi-species model or some other model for which setup functions may be developed in the future. We could store this info in two ways: either in a setup_fun slot that holds the name of the setup function, or in the class of the object, if in addition to the MizerParams class we introduce subclasses CommunityModel, TraitModel, MultispeciesModel, etc. I don't see a benefit of the latter approach, so I think I am leaning towards the more straightforward former approach.

The fact that the params object has a slot that contains a function creates a problem when we save params objects to disk and then share them with each other. To see this, take a look at the stock recruitment function stored in the @srr slot of one of the params objects that Ken recently added to mizer:

Barents_params@srr

This produces:

<srcref: file "/Users/ken/Documents/Source/Mizer/R/MizerParams-class.R" chars 2391:20 to 2396:1>
<environment: namespace:mizer>
Warning messages:
1: In getSrcLines(srcfile, x[7L], x[8L]) :
  restarting interrupted promise evaluation
2: In getSrcLines(srcfile, x[7L], x[8L]) :
  internal error -3 in R_decompress1

Notice that it points to code on Ken's computer and hence R produces some cryptic warning messages.

Note also that the @srr slot contains a closure (http://adv-r.had.co.nz/Functional-programming.html#closures), i.e., it contains its environment, which in this case is the entire content of the mizer namespace. I don't know if that is a problem, because I have not really gotten my head around closures. But it sounds like it might be a problem when sharing params objects. It does for example lead to inflated file sizes when saving the params objects to disk.

The easiest way to avoid these complications is to not store the function in the @srr slot but instead a string with the name of the function. This would be similar to how we currently specify the gear selectivity functions and the predation kernel type.

Another option would be to take the option of using alternative stock-recruitment relationships out of the main mizer package and provide it via a mizer extension and in the process design it more carefully.

Both these options would break backwards compatibility, but I don't think many people have used alternative stock-recruitment relationships yet.

All comments welcome.

Currently there are small discrepancies between the numerical results for the integrals over the feeding kernel using the old method and using the FFT method. As @astaaudzi has noticed, these differences can accumulate quite substantially over longer simulations.

There are two sources of discrepancies that we should remove:

1. The kernels used should be truncated in the same way (thanks to @juliablanchard for pointing this out)
2. Both methods should use the same discretisation scheme for the Riemann sums.

We propose to allow the user to input fishing gears and parameters as an additional input to the MizerParam function. The data frame will contain a species object, where the user specifies one species to be caught by a particular gear; a gear object, where the user specifies a gear name; a catchability function, where the user specifies the function used to generate the catchability curve; maximum catchability, defining the maximum catchability of that species with that gear, and parameters of the catchability curve. The same species can appear in different rows of the data frame with different gears.

The mizerParam class has a slot called "activity". What is its meaning? I can not find a place in the code where it is actually used.

It is often interesting to watch the sizespectrum evolve over time. So the plotSpectra() function should be extended to produce animations. The following code shows how easy it is to create animated size spectra with plotly:

library(mizer)
library(reshape2)
library(plotly)
# Create simulation
NS_params <- set_multispecies_model(NS_species_params)
sim <- project(NS_params, t_max = 20, effort = 1)
# Convert simulation result into a tibble
nf <- melt(sim@n)
# Add Plankton
n_ppf <- melt(sim@n_pp)
n_ppf$sp <- "Plankton"
nf <- rbind(nf, n_ppf)

plot_ly(nf) %>%
    # show only part of plankton spectrum
    filter(w > 10^-6) %>% 
    # calculate biomass density with respect to log size
    mutate(b = value * w^2) %>% 
    # Plot lines
    add_lines(
        x = ~w, y = ~b,
        color = ~sp,
        frame = ~time,
        line = list(simplify = FALSE)
    ) %>% 
    # Use logarithmic axes
    layout(xaxis = list(type = "log", exponentformat = "power"),
           yaxis = list(type = "log", exponentformat = "power"))

There is a desire, expressed by @juliablanchard and @astaaudzi to have convenient access to the full feeding kernel array (predator species x predator size x prey size) even when the model uses a constant predator/prey mass ratio, so that the information in the array is highly redundant. Having this array simplifies the calculation of diets, for example.

The proposal is to introduce a parameter to the functions for setting up MizerParams objects that is a flag determining whether to keep the full feeding kernel around. By default this is set to TRUE when the number of size classes is no more than 200 and to FALSE otherwise. The reason to set the default to FALSE when using more than 200 size classes is that the array then takes up substantial amount of storage.

Currently mizer uses S4 classes for the MizerParams and MizerSim objects, instead of the more usual S3 methods. S4 classes are more strictly enforcing the structure of the object. This has the unfortunate consequence that if a user saves a MizerParams or MizerSim object with one version of mizer, then it can not be loaded with a later version of mizer if in that version the MizerParams class has gained extra slots, for example. We can anticipate having to add extra slots in the future for new mizer features.

Slots in S3 objects are accessed via $ whereas slots in S4 objects are accessed via @. So currently we have syntax like params@species_params$w_mat because the MizerParams object is S4 and the species parameter data frame is S3. If we switch from S4 to S3 for everything, the syntax would change to params$species_params$w_mat. The effect on the user will be that they need to replace all @ signs in their code (that are not preceeded by a whitespace) with $ signs. Luckily we would be able to give them a regular expression with which they can make the change with one global search and replace. The advantage is that from then on the user will never again have to think about whether they need to use @ or $.

The mizer plot functions can be sometimes a bit slow, especially when running a model with a large number of size classes. One way to speed this up would be to skip some size classes to only plot enough points to obtain a smooth plot.

In addition to the multiple unstructured resource components discussed in #46, it would be good to also allow multiple size-structured resource components.

Asta @astaaudzi wrote on this topic in #46: I think that having several background spectra, both size structured and/or unstructured would definitely be useful and needed, especially in multi-species models. A parameter defining availability of these spectra to consumers is also required. How to setup up these parameters in the most efficient and clear way should be discussed. In my fork of mizer I have added two additional size structured background resources (which now has plankton, benthos and algal spectra) and introduced three additional parameters in the species parameter file, defining availability of these resources to consumers (avail_pp, avail_bb, avail_aa). As plankton and benthos resources have different minimum and maximum sizes this produces a natural ontogenetic diet transition from plankton to benthos and then to fish, if a species is a predator. For planktivores and herbivores the availability of benthos is set to 0 and for benthivores the availability of plankton is small, so they switch to benthic food when it fits into the predation kernel.

In this way the availability parameter could be added into the interaction matrix, because it is just a scalar. The rate with which the resource is consumed would then be determined by the standard gamma and h parameters.

Extending the n_pp slot of the MizerParams class to hold multiple spectra is straightforward, as is extending the calculation of the contribution of these to encounter rates and of the mortality rates for them. What needs discussion is how to let the user specify the dynamics for these components. Would it be sensible to follow the same route as proposed for the dynamics of unstructured resources of allowing the user to specify a function for each size-structured resource that implements the dynamics for that resource?

In that case the proposal would be to extend the slot resource_dynamics, which currently holds a list with the dynamics functions for each of the unstructured resources to also hold the dynamics functions for each size-structured resource. The default dynamics function would be the semichemostat function currently used for plankton, but the user could write their own more sophisticated functions. Similarly, the list in the slot resource_params will also have entries for the parameters needed by these functions. The current plankton-specific slots rr_pp and cc_pp would be moved to that list.

However, given that unstructured resources are quite a different beast from size-structured resources, maybe it would be less confusing to have separate slots structured_resource_dynamics and unstructured_resource_dynamics, but perhaps with shorter names if we can think of any.

Since its inception, mizer had three "wrapper functions" -- functions for setting up MizerParams objects with default values: set_community_model(), set_trait_model() and MizerParams(), for consistency also available as set_multispecies_model(). Each of these chooses default values for parameters that the user does not specify. Now we have become unhappy with several of those choices for defaults.

1. We changed the defaults for h, the coefficient of the maximum intake rate, to ensure that fish reach their maturity size at the correct age.
2. We changed the default for ks, the coefficient of basic metabolic rate, to ensure that a species stops being viable at a critical feeding level of 0.2, see #101 .
3. I am proposing to change the default for m, the exponent of the proportion of available energy that is invested into reproduction to better match the mizer growth curve to the von Bertalanffy growth curve, see #108 .
4. @Kenhasteandersen is proposing changes to several other defaults to match those more carefully chosen in his new book.

Currently these changes are turned on only when the user sets options(mizer_new = TRUE). That mechanism however can lead to confusion, as we experienced in Lysekil recently. So please make recommendations for how you think we should keep backwards compatibility for these wrapper functions. Or argue in favour of dropping backwards compatibility for the wrapper functions.

Note that is a completely separate issue from the backwards compatibility of the project() method. The issue discussed here only affects the setting up of MizerParams objects.

There are important cases where species are expected to have variable PPMR with size. This applies to planktivorous fish, that feed on small prey even when they grow large. Same is true for many coastal benthivorous fish, which are dominant species in e.g. Tasmanian coastal model. This point was made by @astaaudzi in a comment in issue #40

A particular choice for a feeding kernel with variable PPMR is published in Heneghan et al. 2016 (http://journal.frontiersin.org/article/10.3389/fmars.2016.00201/abstract) and only requires one extra parameter m. The code is more or less implemented in https://github.com/astaaudzi/mizer-rewiring (although the final version has not been pushed yet). We are just waiting for the resolution on the differences between different integration method for the feeding kernel (issue #39)

The way this should be added to mizer is to provide a new setup function that populates the feeding kernel slot in a MizerParams object with this new feeding kernel. No other changes to the mizer code are needed for this.

When MizerParams() calculates the size-dependent rates, it makes allometric assumptions for many rates. For example the metabolic loss scales as size to the power p, maximum intake rate scales as size to the power n and search volume scales as size to the power q. Currently these parameters are specified as arguments to MizerParams() (and its relevant constituent functions, like setMetab(), setIntakeMax() and setSearchVolume()). Hence the same values are used for all species.

Given that there is some evidence for variation in the scaling exponents across species and that the scaling exponents have quite a large effect on the shape of the within-species size distribution, I propose the option to specify the values for the exponents n, p and q via columns in the species parameter data frame, which would allow different values for different species if desired. Of course the default would not change and still be the same value for all species.

It would be good to have a flexible way to setup multiple resource spectra and availability of these spectra to size structured groups in the model. The numbers in the spectra could be added to the numbers in the simulated size structured groups, turning the number-at-size arrays (@n or @n_pp) into a list. Each resource spectrum can have its own parameters of kappa, lambda and regeneration rate (r_pp in case of current setup). The availability of these spectra to size structured groups could be set in a separate dataframe, in a similar way as the current interaction matrix

This is just to let the mizer community know that we have been awarded a grant by the EU to study fisheries discards and as part of the MINOUW project (http://minouw.icm.csic.es/). In particular it says in the grant agreement:

"Development of size-spectrum models to describe the effect of banning discards on structure and resilience of marine ecosystems. These models do the biomass bookkeeping of fish as they eat and grow. They will allow the mass balance of exploited and unexploited fish species to be computed, in the presence and absence of discarding, showing the changes in structure and resilience that emerge. This work will be done on the platform of 'mizer' (publicly available R code for modelling multi-species, size-spectrum dynamics in the marine environment) (Scott et al. 2013), modified for MINOUW to examine effects of discards. The computational task will be to redevelop mizer currently at TRL 4 for
high performance at TRL 7, demonstrated in the case study areas."

The details of exactly what we will do are yet to be developed. The work will be led by Richard Law and myself at the University of York and we are looking for a two-year postdoc to do the modelling work.

The mizer code is currently differing from the mizer vignette is some places. One of them is that the mizer code subtracts both the energy needed for metabolic activity $k_{s.i} w^p$ and the energy needed for other activity (movement) $k w$, whereas the model presented in chapter 3 of the vignette only has the former. Is there an argument why the energy cost of movement scales linearly with weight?

Asta (@astaaudzi) suggested in #46 that we might like to have different assimilation efficiencies for different food sources. Currently mizer does not allow that, having a species parameter alpha[i] that depends on the predator species i only and does not depend on the prey species. The proposal is to introduce a slot alpha in the MizerParams class for a matrix alpha[i, j] depending on both predator and prey species.

Unfortunately there is an efficiency problem with this idea: due to the way satiation is handled in mizer it would require summing over prey species twice. Mizer first calculates the encountered food. This involves summing over all prey species, independently of how well the prey species biomass can be assimilated. That rate of food encounter is then used to calculate the feeding level that determines what proportion of the encountered food is actually consumed by a predator. If we now wanted to implement a prey-specific assimilation rate we would have to re-do the sum over all prey species j and this time multiply by alpha[i, j]. Not horrible, but a slightly larger change than we might be prepared to accept.

If we decided to go ahead with this we would still want to make it convenient for the user to set up the assimilation efficiency matrix via the species parameter data frame if they want to. One proposal would be to allow two species parameters alpha_prey and alpha_pred and then to set alpha[i, j] <- alpha_pred[i] * alpha_prey[j].

I suspect that in many cases the assimilation efficiency will depend more on the prey species than on the predator species. in which case the user would only specify alpha_prey and leave alpha_pred at its default value. So the user would not have to specify more parameters but rather would be allowed to specify the parameters they actually care about.

There are already many extensions to mizer out there. @astaaudzi has an extension to use temperature-dependent rates and an extension to use multiple resource spectra. @camillanovaglio has a multi-fleet extension. @baldrech has an evolution extension. @reumj has an extension with size-dependent interaction matrix. I am sure there are many others.

It is therefore time that we developed an official mizer extension mechanism, so that extensions can be packaged up as R packages and will plug in flawlessly into mizer. So instead of adding more and more features to core mizer, there would be a growing number of extension packages and users would install only the packages relevant to them.

These extensions should be able to

• add more dynamic variables in addition to the currently existing n and n_pp slots. To avoid any restriction on the number and types of these variables, they will be passed around in a list.
• register call-back functions with project() to update any new variables inside the time loop.
• add more parameters to the MizerParams object. To avoid restrictions, these will be put in a list slot.
• register call-back functions with MizerParams() to set the new parameters.
• overwrite any of the rate functions called by getRates() to change how the rates are calculated.

We will flesh out this extension framework over the coming months as we work on converting the existing extensions into pluggable extensions.

We have two sources of documentation for the user.

1. There are the help pages that are automatically produced from the roxygen comments in the code, that the user can find on the "Reference" pages on the website (https://sizespectrum.org/mizer/dev/reference/index.html)
2. There are the articles on the website that are built from the .Rmd files in the vignettes subdirectory.

As we work on the code, it is relatively easy to remember to update the code documentation, but it is also easy to miss some sections in the vignettes that need updating. So any help with checking over the articles from time to time to update them would be much appreciated.

There is a vignette entitled "The general mizer size-spectrum model", that is currently an incomplete compilation of the model description that is scattered throughout the code documentation. For example its "Fishing" section currently says: "still to be written", even though the content that should go there is already contained in the help files of setFishing() and getFMort().

Maybe the "The general mizer size-spectrum model" vignette should be replaced by a help page that pulls together all the relevant information directly from the code documentation. This can be done with the roxygen @inheritSection tag, provided there are sections in the code documentation that are written in such a way as to be suitable for this purpose. Any help towards improving the code documentation to facilitate this would be much appreciated.

Generally in the code documentation there are two places where the model should be described:

1. in the functions for setting up a model, see for example the help for setParams(),
2. in the functions calculating the rates (you can find a list of these on the help page for getRates().

You will see that the help page for setParams() is much more complete than that for getRates(). That is because the former includes sections from the help pages of its constituent functions. Probably getRates() should do the same.

Currently the appropriate default value for gamma, the coefficient in the search volume, can only be determined when the predation kernel is log-normal. The historic reason for this is that for a lognormal kernel, the default value can be calculated analytically. For arbitrary feeding kernels it needs to be calculated numerically, which however is no problem. The code to do this should go into the non-exported function get_gamma_default().

By default mizer uses an allometric function for the proportion of the available energy that a mature individual invests into reproduction:

repro_prop(w) = (w / w_inf) ^ (m - n)

where n is the scaling exponent of the energy income rate. The exponent m is the scaling exponent of the rate at which energy is invested into reproduction. We currently choose the default m = 1, for no particular reason as far as I know. I suspect that different species will use different strategies, described by different values for m.

We could choose the default value for m so that the mizer growth curve and the von Bertalanffy growth curve agree also at some point between maturity size and asymptotic size (in the new version of mizer we already choose h so that the two growth curves agree at maturity size, and the above form for repro_prop guarantees that they agree at asymptotic size).

Or we could set the default for m so that the slope of the growth curve for mature individuals agrees with the slope of the von Bertalanffy curve at maturity size. I suspect this one @Kenhasteandersen might prefer, because in the old version of mizer he used this criterion to set the default for h.

Please comment.