1.1 What is aLBI?

aLBI (Assessment of Length-Based Indicators) is an R package designed for data-limited fish stock assessment using only catch length-frequency data. The package operationalises the Froese (2004) length-based sustainability indicators and the Cope & Punt (2009) decision framework for estimating the probability of a stock falling below target and limit spawning biomass reference points.

Most fish stock assessments in tropical and coastal developing-country fisheries suffer from the absence of the age, tagging, or time-series abundance data required by classical assessment models. aLBI was developed to address this gap: given a single snapshot of the catch length distribution, it estimates key biological reference lengths, computes three sustainability indicators, and quantifies their uncertainty through a novel three-tier simulation framework.

1.2 What’s New in the Current Version

This version introduces several major methodological improvements to FishPar and FishSS, and adds two new functions (FreqTM, LWR):

1.3 Package Citation

If you use aLBI in published research, please cite:

Ali, A., Sarker, M. R., & Alam, M. S. (2025). Development of a simple R package (aLBI) for the estimation of stock status from the length frequency data. Fisheries Research, 288, 107467. https://doi.org/10.1016/j.fishres.2025.107467

2.0.1 Required Dependencies

aLBI depends on base R packages (stats, graphics, grDevices, utils). The openxlsx package is required only when save_output = TRUE in FishPar.

required_packages <- c("aLBI", "readxl", "openxlsx", "dplyr", "ggplot2")

for (pkg in required_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    warning(paste("Package", pkg, "is required but not installed.",
                  "Install with: install.packages('", pkg, "')"))
  } else {
    suppressPackageStartupMessages(library(pkg, character.only = TRUE))
  }
}

4.1 Input Data Formats

aLBI functions accept data in two general formats:

Raw individual measurements (for FrequencyTable, FreqTM, LWR):

Aggregated length-frequency data (for FishPar):

Data may be loaded from .xlsx, .csv, or created directly in R. The package ships with example datasets in inst/exdata/.

# List all bundled example datasets
list.files(system.file("exdata", package = "aLBI"))
# ExData.xlsx    — raw individual fish lengths
# LC.xlsx        — length-frequency table (for FishPar)
# lenfreqM.xlsx  — multi-month individual lengths (for FreqTM)
# cpdata.xlsx    — Cope & Punt (2009) look-up table
# LWdata.xlsx    — paired length-weight measurements

5.1 Methodological Background

Before any length-based analysis can be performed, individual fish length measurements must be grouped into discrete length classes (bins). The choice of bin width directly affects the shape of the apparent length-frequency distribution, the number of identifiable cohort modes, and the precision of derived indicators.

FrequencyTable uses the Optimum Bin Size (OBS) formula of Wang et al. (2020) to automatically select a biologically appropriate bin width when the user does not specify one. The OBS formula is:

\[h = 0.9 \cdot \min(s,\ IQR/1.34) \cdot n^{-1/5}\]

where s is the standard deviation of lengths, IQR is the interquartile range, and n is the sample size. The resulting h is rounded to the nearest integer for practical use.

The upper bound of each class is used as the class representative value (e.g., the class [9, 10) is represented by 10), consistent with common practice in fisheries length-frequency analysis.

5.2 Function Arguments

5.3 Usage Example

library(readxl)

# Load the bundled example raw-length dataset
lenfreq_path <- system.file("exdata", "ExData.xlsx", package = "aLBI")

if (lenfreq_path == "") {
  stop("ExData.xlsx not found. Ensure aLBI is correctly installed.")
}

length_data <- readxl::read_excel(lenfreq_path)
cat("Dataset dimensions:", nrow(length_data), "rows x", ncol(length_data), "columns\n")
#> Dataset dimensions: 1177 rows x 1 columns
head(length_data, 8)
#> # A tibble: 8 × 1
#>   Length
#>    <dbl>
#> 1   75  
#> 2   66.7
#> 3   66  
#> 4   64  
#> 5   63.3
#> 6   63.2
#> 7   62.3
#> 8   62.2

# Run FrequencyTable with automatic bin width (Wang's OBS formula)
freq_result <- FrequencyTable(
  data        = length_data,
  bin_width   = NULL,     # auto-calculated via Wang's formula
  Lmax        = NULL,     # use observed maximum
  output_file = "FrequencyTable_Output.xlsx"
)

