Selectivity Options for Age-Structured Stock Assessments

A comparative framework using RTMB

Author
Affiliation

James Ianelli

NOAA/NMFS Alaska Fisheries Science Center

Published

April 7, 2026

Abstract

Selectivity is a fundamental component of age-structured stock assessment models, governing how fishing mortality and survey catchability vary with age or size. The choice of selectivity parameterization affects estimates of population abundance, fishing mortality, and management reference points. This manuscript compares several selectivity formulations—including standard logistic, double logistic, spline-based, and time-varying 2D autoregressive approaches—within a unified RTMB framework. We examine parameter identifiability, prior sensitivity, and the practical consequences of selectivity misspecification using simulation and application to eastern Bering Sea walleye pollock (Gadus chalcogrammus). Results highlight trade-offs between flexibility and estimability and provide guidance for practitioners choosing among selectivity options in operational assessments.

Keywords

selectivity, stock assessment, RTMB, fisheries, walleye pollock

1 Introduction

Selectivity functions describe the relative vulnerability of fish to capture as a function of age (or size) and are among the most influential—yet least observable—components of stock assessment models. The assumed shape of the selectivity curve directly affects estimates of spawning biomass, recruitment, and fishing mortality, and misspecification can propagate into biased reference points and management advice (Punt, Hurtado-Ferro, and Whitten 2013; Thompson 1994).

Despite its importance, selectivity is often treated as a modelling convenience rather than a biological or technological quantity to be estimated carefully. Most operational assessments adopt a logistic or double-logistic functional form, chosen for parsimony rather than fidelity to the underlying catch process. More flexible alternatives—such as penalized splines or non-parametric approaches—can reduce structural bias but may introduce identifiability problems, particularly when data are sparse or conflicting.

This manuscript develops a comparative framework for evaluating selectivity options within the R Template Model Builder (RTMB) environment. RTMB provides automatic differentiation, Laplace approximation for random effects, and integration with Bayesian sampling via SparseNUTS, making it a natural platform for exploring both maximum likelihood and posterior-based inference on selectivity parameters.

We organize the comparison around five selectivity families:

  1. Standard logistic (Section): the two-parameter ascending logistic, widely used for fisheries where retention is assumed to be monotonically increasing with age.
  2. Double logistic (Section): a dome-shaped three-parameter formulation allowing selectivity to decline at older ages, motivated by gear avoidance, ontogenetic habitat shifts, or differential availability.
  3. Generalized logistic (Richards) (Section): a flexible extension of logistic selectivity that can represent asymmetric transitions and plateau-like behavior while remaining smooth and differentiable.
  4. Spline-based (Section): a penalized B-spline approach offering flexible, data-driven selectivity shapes with smoothness controlled by a penalty parameter.
  5. Time-varying 2D AR1 (Section): a separable autoregressive structure over year and age dimensions, allowing selectivity to evolve smoothly through time while maintaining age-specific coherence via a Kronecker-structured precision matrix.

For each family, we present the mathematical formulation, an RTMB implementation with prior specification, parameter sensitivity analysis, and MCMC-based posterior evaluation. We then apply all to eastern Bering Sea walleye pollock data to illustrate practical trade-offs.

1.1 Review Scope and Literature Context

This review focuses on methodological literature that is directly relevant to three practical assessment decisions: (i) choice of selectivity functional form, (ii) model selection among competing selectivity hypotheses, and (iii) construction of time-varying selectivity processes with explicit regularization.

1.1.1 Selectivity form choice in integrated assessments

Operationally, selectivity in integrated assessments is best interpreted as the combined effect of gear contact/retention and fish availability, not as a purely mechanistic gear curve (Punt, Hurtado-Ferro, and Whitten 2013; Privitera-Johnson, Methot, and Punt 2022). In that setting, the choice of selectivity parameterization is consequential because it changes inference on spawning biomass, recruitment, and fishing mortality.

Selectivity-form selection therefore requires both statistical and biological diagnostics, including residual pattern checks, sensitivity analyses, and information-based comparisons among plausible alternatives (Punt, Hurtado-Ferro, and Whitten 2013). A key warning is that non-monotone selectivity can be strongly confounded with natural mortality, particularly when older ages are weakly informed by composition data (Thompson 1994).

1.1.2 Gear-selectivity estimation foundations

Classical selectivity-estimation literature remains important for assessment modeling because it clarifies what selectivity data can identify and where latent availability effects remain unresolved (Millar and Fryer 1999; Wileman et al. 1996). This evidence base supports explicit regularization whenever flexible age-specific selectivity surfaces are estimated.

1.1.3 Time-varying selectivity literature

Simulation studies show that time-varying selectivity can materially improve or degrade assessment performance depending on how process flexibility matches the data-generating mechanism (Linton and Bence 2011). Practitioner guidance also emphasizes that subjective choices such as block boundaries, smoothing penalties, and process assumptions can dominate outcomes (Martell and Stewart 2014).

Recent semi-parametric approaches formalize autocorrelation across both age and time, motivating separable latent-process structures that can be estimated within modern AD frameworks (Xu et al. 2019). A foundational applied example is the likelihood-based Southern Bluefin Tuna assessment of (Butterworth, Ianelli, and Hilborn 2003). This method allows temporally structured selectivity changes and treats catch-at-age data as error-prone rather than fixed. Broader SCA model-selection and time-varying process work reinforces the same bias-variance trade-off (Wilberg and Bence 2006, 2008).

1.2 Statistical Framing as a Latent Process

Although selectivity is often presented as a deterministic curve, it is more usefully viewed as a latent process that allocates fishing intensity over age and time. Under this framing, selecting a functional form is equivalent to selecting a prior precision structure on an unobserved selectivity surface.

