diff --git a/NEWS.md b/NEWS.md index a3843411..1d58dff6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,16 @@ # mizer 3.1.0.9000 (development version) +- The "Getting started" guide now includes a self-contained "A worked example: + the Celtic Sea" section that takes a real ecosystem from raw species parameters + through building, finding the steady state, calibrating to observed biomasses + and growth, checking against observed yields, setting the resilience to + fishing, and projecting a fishing scenario whose sustainable-yield curve is + interpreted — the whole mizer workflow in one place (#450). + +- The first example in the "Getting started" guide no longer prints the + parameter-default notes from `newMultispeciesParams()`, which were alarming to + new users out of context (#450). + - The `p` argument of `setMetabolicRate()` is deprecated (#459). It never had any effect on a `MizerParams` object: such an object always has a `p` column already, and the argument was only ever used to fill in a missing one, so it diff --git a/vignettes/mizer.Rmd b/vignettes/mizer.Rmd index da50e8bc..f9495531 100644 --- a/vignettes/mizer.Rmd +++ b/vignettes/mizer.Rmd @@ -30,6 +30,11 @@ described in more detail in sections below. 4. [Exploring the results]. +Once you have seen these stages individually, the +[worked example](#a-worked-example-the-celtic-sea) near the end of this page +puts them together into a complete, self-contained study of a real ecosystem — +building, calibrating and projecting a Celtic Sea model. + If you run into any difficulties or have any questions or suggestions, let us know about it by posting about it on our [issue tracker](https://github.com/sizespectrum/mizer/issues/new). You can also twitter to @[mizer_model](https://twitter.com/mizer_model). We love to hear from @@ -145,11 +150,17 @@ For a more complicated multi-species model you need to provide a data frame with some species parameters. An example of a North Sea model is included with the package. Here we also use a species interaction matrix for the North Sea species. -```{r} +```{r message=FALSE, warning=FALSE} params <- newMultispeciesParams(NS_species_params, NS_interaction) ``` -The notes printed out by the function show us that mizer calculated -default values for many parameters that were not provided. +When you run this yourself, the function prints some notes explaining that, +because only a few species parameters were supplied, mizer has calculated +default values for many others (for example the maximum intake rate `h` and the +search volume `gamma`) using size-based theory. These notes are informational, +not errors: they simply flag the parameters you may want to refine later once +you calibrate the model to data. We have suppressed them here to keep the output +short. Calibrating such a model to observed biomasses and yields is demonstrated +in the [worked example](#a-worked-example-the-celtic-sea) below. @@ -187,9 +198,224 @@ getProportionOfLargeFish(sim, We can then use the full power of R to work with these results. The functionality provided by mizer to explore the simulation results is more -fully described in +fully described in [the section on exploring the simulation results.](https://sizespectrum.org/mizer/articles/exploring_the_simulation_results.html) +## A worked example: the Celtic Sea + +The four stages above introduced the mizer API one command at a time. This +section puts them together into a single, self-contained study of a real +ecosystem — the Celtic Sea — taking a model from raw species parameters all the +way to a fishing-scenario result: + +1. **Build** a multi-species model from species parameters, an interaction + matrix and fishing-gear definitions. +2. **Bring it to steady state** so the size spectra are self-consistent. +3. **Calibrate** it so that its steady state reproduces the *observed* biomasses + and growth rates. +4. **Check** it against an independent quantity — the observed fisheries yields. +5. **Set its resilience** to fishing. +6. **Project a fishing scenario** and **interpret an indicator** — here, the + total sustainable yield as a function of fishing effort. + +The model is the one built step-by-step, with much more explanation of every +choice, in the [Build part of the mizer +course](https://mizer.course.sizespectrum.org/build/). It is based on the Celtic +Sea model of @spence_parameterizing_2021, with species parameters and +observations drawn from FishBase and the ICES stock-assessment database. Here we +concentrate on seeing the whole workflow run in one place. + +```{r include=FALSE} +knitr::opts_chunk$set(message = FALSE, warning = FALSE) +``` + +### Getting the data + +The model needs three inputs: a species-parameter data frame, a species +interaction matrix, and a table of fishing-gear parameters. We also load a set +of observed commercial yields to check the model against later. These small +files live in the course repository; the code below downloads them into your +working directory the first time you run it. + +```{r} +base_url <- "https://github.com/gustavdelius/mizerCourse/raw/master/build/" +files <- c("celtic_species_params.rds", "celtic_gear_params.csv", + "celtic_interaction.csv", "celtic_yields.rds") +for (f in files) { + if (!file.exists(f)) download.file(paste0(base_url, f), destfile = f) +} + +celtic_species_params <- readRDS("celtic_species_params.rds") +celtic_gear_params <- read.csv("celtic_gear_params.csv") +celtic_interaction <- as.matrix(read.csv("celtic_interaction.csv", + row.names = 1)) +``` + +The species-parameter data frame has one row per species. Besides the required +maximum size `w_max` it carries a few life-history parameters and, crucially, +the **observed** average biomass of each species (`biomass_observed`, in grams +per square metre) that we will calibrate to. Some species have no observation +and are left as `NA`. + +```{r} +celtic_species_params[, c("species", "w_max", "biomass_observed")] +``` + +### Building the model + +`newMultispeciesParams()` turns these inputs into a `MizerParams` object, +filling in with size-based defaults every parameter we have not supplied. We +set the allometric exponents `n` and `p` both to 3/4 and switch the commercial +gear on at unit effort. + +```{r} +cel <- newMultispeciesParams( + species_params = celtic_species_params, + gear_params = celtic_gear_params, + interaction = celtic_interaction, + n = 3/4, p = 3/4, + initial_effort = 1 +) +``` + +As before, mizer prints notes about the defaults it has chosen; we have +suppressed them here to keep the output short. + +### Finding the steady state + +A freshly built model has only a rough spectrum. `steady()` runs the +size-spectrum dynamics, holding reproduction and the resource fixed, until the +community settles onto a steady state. + +```{r} +cel <- steady(cel) +plotSpectra(cel, power = 2) +``` + +Each curve is now a smooth species spectrum, and together they line up along the +background resource — the hallmark of a self-consistent size-spectrum model. + +### Calibrating to observed biomass + +The steady state is self-consistent, but its species abundances are still +arbitrary. Calibration rescales them to match observation. We do this in stages, +**re-running `steady()` after every adjustment** because each change pushes the +community off its steady state: + +- `calibrateBiomass()` sets the overall abundance scale (the resource level + `kappa`) so that the *total* community biomass matches the sum of the observed + biomasses. +- `matchBiomasses()` then rescales each species individually to its own + `biomass_observed`. +- `matchGrowth()` adjusts intake, search volume and metabolism so the fish reach + maturity at the right age. + +`matchBiomasses()` and `matchGrowth()` pull on different parameters, so we +alternate them, re-converging each time, until both are satisfied. + +```{r} +cel <- calibrateBiomass(cel) +for (i in 1:4) { + cel <- matchBiomasses(cel) + cel <- matchGrowth(cel) + cel <- steady(cel) +} +``` + +`plotBiomassObservedVsModel()` shows how well the calibrated steady state +reproduces the observations: points on the diagonal are a perfect match. + +```{r} +plotBiomassObservedVsModel(cel) +``` + +Most species sit close to the 1:1 line. A few remain off — the worst by nearly +an order of magnitude — and these would need more iteration, or a look at their +predation-kernel parameters, before you trusted the model for them. Diagnosing +and fixing such residuals is exactly what the [refinement +tutorial](https://mizer.course.sizespectrum.org/build/refine.html) of the course +is about. + +### Checking against observed yield + +Biomass was our calibration target. The fisheries **yield** is an independent +quantity: we did not tune to it, so comparing modelled and observed yields is a +genuine test of the model. We attach the observed yields and plot them against +the model. + +```{r} +celtic_yields <- readRDS("celtic_yields.rds") +species_params(cel)$yield_observed <- as.numeric(celtic_yields) +plotYieldObservedVsModel(cel) +``` + +The modelled yields come out at the right overall level and mostly within a +factor of a few of the observations — reassuring, given that yield was not a +calibration target. Where a species' yield is badly off, the usual culprit is +its gear selectivity or catchability, which you would adjust next (see the +[landings tutorial](https://mizer.course.sizespectrum.org/build/landings.html)). + +### Setting the resilience to fishing + +Calibrating the steady state fixes *where* the community sits, but not *how +sensitively* it responds when we change fishing. That sensitivity is governed by +the strength of density dependence in reproduction. `setBevertonHolt()` sets it +through the `reproduction_level` — the fraction of the maximum recruitment that +is realised at the steady state. A value of 0.5 gives moderate density +dependence, a common default in the absence of stock-specific information. + +```{r} +cel <- setBevertonHolt(cel, reproduction_level = 0.5) +``` + +This changes the reproduction parameters without moving the steady state, so the +biomass and yield fits above are preserved. + +### Projecting a fishing scenario + +We now have a calibrated, resilient model. Let us use it to ask a management +question: **is the fishery being fished at the effort that maximises the total +sustainable yield?** + +We sweep fishing effort from zero to twice its current value. For each level we +project to the new steady state and record the total community yield — the yield +that could be taken indefinitely at that effort. + +```{r} +efforts <- c(0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2) +sustainable_yield <- sapply(efforts, function(e) { + p <- projectToSteady(cel, effort = e, t_max = 100, + return_sim = FALSE, progress_bar = FALSE) + sum(getYield(p)) +}) + +plot(efforts, sustainable_yield, type = "b", pch = 19, + xlab = "Fishing effort (relative to current)", + ylab = "Total sustainable yield (g/m²/year)") +abline(v = 1, lty = 2) +``` + +### Interpreting the result + +The curve is the community-level analogue of a single-stock yield curve. It +rises, peaks, and then falls as heavier fishing starts to erode the stocks +faster than they can rebuild. The dashed line marks the **current** effort. + +The peak — the effort giving maximum sustainable yield for the community as a +whole — lies a little *below* the current effort. In other words, this fishery +is already working slightly harder than the level that would maximise its total +long-term catch: pushing effort higher buys no extra yield and depletes the +community, while easing off a little would *increase* the sustainable catch. That +is a textbook signature of a lightly overfished community, recovered here purely +from the model we calibrated to biomass and growth. + +This is only the beginning of what the calibrated model can tell you. From here +you could look at individual species' yield curves to find each stock's own +$F_\text{MSY}$ (`plotYieldVsF()`), change the gear selectivity to protect +juveniles, or explore how the community size spectrum flattens under fishing. +The [Use part of the course](https://mizer.course.sizespectrum.org/use/) +develops several such scenarios. + ## Size-spectrum models Size spectrum models have emerged as a conceptually simple way to model a large community of individuals which grow and change trophic level during life. There diff --git a/vignettes/reflib.bib b/vignettes/reflib.bib index a01b59bb..0356242c 100644 --- a/vignettes/reflib.bib +++ b/vignettes/reflib.bib @@ -1251,3 +1251,15 @@ @article{jacobsen_consequences_2014 file = {Snapshot:/home/finlay/.mozilla/firefox/99s2uy74.default/zotero/storage/BUSIVZU6/20132701.html:text/html} } + + +@article{spence_parameterizing_2021, + title = {Quantifying uncertainty and dynamical changes in multi-species fishing mortality rates, catches and biomass by combining state-space and size-based multi-species models}, + volume = {22}, + number = {4}, + journal = {Fish and Fisheries}, + author = {Spence, Michael A. and Thorpe, Robert B. and Blackwell, Paul G. and Scott, Finlay and Southwell, Richard and Blanchard, Julia L.}, + year = {2021}, + pages = {667--681}, + doi = {10.1111/faf.12543} +}