--- title: "aLBI: Length-Based Indicators and Fish Stock Assessment in R" subtitle: "A Comprehensive Guide to the aLBI Package" author: "Ataher Ali" date: "`r Sys.Date()`" vignette: > %\VignetteIndexEntry{aLBI: Length-Based Indicators and Fish Stock Assessment in R} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} output: rmarkdown::html_vignette: toc: true toc_depth: 4 number_sections: true --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, warning = FALSE, message = FALSE, cache = FALSE ) options(repos = c(CRAN = "https://cran.rstudio.com")) ``` *** # Introduction ## What is aLBI? `aLBI` (**A**ssessment of **L**ength-**B**ased **I**ndicators) 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**. ## 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`): | Component | Change | |---|---| | `FishPar` — Three-tier uncertainty | **Novel** propagation of L_mat / L_opt uncertainty to P_mat, P_opt, P_mega; separate Bootstrap-only CI also reported | | `FishPar` — `save_output` | Single toggle (default `FALSE`) replaces previous auto-saving | | `FishPar` — Excel output | Six sheets now, including a dedicated `FishSS_Inputs` sheet | | `FishPar` — Visualisation | Zone-shaded 6-panel annotation; dumbbell Target vs. Catch; 9-panel three-tier histograms | | `FishSS` — Trigger logic | Robust `if (Pobj >= 200) Popt else Pmat` replaces fragile table-row check | | `FishSS` — Return value | Now returns `Pobj`, `Px_trigger`, `Px_value`, `StockStatus`, `Selectivity` | | `FreqTM` | Multi-month frequency tables with consistent bin structure | | `LWR` | Length–weight relationship fitting and visualisation | ## 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. > *** # Installation ```{r install, eval=FALSE} # Install from CRAN (stable release) install.packages("aLBI") # Install the latest development version from GitHub # install.packages("devtools") # required devtools::install_github("Ataher76/aLBI") ``` ### Required Dependencies `aLBI` depends on base R packages (`stats`, `graphics`, `grDevices`, `utils`). The `openxlsx` package is required only when `save_output = TRUE` in `FishPar`. ```{r load_packages, message=FALSE, warning=FALSE} 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)) } } ``` *** # Package Functions Overview | Function | Purpose | Key Input | Key Output | |---|---|---|---| | `FrequencyTable` | Length-frequency table from raw measurements | Individual lengths | Binned frequency table | | `FreqTM` | Length-frequency tables across multiple months | Monthly length data | Per-month frequency tables | | `FishPar` | Length parameters + Froese indicators with three-tier uncertainty | Length-frequency table | Parameters, three-tier CIs, plots | | `FishSS` | Stock status assessment via Cope & Punt (2009) | FishPar outputs + CPdata | Stock status probabilities | | `LWR` | Length–weight relationship | Length & weight data | Regression model + plot | The typical workflow is sequential: `FrequencyTable` → `FishPar` → `FishSS` *** # Data Requirements and Preparation ## Input Data Formats `aLBI` functions accept data in two general formats: **Raw individual measurements** (for `FrequencyTable`, `FreqTM`, `LWR`): | Length | (Weight) | |---|---| | 12.3 | 24.1 | | 11.8 | 21.7 | | ... | ... | **Aggregated length-frequency data** (for `FishPar`): | Length | Frequency | |---|---| | 10 | 45 | | 11 | 102 | | 12 | 138 | | ... | ... | Data may be loaded from `.xlsx`, `.csv`, or created directly in R. The package ships with example datasets in `inst/exdata/`. ```{r list_data_files, eval=FALSE} # 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 ``` *** # FrequencyTable: Generating Length-Frequency Distributions ## 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. ## Function Arguments | Argument | Type | Default | Description | |---|---|---|---| | `data` | Data frame or vector | — | Individual fish length measurements | | `bin_width` | Numeric | `NULL` | Fixed bin width (cm). If `NULL`, Wang's OBS formula is used | | `Lmax` | Numeric | `NULL` | Maximum expected length. If `NULL`, observed maximum is used | | `output_file` | Character | `"FrequencyTable_Output.xlsx"` | Name of saved Excel file | ## Usage Example ```{r FrequencyTable_load, eval=TRUE} 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") head(length_data, 8) ``` ```{r FrequencyTable_run, eval=FALSE} # 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 ``` ```{r FrequencyTable_manual_example, eval=TRUE} # 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") ``` ## 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. *** # FreqTM: Multi-Month Length-Frequency Tables ## 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. ## Function Arguments | Argument | Type | Default | Description | |---|---|---|---| | `data` | Data frame | — | Columns for month and individual fish lengths (any column names) | | `bin_width` | Numeric | `NULL` | Bin width. If `NULL`, Wang's OBS formula is applied to the full dataset | | `Lmax` | Numeric | `NULL` | Maximum expected length | | `date_config` | List | `list(day=1, year=2025)` | Sets day and year for converting month names to dates | | `output_file` | Character | `"FreqTM_Output.xlsx"` | Excel output file name | ## Usage Example ```{r FreqTM_load, eval=TRUE} # 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") cat("Total observations:", nrow(monthly_data), "\n") head(monthly_data, 6) ``` ```{r FreqTM_run, eval=FALSE} # 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]] ``` ## 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. *** # FishPar: Length-Based Indicators with Three-Tier Uncertainty ## 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. ## Monte Carlo Simulation for Length Parameters ### From L_max to L_inf In data-limited contexts, the asymptotic growth length *L*_∞ cannot be estimated from age-at-length or tagging data. Instead, it is approximated from the maximum observed length in the catch via: $$L_{\infty} = L_{\max} / 0.95$$ This reflects the empirical observation that the largest fish in a well-sampled catch attains approximately 95% of the asymptotic length (Froese & Pauly 2020). To propagate uncertainty from incomplete sampling of the largest size classes, *L*_max at each MC iteration is drawn from a **truncated Gaussian distribution**: $$L_{\max,i} \sim \mathcal{N}(\bar{L}_{\max},\ 0.05 \cdot \bar{L}_{\max}), \quad \text{bounded to } [0.90 \cdot \bar{L}_{\max},\ 1.10 \cdot \bar{L}_{\max}]$$ ### From L_inf to L_mat Length at first maturity *L*_mat is estimated using the cross-species allometric regression of Froese & Binohlan (2000): $$\log_{10}(L_{\text{mat},i}) = 0.8979 \cdot \log_{10}(L_{\infty,i}) - 0.0782 + \varepsilon_{1,i}$$ where $\varepsilon_{1,i} \sim \mathcal{N}(0,\ 0.015)$ represents regression residual uncertainty in log₁₀ space. ### From L_mat to L_opt The optimal fishing length *L*_opt — the body length at which yield per recruit is maximised — is estimated as: $$\log_{10}(L_{\text{opt},i}) = 1.053 \cdot \log_{10}(L_{\text{mat},i}) - 0.0565 + \varepsilon_{2,i}$$ where $\varepsilon_{2,i} \sim \mathcal{N}(0,\ 0.015)$. Both regressions were derived from 248 fish species and exhibit $r^2 \geq 0.98$ (Froese & Binohlan 2000). The boundaries of the **optimal size range** follow Froese (2004): $$L_{\text{opt\_m10}} = 0.90 \cdot L_{\text{opt}}; \qquad L_{\text{opt\_p10}} = 1.10 \cdot L_{\text{opt}}$$ Point estimates for all six parameters are the means of the MC distributions; **95% confidence intervals** are the 2.5th and 97.5th percentiles. ## 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. ### The Composite Objective Index (P_obj) Following Cope & Punt (2009): $$P_{\text{obj}} = P_{\text{mat}} + P_{\text{opt}} + P_{\text{mega}}$$ This index ranges from 0 to 300 and governs the trigger indicator and look-up table selection in `FishSS`: | P_obj range | Trigger indicator (P_x) | Interpretation | |---|---|---| | < 100 | P_mat | Predominantly juvenile/sub-optimal catch | | 100 – 200 | P_mat | Catch spans the maturity ogive | | ≥ 200 | P_opt | Predominantly mature and optimally-sized catch | ## 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. ### The Step-Function Problem P_mat is not a smooth function of L_mat — it is a **step function** that jumps discontinuously whenever L_mat crosses a length class boundary. For example, if the length class width is 1 cm and there are 158 fish in the 12-cm class, then: - L_mat = 11.9 cm → P_mat = 61.9% - L_mat = 12.1 cm → P_mat = 43.1% A 0.2 cm shift in L_mat causes an 18.7 percentage-point jump in P_mat. The MC distribution of L_mat will inevitably straddle such boundaries, so the expected P_mat across all MC iterations will differ from the P_mat computed at the mean L_mat alone. This is why **model outputs systematically differ from manual calculations**, and why the difference is methodologically correct, not an error. ### The Three Tiers | Tier | Parameters | Catch data | Uncertainty captured | |---|---|---|---| | **Tier 1** — Parameter only | MC-varying (L_mat_i, L_opt_i) | Observed (fixed) | Parameter estimation uncertainty only | | **Tier 2** — Data only | Fixed at MC means | Bootstrap resampled | Catch-data sampling variability only | | **Tier 3** — Total (novel) | MC-varying | Bootstrap resampled | **Total = Parameter + Data** | At each of the B iterations, the loop computes all three tiers simultaneously: **Tier 1** (i-th MC parameters × observed catch): $$I_i^{(1)} = f(\mathbf{x}_{\text{obs}},\ \theta_i)$$ **Tier 2** (fixed mean parameters × bootstrap resample): $$I_i^{(2)} = f(\mathbf{x}_i^*,\ \bar\theta)$$ **Tier 3** (i-th MC parameters × bootstrap resample): $$I_i^{(3)} = f(\mathbf{x}_i^*,\ \theta_i)$$ Under approximate parameter–data independence, the total variance satisfies: $$\text{Var}(I^{(3)}) \approx \text{Var}(I^{(1)}) + \text{Var}(I^{(2)})$$ All 95% CIs are computed as the 2.5th–97.5th percentiles of the B-iteration distributions — no normality assumption is imposed. > **Reporting recommendation:** Report Tier 3 as the primary CI in manuscripts > (total uncertainty, most conservative). Tier 2 is appropriate only when length > parameters are independently known. Tier 1 is most useful as a diagnostic showing > the isolated sensitivity of indicators to allometric regression uncertainty. ## Function Arguments | Argument | Type | Default | Description | |---|---|---|---| | `data` | Data frame | — | Two columns: `Length` and `Frequency` | | `resample` | Integer | `1000` | Number of MC / bootstrap iterations. Use ≥ 5000 for publication | | `progress` | Logical | `FALSE` | Display text progress bar | | `Linf` | Numeric | `NULL` | Known L_inf (overrides L_max/0.95 estimate) | | `Linf_sd` | Numeric | `0.5` | Gaussian noise SD added to each L_inf sample (cm) | | `Lmat` | Numeric | `NULL` | Known L_mat (overrides regression estimate) | | `Lmat_sd` | Numeric | `0.5` | Gaussian noise SD added to each L_mat sample (cm) | | `save_output` | Logical | `FALSE` | Write all PDFs and Excel to `getwd()`. Set `TRUE` to save | ## Usage Example ```{r FishPar_load_data, eval=TRUE} 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") print(lf_data) cat("\nTotal individuals:", sum(lf_data[[2]]), "\n") ``` ```{r FishPar_run, eval=FALSE} # ── 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 ) ``` ## Accessing Results ```{r FishPar_access, eval=FALSE} # ── 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 ``` ### Expected Output Structure ``` $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 ``` ## Generated Plots When `save_output = TRUE`, `FishPar` writes ten PDF files to `getwd()`: | File | Description | |---|---| | `Length_Frequency_Plot.pdf` | Bar histogram + Gaussian KDE smooth + weighted mean length | | `Main_Graph_Annotations.pdf` | 6-panel zone-shaded LFD with each parameter as a vertical line | | `Length_Parameters_Histograms.pdf` | 6-panel MC distributions with density overlay | | `Length_Parameters_Density.pdf` | 6-panel MC kernel densities with shaded fills | | `Length_Parameters_CI.pdf` | Forest plot: mean ± 95% CI for all six parameters | | `Froese_Indicators_Histograms.pdf` | 9-panel histograms (3 tiers × 3 indicators) | | `Froese_Indicators_Density.pdf` | 9-panel kernel densities (3 tiers × 3 indicators) | | `Froese_Indicators_CI.pdf` | Forest plot: Tier 3 CI with Froese (2004) target diamonds | | `Froese_Uncertainty_3Tier.pdf` | Three-panel comparative forest plot (one panel per tier) | | `Target_vs_Catch_Dumbbell.pdf` | Dumbbell chart: observed vs. target with management status colours | ### Reading the Three-Tier Uncertainty Plot The `Froese_Uncertainty_3Tier.pdf` is the most informative diagnostic plot. It shows three panels side by side, one per tier, each displaying all three indicators (P_mat, P_opt, P_mega) on the x-axis. Each panel has an **independently scaled y-axis**, which is crucial for readability: - **Tier 1 (red)** and **Tier 3 (green)**: y-axis spans the full wide CI range (often nearly 0–100%), reflecting substantial parameter uncertainty. - **Tier 2 (blue)**: y-axis is **zoomed in** to the narrow CI range (typically ±1–3 percentage points), reflecting only the sampling precision of the catch data. The contrast between Tier 2 (tight) and Tier 3 (wide) visually communicates the dominant role of allometric regression uncertainty in total indicator uncertainty. ## The Excel Workbook (`FishPar_Results.xlsx`) The Excel workbook contains six sheets: | Sheet | Contents | |---|---| | `Length_Parameters` | MC mean ± 95% CI for all six length parameters | | `Froese_Tier1_ParOnly` | Tier 1 CI — parameter uncertainty only | | `Froese_Tier2_DataOnly` | Tier 2 CI — data sampling uncertainty (classical bootstrap) | | `Froese_Tier3_Joint` | Tier 3 CI — total uncertainty (recommended for reporting) | | `Target_vs_Catch` | Observed P_mat, P_opt, P_mega vs. Froese (2004) targets | | `FishSS_Inputs` | LM_ratio, P_mat, P_opt, P_mega, P_obj, Total_ind — ready for `FishSS()` | > **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. *** # FishSS: Stock Status Assessment ## 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) ### Column Selection Logic The look-up table has 10 probability columns (A–J) corresponding to different combinations of P_obj regime and LM_ratio category. In `aLBI`, the trigger indicator column (`Tx`) is converted to a percentage scale (0–100%, in steps of 5%) to maintain exact numerical alignment with `FishPar` indicator outputs: | P_obj | LM_ratio | Columns used | Trigger P_x | |---|---|---|---| | < 100 | ≤ 0.75 | A (target), C (limit) | P_mat | | < 100 | ≥ 0.90 | B (target), D (limit) | P_mat | | 100–200 | ≤ 0.75 | E (target), G (limit) | P_mat | | 100–200 | ≥ 0.90 | F (target), H (limit) | P_mat | | ≥ 200 | Any | I (target), J (limit) | P_opt | The trigger indicator P_x is the catch proportion used as the "trigger value" to locate the appropriate row in the selected columns: $$P_x = \begin{cases} P_{\text{mat}} & \text{if } P_{\text{obj}} < 200 \\ P_{\text{opt}} & \text{if } P_{\text{obj}} \geq 200 \end{cases}$$ > **Important fix in the current version:** Earlier versions used a fragile > `if (p[[1,2]] > 0)` check to determine the trigger indicator, which silently failed > when CPdata was not sorted in descending order by Tx. The current version uses the > direct rule above, which is robust to any row ordering of CPdata. ### Selectivity Classification Based on P_obj, `FishSS` classifies the catch selectivity pattern: | P_obj | Pattern | |---|---| | < 100, P_opt = P_mega = 0 | Fish small, immature | | < 100, P_opt or P_mega > 0 | Fish small and optimally-sized or all but biggest | | 100–200 | Fish maturity ogive | | ≥ 200, P_opt < 100 | Fish optimally-sized and bigger | | ≥ 200, P_opt = 100 | Fish optimally-sized | ## Function Arguments | Argument | Type | Description | |---|---|---| | `data` | Data frame | Cope & Punt (2009) look-up table (CPdata, columns: Tx, A–J) | | `LM_ratio` | Numeric | L_mat / L_opt ratio (from `FishPar` results) | | `Pmat` | Numeric | Proportion of mature fish, % (from `FishPar`) | | `Popt` | Numeric | Proportion at optimal size, % (from `FishPar`) | | `Pmega` | Numeric | Proportion of mega-spawners, % (from `FishPar`) | ## Usage Example ```{r FishSS_load_cpdata, eval=TRUE} # 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") cat("Columns:", paste(colnames(cpdata), collapse = ", "), "\n") print(head(cpdata, 6)) ``` ```{r FishSS_run, eval=FALSE} # ── 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 ``` ### Expected Return Value ```{r FishSS_output_example, eval=FALSE} 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 ``` ## Interpreting Stock Status The two probability values returned in `$StockStatus` should be interpreted together: | `p_below_target` | `p_below_limit` | Likely stock status | |---|---|---| | < 50% | < 10% | **Healthy** — stock likely above target reference point | | 50–80% | 10–40% | **Cautionary** — stock may be below target; increase monitoring | | > 80% | > 40% | **Overfished** — stock likely below target; reduce fishing pressure | | > 95% | > 70% | **Depleted** — stock likely at or below limit reference point | > **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. ### Feeding Uncertainty into FishSS The current `FishSS` uses **point estimates** of P_mat, P_opt, and P_mega. To propagate the full three-tier uncertainty through to the stock status probabilities, one can run `FishSS` multiple times using values sampled from the Tier 3 bootstrap distribution. This is facilitated by the `froese_par_joint` matrix returned by `FishPar`. *** # LWR: Length–Weight Relationship ## 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. ## Function Arguments | Argument | Type | Default | Description | |---|---|---|---| | `data` | Data frame | — | Two columns: length (col 1) and weight (col 2) | | `log_transform` | Logical | `TRUE` | Apply log–log transformation for linearisation | | `point_col` | Character | `"black"` | Data point colour | | `line_col` | Character | `"red"` | Regression line colour | | `shade_col` | Character | `"red"` | Confidence interval ribbon colour | | `point_size` | Numeric | `2` | Data point size | | `line_size` | Numeric | `1` | Regression line thickness | | `alpha` | Numeric | `0.2` | Confidence interval ribbon transparency | | `main` | Character | `"Length-Weight Relationship"` | Plot title | | `xlab` | Character | `NULL` | Custom x-axis label (auto-generated if NULL) | | `ylab` | Character | `NULL` | Custom y-axis label (auto-generated if NULL) | | `save_output` | Logical | `FALSE` | Save plot (PDF) and model summary (TXT) to `getwd()` | ## Usage Example ```{r LWR_load, eval=TRUE} # 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") cat("Length range:", round(min(LWdata[[1]]), 1), "–", round(max(LWdata[[1]]), 1), "cm\n") cat("Weight range:", round(min(LWdata[[2]]), 2), "–", round(max(LWdata[[2]]), 2), "g\n") head(LWdata) ``` ```{r LWR_run, eval=FALSE} # 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 ``` ### Output Annotation The LWR plot is annotated with: - **Regression equation** in back-transformed form: W = a · L^b - **R²** — proportion of variance in weight explained by length - **p-value** — significance of the relationship - **95% confidence interval ribbon** around the fitted line *** # Integrated Workflow: From Raw Data to Stock Status This section demonstrates the complete `aLBI` workflow from raw length measurements to a stock status assessment. ```{r workflow_overview, eval=FALSE} # ───────────────────────────────────────────────────────────────────────── # STEP 1: Build the length-frequency table from raw measurements # ───────────────────────────────────────────────────────────────────────── length_data <- readxl::read_excel("my_length_data.xlsx") freq_result <- FrequencyTable( data = length_data, bin_width = 1, # 1-cm bins (or NULL for Wang's OBS formula) save_output = FALSE ) # Extract the frequency table for FishPar lf_table <- freq_result$lfreq print(lf_table) # ───────────────────────────────────────────────────────────────────────── # STEP 2: Estimate length parameters and Froese indicators # ───────────────────────────────────────────────────────────────────────── results <- FishPar( data = lf_table, resample = 5000, # 5000 iterations for stable CIs save_output = FALSE # writes FishPar_Results.xlsx + all PDFs ) # Review key outputs cat("=== Length Parameters ===\n"); print(results$estimated_length_par) cat("\n=== Froese Indicators (Tier 3 — Total Uncertainty) ===\n") print(results$froese_par_joint) cat("\nPobj:", round(results$Pobj, 2), "(0–300 scale)\n") cat("LM_ratio:", round(results$LM_ratio, 3), "\n") # ───────────────────────────────────────────────────────────────────────── # STEP 3: Assess stock status # ───────────────────────────────────────────────────────────────────────── cpdata <- readxl::read_excel( system.file("exdata", "cpdata.xlsx", package = "aLBI") ) # Use Tier 3 means as point estimates for FishSS stock_status <- FishSS( data = cpdata, LM_ratio = results$LM_ratio, Pmat = results$froese_par_joint$Mean[1], Popt = results$froese_par_joint$Mean[2], Pmega = results$froese_par_joint$Mean[3] ) cat("\n=== Stock Status Assessment ===\n") cat("P_obj:", stock_status$Pobj, "\n") cat("Trigger indicator:", stock_status$Px_trigger, "(value =", round(stock_status$Px_value, 1), "%)\n") cat("Selectivity pattern:", stock_status$Selectivity, "\n") cat("P(biomass < 0.40 SB0):", stock_status$StockStatus["p_below_target"], "%\n") cat("P(biomass < 0.25 SB0):", stock_status$StockStatus["p_below_limit"], "%\n") # ───────────────────────────────────────────────────────────────────────── # STEP 4 (optional): Length-weight relationship # ───────────────────────────────────────────────────────────────────────── lw_data <- readxl::read_excel("my_lw_data.xlsx") lwr_result <- LWR(data = lw_data, log_transform = TRUE, save_output = FALSE) ``` *** # Why Model Outputs Differ from Manual Calculations ## 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. ## Worked Example Consider the following dataset (N = 844 fish, bin width = 1 cm): ```{r manual_example_data, eval=TRUE} # 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") print(example_lf) ``` ### Manual Calculation (Fixed Parameters) ```{r manual_calc, eval=TRUE} # 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 )) # 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 )) ``` ### Why the Step-Function Matters ```{r step_function_demo, eval=TRUE} # 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") 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) # 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 )) cat("A 0.2 cm shift in Lmat (11.9 → 12.1) changes Pmat by", jump_12cm, "pp\n") ``` ### Comparison: Manual vs. Model (Tier 3) ```{r comparison_table, eval=TRUE} # 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) cat("\nKey insight: Tier 2 ≈ Manual because both use the same fixed Lmat/Lopt.\n") cat("Tier 3 differs because it integrates over the MC distribution of Lmat,\n") cat("propagating the 2.64 cm uncertainty range [10.45, 13.09] through the\n") cat("step function, which shifts the expected Pmat from 61.85% to 56.47%.\n") ``` ## 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. *** # Reproducibility and Best Practices ## 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`: ```{r seed_example, eval=FALSE} set.seed(123) results <- FishPar(data = lf_data, resample = 1000) ``` ## Choosing the Number of Iterations The default `resample = 1000` is suitable for exploratory analysis and vignette building. For manuscript submission: | Purpose | Recommended `resample` | |---|---| | Exploration and diagnostics | 1 000 | | Preliminary results | 2 000 | | Manuscript submission | **5 000** | | High-precision intervals | 10 000 | ## 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: ```{r user_params, eval=FALSE} # 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. ## 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. *** # Conclusion The `aLBI` package provides a rigorous, reproducible framework for data-limited fish stock assessment from catch length-frequency data alone. Its key contributions are: - **Accessible implementation** of the Froese (2004) and Cope & Punt (2009) frameworks in a single, well-documented R package. - **Monte Carlo uncertainty propagation** for all six biological reference lengths, with percentile-based 95% CIs that make no distributional assumptions. - **Three-tier uncertainty decomposition** (Tier 1: parameter only; Tier 2: data only; Tier 3: total/joint) — the first implementation to explicitly propagate allometric regression uncertainty through to the Froese sustainability indicators. - **Comprehensive visualisation suite** (10 publication-ready PDFs) with coherent colour coding, zone annotations, and independent y-axis scaling in comparative plots. - **Seamless workflow integration**: the `FishSS_Inputs` Excel sheet eliminates manual value transfer between `FishPar` and `FishSS`. *** # Acknowledgements The author expresses sincere gratitude to **Dr. Mohammed Shahidul Alam** for his unwavering guidance, expert supervision, and continuous encouragement throughout the development of this package. Heartfelt thanks are also extended to the reviewers, collaborators, and open-source R community members whose feedback has substantially improved both the methodology and the implementation. *** # References 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. Cope, J. M., & Punt, A. E. (2009). Length-based reference points for data-limited situations: Applications and restrictions. *Marine and Coastal Fisheries*, 1(1), 169–186. Efron, B., & Tibshirani, R. J. (1993). *An Introduction to the Bootstrap*. Chapman & Hall, New York. Froese, R. (2004). Keep it simple: Three indicators to deal with overfishing. *Fish and Fisheries*, 5(1), 86–91. Froese, R., & Binohlan, C. (2000). Empirical relationships to estimate asymptotic length, length at first maturity and length at maximum yield per recruit in fishes. *Journal of Fish Biology*, 56, 758–773. Froese, R., & Pauly, D. (Eds.) (2020). *FishBase*. World Wide Web electronic publication. Silverman, B. W. (1986). *Density Estimation for Statistics and Data Analysis*. Chapman & Hall, London. Wang, K., Zhang, C., Xu, B., Xue, Y., & Ren, Y. (2020). Selecting optimal bin size to account for growth variability in Electronic LEngth Frequency ANalysis (ELEFAN). *Fisheries Research*, 225, 105474. *** # Contact | | | |---|---| | **Author** | Ataher Ali | | **E-mail** | ataher.cu.ms@gmail.com | | **Package** | [github.com/Ataher76/aLBI](https://github.com/Ataher76/aLBI) | | **CRAN** | [cran.r-project.org/package=aLBI](https://cran.r-project.org/package=aLBI) | *Please report bugs and feature requests via the GitHub Issues page.*