Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merging corrections from dev branch #16

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions build/create-first-model.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ curonian_interaction <- read.csv("curonian_interaction.csv", row.names = 1)

We use `rownames = 1` to let `read.csv` know that the first column in the spreadsheet, which contains the predator names, should be used as the names of the rows of the interaction matrix.

When you repeat this work in your worksheet, you can check that the data was read in correctly by clicking on `curonian_species_params`, `curonian_species_params` and `curonian_interaction` in the "Environment" tab. That will open the data frames in your editor window for you to inspect.
When you repeat this work in your worksheet, you can check that the data was read in correctly by clicking on `curonian_species_params`, `curonian_gear_params` and `curonian_interaction` in the "Environment" tab. That will open the data frames in your editor window for you to inspect.

We will now set up a multi-species mizer model using the function `newMultispeciesParams()`. Besides the species parameters, the gear parameters and the interaction matrix, the other information that flows into a multi-species model are the resource parameters, the allometric exponents `n` and `p` and the fishing effort.

Expand Down Expand Up @@ -84,7 +84,7 @@ We will now project to the steady state, which will finally give us realistic sp
cur_model <- steady(cur_model)
```

We can ignore the warning from the `setBevertonHolt()` function about unrealistic reproductive efficiencies. Those warnings are an artefact of how the reproduction level is set by default by the `matchBiomasses()` function. We could fix those defaults, but we are not yet concerned with the reproduction dynamics so we don't have to do that and just ignore the warnings.
We can ignore the warning from the `setBevertonHolt()` function about unrealistic reproductive efficiencies. Those warnings are an artefact of how the reproduction level is set by default. We could fix those defaults, but we are not yet concerned with the reproduction dynamics so we don't have to do that and just ignore the warnings.

::: {.callout-important collapse="false"}
Remember, the steady() function does not yet calibrate your model to any observed biomasses or catches. It simply runs the model to find a state where species spectra stop changing. If you start a simulation from initial biomasses that the set of current species and resource parameters cannot support, species biomasses will change drastically at the start. Perhaps there are too many predators at the start, or too many fish for the resource to support it. After a while these fluctuations will settle down (hopefully) and species biomasses will be stable for a given set of parameters. Of course, since we are keeping reproduction and resource levels constant, this is not yet a final model.
Expand Down
2 changes: 1 addition & 1 deletion build/set-multiple-resources.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ title: "Setting up multiple resources"

Most mizer models simulate multiple size structured species and one background resource. But what if you want to account for different types of 'resources' with different dynamics? The most obvious case is the separation of pelagic and benthic resource. Benthic resource may not be important for offshore or pelagic ecosystems, but if you want to model coastal or freshwater ecosystems this may become essential.

Adding benthic resource into multi-species and size based models is not new. Some researchers incorporate benthic food sources by explicitly modelling key benthic species as size-structured groups. This is fine if you have data on their biomases for calibration and if you are interested in their dynamics, but modelling many benthic species separately might be overwhelming. Instead, just like with plankton species, we might want to model all benthic organisms as a separate resource spectrum. That is what the mizerMR extension does.
Adding benthic resource into multi-species and size based models is not new. Some researchers incorporate benthic food sources by explicitly modelling key benthic species as size-structured groups. This is fine if you have data on their biomasses for calibration and if you are interested in their dynamics, but modelling many benthic species separately might be overwhelming. Instead, just like with plankton species, we might want to model all benthic organisms as a separate resource spectrum. That is what the mizerMR extension does.

The origin of mizerMR is a [model for Tasmanian coastal rocky reef ecosystem](https://doi.org/10.1101/2022.06.20.496925), where most large fish are not predators (like in a classical size based marine ecosystem), but feed on benthic invertebrates. In such reefs, benthic production pathway is likely to be a major source of energy. The mizerMR framework has also been used in a [Baltic Sea mizer model](https://onlinelibrary.wiley.com/doi/full/10.1111/gcb.16341) and a [south-east Australian mizer model](https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/1365-2664.14086). There are other ways to incorporate benthic resources, and you can find some examples [here](https://onlinelibrary.wiley.com/doi/epdf/10.1111/geb.13348), [here](https://link.springer.com/article/10.1007/s12080-010-0078-9) and [here](https://onlinelibrary.wiley.com/doi/abs/10.1111/oik.05630), but none of these approaches are currently implemented in the main mizer code.

Expand Down
10 changes: 6 additions & 4 deletions understand/dynamics-of-spectra.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ library(tidyverse)

In the previous tutorial, in the section on [trophic cascades](species-interactions.qmd#trophic-cascades), we already simulated the size-spectrum dynamics to find the new steady state. But we only looked at the final outcome once the dynamics had settled down to the new steady state. We reproduce the code here:

```{r}
```{r message=FALSE}
# Create trait-based model
mp <- newTraitParams() |>
# run to steady state with constant reproduction rate
Expand All @@ -60,7 +60,8 @@ species_params(mp_lessRes)$interaction_resource[8:11] <- 0.8
mp_lessRes_steady <- projectToSteady(mp_lessRes)

# We compare the steady states
plotSpectra2(mp_lessRes_steady, "less resource", mp, "original",
plotSpectra2(mp_lessRes_steady, name1 = "less resource",
mp, name2 = "original",
total = TRUE, power = 2,
ylim = c(1e-8, NA), wlim = c(1e-3, NA))
```
Expand Down Expand Up @@ -114,7 +115,7 @@ The above simulation was run with constant abundance in the smallest size class

