5. Single- vs. multi-species models

2026-05-11

The key switch between model types is msmMode in fit_mod():

msmMode	Predation mortality (M2)
0	Off — single-species; total M = M1
1	Type II MSVPA (Holsman et al. 2015)
2	Type III MSVPA

In single-species mode the model is a standard statistical catch-at-age model. In multi-species mode, a bioenergetics-based predation component (M2) is added to natural mortality so that predator abundance directly affects prey survival. The two modes share the same data format and most fit_mod() arguments — switching between them is a one-line change.

When to use each mode

Single-species is appropriate when:

• Predation interactions are unknown or negligible
• The species of interest has no major predators or prey in the system
• You are building an initial model or debugging data inputs

Multi-species is appropriate when:

• Predation is an important source of mortality (e.g., groundfish in the Bering Sea or Gulf of Alaska)
• You want to evaluate how climate-driven changes in predator or prey abundance propagate through the food web
• Reference points should account for ecosystem effects

A common workflow is to fit single-species models first and use those estimates as starting values (inits) for the multi-species run.

Bering Sea example

library(Rceattle)

# BS2017SS and BS2017MS use the same data structure but differ in M1_base:
# the multi-species dataset has lower residual mortality because some of
# that mortality is explained by predation (M2) in multi-species mode.
data(BS2017SS)
data(BS2017MS)
BS2017SS$projyr <- 2060
BS2017MS$projyr <- 2060

plot_data(BS2017SS)

# Single-species
ss_run <- fit_mod(
  data_list    = BS2017SS,
  inits        = NULL,
  file         = NULL,
  estimateMode = 0,
  random_rec   = FALSE,
  msmMode      = 0, # <-- single-species
  fit_control = fit_control(
    phase        = TRUE,
    verbose      = 1))

# Multi-species — start from single-species MLEs for stability
ms_run <- fit_mod(
  data_list    = BS2017MS,
  inits        = ss_run$estimated_params,  # warm start
  file         = NULL,
  estimateMode = 0,
  niter        = 5,        # iterations between population and predation dynamics
  random_rec   = FALSE,
  msmMode      = 1,        # <-- Type II MSVPA
  suitMode     = 0, # empirical suitability from diet data
  fit_control = fit_control(
    verbose      = 1))

The niter argument controls how many times the population dynamics and predation mortality equations are iterated to convergence within each gradient step. Five iterations is typically sufficient.

Gulf of Alaska example

GOA2018SS is a combined three-species data list (pollock, arrowtooth flounder, Pacific cod) for the 2018 GOA assessment. The same dataset is used for both single- and multi-species runs.

data("GOA2018SS")
plot_data(GOA2018SS)

# Single-species
goa_ss <- fit_mod(
  data_list    = GOA2018SS,
  inits        = NULL,
  file         = NULL,
  estimateMode = 0,
  random_rec   = FALSE,
  msmMode      = 0,
  fit_control = fit_control(
    phase        = TRUE,
    verbose      = 1))

# Single-species with estimated M
goa_ss_M <- fit_mod(
  data_list    = GOA2018SS,
  inits        = goa_ss$estimated_params,
  file         = NULL,
  estimateMode = 0,
  M1Fun        = build_M1(M1_model = c(1, 2, 1)),  # Estimate M for spp 1 & 3
  random_rec   = FALSE,
  msmMode      = 0,
  fit_control = fit_control(
    phase        = TRUE,
    verbose      = 1))

# Multi-species
goa_ms <- fit_mod(
  data_list    = GOA2018SS,
  inits        = goa_ss$estimated_params,
  file         = NULL,
  estimateMode = 0,
  M1Fun        = build_M1(M1_model = c(1, 2, 1)),
  niter        = 3,
  random_rec   = FALSE,
  msmMode      = 1,
  suitMode     = 0,
  fit_control = fit_control(
    phase        = TRUE,
    verbose      = 1))

Predator-prey suitability (suitMode)

suitMode determines how prey preference is estimated for each predator species. It can be a single value (same for all predators) or a vector of length nspp.

suitMode	String	Description	Available?
0	"Empirical"	Empirical proportions from observed stomach content data (MSVPA-style; Holsman et al. 2015)	Yes
1	"GammaLength"	Gamma suitability on prey length-at-age	Not yet
2	"GammaWeight"	Gamma suitability on prey weight-at-age	Yes
3	"LognormalLength"	Log-normal suitability on prey length-at-age	Not yet
4	"LognormalWeight"	Log-normal suitability on prey weight-at-age	Yes
5	"NormalLength"	Normal suitability on prey length-at-age	Yes
6	"NormalWeight"	Normal suitability on prey weight-at-age	Yes

Values 1 and 3 are blocked at runtime pending validation of the growth-model integration. Use suitMode = 0 (empirical) when diet composition data are available, or suitMode = 2/4/5/6 (parametric) when they are not.

Comparing outputs

mod_list  <- list(ss_run, ms_run)
mod_names <- c("Single-species", "Multi-species")

# Population trajectories
plot_biomass(Rceattle = mod_list, model_names = mod_names)
plot_ssb(Rceattle = mod_list, model_names = mod_names)
plot_recruitment(Rceattle = mod_list, model_names = mod_names, add_ci = TRUE)
plot_depletionSSB(Rceattle = mod_list, model_names = mod_names)

# Mortality breakdown (M1 vs M2) — these accept model lists
plot_m_at_age(Rceattle = mod_list, model_names = mod_names, age = 2)
plot_m2_at_age_prop(Rceattle = mod_list, model_names = mod_names)

# plot_mortality() accepts only a single model; call separately per model
plot_mortality(ms_run)   # tile plot of total mortality-at-age across years
plot_mortality(ms_run, type = 2)  # facet line version
plot_mortality(ms_run, M2 = FALSE) # total M (M1 + M2); M2 = TRUE plots M2 only

# Predation (multi-species only) — these also accept single models only
plot_b_eaten(Rceattle = mod_list, model_names = mod_names)
plot_b_eaten_prop(Rceattle = mod_list, model_names = mod_names)

For model averaging across single- and multi-species variants, see ?model_average.