inst/skills/analyse-and-plot/SKILL.md

In mizer: Dynamic Multi-Species Size Spectrum Modelling

name: analyse-and-plot description: >- Analyse and visualise the results of a mizer simulation or the state of a MizerParams object. Use whenever the user wants to extract, summarise or plot size spectra, biomass, numbers, yield, SSB, feeding level, mortality, diet, trophic level, community indicators, growth curves or the resource — including comparing two models and animating spectra through time. Also covers choosing what a density plot shows (biomass, per_log_size, size_axis), the ArraySpeciesBySize/ArrayTimeBySpecies wrappers the getters return, and writing a custom indicator with sizeIntegral(). Prefer these functions over hand-written array wrangling or custom ggplot code. If a plotting argument that used to work now errors (power=), see the upgrade-mizer-code skill.

Analysing and plotting mizer results

Mizer ships a large family of extraction, summary, and plotting functions. Always prefer these over hand-written array wrangling or custom ggplot code — they handle size-range integration, species colours/linetypes, and units for you.

Most functions accept either a MizerSim object (returning a time series) or a MizerParams object (returning a single value from the initial state). So getBiomass(sim) gives biomass over time, getBiomass(params) gives biomass now.

To get the single value at one time step of a simulation, extract a MizerParams snapshot with finalParams(sim) (last step), initialParams(sim) (first step), or getParams(sim, time_range = ...) (averaged over a range) and pass that in:

getMeanMaxWeight(finalParams(sim))             # value at the last time step
getSSB(getParams(sim, time_range = 1990:2000)) # averaged over a period

If you need a function you don't see here, grep for "plot" or the specific name in the bundled API index (path at the end of MIZER-AGENTS.md) before writing custom code — don't read the whole file. The index gives you the name; read the help page for the arguments.

Accessing simulation arrays

These extract raw arrays from a MizerSim object.

| Function | Returns | Dimensions | |---|---|---| | N(sim) | species abundance density | time × species × size | | NResource(sim) | resource abundance density | time × size | | finalN(sim) | species abundance at last time | species × size | | finalNResource(sim) | resource abundance at last time | size | | getEffort(sim) | fishing effort | time × gear | | getTimes(sim) | saved time steps | time |

N(sim)[, , 1]              # time × species in smallest size class
N(sim)["2010", "Cod", ]    # size vector for Cod in year 2010
finalN(sim)["Cod", ]       # size vector for Cod at the final time step

Summary functions

These functions compute derived quantities from abundances. All accept MizerSim or MizerParams. The result is a classed array that can be plotted directly with plot() — see below.

| Function | Returns | Dimensions | |---|---|---| | getBiomass(sim, min_w, max_w) | total biomass | time × species | | getSSB(sim) | spawning stock biomass | time × species | | getN(sim, min_w, max_w) | total abundance | time × species | | getYield(sim) | total yield across gears | time × species | | getYieldGear(sim) | yield by gear | time × gear × species | | getFeedingLevel(sim) | feeding level at size | time × species × size | | getPredMort(sim) | predation mortality at size | time × species × size | | getFMort(sim) | fishing mortality at size | time × species × size | | getFMortGear(sim) | fishing mortality by gear | time × gear × species × size | | getDiet(params) | diet resolved by prey at size | predator × size × prey | | getTrophicLevel(params) | trophic level at size | species × size | | getTrophicLevelBySpecies(params) | mean trophic level per species | species |

Size range: getBiomass() and getN() accept min_w, max_w, min_l, max_l to restrict the calculation to a size range.

getSSB(sim)                              # SSB of all species over time
getBiomass(sim, min_w = 10, max_w = 1e4) # biomass of 10g–10kg fish
getYield(sim)["2010", ]                  # yield in year 2010

Indicator functions

These compute community-level indicators. All accept MizerSim (time series) or MizerParams (single value from the initial state). See ?indicator_functions.

| Function | Key arguments | Returns | |---|---|---| | getProportionOfLargeFish(sim) | threshold_w = 100, biomass_proportion | proportion of large fish through time | | getMeanWeight(sim) | min_w, max_w, species | mean community weight through time | | getMeanLength(sim) | min_w, max_w, species | mean community length through time | | getMeanMaxWeight(sim) | measure = "both"/"numbers"/"biomass" | mean asymptotic weight through time | | getCommunitySlope(sim) | min_w, max_w, species | slope, intercept, R² through time |

lfi <- getProportionOfLargeFish(sim, min_w = 10, max_w = 5000, threshold_w = 500)
slope <- getCommunitySlope(sim, min_w = 10, max_w = 5000)

Writing your own indicator

