Fish growth (von Bertalanffy) and length-weight condition

1. Linf needs old fish. With few fish near the asymptote (hard_nls, hard_nimble_fixed) Linf, and so age at 90% Linf, is not identified: nls occasionally returns absurd Linf with wide intervals (age at 90% Linf fails the gate); the Bayesian model returns Linf biased +16 to +26% toward its prior. The module notes when < 5 fish are at or beyond the estimated age at 90% Linf; treat Linf as unreliable then and compare groups by length at ages that were sampled.
2. Length-selective gear biases growth (section 4): K -14%, t0 -0.36 years, coverage down to 0.58. Selectivity is not modelled; drop partially selected ages.
3. Single-parameter group differences have low power (section 6): 0.59 for a 50 mm Linf difference with 250 fish per group. Use Kimura's joint test or whole-curve comparisons.
4. Hierarchical SDs with few groups. At 10 groups the posterior median of the among-group SDs is biased +6 to +7% (passes); the posterior mean was +14% (failed). With 6 groups x 40 fish sd_log_Linf is biased +24.5% and FAILS the gate (section 5): with fewer than 10 groups (independent verification: 8 groups x 70 fish also FAILS, +14.7% and +16.2%), interpret the among-group SD of Linf through its upper 95% limit ("variation in Linf is at most about N%"), not as a point value. nlme's REML gives smaller SDs than the Bayesian median on the example data.
5. Hierarchical convergence needs longer chains than the fixed-group model: on the example data (12 groups, 2106 fish) 10,000 iterations gave max R-hat 1.010 and min bulk ESS 335 (t0), short of the real-data standard, so the hierarchical default is 20,000 iterations with 5,000 burn-in.
6. Not modelled: ageing error, individual growth variation beyond the residual, Linf-K correlation among groups, growth models other than von Bertalanffy, length-at-capture date beyond age_offset.

Date: 2026-10-19. R 4.5.3, FSA 0.10.1, nimble (ecotools Rlib), fishmethods and nlme for cross-checks. Gate: recovery_gate() with its fixed thresholds (|bias| <= max(10% of truth, 3 MC SE), or 0.1 x empirical SD when the truth is 0; coverage >= 0.95 - 3 MC SE, i.e. 0.921 for 500 reps and about 0.884 for 100 reps; >= 90% of reps converged). Recovery uses fast_settings (Nimble: 3 chains x 3,000 iterations, 1,000 burn-in) unless a config overrides it; a Bayesian rep counts as converged when max R-hat < 1.05.

Code state: the nls runs and default_nimble_fixed used the files hashed in dev/fisheries-growth/code_md5_validation.txt. On 2026-10-19 the hierarchical engine was changed to report the posterior median of sd_log_Linf and sd_log_K (section 3), and hier_nimble, null_hier, hard_nimble_fixed and hard_hier were rerun with dev/fisheries-growth/code_md5_rerun.txt. That change touches only the hierarchical point estimates of the two SDs, so the nls and fixed-group runs are unaffected. Later edits touch only text in output.R and default settings in module.R (final hashes: dev/fisheries-growth/code_md5_final.txt).

Reproduce: bash dev/fisheries-growth/run_validation.sh <run> <nsim> with configs in modules/fisheries-growth/validation/<run>.yaml (outputs in results/fisheries-growth/validation/<run>/), bin/Rscript dev/fisheries-growth/kimura_fpr.R, bin/Rscript dev/fisheries-growth/crosscheck.R. Tables come from dev/fisheries-growth/validation_tables.R. Times per fit: nls 0.3 s; Bayesian fixed groups (3 groups x 150 fish) 12 s; hierarchical (10 groups x 80 fish) 35 s, all at fast settings.

Default design (module defaults): 3 groups x 150 aged fish; catch-at-age P(a) ∝ exp(-0.35(a-1)), ages 1-15 (about 13 fish per group at age >= 8); Linf 600/650/600 mm, K 0.30/0.25/0.30, t0 -0.5; multiplicative error sigma 0.08 (log scale); every fish weighed, W = 3.5e-6 L^3.18 with condition 1/1.05/1 and log-scale SD 0.08.

