3. Model diagnostics

Naming the parameter, not just the block

TMB names every element of obj$par after its parameter block, so on their own max_gradient, estimability, the two bounds checks and hessian_conditioning could say sel_coff_dev was flagged but not which of its 392 elements. Each now reports the fleet, sex, age, bin and year instead, from [parameter_index()], and follows the block name with the quantity it estimates, from [parameter_dictionary()]:

  [WARN] max_gradient     Maximum absolute marginal gradient = 0.00222 on log_M1 (M1): Pollock, age 1.
  [WARN] hessian_conditioning Condition number = 2.4e+06; loads on: log_F (F): 67% + rec_dev (recruitment deviations): 26%.
  log_F (F)                        Hake_fishery, 38 years in 1980-2021  (38 of 44 parameters; 67% of the direction)
  rec_dev (recruitment deviations) 24 years in 1981-2013  (24 of 44 parameters; 26% of the direction)

A scattered set of years or ages is counted against its span (“38 years in 1980-2021”); a consecutive run prints as a range. The count in parentheses is how many of the block’s parameters carry 90% of its loading.

The checks do not all count positions in the same vector. max_gradient, estimability and the two bounds checks are captured from the hindcast optimization; hessian_conditioning reads the sdreport, which under an estimating HCR belongs to the projection fit and holds only log_Ftarget / log_Flimit. Each check therefore verifies that the index describes the vector it is about to label, and prints block names alone when it does not: a coordinate naming the wrong parameter is worse than no coordinate.

parameter_index() itself describes fit$obj, which under any HCR but "NoFishing" is that projection object. Pass a fit run with estimateMode = "Hindcast" to browse the hindcast parameters directly.

Composition data

plot_comp() draws composition Pearson-residual bubbles and overlays the observed (shaded area) and predicted (line) age or length compositions for every fleet, both annually and aggregated across years. Joint-sex data are drawn on a single bin axis with females above and males below zero.

plot_comp(model_1)

Survey indices

plot_index() shows observed vs. predicted survey biomass; pass log = TRUE for the log scale. The predicted line is the median of the fitted observation distribution, which is where the lognormal likelihood is centred, so about half the observations should fall below it. plot_index_residual() plots the Pearson residual – observed minus predicted, standardized by the sd the row was fitted with – a useful first check for time trends or heteroscedasticity. A positive residual means the survey saw more than the model predicted. These return ggplot objects.

plot_index(model_1)
plot_index(model_1, log = TRUE)
plot_index_residual(model_1)

plot_catchability() draws the other half of the same fit: the catchability the model scaled each survey by, faceted by fleet. It reads the realized q, so one figure covers an estimated or fixed mean, annual deviations under Time_varying_q, an environmental regression, a q linkage such as rw(1 | Year), and the closed-form "Analytical" value. A fleet with a time-invariant q draws a flat line. Hindcast years only, the model does not project catchability.

plot_catchability(model_1)
plot_catchability(model_1, log = TRUE)

Catch

plot_catch(model_1)

Choosing the one-step-ahead method

method selects how the conditional distribution is read. The default, "oneStepGaussianOffMode" (WHAM’s and SAM’s), approximates it as Gaussian: the observation is treated as a free variable, the mode of the conditional density is found in that coordinate, and the observation is standardized against it.

method = "cdf" instead asks the model for the conditional CDF and returns qnorm(F(x)). For a composition bin the conditional is a binomial, written through the regularized incomplete beta so it is defined at the fractional counts composition data hold; for the aggregate series it is pnorm. On a fixed-effect model nothing is approximated; with random effects see below.

The two do not agree on composition data, and the disagreement is measurable. Simulating observations from a fitted model and residualizing them at the parameters that generated them makes the right answer exactly standard normal, so the methods can be scored (tools/verify/verify-osa-cdf.R, BS2017SS, 20 replicates of 4538 composition bins):