# Inspect outputs
freq_result$lfqTable   # Full frequency distribution table with class intervals
freq_result$lfreq      # Condensed table: upper class boundary and frequency

# Inline reproducible example with generated data
set.seed(42)
example_lengths <- data.frame(
  Length = round(rnorm(300, mean = 28, sd = 6), 1)
)
example_lengths$Length <- pmax(10, pmin(50, example_lengths$Length))

cat("Sample: n =", nrow(example_lengths),
    "| mean =", round(mean(example_lengths$Length), 2),
    "| range:", round(min(example_lengths$Length), 1),
    "–", round(max(example_lengths$Length), 1), "cm\n")
#> Sample: n = 300 | mean = 27.87 | range: 10 – 44.2 cm

5.4 Interpreting the Output

FrequencyTable returns a list with two elements:

• $lfqTable: The complete frequency distribution table with class intervals (lower bound, upper bound), midpoint, frequency, relative frequency (%), and cumulative frequency.
• $lfreq: A condensed two-column data frame (Length, Frequency) suitable for direct use in FishPar.

Tip: The $lfreq component is the recommended input format for FishPar. Save it as an Excel file using the output_file argument for reproducibility.

6.1 Methodological Background

Seasonal variation in catch length composition is ecologically informative — it can reflect recruitment pulses, migratory behaviour, or gear-based seasonal selectivity. FreqTM constructs length-frequency tables for each month in a dataset using a consistent bin structure derived once from the full dataset, ensuring that length classes are directly comparable across months.

This is critical for analyses such as ELEFAN (Electronic Length Frequency Analysis) where the same class structure must be maintained across all time periods.

6.2 Function Arguments

6.3 Usage Example

# Load multi-month length data
lenfreqM_path <- system.file("exdata", "lenfreqM.xlsx", package = "aLBI")

if (lenfreqM_path == "") {
  stop("lenfreqM.xlsx not found. Ensure aLBI is correctly installed.")
}

monthly_data <- readxl::read_excel(lenfreqM_path)
cat("Months present:", paste(unique(monthly_data[[1]]), collapse = ", "), "\n")
#> Months present: Jan, Feb, Mar, Apr, May, Jun, Jul, Aug
cat("Total observations:", nrow(monthly_data), "\n")
#> Total observations: 1500
head(monthly_data, 6)
#> # A tibble: 6 × 2
#>   Months Length
#>   <chr>   <dbl>
#> 1 Jan       9.3
#> 2 Jan      11.5
#> 3 Jan       8.7
#> 4 Jan       7.4
#> 5 Jan       8.6
#> 6 Jan       9

# Generate monthly frequency tables with consistent bin structure
monthly_freq <- FreqTM(
  data        = monthly_data,
  bin_width   = NULL,
  Lmax        = NULL,
  date_config = list(day = 15, year = 2024),
  output_file = "FreqTM_Output.xlsx"
)

# View results for the first month
monthly_freq[[1]]

6.4 Notes on Date Configuration

The date_config argument converts month names (e.g., "January", "February") into date objects formatted as day.month.year (e.g., 15.01.2024). This facilitates direct use in time-series growth analyses that require date-indexed length data.

7.1 Methodological Framework

FishPar is the core function of the aLBI package. It implements a rigorous two-stage simulation pipeline:

1. Stage 1 — Monte Carlo (MC) simulation for biological reference lengths (L_max, L_inf, L_mat, L_opt).
2. Stage 2 — Three-tier uncertainty propagation for Froese (2004) sustainability indicators (P_mat, P_opt, P_mega).

The separation of these two stages is methodologically important and reflects a key advance over both the classical manual calculation approach and earlier versions of the package.

7.2 Monte Carlo Simulation for Length Parameters

7.3 Froese (2004) Sustainability Indicators

Three length-based sustainability indicators are computed from the catch composition and the estimated reference lengths.

P_mat — proportion of mature fish in the catch:

\[P_{\text{mat}} = \frac{\sum_{k: L_k \geq L_{\text{mat}}} f_k}{N} \times 100\%\]