core/run.R selftest fisheries-growth (4 reps): contract checks pass; with 4 reps the coverage bound is 0.80 and some estimands fail on Monte Carlo noise alone, as expected (results/fisheries-growth/selftest/).

500/500 fits, no errors.

Other default-design runs (full tables in the appendix):

dev/fisheries-growth/kimura_fpr.R; reject at p < 0.05. Groups with identical growth (Linf 600, K 0.30, t0 -0.5 for all three groups):

All within 2 MC SE of 0.05. Power at the default design (Linf 600 vs 650, K 0.30 vs 0.25; 500 reps): all common 0.996, common Linf 0.440, common K 0.362, common t0 0.058 (t0 is truly equal).

Design hier_nimble: 10 groups (lakes/years) x 80 fish, log Linf_g ~ N(log 600 - 0.05 z_g, 0.08), log K_g ~ N(log 0.30 + 0.10 z_g, 0.15), common t0 -0.5, group covariate "temperature" (mean 20, SD 2, standardised). nsim = 100, 96/100 reps with R-hat < 1.05.

PASS on all 69 estimands (group-level Linf_g, K_g, age at 90% Linf, differences, condition ratios included; appendix). The hyper-parameters:

History, reported so the change is visible: the first run (same seeds) reported the posterior mean of the SDs and FAILED them on bias: sd_log_Linf +14.6% (MCSE 3.2%), sd_log_K +13.5% (MCSE 3.0%), coverage 0.969/0.969. With 10 groups the SD posterior is right-skewed and its mean sits above the truth. We now report the posterior median for these two SDs only (README, methods text). The thresholds were not changed. With 6 groups the median is still biased (section 5).

Negative control null_hier (same design, Linf_beta = K_beta = 0; 98/100 converged): the 95% credible interval excluded 0 in 4.1% (Linf_beta) and 3.1% (K_beta) of reps; all 69 estimands PASS (sd_log_Linf +5.5%, sd_log_K +5.6%).

misspec_selectivity: the same default design, but the gear never catches fish < 250 mm, while the model assumes each age class is sampled at random. nls engine, 500 reps. 12 of 25 estimands FAIL:

Length-weight and condition still PASS (the selectivity does not bias W given L). This is the failure mode the README warns about.

hard_nls: few old fish (60 fish/group, Z = 0.9; about 0.1 fish per group at age >= 8). 494/500 fits (6 errors: "nls fit failed from every starting value", reported to the user as such). 23/25 estimands pass the gate, but this hides the real problem: with almost no fish near the asymptote, Linf is not identified in a minority of data sets. Median Linf_A is 596 (truth 600) but 6.3-7.1% of fits put Linf above 1.5 x truth (mean Linf_A 2,380 mm; the bias passes only because 3 MC SE is then enormous), the median 95% CI is 56-68% of Linf wide, and age at 90% Linf FAILS (coverage 0.907 and 0.909 < 0.921; mean 55 vs 7.2 years, driven by those fits). Condition, length-weight and the K/t0 estimates pass. The module flags this situation in every fit (note when < 5 fish are at or beyond the estimated age at 90% Linf).

hard_nimble_fixed: same design, Bayesian fixed groups (97/100 converged). 9 FAIL: Linf biased +16 to +26% (posterior mean 722-757 vs 600-650), K_C -12.5%, t0 -21 to -31%, age at 90% Linf +33 to +59%; coverage stays 0.92-0.97. When the data do not identify Linf the posterior follows the prior N(log max length, 1), which is wide and skewed upward. Condition, length-weight, differences between groups and the residual SD pass.