The null standard error on the autocorrelation is 0.015, so only the last row is consistent with the residuals being iid standard normal. A composition bin holds a count, so its conditional CDF is a step function and qnorm(F(x)) inherits the step; the randomized quantile residual qnorm(F(x) − U f(x)) (Dunn and Smyth 1996, as Trijoulet et al. 2023 prescribe) removes it. That is why discrete defaults to TRUE under "cdf", and why "cdf" with discrete = FALSE is the worst of the three, not the best. On the aggregate index and catch series, which really are Gaussian, all three agree to 5e-5 wherever the 8.04 ceiling below does not bind.

The step it randomizes over is the model’s, not the data’s: the counts the likelihood reads are proportion × Sample_size, so on a fleet whose sample size has been tuned they are not whole fish. That is the right choice: an OSA residual asks whether the observation is consistent with the model as specified, but it is worth knowing that a heavily down-weighted fleet is being residualized against a multinomial it was never sampled from.

osa_cdf <- osa_residuals(model_1, source = c("index", "comp"), method = "cdf")
osa_diagnostics(osa_cdf)

Three things to know before switching:

• predicted is NA on every row. The method never forms a conditional mean, which is the point, under the Gaussian methods that mean is an expected bin count that can come back negative on a sparse bin, and the residual on those rows is then biased positive. Fitted values are in model_1$quantities and in residuals(model_1, type = "pearson"). (sd is NA under the default method as well; only method = "oneStepGaussian" returns one.)

• A Dirichlet-multinomial fleet cannot use it. Its conditional is a beta-binomial, which has no closed-form CDF and cannot be summed at a fractional count. Those fleets are residualized with "oneStepGaussianOffMode", announced in a message, and recorded in the object’s method attribute. Their predicted is blanked too, so the column keeps one meaning across the object.

• |residual| is censored at 8.04, in both directions. The upper end is forced by double precision, F is recovered as 1/(1 + exp(.)), which saturates at the last double below one. The lower end is not: that same expression takes a small F down to about 1e-308. It is censored to match anyway, because an asymmetric ceiling would show as a long left tail against a wall on the right, which is what skewness in the residuals looks like. The cost is real, osa_diagnostics() computes SDNR and the tail statistics on the censored values, so it bites hardest on a short series where one observation holds the statistic. osa_residuals() warns when any residual sits there. Which method to reach for then is measured, not obvious (tools/verify/verify-osa-cdf-accuracy.R, a 12-year survey with one observation multiplied by 200):

  method	min	max	non-finite	SDNR
  oneStepGaussianOffMode	−11.77	0.19	1	NA
  oneStepGaussian	−11.77	38.98	0	12.89
  oneStepGeneric	−3.95	3.33	0	2.21
  cdf	−8.04	8.04	0	4.59

  Only "oneStepGaussian" reports the uncensored value here: the package default returns NaN on that row, so its SDNR is unusable rather than merely large, and "oneStepGeneric" compresses the outlier harder than the ceiling does. The default’s failure is magnitude-dependent (at a ×20 outlier it is finite and matches "oneStepGaussian") so reach for that method specifically, and on the fleet in question: it costs an nlminb per observation.

With random effects, the choice reverses for the aggregate series

"cdf" is exact only on a fixed-effect model. With random effects oneStepPredict() integrates the CDF over the latent states by Laplace, and that integrand is a Gaussian times a sigmoid rather than a density, where, for a Gaussian observation, the Gaussian methods integrate a density and are exact. Measured against the exact Kalman innovations of a linear-Gaussian state space model, as the latent state becomes more informative relative to the observation (reproduced by tools/verify/verify-osa-cdf-accuracy.R):

That is a result about a linear-Gaussian model, and it does not transfer wholesale. The Gaussian methods are exact when the one-step-ahead predictive is Gaussian, which needs the model to be linear in the random effects, and Rceattle’s index and catch are exp() of cumulated log recruitment deviations pushed through the population dynamics. fullGaussian and oneStepGaussian cannot differ for a Gaussian conditional, so their disagreement measures the departure: on a 17-deviation fixture they differ by 0.091 on both index and catch, while "cdf" differs from oneStepGaussian by 0.017 on index. No method is exact for index or catch under random effects.