First check that a built-in does not already cover it: most custom indicators turn out to be getBiomass()/getN() over a size range, or one of the four above with different arguments. If none fits, an indicator is an integral over the size spectrum, $\int N_i(w)\, K_i(w)\, dw$, where $K_i(w)$ is a weighting factor (supplied to the weighting argument of sizeIntegral()). sizeIntegral() does that integral for you:

# Abundance between 10g and 5kg (default weighting factor weighting = 1)
sizeIntegral(params, min_w = 10, max_w = 5000)

# Biomass between 10g and 5kg (weighting factor is body weight)
sizeIntegral(params, weighting = w(params), min_w = 10, max_w = 5000)

# Biomass through time, wrapped ready to plot
sizeIntegral(sim, weighting = w(params), value_name = "Biomass", units = "g")

Give it the object and the weighting factor $K$ (e.g. body weight w(params) for biomass, or 1 for numbers); it selects the size range, uses the quadrature scheme the model is actually on and wraps the result in the appropriate mizer array class. Doing the sum by hand instead means getting all of that right yourself, silently and only for some users. Three things are worth knowing:

• The size range is an argument, not a subsetting step. min_w/max_w or min_l/max_l are passed through to get_size_range_array(), which does the length-weight conversion per species and accepts either a single number or one value per species. Never subset the size grid by hand.
• Pass the whole product to weighting. If $K$ is a product of several size-dependent terms, build the combined product first — SSB uses sweep(params@maturity, 2, params@w, "*") (maturity $\times$ body weight), yield uses fishing mortality $\times$ body weight. Bin-averaging happens inside on the weighting array as a whole, and the average of a product is not the product of the averages. Do not include params@dw: sizeIntegral() handles bin widths and bin-averaging automatically.
• Extra dimensions in weighting are kept. A gear × species × size weighting array (like getFMortGear()) gives a gear × species result; a weighting array whose first dimension is named "time" is lined up with the times of the simulation rather than multiplied out against them.

The result is already an ArrayTimeBySpecies when it is one, so you inherit the whole toolkit described below — plot(), plot2(), plotRelative(), addPlot() — for free. For a quantity that keeps the size dimension, and so is not an integral over sizes, wrap it yourself:

ArraySpeciesBySize(my_size_resolved, value_name = "My index", params = params,
                   representation = "average")

Use representation = "average" for a quantity that is a bin average (anything integrated over a bin) and "point" for one sampled at the bin boundary, such as a growth rate; the tag drives the half-bin plotting shift.

If your indicator decomposes the encounter rate — a diet or trophic-level style quantity — see the note on encounter_kernel() in the extend-mizer skill before pairing pred_kernel() with getEncounter().

Plotting mizer arrays

The arrays returned by the summary and rate functions carry a mizer array class and have their own plot() method, so you can visualise any quantity without a dedicated plot function or custom ggplot code. They also carry a value_name, type, units and their params, and have print(), summary() and as.data.frame() methods.

| Class | Typical source | plot() shows | |---|---|---| | ArrayTimeBySpecies | getBiomass(sim), getSSB(sim), getYield(sim), getN(sim) | value vs time, one line per species | | ArraySpeciesBySize | getFeedingLevel(params), getPredMort(params), getEncounter(params) | value vs size, one line per species | | ArrayTimeBySpeciesBySize | getFMort(sim), getPredMort(sim) | one time slice vs size (set with time) | | ArrayResourceBySize | NResource(params), finalNResource(sim), getResourceMort(params), resource_rate(params), resource_capacity(params), resource_level(params) | resource quantity vs size | | ArrayTimeByResourceBySize | NResource(sim) | one time slice vs size (set with time) |

plot(getBiomass(sim))          # value vs time, one line per species
plot(getFeedingLevel(params))  # value vs size, one line per species
plot(getResourceMort(params))  # plankton resource mortality vs size

The array plots come with a small toolkit for combining and comparing them. Every one of these has a method for every array class in the table above:

| Function | What it does | |---|---| | addPlot() | adds a compatible array as extra lines on an existing plot | | plot2() | compares two compatible arrays (colour = species, linetype = which object) | | plotRelative() | shows the relative difference 2 (y - x) / (x + y) between two compatible arrays | | plotHover() | turns any of these ggplots into a hover-enabled plotly plot |

# Add another compatible array as extra lines on an existing plot
p <- plot(getBiomass(sim), species = "Cod")
addPlot(p, getBiomass(sim), species = "Herring", linetype = "dashed")

# Compare two compatible arrays
plot2(getFMort(params), getFMort(params2), "Before", "After")
plotRelative(getEGrowth(params), getEGrowth(params2))  # relative difference