hard_hier: hierarchical with 6 groups x 40 fish, Z = 0.6 (15,000 iterations, 5,000 burn-in). 97/100 reps converged (the first attempt at the fast setting of 3,000 iterations converged in only 44/100, so this design needs longer chains; 57 s per fit). 44 of 45 estimands PASS; sd_log_Linf FAILS: posterior median 0.0996 vs truth 0.08 (+24.5%, MCSE 5.8%), coverage 0.948. sd_log_K passes (+5.7%, coverage 0.959). Linf_typical +1.4% (coverage 0.969), K_typical -0.1%, t0 -1.1%, Linf_beta -0.8% (power 0.04), K_beta -6.0% (power 0.13); group-level Linf_g, K_g, age at 90% Linf all pass. With 6 groups the data barely inform the among-group SD of Linf, and the half-normal(0, 0.5) prior pulls it upward; the 95% interval is still honest. With 6 groups the covariate effects have power <= 0.13 for effects of this size.

--vary n_per_group=30,60,100,150,250 at the default parameters (results/fisheries-growth/power/). Proportion of reps whose 95% CI excludes the null:

Coverage stayed 0.91-0.98 at every sample size. Single-parameter differences need many fish because Linf and K trade off; Kimura's joint test that the curves differ has power 0.996 at 150/group (section 2).

dev/fisheries-growth/crosscheck.R on the example data (Lake Erie walleye 2014, female vs male, n = 1006; and 2003-2014 for the hierarchical model, n = 2106):

• example/config.yaml (nls, sexes, 2014): report results/fisheries-growth/example/. Female Linf 680 mm (661-699), K 0.312; male Linf 559 mm (548-569), K 0.390; males reach 90% of Linf at 4.8 years vs 6.4. Kimura: Linf and K differ (chi-square 113.4 and 8.8), t0 does not (p = 0.21). b = 3.34; males 1.8% heavier at the same length. Every sentence in the report was checked against the estimates table; figures inspected.
• example/config_hierarchical.yaml (Nimble, year as random effect, calendar year as covariate): 20,000 iterations, 5,000 burn-in, 3 chains: max R-hat 1.002, min bulk ESS 581 (meets the real-data standard; at 10,000 iterations it did not: R-hat 1.010, ESS 335). Typical year: Linf 584 mm (565-605), K 0.402 (0.362-0.445); among-year SD 0.049 on log Linf (about 5%) and 0.126 on log K (about 13%); calendar year: +0.3% Linf per SD (CI -3.2 to +3.8%), +5.7% K per SD (CI -3.1 to +15.3%, P(positive) 0.91), i.e. no clear trend. 2004 stands out (Linf 636 mm, K 0.294). Condition fell: 2014 fish weigh 7.5% less at the same length than 2003 fish (CI -9.1 to -5.8%). Report sentences checked against the table; conditional-effects plots (Linf and K vs year, response scale with 50/95% bands and shrunken group estimates) inspected.

• Both example reports in results/fisheries-growth/example*/ were regenerated on 2026-10-19 with the final verified code (hashes in dev/fisheries-growth/verify/code_md5_verified.txt); estimates are unchanged (hierarchical: max R-hat 1.002, min bulk ESS 581). The final hierarchical report also notes that year 2005 has only 4 fish at or beyond its estimated age at 90% of Linf.

Script: dev/fisheries-growth/verify/verify.R (own simulator, no module code used to simulate; module fit() + core summarise_recovery()/recovery_gate()). Outputs in dev/fisheries-growth/verify/out/; logs logs/fisheries-growth_verify_*.log.

Literature / code check. VB typical parameterisation L = Linf(1 - exp(-K(age - t0))) matches FSA (makeGrowthFun, param 1); age at 90% Linf = t0 + ln10/K is correct. Kimura LRT -2 ln Lambda = N ln(RSS_H/RSS_general) with (G-1) df per constrained parameter matches fishmethods::growthlrt source (-N*log(RSS1/RSS_k)); rerun on the 2014 walleye data with additive error: 132.448 / 4.968 / 5.644 / 638.565 vs fishmethods 132.45 / 4.97 / 5.64 / 638.56; parameter estimates and SEs identical (Linf_female 684.05, SE 7.575). Walleye Ws from FSA::wsVal: log10 Ws = -5.453 + 3.18 log10 TL (Murphy et al. 1990, RLP, min 150 mm); hand Wr means 88.77 / 90.30 = report. Length-weight bias correction exp(s^2/2) and condition ratio exp(group coefficient) are correct. simulate() and fit() share no model code (fit uses formula strings / Nimble code); the independent simulator below rules out cancelling errors.

