Fit a Tiny Assessment Model (TAM)

Estimates abundance, recruitment, and mortality from catch-at-age and survey observations. Returns population trends, observation diagnostics, parameter estimates, and uncertainty tables. Use N_settings, F_settings, and M_settings to specify where unexplained biological variation is allowed. See prepare_tam() for options and tinyAM-model for the model equations.

Usage

fit_tam(
  data,
  years = NULL,
  ages = NULL,
  N_settings = list(process = "iid", init = "exp", mu_form = NULL, sd_form = ~1),
  rec_settings = list(process = "rw", init = "free", mu_form = NULL, sd = NULL, phi = NULL),
  F_settings = list(process = "rw", mu_form = NULL, sd_form = ~1, age_breaks = NULL, cor_form = ~1,
        mean_ages = NULL),
  M_settings = list(process = "off", mu_form = NULL, mu_supplied = ~I(0.2), sd_form = ~1,
        age_breaks = NULL, first_dev_year = NULL, mean_ages = NULL),
  catch_settings = list(sd_form = ~1, sd_supplied = NULL),
  index_settings = list(sd_form = ~1, sd_supplied = NULL, q_form = ~q_block, q_link = "log"),
  proj_settings = NULL,
  interval = 0.95,
  add_osa_res = FALSE,
  silent = FALSE,
  start_par = NULL,
  grad_tol = 0.01
)

# Optional projections:
# proj_settings = list(n_proj = 3, n_mean = 3, F_mult = 1, catch = NULL)

Arguments

data

A named list of tidy observation tables (e.g., catch, index, weight, maturity) and optional pop_covariates. A tam_data object or plain list can be used. The weight table may distinguish complete stock and catch weights with type; stock weights supply both purposes when catch weights are absent. See cod_obs for an example.

years

Consecutive historical model years (at least two). Inferred from observed data (non-projection) if NULL.

ages

Consecutive model ages (at least two). Inferred from observed data (non-projection) if NULL. If the ages in the data extend beyond max(ages), the "obs" column is summed for catch and index data within each year (and survey for indices). Older-age weight and maturity are retained for abundance-based aggregation inside nll_fun(), preserving biomass and mature biomass in the plus group.

N_settings

A list with elements:

process: "off" for deterministic cohort survival, "iid" for independent cohort residuals, "rw" for residuals that accumulate through time, or "ar1" for residuals correlated between years and ages. See tinyAM-model.
init: "exp" (default), "free", "latent", "mono" or "capped", independently of process. All use first-year recruitment log_r0 as the starting anchor. "exp" uses deterministic survivorship; "free" estimates fixed older-age log_n0 states; "latent" integrates them without an initial density. "mono" constrains abundance to decline from r0; "capped" bounds each older age by r0 without imposing age ordering. See Details.
mu_form: ordinary fixed formula for proportional cohort gains/losses, default NULL. Covariates come from pop_covariates at the destination year and older age, after survivors are combined (including the plus group). Recruitment and initialization are unaffected. These effects do not identify immigration separately from mortality or observation changes. RW residuals require independent temporal mean contrasts; constant directions cancel and are rejected. See tinyAM-model.
sd_form: ordinary formula for log SD of older-age survival residuals, default ~ 1. Use ~ factor(age) or a supplied age-sharing column in data$pop_covariates. Its design must be constant across years. Recruitment SD is controlled separately by rec_settings; initial abundance has no separate SD. Random-effect SD terms are not enabled.

rec_settings

A list controlling recruitment separately from cohort survival:

process: "off" for recruitment exactly following its mean formula, or "rw" (default), "iid", or "ar1" for annual log-recruitment residuals. Boundary recruitment states are controlled by init.
init: "free" (default) estimates first-year recruitment and additional pre-parent-SSB boundary states as fixed parameters; "latent" integrates them without an initial density or separate SD. Independent of process, including "off", and of older-age N_settings$init. Simulations condition on these states; joint uncertainty draws can vary them. See recruitment_formulas.
mu_form: one-sided formula for the expected log recruitment. NULL is allowed for RW; off/IID/AR1 require an explicit formula, such as ~ 1 for a constant mean. Use bh() or ricker() for a median stock-recruit curve with no residuals or IID/AR1 residuals. Ordinary covariates come from data$pop_covariates at the youngest modeled age. See recruitment_formulas.
sd: positive fixed recruitment SD, or NULL to estimate it. This is an innovation SD for RW/AR1 and a marginal SD for IID. Must be NULL when process is off.
phi: fixed AR1 correlation in [0, 1), or NULL to estimate it. Must be NULL for other processes. Estimated correlations are positive.

