---
title: "5. Single- vs. multi-species models"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{5. Single- vs. multi-species models}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE,
  fig.width = 7,
  fig.height = 5
)
```

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

```{r}
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)
```

```{r}
# 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.

```{r}
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

```{r}
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`.