To see the effect we run the same code as above after deleting the two lines that turned off the reproduction dynamics. We also specify with `t_save = 2` that we want to save the spectrum only every second year, which speeds up the display of the animation.

```{r}
```{r message=FALSE}
# Create trait-based model
mp <- newTraitParams() |>
# run to steady state with constant reproduction rate
Expand Down Expand Up @@ -241,7 +242,8 @@ Changing the reproduction level has no effect on the steady state, because that
```{r}
mp9 <- projectToSteady(mp9)

plotSpectra2(mp, "reproduction_level = 0.25", mp9, "reproduction_level = 0.9",
plotSpectra2(mp, name1 = "reproduction_level = 0.25",
mp9, name2 = "reproduction_level = 0.9",
total = TRUE, power = 2, ylim = c(1e-8, NA), wlim = c(1e-3, NA))
```

Expand Down
3 changes: 2 additions & 1 deletion understand/predation-growth-and-mortality.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,8 @@ We can now see how fishing affects the adult size spectrum.
# We find the new steady state with fishing mortality imposed
params_fishing <- singleSpeciesSteady(params_fishing)

plotSpectra2(params, "No Fishing", params_fishing, "Knife-edge",
plotSpectra2(params, name1 = "No Fishing",
params_fishing, name2 = "Knife-edge",
power = 2, wlim = c(10, NA))

```
Expand Down
4 changes: 2 additions & 2 deletions understand/single-species-spectra.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -447,8 +447,8 @@ params_changed_maturity <- singleSpeciesSteady(params_changed_maturity)
Then we can use the function `plotSpectra2()` to plot the old size spectrum and the new size spectrum on the same graph.

```{r}
plotSpectra2(params, "Early maturity",
params_changed_maturity, "Late maturity",
plotSpectra2(params, name1 = "Early maturity",
params_changed_maturity, name2 = "Late maturity",
power = 2, resource = FALSE, wlim = c(10, NA))
```

Expand Down
8 changes: 5 additions & 3 deletions understand/species-interactions.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ In this first week we aim for understanding, not realism. So in this tutorial we

We use the `newTraitParams()` function to create our idealised trait-based multi-species model. The function has many parameters, but we will just keep the defaults. Unlike the `newSingleSpeciesParams()` function, the `newTraitParams()` function does not set the initial spectra to their steady state values. We thus need to run the result through the `steady()` function (we'll discuss that function more next week). We assign the resulting MizerParams object to the variable `mp`.

```{r}
```{r message=FALSE}
mp <- newTraitParams() |> steady()
```

Expand Down Expand Up @@ -231,7 +231,8 @@ mp_lessRes_steady <- projectToSteady(mp_lessRes)
```

```{r}
plotSpectra2(mp_lessRes_steady, "less resource", mp, "original",
plotSpectra2(mp_lessRes_steady, name1 = "less resource",
mp, name2 = "original",
total = TRUE, power = 2, ylim = c(1e-8, NA), wlim = c(1e-3, NA))
```

Expand All @@ -240,7 +241,8 @@ There is much to see in this graph. We can see how the reduction in the abundanc
Perhaps you would prefer to plot the above graph with `power = 1`, which will show the biomass density in weight instead of the Sheldon spectrum. Different people find different options more intuitive.

```{r}
plotSpectra2(mp_lessRes_steady, "less resource", mp, "original",
plotSpectra2(mp_lessRes_steady, name1 = "less resource",
mp, name2 = "original",
total = TRUE, power = 1, ylim = c(1e-8, NA), wlim = c(1e-3, NA))
```

Expand Down
11 changes: 6 additions & 5 deletions use/change-resources.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -76,14 +76,14 @@ c_R <- (mu_bar + r_R) / r_R * N_bar
mu_R <- seq(0, 2, by = 0.05)
N_R <- c_R * r_R / (mu_R + r_R)
df <- melt(data.frame(N_R, mu_R, c_R), id.vars = "mu_R")
df$dd <- "Low"
df$dd <- "High"

r_R <- 2
c_R <- (mu_bar + r_R) / r_R * N_bar
mu_R <- seq(0, 2, by = 0.05)
N_R <- c_R * r_R / (mu_R + r_R)
df2 <- melt(data.frame(N_R, mu_R, c_R), id.vars = "mu_R")
df2$dd <- "High"
df2$dd <- "Low"

df <- rbind(df, df2)
df <- df[df$value < 2.1, ]
Expand All @@ -97,7 +97,7 @@ ggplot(rbind(df, df2)) +
breaks = N_bar, labels = c("N_R")) +
scale_x_continuous(name = "Resource mortality rate [1/year]",
breaks = mu_bar, labels = c("mu_R")) +
scale_colour_manual(values = c("black", "blue")) +
scale_colour_manual(values = c("blue","black")) +
scale_linetype_manual(values = c("solid", "dotted"))
```

Expand Down Expand Up @@ -210,7 +210,8 @@ cm_less_be <- setResourceSemichemostat(cm, rp_less_be)
Let's check that resource abundances really changed by plotting the spectra

```{r}
plotSpectra2(cm_more_be, "More", cm_less_be, "Less", power = 2)
plotSpectra2(cm_more_be, name1 = "More",
cm_less_be, name2 = "Less", power = 2)
```

As we discussed before, while the `setResourceSemichemostat()` has set the resource carrying capacity so that initially the new resource abundance is steady, the fish abundances will change over time in response to the changed resource abundance and as a result also the resource abundance will start to change. Let us investigate that by projecting for 50 years with the new level and compare to the unchanged conditions.
Expand Down Expand Up @@ -264,7 +265,7 @@ cm_more_be <- setResourceSemichemostat(cm, rp_more_be)
cm_less_be <- setResourceSemichemostat(cm, rp_less_be)

# Check on the spectra plot that the changes have taken place
plotSpectra2(cm_more_be, "More", cm_less_be, "Less", power = 2)
plotSpectra2(cm_more_be, name1 = "More", cm_less_be, name2 = "Less", power = 2)
```

Now we will project for 50 years with the new level and see how the biomasses are changed
Expand Down
2 changes: 1 addition & 1 deletion use/fishing-scenarios.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ We can see that, surprisingly, after doubling fishing effort on all species, the
To understand how that is possible, we have to look at how the spectra have changed. After all mizer is all about changes at size. Looking at overall biomasses does not tell use about size dynamics. We can compare the spectra using the `plotSpectra2()` function:

```{r}
plotSpectra2(sim_double_effort, 'double', cm, 'same', power = 2)
plotSpectra2(sim_double_effort, name1 = 'double', cm, name2 = 'same', power = 2)
```

Again, the changes are more visible if we plot relative change with respect to the initial spectra:
Expand Down
5 changes: 2 additions & 3 deletions use/further-scenarios.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ The abundance of some (but not all) of the smaller to medium sizes of prey are a
Mizer can also be used to carry out projections with changes in fishing effort. We will start by reading in a time series of catch and effort. We will read in the effort data for this fishery. Note: these are not the correct effort values, they are only for illustration in this example. Then we plot effort and catch through time.

```{r, code_folding=TRUE}
dat <- readRDS("toothfish/longline.rds")
dat <- readRDS("longline.rds")
plot_dat <- melt(dat, "Year")

ggplot(plot_dat, aes(x = Year, y = value)) +
Expand Down Expand Up @@ -281,8 +281,7 @@ Let's take a look at the relative exploitation status of the stocks using the pr

```{r}
# plot change in biomass under each scenario relative to unfished values
# get saved values from steady state without fishing that we generated earlier
sim0<-readRDS("sim0.rds")
# using the sim0 values from steady state without fishing that we generated earlier
# get the unfished biomasses
B_unfished <- getBiomass(sim0)[1, ]
#scen 1
Expand Down
File renamed without changes.
2 changes: 1 addition & 1 deletion use/tune-resilience.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ To measure the resilience of our species to fishing, we will change the fishing
plotYieldCurve(cm_generic_gear, species = "smallfish", F_max = 2)
```

We refer to the yield achieved in the steady state as the *sustainable yield* (SY) because because it can be sustained indefinitely. The sustainable yield for this species has a maximum (MSY) at a fishing mortality of about 0.8. We refer to this value as $F_{MSY}$. The small fish have a high resilience to fishing. This is close enough to expectations, so we judge that it is not necessar to make changes to the resource level of "smallfish".
We refer to the yield achieved in the steady state as the *sustainable yield* (SY) because because it can be sustained indefinitely. The sustainable yield for this species has a maximum (MSY) at a fishing mortality of about 0.8. We refer to this value as $F_{MSY}$. The small fish have a high resilience to fishing. This is close enough to expectations, so we judge that it is not necessary to make changes to the resource level of "smallfish".

Next we look at ruffe. This is classified as a medium resilience species, so we would like the $F_{MSY}$ to be between 0.1 and 0.5.

Expand Down