F_settings

A list with elements:

process: "off" makes historical F equal its explicit mean formula, without F latent states, residual SDs or correlations. "iid" gives independent departures from the mean log F, "rw" for departures that accumulate over time, or "ar1" for departures correlated between years and ages. "cor_rw" accumulates changes over time with AR1 correlation between annual increments at neighbouring states within each fleet. It estimates one signed age correlation per active correlation group; its SD is the marginal SD of annual log-F increments.
sd_form: ordinary age-based formula for log process SD, default ~ 1, evaluated on data$catch. Sharing an SD does not share latent F states. IID SDs describe residuals; RW/correlated-RW SDs describe annual increments; stationary AR1 retains its innovation-scale convention.
age_breaks: optional increasing integer interval boundaries for shared actual F. Intervals include their lower boundary and exclude their upper boundary; the final boundary must be max(ages) + 1. For ages 1:6, c(1, 2, 3, 4, 5, 7) shares F at ages 5 and 6. Ages below the first boundary have exactly zero F in all years, with no F state or process density. Positive catch observations at those ages are rejected. NULL retains separate F states at every modeled age. All mean components and process SDs must agree within a shared block; refine blocks for age-varying means. Process correlation uses ordered blocks, not replicated fish ages. With fleet labels, supply a common vector or a named list covering all fleets; NULL entries retain separate states by age within that fleet.
cor_form: categorical correlation groups constant within each fleet. Default ~ 1 shares correlations; ~ fleet estimates separate values. For AR1, list(age = ~ fleet, year = ~ 1) allows different sharing by dimension. A single formula applies to both dimensions. Sharing parameters never correlates innovations between fleets. Unused groups are omitted.
mu_form: an optional formula for mean-$F$ (coefficients estimated as log-scale parameters log_mu_f), evaluated on data$catch. Required for "off", "iid" and "ar1"; use ~ 1 for one estimated mean level or ~ factor(age) for age-specific means. For "rw", NULL is allowed. Mean columns constant through time cancel from RW increments and are removed with a warning. Remaining columns must be identifiable from historical year-to-year differences. Ordinary F mean designs use historical support only; projection rows carry the terminal historical design by age. No future-year coefficients are fitted. With IID residuals, iid(), rw() and ar1() can add random mean effects, e.g. ~ factor(age) + rw(year). Temporal mean terms must have disjoint support, for example one per fleet; overlapping temporal mean terms and structured means with temporal residual processes are not enabled. mono() and logistic() can instead describe a rising mean age curve with no residual process or with IID/AR1 residuals, e.g. ~ logistic(age). Ordinary terms supply baseline F; realized F can depart from the rising mean. Static curves are not identifiable with RW/correlated-RW residuals. See catchability_curves. double_logistic() supplies WHAM's unnormalized four-parameter dome. Catch covariates must be available even where observations are missing. They must agree across biological ages pooled into a modeled plus group. The catch response obs cannot enter an F formula. With fleet labels, each fleet has its own F states; total F is their sum for survival and Z. for_fleet() gives one fleet its own formula, e.g. ~ 0 + for_fleet("Trawl", ~ factor(year) + logistic(age)) + for_fleet("Gillnet", ~ factor(year) + factor(age)). Ordinary terms outside scopes share their coefficients. Each inner formula includes an intercept unless removed explicitly; avoid redundant global and scoped baselines. smooth_formulas adds penalized s() and age-year t2() relationships. Smooth means support "off", "iid" and "ar1" F processes; the latter combinations are experimental and require mean/process separation checks.
mean_ages: optional vector of ages to include in population weighted average F (F_bar) calculations. All ages used if absent.

M_settings

A list with elements:

process: "off" uses only supplied/mean M. Otherwise "iid", "rw", and "ar1" have the same meanings as for F, on selected years/age blocks.
sd_form: ordinary age-based formula for log process SD, default ~ 1. Evaluated on data$pop_covariates and required to be constant within each fitted M age block and across years. Use finer age_breaks for distinct state SDs.
mu_form: optional formula for mean-$M$ (applied on the log scale) built on data$pop_covariates, yielding coefficients mu_m. These enter the log-$M$ surface directly and may therefore be positive or negative; they intentionally avoid a log_ prefix to prevent automatic back-transformation when tidied. If provided together with mu_supplied, the intercept in mu_form is dropped (warning) so supplied levels act as fixed offsets. iid(), rw() and ar1() support random mean effects with IID residuals or the residual process off. Use ~ 0 + rw(year) to vary supplied M through time. All mean components must be constant within fitted M age blocks. With AR1 residuals, one age-defined IID random slope is also supported, e.g. ~ 0 + capelin + iid(M_age_group, by = capelin). The grouping column must be constant across years at each age and the by column numeric. This experimental combination retains shared actual M within blocks; see formula_effects for the model and statistical cautions.
mu_supplied: optional one-sided formula giving supplied (non-estimated) $M$, e.g. ~ I(0.2) or a column reference such as ~ M_assumption stored in the data$pop_covariates data.frame.
age_breaks: optional integer break points used by cut_ages() to define age_blocks sharing absolute latent M states across ages. Existing specifications retain their original interpretation. Alternatively, end at max(ages) + 1 to use explicit half-open intervals as for F: with ages 1:6, c(2, 5, 7) makes blocks 2–4 and 5–6. Ages outside M process blocks retain their existing mean M; they are not set to zero. When a narrower set of age_breaks than modeled ages is provided, process states are only estimated for ages within the specified range; ages outside this range are fixed to their mean or assumed levels. By default, blocks use cut_ages(ages[-1], unique(range(ages[-1]))), usually one block excluding the youngest age; two older ages form separate singleton blocks. This restricts how M can trade off against recruitment; it does not guarantee identifiability.
first_dev_year: one historical modeled year at which M process states begin. Defaults to the second year if NULL. Earlier M stays at its mean or supplied value; earlier years are not coupled to the first latent state. Ignored when the M process is off.
mean_ages: optional vector of ages to include in population weighted average M (M_bar) calculations. All ages used if absent.

catch_settings

A list with elements:

sd_form: formula for log observation SD for catch-at-age data. Ordinary terms and Gaussian formula_effects are supported, e.g. ~ iid(age).
sd_supplied: optional one-sided formula giving supplied SDs (on the natural scale of the log-observation residuals) for catch-at-age data. When provided, the intercept is removed from sd_form so supplied SDs act as offsets.

index_settings

A list with elements:

sd_form: formula for log observation SD for index-at-age data. Ordinary terms and Gaussian formula_effects are supported, e.g. ~ survey + rw(age, by = survey).
sd_supplied: optional one-sided formula giving supplied SDs (on the natural scale of the log-observation residuals) for index-at-age data. When provided, the intercept is removed from sd_form so supplied SDs act as offsets.
q_form: formula for catchability, evaluated on the index table. Ordinary ~ q_block allows any pattern across blocks; ~ mono(q_block) is non-decreasing across ordered blocks. See mono() for independent non-decreasing curves by survey. Structured effects such as ~ survey + ar1(year, by = survey) allow shared annual changes across ages. ~ survey + logistic(age, by = survey) fits rising selectivity curves. ~ q_key + power(N, by = q_power_group) allows abundance-dependent catchability using internally calculated abundance at survey time. Missing power-group labels mean p = 1, not missing observations. Powers require the log q link. See power() and formula_effects. double_logistic() allows dome-shaped selectivity and multiplies baseline q after its inverse link. Penalized s() and t2(age, year) terms instead modify the selected q-link predictor. See smooth_formulas for supported bases, grouping and penalties. for_survey() applies a formula to one named survey, e.g. ~ 0 + for_survey("RV", ~ logistic(age)) + for_survey("Sentinel", ~ q_block). Effects are built on each survey's own reporting ages, then contribute zero to the predictor elsewhere.
q_link: "log" (default) allows any positive catchability; "logit" restricts catchability to between zero and one. Use the logit link only when the survey index has an absolute abundance scale and this restriction is scientifically appropriate. Relative indices or surveys affected by herding may need the log link. The same formula and optional mono() terms work with either link. With "logit", coefficients and dq act on logit-q rather than log-q; covariate effects are no longer constant multipliers of q. This is a restriction, not a prior favouring a value near one. q can approach a boundary with large coefficients and SEs.

