0. Model options and functionality

1. Model dimensions and operating mode

Equations: see vignette 8 §Population dynamics.

Dimensions and the time horizon are set on the data_list; the population and estimation switches are arguments to fit_mod():

2. Predation: suitMode (per predator)

Equations: see vignette 8 §Predation sub-model (suitability, predation mortality, consumption / bioenergetics).

suitMode sets how prey preference is modelled; length 1 (all predators) or nspp.

Where year-specific diet data are not supplied (Year = 0 in diet_data), suitability is averaged across suit_styr:suit_endyr, the full hindcast by default.

A predicted diet needs estimated suitability and a predation model: suitMode 2, 4 or 6 with msmMode > 0. Under suitMode = 0 the model reads suitability from diet_data rather than predicting stomach proportions, and under msmMode = 0 no predation is computed at all, so in both cases there is no fitted diet to compare against. residuals(source = "diet"), plot_diet_comp(), plot_diet_comp1() and plot_diet_comp2() refuse on such a fit instead of reporting every observation as a complete misfit.

3. Selectivity (per fleet)

Equations: see vignette 8 §Selectivity parameterizations.

Set the Selectivity column of fleet_control; integers and strings both work:

There is no 10: the integer codes are frozen, because renumbering one would silently reinterpret every saved config and every fitted .rds that holds it.

The AR(1) forms reuse the Sel_curve_pen columns as logit-scale correlations, one per dimension: Sel_curve_pen1 across selectivity bins (ages or length bins, per Selectivity_dimension), Sel_curve_pen2 across years, and Sel_curve_pen3 across cohorts for 3DAR1. Before 5.9.0 the 2DAR1 density applied these two the wrong way round, so a fit that started from unequal Sel_curve_pen1 / Sel_curve_pen2 (or that warm-starts from inits produced by an earlier version) begins from a mirrored point. Refit from inits = NULL, or transpose the two, when moving a 2DAR1 model across that boundary. Fits that start both correlations at 0 are unaffected: TMBphase() holds them there, and only the two reported estimates exchange names.

Time-varying behavior is set by Time_varying_sel (per fleet):

Not every mode applies to every selectivity form, because the model only scores a deviation where it defines a density for it:

The two non-parametric modes differ in what they are a deviation of. "IID" scores each annual coefficient deviation about an estimated base curve; "RandomWalk" scores the year-to-year change in the realized selectivity, which is renormalized to mean 1 within each year.

Neither is integrated out, and fit_mod() refuses random_sel = TRUE on both, fit them with random_sel = FALSE, the penalized formulation the AMAK model intends. The walk’s density says nothing about the level of a year’s coefficients, so those directions are improper and the deviation standard deviation collapses to zero. The IID density is proper, but the AMAK shape penalties (Sel_curve_pen1, Sel_curve_pen2) and the always-on average-selectivity term are charged on each year’s realized curve, so the prior actually integrated is the IID normal times those penalties. None of them scales with Time_varying_sel_sd, so the estimated sd is biased low whatever the penalty weights are (measured at about 5% of the precision at sd 0.35 from the average-selectivity term alone, more as the sd grows), and Sel_curve_pen1 is one-sided, which leaves the Laplace objective only piecewise smooth. NonParametricIntegrable (5.40.0) is the form to fit with random_sel = TRUE: charging the shape penalties once on the base curve makes the mode integrable without moving any penalized fit.

Even there, the Laplace marginal likelihood underestimates a variance component on small multinomial samples: on Atka2022’s fishery (input sample sizes 2 to 236) a true deviation SD of 0.35 comes back at 0.24 to 0.26, 0.70 near 0.61, and 0.35 near 0.32 with ten times the sample size (tools/verify/verify-sim-recovery-np-integrable.R). Read the reported SD as a lower bound unless the compositions are well sampled.