Low-dimensional parametric forms (for example, logistic and double logistic) impose strong structural priors and are robust when data are sparse. Flexible forms (for example, spline and time-varying AR structures) reduce structural bias risk, but require explicit regularization to prevent confounding and overfitting (Martell and Stewart 2014; Xu et al. 2019).

This perspective motivates three principles for applied assessments:

  1. Match effective selectivity complexity to data information content.
  2. Treat regularization (penalties/priors/precision structures) as essential, not optional.
  3. Use model-selection tools (AIC, DIC, cross-validation, posterior predictive checks) to evaluate bias-variance trade-offs under explicit prior structure (Punt, Hurtado-Ferro, and Whitten 2013; Wilberg and Bence 2008).

Regularization can also be framed as numerical stabilization: penalties such as ridge/Tikhonov terms or spline roughness controls impose smoothness that improves conditioning and mitigates overfitting, echoing the numerical recipes literature on stabilization and controlled flexibility in nonlinear optimization (Press et al. 1992).

Within this unifying view, logistic, dome-shaped, spline, blocked, and AR(1) approaches differ primarily in the precision matrix implied for latent selectivity deviations, rather than only in curve geometry.

Standard Logistic Selectivity

2 Standard Logistic Selectivity

2.1 Overview

The standard (ascending) logistic selectivity is the simplest and most commonly used parameterization in age-structured stock assessments. It assumes that selectivity increases monotonically with age and asymptotes at full selectivity. Its main advantage is parsimony: a small number of parameters can often stabilize inference when composition data are limited [Methot and Wetzel (2013); Punt, Hurtado-Ferro, and Whitten (2013)].

As a consequence, logistic selectivity is frequently used as a baseline candidate in model-comparison workflows. However, when older fish exhibit reduced vulnerability, strict monotonicity can induce structural bias and shift variation into other parameters [Thompson (1994)].

2.2 Model definition

With parameters \(a_{50}\) (age at 50% selectivity) and \(\delta\) (slope width from 50% to 95% selectivity):

\[ \text{sel}(a) = \frac{1}{1 + \exp\left[-\log(19)\,\frac{a - a_{50}}{\delta}\right]} \]

where \(\delta > 0\) controls how quickly selectivity transitions from low to high values.

2.3 Scenario exploration

2.4 RTMB implementation

Placeholder: RTMB objective function with priors on \(a_{50}\) and \(\delta\), optimization, and MCMC sampling.

Double Logistic Selectivity

2.5 Overview

This notebook documents the 3-parameter double logistic selectivity curve used for prior elicitation in the EBS pollock assessment workflow. The formulation allows dome-shaped selectivity, where vulnerability increases with age to a peak and then declines. This behavior is motivated by gear avoidance at older ages, ontogenetic habitat shifts, or differential availability between age classes.

Within integrated assessments, double-logistic forms are often a practical extension of monotone logistic selectivity because they can represent declining vulnerability at older ages while remaining low-dimensional (Punt, Hurtado-Ferro, and Whitten 2013). Interpretation still requires care: dome shape can be confounded with natural mortality and cohort processes if independent information is weak (Thompson 1994).

2.6 Model definition

Let age be \(a\), with parameters \(p_1\), \(p_2\), and \(p_3\) and derived inflection points \(\gamma_1\) and \(\gamma_2\):

\[ \gamma_1 = p_1 + p_2, \qquad \gamma_2 = 2p_1 + p_2 + p_3. \]

The ascending limb is a standard logistic and the descending limb is one minus a logistic:

\[ \text{asc}(a) = \frac{1}{1 + \exp\left[-\log(19)\,\frac{a - \gamma_1}{p_1}\right]}, \]

\[ \text{desc}(a) = 1 - \frac{1}{1 + \exp\left[-\log(19)\,\frac{a - \gamma_2}{p_3}\right]}. \]

The full double logistic selectivity is

\[ \text{sel}(a) = \min\left(1,\ \frac{\text{asc}(a)\,\text{desc}(a)}{0.95^2}\right). \]

2.6.1 Interpretation

  • \(p_1 > 0\) controls the ascending slope and is the distance from \(\gamma_1\) (50% selectivity) to the 95% point.
  • \(p_2\) shifts the ascending limb through \(\gamma_1 = p_1 + p_2\).
  • \(p_3 > 0\) controls the descending slope and is the distance from \(\gamma_2\) (50% on the descending limb) to the 5% point.
  • The dome width is \(\gamma_2 - \gamma_1 = p_1 + p_3\).
  • The factor \(0.95^{-2}\) normalizes the curve so that when both limbs are at 0.95 the product is near 1 (then capped at 1.0).

2.7 Scenarios

Show R code 
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)

dbl_logistic_sel <- function(age, p1, p2, p3) {
  gamma1 <- p1 + p2
  gamma2 <- 2 * p1 + p2 + p3

  asc <- 1 / (1 + exp(-log(19) * (age - gamma1) / p1))
  desc <- 1 - 1 / (1 + exp(-log(19) * (age - gamma2) / p3))

  sel <- asc * desc * 0.95^(-2)
  pmin(sel, 1.0)
}

ages <- seq(1, 15, by = 0.1)

scenarios <- tibble(
  scenario = c(
    "Wide dome",
    "Asymmetric\n(fast rise, slow decline)",
    "Asymmetric\n(slow rise, fast decline)",
    "Nearly asymptotic"
  ),
  p1 = c(2.0, 0.8, 2.0, 1.5),
  p2 = c(4.0, 4.0, 3.0, 4.0),
  p3 = c(2.5, 3.0, 0.8, 8.0)
) %>%
  mutate(
    gamma1 = p1 + p2,
    gamma2 = 2 * p1 + p2 + p3,
    dome_width = gamma2 - gamma1
  )

