library(mizer) gear_names <- c("Industrial", "Pelagic", "Beam", "Otter") times <- seq(from = 1, to = 10, by = 1) effort_array <- array(NA, dim = c(length(times), length(gear_names)), dimnames = list(time = times, gear = gear_names)) effort_array[, "Industrial"] <- 0.5 effort_array[, "Pelagic"] <- seq(from = 1, to = 2, length = length(times)) effort_array[, "Beam"] <- seq(from = 1, to = 0, length = length(times)) effort_array[, "Otter"] <- seq(from = 1, to = 0.5, length = length(times))
In the sections on the multispecies model and on
running a simulation we saw how to set up a model
and project it forward through time under our desired fishing scenario. The
result of running a projection is an object of class MizerSim
. What do we then
do? How can we explore the results of the simulation? In this section we
introduce a range of summaries, plots and indicators that can be easily produced
using functions included in mizer.
We will use the following MizerSim
object for these examples, where the
effort array is the one we created in the
previous section on running a simulation:
sim <- project(NS_params, effort = effort_array, dt = 0.1, t_save = 1)
The projected species abundances at size through time can be obtained with
N(sim)
. This returns a three-dimensional array (time x species x size). Consequently,
this array can get very big so inspecting it can be difficult. In the example we
have just run, the time dimension of n
has r dim(N(sim))[1]
rows
(one for the initial population and then one for each of the saved time steps).
There are also 12 species each with 100 sizes. We can check this by running the
dim()
function and looking at the dimensions of the n
array:
dim(N(sim))
To pull out the abundances of a particular species through time at size you can subset the array. For example to look at Cod through time you can use:
N(sim)[, "Cod", ]
This returns a two-dimensional array: time x size, containing the cod
abundances. The time dimension depends on the value of the argument t_save
when project()
was run. You can see that even though we specified dt
to be
0.1 when we called project()
, the t_save = 1
argument has meant that the
output is only saved every year.
Often we are particularly interested in the results at the final time-step. These we can access with
finalN(sim)
which is a two dimensional array (species x size).
The projected resource abundances can be accesses similarly with
NResource(sim)
This returns a two-dimensional array (time x size). And if we are only interested in the final time step
finalNResource(sim)
returns a vector with one entry for each size class.
As well as the summary()
methods that are available for both MizerParams
and
MizerSim
objects, there are other useful summary functions to pull information
out of a MizerSim
object. A description of the different summary functions
available is given in
the summary functions help page.
All of these functions have help files to explain how they are used. (It is also
possible to use most of these functions with a MizerParams
object if you also
supply the population abundance as an argument. This can be useful for exploring
how changes in parameter value or abundance can affect summary statistics and
indicators. We won't explore this here but you can see their help files for more
details.)
The functions getBiomass()
and getN()
have additional arguments that allow the
user to set the size range over which to calculate the summary statistic. This
is done by passing in a combination of the arguments min_l
, min_w
, max_l
and max_w
for the minimum and maximum length or weight.
If min_l
is specified there is no need to specify min_w
and so on. However,
if a length is specified (minimum or maximum) then it is necessary for the
species parameter data.frame (see
the species parameters section)
to include
the parameters a
and b
for length-weight conversion. It is possible to mix
length and weight constraints, e.g. by supplying a minimum weight and a maximum
length. The default values are the minimum and maximum weights of the spectrum,
i.e. the full range of the size spectrum is used.
Here we show a simple demonstration of using a summary function using the sim
object we created earlier. Here, we use getSSB()
to calculate the SSB of each
species through time (note the use of the head()
function to only display the
first few rows).
ssb <- getSSB(sim) dim(ssb) head(ssb)
As mentioned above, we can specify the size range for the getsummaryBiomass()
and
getN()
functions. For example, here we calculate the total biomass of each
species but only include individuals that are larger than 10 g and smaller than
1000 g.
biomass <- getBiomass(sim, min_w = 10, max_w = 1000) head(biomass)
Functions are available to calculate a range of indicators from a MizerSim
object after a projection. A description of the different indicator functions
available is given in the indicator functions help page..
You can read the help pages for each
of the functions for full instructions on how to use them, along with examples.
With all of the functions in the table it is possible to specify the size range of
the community to be used in the calculation (e.g. to exclude very small or very
large individuals) so that the calculated metrics can be compared to empirical
data. This is used in the same way that we saw with the function getBiomass()
in
the section on summary functions for MizerSim objects..
It is also
possible to specify which species to include in the calculation. See the help
files for more details.
For these examples we use the sim
object we created earlier.
The slope of the community can be calculated using the getCommunitySlope()
function. Initially we include all species and all sizes in the calculation (only
the first five rows are shown):
slope <- getCommunitySlope(sim) head(slope)
This gives the slope, intercept and $R^2$ value through time (see the help file
for getCommunitySlope
for more details).
We can include only the species we want with the species
argument. Below we
only include demersal species. We also restrict the size range of the community
that is used in the calculation to between 10 g and 5 kg. The species
argument is a
character vector of the names of the species that we want to include in the
calculation.
dem_species <- c("Dab", "Whiting", "Sole", "Gurnard", "Plaice", "Haddock", "Cod", "Saithe") slope <- getCommunitySlope(sim, min_w = 10, max_w = 5000, species = dem_species) head(slope)
R is very powerful when it comes to exploring data through plots. Two useful
packages for plotting are ggplot2
and plotly
. These use data.frames for input data
whereas many of the mizer functions return arrays or
matrices. Fortunately it is straightforward to turn arrays and matrices into
data.frames using the melt()
function from the reshape2
package that mizer
makes available to you.
Although mizer
does include some dedicated plots, it is definitely worth your
time getting to grips with these other plotting packages. This
will make it possible for you to make your own plots. We provide some details
in the section on using ggplot2 and plotly with mizer.
Included in mizer
are several dedicated plots that use MizerSim
objects as
inputs (see the plots help page.). As
well as displaying the plots, these functions all return objects of type ggplot
from the ggplot2
package, meaning that they can be further modified by the user
(e.g. by changing the plotting theme). See the help page of the individual plot
functions for more details. The generic plot()
method has also been overloaded
for MizerSim
objects. This produces several plots in the same window to
provide a snapshot of the results of the simulation.
Some of the plots plot values by size (for example plotFeedingLevel()
and
plotSpectra()
). For these plots, the default is to use the data at the final
time step of the projection. With these plotting functions, it is also possible
to specify a different time, or a time range to average the values over before
plotting.
Using the plotting functions is straightforward. For example, to plot the total
biomass of each species against time you use the plotBiomass()
function:
plotBiomass(sim)
As mentioned above, some of the plot functions plot values against size at a point in time (or averaged over a time period). For these plots it is possible to specify the time step to plot, or the time period to average the values over. The default is to use the final time step. Here we plot the abundance spectra (biomass), averaged over time = 5 to 10:
plotSpectra(sim, time_range = 5:10)
As mentioned above, and as we have seen several times in this guide, the generic
plot()
method has also been overloaded. This produces 5 plots in the same
window (plotFeedingLevel()
, plotBiomass()
, plotPredMort()
, plotFMort()
and plotSpectra()
). It is possible to pass in the same arguments that these
individual plots use, e.g. arguments to change the time period over which the
data is averaged.
plot(sim)
The next section describes how to use what we have learned to model the North Sea.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.