Target: P_mat = 100% (all caught fish should be mature).

P_opt — proportion of fish at optimal harvest size:

\[P_{\text{opt}} = \frac{\sum_{k: L_{\text{opt\_m10}} \leq L_k \leq L_{\text{opt\_p10}}} f_k}{N} \times 100\%\]

Target: P_opt = 100% (ideally all fish are in the optimal size window).

P_mega — proportion of mega-spawners (large, highly fecund fish):

\[P_{\text{mega}} = \frac{\sum_{k: L_k > L_{\text{opt\_p10}}} f_k}{N} \times 100\%\]

Target: P_mega ≥ 20% (at least 20% of the catch should be large, reproductively valuable fish).

where \(L_k\) is the midpoint of length class \(k\), \(f_k\) is its catch frequency, and \(N = \sum_k f_k\) is the total sample size.

7.4 Three-Tier Uncertainty Propagation (Novel)

The critical methodological advance of the current FishPar is its three-tier uncertainty decomposition. This directly addresses the fundamental limitation of manual and classical bootstrap approaches: both compute P_mat, P_opt, and P_mega using fixed point-estimate values of L_mat and L_opt, thereby ignoring the substantial uncertainty in those parameters.

7.5 Function Arguments

7.6 Usage Example

library(readxl)

# Load the bundled length-frequency dataset
lf_path <- system.file("exdata", "LC.xlsx", package = "aLBI")
if (lf_path == "") stop("LC.xlsx not found. Please reinstall aLBI.")

lf_data <- readxl::read_excel(lf_path)
cat("Length-frequency data:\n")
#> Length-frequency data:
print(lf_data)
#> # A tibble: 15 × 2
#>    LengthClass Frequency
#>          <dbl>     <dbl>
#>  1          12        26
#>  2          15       166
#>  3          18       244
#>  4          21       582
#>  5          24       973
#>  6          27      1067
#>  7          30       963
#>  8          33       511
#>  9          36       472
#> 10          39       286
#> 11          42       173
#> 12          45       171
#> 13          48        83
#> 14          51        36
#> 15          54        36
cat("\nTotal individuals:", sum(lf_data[[2]]), "\n")
#> 
#> Total individuals: 5789

# ── Basic usage (plots shown on screen; no files written) ──────────────────
results <- FishPar(
  data        = lf_data,
  resample    = 1000,
  progress    = FALSE,
  Linf        = NULL,    # derive L_inf from L_max
  Linf_sd     = 0.5,
  Lmat        = NULL,    # derive L_mat from allometric regression
  Lmat_sd     = 0.5,
  save_output = FALSE    # set TRUE to write PDFs + Excel to disk
)

# ── With known L_inf from FishBase or literature ──────────────────────────
results_linf <- FishPar(
  data        = lf_data,
  resample    = 1000,
  Linf        = 18.5,    # species-specific L_inf (cm)
  Linf_sd     = 0.5,
  save_output = FALSE
)

# ── Save all outputs to working directory ─────────────────────────────────
results_saved <- FishPar(
  data        = lf_data,
  resample    = 5000,    # increase for publication-quality CIs
  save_output = TRUE     # writes PDFs + FishPar_Results.xlsx
)

7.7 Accessing Results

# ── Length parameters (MC mean ± 95% CI) ─────────────────────────────────
results$estimated_length_par

# ── Three-tier Froese indicator CIs ──────────────────────────────────────
results$froese_par_mc         # Tier 1: parameter uncertainty only
results$froese_par_bootstrap  # Tier 2: data-sampling uncertainty only
results$froese_par_joint      # Tier 3: total uncertainty (recommended)

# ── Derived scalars ───────────────────────────────────────────────────────
results$LM_ratio    # L_mat / L_opt
results$Pobj        # P_mat + P_opt + P_mega (0–300 scale, Cope & Punt 2009)
results$Total_ind   # total observed individuals

# ── Target comparison ─────────────────────────────────────────────────────
results$froese_ind_vs_target

$estimated_length_par
  Parameters Mean_estimate Lower_CI Upper_CI