Time_varying_sel is soft-deprecated as a whole, and data_check() says so on use. The linkage grammar expresses the same IID, random-walk and AR(1) deviations, (1 | Year), rw(1 | Year), ar1(1 | Year), and time blocks as ~ cut(Year, ...), per selectivity parameter rather than per fleet, and additionally lets the deviation standard deviation be estimated, given a prior, or held fixed:

build_selectivity(linkages = list(
  inf_asc = linkage_spec(~ ar1(1 | Year), by = ~ fleet)))

Name whichever parameter should vary: slp_asc, inf_asc, slp_desc, inf_desc, or a DoubleNormal alias. The GOA pollock 2025 assessment runs that way. apical is the one linkage parameter every estimated form has: the height of one sex’s curve relative to the other’s, see Environmental linkages and priors.

Bin indexing: with non-parametric and AR(1) forms, selectivity is indexed over N_sel_bins.

Sel_norm_bin says where selectivity is normalized to 1. Set it to "Max" to normalize by the largest value (non-differentiable), to "Off" to leave the curve unnormalized, or to a bin to normalize there. Give Sel_norm_bin_upper as well and the reference is mean selectivity over that range instead.

A numeric Sel_norm_bin is read on the fleet’s own Selectivity_dimension: an absolute age for an age-based fleet (6 means age 6, not the sixth bin) or a 1-based length-bin ordinal for a length-based one. Bin_first_selected beside it uses the opposite convention and is always a bin ordinal. Below Bin_first_selected the curve is zeroed, so a reference taken there would divide by nothing; such a value is read as "Max", which is what a negative has always meant. Legacy numbers still work: blank is "Off", and 0, -1 or -999 are "Max".

mirror <- fleet_control[fleet_control$Fleet_name == "Fishery", ]
mirror$Fleet_name         <- "Fishery_CPUE"
mirror$Fleet_code         <- nrow(fleet_control) + 1L
mirror$Fleet_type         <- "Survey"
mirror$Catchability_index <- mirror$Fleet_code
fleet_control <- rbind(fleet_control, mirror)

4. Catchability (per fleet)

Equations: see vignette 8 §Catchability parameterizations.

Set the Catchability column of fleet_control:

build_catchability(linkages = list(
  q = linkage_spec(~ ar1(1 | Year), by = ~ fleet)))

build_catchability(linkages = list(
  q = linkage_spec(~ ar1(1 | Year), by = ~ fleet, fleet = <Fleet_code>,
                   observe = "<that fleet's env_data column>",
                   obs_sd  = <its measurement SD>)))

[2] Fishery │ Fishery │ sel: Logistic │ q: Estimated (fixed: no index data)

5. Composition likelihoods

Set Comp_distribution (or CAAL_distribution) per fleet:

Diet composition uses an analogous switch on the bioenergetics control sheet: Diet_distribution = 0 (multinomial) or 1 (Dirichlet-multinomial). Conditional age-at-length (CAAL) data go through their own data path and are weighted via CAAL_weights per fleet.

6. Recruitment / stock-recruit (recFun = build_srr())

Priors and environmental relationships go on the curve through build_srr(); see ?build_srr and vignette("environmental-linkages-and-priors").

Under multispecies mode (msmMode > 0) spawning biomass per recruit is undefined, because mortality includes predation. A curve fitted in the hindcast then estimates its initial recruitment level instead of deriving it, and the choice of initMode can move the fitted stock substantially; the curve can also enter as a recruitment penalty (srr_fun = "mean" with the curve as srr_pred_fun). A prior on Beverton-Holt steepness is refused.

7. Natural mortality (M1Fun = build_M1())

8. Growth and weight-at-age (growthFun = build_growth())

9. Harvest control rules (HCR = build_hcr())

DynamicHCR = TRUE switches from static to dynamic SB₀ reference points. HCRorder controls the order in which species F is solved for iterative multi-species HCRs (e.g. predators before prey). See ?build_hcr for more details.

10. Projection and MSE

run_mse() parallelizes over simulations, on a FORK cluster where the platform allows one and PSOCK on Windows. The MSE vignette has the seed plumbing and the reproducibility guarantees.