sel_data <- scenarios %>%
  rowwise() %>%
  do({
    scen <- .
    tibble(
      scenario = scen$scenario,
      age = ages,
      selectivity = dbl_logistic_sel(ages, scen$p1, scen$p2, scen$p3),
      p1 = scen$p1,
      p2 = scen$p2,
      p3 = scen$p3,
      gamma1 = scen$gamma1,
      gamma2 = scen$gamma2
    )
  }) %>%
  ungroup() %>%
  mutate(
    label = sprintf("%s\n(p1=%.1f, p2=%.1f, p3=%.1f)", scenario, p1, p2, p3)
  )
Show R code 
knitr::kable(
  scenarios,
  digits = 2
)
Table 1: Scenario parameters and derived inflection points
scenario p1 p2 p3 gamma1 gamma2 dome_width
Wide dome 2.0 4 2.5 6.0 10.5 4.5
Asymmetric
(fast rise, slow decline) 0.8 4 3.0 4.8 8.6 3.8
Asymmetric
(slow rise, fast decline) 2.0 3 0.8 5.0 7.8 2.8
Nearly asymptotic 1.5 4 8.0 5.5 15.0 9.5

The scenario summary in Table lists the four parameter sets and the derived inflection points (\(\gamma_1\), \(\gamma_2\)) used in the selectivity scenarios.

2.8 Faceted scenario curves

Each scenario appears on its own panel in Figure so the ascent, descent, and dome width can be compared against the 50% and 95% reference lines.

