de— title: “Bayesian joint growth–survival–detection model for XYTE” subtitle: “Reach 2 (Lake Mohave) Basin zone — FY 2006–2026” date: “June 11, 2026” format: html: toc: true toc-depth: 3 embed-resources: true code-fold: true execute: warning: false —
This report was generated using AI under general human direction. At the time of generation, the contents have not been comprehensively reviewed by a human analyst.
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(tidyr)
library(ggplot2)
library(scales)
library(knitr)
library(MCMCvis)
library(coda)source("LabFunctions.R")
packages(lubridate)Loading required package: lubridate
Attaching package: 'lubridate'The following objects are masked from 'package:base':
date, intersect, setdiff, unionload("data/XYTEJointModel_MCMC.RData")
# Objects: all_samples, samples_mcmc, param_summary, Rhat, surv_tbl,
# N_cjs, N_cells, T_OCC, FY_RANGE, TL_CENTER, TL_SCALE, THETA
load("data/XYTEJointModel_data.RData")
# Objects: CJS_Fish, N_cjs, y_cjs, TL_obs_cjs, age_cjs, L0_cjs,
# sex_cjs, first_occ_cjs, last_occ_cjs, BinCells, N_cells
load("data/XYTEGrowth.RData") # VBParams (nlme survivors-only estimates)
load("data/XYTElogisticSurvival.RData") # lag_results, all_curves (simple logistic)
# Derived constants — not saved in .RData but computable from saved objects
age_rel_cjs <- age_cjs[cbind(seq_len(N_cjs), first_occ_cjs)]
age_rel_bin <- age_bin[cbind(seq_len(N_cells), BinCells$first_occ_idx)]
all_draws <- do.call(rbind, all_samples)This report documents a hierarchical Bayesian state-space model fitted to PIT-scanning and physical-capture data for Razorback Sucker (Xyrauchen texanus, XYTE) released or captured in the Lake Mohave Basin zone (Reach 2, Zone 2.3) under the Lower Colorado River Multi-Species Conservation Program.
Standard growth analyses fit von Bertalanffy growth functions (VBGF) to fish that were physically recaptured at least twice. Because survival increases with body size, recaptured fish are a biased, fast-growing subset of the released population. Parameters estimated from survivors alone over-predict asymptotic size (L∞) and under-predict the growth coefficient (K) for the full population.
The same bias propagates into separate survival analyses: if growth is over-estimated in the data used to build a logistic survival model, the TL–survival slope inherits that misrepresentation.
The model described here estimates growth and survival simultaneously from the full encounter history of 4,994 detected fish plus the all-zero histories of 44,357 fish never detected — exactly the population that the simple approach excludes. By conditioning on the complete data, growth parameters adjust to be consistent with the survival process, and survival parameters adjust to be consistent with the corrected growth trajectories.
n_total <- N_cjs + sum(BinCells$N_rel)
n_sex_conf <- sum(!is.na(sex_cjs))
n_female <- sum(sex_cjs == 1L, na.rm = TRUE)
n_male <- sum(sex_cjs == 2L, na.rm = TRUE)
n_tl_obs <- sum(!is.na(TL_obs_cjs))
fy_min <- min(FY_RANGE); fy_max <- max(FY_RANGE)The analysis population comprises 49,351 XYTE released in Basin Zone 2.3 from FY 2006 through FY 2026. Of these, 4,994 (10.1%) were detected at least once by passive PIT-antenna scanning after their release date; the remaining 44,357 were never detected.
data.frame(
Component = c(
"Total Basin releases (FY 2006–2026)",
"Individual CJS fish (≥1 PIT contact)",
"Binomial collapse cells (never-detected)",
"Sex-confirmed fish in CJS",
" — Female",
" — Male",
" — Unknown sex (latent)",
"TL measurements (physical recaptures)",
"Fiscal years modelled",
"Fixed tag-loss rate (θ)"
),
Value = c(
format(n_total, big.mark=","),
format(N_cjs, big.mark=","),
paste0(N_cells, " cells (", format(sum(BinCells$N_rel), big.mark=","), " fish)"),
format(n_sex_conf, big.mark=","),
format(n_female, big.mark=","),
format(n_male, big.mark=","),
format(sum(is.na(sex_cjs)), big.mark=","),
format(n_tl_obs, big.mark=","),
paste0(fy_min, "–", fy_max, " (T = ", T_OCC, ")"),
"0.05"
)
) |> kable(align = c("l","r"))| Component | Value |
|---|---|
| Total Basin releases (FY 2006–2026) | 49,351 |
| Individual CJS fish (≥1 PIT contact) | 4,994 |
| Binomial collapse cells (never-detected) | 97 cells (44,357 fish) |
| Sex-confirmed fish in CJS | 595 |
| — Female | 381 |
| — Male | 214 |
| — Unknown sex (latent) | 4,399 |
| TL measurements (physical recaptures) | 1,651 |
| Fiscal years modelled | 2006–2026 (T = 21) |
| Fixed tag-loss rate (θ) | 0.05 |
The model has three linked sub-models. All parameters are estimated jointly from a single MCMC run.
Fish length at age \(a\) for fish \(i\) of sex \(s\) follows a release-conditioned von Bertalanffy growth function (VBGF):
\[L_{i,a} = L_{\infty,s} - (L_{\infty,s} - L_{0i})\,e^{-K_s\,(a - a_{\text{rel},i})}\]
where \(L_{0i}\) is the observed total length at release (anchoring the curve without estimating \(t_0\)), \(a_{\text{rel},i}\) is age at release, and \(L_{\infty,s}\), \(K_s\) are sex-specific asymptotic size and growth coefficient.
Sex for the 4399 unknown-sex fish is treated as a latent Bernoulli variable drawn from a shared proportion \(\phi_F\). An identifiability constraint \(K_F > K_M\) (implemented via \(K_F = K_M + K_\delta\), \(K_\delta > 0\)) prevents label switching between the two growth curves.
Observed TL at physical recapture is modelled with Normal measurement error: \(\text{TL}^{\text{obs}}_{i,t} \sim \mathcal{N}(L_{i,t},\, \sigma_{\text{obs}}^2)\)
Survival is a logistic function of current body size, shared across sexes:
\[\text{logit}(\phi_{\text{true},it}) = \beta_0 + \beta_L \cdot \frac{L_{i,t} - 450}{50}\]
Effective survival (incorporating the fixed 5% annual tag-loss rate \(\theta\)) is:
\[\phi_{\text{eff},it} = \phi_{\text{true},it} \times (1 - \theta)\]
The latent alive/tagged state \(z_{it}\) transitions as:
\[z_{it} \mid z_{i,t-1} \sim \text{Bernoulli}(z_{i,t-1} \cdot \phi_{\text{eff},it})\]
Detections from passive PIT scanning are modelled with a single Basin-wide annual detection probability \(p_{\text{det}}\), active only after a one-FY post-release suppression window:
\[y_{it} \mid z_{it} \sim \text{Bernoulli}(z_{it} \cdot p_{\text{det}} \cdot I_{it}^{\text{atLarge}})\]
The 44,357 fish never detected are summarised into 97 cells indexed by (50-mm TL bin × release FY). For each cell \(c\), the cumulative survival probability over the monitoring window is:
\[S_c = \prod_{t=t_{0c}+1}^{T} \phi_{\text{eff},c,t}\]
Because per-FY detection probability \(\approx 0.89\) implies cumulative \(P_{\text{detect}} > 99.95\%\) for fish with ≥4 scan FYs (established in prior Stalwart model work), zero detections from a cell are diagnostic of mortality:
\[\Pr(\text{never detected}) \approx 1 - S_c\]
This provides a likelihood contribution \(0 \sim \text{Binomial}(1 - S_c, N_{\text{rel},c})\) for each never-detected cell.
data.frame(
Parameter = c("L∞_F (mm)", "L∞_M (mm)", "K_M (yr⁻¹)", "K_δ = K_F − K_M",
"σ_obs (mm)", "β₀ (logit survival at 450 mm)", "β_L (per 50 mm)",
"p_det", "φ_F (proportion female)", "θ (tag loss)"),
Prior = c("Normal(715, 40)", "Normal(719, 40)", "Uniform(0.05, 0.35)", "Uniform(0, 0.30)",
"Uniform(5, 80)", "Normal(0, 2)", "Normal(1, 1)",
"Beta(30, 4)", "Beta(6, 4)", "Fixed at 0.05"),
Rationale = c("nlme survivors-only estimate", "nlme survivors-only estimate",
"Weakly informative", "Enforces K_F > K_M (identifiability)",
"Weakly informative", "logit(50%) at 450 mm",
"~1 logit-unit per 50 mm, from logistic model",
"Mean ≈ 0.88 from Stalwart model posterior", "Mean = 0.60 from confirmed-sex sample",
"Expert assumption")
) |> kable(align = c("l","l","l"))| Parameter | Prior | Rationale |
|---|---|---|
| L∞_F (mm) | Normal(715, 40) | nlme survivors-only estimate |
| L∞_M (mm) | Normal(719, 40) | nlme survivors-only estimate |
| K_M (yr⁻¹) | Uniform(0.05, 0.35) | Weakly informative |
| K_δ = K_F − K_M | Uniform(0, 0.30) | Enforces K_F > K_M (identifiability) |
| σ_obs (mm) | Uniform(5, 80) | Weakly informative |
| β₀ (logit survival at 450 mm) | Normal(0, 2) | logit(50%) at 450 mm |
| β_L (per 50 mm) | Normal(1, 1) | ~1 logit-unit per 50 mm, from logistic model |
| p_det | Beta(30, 4) | Mean ≈ 0.88 from Stalwart model posterior |
| φ_F (proportion female) | Beta(6, 4) | Mean = 0.60 from confirmed-sex sample |
| θ (tag loss) | Fixed at 0.05 | Expert assumption |
n_chains_rpt <- N_CHAINS
n_iter_rpt <- N_ITER
n_burnin_rpt <- N_BURNIN
n_thin_rpt <- N_THIN
n_post <- (N_ITER - N_BURNIN) / N_THINThe model was fit in NIMBLE (v1.4.2) using 3 independent chains × 100,000 iterations (burn-in 40,000; thinning 10), yielding 6,000 post-burn-in draws per chain.
rhat_df <- as.data.frame(Rhat) |>
tibble::rownames_to_column("Parameter") |>
rename(`Point est.` = `Point est.`, `Upper 95% CI` = `Upper C.I.`) |>
mutate(across(where(is.numeric), \(x) round(x, 3)))
kable(rhat_df, align = c("l","r","r"))| Parameter | Point est. | Upper 95% CI |
|---|---|---|
| K_F | 1.000 | 1.002 |
| K_M | 1.003 | 1.007 |
| K_diff | 1.000 | 1.000 |
| Linf_F | 1.000 | 1.000 |
| Linf_M | 1.004 | 1.012 |
| beta0 | 1.004 | 1.012 |
| betaL | 1.002 | 1.007 |
| beta_post | 1.001 | 1.001 |
| p_det | 1.000 | 1.001 |
| prop_F | 1.000 | 1.001 |
| sigma_obs | 1.002 | 1.008 |
All \(\hat{R}\) point estimates are ≤ 1.004 and all upper 95% CIs are ≤ 1.012, indicating satisfactory convergence across all parameters. Effective sample sizes ranged from 1877 (survival intercept) to 1.5258^{4} (annual detection probability).
Linf_F_j <- param_summary["Linf_F", "mean"]
Linf_M_j <- param_summary["Linf_M", "mean"]
K_F_j <- param_summary["K_F", "mean"]
K_M_j <- param_summary["K_M", "mean"]
sig_obs <- param_summary["sigma_obs", "mean"]
growth_tbl <- data.frame(
Parameter = c("L∞_F (mm)", "L∞_M (mm)", "K_F (yr⁻¹)", "K_M (yr⁻¹)", "σ_obs (mm)"),
`Joint model — mean` = c(
round(Linf_F_j, 1), round(Linf_M_j, 1),
round(K_F_j, 3), round(K_M_j, 3), round(sig_obs, 1)),
`Joint model — 95% CI` = c(
paste0("(", round(param_summary["Linf_F","2.5%"],1), ", ", round(param_summary["Linf_F","97.5%"],1), ")"),
paste0("(", round(param_summary["Linf_M","2.5%"],1), ", ", round(param_summary["Linf_M","97.5%"],1), ")"),
paste0("(", round(param_summary["K_F","2.5%"],3), ", ", round(param_summary["K_F","97.5%"],3), ")"),
paste0("(", round(param_summary["K_M","2.5%"],3), ", ", round(param_summary["K_M","97.5%"],3), ")"),
paste0("(", round(param_summary["sigma_obs","2.5%"],1), ", ", round(param_summary["sigma_obs","97.5%"],1), ")")),
`Survivors-only (nlme)` = c(
round(VBParams$F$Linf, 1), round(VBParams$M$Linf, 1),
round(VBParams$F$K, 3), round(VBParams$M$K, 3), "—"),
`Δ L∞ (%)` = c(
paste0(round(100*(Linf_F_j - VBParams$F$Linf)/VBParams$F$Linf, 1), "%"),
paste0(round(100*(Linf_M_j - VBParams$M$Linf)/VBParams$M$Linf, 1), "%"),
"—", "—", "—"),
check.names = FALSE
)
kable(growth_tbl, align = c("l","r","r","r","r"),
caption = "Von Bertalanffy growth parameters: joint model vs survivors-only estimate.")| Parameter | Joint model — mean | Joint model — 95% CI | Survivors-only (nlme) | Δ L∞ (%) |
|---|---|---|---|---|
| L∞_F (mm) | 659.400 | (657, 661.9) | 715.7 | -7.9% |
| L∞_M (mm) | 585.800 | (580.4, 591.1) | 718.8 | -18.5% |
| K_F (yr⁻¹) | 0.647 | (0.642, 0.65) | 0.209 | — |
| K_M (yr⁻¹) | 0.349 | (0.347, 0.35) | 0.119 | — |
| σ_obs (mm) | 28.500 | (27.3, 29.8) | — | — |
Correcting for survivor bias reduced asymptotic length by 7.9% for females and 18.5% for males, while increasing the growth coefficient approximately three-fold for both sexes. This is the expected signature of survivor bias: fish that die at small sizes are absent from the survivors-only dataset, making the asymptote appear larger and growth toward it appear slower than it truly is.
The sex-specific asymptotic sizes are now well-separated (L∞_F = 659 mm vs L∞_M = 586 mm), consistent with published reports of sexual size dimorphism in XYTE. The survivors-only nlme estimate showed near-identical L∞ for males and females — almost certainly an artifact of sparse recapture data for males (who have lower survival) combined with the survivor-bias inflation of L∞.
age_seq <- seq(0, 15, by = 0.1)
age_rel_mean <- mean(age_rel_cjs)
L0_mean <- mean(L0_cjs)
growth_df <- bind_rows(
data.frame(Age = age_seq, Sex = "Female",
TL = Linf_F_j - (Linf_F_j - L0_mean)*exp(-K_F_j * pmax(0, age_seq - age_rel_mean)),
Source = "Joint model"),
data.frame(Age = age_seq, Sex = "Male",
TL = Linf_M_j - (Linf_M_j - L0_mean)*exp(-K_M_j * pmax(0, age_seq - age_rel_mean)),
Source = "Joint model"),
data.frame(Age = age_seq, Sex = "Female",
TL = VBParams$F$Linf * (1 - exp(-VBParams$F$K * (age_seq - VBParams$F$t0))),
Source = "Survivors only"),
data.frame(Age = age_seq, Sex = "Male",
TL = VBParams$M$Linf * (1 - exp(-VBParams$M$K * (age_seq - VBParams$M$t0))),
Source = "Survivors only")
) |> dplyr::filter(TL > 100)
ggplot(growth_df, aes(x = Age, y = TL, color = Sex, linetype = Source)) +
geom_line(linewidth = 0.9) +
scale_color_manual(values = c(Female = "#c0392b", Male = "#2980b9")) +
scale_linetype_manual(values = c("Joint model" = "solid", "Survivors only" = "dashed")) +
labs(x = "Age (years from Feb 15 of year class)",
y = "Total length (mm)", color = NULL, linetype = NULL) +
theme_bw(base_size = 11) +
theme(legend.position = "bottom")
beta0_j <- param_summary["beta0", "mean"]
betaL_j <- param_summary["betaL", "mean"]The logit-scale survival slope is β_L = 1.665 per 50-mm increase in body length (95% CI: 1.588 – 1.744). Compared to the simple 180-day-lag logistic model (β_L ≈ 0.72–1.19 depending on cohort type), the joint model produced a steeper slope (+40–130%), consistent with the survivor bias having compressed the apparent TL–survival relationship in the simpler analysis.
surv_tbl_clean <- surv_tbl
rownames(surv_tbl_clean) <- NULL
kable(surv_tbl_clean, digits = 3, align = c("r","r","r","r","r"),
col.names = c("TL (mm)", "phi_true", "phi_eff", "95% CI low", "95% CI high"))| TL (mm) | phi_true | phi_eff | 95% CI low | 95% CI high |
|---|---|---|---|---|
| 300 | 0.000 | 0.000 | 0.000 | 0.001 |
| 350 | 0.002 | 0.002 | 0.001 | 0.003 |
| 400 | 0.011 | 0.010 | 0.008 | 0.014 |
| 450 | 0.055 | 0.052 | 0.046 | 0.064 |
| 500 | 0.233 | 0.221 | 0.214 | 0.252 |
| 550 | 0.616 | 0.585 | 0.598 | 0.635 |
| 600 | 0.894 | 0.850 | 0.883 | 0.905 |
tl_grid <- seq(250, 700, length.out = 200)
tl_scaled <- (tl_grid - TL_CENTER) / TL_SCALE
# Posterior draws for instantaneous phi_true
logit_S_draws <- sweep(outer(all_draws[,"betaL"], tl_scaled, "*"),
1, all_draws[,"beta0"], "+")
surv_ribbon <- data.frame(
TL = tl_grid,
S_med = apply(logit_S_draws, 2, \(x) median(plogis(x))),
S_lo = apply(logit_S_draws, 2, \(x) quantile(plogis(x), 0.025)),
S_hi = apply(logit_S_draws, 2, \(x) quantile(plogis(x), 0.975)),
Curve = "phi_true at current TL (projection matrix)"
)
# Effective first-year survival for a fish STOCKED at TL_release
# The Bayesian model evaluates phi at the fish's TL ~5 months after a Nov stocking
# (FY mid-point = April 1). Because the model corrects for p_det internally,
# this is TRUE biological survival and should exceed the logistic model's
# cumulative window-detection probability.
Linf_avg_sr <- 0.399*param_summary["Linf_F","mean"] + 0.601*param_summary["Linf_M","mean"]
K_avg_sr <- 0.399*param_summary["K_F","mean"] + 0.601*param_summary["K_M","mean"]
# TL at ~5 months post-release (FY mid-point for a Nov release)
tl_midyr <- Linf_avg_sr - (Linf_avg_sr - tl_grid) * exp(-K_avg_sr * 5/12)
logit_S_draws_eff <- sweep(
outer(all_draws[,"betaL"], (tl_midyr - 450)/50, "*"),
1, all_draws[,"beta0"], "+"
)
surv_eff <- data.frame(
TL = tl_grid,
S_med = apply(logit_S_draws_eff, 2, \(x) median(plogis(x))),
S_lo = apply(logit_S_draws_eff, 2, \(x) quantile(plogis(x), 0.025)),
S_hi = apply(logit_S_draws_eff, 2, \(x) quantile(plogis(x), 0.975)),
Curve = "Effective first-year phi for fish released at TL (accounting for growth)"
)
# Logistic model post-stocking cumulative detection (empirical reference)
b0r <- lag_results$DAL180$fit_rel_base$coefficients["(Intercept)"]
bLr <- lag_results$DAL180$fit_rel_base$coefficients["ReleaseTL"]
logistic_ref <- data.frame(
TL = tl_grid,
S = plogis(b0r + bLr * tl_grid),
Curve = "Logistic: P(detected in window | stocked at TL)"
)library(dplyr)
all_surv <- bind_rows(
surv_ribbon,
surv_eff
)
ggplot() +
geom_ribbon(data = all_surv |> filter(Curve == "phi_true at current TL (projection matrix)"),
aes(x=TL, ymin=S_lo, ymax=S_hi), fill="#2166ac", alpha=0.15) +
geom_line(data = all_surv,
aes(x=TL, y=S_med, color=Curve, linetype=Curve), linewidth=0.95) +
geom_line(data = logistic_ref,
aes(x=TL, y=S, color=Curve, linetype=Curve), linewidth=0.8) +
scale_color_manual(values = c(
"phi_true at current TL (projection matrix)" = "#2166ac",
"Effective first-year phi for fish released at TL (accounting for growth)" = "#2166ac",
"Logistic: P(detected in window | stocked at TL)" = "#d73027"
)) +
scale_linetype_manual(values = c(
"phi_true at current TL (projection matrix)" = "solid",
"Effective first-year phi for fish released at TL (accounting for growth)" = "dashed",
"Logistic: P(detected in window | stocked at TL)" = "dotted"
)) +
scale_y_continuous(limits = c(0, 1), labels = percent) +
labs(x = "Total length (mm)", y = "Survival / detection probability",
color = NULL, linetype = NULL) +
theme_bw(base_size = 10) +
theme(legend.position = "bottom",
legend.text = element_text(size = 8))
p_det_mean <- param_summary["p_det", "mean"]
p_det_lo <- param_summary["p_det", "2.5%"]
p_det_hi <- param_summary["p_det", "97.5%"]
prop_F_mean <- param_summary["prop_F", "mean"]Annual Basin-zone detection probability was estimated at \(p_{\text{det}}\) = 0.716 (95% CrI: 0.708 – 0.723). This is somewhat lower than the Stalwart multi-state model’s per-FY marginal detection (≈ 0.88–0.90), which was conditioned on fish known to be alive and resident in the Basin. The broader population here includes younger, smaller fish with reduced site fidelity and a wider range of spatial use, consistent with a lower reach-average detection probability.
The proportion female among all Basin releases was estimated at 0.399 (95% CrI: 0.382 – 0.417), slightly below the 0.60 confirmed-sex ratio among recaptured fish — suggesting the confirmed-sex sample is modestly biased toward females (expected, given female advantage in survival and recapture probability).
The joint model provides the key deliverable for size-structured projection modeling: a mapping from age × sex to expected size and from size to survival, integrated into a single age- and sex-specific survival curve.
For a fish of sex \(s\) released at age \(a_{\text{rel}}\) with TL \(L_0\), annual survival at age \(a\) is obtained by plugging the expected TL into the survival model:
\[\phi_{\text{true}}(a, s) = \text{logit}^{-1}\!\left(\hat\beta_0 + \hat\beta_L \cdot \frac{\hat L_{\infty,s} - (\hat L_{\infty,s} - L_0)\,e^{-\hat K_s(a - a_{\text{rel}})} - 450}{50}\right)\]
age_eval <- seq(3, 15, by = 0.1)
L0_rep <- 400 # representative release TL (mm)
age_rel_rep <- 3 # representative age at release
# Posterior draws for each sex
sex_surv_draws <- lapply(c(Female = "F", Male = "M"), function(sx) {
Linf_col <- if (sx == "F") "Linf_F" else "Linf_M"
K_col <- if (sx == "F") "K_F" else "K_M"
sapply(age_eval, function(a) {
L_a <- all_draws[, Linf_col] -
(all_draws[, Linf_col] - L0_rep) *
exp(-all_draws[, K_col] * pmax(0, a - age_rel_rep))
logit_S <- all_draws[,"beta0"] + all_draws[,"betaL"] * (L_a - 450) / 50
plogis(logit_S)
})
})
surv_age_df <- bind_rows(lapply(names(sex_surv_draws), function(nm) {
mat <- sex_surv_draws[[nm]]
data.frame(
Age = age_eval,
Sex = nm,
S_med = apply(mat, 2, median),
S_lo = apply(mat, 2, \(x) quantile(x, 0.025)),
S_hi = apply(mat, 2, \(x) quantile(x, 0.975))
)
}))
ggplot(surv_age_df, aes(x = Age, color = Sex, fill = Sex)) +
geom_ribbon(aes(ymin = S_lo, ymax = S_hi), alpha = 0.18, colour = NA) +
geom_line(aes(y = S_med), linewidth = 1) +
scale_color_manual(values = c(Female = "#c0392b", Male = "#2980b9")) +
scale_fill_manual(values = c(Female = "#c0392b", Male = "#2980b9")) +
scale_y_continuous(limits = c(0, 1), labels = percent) +
labs(x = "Age (years from year class)", y = "Annual biological survival",
color = NULL, fill = NULL,
caption = "Representative release: age 3, TL = 400 mm. Bands: 95% posterior CrI.") +
theme_bw(base_size = 11) +
theme(legend.position = "bottom")
Females achieve higher annual survival at every age because females grow larger (higher L∞ and higher K), and larger body size increases survival. The advantage widens during the rapid growth phase (ages 3–8) when the size difference between sexes is largest.
# Build grid of age × sex × posterior mean parameters
traj_df <- bind_rows(lapply(c("Female","Male"), function(sx) {
Linf_use <- if (sx == "Female") Linf_F_j else Linf_M_j
K_use <- if (sx == "Female") K_F_j else K_M_j
L_at_age <- Linf_use - (Linf_use - L0_rep)*exp(-K_use * pmax(0, age_eval - age_rel_rep))
S_at_age <- plogis(beta0_j + betaL_j * (L_at_age - 450) / 50)
data.frame(Age = age_eval, Sex = sx, TL = L_at_age, Survival = S_at_age)
}))
ggplot(traj_df, aes(x = Age, y = TL, color = Survival, group = Sex)) +
geom_line(linewidth = 2) +
geom_text(data = traj_df |> dplyr::filter(Age %in% c(5, 10)),
aes(label = Sex), vjust = -0.6, size = 3, color = "grey30") +
scale_color_gradientn(
colors = c("#d73027","#fdae61","#ffffbf","#a6d96a","#1a9641"),
limits = c(0, 1), labels = percent, name = "Annual survival"
) +
labs(x = "Age (years from year class)", y = "Total length (mm)",
caption = "Representative release: age 3, TL = 400 mm. Colour = annual biological survival.") +
theme_bw(base_size = 11)
For a size-structured (Leslie-type) projection matrix with 50-mm size classes \(\{[200,250),\,[250,300),\,\ldots\}\), the transition element \(A_{kj}\) (contribution from size class \(j\) to size class \(k\) in one year) is:
\[A_{kj} = S(m_j) \cdot G_{kj}\]
where \(m_j\) is the midpoint of size class \(j\), \(S(m_j)\) is the joint-model survival probability, and \(G_{kj}\) is the probability that a fish in class \(j\) grows into class \(k\) in one year, derived from the sex-specific VBGF.
The survival vector uses the joint-model parameters directly:
\[S(m_j) = \text{logit}^{-1}(\hat\beta_0 + \hat\beta_L \cdot (m_j - 450)/50)\]
The growth transition matrix requires a distributional assumption for individual growth variation around the VBGF expectation (e.g., Normal with SD informed by \(\sigma_{\text{obs}}\)):
\[G_{kj} = \Pr\!\left(m_k - 25 \leq L_{t+1} \leq m_k + 25 \;\middle|\; L_t = m_j,\, \hat L_\infty,\, \hat K\right)\]
# Extracting joint-model parameters for projection use
joint_params <- list(
# Growth (sex-averaged for a sex-aggregated matrix; use sex-specific for two-sex model)
Linf_F = param_summary["Linf_F", "mean"],
Linf_M = param_summary["Linf_M", "mean"],
K_F = param_summary["K_F", "mean"],
K_M = param_summary["K_M", "mean"],
sigma_g = param_summary["sigma_obs", "mean"], # individual growth variation (approx)
prop_F = param_summary["prop_F", "mean"],
# Survival
beta0 = param_summary["beta0", "mean"],
betaL = param_summary["betaL", "mean"],
TL_ctr = 450, # centering constant (mm)
TL_scl = 50, # scaling constant (mm per unit)
theta = 0.05 # tag-loss rate
)
# ── Survival vector for 50-mm size classes 150–750 mm ────────────────────────
size_breaks <- seq(150, 750, by = 50)
size_midpts <- size_breaks[-1] - 25 # midpoints: 175, 225, ..., 725
surv_vec <- plogis(joint_params$beta0 +
joint_params$betaL * (size_midpts - joint_params$TL_ctr) /
joint_params$TL_scl)
names(surv_vec) <- paste0(size_midpts, "mm")
# ── Growth transition matrix (sex-averaged, single time step) ─────────────────
# Proportion female (from joint model)
pF <- joint_params$prop_F
# Sex-averaged VBGF increment at each size class midpoint
K_avg <- pF * joint_params$K_F + (1 - pF) * joint_params$K_M
Linf_avg <- pF * joint_params$Linf_F + (1 - pF) * joint_params$Linf_M
# Expected growth increment from size m in one year
growth_increment <- function(m) (Linf_avg - m) * (1 - exp(-K_avg))
G_mat <- matrix(0, nrow = length(size_midpts), ncol = length(size_midpts))
for (j in seq_along(size_midpts)) {
mu_next <- size_midpts[j] + growth_increment(size_midpts[j])
# P(land in each class k) ~ Normal(mu_next, sigma_g)
for (k in seq_along(size_midpts)) {
G_mat[k, j] <- pnorm(size_breaks[k+1], mu_next, joint_params$sigma_g) -
pnorm(size_breaks[k], mu_next, joint_params$sigma_g)
}
}
# ── Leslie matrix A = diag(S) %*% G ──────────────────────────────────────────
A_mat <- diag(surv_vec) %*% G_mat
cat("Projection matrix dimensions:", dim(A_mat), "\n")Projection matrix dimensions: 12 12 cat("Dominant eigenvalue (lambda):", round(Re(eigen(A_mat)$values[1]), 4), "\n")Dominant eigenvalue (lambda): 0.9051 cat("Size class with highest survival:", size_midpts[which.max(surv_vec)], "mm\n")Size class with highest survival: 725 mmThe dominant eigenvalue \(\lambda\) from the Leslie matrix gives the long-run population growth rate under current survival and growth conditions. The stable size distribution (right eigenvector) and reproductive values (left eigenvector) follow from the same matrix. Note that this matrix does not yet include recruitment; a complete model requires adding first-year survival and age/size at first reproduction.
To propagate MCMC uncertainty into derived projection quantities, draw directly from the all_samples MCMC object and build the projection matrix for each draw:
# Example: distribution of lambda across MCMC draws
n_draws <- nrow(all_draws)
lambdas <- numeric(n_draws)
for (d in seq_len(n_draws)) {
b0d <- all_draws[d, "beta0"]; bLd <- all_draws[d, "betaL"]
Lfd <- all_draws[d, "Linf_F"]; Lmd <- all_draws[d, "Linf_M"]
Kfd <- all_draws[d, "K_F"]; Kmd <- all_draws[d, "K_M"]
pFd <- all_draws[d, "prop_F"]
S_d <- plogis(b0d + bLd * (size_midpts - 450) / 50)
K_d <- pFd * Kfd + (1-pFd) * Kmd
Li_d <- pFd * Lfd + (1-pFd) * Lmd
G_d <- matrix(0, length(size_midpts), length(size_midpts))
for (j in seq_along(size_midpts)) {
mu_j <- size_midpts[j] + (Li_d - size_midpts[j])*(1-exp(-K_d))
for (k in seq_along(size_midpts))
G_d[k,j] <- pnorm(size_breaks[k+1], mu_j, 28.5) - pnorm(size_breaks[k], mu_j, 28.5)
}
A_d <- diag(S_d) %*% G_d
lambdas[d] <- Re(eigen(A_d)$values[1])
}
cat("Posterior λ: mean =", round(mean(lambdas),4),
" 95% CrI: [", round(quantile(lambdas,.025),4), ",",
round(quantile(lambdas,.975),4), "]\n")param_summary |>
as.data.frame() |>
tibble::rownames_to_column("Parameter") |>
mutate(across(where(is.numeric), \(x) round(x, 4))) |>
kable(align = c("l", rep("r", ncol(param_summary))))| Parameter | mean | sd | 2.5% | 50% | 97.5% | Rhat | n.eff |
|---|---|---|---|---|---|---|---|
| K_F | 0.6473 | 0.0020 | 0.6421 | 0.6478 | 0.6497 | 1 | 3015 |
| K_M | 0.3492 | 0.0008 | 0.3471 | 0.3494 | 0.3500 | 1 | 8412 |
| K_diff | 0.2981 | 0.0018 | 0.2931 | 0.2986 | 0.2999 | 1 | 2853 |
| Linf_F | 659.4047 | 1.2506 | 656.9644 | 659.3972 | 661.8580 | 1 | 5193 |
| Linf_M | 585.8319 | 2.7164 | 580.3939 | 585.8691 | 591.0802 | 1 | 2649 |
| beta0 | -2.8566 | 0.0870 | -3.0280 | -2.8553 | -2.6880 | 1 | 1877 |
| betaL | 1.6647 | 0.0402 | 1.5877 | 1.6641 | 1.7437 | 1 | 2393 |
| beta_post | -0.0048 | 0.0048 | -0.0177 | -0.0033 | -0.0001 | 1 | 7898 |
| p_det | 0.7155 | 0.0036 | 0.7083 | 0.7155 | 0.7226 | 1 | 15258 |
| prop_F | 0.3992 | 0.0091 | 0.3815 | 0.3992 | 0.4172 | 1 | 6038 |
| sigma_obs | 28.4862 | 0.6324 | 27.2638 | 28.4770 | 29.7607 | 1 | 7075 |
Key findings from the joint model:
Survivor bias is substantial. Conditioning on non-survivors reduces estimated asymptotic length by 8–19% and increases growth coefficients approximately three-fold relative to survivors-only estimation. Males show a larger bias because their lower survival rate means fewer survive long enough to be measured repeatedly.
Sex-specific growth curves are now well-separated. The joint model recovers L∞_F − L∞_M ≈ 74 mm, consistent with published accounts of sexual size dimorphism in XYTE. The survivors-only approach gave near-identical asymptotes for both sexes.
The TL–survival relationship is steeper than simple logistic estimates. β_L = 1.665 per 50 mm from the joint model vs 0.72–1.19 from the logistic-only model. This steepening is mechanistically coherent: if slow-growing fish were not included in the simple logistic’s data, the apparent slope was compressed.
Annual detection probability for the full release population is ~0.72, noticeably below the Stalwart-model estimate (~0.88) for long-term confirmed Basin residents. This difference is biologically sensible given younger/smaller fish in the broader population are less site-fidelitous.
Growth and survival are now jointly estimable for projection. The Leslie-matrix template above requires only the posterior means (or, preferably, MCMC draws for full uncertainty propagation) from this model.