11. Diagnostics and inference

12. Plotting

Pop / mortality / fit / predation / biology plots, all named after the quantity they show: plot_biomass, plot_ssb, plot_recruitment, plot_depletion, plot_ssb_depletion, plot_exploitable_biomass, plot_catch, plot_f, plot_mortality, plot_m_at_age, plot_m2_at_age_prop, plot_index, plot_index_residual, plot_comp, plot_data, plot_form, plot_b_eaten, plot_b_eaten_prop, plot_diet_comp, plot_ration, plot_selectivity, plot_selectivity_vs_maturity, plot_maturity, plot_stock_recruit, plot_timeseries. Most accept a single fitted model or a list of models with model_names.

These are built with ggplot2 (colorblind-safe Okabe-Ito palette for discrete series, viridis for continuous ones; theme_classic with a panel border for the trajectory figures, theme_bw for the OSA and composition panels) and return the ggplot object, so they print when called interactively and take further ggplot2 layers (e.g. plot_biomass(fit) + ggplot2::ggtitle("My title")). Pass file = "stem" to also write the figure to a PNG. The suffix is per-plotter and not derived from the function name – plot_recruitment() writes stem_R_trajectory.png, plot_catch() stem_fishery_catch.png, plot_data() stem_data_plot.png – and the per-fleet plotters write one file each: plot_comp() writes stem_comp_pearson.png plus stem_comp_annual_fleet<N>_<age|length>.png and stem_comp_aggregated_fleet<N>_<age|length>.png. tests/testthat/test-plot-save-paths.R records every name.

Shared arguments. The timeseries plotters (plot_biomass, plot_ssb, plot_recruitment, the depletions, plot_exploitable_biomass, plot_f), the predation plotters (plot_b_eaten, plot_b_eaten_prop, plot_m_at_age, plot_m2_at_age_prop, plot_ration) and plot_selectivity take one common vocabulary, documented together in ?"rceattle-plot-args": line_col, lwd, lty, alpha, species, spnames, minyr, maxyr, incl_proj, incl_mean, add_ci, model_names, file, width and height. Each word means the same thing wherever it appears, but not every plotter takes every one: incl_mean is on the predation plotters, add_ci only where the quantity holds standard errors, and alpha only where the figure has a ribbon or a fan. Check args(). The other plotters still take their own arguments.

• line_col takes colour names, hex codes, or base-graphics palette indices (line_col = 1), and lwd keeps the base-graphics scale where the default 3 is a standard-weight line.
• line_col and lty supply values for whichever variable the figure separates by colour or line type, which is not always the model. In plot_b_eaten_prop() colour separates predators; in plot_ration() line type separates the sexes. Each function’s help says which. Too few colours are recycled, with a warning naming what they coloured. Where colour is a magnitude (the year fan in plot_selectivity()) line_col gives the ramp anchors, so line_col = "black" draws the fan in black.
• species selects (by index, name, logical mask, or "all", in the order given) and spnames labels.
• minyr and maxyr narrow the data, not just the axis, so a panel with a free y scale rescales to the window.
• Arguments left over from the base-graphics version, right_adj, top_adj, mod_cex, legend.pos, single.plots, theta, are still accepted, and ignored, by the plotters that had them; they were never added to the ones that did not. Set the equivalent on the returned ggplot instead, e.g. plot_biomass(fit) + ggplot2::theme(legend.position = "right").

Display units. The model works in mt and thousands of fish: numbers-at-age are in thousands and weight-at-age in kg, so biomass and SSB come out in mt and recruitment in thousands of fish. The timeseries plotters rescale to million mt and millions for the axis, and derive the recruitment age from minage rather than assuming age 1. Supply the model’s inputs on that convention, catch and index in mt, weight-at-age in kg, or the axis labels will not describe what is plotted. Pass ylab to override the derived label.