Show R code 
p1_plot <- ggplot(sel_data, aes(x = age, y = selectivity)) +
  geom_line(linewidth = 1.2, color = "steelblue") +
  geom_hline(yintercept = c(0.5, 0.95), linetype = "dashed",
             color = "gray50", alpha = 0.7) +
  geom_vline(aes(xintercept = gamma1), linetype = "dotted",
             color = "red", alpha = 0.5) +
  geom_vline(aes(xintercept = gamma2), linetype = "dotted",
             color = "orange", alpha = 0.5) +
  facet_wrap(~label, ncol = 2) +
  scale_x_continuous(breaks = seq(2, 14, by = 2)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = seq(0, 1, by = 0.25)) +
  labs(
    title = "3-parameter double logistic selectivity curves",
    subtitle = "Dashed lines show 50% and 95% selectivity",
    x = "Age",
    y = "Selectivity",
    caption = "Red dotted = gamma1, orange dotted = gamma2"
  ) +
  theme_bw(base_size = 11) +
  theme(
    strip.text = element_text(size = 9),
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

p1_plot
Figure 1: Double logistic selectivity for each parameter scenario.

2.9 Overlay comparison

All scenarios are overlaid in Figure to highlight differences in peak age and decline shape when the curves are viewed on a common axis.

Show R code 
p2_plot <- ggplot(sel_data, aes(x = age, y = selectivity, color = scenario)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = c(0.5, 0.95), linetype = "dashed",
             color = "gray50", alpha = 0.5) +
  scale_x_continuous(breaks = seq(2, 14, by = 2)) +
  scale_y_continuous(limits = c(0, 1.05), breaks = seq(0, 1, by = 0.25)) +
  scale_color_brewer(palette = "Set1", name = "Scenario") +
  labs(
    title = "Comparison of double logistic selectivity curves",
    subtitle = "All scenarios overlaid",
    x = "Age",
    y = "Selectivity"
  ) +
  theme_bw(base_size = 12) +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

p2_plot
Figure 2: All scenarios overlaid for direct comparison.

2.10 Parameter sensitivity analysis

The three panels in Figure show how each parameter changes the curve when the other two are held fixed.

Show R code 
p1_vary <- expand_grid(
  p1 = seq(0.5, 3, by = 0.5),
  p2 = 5,
  p3 = 2,
  age = ages
) %>%
  rowwise() %>%
  mutate(selectivity = dbl_logistic_sel(age, p1, p2, p3)) %>%
  ungroup()

p3a <- ggplot(p1_vary, aes(x = age, y = selectivity, color = factor(p1))) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d(name = "p1", option = "C") +
  scale_y_continuous(limits = c(0, 1.05)) +
  labs(title = "Sensitivity to p1 (ascending slope)",
       subtitle = "p2 = 5, p3 = 2 fixed",
       x = "Age", y = "Selectivity") +
  theme_bw(base_size = 9) +
  theme(legend.position = "right", aspect.ratio = 1)

p2_vary <- expand_grid(
  p1 = 1.5,
  p2 = seq(2, 8, by = 1),
  p3 = 2,
  age = ages
) %>%
  rowwise() %>%
  mutate(selectivity = dbl_logistic_sel(age, p1, p2, p3)) %>%
  ungroup()

p3b <- ggplot(p2_vary, aes(x = age, y = selectivity, color = factor(p2))) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d(name = "p2", option = "D") +
  scale_y_continuous(limits = c(0, 1.05)) +
  labs(title = "Sensitivity to p2 (horizontal shift)",
       subtitle = "p1 = 1.5, p3 = 2 fixed",
       x = "Age", y = "Selectivity") +
  theme_bw(base_size = 9) +
  theme(legend.position = "right", aspect.ratio = 1)

p3_vary <- expand_grid(
  p1 = 1.5,
  p2 = 4,
  p3 = seq(0.5, 6.5, by = 1),
  age = ages
) %>%
  rowwise() %>%
  mutate(selectivity = dbl_logistic_sel(age, p1, p2, p3)) %>%
  ungroup()

p3c <- ggplot(p3_vary, aes(x = age, y = selectivity, color = factor(p3))) +
  geom_line(linewidth = 1) +
  scale_color_viridis_d(name = "p3", option = "E") +
  scale_y_continuous(limits = c(0, 1.05)) +
  labs(title = "Sensitivity to p3 (descending slope)",
       subtitle = "p1 = 1.5, p2 = 4 fixed",
       x = "Age", y = "Selectivity") +
  theme_bw(base_size = 9) +
  theme(legend.position = "right", aspect.ratio = 1)

p_sensitivity <- (p3a | p3b) / (p3c | plot_spacer()) +
  plot_annotation(
    title = "Parameter sensitivity analysis",
    subtitle = "Effect of varying each parameter individually"
  )

p_sensitivity
Figure 3: Sensitivity to p1, p2, and p3 with other parameters fixed.

2.11 RTMB implementation with priors

A common approach is to estimate on the log scale for the positive parameters and apply lognormal priors on \(p_1\) and \(p_3\):

\[ \log(p_k) \sim \mathcal{N}(\mu_k, \sigma_k^2), \qquad k \in \{1, 3\}. \]

If using a lognormal prior specified by median \(m\) and coefficient of variation \(\text{CV}\) on the natural scale, a convenient conversion is

\[ \mu = \log(m), \qquad \sigma = \sqrt{\log(\text{CV}^2 + 1)}. \]

Parameter \(p_2\) can remain unconstrained (normal prior) or be modeled on the log scale if you want to enforce \(p_2 > 0\).

outer mgc:  0 
outer mgc:  0 
outer mgc:  0.006737636 
outer mgc:  0.006737636 
outer mgc:  0.001 
outer mgc:  0.001 
outer mgc:  0.003252194 
outer mgc:  0.003252194 
outer mgc:  2.5

2.12 MCMC evaluation



Gradient evaluation took 4.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.018 seconds (Warm-up)
               0.208 seconds (Sampling)
               0.226 seconds (Total)



Model 'selectivity_set1' has 3 pars, and was fit using NUTS with a 'diag' metric
1 chain(s) of 1300 total iterations (150 warmup) were used
Average run time per chain was 0.23 seconds 
Minimum ESS=360.7 (31.37%), and maximum Rhat=1.015
There were 0 divergences after warmup
Show R code 
SparseNUTS::plot_marginals(fit_set1, pars = 1:3)
Figure 4: Marginal posterior distributions for selectivity parameters (SparseNUTS, prior set 1).

The marginal posterior densities in Figure summarize \((p_1, p_2, p_3)\) under prior set 1.

Show R code 
pairs(fit_set1)
Figure 5: Pairwise posterior scatterplots for selectivity parameters.
Show R code 
df_set1 <- mcmc_selectivity_df(fit_set1, age)
plot_mcmc_selectivity(df_set1)
Figure 6: MCMC selectivity draws for prior set 1 (p1 median 1.5, p3 median 2.5, p2 sd 1).

2.12.1 Additional prior sets



Gradient evaluation took 2.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.018 seconds (Warm-up)
               0.122 seconds (Sampling)
               0.14 seconds (Total)



Model 'selectivity_set2' has 3 pars, and was fit using NUTS with a 'diag' metric
1 chain(s) of 1300 total iterations (150 warmup) were used
Average run time per chain was 0.14 seconds 
Minimum ESS=393.8 (34.24%), and maximum Rhat=1.003
There were 0 divergences after warmup


Gradient evaluation took 2.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.017 seconds (Warm-up)
               0.129 seconds (Sampling)
               0.146 seconds (Total)



Model 'selectivity_set3' has 3 pars, and was fit using NUTS with a 'diag' metric
1 chain(s) of 1300 total iterations (150 warmup) were used
Average run time per chain was 0.15 seconds 
Minimum ESS=427 (37.13%), and maximum Rhat=1.005
There were 0 divergences after warmup
Show R code 
df_set2 <- mcmc_selectivity_df(fit_set2, age)
plot_mcmc_selectivity(df_set2)
Figure 7: MCMC selectivity draws for prior set 2 (p1 median 1.5, p3 median 8, p2 sd 1).
Show R code 
df_set3 <- mcmc_selectivity_df(fit_set3, age)
plot_mcmc_selectivity(df_set3)
Figure 8: MCMC selectivity draws for prior set 3 (p1 median 1.5, p3 median 8, p2 sd 2).

2.13 Generalized Logistic (Richards)

The generalized logistic (Richards) form is a practical extension of the standard and double logistic families when additional curvature is required, particularly for asymmetric transitions and near-plateau behavior. The interactive tool linked from Maunder’s selectivity page (Fitting Double Richards to Double Normal Selectivity) shows that a double generalized logistic can closely approximate a split-normal-with-plateau target while preserving smooth derivatives useful for optimization.

2.13.1 Model definition

For an ascending Richards curve at age \(a\):

\[ \text{Richards}_{\uparrow}(a;\,a_{50},k,\nu)= \left(1+\exp\left[-k(a-a_{50})\right]\right)^{-1/\nu}, \]

where \(a_{50}\) is an inflection-location parameter, \(k>0\) controls transition rate, and \(\nu>0\) controls asymmetry (shape).

A descending limb can be written as:

\[ \text{Richards}_{\downarrow}(a;\,a_{50},k,\nu)= 1-\left(1+\exp\left[-k(a-a_{50})\right]\right)^{-1/\nu}. \]

A convenient double-Richards selectivity is then:

\[ \text{sel}(a)= \frac{ \text{Richards}_{\uparrow}(a;\,a_{50,1},k_1,\nu_1)\, \text{Richards}_{\downarrow}(a;\,a_{50,2},k_2,\nu_2) }{ \max_a\left[ \text{Richards}_{\uparrow}(a;\,a_{50,1},k_1,\nu_1)\, \text{Richards}_{\downarrow}(a;\,a_{50,2},k_2,\nu_2) \right] }, \]

which scales peak selectivity to 1.0.

2.13.2 Optional parameter settings

Show R code 
richards_scenarios <- tibble::tibble(
  scenario = c(
    "Near-symmetric dome",
    "Fast rise, slow decline",
    "Slow rise, fast decline",
    "Broad plateau-like dome"
  ),
  a50_up = c(4.5, 4.0, 5.2, 4.8),
  k_up = c(1.2, 2.0, 0.8, 1.1),
  nu_up = c(1.0, 0.7, 1.4, 0.8),
  a50_down = c(9.8, 10.5, 9.5, 11.0),
  k_down = c(1.1, 0.8, 1.9, 0.7),
  nu_down = c(1.0, 1.4, 0.7, 1.3)
)

knitr::kable(richards_scenarios, digits = 2)
Table 2: Optional parameter settings for double-Richards selectivity examples.
scenario a50_up k_up nu_up a50_down k_down nu_down
Near-symmetric dome 4.5 1.2 1.0 9.8 1.1 1.0
Fast rise, slow decline 4.0 2.0 0.7 10.5 0.8 1.4
Slow rise, fast decline 5.2 0.8 1.4 9.5 1.9 0.7
Broad plateau-like dome 4.8 1.1 0.8 11.0 0.7 1.3
Show R code 
richards_up <- function(age, a50, k, nu) {
  (1 + exp(-k * (age - a50)))^(-1 / nu)
}

richards_down <- function(age, a50, k, nu) {
  1 - richards_up(age, a50, k, nu)
}

double_richards_sel <- function(age, a50_up, k_up, nu_up, a50_down, k_down, nu_down) {
  raw <- richards_up(age, a50_up, k_up, nu_up) *
    richards_down(age, a50_down, k_down, nu_down)
  raw / max(raw)
}

ages <- seq(1, 15, by = 0.1)

richards_df <- richards_scenarios |>
  dplyr::rowwise() |>
  dplyr::do({
    s <- .
    tibble::tibble(
      scenario = s$scenario,
      age = ages,
      selectivity = double_richards_sel(
        age = ages,
        a50_up = s$a50_up, k_up = s$k_up, nu_up = s$nu_up,
        a50_down = s$a50_down, k_down = s$k_down, nu_down = s$nu_down
      )
    )
  }) |>
  dplyr::ungroup()

ggplot2::ggplot(richards_df, ggplot2::aes(x = age, y = selectivity, color = scenario)) +
  ggplot2::geom_line(linewidth = 1.1) +
  ggplot2::scale_y_continuous(limits = c(0, 1.05), breaks = seq(0, 1, by = 0.25)) +
  ggplot2::labs(
    x = "Age",
    y = "Selectivity",
    color = "Scenario"
  ) +
  ggplot2::theme_bw(base_size = 11) +
  ggplot2::theme(panel.grid.minor = ggplot2::element_blank())
Figure 9: Example double-Richards selectivity curves under optional parameter settings.

In assessment practice, this flexibility can reduce structural bias relative to overly rigid parametric shapes, but it also introduces additional parameters and therefore greater identifiability risk. With limited composition information, shape parameters can become weakly determined and confound with other model components such as natural mortality, catchability, and recruitment patterns.

A defensible use case is as a candidate or sensitivity model, with regularization or parameter bounds and explicit comparison against simpler alternatives (logistic/double logistic) and semi-parametric alternatives (for example, spline selectivity). Diagnostics should include residual structure, parameter-profile behavior, and sensitivity of management quantities to the added flexibility.

Spline-Based Selectivity

3 Spline-Based Selectivity

Penalized B-spline selectivity with smoothness control.

3.1 Overview

Spline-based selectivity provides a flexible, semi-parametric alternative to logistic functional forms. Rather than imposing a specific shape (monotone or dome-shaped), a B-spline basis is used to represent selectivity as a smooth function of age, with the degree of smoothness controlled by a penalty on the second differences of the spline coefficients.

This formulation is consistent with a broader selectivity-estimation literature that emphasizes flexible shape representation together with explicit regularization [Millar and Fryer (1999); Wileman et al. (1996)]. In assessment applications, spline penalties act as the key control on effective complexity and therefore on bias-variance trade-offs [Martell and Stewart (2014)].

3.2 Model definition

Let \(\mathbf{B}(a)\) be a B-spline basis matrix evaluated at ages \(a\), with \(K\) basis functions and coefficient vector \(\boldsymbol{\beta}\). The log-selectivity is

\[ \log\,\text{sel}(a) = \mathbf{B}(a)\,\boldsymbol{\beta} \]

with selectivity normalized so the mean over a reference age range equals 1:

\[ \text{sel}(a) = \frac{\exp(\mathbf{B}(a)\,\boldsymbol{\beta})} {\frac{1}{|\mathcal{A}_{\text{ref}}|}\sum_{a' \in \mathcal{A}_{\text{ref}}} \exp(\mathbf{B}(a')\,\boldsymbol{\beta})}. \]