1       Lmax         18.00    16.20    19.80
2       Linf         18.95    17.05    20.84
3       Lmat         11.71    10.45    13.09
4       Lopt         11.72    10.46    13.10
5   Lopt_p10         12.89    11.51    14.41
6   Lopt_m10         10.55     9.41    11.79

$froese_par_joint                 # Tier 3 — primary reporting CI
  Parameters  Mean Lower_CI Upper_CI          Source
1       Pmat 56.47    28.78    77.13  Tier3_Joint_total
2       Popt 38.38    26.78    52.96  Tier3_Joint_total
3      Pmega 36.40    14.22    62.68  Tier3_Joint_total

7.8 Generated Plots

When save_output = TRUE, FishPar writes ten PDF files to getwd():

7.9 The Excel Workbook (FishPar_Results.xlsx)

The Excel workbook contains six sheets:

The FishSS_Inputs sheet is designed to streamline the workflow: all values required by FishSS() can be read directly from this sheet without manual extraction.

8.1 Methodological Background

FishSS implements the Cope & Punt (2009) decision framework for length-based stock status assessment. It uses the composite objective index P_obj and the length-maturity ratio LM_ratio to select the appropriate columns from the published look-up table (Table 5, Cope & Punt 2009), then interpolates to find:

• p_below_target: Probability that current spawning biomass is below 0.4 × SB₀ (target reference point)
• p_below_limit: Probability that current spawning biomass is below 0.25 × SB₀ (limit reference point)

8.2 Function Arguments

8.3 Usage Example

# Load the Cope & Punt (2009) look-up table
cpdata_path <- system.file("exdata", "cpdata.xlsx", package = "aLBI")
if (cpdata_path == "") stop("cpdata.xlsx not found. Reinstall aLBI.")

cpdata <- readxl::read_excel(cpdata_path)
cat("CPdata structure:", nrow(cpdata), "rows ×", ncol(cpdata), "columns\n")
#> CPdata structure: 21 rows × 11 columns
cat("Columns:", paste(colnames(cpdata), collapse = ", "), "\n")
#> Columns: Tx, A, B, C, D, E, F, G, H, I, J
print(head(cpdata, 6))
#> # A tibble: 6 × 11
#>      Tx     A     B     C     D     E     F     G     H     I     J
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1   100     0     0     0     0     0     0     0     0   100   100
#> 2    95     0     0     0     0    22     0    11     0   100    93
#> 3    90     0     0     0     0   100    44    83    22   100    74
#> 4    85     0     0     0     0   100   100   100    67   100    63
#> 5    80     0     0     0     0   100   100   100   100    89    52
#> 6    75     0     0     0     0   100   100   100   100    74    37

# ── Method 1: Using FishPar results directly ──────────────────────────────
# (Recommended — uses Tier 3 joint means as point estimates)
stock_status <- FishSS(
  data     = cpdata,
  LM_ratio = results$LM_ratio,
  Pmat     = results$froese_par_joint$Mean[1],  # Tier 3 mean P_mat
  Popt     = results$froese_par_joint$Mean[2],  # Tier 3 mean P_opt
  Pmega    = results$froese_par_joint$Mean[3]   # Tier 3 mean P_mega
)

# ── Method 2: Using FishSS_Inputs sheet from saved Excel ─────────────────
library(openxlsx)
ss_inputs <- openxlsx::read.xlsx("FishPar_Results.xlsx",
                                  sheet = "FishSS_Inputs")

stock_status <- FishSS(
  data     = cpdata,
  LM_ratio = ss_inputs$Value[ss_inputs$Parameter == "LM_ratio"],
  Pmat     = ss_inputs$Value[ss_inputs$Parameter == "Pmat"],
  Popt     = ss_inputs$Value[ss_inputs$Parameter == "Popt"],
  Pmega    = ss_inputs$Value[ss_inputs$Parameter == "Pmega"]
)

# View the full stock status assessment
stock_status

stock_status

# $Pobj
# [1] 131.25     <- P_mat + P_opt + P_mega (0–300 scale)

# $Px_trigger
# [1] "Pmat"     <- trigger indicator (P_mat when Pobj < 200)

# $Px_value
# [1] 56.47      <- value of trigger indicator (%)

# $StockStatus
# p_below_target  p_below_limit
#            100            100   <- probabilities (%) from look-up table

