This vignette fits a sequence of five single- and multi-species
models with Rceattle, applied to walleye pollock,
arrowtooth flounder, and Pacific cod in the Gulf of Alaska and eastern
Bering Sea. The model follows Holsman et al. (2016) and Adams et
al. (2022). Every function has help available via ?. The
workflow is:
Install Rceattle
Load Rceattle and associated excel data
Fit assessment models:
model 1: a single-species assessment model estimated with recruitment estimated as penalized effects;
model 2: three single-species assessment models jointly estimated with recruitment estimated as penalized effects;
model 3: as model 2, but with natural mortality estimated for each species;
model 4: as model 3, but with recruitment estimated as random effects;
model 5: full multi-species assessment models with time- and age- varying natural mortality due to predation.
Compare models by AIC and Mohn’s rho (retrospective analysis)
Review plots of input data, diagnostics, and results.
Rceattle and associated dependencies and can be
installed as follows:
# Not run: this installs the package the vignette is already using, so
# executing it in CI would reinstall mid-render.
# Rceattle (pulls CRAN dependencies automatically; e.g. TMB, Matrix, dplyr)
install.packages("remotes")
remotes::install_github("afsc-assessments/Rceattle")
# Optional: TMBhelper provides richer optimization diagnostics.
# Rceattle falls back to plain nlminb + sdreport if it's not installed.
# devtools::install_github("kaskr/TMB_contrib_R/TMBhelper")Rceattle needs a data object to run. Data from the 2018
Gulf of Alaska and 2017 eastern Bering Sea single- and multi-species
assessments are bundled with the package:
library(Rceattle)
data("GOApollock") # Data for GOA-pollock (1-species) ?GOApollock for more information on the data
data("BS2017SS") # Combined data for 3 species. ?BS2017SS for more information on the data
data("BS2017MS") # Multi-species data. Note: the only difference is the residual mortality (M1_base) is lowerData can be visualized using the plot_data function:
The data files can be written to excel, modified within excel, and
read back into R as follows:
Rceattle::write_data(data_list = BS2017SS, file = "BS2017SS.xlsx")
mydata <- Rceattle::read_data( file = "BS2017SS.xlsx")A description of the inputs and file structure can be found on the
first sheet (meta_data) and the associated help function
?BS2017SS.
Models are fit with the fit_mod function. The first
model is a single-species fit with no stock–recruit relationship to Gulf
of Alaska pollock. It has 8 fleets: 1 fishery, 6 surveys, and 1 fleet
turned off (its data excluded from the likelihood; see
model_1$fleet_control).
model_1 <- Rceattle::fit_mod(data_list = GOApollock,
inits = NULL, # Initial parameters = 0
estimateMode = 0, # Estimate
random_rec = FALSE, # No random recruitment
fit_control = fit_control(
phase = TRUE,
verbose = 1))
summary(model_1)The second model fits three single-species models jointly, with no
stock–recruit curve and natural mortality fixed at the input values
(mydata$M1_base). It has 7 fleets: 3 fisheries (one per
species) and 4 surveys (see mydata$fleet_control).
Single-species mode is the default (msmMode = 0):
model_2 <- Rceattle::fit_mod(data_list = mydata,
inits = NULL, # Initial parameters = 0
estimateMode = 0, # Estimate
random_rec = FALSE, # No random recruitment
msmMode = 0, # Single-species
fit_control = fit_control(
phase = TRUE,
verbose = 1))
summary(model_2)Alternatively, natural mortality can be estimated with the
build_M1 function;
M1_model = "sex_age_invariant" estimates a single M per
species, constant across sex and age:
model_3 <- Rceattle::fit_mod(data_list = mydata,
estimateMode = 0, # Estimate
M1Fun = build_M1(M1_model = "sex_age_invariant"), # Estimate sex-invariant M)
random_rec = FALSE, # Penalized likelihood
msmMode = 0, # Single species mode
fit_control = fit_control(
phase = TRUE,
verbose = 1))
summary(model_3)model_4 treats recruitment as random effects via
random_rec. It starts from the previous model’s MLEs and
skips phasing to reduce run-time:
model_4 <- Rceattle::fit_mod(data_list = mydata,
inits = model_3$estimated_params, # Start estimation at model 3's MLEs
estimateMode = 0, # Estimate
M1Fun = build_M1(M1_model = "sex_age_invariant"),
random_rec = TRUE, # Random recruitment
msmMode = 0, # Single species mode
fit_control = fit_control(
phase = FALSE,
verbose = 1))
summary(model_4)Setting msmMode = 1 estimates time- and age-varying
natural mortality from predation, running the model in multi-species
mode. This follows the MSVPA parameterization of Magnusson
(1995).
model_5 <- Rceattle::fit_mod(data_list = BS2017MS,
inits = model_3$estimated_params, # Initial parameters from single species estimates
M1Fun = build_M1(M1_model = "sex_age_invariant"),
estimateMode = 0, # Estimate
niter = 3, # 3 iterations around population and predation dynamics
random_rec = FALSE, # No random recruitment
msmMode = 1, # MSVPA based
suitMode = 0, # Empirical suitability
fit_control = fit_control(
verbose = 1))
summary(model_5)fit_mod() returns an object of class
Rceattle. Rather than reaching into it by hand, most of
what you need is reachable through standard S3 methods:
summary(model_1) # reference points, likelihood, convergence
coef(model_1) # estimated parameters
head(vcov(model_1)) # covariance matrix (NULL when getsd = FALSE)
logLik(model_1) # log-likelihood, holding its degrees of freedom
AIC(model_1) # so stats::AIC() works
head(residuals(model_1))
head(as.data.frame(model_1))plot(), simulate() and
profile() dispatch on the fit too, and the diagnostics
return their own classes with print() methods that
summarize the result, retrospective() gives an
Rceattle_retro, jitter() an
Rceattle_jitter, self_test() an
Rceattle_selftest, profile() an
Rceattle_profile, and osa_residuals() an
rceattle_osa. A data list read by read_data()
is an Rceattle_data, which also prints and summarizes.
The numbers themselves live in two places. Derived quantities on
their natural scale are in fit$quantities, biomass,
ssb, N_at_age, M_at_age,
sel_at_age, reference points and the likelihood components.
The parameters the optimizer actually estimated are in
fit$estimated_params, mostly on a log scale, under
abbreviated names. parameter_dictionary() explains those
names:
parameter_dictionary("R_log_sd") # log_sd of recruitment -> sigma_R
parameter_dictionary(process = "selectivity")Anything beginning log_, or ending _log_sd,
needs exp() to reach the quantity named in the dictionary’s
natural column.
quantity_dictionary() does the same job for the reported
quantities, what each one means, the units it is in, its dimensions, and
whether sdreport() gives it a standard error:
quantity_dictionary("ssb_depletion")
quantity_dictionary(process = "reference_points")
# Everything that holds a standard error
dict <- quantity_dictionary()
dict[dict$se, c("quantity", "units")]Units follow the model throughout: numbers-at-age are
thousands of fish and weight-at-age is
kg, so their product is mt. Two things
worth reading off the dictionary before using a quantity.
ssb is female spawning-stock biomass, and every
per-recruit reference point (SPR0, SPRlimit,
SPRtarget) is defined only in single-species mode, on a
multispecies fit the model leaves them at zero, which is not an estimate
of zero.
Models can be compared by marginal or joint negative log-likelihood, AIC, and retrospective analysis. AIC cannot compare a penalized-likelihood fit against a marginal-likelihood fit. The likelihood components are stored in the model objects built above:
# Evaluate 1 model
AIC(model_1) # AIC, via the logLik() method below
model_1$opt$objective # Negative log-likelihood
model_1$quantities$jnll_comp # Negative log-likelihood components
# Compare AIC across 3-species penalized likelihood models
sapply(list(model_2, model_3, model_5), AIC)report_tables() collects the quantities an assessment
reports, the fit and its likelihood decomposition, the time series with
uncertainty, the reference points, and the fits to index and catch, into
one set of tidy tables. Every table holds a model column,
so several fits stack for a like-for-like comparison:
Read the basis column of reference_points
before quoting anything from it. A CEATTLE fit leaves a number
behind for reference points it never estimated, and each one looks like
an estimate: Ftarget and Flimit sit at 1.0/yr
unless the harvest control rule estimates them, and in multispecies mode
SB0 is a placeholder until a no-fishing projection derives
it, which would make B40% a few hundred tonnes. Those cells are returned
as NA with the reason in basis, rather than as
the number the model happens to hold.
Diagnostics are attached rather than computed: a retrospective or a
jitter is tens to hundreds of refits, so report_tables()
never runs one on your behalf. Pass them in once you have them, and the
matching sections appear:
standard_output() relabels the result into the NOAA
standardized assessment output that the stockplotr and
asar packages read, so Rceattle results can be plotted and
written into a report by the same tooling used for SS3, BAM, WHAM and
FIMS. That standard describes one stock and has no species dimension, so
a fit with more than one needs species to pick the stock to
emit:
head(standard_output(tabs, model = "pollock"))
head(standard_output(tabs, model = "three_species", species = 1))The retrospective function measures retrospective bias,
returning a list of refit models and a data.frame of Mohn’s rho
values.
model_1_retro <- retrospective(model_1, peels = 7)
# `mohns` is long: columns Object, `Forecast year` (0 = terminal peel), N, then one
# column per species. Read the rho by name -- column 2 is the peel year, not a rho.
rho <- subset(model_1_retro$mohns, Object == "biomass" & `Forecast year` == 0)
plot_biomass(model_1_retro$Rceattle_list,
model_names = paste("Pollock; mohns =", round(rho[["Pollock"]], 3)))Rceattle includes plots for hindcast and forecast time
series of biomass, SSB, recruitment, biomass consumed as prey,
mortality-at-age, and depletion, with uncertainty. Time-series plots
accept a list of models (see above).
model_list <- list(model_2, model_3, model_4, model_5)
model_names <- list("single-spp", "Single-spp M", "Single-spp M - RE", "Multi-spp")
# Biomass
plot_biomass(model_list, model_names = model_names, incl_proj = TRUE, add_ci = TRUE)# Recruitment
plot_recruitment(model_list, model_names = model_names, incl_proj = TRUE, add_ci = TRUE)
# SSB
plot_ssb(model_list, model_names = model_names, incl_proj = TRUE, add_ci = TRUE)
# Biomass depletion (total biomass / B0)
plot_depletion(model_list, model_names = model_names, incl_proj = FALSE, add_ci = FALSE)
# SSB depletion (female spawning biomass / SB0), the quantity a Tier 3 HCR
# compares against B40%
plot_ssb_depletion(model_list, model_names = model_names, incl_proj = FALSE, add_ci = FALSE)Diagnostic plots cover selectivity, composition fits, index fits, and catch fits. Composition plots take a single model; selectivity, index and catch plots take one or more.
Each fleet is drawn on the dimension its selectivity was estimated
on, age, or length bin for a fleet whose
Selectivity_dimension is "Length". With one
model, colour is the year, so time-varying selectivity reads as a fan.
With several, colour separates the models and the fan moves to
transparency.
# Selectivity
plot_selectivity(model_2)
# Overlay models; `line_col` and `lwd` apply as in the other plotters
plot_selectivity(model_list, model_names = model_names, lwd = 2)
# Plot composition fit
plot_comp(model_2)
# Plot index
plot_index(model_list) # - All models
# Plot catch data
plot_catch(model_list)Outputs of the model can queried using the $ sign:
model_1$initial_params - list of initial parameter
values prior to estimation
model_1$bounds - list of parameter bounds
model_1$map - list of parameter map in factor form
(model_1$map$mapFactor) or list form
(model_1$map$mapList)
model_1$estimated_params - list of estimated
parameters post estimation
model_1$quantities - list of derived quantities
(SSB, biomass, selectivity, etc)
model_1$data_list - the input data
model_1$run_time - model run time
model_1$obj - the TMB makeadfun object.
model_1$opt - the fit_tmb optimization
object.
model_1$sdrep - the TMB sdreport object.