Skip to contents

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 means

Note 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:

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:

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

  1. Declare the PARAMETER in ceattle.cpp and use it in the relevant .hpp module.
  2. Add a default starting value in R/2-build_params.R.
  3. Add mapping logic in R/3-build_map.R (and the matching build_map_*() helper) so it is turned on/off for the right configurations.
  4. Add it to the phasing order in R/6-phaser.R.
  5. 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.
  6. If the parameter is a random effect, give it a SIMULATE draw beside its density (see Change the likelihood), or sim_mod(process = ) will hand back its fitted values as though they were simulated.
  7. (Optional) expose an on/off or random-effect switch as a fit_mod() argument.
  8. Document it in inst/extdata/meta_data_names.xlsx and/or the fit_mod() roxygen.

Add a new reported quantity

  1. REPORT() it in ceattle.cpp (and ADREPORT() it if it needs a standard error). A REPORT inside a comment block is not reported, enumerate from a real fit, not from a grep, when checking what a model actually returns.
  2. Give it dimension names in R/6-rename_output.R if it is an array.
  3. 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.R reads the template and a fit and asserts the two agree, so a quantity added here and not there fails a test.
  4. If a quantity is only defined for some configurations, the per-recruit reference points are computed only under msmMode = 0, say so in its meaning, and gate it in report_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:

  1. Add a row to .rce_column_schema() with the column’s canonical name, type, default, any back-compat aliases, 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.
  2. Consume the column where it matters (data_check(), rearrange_data(), or the cpp) by its canonical name, aliases are already upgraded on entry.
  3. If it is a rename of an existing column, put the old spelling in aliases and pin it in test-schema-canonical.R.

A whole new data_list element (a new sheet / array, e.g. an index covariate):

  1. Add the DATA_* object in ceattle.cpp and consume it in the module.
  2. Add read and write support in R/0-read_write_excel_data.R. A new element with no write_data()/read_data() handling round-trips to nothing through the standard xlsx format; the feature is silently lossy. (This is exactly how index_cov was lost; it is the single most common data-input bug.)
  3. Add validation in R/1-data_check.R and, if the model needs it, a row in the data_requirements table (R/1-data_requirements_table.R).
  4. Add reshaping in R/5-rearrange_data.R if TMB needs a different layout.
  5. Document it in meta_data_names.xlsx and/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:

  1. Implement the case in the module .hpp (e.g. selectivity.hpp case N, predation.hpp smode N).
  2. Register the string ↔︎ integer mapping in the relevant map in R/0-switches.R (e.g. sel_map, suitMode_map); otherwise validate_switches() rejects the new name.
  3. 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)).
  4. Add bounds in R/4-build_parameter_bounds.R if the option introduces bounded parameters (see the suitMode %in% c(1:2) branch).
  5. Add a test named test-<area>-<topic>.R in tests/testthat/.
  6. 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 *_sim name. TMB never clears its report environment, so anything reported inside a SIMULATE block stays visible in obj$report() afterwards. A draw written into the observed-data object and reported under its name is readable as the data it replaced, which is why comp_sim, caal_sim and diet_sim are 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.rds and the golden-jnll regression test pin the objective-function value so unintended numerical changes surface.
  • test-tmb-makeadfun_smoke.R checks the model builds and evaluates.
  • test-parameter-recovery.R fits simulated data and checks estimates.
  • helpers-make-msm-data.R builds 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 NByageFixed sheet and selectivity via emp_sel with Selectivity = "Fixed", set estDynamics = 1, and run fit_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:

  • estimateMode 0 = “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_map skips 10, 12 and 14.
  • initMode 1 (Equilibrium) and 5 (OffsetEquilibrium) differ by one term.
  • msmMode accepts 3–9, then data_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.