Independent planted truth (design differs from the builder's: 4 unequal groups north/south/ east/west, n = 100/120/80/150, ages 2-18 drawn as recruitment age + geometric survival (Z 0.25) plus a 0.4-yr capture offset, Linf 450/480/450/520, K 0.20/0.20/0.26/0.17, t0 -1.0/-0.6/-1.0/-0.8, W = 10^-5.2 * cond * L^3.05, lw sigma 0.10, condition 1/0.94/1.06/1, 70% of fish weighed, reference length 400):

Builder's runs rerun. recover with validation/default_nls.yaml (500 reps) and validation/hier_nimble.yaml (100 reps) reproduced the recorded summaries exactly (max abs difference in mean estimates 0; PASS 25/25 and 69/69).

Controls.

• Identical growth and condition in 3 groups (nls, 1,000 reps): Wald false-positive rates 0.043-0.052 for every Linf/K/t0 difference and condition ratio; Kimura tests reject 0.055 (Linf), 0.061 (K), 0.056 (t0), 0.053 (all) (binomial 95% range at n = 1000: 0.037-0.064).
• Hierarchical covariate effects = 0 (8 lakes, 97 converged): false-positive rate 0.031 (Linf), 0.062 (K).
• Negative control, ageing error (fish older than 6.4 read 1 year young, older than 9.4 read 1-3 years young; not modelled): 13/33 FAIL (Linf up to +9.6% with coverage 0.55, age at 90% Linf +10 to +15%, K coverage 0.85, residual SD coverage 0.76). The gate can fail.

Usability. fit on example/config.yaml reproduces the stored estimates exactly. Every sentence of both example reports was checked against the estimates (Linf, 90% Linf, differences, Kimura, b, condition %, Wr, covariate % effects = 100(exp(beta)-1), SD of calendar year 3.61); all correct, association (not causal) wording for covariates. Figures viewed (growth curves with CI and PI bands, conditional effects of calendar year on Linf/K on the response scale, condition). Messy-data test (dev/fisheries-growth/verify/messy/): wrong length column and duplicated id already gave clear messages. Fixed in module.R (read_data/check_data only; fit and simulate unchanged):

• a missing filter: column now stops with the column list (was "no rows left");
• group values not in group_levels are listed with a "did you mean" hint, and a rare group that looks like a typo of another (femal, Male) gets the same hint;
• NA/non-numeric age and length messages list fish ids (row numbers were post-filter and misleading);
• ages > 150 are flagged as probably months/days;
• weights more than 2.5x off the length-weight line (e.g. some weights in kg) are now a data problem (they were only a note after fitting). Both example datasets and 20 default simulations raise no problems. Also changed: output.R adds a sentence when a hierarchical fit has < 10 groups saying the among-group SDs tend to be overestimated; README/Known limit 4 corrected from "fewer than about 8 groups" to "fewer than 10 groups". Selftest (4 reps) runs. Code hashes after verification: dev/fisheries-growth/verify/code_md5_verified.txt.

Verdict: VERIFIED_WITH_LIMITS. The default design and an independent design pass for nls, additive nls and Bayesian fixed groups; the hierarchical model recovers everything except the among-group SD of log Linf below 10 groups (documented limit 4).

Script: dev/fisheries-growth/verify/v2/verify2.R (own simulator written from the VB equations; no module code used to simulate; module fit() + core summarise_recovery()/recovery_gate()). Outputs: dev/fisheries-growth/verify/v2/out/; logs logs/fisheries-growth_verify2_*.log. Code state: dev/fisheries-growth/verify/v2/code_md5_verify2.txt. Note: output.R differs from verify/code_md5_verified.txt because a figures session edited the figure code (fg_figures only) at 08:15 and 08:23 after the first verification; engines, lw.R, simulate.R and module.R are unchanged, so recovery results are unaffected. The example fit below was run with the 08:23 output.R.

Literature: VB form matches FSA::makeGrowthFun(param = 1); Kimura statistic matches the fishmethods::growthlrt source (-N*log(RSS1/RSSk)); walleye Ws (-5.453, 3.18, min 150 mm, Murphy et al. 1990 RLP) matches FSA::wsVal; hierarchical lognormal Linf/K follows Helser & Lai (2004).

Design (different from the builder's and verifier 1's): lengths in cm, weights in kg, 2 groups (inlet n = 180, outlet n = 90), dome-shaped catch-at-age P(a) ∝ a^2 exp(-0.25a), ages 3-28 + 0.5, Linf 92/80 cm, K 0.13/0.17, t0 -1.4/-0.9, W = 7e-6 cond L^3.1 (cond 1/0.93), lw sigma 0.12, 60% weighed.

With 14 groups the among-group SDs pass, consistent with Known limit 4 (problem below 10 groups).

Controls:

• Null (identical growth and condition, 1,000 reps, nls): Wald false-positive rates 0.049 (condition), 0.054 (Linf), 0.067 (K), 0.061 (t0); Kimura rejection 0.057 (Linf), 0.070 (K), 0.068 (t0), 0.070 (all). Binomial 95% range at 1,000 reps is 0.037-0.064, so with a small second group (90 fish) the nominal 5% tests are mildly liberal (about 6-7%). Treat p-values near 0.05 with caution when a group has fewer than about 100 fish.
• Negative control (data from a Gompertz curve, fitted as VB, 300 reps): 4/12 FAIL (K -21% with coverage 0.10 and -15% with coverage 0.73; K difference +13%). The gate can fail.

Usability: fit on example/config.yaml reproduces the stored estimates exactly. Every number in the report's interpretation was checked against the estimates table (Linf 680/559, 90% Linf 612/503 mm at 6.4/4.8 yr, Linf difference 121 mm (99-143), K difference 0.078, condition +1.8% (+0.7 to +2.9%), mean weight 1223 g, Wr 88.8/90.3); wording is correct and non-causal. Figures (growth curves with confidence and prediction bands drawn over observed ages, parameter comparison, condition) were viewed and agree with the table. Messy data (dev/fisheries-growth/verify/v2/messy/): NA, blank and "4+" ages; misspelled/case/space group values ("Femal", "Female", " male"); weights in kg; a renamed length column; a misspelled filter column; "512mm" and negative lengths all gave clear messages that name fish ids and suggest fixes. One cosmetic issue: a negative length also prints an R warning "NaNs produced" from the weight-outlier check.

Changes made: none to module code (no blocking problems found).

Verdict: VERIFIED_WITH_LIMITS (agrees with verifier 1; limits as listed above, plus mild liberalness of the group-difference tests with small groups).

default_nls: nsim 500, fit errors 0, converged 500 (nls, no R-hat), coverage bound 0.921

additive_nls: nsim 500, fit errors 0, converged 500 (nls, no R-hat), coverage bound 0.921

default_nimble_fixed: nsim 100, fit errors 0, converged 100 (R-hat < 1.05 share 1.00), coverage bound 0.885

hard_nls: nsim 500, fit errors 6, converged 494 (nls, no R-hat), coverage bound 0.921

misspec_selectivity: nsim 500, fit errors 0, converged 500 (nls, no R-hat), coverage bound 0.921

null_nls: nsim 500, fit errors 0, converged 500 (nls, no R-hat), coverage bound 0.921

hier_nimble: nsim 100, fit errors 0, converged 96 (R-hat < 1.05 share 0.96), coverage bound 0.883

null_hier: nsim 100, fit errors 0, converged 98 (R-hat < 1.05 share 0.98), coverage bound 0.884

hard_nimble_fixed: nsim 100, fit errors 0, converged 97 (R-hat < 1.05 share 0.97), coverage bound 0.884

hard_hier: nsim 100, fit errors 0, converged 97 (R-hat < 1.05 share 0.97), coverage bound 0.884