proj_settings

Optional list with elements:

n_proj: number of years to project (default NULL disables projections).
n_mean: number of recent historical years (default 3) used to average each weight type and maturity by age. Catch/index observations in projection years are missing. Other columns are copied from the last available year in each table.
F_mult: multiplier of terminal historical F (default 1, status quo). Can be a value of length 1 or length = n_proj. When it is a vector of length one, that multiplier is recycled across all projection years.
catch: total annual catch weight, in catch-number times mean-weight units (e.g. tonnes for thousand fish and kg). Supply one value recycled across years, or n_proj values. Specify either catch or F_mult, not both. catch = "observed" is restricted to fit_hindcast() evaluation; it is not a forecast input for an ordinary fit. See tinyAM-model for the catch solver and feasibility safeguards.

interval

Level in (0, 1) to use to generate confidence intervals, where applicable; default 0.95.

add_osa_res

Logical; add one-step-ahead residuals using "oneStepGaussianOffMode"? To select another method after fitting and retain its results, use add_osa_residuals().

silent

Logical; if TRUE, disables tracing information.

start_par

Optional fitted tam_fit or list of starting parameter values, matching the structure returned by make_par(). Useful for warm starts or retrospective runs where parameter estimates from a previous fit are reused to improve convergence and speed. Non-matching or missing entries are ignored. Fitted-object starts require the current stored parameter modes; regenerate older saved fits. Catchability coefficients are omitted if q_link changes.

grad_tol

Numeric tolerance passed to check_tam() when evaluating the fitted object's gradients. Defaults to 1e-2.

Value

A tam_fit list with components:

call: matched call.
dat: data list returned by prepare_tam().
obj: RTMB ADFun object.
opt: [stats::nlminb()] optimization result.
rep: list from obj$report().
sdrep: RTMB::sdreport() result.
fixed_par: fixed parameter estimates in a tidy format (see tidy_par()).
random_par: list of random parameter estimates in a tidy format (see tidy_par()).
is_converged: Did numerical checks pass? (see check_tam())
diagnostics: Numerical, structural and advisory checks from check_tam().
obs_pred: obs$catch and obs$index data augmented with predicted values, parameter estimates, and standardized residuals (see tidy_obs_pred()).
pop: A collection of population summaries in tidy format (see tidy_pop()).
rec: Recruitment predictions and residuals, with parent pairs and curve tables for stock-recruit models (see tidy_recruitment()).
formula_effects: When Gaussian q, F/M mean or observation SD formula terms are present, signed effect levels, RW increments and numeric-by contributions for reporting and dashboards.

Details

Builds data with prepare_tam() and initial parameters with make_par(). RTMB integrates random effects using the Laplace approximation; nlminb optimizes the remaining fixed effects. mono() increments dq and dmu_f have a zero lower bound, as do initial-abundance decrements d_n0. Their SEs remain on the increment scale, but symmetric Wald intervals are only local approximations at a boundary. Random-effect blocks are chosen from the active processes and Gaussian formula effects. Initial abundance and recruitment boundaries are integrated only when their respective init setting is "latent"; see tinyAM-model and recruitment_formulas. make_par() describes the parameter layout.

A warning is issued if the number of random effects exceeds 1.5 times the number of observed data points (rough identifiability check).

Parameters without likelihood support are removed or fixed: examples include constant mean directions canceled by RW increments, RW SDs with no increments, and AR1 correlations on singleton axes. The latent states themselves retain the conventions in tinyAM-model.

is_converged summarizes optimizer, gradient and uncertainty checks. Also inspect opt$convergence and opt$message; a passing numerical check does not establish that the data identify every biological component. Fitting also warns about sparse catchability effects and wide formula-effect uncertainty. These advisories are retained by check_tam() and do not change the likelihood, parameter constraints or numerical convergence criteria.

Omitted entries in a settings list use package defaults; update() instead retains the original fit's entries. Use explicit NULL to remove an optional formula or supplied offset. Unknown or duplicate entry names are errors. Supply spawning timing in data$maturity$spawn_time; see prepare_tam().

Examples