"ecov" is the exception, and there a Gaussian method is right. Its conditional (a Gaussian measurement of an AR1 latent, every other data term unconditional) genuinely is linear-Gaussian: the two Gaussian methods agree to 4e-14 on the QAR1 fixture, while "cdf" sits 0.139 away, about a quarter of the residual sd. It does not reverse for compositions: their conditional is discrete and skewed, which is what the Gaussian methods get wrong, and by much more. Simulating from a 22-random-effect model with the recruitment deviations redrawn and residualizing at the generating parameters (1680 residuals), oneStepGaussianOffMode gives mean +0.513, sd 0.404 and KS rejection in every one of 120 replicates, against +0.006, sd 1.002 and 6 of 120, the nominal 5%, for "cdf" with discrete = TRUE.

What limits "cdf" on compositions is the depth of the conditioning rather than the random effects themselves. On BS2017SS with random_rec = TRUE (159 random effects, 4538 composition bins) it returns 1879 non-finite residuals, against 0 for "oneStepGaussianOffMode" on the same fit, and takes hours; the same rows residualized on their own come back clean, and the 22-random-effect model above is residualized correctly. So try it, read the warning it issues, and where the residuals come back non-finite in bulk fall back to a Gaussian method for that source, remembering that its composition residuals are the under-dispersed ones in the comparison above.

Both halves of that split are what the SAM authors do, in two packages. stockassessment::residuals.sam() names no method, so it takes oneStepPredict()’s "oneStepGaussianOffMode" default with discrete = FALSE, on a model that is random effects throughout, whose observations are Gaussian on the log scale. Their composition package compResidual::resMulti(), the implementation Trijoulet et al. (2023) cite, instead hardcodes method = "cdf", discrete = TRUE. Rceattle fits both kinds of data, so it makes the same choice per observation type inside one call rather than across two packages.

One caveat that applies to every method: Estimate_index_sd = "Analytical", Estimate_catch_sd = "Analytical" and the analytical Catchability forms concentrate their parameter out using all the observations, including the one being residualized, so the conditioning is not strictly one-step-ahead and the residuals are approximate – by an amount that shrinks as the series lengthens.

Fleets with tail accumulation

Where a fleet folds its composition tails into a boundary bin (Comp_accum_young / Comp_accum_old), both the OSA and the Pearson residuals describe the bins the likelihood actually fit, not the raw bins in comp_data. A fleet on 10 ages with Comp_accum_young = 3 reports bins 3 to 10: bin 3 holds ages 1 to 3 combined, and the accumulated column (Accumulated from residuals()) marks it, since it stands for a range rather than the single age it is named for. Expect fewer residual bins than raw age bins for such a fleet.

The two differ in one place. The one-step-ahead decomposition drops each composition’s last bin, which is fixed by the sum-to-sample-size constraint. Under an old-tail accumulation that dropped bin is the upper accumulated one, so the fleet’s upper boundary bin appears in the Pearson residuals and not in the OSA residuals.

Following Stewart and Monnahan (2025), a practical workflow is:

1. Inspect the aggregate fit across years (plot_comp(model_1)): systematic lack of fit here points to structural misspecification (e.g. the wrong selectivity shape). Counts are pooled across the fitted hindcast years, and the interval holds 95% of the data the model predicts, not a confidence interval on the mean, so observations outside it are the ones to explain. The panel is annotated with the input sample size, the effective sample size the likelihood assumed (ESS (likelihood)), and the one this fit’s residuals imply (ESS (McAllister-Ianelli)). The gap between those two is the reweighting decision: they agree when the fleet is already weighted the way its residuals suggest.
2. Compare the SDNR to its null interval (from osa_diagnostics()): a value outside the interval flags over- or under-dispersion, often a data-weighting or effective-sample-size issue.
3. Compare the lower/upper tail statistics to their null intervals: these localize departures to the tails of the residual distribution. They are reported as order statistics, lower_r / upper_r name which one, and lower_p / upper_p the nominal probability it sits at, which is near but not exactly 2.5%, against the exact Beta(r, n - r + 1) null. The statistic and its interval must come from the same estimator, so osa_diagnostics() no longer simulates and its nsim / seed arguments are ignored. The residuals themselves remain randomized-quantile residuals: pass a seed to osa_residuals() to control those.
4. Inspect Pearson bubble plots (plot_comp(model_1), the default): patterns in the sign of residuals across age/year can identify where lack of fit occurs, Stewart and Monnahan recommend keeping these alongside OSA residuals.
5. Revise the model or data weighting where misspecification is found, using this together with the other diagnostics rather than as an accept/reject test.