Here, \(\mathcal{A}_{\text{ref}}\) is a user-specified reference age range (typically ages expected to be most highly selected). This normalization keeps the objective smooth and differentiable, avoiding the non-differentiability introduced by a \(\max(\cdot)\) denominator.

Smoothness is enforced via a penalty on second differences:

\[ \text{penalty} = \frac{\lambda}{2} \sum_{k=3}^{K} (\beta_k - 2\beta_{k-1} + \beta_{k-2})^2 \]

where \(\lambda\) controls the trade-off between fit and smoothness.

3.3 Basis visualization

3.4 Example curves

3.5 RTMB implementation

Placeholder: RTMB objective with spline coefficients as parameters, smoothness penalty as a tuning parameter or estimated, and MCMC sampling.

Time-Varying Selectivity with 2D Autoregressive Structure

4 Time-Varying Selectivity with 2D Autoregressive Structure

Separable AR1 processes in year and age via RTMB

4.1 Overview

In many fisheries, selectivity is not constant over time. Changes in fishing gear, spatial effort allocation, and fish distribution can all shift the age-specific vulnerability pattern across years. A time-varying selectivity model captures these dynamics by allowing log-selectivity deviations to vary over a year \(\times\) age matrix, with temporal and ontogenetic smoothness enforced through a separable two-dimensional autoregressive prior.