# $Selectivity
# [1] "Fish maturity ogive"   <- catch selectivity pattern

8.4 Interpreting Stock Status

The two probability values returned in $StockStatus should be interpreted together:

Important: These probabilities are based on a simulation study integrating steepness values and assume that the catch length distribution is representative of the exploited population. They should be interpreted alongside biological knowledge of the target species.

9.1 Methodological Background

The length–weight relationship (LWR) describes the power-law relationship between fish body length (L) and weight (W):

\[W = a \cdot L^b\]

Taking logarithms linearises this to:

\[\ln(W) = \ln(a) + b \cdot \ln(L)\]

The parameter b is biologically interpretable:

• b ≈ 3.0: isometric growth — body shape is constant as fish grow.
• b < 3.0: negative allometry — fish become more slender with length.
• b > 3.0: positive allometry — fish become more robust (heavier) with length.

9.2 Function Arguments

9.3 Usage Example

# Load the bundled length-weight dataset
lw_path <- system.file("exdata", "LWdata.xlsx", package = "aLBI")
if (lw_path == "") stop("LWdata.xlsx not found. Reinstall aLBI.")

LWdata <- readxl::read_excel(lw_path)
cat("LW dataset: n =", nrow(LWdata), "fish\n")
#> LW dataset: n = 554 fish
cat("Length range:", round(min(LWdata[[1]]), 1), "–",
    round(max(LWdata[[1]]), 1), "cm\n")
#> Length range: 11.4 – 58.7 cm
cat("Weight range:", round(min(LWdata[[2]]), 2), "–",
    round(max(LWdata[[2]]), 2), "g\n")
#> Weight range: 11 – 1120 g
head(LWdata)
#> # A tibble: 6 × 2
#>   Length Weight
#>    <dbl>  <dbl>
#> 1   58.7   1025
#> 2   58.6   1013
#> 3   54.6   1011
#> 4   54.3    941
#> 5   54.2    946
#> 6   53     1120

# Fit and plot the length-weight relationship
lwr_result <- LWR(
  data          = LWdata,
  log_transform = TRUE,      # log-log linearisation (recommended)
  point_col     = "black",
  line_col      = "#C0392B", # red regression line
  shade_col     = "#C0392B",
  point_size    = 2,
  line_size     = 1,
  alpha         = 0.2,
  main          = "Length-Weight Relationship",
  xlab          = NULL,      # auto-label: "log(Length)" or "Length"
  ylab          = NULL,      # auto-label: "log(Weight)" or "Weight"
  save_output   = FALSE      # set TRUE to write PDF + TXT to disk
)

# Model summary
lwr_result$model_summary   # intercept, slope, r², p-value
lwr_result$plot            # ggplot2 object

11.1 The Core Problem: Fixed vs. Distributed Parameters

A common question from users (and a valid peer-review concern) is why the P_mat, P_opt, and P_mega values from FishPar differ from values computed by hand using the same Froese & Binohlan (2000) regressions. This section explains the discrepancy with a worked numerical example.

11.2 Worked Example

Consider the following dataset (N = 844 fish, bin width = 1 cm):

# Example length-frequency data
example_lf <- data.frame(
  Length    = 6:18,
  Frequency = c(2, 3, 22, 54, 110, 131, 158, 141, 90, 63, 45, 23, 2)
)
N_total <- sum(example_lf$Frequency)
cat("N =", N_total, "fish; length range:", min(example_lf$Length),
    "–", max(example_lf$Length), "cm\n")
#> N = 844 fish; length range: 6 – 18 cm
print(example_lf)
#>    Length Frequency
#> 1       6         2
#> 2       7         3
#> 3       8        22
#> 4       9        54
#> 5      10       110
#> 6      11       131
#> 7      12       158
#> 8      13       141
#> 9      14        90
#> 10     15        63
#> 11     16        45
#> 12     17        23
#> 13     18         2

