---
title: "Jack mackerel base model: MCMC uncertainty"
subtitle: "ADNUTS default run for the single-stock model h1_1.06"
author: "Jim Ianelli"
date: 2026-09-10
lang: en
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
    lightbox: true
    code-fold: true
    theme: cosmo
execute:
  echo: false
  warning: true
  message: false
---

[Back to the jack mackerel wiki](../../assessments.html) · [Current SC14 assessment](../../assessment/SC14.html)

```{r}
#| label: read-saved-results
library(jsonlite)
library(knitr)
manifest <- fromJSON("run-manifest.json", simplifyVector = TRUE)
diag <- fromJSON("diagnostic-summary.json", simplifyVector = TRUE)
parameters <- read.csv("parameter-summary.csv", check.names = FALSE)
derived <- read.csv("derived-summary.csv", check.names = FALSE)
derived_diag <- fromJSON("derived-status.json", simplifyVector = TRUE)
terminal <- derived[derived$year == 2026 & derived$type %in% c("SSB", "Recruits"), ]
chains <- read.csv("chain-diagnostics.csv", check.names = FALSE)
green <- read.csv("kobe-green-probability-2007-2026.csv")
stopifnot(nrow(green)==20L, identical(green$year,2007:2026), all(green$retained_draws==3000L),
          all(abs(green$probability_green-green$green_draws/green$retained_draws)<1e-12))
fmt <- function(x, digits=2) ifelse(is.finite(x), formatC(x, format="f", digits=digits, big.mark=","), "Not available")
parameter_rows <- parameters[parameters$is_parameter, ]
passes <- identical(diag$screening_status, "passes reported screening") && isTRUE(derived_diag$selected_quantities_pass_rhat_ess_screen)
stopifnot(identical(as.integer(unlist(manifest$actual_dimensions)), c(2000L,3L,1594L)),
          diag$dimensions$retained_total == 3000L,
          diag$dimensions$active_parameters == 1593L,
          nrow(parameter_rows) == 1593L,
          sum(chains$retained_draws) == 3000L)
```

## What the run shows

The default ADNUTS run sampled **1,593 active parameters** in the single-stock base model **h1_1.06**. Three chains each completed 2,000 iterations; the first 1,000 iterations of each chain were warmup. The summaries below use the remaining **3,000 draws**.

```{r}
#| results: asis
#| label: interpretation-status
if (passes) {
  cat('::: {.callout-note title="Sampling checks passed"}\nThe run passes the parameter and sampler screening checks reported below. Its uncertainty estimates are conditional on this model, its data, its fixed parameter values, and its likelihood and penalty terms. Model adequacy also requires evaluation of its assumptions and fit to the observations.\n:::\n')
} else {
  cat('::: {.callout-warning title="The default run needs further sampling checks"}\nThe run completed with zero retained divergences and zero tree-depth-limit events. Some parameters have slow mixing and fail the screening criteria below. The figures and intervals provide an exploratory description of these retained draws. Further sampling and validation are needed before using these uncertainty estimates for management advice.\n:::\n')
}
```

The largest finite parameter R-hat is **`r fmt(diag$parameter_checks$maximum_finite_rhat, 3)`**. The smallest finite bulk and tail effective sample sizes are **`r fmt(diag$parameter_checks$minimum_finite_bulk_ess, 0)`** and **`r fmt(diag$parameter_checks$minimum_finite_tail_ess, 0)`**, respectively. There are **`r diag$sampler_checks$divergent_transitions` divergent transitions** and **`r diag$sampler_checks$max_treedepth_hits` maximum-tree-depth events** among 3,000 retained iterations.

The 2026 spawning-biomass sample median is **6.14 million tonnes**, with an exploratory central 95% interval of **4.47–8.43 million tonnes**. Its R-hat is **1.018**, so this quantity also needs further between-chain agreement checks. Recruitment at age 1 has a sample median of **5.38 billion fish**; detailed units, intervals and diagnostics are given below.

