Add one-step-ahead & process residuals; rework residual/plot API (v4.6.0)#97
Open
grantdadams wants to merge 23 commits into
Open
Add one-step-ahead & process residuals; rework residual/plot API (v4.6.0)#97grantdadams wants to merge 23 commits into
grantdadams wants to merge 23 commits into
Conversation
First phase of one-step-ahead (OSA) residual support via TMB::oneStepPredict, covering the aggregate catch and index series. C++ (ceattle_v01_11.cpp): - Add DATA_VECTOR(obsvec) + DATA_VECTOR_INDICATOR(keep, obsvec), per-row index maps (catch_obsvec_idx / index_obsvec_idx), and an osa_mode flag (unused this phase; reserved for composition/diet). - Index (slot 0) and catch (slot 1) likelihoods read the observation from obsvec and gate it with keep. With keep == 1 (normal fitting) the joint nll is bit-identical to the pre-change model (verified against a baseline DLL: max |jnll_comp| diff = 0.0e+00). R: - build_osa_data() (R/0-osa_data.R) assembles obsvec, the obs_ctl position-map metadata, and the per-type obsvec index vectors in rearrange_data(), after the data.frame->matrix coercion so obs_ctl stays a data frame. - osa_residuals() + osa_diagnostics() (SDNR + simulated tail null intervals, Stewart & Monnahan 2025) + rceattle_osa S3 class with a plot method (Q-Q and residual-by-year, afscOSA styling). fit_mod() carries obs_ctl onto the fit. - residuals.Rceattle() gains type = "osa". Tests (tests-Likelihoods/test-osa-residuals.R): obsvec construction, obs_ctl-is-data-frame guard, finite-nll smoke, and an end-to-end osa_residuals() check on GOApollock.
Extends one-step-ahead residuals to age/length composition (comp) and conditional age-at-length (caal) data via a keep-aware conditional decomposition. C++: - New src/TMB/comp_osa.hpp: keep-aware dmultinom_osa / ddirmultinom_osa that decompose a composition into conditional binomials / beta-binomials (Trijoulet et al. 2023; ported from WHAM age_comp_osa.hpp), with osa_order() and osa_squeeze() helpers. With keep == 1 the conditional log-densities sum to the joint density. - comp (slot 2) and caal (slot 3) gain an osa_mode branch. osa_mode == 0 is the original weighted fitting code, byte-identical to before (verified: reported jnll_comp diff = 0.0e+00 vs a baseline DLL, with composition data). osa_mode == 1 reads bin counts from obsvec and uses the unweighted keep-gated density; AFSC MultinomialAFSC fleets are residualized under the full multinomial. - Adds comp_obsvec_idx / caal_obsvec_idx start-position maps. - Fixes a pre-existing CAAL bug: the unweighted slot-3 likelihood was written into slot 2 (comp); corrected to slot 3. Affects unweighted_jnll_comp only, not the fitted objective. R: - build_osa_data() appends comp and caal bin counts ((proportion + 1e-5) * N) to obsvec, flags each composition's final (sum-to-N) bin, and records the conditioning length for caal. - osa_residuals() supports types "comp"/"caal", rebuilds the model in osa_mode=1 at the fitted parameters, residualizes only non-last bins, defaults to all available types, propagates the caal conditioning length, and warns on non-finite residuals (poor convergence / sparse compositions). - plot.rceattle_osa() gains an afscOSA-style composition bubble panel. Tests: comp and caal obsvec construction, plus comp OSA end-to-end on GOApollock (837 residuals from 930 bins; deterministic; aggregate unaffected).
Extends one-step-ahead residuals to predator diet (stomach-content) composition for multispecies models with estimated suitability. C++ (ceattle_v01_11.cpp): - The diet likelihood (slot 18) gains an osa_mode branch. osa_mode == 0 is the original weighted per-stomach code, byte-identical to before (verified: reported jnll_comp diff = 0.0e+00 vs a baseline DLL, with the diet likelihood ACTIVE via suitMode = 4). osa_mode == 1 reads each stomach's prey counts from obsvec and uses the unweighted keep-gated dmultinom_osa / ddirmultinom_osa decomposition (shared with comp/caal). The "other prey" category is the last bin, fixed by the sum-to-1 constraint and dropped. - Adds the per-stomach diet_obsvec_idx start-position map (length n_stomach_obs). R: - build_osa_data() appends, per stomach, the prey counts plus the "other prey" residual category ((proportion + 1e-5) normalized, scaled to the stomach sample size), matching the TMB likelihood. Inclusion mirrors the template guard: predator suitMode > 0 and >= 1 prey item. suitMode is set from the fit_mod argument before rearrange_data, so the R and C++ inclusion agree. - osa_residuals() supports types = "diet" (opt-in; not in the default set since it applies only to multispecies models and is expensive). Tests: a direct build_osa_data() diet unit test (stomach grouping, "other prey" last bin, count values). End-to-end diet OSA verified on a multispecies model (suitMode = 4): 9300 diet positions -> 9000 finite residuals, deterministic.
Adds process_residuals(), a complement to the observation-based osa_residuals() that validates the model's random-effect process deviations (Nielsen and Berg 2014). Each deviation set carries a Gaussian process prior. The posterior mode is shrunk toward that prior, so a single draw is taken from the joint posterior of the deviations (from the joint precision when they are random effects, or the fixed-effect covariance when they are penalized fixed effects) and standardized by the process SD; under a correct process these are approximately iid N(0,1). A random draw is used, so results are stochastic -- a seed is exposed. A generic registry maps any supported deviation parameter to its estimated elements (via the map and arrayInd), prior mean/SD, and species/fleet/year/age labels. Supported processes: "recruitment" (rec_dev), "initial" (init_dev), and "catchability" (index_q_dev); process = "all" returns every supported process present. Selectivity and natural-mortality deviations use random-walk / 2D-AR1 priors and are not yet supported (they need full prior-covariance decorrelation). - residuals.Rceattle() gains type = "process". - Results use the rceattle_osa class, so osa_diagnostics() and the plot method apply directly. - Adds Matrix to Imports (Matrix::solve on the joint precision). Verified on GOApollock: 49 recruitment residuals (one per hindcast year, correctly mapped, deterministic, SDNR within the null interval); process = "all" returns recruitment + initial + catchability (154 finite residuals).
Adds a residual_type argument to the existing diagnostic plotters so OSA residuals are available through the familiar functions (Stewart & Monnahan 2025), not only via plot() on an rceattle_osa object: - plot_comp(fit, residual_type = "osa") draws the OSA Q-Q plot (with SDNR and tail-statistic annotation) and signed composition bubbles. - plot_indexresidual(fit, residual_type = "osa") draws the OSA Q-Q and residual-by-year panels for the aggregate catch/index series. The default residual_type = "pearson" keeps the legacy Pearson-residual and aggregate-fit plots unchanged; Stewart & Monnahan recommend retaining Pearson bubbles alongside OSA residuals. All four plot paths (aggregate, comp, index residual, process) verified to render.
Review pass over the OSA feature (correctness + readability): - process_residuals(): replace a silent fallback that could mislabel deviations (assigning the wrong species/year/age) when the estimated set is not a contiguous leading block -- e.g. partially time-varying catchability -- with a clear error, so residuals are never silently mislabeled. - build_osa_data(): factor the near-identical comp and caal row loops through a shared append_composition() helper, keeping each type's column extraction explicit. No behavior change. Docs: - NEWS.md: document osa_residuals(), osa_diagnostics(), process_residuals(), the residuals() type = "osa"/"process" options, the plot_comp()/plot_indexresidual() residual_type = "osa" options, and the CAAL unweighted-slot bug fix. - model-diagnostics vignette: a new OSA-residuals section with the Stewart & Monnahan (2025) 5-step diagnostic workflow, plus process residuals. - _pkgdown.yml: index the new diagnostic and plot functions.
- Bump version 4.5.0 -> 4.6.0 - Split NEWS: OSA residuals under a new 4.6.0 section; the 4.5.0 convergence diagnostics feature stays under 4.5.0 - Fix stale "later phase" comments left over from the incremental development (the obsvec DATA block in ceattle_v01_11.cpp, build_osa_data(), and plot.rceattle_osa()) so they describe the finished implementation (comp / caal / diet all supported). Comment-only; no change to results.
Consolidates projection_uncertainty with the other optimizer / reporting controls in fit_control(). Passing it directly to fit_mod() still works via the existing deprecation path (warns and forwards into fit_control), the same backward-compatible mechanism used for phase, getsd, and the other former fit_mod() control arguments.
Driven by simulation-testing performance and a residual-math review: - fit_control(osa = FALSE) (new, default off): skip building the composition OSA observation data during a fit. The aggregate index/catch obsvec the TMB template always reads is still built; only the costly comp/caal/diet metadata is gated, so the fitted objective is unchanged. Simulation drivers (run_mse, self_test, jitter, retrospective) inherit the fast default; set osa = TRUE to compute composition OSA residuals. build_osa_data() also assembles obsvec / obs_ctl in a single pass instead of growing them in a loop. - residuals() now follows stats::residuals.glm(): type selects the residual kind (response / pearson / osa / process) and a new source argument selects the data source(s); all applicable sources are returned by default. Pearson residuals are now available for index/catch (standardized by the realized observation log-SD). index/catch residuals are restricted to genuine observations (projection / NA / non-positive rows dropped). The legacy type = <source> form still works with a deprecation warning. - process_residuals() warns when catchability uses a correlated deviate prior (random walk or AR1), where standardizing by the marginal SD is only approximate; it is exact for the iid prior and for recruitment / initial. NEWS, model-diagnostics vignette, and tests updated; R CMD check clean.
plot.rceattle_osa() now draws up to two separate figures: - Aggregate (index / catch): a Q-Q panel only (no age/length bin -> no bubbles). - Composition (comp / caal): a Q-Q panel plus signed OSA-residual and Pearson-residual bubble panels, with age-based bins (age composition and conditional age-at-length) in the left column and length-based bins in the right column, each on its own bin axis. Panel headers are two rows -- fleet name (from fleet_control) over the data type. osa_residuals() now carries a fleet_name column and, for composition types, attaches the matching Pearson residuals (the "pearson" attribute) so plot() can show OSA and Pearson bubbles together. Adds a plot smoke test; NEWS and docs updated; R CMD check clean.
- plot.rceattle_osa() gains source / species / combine arguments: source and species subset the figures (mirroring residuals()), and combine = FALSE draws the age and length composition as separate figures instead of side-by-side columns. - osa_residuals(): rename the types argument to source (now also accepts "all" and "diet"), so residuals(), osa_residuals(), and plot() select data sources with one consistent vocabulary. type stays reserved for the residual kind (response / pearson / osa / process), per stats::residuals.glm(). - residuals.Rceattle() gains a species filter. - Update internal callers (plot_comp, plot_indexresidual, dev scripts), tests, the model-diagnostics vignette, NEWS, and docs. R CMD check clean.
…siduals - plot_comp() rewritten in ggplot2 (consistent with plot.rceattle_osa): composition Pearson bubbles plus observed-vs-fitted annual and aggregated figures. geom_area spans only the observed bins (no boundary overshoot), zero observed proportions are kept (only NA dropped), joint-sex is drawn on one age/length axis (females above, males below zero), and a species filter is added. .comp_resid_long() reuses residuals(type = "pearson"). - plot.rceattle_osa(): joint-sex composition bins are re-based onto a shared age/length axis and faceted by sex (OSA and Pearson bubbles); osa_residuals() now carries nages/nlengths for this. - osa_residuals(): split the oneStepPredict() call by observation type so 'discrete' is applied per type -- aggregate index/catch stay continuous, while composition can be discrete (randomized quantile via the generic method). - residuals(): add source = "diet" (predator/prey schema, returned on its own) and extract .pearson_proportion() shared by the comp/caal/diet branches; plot_diet_comp() now draws its diet residuals from residuals(source = "diet"). - Add diet residual tests (mock + end-to-end make_msm_test_data fit) and a plot_comp smoke test. NEWS and docs updated; R CMD check clean.
Reflects the ggplot2 rewrite: Pearson-residual bubbles plus observed (shaded area) vs predicted (line) composition, annual and aggregated, with joint-sex drawn on one bin axis.
…rison - tests: add a deterministic 'fitted objective unchanged' golden test (BS2017SS total jnll at the initial parameters, estimateMode = 3) that guards the OSA obsvec/keep machinery leaving normal fitting byte-identical. - R/dev: remove the now-redundant OSA verification scratch scripts (osa_phase4_check, osa_plot_check, osa_caal_construct) -- their checks are covered by the testthat suite. The jnll-identical regression (osa_phase1-3_check) and the test runner are kept. - tests/comparison/WHAM-growth-comparison.R: add a Rceattle-vs-WHAM OSA residual comparison (aggregate index/catch, plus a composition template) using the growth fork's make_osa_residuals(); the fork's age_comp_osa.hpp decomposition matches Rceattle's comp_osa.hpp.
…s, OSA speedup, diet tests comp_offset is configurable via fit_control(comp_offset=) and stored on data_list (default 1e-5; 0 reproduces a WHAM-style multinomial). It lives on data_list so internal re-fits inherit it, is applied identically in the fitting likelihood and the OSA observation vector, and is read by the TMB template as a DATA_SCALAR. switch_check() fills it silently (set_default gains msg = NULL); build_osa_data() and rearrange_data() read it from data_list. osa_residuals() gains parallel (default TRUE): the per-observation one-step-ahead loop runs via parallel::mclapply(), a near-linear speedup for random-effects models. It probes TMB::openmp() and falls back to serial when several TMB models are loaded (e.g. alongside WHAM, or the full test suite). The discrete randomized-quantile path stays serial for reproducibility. Likelihood / growth fixes validated against the WHAM growth fork: - recruitment / initial-abundance bias correction centered at -sigma^2/2 - plus-group growth SD pinned in the second growth builder - within-year plus-group mean length now pinned at Jan-1 for Richards too (was von-Bertalanffy-only), so the plus-group convention is consistent across growth families - est_phase = 0 intercept linkages map the base parameter off (fix-at-init) - init_dev added to the random vector only when actually estimated Tests and docs: golden jnll = 1537036.2876293703 (default comp_offset); new diet OSA tests (decomposition invariant + end-to-end); plus-group within-year test now asserts pinning for VB and Richards; recruitment test uses the -sigma^2/2 mean; composition-likelihood test uses the renormalized multinomial; WHAM comparison scripts pinned to comp_offset = 0; NEWS, model-diagnostics vignette, Rd docs.
…s, OSA speedup, diet tests comp_offset is configurable via fit_control(comp_offset=) and stored on data_list (default 1e-5; 0 reproduces a WHAM-style multinomial). It lives on data_list so internal re-fits inherit it, is applied identically in the fitting likelihood and the OSA observation vector, and is read by the TMB template as a DATA_SCALAR. switch_check() fills it silently (set_default gains msg = NULL); build_osa_data() and rearrange_data() read it from data_list. osa_residuals() gains parallel (default TRUE): the per-observation one-step-ahead loop runs via parallel::mclapply(), a near-linear speedup for random-effects models. It probes TMB::openmp() and falls back to serial when several TMB models are loaded (e.g. alongside WHAM, or the full test suite). The discrete randomized-quantile path stays serial for reproducibility. Likelihood / growth fixes validated against the WHAM growth fork: - recruitment / initial-abundance bias correction centered at -sigma^2/2 - plus-group growth SD pinned in the second growth builder - within-year plus-group mean length now pinned at Jan-1 for Richards too (was von-Bertalanffy-only), so the plus-group convention is consistent across growth families - est_phase = 0 intercept linkages map the base parameter off (fix-at-init) - init_dev added to the random vector only when actually estimated Tests and docs: golden jnll = 1537036.2876293703 (default comp_offset); new diet OSA tests (decomposition invariant + end-to-end); plus-group within-year test now asserts pinning for VB and Richards; recruitment test uses the -sigma^2/2 mean; composition-likelihood test uses the renormalized multinomial; WHAM comparison scripts pinned to comp_offset = 0; NEWS, model-diagnostics vignette, Rd docs.
`residuals()` now follows stats::residuals.glm() strictly: `type` selects the
residual kind only ("response"/"pearson"/"osa"/"process"); data sources are
selected with `source`. Passing a source name via `type` previously rerouted it
to `source` with a one-time warning; that back-compat path and its warning are
removed (match.arg now errors on an unknown `type`). Test and NEWS updated.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Releases v4.6.0. Adds one-step-ahead (OSA) and process residuals for model validation, reworks the residual and composition-plotting API into a consistent,
glm-style interface, and makes OSA composition diagnostics fast enough to leave on during simulation testing. This release also adds likelihood and growth corrections that intentionally change the fitted objective for affected models.The new OSA machinery itself (the
obsvec/keep/osa_modedata path) leaves normal fitting byte-identical; the objective changes come only from the explicit corrections below.22 commits; 61 files changed.
New features — model validation residuals
osa_residuals()— one-step-ahead residuals (Thygesen et al. 2017; Trijoulet et al. 2023) for all fitted observation types: survey index, fishery catch, age/length composition, conditional age-at-length, and predator diet. Composition uses the conditional binomial / beta-binomial decomposition (a direct port of WHAM'sage_comp_osa.hpp); theoneStepPredict()call is split by observation type so continuous (lognormal) index/catch and (optionally discrete) composition are residualized correctly. Aparallelargument (defaultTRUE) runs the per-observation loop viaparallel::mclapply— a near-linear speedup for random-effects models.process_residuals()— SAMprocres-style residuals on the model's random-effect deviations (recruitment, initial abundance, catchability). Exact for iid priors; warns when a catchability deviate prior is correlated (random-walk/AR1).osa_diagnostics()— SDNR and lower/upper tail statistics with their standard-normal null intervals (Stewart & Monnahan 2025; Francis 2014).plot()forrceattle_osa— Q-Q (with SDNR/tail annotation) plus signed OSA and Pearson bubble panels, styled after NOAA-AFSCafscOSA. Aggregate and composition are separate figures; ages and lengths are separate columns; joint-sex bins share one axis (females above, males below); panel headers use fleet names.Residual / plotting API (glm-style)
residuals()followsstats::residuals.glm():typeselects the residual kind ("response","pearson","osa","process") andsourceselects the data source ("index"/"catch"/"comp"/"caal"/"diet"/"all"), plus aspeciesfilter. Pearson residuals are now available for index/catch and for diet. .plot_comp()re-implemented in ggplot2: composition Pearson bubbles plus observed-vs-fitted annual and aggregated figures. Fixes a histogram bug where the polygon extended past the first/last bin; zero observed proportions are retained (onlyNAdropped); joint-sex is drawn on one axis.residual_type = "osa"routes to the OSA diagnostics.plot_diet_comp()now draws its diet residuals fromresiduals(source = "diet")— a single source of truth for residual math.fit_control / performance
fit_control(osa = FALSE)(new, default off): a fit skips assembling the composition OSA observation data, so repeated refits stay fast (run_mse(),self_test(),jitter(),retrospective()). Setosa = TRUEto compute composition OSA residuals (~150 ms → ~5 ms per fit TRUE -> FALSE).fit_control(comp_offset = NULL)(new): the small proportion offset added to age/length-composition and CAAL bins before the multinomial / Dirichlet-multinomial.NULLinheritsdata_list$comp_offset(default1e-5). Stored ondata_listso internal re-fits inherit it, and applied identically in the fitting likelihood and the OSA observation vector.projection_uncertaintymoved fromfit_mod()intofit_control()(deprecation-forwarded).Likelihood & growth corrections (change results for affected models)
These bring Rceattle into exact agreement with the WHAM growth fork at the MLE and fix statistically incorrect behavior.
−σ²/2(was+σ²/2), soE[R] = R0(mean-unbiased). Woops!Comp_loglike = "Multinomial") is now evaluated through the OSA conditional-binomial decomposition, which normalizes the predicted proportions to sum to 1 (matching WHAM); the previous plaindmultinomused un-normalized proportions after the offset. OJO We may want to revisit this....exp(sd_Linf), and the within-year plus-group mean length is pinned at its Jan-1 value for Richards as well as von Bertalanffy — so the plus-group convention is consistent across growth families.Bug fixes
unweighted_jnll_compinstead of the CAAL slot (diagnostic only; not the fitted objective).init_devwas added to the random-effect vector even when fully mapped off (initMode = "Equilibrium"), producing an NA/NaN gradient underrandom_rec; it is now treated as random only when at least one element is estimated.est_phase = 0intercept linkages now map the base parameter off as well (honoring the documented "fix at init" contract), so a fixed-at-init linkage truly fixes the parameter.Backward compatibility & verification
obsvec/keep/osa_modemachinery feedsTMB::oneStepPredict()while leaving normal fitting (osa_mode == 0) byte-identical — confirmed by a jnll-identical regression against the pre-change template.GiancarloMCorrea/wham, refgrowth) to ~1e-10, which in turn enables exact OSA cross-validation (aggregate residuals identical under fixed effects; order-dependent once random effects couple the observations).fishfollower/SAM)procres.make_msm_test_data) models.Testing
Full
testthatsuite passes, with new tests for OSA construction, end-to-end OSA / process / diet residuals (including a diet decomposition-invariant check), the residual-kind API, the composition/diet plots, and growth plus-group pinning for VB and Richards. Golden-jnll regression updated to the corrected objective (1537036.2876293703on BS2017SS). (Recommend re-runningR CMD check --as-cranon the final branch before merge.)Known follow-ups (documented, not blocking)
ADREPORTof standardized residuals.eval = FALSEOSA chunks in the model-diagnostics vignette once a converged CAAL example is available.