Simulation and applied studies show that assumptions about time-varying selectivity can strongly influence stock-status inference, particularly when the estimation model cannot represent the data-generating dynamics [Linton and Bence (2011); Martell and Stewart (2014)]. Recent semi-parametric methods have therefore emphasized autocorrelation across both age and year dimensions, directly motivating separable latent-process formulations [Xu et al. (2019)].

This formulation uses RTMB’s dautoreg and dseparable functions to construct a Kronecker-structured precision matrix, providing a computationally efficient and differentiable prior for the year–age deviation surface.

4.2 Model definition

4.2.1 Selectivity-at-age by year

Let \(f\) index fishery, \(y\) index year, and \(a\) index age. For each fishery, a matrix of log-selectivity deviations \(\boldsymbol{\epsilon}_f\) is defined over years \(y = 1, \ldots, Y\) and estimated ages \(a = a_{\min}, \ldots, a_{\max}\):

\[ \log s_{f,y,a} = \epsilon_{f,y,a} \]

The realized selectivity is obtained by exponentiating and normalizing so that the mean across ages within each year equals 1:

\[ s_{f,y,a} = \frac{\exp(\epsilon_{f,y,a})} {\frac{1}{a_{\max} - a_{\min} + 1}\sum_{a'=a_{\min}}^{a_{\max}} \exp(\epsilon_{f,y,a'})} \] For ages beyond \(a_{\max}\), selectivity may be held constant at the value for \(a_{\max}\) (a “plus-group” extension):

\[ s_{f,y,a} = s_{f,y,a_{\max}}, \qquad a > a_{\max}. \]

4.2.2 Block structure for change years

In practice, selectivity is not re-estimated in every year. Instead, change years define blocks within which selectivity is constant. Let \(\mathbf{C}_f\) be a binary indicator matrix where \(C_{f,y} = 1\) if year \(y\) starts a new selectivity block for fishery \(f\). Within a block, selectivity repeats the values from the block’s initial year:

\[ s_{f,y,a} = \begin{cases} s_{f,y,a} & \text{if } C_{f,y} = 1 \\ s_{f,y-1,a} & \text{if } C_{f,y} = 0 \end{cases} \]

This reduces the number of free parameters from \(Y \times (a_{\max} - a_{\min} + 1)\) to \(B_f \times (a_{\max} - a_{\min} + 1)\) where \(B_f\) is the number of selectivity blocks for fishery \(f\).

4.3 Separable 2D AR1 prior

4.3.1 Marginal AR1 processes

The log-selectivity deviations for each fishery are assumed to follow a separable two-dimensional Gaussian Markov random field (GMRF). The key idea is that temporal (year) and ontogenetic (age) correlation structures are modeled independently, and their joint distribution is constructed via a Kronecker product.

Define two stationary AR1 processes:

Year dimension. For a sequence \(x_1, \ldots, x_Y\):

\[ x_y = \rho_y\, x_{y-1} + \sigma_y\,\eta_y, \qquad \eta_y \sim \mathcal{N}(0,1) \]

where \(\rho_y \in (-1, 1)\) is the year-to-year autocorrelation. The marginal variance is \(\sigma_y^2 / (1 - \rho_y^2)\) and the precision (inverse covariance) matrix for a length-\(Y\) sequence is the tridiagonal matrix \(\mathbf{Q}_y\) with:

\[ [\mathbf{Q}_y]_{ij} = \begin{cases} 1 + \rho_y^2 & i = j,\ 1 < i < Y \\ 1 & i = j = 1 \text{ or } i = j = Y \\ -\rho_y & |i - j| = 1 \\ 0 & \text{otherwise} \end{cases} \]

scaled by \(1 / \sigma_y^2\).

Age dimension. Analogously, for ages \(a_{\min}, \ldots, a_{\max}\):

\[ z_a = \rho_a\, z_{a-1} + \sigma_a\,\nu_a, \qquad \nu_a \sim \mathcal{N}(0,1) \]

with autocorrelation \(\rho_a\) and precision matrix \(\mathbf{Q}_a\) of the same tridiagonal form.

4.3.2 Kronecker product structure

The joint precision matrix for the vectorized year–age deviation \(\text{vec}(\boldsymbol{\epsilon}_f)\) is the Kronecker product:

\[ \mathbf{Q} = \mathbf{Q}_y \otimes \mathbf{Q}_a \]

This encodes the assumption that the year and age correlation structures are separable: the correlation between \(\epsilon_{y,a}\) and \(\epsilon_{y',a'}\) factorizes as