plotHover(getBiomass(sim))     # interactive (hover) version of any array plot

plot2() and plotRelative() prepare each of their two arrays separately, using the model attached to that array. This matters when the two models differ in the length-weight relationship w = a l^b: with size_axis = "l" each spectrum is then drawn at its own lengths, and a density is rescaled by its own Jacobian. The two length grids no longer coincide when that happens, so plotRelative() interpolates both series (linearly in the logarithm of size) onto the union of their coordinates, restricted to the range both cover. On a weight axis, and whenever the two models agree, the grids coincide and nothing is approximated.

The two arrays must hold the same kind of value: comparing a "density" with a "value" is an error, because the type decides both the Jacobian and the y-axis scaling and there is no pair of axes that carries both. A differing value_name or differing units only warn.

Common arguments

Most analysis and plotting functions — including plot() on an array and the dedicated plot…() functions below — share these optional arguments:

| Argument | Effect | |---|---| | species | character vector — restrict to a subset of species | | time_range | numeric vector — average over this time period (plots against size) | | tlim | numeric vector c(min, max) — restrict the time axis (plots against time) | | wlim/llim | numeric vector c(min, max) — restrict the size (x) axis | | ylim | numeric vector c(min, max) — restrict the value (y) axis | | highlight | character vector — draw named species with thicker lines | | total | logical — add a line for the community total. The total of everything the object holds, so it does not change when you select species or hide the resource; on a length axis it is summed at equal length | | background | logical — whether species marked with markBackground() are drawn. They are drawn only when the selection asks for them, always under a single grey "Background" legend entry; background = FALSE removes them | | log_x, log_y | logical — log-scale the x or y axis | | size_axis | "w" (default) or "l" — plot against weight or against length |

wlim/llim (size axis) and ylim (value axis) only set the visible window: data outside the range is hidden but nothing is recomputed. To change the underlying numbers — for example the size range that a biomass is summed over — pass min_w/max_w (or min_l/max_l) to the get…() function instead, e.g. plot(getBiomass(sim, min_w = 10)).

Which arguments apply depends on the array's shape:

• plot(<ArrayTimeBySpecies>) accepts species, tlim, total, background, highlight, log_x, log_y, ylim.
• plot(<ArraySpeciesBySize>) accepts species, highlight, total, background, log_x, log_y, wlim, llim, ylim, size_axis, per_log_size, all.sizes. size_axis, llim and per_log_size belong to the size shapes only — a plot against time has no size axis to convert.
• plot(<ArrayTimeBySpeciesBySize>) takes one time slice and hands it to the ArraySpeciesBySize method, so it accepts everything that method does plus time (default: the last time step). It has no tlim: only one time is shown.
• plot(<ArrayResourceBySize>) accepts log_x, log_y, wlim, llim, ylim, size_axis, per_log_size. The resource is a single spectrum, so there is nothing for species, highlight, total or background to select.
• plot(<ArrayTimeByResourceBySize>) accepts the same as ArrayResourceBySize plus time.

All five also accept return_data = TRUE, which returns the data frame behind the plot instead of the plot, and y_ticks to set the number of y-axis ticks.

What kind of value an array holds

Every mizer array declares what kind of quantity it holds, in its type attribute, because two kinds need handling that the numbers alone do not reveal:

| type | Meaning | What the plots do with it | |---|---|---| | "value" | a rate, an amount — the default | nothing special | | "density" | an amount per gram of body weight | converts the values, not just the axis, when plotted against length | | "proportion" | a fraction | shows the whole of the interval from 0 to 1 on a linear y axis |

Read it with array_type(x), and set it when you build an array of your own:

```{r eval=FALSE} ArraySpeciesBySize(x, value_name = "Number density", units = "1/g", type = "density", params = params)

### Plotting densities

A density is an amount *per unit size*, so its numerical value depends on which
size variable it is a density in. Changing that variable — weight to length, or
size to log size — therefore changes the plotted **values**, not just the axis:
it needs a Jacobian factor. The plot functions apply it for you, for the arrays
that declare themselves densities:

| Source | Density |
|---|---|
| `initialN(params)`, `finalN(sim)`, `N(sim)`, `get_initial_n(params)` | consumer number density, per gram |
| `initialNResource(params)`, `finalNResource(sim)`, `NResource(sim)` | resource number density, per gram |
| `resource_capacity(params)` | resource carrying capacity, per gram |
| `getFluxGradient(params)` | rate of change of the flux, per gram per year |

Which density you get is set by two independent arguments:
`size_axis` chooses the size variable and `per_log_size` chooses whether the
values are per size or per logarithmic size. The factors are built from the
length-weight relationship $w = a\, l^b$ of each species, taken from the `a` and
`b` columns of `species_params`:

| Argument | Factor |
|---|---|---|
| `size_axis = "w"`, `per_log_size = TRUE` | $dw/d\log w = w$ |
| `size_axis = "l"`, `per_log_size = FALSE` | $dw/dl = b\, w / l$ |
| `size_axis = "l"`, `per_log_size = TRUE` | $dw / d\log l = b\, w$ |

**`log_x` does not change the y-axis.** Showing size on a logarithmic axis is
a display choice; you need to use `per_log_size` to convert a density per unit
size into a density per logarithmic size interval. Conflating the two is the
usual reason a spectrum looks like it has the wrong slope.

```r
plot(initialN(params), per_log_size = TRUE)                 # per log weight
plot(initialN(params), size_axis = "l", per_log_size = TRUE) # per log length
plot(initialNResource(params), per_log_size = TRUE)          # resource too

Plotting size spectra

plotSpectra() is the function you want for plots of the abundance or biomass density against size, one line per species. Unlike a plain plot() of a species density array it also overlays the resource spectrum (resource = TRUE, the default). Which density it shows is set by biomass and per_log_size, described below.

By default it shows the final time step of a simulation; pass time_range to average over a period, or give it a MizerParams object to see the current state. The common arguments above all apply, and plotlySpectra() is the interactive twin.

plotSpectra(params)                                   # spectra of the current state
plotSpectra(sim, per_log_size = TRUE, time_range = 1990:2000)
plotSpectra(sim, species = c("Cod", "Herring"), resource = FALSE)
plotSpectra(sim, biomass = TRUE, size_axis = "l")     # biomass density against length

The resource has its own length convention. It is a composite of many taxa, so instead of a taxonomic weight-length relationship it uses the equivalent spherical diameter of an organism with the density of water (a = pi/6, b = 3, in resource_params()). It therefore appears on a length axis, but measured differently from the fish: a fish of a given weight is about 3.7 times longer than a sphere of that weight. That gap at the resource-consumer boundary is real biology, not an artefact.

Which density a spectrum plot shows

plotSpectra(), plotSpectra2() and animate() describe the plotted quantity with two independent logical arguments:

| | per_log_size = FALSE | per_log_size = TRUE | |---|---|---| | biomass = FALSE | number density | number density per log size | | biomass = TRUE | biomass density | biomass density per log size |

The older single power argument is the sum of the two (0, 1, 1, 2 across that table) and is still accepted.

Cumulative distributions

plotCDF(object, species, biomass, normalise) plots cumulative abundance or biomass over size — steadier than a density spectrum for eyeballing where biomass sits. biomass = TRUE (default) accumulates biomass, biomass = FALSE accumulates numbers; normalise = FALSE plots the cumulative total rather than the proportion. The per_log_size argument is not used: a cumulative total does not depend on it.

plotCDF(NS_params, species = c("Cod", "Herring"))
plotCDF(NS_sim, biomass = FALSE, normalise = FALSE)

Comparing two size distributions

| Function | Shows | |---|---| | plotSpectra2(object1, object2, name1, name2) | two abundance spectra overlaid | | plotSpectraRelative(object1, object2) | relative difference of two spectra | | plotCDF2(object1, object2, name1, name2) | two cumulative distributions overlaid |

plotSpectra2(params, params2, "Before", "After")
plotSpectraRelative(params, params2)         # 2 (N2 - N1) / (N1 + N2)
plotCDF2(sim, sim2, "Unfished", "Fished")

Animating through time

animate() plays a spectrum or array through the course of a simulation (animateSpectra() is a retained alias).

animate(sim)                 # abundance spectra over time
animate(getFMort(sim))       # an ArrayTimeBySpeciesBySize over time
animate(NResource(sim))      # an ArrayTimeByResourceBySize over time

animate() accepts most of the common arguments from plot().

Scanning a model over a range — scanModel()

Everything above measures one model. scanModel() measures a family of them: it varies one aspect of the model over a range of values and, at each value, projects until the model settles and measures a quantity on the attractor it settled on. The result is a MizerScan, a data frame that knows how to plot itself.

scan <- scanModel(params,
                  scan_values = seq(0, 1.5, 0.1),
                  set_func    = scanFishingMortality("Cod"),  # what to vary
                  value_func  = getYield,                      # what to measure
                  species     = "Cod")
plot(scan)
attr(scan, "at_max")      # the F at which the yield is largest, i.e. F_MSY