This separate working document accompanies the assessment and provides an analysis for Scientific Committee review.

## Model and sampler

The [SC14 assessment source](https://github.com/SPRFMO/jjm/blob/main/doc/SC14.qmd) identifies model 1.06 as the final model in the working comparison. This run uses its **single-stock hypothesis**, with data from 1970 through 2026, control file `h1_1.06.ctl`, data file `1.06.dat`, and the saved fitted parameter file. The run retains the working assessment’s 2026 catch projections and other model inputs. This analysis covers the single-stock hypothesis.

The fitted objective was `r sprintf('%.11f', as.numeric(sub('.*Objective function value = ([^ ]+).*', '\\1', manifest$original_parameter_header)))`. An isolated evaluation at the saved parameters reproduced it to within $2\times10^{-9}$. The original fitted maximum gradient component was approximately 0.000167. The preparation evaluated the saved parameter values using zero optimization iterations. The original fitted gradient supplies the convergence evidence for that mode.

The target is proportional to the exponential of the negative complete ADMB objective over the active parameters and their existing bounds. Existing likelihoods, penalties and priors remain in place; ADMB applies the required Jacobian for its sampling transformations. Parameters fixed in the control file remain fixed, including natural mortality (0.28 per year), steepness (0.65), and recruitment deviation standard deviation (0.6). The run uses the existing prior specification.

[ADNUTS](https://cole-monnahan-noaa.github.io/adnuts/articles/adnuts.html) provides NUTS sampling for ADMB models. The installed package's `sample_nuts()` defaults were used directly.

| Setting | Value |
|:--|:--|
| Model executable | ADMB 13.2, `jjm2` |
| R package | adnuts `r manifest$version` |
| Chains | 3 |
| Iterations per chain | 2,000 total; 1,000 warmup; 1,000 retained |
| Total retained draws | 3,000 |
| Initial mass matrix | Unit diagonal |
| Mass matrix adaptation | Diagonal adaptation during warmup |
| Target acceptance probability | 0.8 |
| Maximum tree depth | 12 |
| Step size | Adapted |
| Thinning | 1 |
| Initialization | Package default: all chains start from the saved fitted mode |
| Chain seeds | `r paste(manifest$package_generated_chain_seeds, collapse=', ')` |
| Sampling wall time | `r fmt(manifest$wall_seconds / 3600, 2)` hours |

The root R seed was 20260910; the package generated the individual chain seeds. Sampler controls retained their default values throughout this run. The `-fut_sel 3` model argument preserves the assessment's future-selectivity convention. Interpretation focuses on historical fitted quantities. Derived quantities were evaluated from retained draws after sampling, keeping the sampler's default `mceval = FALSE` unchanged.

## Do the chains agree?

R-hat measures agreement between chains, while effective sample size (ESS) accounts for autocorrelation. The screening used here requires rank-normalized split R-hat no greater than 1.01 and bulk and tail ESS of at least 400. The ESS threshold of 400 is the screening choice for this report. These checks follow the definitions in the [Stan diagnostic guidance](https://mc-stan.org/learn-stan/diagnostics-warnings.html). Convergence assessment combines these numerical screens with chain plots and checks of model exploration. All chains start at the same mode under the requested defaults; suitably dispersed starts would provide an additional check of exploration across modes.

```{r}
#| label: tbl-parameter-diagnostics
#| tbl-cap: "Parameter checks across 1,593 active parameters. Log posterior is assessed separately. Failure counts overlap."
p <- diag$parameter_checks
kable(data.frame(Check=c("Pass all parameter screens", "R-hat > 1.01", "Bulk ESS < 400", "Tail ESS < 400", "Undefined diagnostics", "Nonfinite parameter draws", "Constant parameter in any chain"),
                 Parameters=c(p$passing,p$rhat_above_1_01,p$bulk_ess_below_400,p$tail_ess_below_400,p$undefined_diagnostics,p$parameters_with_nonfinite_draws,p$parameters_constant_in_any_chain)))
```

```{r}
#| label: tbl-chain-diagnostics
#| tbl-cap: "Sampler checks use only the 1,000 retained iterations in each chain."
x <- chains[, c("chain", "retained_draws", "divergent_transitions", "max_treedepth_hits", "ebfmi", "mean_accept_stat", "mean_leapfrog_steps")]
names(x) <- c("Chain", "Retained", "Divergences", "Depth-12 events", "E-BFMI", "Mean acceptance", "Mean leapfrog steps")
kable(x, digits=c(0,0,0,0,3,3,0))
```

Divergences can signal biased exploration. Maximum-tree-depth events indicate that the sampler reached its trajectory-length limit. Energy Bayesian fraction of missing information (E-BFMI) is computed within each retained chain; values below 0.3 are flagged. The full screening requires zero divergences, zero depth-limit events, finite diagnostics, and E-BFMI of at least 0.3 in every chain.

```{r}
#| label: tbl-slow-parameters
#| tbl-cap: "Ten parameters with the weakest between-chain agreement, followed by lowest bulk ESS. Native ADMB parameter scales are retained."
o <- order(parameter_rows$passes_parameter_checks,
           -ifelse(is.finite(parameter_rows$rhat), parameter_rows$rhat, Inf),
           ifelse(is.finite(parameter_rows$ess_bulk), parameter_rows$ess_bulk, -Inf))
x <- head(parameter_rows[o,c("variable","rhat","ess_bulk","ess_tail","mcse_mean")],10)
names(x) <- c("Parameter","R-hat","Bulk ESS","Tail ESS","MCSE of mean")
kable(x, digits=c(0,3,0,0,4), row.names=FALSE)
```

![Six parameter traces, including the weakest parameter diagnostics, recruitment scale and log posterior. Each chain is distinguished by color and line type; warmup is shaded.](trace-selected.png){#fig-traces fig-alt="Separate chain traces for selected model parameters and log posterior. A shaded region marks the first 1000 warmup iterations, which are excluded from all reported summaries."}

![Retained-draw marginal distributions for the same parameters, shown separately by chain. Agreement in these plots is assessed alongside the numerical diagnostics.](marginals-selected.png){#fig-marginals fig-alt="Chain-specific density curves for six selected quantities using only retained draws. Different line types and colors identify the three chains."}

## Biomass and recruitment uncertainty

These summaries use all 3,000 retained draws, with chain membership preserved for diagnostics. Intervals are the equal-tailed 2.5th and 97.5th percentiles. These intervals describe the sampled draws. Further sampling is needed to support reliable posterior credible intervals.

```{r}
#| label: tbl-terminal-results
#| tbl-cap: "2026 summaries from all 3,000 retained draws. Intervals are exploratory because the full-model sampling checks fail."
x <- terminal[match(c("SSB", "Recruits"), terminal$type), ]
kable(data.frame(Quantity=x$quantity, Units=x$units,
                 "Fitted mode"=fmt(x$mode, 0),
                 "Sample median"=fmt(x$median, 0),
                 "2.5th–97.5th percentiles"=paste0(fmt(x$lower95,0), "–", fmt(x$upper95,0)),
                 "R-hat"=fmt(x$rhat,3), "Bulk ESS"=fmt(x$ess_bulk,0),
                 "Tail ESS"=fmt(x$ess_tail,0), check.names=FALSE))
```

The original fitted values are 6,108.87 thousand tonnes of spawning biomass and 5,559.86 million recruits at age 1 in 2026. The mode curves below were regenerated at the saved fitted parameters and checked against the original objective. The central 50% band shows the interquartile range; the wider band contains the central 95% of retained values.

The derived-quantity screen flags **`r derived_diag$flagged_quantity_years` of `r derived_diag$assessed_quantity_years` quantity-year combinations** (spawning biomass and recruitment over 1970–2026). Its maximum R-hat is `r fmt(derived_diag$max_rhat, 3)`, and minimum bulk/tail ESS are `r fmt(derived_diag$min_ess_bulk, 0)`/`r fmt(derived_diag$min_ess_tail, 0)`. Overall sampling assessment combines these derived-quantity diagnostics with the full set of parameter and sampler checks.

![Spawning biomass over 1970–2026: retained-sample median, central 50% and 95% intervals, and the original fitted mode.](ssb-trajectory.png){#fig-ssb fig-alt="Historical spawning biomass in thousand tonnes, with sample median and shaded central 50 and 95 percent intervals. A dashed line marks the original fitted mode; the intervals are exploratory given incomplete chain agreement."}

![Recruitment at age 1 over 1970–2026: retained-sample median, central 50% and 95% intervals, and the original fitted mode.](recruitment-trajectory.png){#fig-recruitment fig-alt="Age-one recruitment in million fish by year, with retained-sample median, central 50 and 95 percent intervals, and dashed fitted-mode estimates."}

![The three retained chains for 2026 spawning biomass. The horizontal reference is the original fitted mode.](terminal-ssb-trace.png){#fig-terminal-trace fig-alt="2026 spawning biomass from each of the three chains across 1000 retained iterations, with a horizontal fitted-mode reference. Chain line styles and colors allow assessment of overlap and drift."}

![Chain-specific distributions for 2026 spawning biomass, using retained draws only.](terminal-ssb-density.png){#fig-terminal-density fig-alt="Three separate density curves for retained 2026 spawning-biomass draws, with a fitted-mode vertical reference, allowing visual comparison of chain locations and spreads."}



## 100 Kobe trajectories, 2007–2026 {#kobe-trajectories}

The Kobe display follows **100 retained draws through the last 20 model years, 2007–2026**. Draws were sampled without replacement using seed 20260910: 33 from chain 1, 33 from chain 2 and 34 from chain 3. Each trajectory follows the same parameter draw across all 20 years. These are historical fitted trajectories conditional on the base model; the full-run mixing limitations also apply here.

[Explore the interactive Kobe plot](kobe-100-trajectories-2007-2026.html) to select a trajectory and inspect its year-specific values. Download [compact Kobe-format data](kobe-format.csv) (`iter`, `year`, `stock`, `harvest`) or [all quantities and draw identities](kobe-trajectories-100-draws-2007-2026.csv). Each CSV has **2,000 rows: 100 draws × 20 years**.

The horizontal coordinate is spawning biomass divided by the draw's mean annual MSY spawning biomass over **2017–2026**, following the assessment's fixed-biomass-reference convention. The vertical coordinate is fishing mortality divided by the **matching annual FMSY**. Fishing mortality is the unweighted mean across all 12 model ages after summing the four fleets:

$$
\mathrm{stock}_{d,y}=\frac{SSB_{d,y}}{\operatorname{mean}_{t=2017}^{2026}(SSBMSY_{d,t})},
\qquad
\mathrm{harvest}_{d,y}=\frac{\frac{1}{12}\sum_{f=1}^{4}\sum_{a=1}^{12}F_{d,f,y,a}}{FMSY_{d,y}}.
$$

Both coordinates are dimensionless. Reference lines at one mark each denominator. The thick curve joins the coordinate-wise medians of the 100 selected draws for each year; each coordinate is summarized separately to provide a visual guide to the selected draws.

![One hundred matched Kobe trajectories for 2007–2026. Thin lines show selected draws; the thick line joins their annual coordinate-wise medians. The open circle marks 2007 and the filled diamond marks 2026.](kobe-100-trajectories-2007-2026.png){#fig-kobe fig-alt="Kobe plot of 100 retained posterior trajectories over 2007 to 2026. Horizontal axis is spawning biomass relative to each draw's mean 2017 to 2026 MSY biomass; vertical axis is fishing mortality relative to annual FMSY. The median moves from approximately 0.54 and 3.50 in 2007 to 0.90 and 0.86 in 2026. Reference lines divide the plane at one; these are exploratory results because full-run chain screening failed."}

Annual reference points come from the frozen model's existing annual MSY routine. The original evaluation omitted annual FMSY, so a reporting-only copy emitted that quantity and the native F/FMSY check for the selected draws. The added output leaves their spawning biomass, fishing mortality and annual MSY biomass identical to the original evaluation. All annual reference points are positive and none reaches the routine's fishing-mortality grid limits. See the [Kobe methods and data dictionary](kobe-README.md), [selection manifest](selection-manifest.csv), [validation record](kobe-status.json), and [supplementary reproduction files](kobe-reproducibility.zip).

## Probability of Kobe green through time {#kobe-green-probability}

For each year from **2007 through 2026**, the estimated posterior probability of being in the green quadrant is the proportion of **all 3,000 retained draws** with spawning biomass above its MSY reference **and** fishing mortality below FMSY:

$$
\widehat{P}(\mathrm{green}_y\mid\mathrm{data},\mathrm{model})=
\frac{1}{3000}\sum_{d=1}^{3000}
\mathbf{1}\!\left[\mathrm{stock}_{d,y}>1\;\text{and}\;\mathrm{harvest}_{d,y}<1\right].
$$

The two conditions are evaluated jointly within each draw, preserving their posterior dependence. The reference definitions match the Kobe trajectories above: each draw's mean SSBMSY for 2017–2026 and annual FMSY. Both inequalities are strict. Each chain contributes 1,000 retained draws, giving equal chain weights.

The **2026 estimate is `r fmt(100*green$probability_green[green$year==2026],1)`%**, corresponding to **`r green$green_draws[green$year==2026]` of 3,000 draws** in green. These are exploratory posterior estimates conditional on this base model. Further sampling is needed to improve the full model's chain agreement and strengthen the probability estimates.

![Estimated posterior probability of the Kobe green quadrant for each year, using all 3,000 retained draws. The thick green line and points give the pooled estimate; thin patterned lines give estimates from each chain.](kobe-green-probability-2007-2026.png){#fig-kobe-green fig-alt="Annual estimated posterior probability of simultaneously exceeding the MSY spawning-biomass reference and having fishing mortality below annual FMSY over 2007 to 2026. A thick green line shows the proportion across all 3000 retained draws, and three patterned lines show individual-chain proportions. The downloadable table gives every annual value and Monte Carlo precision estimate."}

```{r}
#| label: tbl-kobe-green
#| tbl-cap: "Annual estimated posterior probabilities of Kobe green. Each year uses 3,000 retained draws. Monte Carlo standard errors (MCSE) are expressed in percentage points."
gf <- function(z) ifelse(is.finite(z), formatC(z,format="f",digits=1), "—")
kable(data.frame(Year=green$year,
  "Green draws"=green$green_draws,
  "Probability (%)"=gf(100*green$probability_green),
  "Chain range (%)"=paste0(gf(100*pmin(green$chain1_probability,green$chain2_probability,green$chain3_probability)),"–",
                            gf(100*pmax(green$chain1_probability,green$chain2_probability,green$chain3_probability))),
  "MCSE (percentage points)"=gf(100*green$probability_mcse),check.names=FALSE),row.names=FALSE)
```

MCSE measures Monte Carlo precision of the estimated probability while accounting for correlation within chains. The calculation uses the binary green-quadrant indicator and preserves chain identity. The chain range shows the smallest and largest of the three chain-specific estimates. For years with a constant indicator across retained draws, the table reports the observed proportion and marks MCSE with a dash. A zero count describes the finite set of sampled draws. The full-model diagnostics above provide the context for interpreting these probabilities.

Download the [annual probability table](kobe-green-probability-2007-2026.csv), [chain-specific estimates](kobe-green-by-chain-2007-2026.csv), [all draw-level coordinates and green indicators](kobe-green-draws-2007-2026.csv.gz), or [vector figure](kobe-green-probability-2007-2026.svg). The [calculation and evaluation scripts](kobe-green-reproducibility.zip) and [validation record](kobe-green-status.json) document the derivation. The previously displayed 100 trajectories retain their original draw selection; the probability calculation uses the full retained sample.

## Interpretation and next steps

The chain-specific marginal plots show a shift in chain 1 for several of the slowest fishing-mortality and survey-catchability parameters. This agrees with the elevated R-hat and low ESS. Increasing the retained run length and checking suitably dispersed starting points are reasonable next checks; either would be a separately labelled follow-up run. The current analysis preserves the requested defaults.

## What uncertainty is covered?

MCMC explores parameter uncertainty conditional on this single-stock model. Alternative stock hypotheses, additional processes, different data treatments and fixed parameter choices can be assessed through separate sensitivity analyses. The defaults initialize every chain at the same mode. A run that needs further work should be followed by a separately labelled analysis with justified changes to initialization, run length or sampler controls.

The Kobe calculation constructs stock-wide fishing mortality from F at age across all four fleets. The model's native `Fbar` field describes a fishery-specific parameter, and its depletion fields describe total-biomass ratios. The source notes record these definitions so each reported quantity retains its biological meaning.

## Files and reproducibility

- [Editable Quarto source](base-model-adnuts.qmd)
- [Recorded run settings](run-manifest.json) and [sampling diagnostic summary](diagnostic-summary.json)
- [All parameter summaries](parameter-summary.csv) and [chain diagnostics](chain-diagnostics.csv)
- [Derived quantity summaries](derived-summary.csv) and [source notes](source-notes.md)
- [Kobe data in R format](kobe-trajectories.rds), [vector figure](kobe-100-trajectories-2007-2026.svg), and [trajectory extraction script](build-kobe-trajectories.R)
- [Software environment](session-info.txt) and [source/output checksums](release-manifest.json)
- [Reproduction instructions](reproduction-README.md) and [reproducibility files](reproducibility.zip), including the exact model source, control, data, saved parameter input, and analysis scripts

The raw fit and full ADMB evaluation output are retained in the local run archive. Their checksums identify the draws used here. Report rendering uses the saved summaries and figures from the completed run.

## Validation record

| Check | Outcome | Evidence |
|:--|:--|:--|
| Base-model identity and inputs | Pass | Frozen file hashes; original and prepared objective agree within 2 × 10⁻⁹ |
| Requested sampler defaults | Pass | Recorded package formals, resolved controls and all three executed commands |
| Sampling completion and draw order | Pass | 2,000 × 3 iterations; all 4,779,000 pooled parameter values match the retained fit bit for bit |
| Numerical sampler diagnostics | Pass | Zero retained divergences and zero depth-limit events; E-BFMI ≥0.984 |
| Parameter mixing and precision | Fail | 366/1,593 active parameters fail one or more reported screens |
| Native derived-output completeness | Pass | All 3,000 draws and 57 years verified; all 20 native output types inventoried |
| Kobe draw selection and coordinates | Pass | 100 distinct retained draws; 20 years each; unchanged native quantities and matching annual reference points |
| Annual Kobe green probability | Pass | All 3,000 retained draws and 20 years; joint event counts, chain weights and original 100-draw coordinates verified |
| Reported results | Pass | Tables and figures read the saved, checked CSVs; independent parameter diagnostic calculation agrees |
| HTML, navigation, figure descriptions and file links | Pass | Rendered artifact, figure review and scoped local-link validation recorded in the release checks |
| Manual screen-reader review | Pending | A manual review remains an additional validation step |

These checks document the numerical, file and visual validation completed for this analysis.