# Step 1: derive fixed point-estimate parameters
Lmax_obs <- max(example_lf$Length)       # 18 cm
Linf_pt  <- Lmax_obs / 0.95             # 18.947 cm
Lmat_pt  <- 10^(0.8979 * log10(Linf_pt) - 0.0782)   # 11.71 cm
Lopt_pt  <- 10^(1.053  * log10(Lmat_pt) - 0.0565)   # 11.72 cm
Lopt_p10_pt <- Lopt_pt * 1.10                         # 12.89 cm
Lopt_m10_pt <- Lopt_pt * 0.90                         # 10.55 cm

cat(sprintf(
  "Fixed parameters:\n  Lmax = %.2f, Linf = %.3f\n  Lmat = %.2f, Lopt = %.2f\n  Lopt_m10 = %.2f, Lopt_p10 = %.2f\n",
  Lmax_obs, Linf_pt, Lmat_pt, Lopt_pt, Lopt_m10_pt, Lopt_p10_pt
))
#> Fixed parameters:
#>   Lmax = 18.00, Linf = 18.947
#>   Lmat = 11.72, Lopt = 11.72
#>   Lopt_m10 = 10.55, Lopt_p10 = 12.90

# Step 2: classify each length class
example_lf$Mature  <- example_lf$Length >= Lmat_pt
example_lf$Optimal <- example_lf$Length >= Lopt_m10_pt &
                      example_lf$Length <= Lopt_p10_pt
example_lf$Mega    <- example_lf$Length >  Lopt_p10_pt

# Step 3: compute manual Froese indicators
Pmat_manual  <- 100 * sum(example_lf$Frequency[example_lf$Mature])  / N_total
Popt_manual  <- 100 * sum(example_lf$Frequency[example_lf$Optimal]) / N_total
Pmega_manual <- 100 * sum(example_lf$Frequency[example_lf$Mega])    / N_total

cat(sprintf(
  "\nManual Froese Indicators:\n  Pmat  = %.2f%%\n  Popt  = %.2f%%\n  Pmega = %.2f%%\n  (No CIs possible)\n",
  Pmat_manual, Popt_manual, Pmega_manual
))
#> 
#> Manual Froese Indicators:
#>   Pmat  = 61.85%
#>   Popt  = 34.24%
#>   Pmega = 43.13%
#>   (No CIs possible)

# Demonstrate the step-function effect of Lmat on Pmat
Lmat_values <- seq(9, 17, by = 0.1)

compute_Pmat <- function(Lmat, lf) {
  100 * sum(lf$Frequency[lf$Length >= Lmat]) / sum(lf$Frequency)
}

Pmat_curve <- sapply(Lmat_values, compute_Pmat, lf = example_lf)

cat("How Pmat changes as Lmat varies (step-function):\n")
#> How Pmat changes as Lmat varies (step-function):
df_steps <- data.frame(
  Lmat_range  = c("≤ 10", "(10,11]", "(11,12]", "(12,13]", "(13,14]", "> 14"),
  Pmat        = c(
    compute_Pmat(10.0, example_lf),
    compute_Pmat(10.5, example_lf),
    compute_Pmat(11.5, example_lf),  # ← Manual falls here
    compute_Pmat(12.5, example_lf),
    compute_Pmat(13.5, example_lf),
    compute_Pmat(14.5, example_lf)
  ),
  Notes = c("", "", "← Manual (Lmat=11.71)", "", "", "")
)
print(df_steps)
#>   Lmat_range     Pmat                 Notes
#> 1       ≤ 10 90.40284                      
#> 2    (10,11] 77.36967                      
#> 3    (11,12] 61.84834 ← Manual (Lmat=11.71)
#> 4    (12,13] 43.12796                      
#> 5    (13,14] 26.42180                      
#> 6       > 14 15.75829

# The jump at the 12-cm class boundary
jump_12cm <- 100 * example_lf$Frequency[example_lf$Length == 12] / N_total
cat(sprintf(
  "\nJump at Lmat = 12.0 cm: %.1f pp (= 158 fish / 844 total)\n",
  jump_12cm
))
#> 
#> Jump at Lmat = 12.0 cm: 18.7 pp (= 158 fish / 844 total)
cat("A 0.2 cm shift in Lmat (11.9 → 12.1) changes Pmat by", jump_12cm, "pp\n")
#> A 0.2 cm shift in Lmat (11.9 → 12.1) changes Pmat by 18.72038 pp