\[ \text{Cor}(\epsilon_{y,a},\, \epsilon_{y',a'}) = \rho_y^{|y - y'|}\,\rho_a^{|a - a'|} \]

The corresponding joint covariance matrix is:

\[ \boldsymbol{\Sigma} = \boldsymbol{\Sigma}_y \otimes \boldsymbol{\Sigma}_a \]

where \(\boldsymbol{\Sigma}_y\) and \(\boldsymbol{\Sigma}_a\) are the marginal AR1 covariance matrices for year and age, respectively.

4.3.3 Marginal scale

The overall marginal standard deviation of any single element \(\epsilon_{f,y,a}\) is

\[ \text{SD}(\epsilon_{f,y,a}) = \frac{\sigma_f} {\sqrt{1 - \rho_{y,f}^2}\;\sqrt{1 - \rho_{a,f}^2}} \]

where \(\sigma_f\) is the innovation scale for fishery \(f\). This quantity (scale in the RTMB code) ensures that the specified \(\sigma_f\) represents the innovation standard deviation while the marginal variance accounts for both autocorrelation parameters.

4.3.4 Log-density

The log-density of the separable GMRF, evaluated by RTMB’s dseparable, is:

\[ \log p(\boldsymbol{\epsilon}_f \mid \rho_y, \rho_a, \sigma) = -\frac{YA}{2}\log(2\pi) + \frac{A}{2}\log|\mathbf{Q}_y| + \frac{Y}{2}\log|\mathbf{Q}_a| - \frac{1}{2}\,\text{vec}(\boldsymbol{\epsilon}_f)^\top \left(\mathbf{Q}_y \otimes \mathbf{Q}_a\right) \text{vec}(\boldsymbol{\epsilon}_f) \]

where \(A = a_{\max} - a_{\min} + 1\) is the number of estimated ages. The Kronecker structure means the log-determinant decomposes as a sum of the individual log-determinants (scaled by dimension), and the quadratic form can be evaluated without materializing the full \(YA \times YA\) precision matrix.

4.4 RTMB implementation

The implementation uses dseparable to compose two dautoreg densities. The dautoreg(x, phi) function evaluates the log-density of an AR1 process with autocorrelation phi, and dseparable(f1, f2) constructs their Kronecker product density.

4.4.1 dseparable mechanics

The call dseparable(f1, f2)(X, scale) performs the following:

  1. Treats the matrix X (dimension \(B_f \times A\)) as a 2D random field.
  2. Evaluates f1 along each row (year marginal) and f2 along each column (age marginal).
  3. Combines them via the Kronecker identity to compute the joint log-density without forming the full \(B_f A \times B_f A\) precision matrix.
  4. The scale parameter rescales the entire field so that the marginal standard deviation matches the specified value.

This is efficient because the computational cost scales as \(\mathcal{O}(B_f \cdot A)\) rather than \(\mathcal{O}((B_f \cdot A)^3)\).

4.4.2 Selectivity construction

The get_selectivity function builds the full 3D array (fishery \(\times\) year \(\times\) age) from the estimated log-selectivity blocks:

4.5 Visualization

4.6 Properties and considerations

The separable 2D AR1 formulation has several properties that make it attractive for operational stock assessments:

Parsimony. Three hyperparameters per fishery (\(\rho_y\), \(\rho_a\), \(\sigma\)) control the smoothness of an arbitrarily large year–age surface. The autocorrelation parameters shrink the effective number of free parameters well below the nominal \(B_f \times A\).

Computational efficiency. The Kronecker structure avoids forming or factorizing the full joint precision matrix. Combined with RTMB’s automatic differentiation, gradients of the log-prior are computed at the same \(\mathcal{O}(B_f \cdot A)\) cost as the density itself.

Interpretability. The parameter \(\rho_y\) directly controls how quickly selectivity can change between years (high values produce slow drift, low values allow abrupt shifts), while \(\rho_a\) controls the smoothness of the selectivity curve within any single year.

Integration with block structure. The change-year mechanism can be layered on top: the 2D AR1 prior applies to the matrix of block-specific selectivity patterns, not to every individual year. This reduces dimensionality further while still allowing the prior to borrow strength across blocks.

4.6.1 Comparison of time-varying selectivity approaches

Show R code 
compare_tv <- data.frame(
  Approach = c(
    "Structured temporal selectivity (Butterworth et al. 2003)",
    "Separable AR1 surface (this manuscript)",
    "Block/time‑period selectivity"
  ),
  Structure = c(
    "Pre‑specified temporal structure on selectivity changes",
    "Autocorrelated deviations across age and time",
    "Piecewise‑constant selectivity blocks"
  ),
  Strengths = c(
    "Historical applied example; likelihood‑based CAA with error",
    "Smooth evolution; explicit regularization",
    "Transparent and easy to communicate"
  ),
  Risks = c(
    "Potential confounding if structure is too rigid",
    "Over‑flexibility if penalties are weak",
    "Arbitrary block boundaries"
  ),
  check.names = FALSE
)

knitr::kable(compare_tv)
Table 3: Comparison of time‑varying selectivity approaches in applied assessments.
Approach Structure Strengths Risks
Structured temporal selectivity (Butterworth et al. 2003) Pre‑specified temporal structure on selectivity changes Historical applied example; likelihood‑based CAA with error Potential confounding if structure is too rigid
Separable AR1 surface (this manuscript) Autocorrelated deviations across age and time Smooth evolution; explicit regularization Over‑flexibility if penalties are weak
Block/time‑period selectivity Piecewise‑constant selectivity blocks Transparent and easy to communicate Arbitrary block boundaries

4.7 Comparison and Discussion