Confidence intervals. add_ci = TRUE draws a 95% interval. Strictly positive series take it on the log scale, exp(log(x) ± 1.96 · sd_log), using the model’s log_biomass / log_ssb / log_R where they exist and the delta-method identity sd(log x) = sd(x)/x where they do not. The interval is therefore right-skewed and cannot reach zero for a weak year class or a depleted stock. See ?plot_timeseries for the full derivation.

plot_selectivity() takes it too, but only on a fit run with fit_control(selectivity_se = TRUE) and only for estimated, age-based fleets; otherwise it warns and draws none. Selectivity is identified only up to a scalar, so Sel_norm_bin fixes the curve’s level rather than measuring it: the pinned bin has a standard error of 0 and the band widens away from it. Quote that anchor with the interval.

The error comes from whichever sdreport the fit ends on, and that is not always the one that estimated the curve. Under estimateMode = "Estimate" with an estimating HCR, fit_mod() re-optimizes the projection with every hindcast parameter fixed, so the standard error of every selectivity is exactly 0. estimateMode = "Projection" estimates no selectivity at all, and projection_uncertainty cannot rescue it. Fit the hindcast alone, or ask for uncertainty in the projection:

# Errors from the hindcast fit, the usual choice for a selectivity figure.
mod <- fit_mod(data_list = data, estimateMode = "Hindcast",
               fit_control = fit_control(selectivity_se = TRUE))
plot_selectivity(mod, add_ci = TRUE, minyr = 2015, maxyr = 2015)

# Or keep the HCR projection and reuse the hindcast parameters in it.
mod <- fit_mod(data_list = data, HCR = build_hcr(HCR = 5),
               fit_control = fit_control(selectivity_se = TRUE,
                                         projection_uncertainty = TRUE))

fit_mod() warns when the combination would report zeros, and plot_selectivity() declines to draw a band that is 0 wide at every age rather than render it as certainty. Give minyr / maxyr a single year on a time-varying fleet: every year drawn gets its own band.

Two further conditions belong with any interval you report:

• The normalization. On Atka2022 the survey pins age 4, and its standard deviation runs 0.33, 0.18, 0.12, 0.00, 0.14 over ages 1–5; the fishery pins nothing, normalizes to a mean of one, and is flattest mid-range instead. Same fit, same data, the shape of the band is a property of Sel_norm_bin, not of the information in the data.
• The smoothing. Under random_sel = FALSE the annual deviations are penalized rather than integrated, so these are penalized-likelihood errors conditional on Sel_curve_pen1 / Sel_curve_pen2 and on a fixed Time_varying_sel_sd. Change the smoothing and both the curve and its band move, with nothing in the data to arbitrate.

The band is marginal. Selectivity error is strongly correlated with F and q through the same confounding, so propagating it into a derived quantity needs the covariance (fit_control(getReportCovariance = TRUE)), not this interval alone.

References

• Adams, G.D., Holsman, K.K., Barbeaux, S.J., Dorn, M.W., Ianelli, J.N., Spies, I., Stewart, I.J., Punt, A.E. (2022). An ensemble approach to understand predation mortality for groundfish in the Gulf of Alaska. Fisheries Research, 251, 106303.
• Holsman, K.K., Ianelli, J., Aydin, K., Punt, A.E., Moffitt, E.A. (2016). A comparison of fisheries biological reference points estimated from temperature-specific multi-species and single-species climate-enhanced stock assessment models. Deep Sea Research Part II, 134, 360-378.
• Wassermann, S.N., Adams, G.D., Haltuch, M.A., Kaplan, I.C., Marshall, K.N., Punt, A.E. (2025). Even low levels of cannibalism can bias population estimates for Pacific hake. ICES Journal of Marine Science, 82(1), fsae064.
• Cheng, M.L.H., et al. (2024). 2D / 3D AR(1) selectivity (citation in vignette 8).
• Rogers, L.A., et al. (2024). AR(1) catchability for GOA pollock (citation in vignette 8).