When the diagnosis is a weighting problem rather than a structural one, reweight_comps() runs the McAllister-Ianelli (1997) tuning loop: it refits the model the way it was fitted, updating each fleet’s Comp_weights until the implied weights settle.

tuned <- reweight_comps(model_1, n_iter = 10, tol = 0.01)
tuned$reweight$history    # weight per fleet, per iteration

Restrict it to particular fleets with fleets =, taking either Fleet_code values or names matching fleet_control$Fleet_name, give ids or names, not a mix, since R coerces c(1, "Survey") to character.

tuned <- reweight_comps(model_1, fleets = c("Bottom_trawl", "Fishery"))

Only multinomial fleets are tuned. A Dirichlet-multinomial fleet estimates its own weight inside the likelihood, so tuning it externally would compete with that estimate; those fleets are named and skipped.

The alternatives are Francis (2011) reweighting, applied by hand, and estimating the weights directly with a Dirichlet-multinomial likelihood (Comp_distribution = "DirichletMultinomial"), which avoids the iteration entirely, note that under that likelihood Comp_weights is read as the log of the starting weight, so the model default of 1 is a starting weight of e, and a weight of 1 is written as 0.

Because OSA residuals for discrete composition data are randomized quantile residuals, they are stochastic; osa_residuals() takes a seed for reproducibility.

Process residuals

process_residuals() is the complementary check on the process model. It standardizes the model’s random-effect deviations (recruitment, initial abundance, catchability) against their process prior, drawing once from the joint posterior (SAM’s procres; Nielsen and Berg 2014). Under a correct process model these too are approximately iid standard normal.

pr <- process_residuals(model_1, process = "recruitment")
osa_diagnostics(pr)
plot(pr)

# residuals() exposes them too:
residuals(model_1, type = "process")

Redrawing process error

process widens the question. The default (FALSE) asks whether the estimator recovers its own parameters from new observations; naming a process redraws that process too, and asks whether it recovers something it was not shown.

# Recruitment deviations redrawn as well as the observations
sims <- self_test(model_1, nsim = 10, process = "recruitment")

# "dynamics" covers recruitment, M and growth; "observation" covers
# catchability and selectivity; TRUE covers everything.
sims <- self_test(model_1, nsim = 10, process = "dynamics")

Recruitment covers the initial age structure as well as the hindcast deviations, since the initial numbers-at-age are recruitment from before styr. Natural mortality covers the M1_re random effects, and growth, catchability and selectivity are redrawn wherever they are written as a random linkage (see vignette("environmental-linkages-and-priors")). A process the model gives no distribution to (asking for M when M1_re = 0, say) is left alone and reported in a warning rather than silently returned unchanged.

Recruitment has one such case worth knowing about. Under srr_fun = 0 with srr_pred_fun > 0 (the AMAK/Ianelli configuration, where recruitment is estimated annually and the stock-recruit curve is fitted as a penalty on those same deviations) the deviations are scored by two densities at once.

Two terms on one quantity do not compose into a distribution to draw from. Drawing at R_sd from the first alone ignores the second, and the second is not a density on the deviations alone either: it couples them to the stock-recruit parameters, which are estimated too. Either way the self-test would measure recovery against a process the model does not assume, so the deviations are not redrawn and sim_mod() says why. A random linkage on a recruitment parameter is a separate latent and is still drawn. Set srr_fun to the stock-recruit function itself to redraw the deviations.

The fitted deviations are no longer what generated the data, so the values that did come back attached to each replicate, each with a mask of the cells the draw touched:

truth <- attr(sims, "process_sim")[["Sim_1"]]
rec   <- truth$rec_dev[truth$rec_dev_drawn]