• set_func(params, value) returns a modified MizerParams. Ready-made ones: scanEffort(gear) for fishing effort, scanFishingMortality(species, gear) for the fishing mortality on one species with the rest of the fishing left alone, scanSpeciesParam(species, parameter) for a species parameter. Any function of (params, value) will do, as long as it is idempotent — with continuation = TRUE it is handed the object it returned at the previous scan value, so it must set rather than accumulate.
• value_func(sim) returns a time-by-series object. getBiomass, getYield, getSSB, getN and sizeIntegral all work unchanged, and so does a plain numeric vector over time such as getMeanWeight. Give extra arguments with a closure, not through ...: value_func = function(sim) getBiomass(sim, min_w = 10).

Scanning something that has nothing to do with fishing needs no more than a two-line setter, because project() takes the effort from the params object:

plot(scanModel(params,
               scan_values = 10^seq(10, 12, length.out = 9),
               set_func = function(params, value) {
                   resource_params(params)$kappa <- value
                   params
               },
               scan_name = "Resource capacity", scan_units = "g"),
     log_x = TRUE)

What the band means

How the quantity is measured depends on what the model settled on, which projectUntilSettled() reports:

| Attractor | What is measured | |---|---| | fixed point | read off the settled state, no further projection; ymin == ymax | | limit cycle | averaged over exactly one period; ymin/ymax give the range | | neither | averaged over t_sample years, and the scan values are named in a message |

So the band on the plot is the range of the oscillation, and a Hopf bifurcation appears as the scan value at which it opens up. Averaging over exactly one period is what keeps the curve smooth: a window that is not a whole number of periods leaves a residue of the oscillation in the average. If you need the average more accurately, reduce dt — do not lengthen the window.

Points where the model settled on neither are marked with a cross, and the residual column says how fast the abundances were still changing there, in 1/year. Treat those points as provisional and raise t_max.

plot() takes style = "ribbon" (default: the average as a line inside the band), "envelope" (lines along the edges, no average) or "line" (no band), plus mark_max, reference_lines and the usual log_x/log_y/ylim arguments. A bifurcation diagram over fishing effort is scanEffort() with style = "envelope".

plotYieldVsF()

The one scan common enough to have its own function. plotYieldVsF(params, species) is scanModel() with scanFishingMortality() and getYield(), drawn with the peak marked, so the fishing mortality at the peak is (F_{MSY}):

plotYieldVsF(NS_params, "Cod", F_max = 1.5)
scan <- plotYieldVsF(NS_params, "Cod", F_max = 1.5, return_data = TRUE)
attr(scan, "at_max")      # F_MSY for Cod

The y axis is linear by default, because the yield is exactly zero at F = 0. The current fishing mortality is drawn as a "Current F" reference line. If the species already has an F_MSY species parameter it is also drawn as a reference line, so the value the model gives can be compared with the value assumed.

Dedicated plot functions

Besides the spectrum plots above, mizer has a dedicated plot…() function for each of the common summary quantities. Each is a shortcut for plot() applied to the matching get…() array (e.g. plotBiomass(sim) is plot(getBiomass(sim))). They accept the common arguments above, and each has a plotly…() counterpart (e.g. plotlyBiomass()) for interactive use — the array plot()s use plotHover() instead.

• Against time: plotBiomass(sim), plotYield(sim)
• Against body size (final time step by default, or pass time_range): plotFeedingLevel(sim), plotPredMort(sim), plotFMort(sim)
• Distinct plots:
• plotYieldGear(sim) — yield vs time faceted by gear (one panel per gear)
• plotGrowthCurves(sim) — size at age rather than a size spectrum
• plotDiet(params) — stacked diet composition by prey
• Calibration: plotBiomassObservedVsModel(params) and plotYieldObservedVsModel(params); the latter takes a gear argument that restricts both the modelled and the observed catch to the named gears. See the calibrate-model skill.
• Overview: plot(sim) combines several panels; plot(params) shows the same panels for a model's steady state (without the biomass-through-time panel).

Working with ggplot2

All plotting functions return a ggplot2 object, so you can customise them:

library(ggplot2)
p <- plotBiomass(sim, species = c("Cod", "Herring"))
p + theme_bw() + labs(title = "Biomass through time")
p + geom_hline(aes(yintercept = 1e10), linetype = "dashed")

Species line colours and types come from the linecolour/linetype slots of the MizerParams; change them there for consistent styling across every plot:

params <- setColours(params, c("Cod" = "darkblue"))
params <- setLinetypes(params, c("Cod" = "dashed"))

For interactive exploration prefer the plotly…() twin of a named function, or plotHover() for the compositional array plots.

mizer documentation built on Aug. 31, 2026, 5:08 p.m.