# Model Tier 3 outputs (from FishPar with this dataset)
comparison <- data.frame(
  Indicator = c("Pmat", "Popt", "Pmega"),
  Manual    = c(61.85, 34.24, 43.13),
  Model_T2  = c(61.85, 34.24, 43.13),  # Tier 2 ≈ manual
  Model_T3_Mean = c(56.47, 38.38, 36.40),
  T3_Lower95  = c(28.78, 26.78, 14.22),
  T3_Upper95  = c(77.13, 52.96, 62.68)
)
print(comparison)
#>   Indicator Manual Model_T2 Model_T3_Mean T3_Lower95 T3_Upper95
#> 1      Pmat  61.85    61.85         56.47      28.78      77.13
#> 2      Popt  34.24    34.24         38.38      26.78      52.96
#> 3     Pmega  43.13    43.13         36.40      14.22      62.68

cat("\nKey insight: Tier 2 ≈ Manual because both use the same fixed Lmat/Lopt.\n")
#> 
#> Key insight: Tier 2 ≈ Manual because both use the same fixed Lmat/Lopt.
cat("Tier 3 differs because it integrates over the MC distribution of Lmat,\n")
#> Tier 3 differs because it integrates over the MC distribution of Lmat,
cat("propagating the 2.64 cm uncertainty range [10.45, 13.09] through the\n")
#> propagating the 2.64 cm uncertainty range [10.45, 13.09] through the
cat("step function, which shifts the expected Pmat from 61.85% to 56.47%.\n")
#> step function, which shifts the expected Pmat from 61.85% to 56.47%.

11.3 Why Manual Calculation Cannot Produce Confidence Intervals

The manual approach has four fundamental limitations:

1. Fixed parameters: L_mat = 11.71 cm is treated as the exact true value, ignoring that it is an estimate from a cross-species allometric regression.
2. No resampling framework: Without MC simulation or bootstrapping, there is no distributional information from which a CI can be constructed.
3. Step-function incompatibility: Even applying the delta method is impossible here — the step function is not differentiable at class boundaries, violating the smooth- function assumption that the delta method requires.
4. Two interacting uncertainty sources: Parameter estimation uncertainty and catch-data sampling uncertainty cannot be disentangled analytically; only numerical simulation (the three-tier framework) can quantify both simultaneously.

Summary: The 5.4 pp difference between manual (61.85%) and model (56.47%) P_mat is not a computational error. It is the correct, quantified effect of propagating allometric regression uncertainty through the step-function classification. The three-tier framework in FishPar is the only approach capable of providing honest, defensible confidence intervals for length-based sustainability indicators in data-limited fisheries.

12.1 Setting a Random Seed

All MC and bootstrap simulations in FishPar use R’s default pseudorandom number generator. For exactly reproducible results, set a seed before calling FishPar:

set.seed(123)
results <- FishPar(data = lf_data, resample = 1000)

12.2 Choosing the Number of Iterations

The default resample = 1000 is suitable for exploratory analysis and vignette building. For manuscript submission:

12.3 Providing Species-Specific Parameters

When species-specific L_inf or L_mat values are available from FishBase, tagging studies, or otolith ageing, they should be provided via the Linf and Lmat arguments to override the allometric regression estimates:

# Species: Glossogobius giuris; Linf from FishBase = 21.0 cm
results_specific <- FishPar(
  data        = lf_data,
  resample    = 5000,
  Linf        = 21.0,   # FishBase value
  Linf_sd     = 1.0,    # uncertainty in the FishBase estimate
  save_output = TRUE
)

When Linf is supplied, the function adds Gaussian noise with SD = Linf_sd at each MC iteration, so uncertainty in the literature estimate is still propagated.

12.4 LM_ratio and the Cope & Punt Gap

The Cope & Punt (2009) framework accepts only two LM_ratio categories: ≤ 0.75 or ≥ 0.90. If the estimated LM_ratio falls between 0.75 and 0.90, FishSS returns a warning and NULL. In this case the user should:

1. Check whether a species-specific L_mat from the literature would place LM_ratio clearly in one category.
2. Run sensitivity analyses with LM_ratio set to both 0.75 and 0.90 to bracket the stock status assessment.