Compare estimates against those, not against the operating model. The mask matters: the process draws cover the hindcast only, and beta_linkage_re spans every random linkage whether or not its process was asked for, so a correlation taken over the whole array scores fitted values as though they had been simulated. For the same reason compare_sim(), which compares against the operating model, is not valid for process-drawn replicates.

One caveat before reading a redrawn M as validation. On BS2017SS with a year-varying M random effect, refits recover the simulated deviations (correlation around 0.5, positive in every replicate) but recover their standard deviation poorly, roughly 70% low, and collapsing to zero in a third of replicates. That is the usual behaviour of a variance component estimated from little information, not a fault in the simulation, and a year-varying M is weakly informed because it trades off against selectivity and catchability. Read a self-test that redraws M as a test of deviation recovery rather than of variance recovery.

When the returned list is shorter than nsim, debug = TRUE returns every simulation so the failures can be read directly rather than inferred:

sims <- self_test(model_1, nsim = 10, debug = TRUE)

failed <- names(which(!attr(sims, "converged")))
sims[[failed[1]]]$convergence          # which check failed, and on which parameter

Moving a set of cells together

slots cross, so k slots give a k-dimensional profile. That is wrong for a parameter supplied as a vector, an empirical age-based M profiled over 13 values across ten ages would be 13^10 fits. joint puts every named cell on one grid instead: "multiply" scales each cell’s current value, preserving the shape of the schedule (1 is the fitted model), and "add" shifts it, preserving the increments (0 is the fitted model).

sp <- 1
prof_M_all <- profile(
  fitted = model_1,
  param  = "M1",
  slots  = lapply(seq_len(model_1$data_list$nages[sp]),
                  function(a) c(sp, 1, a)),   # every age bin of species 1
  values = list(seq(0.7, 1.3, by = 0.05)),
  joint  = "multiply"
)

plot_profile(prof_M_all)

Catchability

param = "q" profiles a fleet’s base catchability, and slots takes a fleet name. Fleets sharing a Catchability_index share one q, so naming any member profiles the whole group; the others are added with a message.

prof_q <- profile(
  fitted = model_1,
  param  = "q",
  slots  = list("Bottom_trawl"),
  values = list(seq(0.5, 2.0, by = 0.1))
)

profile() returns Rceattle_list (one fit per grid row, NULL where the fit failed), grid (the user-scale value grid, a multiplier or offset under joint), nll (joint NLL aligned with grid, NA for non-converged fits), alias (the natural-scale name you profiled under, if you used one), joint, and echoes of param and slots. To cross-profile across multiple species, supply one slot per species (e.g. slots = list(1, 2, 3) with param = "sigmaR"); to cross-profile M1 across sex, supply one slot per sex as in the third example above. cores behaves the same way as in the other diagnostic functions.

Where the data disagree

The total says how well the data determine the parameter, not whether the data sources agree about it, a smooth quadratic total can be two surveys pulling in opposite directions. plot_profile() draws the components separately, in the layout r4ss::SSplotProfile() uses.

prof_M <- profile(
  fitted = model_1,
  param  = "M1",
  slots  = list(c(1, 1, 1)),
  values = list(seq(0.20, 0.45, by = 0.025))
)

plot_profile(prof_M, xlab = "M at age 1")

Read it by where each curve bottoms out, not by how deep it is: every series is re-zeroed at its own minimum and a point marks it, so the spread of those points along the x axis is the disagreement. The heavy black line is the total.

When one component dwarfs the rest the others flatten onto the axis and where they prefer the parameter cannot be read. relative = "scaled" puts every curve on 0 to 1 so the minima can be compared directly; it discards magnitude, so raise minfraction alongside it.

profile_components() returns the same numbers, one row per grid point per component, each labelled with the fleet or species it belongs to.

comps <- profile_components(prof_M)
head(comps)

Under random_rec = TRUE the total is the Laplace-approximated marginal likelihood while the components are the inner joint negative log-likelihood, so they will not sum; compare the shapes, not the sums. See ?plot_profile for relative, minfraction, model overlays, and which components weighted = FALSE changes.