This guide is for people who want to modify Rceattle itself: add a parameter, a data input, a new selectivity or suitability form, or a likelihood component. It is not needed to fit models, for that, start with the user articles (Introduction and Model options and functionality).
It consolidates and updates the developer notes on the GitHub wiki.
Repository layout
| Path | Contents |
|---|---|
R/ |
The R side, 67 files. The number prefix tracks pipeline order for
the 1- to 6- chain; everything under
0- is shared machinery. Every file is listed in Every file in R/,
and what it does. |
src/TMB/ |
The C++ TMB template ceattle.cpp and its included
.hpp modules. |
tests/testthat/ |
Unit and regression tests, grouped by area. |
inst/extdata/ |
Excel templates, including meta_data_names.xlsx (the
authoritative column/switch reference). |
data/ |
Built-in example data lists (BS2017SS,
GOA2018SS, NorthernRockfish2022, …). |
vignettes/ |
User articles; vignettes/articles/ holds website-only
articles like this one. |
.github/workflows/ |
R-CMD-check.yaml and pkgdown.yaml CI. |
inst/dev/ |
Committed developer notes that ship with the package:
TRAPS.md (verified traps with their measured numbers — read
it before touching anything it names), SIBLING-REPOS.md
(the consumer repositories an API change is swept against),
SIMPLIFY-LOG.md (design wishes; a defect never goes here),
and two runnable selectivity examples. The working notes — the session
handoff, the triaged backlog and the per-subsystem
TODO-*.md records — live in a separate maintainer
repository, and the open items that code comments used to cite are now
GitHub issues. |
data-raw/ |
Generators for committed artefacts, notably
regenerate_meta_xlsx.R for
meta_data_names.xlsx. |
tools/verify/ |
22 standalone verification harnesses /golden-check does
not cover — the optimizer path, array bounds, the MSE draw digests, the
simulation recoveries. |
.claude/commands/ |
The slash commands referenced below (/golden-check,
/new-column, /verify,
/ecosystem-sweep). Each is a markdown file stating the full
procedure, usable as a checklist by a person. |
The whole model (observation model, population dynamics, and
likelihood) is a single C++ TMB script,
src/TMB/ceattle.cpp, which #includes a set of
topical .hpp modules. fit_mod() is the R entry
point; it calls a chain of helpers that build the TMB inputs, optimize,
and label the output.
I want to change X — where do I start?
The question most people arrive with. Find the row, start at the
file, and read its row in Every file in
R/, which covers all 67. Where a row names a vignette,
read that first — these rows are orientation, and a recipe is the
job.
| I want to … | Start at | Then usually also |
|---|---|---|
| add a selectivity form |
adding-a-selectivity-form.html — read it
first. It names all the files and the two static tests that
catch a half-finished form; this row is orientation, not the job |
the dispatch is switch (sel_case) in
src/TMB/selectivity.hpp, not
switch (sel_type) — sel_case is derived, and
form 13 under a random walk rewrites itself to -13. Five
files are always touched: R/0-switches.R
(sel_map), R/0-column_schema.R,
R/data.R, R/1-data_check.R,
R/3-build_map.R. Codes 1-9, 11, 13 and 15 are taken; the
next free one is 16 |
| add a likelihood component |
src/TMB/ceattle.cpp §13 and the JnllRow
enum |
R/6-rename_output.R (row labels) and
.JNLL_ROW_AXIS in R/9-profile.R — all three
must agree, and test-schema-jnll-rows.R checks it |
| add a workbook column | R/0-column_schema.R |
not only there. switch_check() fills
defaults through 23 hand-written calls, not a loop, so a column with a
default needs a 24th in R/0-switches.R;
meta = TRUE is the constructor default and
R/data.R’s @format is pinned to the schema by
test-schema-canonical.R; and
meta_data_names.xlsx is a committed file regenerated by
data-raw/regenerate_meta_xlsx.R.
.claude/commands/new-column.md has the full eight-step
recipe |
| add a switch value (a new form name) |
R/0-switches.R (the relevant *_map) |
validate_switches() rejects anything not in the
map |
| change what a parameter starts at | R/2-build_params.R |
R/4-build_parameter_bounds.R if it needs bounds |
| change which parameters are estimated | R/3-build_map.R |
R/6-phaser.R if it should be phased |
| add a covariate or time-varying effect to a process |
linkage_spec() in R/0-build_linkage.R
|
the process’s own build_*(). Nothing in the TMB
template if the parameter already accepts a linkage; a new
parameter needs R/0-linkage_encode.R and
src/TMB/linkage.hpp changed in lockstep, and a new
process or link must be added to
LINKAGE_PROCESSES_IMPLEMENTED /
LINKAGE_LINKS_IMPLEMENTED in
R/0-linkage_table.R or it is refused as reserved |
| add a data input |
R/0-column_schema.R +
R/0-read_write_excel_data.R
|
R/1-data_check.R, R/5-rearrange_data.R,
and the DATA_* block in ceattle.cpp. A
data_list element with no read/write support round-trips to
nothing |
| change a reported quantity |
REPORT()/ADREPORT() in
ceattle.cpp §14 |
R/6-rename_output.R and
R/0-quantity_dictionary.R
|
| change predation / suitability | src/TMB/predation.hpp |
R/3-build_map.R (build_map_predation), and
run the hake MSE — /golden-check does not reach estimated
suitability |
| change growth or weight-at-age | src/TMB/growth.hpp |
R/0-build_growth.R |
| change a reference point or the HCR |
src/TMB/spr.hpp and ceattle.cpp §6.2 for
the SPR proxies; §6.6–6.7 for the depletion-based points and the
projection HCR |
R/0-build_hcr.R, and build_hcr_map() in
the same file for its parameters |
| add a plot |
R/7-plot_helpers.R for the shared resolvers, then the
topical R/7-plot_*.R
|
?"rceattle-plot-args" documents the shared argument
vocabulary once; quantity_dictionary() gives the axis label
and units |
| add an observation-error draw | the SIMULATE{} block beside the density that scores it,
in ceattle.cpp
|
sim_mod() implements no observation model in R; a new
likelihood family owes a draw |
| change projection or the MSE | R/10-run_mse.R |
R/10-mse_summary.R; anything touching row counts or
draw order changes every subsequent random draw |
The fit pipeline
fit_mod() (R/6-fit_mod.R) orchestrates the
following stages. The number-prefixed file names mirror this order.
| Stage | File | Function | Role |
|---|---|---|---|
| 1 | R/0-clean_data.R |
clean_data() |
Coerce and clean the incoming data_list. |
| 2 | R/0-switches.R |
switch_check() |
Fill missing switches with defaults; normalize string/integer switches. |
| 3 | R/1-data_check.R |
data_check() |
Validate inputs; error early on unsupported options. |
| 4 | R/2-build_params.R |
build_params() |
Build the starting parameter list from the switches. |
| 5 | R/3-build_map.R |
build_map() |
Build the TMB map (which parameters are estimated
vs. fixed). |
| 6 | R/4-build_parameter_bounds.R |
build_bounds() |
Lower/upper parameter bounds. |
| 7 | R/5-rearrange_data.R |
rearrange_data() |
Reshape the data for TMB (and build the OSA observation vector when requested). |
| 8 | R/6-fit_mod.R |
fit_mod() |
TMB::MakeADFun() → nlminb() →
TMB::sdreport(). |
| 9 | R/6-rename_output.R |
rename_output() |
Label derived quantities on the returned object. |
Phasing (optional staged estimation) is handled by
set_phases() / TMBphase() in
R/6-phaser.R; TMBAIC()
(R/6-tmb_aic.R) computes AIC. Downstream wrappers
(retrospective(), jitter(),
run_mse(), sim_mod(),
self_test(), the plot_*() family) all take the
fitted Rceattle object returned by
fit_mod().
Before that pipeline runs, a user can assemble or edit the
data_list programmatically and see exactly what a given
model needs, see Data assembly and the column schema
below.
Every file in R/, and what it does
All 67, one row each, grouped by what they are for rather
than by their number prefix — the prefix tracks pipeline order only for
the 1- to 6- chain, and the 31 files under
0- are shared machinery that the chain draws on. The nine
files a fit passes through are repeated here so this table is complete
on its own; The fit pipeline gives them
in order.
A note on vocabulary used below, because two words do double duty in
this package. The map means the TMB map,
the mask saying which parameters are estimated and which are held fixed;
a parameter that is “mapped out” is not estimated. A switch
map (sel_map, q_map) is a different
thing: a lookup from the readable name a user types to the integer the
template consumes. Deviates are the annual offsets that
make a parameter time-varying — estimated per year, and penalized or
integrated as random effects depending on the setting.
Apical selectivity is selectivity at the age or length
where the curve peaks.
Specifications and registries — the source of truth for what a setting means. Change a value here and the rest of the package follows.
| File | What it holds |
|---|---|
0-column_schema.R |
Every workbook column: name, sheet, type, default,
allowed values, old spellings, TMB target, and the one-line
doc that is written verbatim into
meta_data_names.xlsx. Note the ?BS2017SS field
text is not generated from it — that lives in
data.R and is only asserted to agree.
has_default is wired for fleet_control only;
for the 34 control and bioenergetics scalars the default
field is documentation. |
0-switches.R |
The 20 switch maps (sel_map, q_map,
initMode_map, …) mapping a readable name to the template’s
integer, plus switch_check() /
validate_switches() / convert_switches() /
revert_switches(). Also the selectivity-penalty sign
registry .RCE_SEL_PEN_POSITIVE and the
Sel_norm_bin codes. switch_check() fills
defaults through 23 hand-written calls, one per column — not a loop over
the schema. |
0-parameter_dictionary.R |
Internal parameter name → natural-scale name, process, meaning,
dimensions. Exported as parameter_dictionary(). |
0-parameter_index.R |
Parameter → model coordinate: which axes a block’s dimensions mean
(species, sex, fleet, year, age, bin), the slot labels for seven
non-selectivity blocks and for each selectivity form, and the
shortened-axis registry. Exported as
parameter_index(). |
0-quantity_dictionary.R |
Reported quantity → meaning, units, dimensions.
Exported as quantity_dictionary(). |
1-data_requirements.R |
data_requirements(): which inputs a given configuration
requires, which are optional, which are ignored. |
1-data_requirements_table.R |
The requirement catalogue that function reads, plus the per-element classification helpers. |
data.R |
The @format field dictionary for the bundled datasets —
hand-maintained, and pinned name-for-name against
0-column_schema.R’s doc by
test-schema-canonical.R. Adding a schema column with
meta = TRUE and not updating this file is a CI
failure. |
Building a model specification — the
build_*() family, called before fitting.
| File | What it holds |
|---|---|
0-build_data.R |
build_data(): assemble or edit a data_list
in R, supplying only the inputs a configuration uses. |
0-model_config.R |
model_config(): the model-level switch bundle, and its
print method. |
0-build_srr.R |
The stock–recruit relationship: form, priors, fixed or estimated, and its linkages. |
0-build_m1.R |
Natural mortality M1: form, priors, linkages, and the
retired-argument refusals. |
0-build_growth.R |
Growth form (von Bertalanffy, Richards), the population length grid, and growth linkages. |
0-build_selectivity.R |
Selectivity linkages, the apical-offset rows, and which forms accept which linkage. |
0-build_catchability.R |
Catchability linkages and the analytical-q forms. |
0-build_composition.R |
Composition weighting, including the Dirichlet-multinomial weight priors. |
0-build_hcr.R |
build_hcr(), the harvest control rule, and
build_hcr_map(), which maps its parameters. |
0-fit_control.R |
The optimizer, sdreport and phasing knobs, bundled for
fit_mod(). |
0-save_config.R |
save_config() / load_config() /
run_config(): a run’s whole configuration as commented
YAML. |
0-priors.R |
prior_normal() and friends, for linkage
coefficients. |
The linkage system — one formula grammar covering every way a parameter can vary (covariate, time block, random effect), for every process.
| File | What it holds |
|---|---|
0-build_linkage.R |
linkage_spec(), the user entry point; materialization;
the random-effect registry. The largest file in this subsystem. |
0-linkage_formula.R |
Splits a formula into its fixed and random parts. |
0-linkage_table.R |
The long-format coefficient registry every process writes into, and the refusals for a reserved-but-unwired process or link. |
0-linkage_encode.R |
Turns that table into the parallel vectors and the dense design
matrix the template reads. Must stay in lockstep with
src/TMB/linkage.hpp. |
Reading and writing data.
| File | What it holds |
|---|---|
0-read_write_excel_data.R |
read_data(), write_data(),
write_template() — the workbook round-trip, including the
legacy sheet-name aliases. |
0-clean_data.R |
Coerce and tidy an incoming data_list; split the
initial equilibrium catch out of catch_data. |
0-combine_data_sets.R |
combine_data(): merge two data lists. |
The fit chain — the nine stages
fit_mod() runs, in order.
| File | What it holds |
|---|---|
1-data_check.R |
data_check(): validate the inputs and refuse
unsupported combinations. Emits conditions only; returns nothing a fit
reads. The per-form Time_varying_sel restrictions live
here. |
2-build_params.R |
build_params(): the starting parameter list, and the
linkage intercept initial values. |
3-build_map.R |
build_map(): the TMB map — which parameters are
estimated and which are held fixed — assembled per process, plus the
shared-parameter-block adjustment. |
4-build_parameter_bounds.R |
build_bounds(): lower and upper bounds, on the scale
each parameter is stored on. |
5-rearrange_data.R |
rearrange_data(): reshape the data into the arrays the
template expects, and build the OSA observation vector. |
6-fit_mod.R |
fit_mod(): the entry point. Assembles the TMB inputs,
optimizes, projects under the HCR, and labels the output. |
6-phaser.R |
set_phases() / TMBphase(): ADMB-style
phased estimation. |
6-rename_output.R |
rename_output(): labels every derived quantity,
including the jnll_comp row names. Also holds the
McAllister–Ianelli weight calculation
(calc_mcall_ianelli()), which 6-reweight.R
iterates. |
0-tmb_helpers.R |
The nlminb + sdreport wrapper, the
estimability check, and the marginal-gap diagnostic. |
Post-fit: refits, reweighting, reporting.
| File | What it holds |
|---|---|
6-refit_like.R |
.refit_like(): re-fit a model reusing a source
data_list’s structural choices. Eight entry points route
through it — every refit diagnostic, plus run_mse(),
reweight_comps(), remove_F() and
sim_mod(). |
6-reweight.R |
reweight_comps(): iterate composition weights toward
the McAllister–Ianelli target. |
6-tmb_aic.R |
TMBAIC(): marginal AIC. |
6-report_tables.R |
report_tables() and standard_output():
assessment tables and the NOAA standard output. |
0-rceattle_class.R |
The nine S3 methods on a fit — print,
summary, plot, coef,
vcov, logLik, residuals,
as.data.frame, simulate — and the Pearson
residual machinery. |
0-rceattle_data_class.R |
print / summary for a
data_list, rendered as an indented spec tree. |
Diagnostics.
| File | What it holds |
|---|---|
0-convergence.R |
convergence_diagnostics() and its thirteen checks, and
the four gates every refit diagnostic shares (converged, timeout,
dropped, errors). |
0-diagnostic_print.R |
The shared header, table and tally the diagnostic print methods use. |
9-diagnostic_helpers.R |
.parallel_lapply() and the argument vocabulary the
refit diagnostics share. |
9-retro_and_jitter.R |
retrospective() (Mohn’s rho) and
jitter(). |
9-self_test.R |
self_test(): simulate from the fit and check the model
recovers its own parameters. |
9-profile.R |
Likelihood profiles, and profile_components() for the
per-component breakdown. |
6-osa_residuals.R |
One-step-ahead residuals — the only residuals with a known null distribution under random effects. |
0-osa_data.R |
Builds the flat observation vector the OSA machinery needs. Called
from 5-rearrange_data.R. |
6-process_residuals.R |
Residuals for the random-effect processes rather than the observations. |
Simulation, projection and MSE.
| File | What it holds |
|---|---|
8-sim_mod.R |
sim_mod(): draws new data by calling the template’s
SIMULATE blocks — no observation model is implemented in R.
Also sample_rec(), compare_sim(), and R
mirrors of the growth equations (get_growth_matrix_r(),
get_weight_at_age_r()). |
10-run_mse.R |
run_mse(): the operating model, the assessment
schedule, and the estimation-model refit. |
10-mse_summary.R |
mse_summary() performance metrics, plus
check_mse() and load_mse() for runs on
disk. |
10-project-no-F.R |
remove_F(): rerun with F set to zero from a given
year. |
11-model_average.R |
model_average(): average derived quantities across
fits. |
Plotting — the shared argument vocabulary lives in
7-plot_helpers.R and is documented once, in
?"rceattle-plot-args".
| File | What it holds |
|---|---|
7-plot_helpers.R |
The resolvers every plotter should use: species, years, colour, line parameters, faceting, age indexing, the projection divider. |
7-plot_ceattle.R |
plot_timeseries() and the wrappers over it (biomass,
SSB, recruitment, depletion, F), plus selectivity, mortality, maturity,
ration and the predation plots. |
7-plot_comp.R |
Composition and diet-composition fits and residuals. |
7-plot_diagnostics.R |
Index and catch fits, catchability, index residuals. |
7-plot_osa.R |
OSA residual diagnostics. |
7-plot_profile.R |
Likelihood profiles. |
7-plot_data.R |
Timeline of which data enter which likelihood. |
7-plot_stock_recruit.R |
The stock–recruit curve with the fitted points. |
7-plot_sel_vs_mat.R |
Selectivity against maturity, for one species. |
Package infrastructure.
| File | What it holds |
|---|---|
Rceattle-package.R |
Package-level roxygen and the import directives. |
0-deprecate.R |
.rce_deprecated_arg(): the warning a renamed argument
raises. |
Three dictionaries you can just run
Worth knowing before reading any code: three of the registries above are exported functions, and between them they answer “what is this variable and what are its units” without opening the template.
quantity_dictionary() # reported quantity -> meaning, units, dimensions
parameter_dictionary() # internal parameter name -> natural-scale name, process
parameter_index() # parameter -> which axis each dimension meansNote the output names changed in the 5.x series: the quantities older
scripts and the wiki call NByage and
biomassSSB are now N_at_age and
ssb. quantity_dictionary() is the authority on
the current name.
The switch system
Model options are supplied either as strings
("Logistic", "NPFMC") or as the integer codes
the TMB template consumes. R/0-switches.R is the single
source of truth for that correspondence, via four functions:
-
switch_check(): fill missing switches with defaults and normalize. -
validate_switches(): error if any switch is not a known code/string. -
convert_switches(): canonical string → integer code for TMB. -
revert_switches(): integer code → canonical string (backward compatibility for older integer-coded data files).
The mappings live in named vectors: sel_map,
tv_sel_map, q_map, tv_q_map,
comp_loglike_map, fleet_map,
initMode_map, suitMode_map, and
hcr_map. Adding a value to a map is what makes a new string
name legal; validate_switches() rejects anything not
listed.
Data assembly and the column schema
A CEATTLE data_list has many possible inputs, but any
one assessment uses only some of them. A few pieces keep that
manageable.
The column schema (R/0-column_schema.R) defines
every fleet_control column in one place.
.rce_column_schema() records each column’s name, type,
default value, older accepted spellings, and a one-line description. The
rest of the package reads the column from there, so adding a column
means editing this one file:
-
switch_check()fills its default, -
write_data()writes it to the workbook andread_data()reads it back, - the meta sheet and the
?BS2017SSfield descriptions document it, -
data_check()andrearrange_data()refer to it by its current name.
If a column has older accepted spellings, list them in the schema’s
aliases field. Whenever a data_list enters the
pipeline, an older spelling is translated to the current name
automatically, so workbooks and scripts written against earlier column
names keep working (test-schema-canonical.R checks
this).
Both spellings at once is the case to know about, because it is what
a script produces when it assigns an older name to a
data_list that read_data() has already
translated: the assignment creates the old column beside the new one.
They must then hold the same setting (an integer switch code and its
string count as the same) and anything else is an error naming both
values. Nothing is merged and no blank is filled in from the other
column, because NA is a real setting in several of these
(Sel_norm_bin and Sel_cap_bin mean “do not
normalize” and “no cap”; Proj_F_proportion means “no F
apportioned”), so filling could not express clearing a value and would
move a number silently. Until 5.25.0 the older spelling was simply
discarded here, which is how several assessment scripts came to be
running settings they did not ask for.
data_requirements() reports, for a
given model setup, which inputs that model requires, which are optional,
and which are ignored, so you can see what a configuration needs before
fitting rather than hitting a validation error partway through.
build_data() with
model_config() assembles or edits a
data_list in R, supplying only the inputs a given model
uses, and prints the model as a readable outline (print() /
summary()). write_template() writes a blank
workbook with the right sheets and columns for a chosen
configuration.
Run configuration (save / load)
save_config() and load_config() write and
read a run’s full configuration, the model setup plus the estimation and
uncertainty settings from fit_control(), as a plain-text
YAML file. Each field is commented with its meaning, and fields left at
their default are omitted, so two saved runs differ in the file only
where they differ in substance, and the file reviews cleanly in version
control.
Apply a saved configuration with
fit_mod(data_list, config = load_config("run.yaml")). It
fills in only the arguments you did not pass yourself, so an explicit
argument always takes precedence. Every fit also stores the
configuration it ran with, available as run_config(fit).
Formulas, priors, and nested build_*() specifications all
save and reload, so reopening a saved configuration reproduces the same
fit.
The TMB template
ceattle.cpp includes these modules (order matches the
top of the file):
| Module | Responsibility |
|---|---|
helper_functions.hpp |
Shared utilities. |
comp_osa.hpp |
Composition likelihoods and the one-step-ahead (OSA) decomposition. |
growth.hpp |
Growth (von Bertalanffy / Richards); weight- and length-at-age. |
selectivity.hpp |
Selectivity forms, dispatched by switch (sel_case) —
sel_case is Selectivity unless form 13 is
under a random walk, which rewrites it to -13: logistic
(1), Ianelli non-parametric (2), double logistic (3), descending
logistic (4), Hake/Taylor non-parametric (5), 2D/3D AR1 (6/7), double
normal (8), NonParametricPM (9), LogisticPM (11),
NonParametricIntegrable (13), DoubleNormalSS3 (15, Stock Synthesis size
pattern 24, with its own sel_dn6 parameter array). Codes
10, 12 and 14 are unused; the next free code is 16. |
recruitment.hpp |
Stock–recruitment functions. |
bioenergetics.hpp |
Temperature-dependent ration / consumption (Ceq). |
predation.hpp |
Prey suitability and predation mortality (MSVPA type-2, estimated iteratively). |
diet_data.hpp |
Diet data for MSVPA suitability, and the predicted stomach contents
(diet_hat) the diet likelihood scores. |
linkage.hpp |
Applies the encoded linkage tables (covariates, time blocks, random
effects) to each process’s parameters. Their priors and random-effect
densities are scored in ceattle.cpp §13. |
likelihood.hpp |
Likelihood math for §13: one composition or CAAL row
(comp_row_nll()), the linkage prior families, and the
composition tail fold. |
comp_sim.hpp |
Composition draws: multinomial (sequential conditional binomial),
Dirichlet, and Dirichlet-multinomial, used by the SIMULATE
blocks. |
The model opens with the DATA_* block (switches such as
estimateMode, msmMode, suitMode,
initMode, srr_fun, HCR, plus the
data objects), then the PARAMETER_* block, then the
population and observation dynamics, and finally the joint negative
log-likelihood. The likelihood is accumulated in the
jnll_comp matrix, rows are likelihood components, columns
are species/fleet. Each row is addressed by a named constant from the
JnllRow enum (JNLL_INDEX,
JNLL_CATCH, JNLL_COMP, …
JNLL_LINKAGE_RE, with JNLL_N_ROWS dimensioning
the matrix), so a row is referred to by name in the C++ rather than by a
bare integer. The enum has two hand-synced partners:
the display names in R/6-rename_output.R, and
.JNLL_ROW_AXIS in R/9-profile.R, which records
whether a row’s columns count fleets or species. Adding or reordering a
row means updating all three by hand, or a row is classified on the
wrong axis and a survey’s likelihood is reported against a species.
tests/testthat/test-schema-jnll-rows.R reads the template
and asserts the three agree.
Building. src/TMB/compile.R compiles
the model (framework = "TMBad"); TMB and
RcppEigen are LinkingTo dependencies. The
spurious Eigen -Warray-bounds warnings are silenced by a
source #pragma in the .cpp rather than a
compiler flag, since CRAN rejects warning-suppressing flags.
The linkage and priors system
Every process’s parameters share one formula-driven linkage system:
R/0-build_linkage.R, R/0-linkage_encode.R,
R/0-linkage_table.R, and R/0-priors.R.
linkage_spec(formula = ~ x, by = ~ ..., init = , priors = )
is the user entry point, attached through the per-process builders:
-
build_srr(): recruitment (R0,alpha,beta), -
build_M1(): natural mortality (M1), -
build_growth(): growth (K,L1,Linf,m) and the SD endpoints (sd_L1,sd_Linf), -
build_catchability(): catchability (q), -
build_selectivity(): selectivity (slp_asc,slp_desc,inf_asc,inf_desc,coff), -
build_composition(): the Dirichlet-multinomial weights (theta_comp,theta_caal,theta_diet), prior-only.
The same formula grammar covers every way a parameter can vary. The
right-hand side can be a covariate (~ temp), a set of time
blocks (~ cut(Year, ...) (the fixed part goes straight to
model.matrix(), so there is no bespoke block helper), or a
year-to-year random effect) independent (~ (1|Year)), a
random walk (rw(1|Year)), or first-order autoregressive
(ar1(1|Year)). ~ 1 with a priors
list leaves the parameter at one value but places a prior on it. Use
by = to say what the effect varies over,
by = ~ species + sex + age_bin for the population
processes, by = ~ fleet for catchability and
selectivity.
Placing a prior on a selectivity, catchability, or
composition-weight parameter is the ~ 1 case, keyed by
`(Intercept)`, e.g.
build_selectivity(linkages = list(
inf_asc = linkage_spec(~ 1, by = ~ fleet,
priors = list(`(Intercept)` = normal(50, 10)))))The Dirichlet-multinomial composition weights take a prior only (they
have no covariate or random-effect form). See
vignette("environmental-linkages-and-priors") for the full
grammar.
Note that inside priors =, the bare constructors
normal(), lognormal(), gamma(),
beta() are not the exported
prior_*() functions. linkage_spec() evaluates
the priors quosure in a data mask
(rlang::eval_tidy(..., data = .prior_dispatch_mask())) that
binds normal → prior_normal, and so on. Both
spellings work inside linkage_spec(), and the bare form is
intentional, do not “fix” it to prior_normal().
Recipes
Add a new estimated parameter
- Declare the
PARAMETERinceattle.cppand use it in the relevant.hppmodule. - Add a default starting value in
R/2-build_params.R. - Add mapping logic in
R/3-build_map.R(and the matchingbuild_map_*()helper) so it is turned on/off for the right configurations. - Add it to the phasing order in
R/6-phaser.R. - Add it to the parameter dictionary in
R/0-parameter_dictionary.R. It is a registry keyed by parameter name, so a new (or retired) block that misses it fails only later. - If the parameter is a random effect, give it a
SIMULATEdraw beside its density (see Change the likelihood), orsim_mod(process = )will hand back its fitted values as though they were simulated. - (Optional) expose an on/off or random-effect switch as a
fit_mod()argument. - Document it in
inst/extdata/meta_data_names.xlsxand/or thefit_mod()roxygen.
Add a new reported quantity
-
REPORT()it inceattle.cpp(andADREPORT()it if it needs a standard error). AREPORTinside a comment block is not reported, enumerate from a real fit, not from a grep, when checking what a model actually returns. - Give it dimension names in
R/6-rename_output.Rif it is an array. - Add it to the quantity dictionary in
R/0-quantity_dictionary.R, with its units and whether it holds a standard error.test-schema-quantity-dictionary.Rreads the template and a fit and asserts the two agree, so a quantity added here and not there fails a test. - If a quantity is only defined for some configurations, the
per-recruit reference points are computed only under
msmMode = 0, say so in itsmeaning, and gate it inreport_tables(). CEATTLE leaves a number behind rather than a gap, so an ungated quantity reports a plausible wrong figure.
Add or change a data input
There are two cases. Follow the one that matches.
A fleet_control column (e.g. a new
per-fleet setting), add it in one place,
R/0-column_schema.R:
- Add a row to
.rce_column_schema()with the column’s canonical name, type, default, any back-compataliases, and its doc string. Everything downstream (switch_check()default-filling,write_data()/read_data()ordering, the meta sheet, the field dictionary) reads from the schema, so no other file needs editing to make the column persist and round-trip. - Consume the column where it matters (
data_check(),rearrange_data(), or the cpp) by its canonical name, aliases are already upgraded on entry. - If it is a rename of an existing column, put the old spelling in
aliasesand pin it intest-schema-canonical.R.
A whole new data_list element (a new
sheet / array, e.g. an index covariate):
- Add the
DATA_*object inceattle.cppand consume it in the module. -
Add read and write support in
R/0-read_write_excel_data.R. A new element with nowrite_data()/read_data()handling round-trips to nothing through the standard xlsx format; the feature is silently lossy. (This is exactly howindex_covwas lost; it is the single most common data-input bug.) - Add validation in
R/1-data_check.Rand, if the model needs it, a row in thedata_requirementstable (R/1-data_requirements_table.R). - Add reshaping in
R/5-rearrange_data.Rif TMB needs a different layout. - Document it in
meta_data_names.xlsxand/or roxygen.
Add a new switch option (e.g. a selectivity form or
suitMode)
For a selectivity form, Adding a selectivity form traces one existing form through every file and names the tests that catch a half-finished addition. The general shape:
- Implement the case in the module
.hpp(e.g.selectivity.hppcase N,predation.hppsmodeN). - Register the string ↔︎ integer mapping in the relevant map in
R/0-switches.R(e.g.sel_map,suitMode_map); otherwisevalidate_switches()rejects the new name. - If the option is not yet ready for use, block it in
R/1-data_check.R, as the length-based suitability modes are (suitMode %in% c(1, 3, 5)). - Add bounds in
R/4-build_parameter_bounds.Rif the option introduces bounded parameters (see thesuitMode %in% c(1:2)branch). - Add a test named
test-<area>-<topic>.Rintests/testthat/. - Document it.
Change the likelihood
The joint negative log-likelihood is assembled in section 13 of
ceattle.cpp, in the jnll_comp matrix,
addressed by the JnllRow row constant and the species/fleet
column: data likelihoods first, then priors and penalties. Add or modify
the relevant term there, with any math it shares with another term in
likelihood.hpp, then recompile and re-run the regression
tests (below) to confirm the objective changes only as intended.
Add or update the matching SIMULATE block at the
same time. Every observation and every process error is drawn
in a SIMULATE block next to the density that scores it, so
the two are edited together. A simulator that has drifted from its
likelihood does not error, it makes self_test() report
recovery against a data-generating process the model never assumed.
Three rules the existing blocks follow:
-
Draw what the density assumes: including the
bias-correction convention and the scale (log vs. natural). Where a
family has no exact sampler (a multivariate normal truncated to the
positive orthant, say) note it in the draw and warn from
sim_mod()rather than substituting a different distribution. -
Report under a
*_simname. TMB never clears its report environment, so anything reported inside aSIMULATEblock stays visible inobj$report()afterwards. A draw written into the observed-data object and reported under its name is readable as the data it replaced, which is whycomp_sim,caal_simanddiet_simare separate copies. - Do not draw what the model does not define. If two densities score one latent, the AMAK/Ianelli stock-recruit penalty is the live example; there is no single distribution to draw from. Leave it and warn.
simulate_state selects which processes are redrawn
(recruitment, M, growth, catchability, selectivity). Its slot order is
the linkage process code, so a random linkage is gated by
simulate_state(linkage_process(i)) with no translation,
change one and change the other (R/0-linkage_encode.R).
simulate_period is a two-slot vector, and only
slot 0 is read. It gates the process draws, whose
densities cover the hindcast alone. Observations are drawn in every
period deliberately, because run_mse() splices the
negative-Year rows back in as the next assessment’s data,
so gating them would hand the estimation model observations that were
never simulated. Slot 1 is reserved for a projection-period process draw
(WHAM uses its equivalent that way) and is inert today; nothing in
sim_mod()’s interface reaches either slot, and the tests
and tools/verify/ scripts that exercise slot 0 do so
through the internal .sim_draw(period = ). Wiring slot 1 up
means giving the projection-period deviations a density to draw from
first.
Both are DATA_IVECTORs, defaulted in
build_osa_data(), the single funnel every data list passes
through on its way to MakeADFun().
Testing
Tests use testthat in a flat layout: every file is a
top-level test-<area>-<topic>.R (for example
test-selectivity-logisticpm.R,
test-growth-sd-plus-group.R,
test-likelihood-index-covariance.R). Shared
helpers-*.R and fixtures/ sit alongside them.
Notable fixtures and guards:
-
tests/testthat/fixtures/fit_baseline.rdsand the golden-jnllregression test pin the objective-function value so unintended numerical changes surface. -
test-tmb-makeadfun_smoke.Rchecks the model builds and evaluates. -
test-parameter-recovery.Rfits simulated data and checks estimates. -
helpers-make-msm-data.Rbuilds small multi-species fixtures.
Run the suite with devtools::test() (or
testthat::test_local()). After any C++ change, recompile
(src/TMB/compile.R or pkgbuild::compile_dll())
before testing, since the tests load the compiled
.so/.dll.
Branches and releases
Feature branches are cut from dev and merged back into
dev; the maintainer merges dev into
main as a release, following
inst/RELEASE-CHECKLIST.md. Operational assessments pin a
tagged release. What a pull request holds (NEWS.md, the
DESCRIPTION version, the affected vignette), the
commit-message convention and what CI runs where are stated once, in CONTRIBUTING.md.
Debugging tips
-
gradient is of length 1(or similar). The population is usually crashing, numbers-at-age going negative or a divide-by-zero. Check starting values and bounds for the offending process. -
Model will not converge. Inspect
mod$identified(when reported) to see which parameters are hard to estimate, and adjust their starting values. -
Reproducing an external model. Supply
numbers-at-age via the
NByageFixedsheet and selectivity viaemp_selwithSelectivity = "Fixed", setestDynamics = 1, and runfit_mod(estimateMode = 3)to evaluate without re-fitting. (See also the Stock Synthesis conversion article.)
The argument names, and what is deliberately left alone
The diagnostics take the fitted model as object
(5.11.0). Rceattle and fit still work and are
silent; they warn from 5.12.0 and go in 6.0.0. To find your own call
sites now:
options(Rceattle.warn_deprecated_args = TRUE)Three keep their own names on purpose. profile() takes
fitted because it is an S3 method for
stats::profile(), and renaming a method’s first formal away
from its generic fails R CMD check.
mse_summary(mse = ) and
compare_sim(operating_mod = ) take an MSE result and a set
of simulations, not a fitted model, so their names say more than
object would.
The rest of the surface below is known-awkward and deliberately unchanged. It is recorded here so it is not reopened each time someone reads the API cold.
Switch codes are frozen. Renumbering any of them silently
reinterprets every saved run_config YAML, every workbook,
and every fitted .rds in ../Rceattle-models
and ../GOA-ATF-ESP, the same integer would mean something
new with nothing to signal it. This outweighs the readability of a
tidier numbering, so these stay as they are:
-
estimateMode0 = “Estimate” but 1 = “Hindcast”, which reads backwards. -
srr_fun’s valid set is non-contiguous: 0/2/4 are live, 1/3/5 soft-deprecated. -
sel_mapskips 10, 12 and 14. -
initMode1 (Equilibrium) and 5 (OffsetEquilibrium) differ by one term. -
msmModeaccepts 3–9, thendata_check()rejects them.
Overloaded columns. Time_varying_q is a mode
code for most fleets but a 1-based env_data column index
(possibly comma-separated, "1,3") when
Catchability is Environmental (5). Code 6 was
the other such form and is removed, but the overload is not academic: it
is why the Time_varying_q mode checks
(validate_switches(), the soft-deprecation, and the 5.16.0
"AR1" refusal) all exempt those fleets. A check that
forgets it rejects a working model for a value the assessor never set as
a mode.
Sel_curve_pen1 is a penalty weight for selectivity
2/9/11/13 and Sel_curve_pen2 a curvature penalty for
2/9/13, while for 6/7 both are logit-scale AR1 correlations. Re-encoding
either would break existing workbooks; they are documented in the column
schema instead.
Argument counts and names. fit_mod() and
run_mse() both take upwards of twenty arguments, and the
iteration counts are spelled peels, njitter,
nsim and n_iter, each idiomatic where it
appears. fit_mod() also still accepts the optimizer
arguments that moved onto fit_control(), with a warning;
the live list is .deprecated_ctl_args in
R/6-fit_mod.R.
One more collision. compare_sim() has an
argument named object that is a character selector for
which part of the model to compare, not a fitted model. It predates this
release and is left alone for the same reason operating_mod
is, but it is the one place where object means something
else.
profile() is the odd diagnostic out. It refits
like retrospective(), jitter() and
self_test() but takes only a bare getsd=, no
phase, no fit_control, so every profile point
refits at .refit_like()’s unphased default. Worth knowing
when reading a profile that looks flatter than expected.
See also
- GitHub wiki, the original developer notes this page consolidates.
- Model parameterizations, equation-level detail on selectivity, catchability, and predation.
- Adams et al. (2022), Fisheries Research 251:106303, the TMB generalisation of CEATTLE. Holsman et al. (2016), Deep-Sea Research II 134:360–378 (the original model. Wassermann et al. (2025), ICES JMS) the hake/cannibalism extension.