The five selectivity families evaluated here span a practical gradient from low-dimensional structural assumptions to high-dimensional latent-process flexibility. The logistic model (Section; Figure) provides an interpretable baseline with strong monotonicity assumptions and minimal parameter burden. The double-logistic model (Section; Table, Figure, Figure) adds biologically plausible dome-shape behavior, but also increases confounding risk with natural mortality and cohort effects when data support is limited (Thompson 1994).

The generalized logistic (Richards) option (Section) extends this parametric family by allowing asymmetric transitions and plateau-like behavior. Its additional flexibility can improve fit when conventional logistic forms are too restrictive, but the associated parameter expansion requires stronger identifiability checks and careful regularization.

Spline selectivity (Section; Figure; Figure) offers useful middle ground: flexible shape estimation with transparent control of roughness via penalization. In review terms, this is a semi-parametric compromise between rigid parametric forms and fully dynamic latent surfaces. The time-varying 2D AR1 approach (Section; Figure; Figure) further extends flexibility by allowing coherent year-by-age evolution with explicit process structure, aligning with recent recommendations for autocorrelated selectivity deviations (Linton and Bence 2011; Xu et al. 2019).

For scientific review and operational assessment use, a defensible workflow is:

  1. Start with a parsimonious baseline (logistic or double logistic) and diagnose fit deficiencies in composition residuals and retrospective patterns.
  2. Introduce additional flexibility (spline or time-varying structures) only when diagnostics indicate systematic lack of fit and data support richer structure.
  3. Compare candidate models with complementary evidence, including information criteria, predictive checks, and sensitivity of management quantities (Punt, Hurtado-Ferro, and Whitten 2013; Wilberg and Bence 2008).
  4. Report regularization choices and identifiability diagnostics explicitly, because inferred selectivity dynamics are conditional on prior/penalty assumptions (Martell and Stewart 2014; Privitera-Johnson, Methot, and Punt 2022).

Overall, the central trade-off is not whether one curve family is universally “best”, but whether the assumed latent structure is commensurate with the available information and management objectives.

References

Butterworth, D. S., J. N. Ianelli, and R. Hilborn. 2003. “A Statistical Model for Stock Assessment of Southern Bluefin Tuna with Temporal Changes in Selectivity.” African Journal of Marine Science 25 (1): 331–61. https://doi.org/10.2989/18142320309504021.
Linton, Brenden C., and James R. Bence. 2011. “Catch-at-Age Assessment in the Face of Time-Varying Selectivity.” ICES Journal of Marine Science 68 (3): 611–25. https://doi.org/10.1093/icesjms/fsq173.
Martell, Steven J. D., and Ian J. Stewart. 2014. “Towards Defining Good Practices for Modeling Time-Varying Selectivity.” Fisheries Research 158: 84–95. https://doi.org/10.1016/j.fishres.2013.11.001.
Methot, Richard D., and Chantell R. Wetzel. 2013. “Stock Synthesis: A Biological and Statistical Framework for Fish Stock Assessment and Fishery Management.” Fisheries Research 142: 86–99. https://doi.org/10.1016/j.fishres.2012.10.012.
Millar, R. B., and R. J. Fryer. 1999. “Estimating the Size-Selection Curves of Towed Gears, Traps, Nets and Hooks.” Reviews in Fish Biology and Fisheries 9: 89–116. https://doi.org/10.1023/A:1008838220001.
Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. 1992. Numerical Recipes in c: The Art of Scientific Computing. 2nd ed. Cambridge University Press. https://cnl.gmu.edu/facility/NumericalRecipesinC.pdf.
Privitera-Johnson, Kristin, Richard D. Methot, and André E. Punt. 2022. “Towards Best Practice for Specifying Selectivity in Age-Structured Integrated Stock Assessments.” Fisheries Research 248: 106248. https://doi.org/10.1016/j.fishres.2022.106248.
Punt, André E., Felipe Hurtado-Ferro, and Athol R. Whitten. 2013. “Model Selection for Selectivity in Fisheries Stock Assessments.” Fisheries Research 141: 1–12. https://doi.org/10.1016/j.fishres.2013.02.008.
Thompson, Grant G. 1994. “Confounding of Gear Selectivity and the Natural Mortality Rate in Cases Where the Former Is a Nonmonotone Function of Age.” Canadian Journal of Fisheries and Aquatic Sciences 51 (12): 2654–65. https://doi.org/10.1139/f94-265.
Wilberg, Michael J., and James R. Bence. 2006. “Performance of Time-Varying Catchability Estimators in Statistical Catch-at-Age Analysis.” Canadian Journal of Fisheries and Aquatic Sciences 63: 2275–85. https://doi.org/10.1139/f06-111.
———. 2008. “Performance of Deviance Information Criterion Model Selection in Statistical Catch-at-Age Analysis.” Fisheries Research 93: 212–21. https://doi.org/10.1016/j.fishres.2008.04.010.
Wileman, D. A., R. S. T. Ferro, R. Fonteyne, and R. B. Millar. 1996. “Manual of Methods of Measuring the Selectivity of Towed Fishing Gears.” ICES Cooperative Research Report 215. Copenhagen: International Council for the Exploration of the Sea (ICES). https://ices-library.figshare.com/articles/report/Manual_of_methods_of_measuring_the_selectivity_of_towed_fishing_gears/18624446/1/files/33403508.pdf.
Xu, Haikun, James T. Thorson, Richard D. Methot, and Ian G. Taylor. 2019. “A New Semi-Parametric Method for Autocorrelated Age- and Time-Varying Selectivity in Age-Structured Assessment Models.” Canadian Journal of Fisheries and Aquatic Sciences 76: 268–85. https://doi.org/10.1139/cjfas-2017-0446.