fit <- fit_tam(
  cod_obs,
  years = 1983:2024,
  ages = 2:14,
  N_settings = list(process = "iid", init = "exp"),
  F_settings = list(process = "rw", mu_form = NULL),
  M_settings = list(process = "off", mu_supplied = ~ I(0.3)),
  catch_settings = list(sd_form = ~ 1),
  index_settings = list(sd_form = ~ 1, q_form = ~ q_block),
  silent = TRUE
)
fit
#> Call:
#> fit_tam(data = cod_obs, years = 1983:2024, ages = 2:14, N_settings = list(process = "iid", 
#>     init = "exp"), F_settings = list(process = "rw", mu_form = NULL), 
#>     M_settings = list(process = "off", mu_supplied = ~I(0.3)), 
#>     catch_settings = list(sd_form = ~1), index_settings = list(sd_form = ~1, 
#>         q_form = ~q_block), silent = TRUE)
#> 
#> Objective: 1138.33
#> ✔ Numerical convergence: passed
#> ✔ Optimizer: 0 (relative convergence (4))
#> ✔ Max |gradient|: 0.00143 adjusted; 0.00143 raw (tolerance 0.01)
#> ✔ Positive-definite Hessian: yes
#> ✔ Model structure: no redundant model terms detected
#> 
#> Advisory findings (2):
#> ℹ Some age/series groups have standardized residual SDs outside 0.5 to 1.5;
#>   fitted-state shrinkage also affects spread.
#> ℹ Some age/series groups have absolute adjacent-year residual correlation above
#>   0.5.
#> 
#> Coefficients:
#>                        Estimate Lower 95% Upper 95% Std. Error SE scale
#> r0 [year=1983] [age=2] 485,911  390,954   603,932   0.111      log     
#> sd_r                   0.497    0.367     0.674     0.155      log     
#> sd_f                   0.493    0.430     0.566     0.070      log     
#> sd_n                   0.595    0.531     0.667     0.059      log     
#> sd_catch ((Intercept)) 0.381    0.301     0.484     0.121      log     
#> sd_index ((Intercept)) 0.428    0.372     0.493     0.071      log     
#> q ((Intercept))        1.751    1.245     2.463     0.174      log     
#> q (q_block3)           1.484    1.142     1.929     0.134      log     
#> q (q_block4)           1.477    1.083     2.013     0.158      log     
#> q (q_block5)           1.284    0.910     1.811     0.176      log     
#> q (q_block6-20)        1.002    0.710     1.416     0.176      log     
#> 
#> Terminal year (2024) estimates:
#>             Estimate Lower 95% Upper 95% Std. Error SE scale
#> abundance   259,302  187,890   357,854   0.164      log     
#> recruitment 39,314   18,416    83,927    0.387      log     
#> ssb         180,686  132,681   246,060   0.158      log     
#> F_bar       0.028    0.018     0.042     0.210      log     
#> M_bar       0.300    0.300     0.300     0.000      log     
#> 
#> SE scale: reported = same scale as estimate; logit = log-odds.
#> A small log-scale SE approximates CV (0.10 is about 10%); use CIs for large SEs.
head(fit$pop$ssb)
#>   year      est      lwr      upr        se se_scale is_proj   unit spawn_time
#> 1 1983 374324.2 300052.6 466980.2 0.1128411      log   FALSE tonnes          0
#> 2 1984 382659.2 308378.0 474832.9 0.1101134      log   FALSE tonnes          0
#> 3 1985 386083.0 306118.9 486935.2 0.1184096      log   FALSE tonnes          0
#> 4 1986 468682.0 364325.2 602930.5 0.1285113      log   FALSE tonnes          0
#> 5 1987 407237.1 319341.0 519325.8 0.1240512      log   FALSE tonnes          0
#> 6 1988 474914.6 367434.9 613833.7 0.1309152      log   FALSE tonnes          0
#>          timing
#> 1 Start of year
#> 2 Start of year
#> 3 Start of year
#> 4 Start of year
#> 5 Start of year
#> 6 Start of year

## Fit with projections (status quo F)
fit2 <- update(fit,
  proj_settings = list(n_proj = 3, n_mean = 3, F_mult = 1)
)

if (interactive()) {
  fits <- list("without_proj" = fit, "with_proj" = fit2)
  vis_tam(fits)
}

Usage

Arguments

Value

Details

See also

Examples