## ----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")) ## ----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") ## ----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)) } } ## ----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_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) ## ----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 ## ----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") ## ----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) ## ----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]] ## ----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") ## ----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 # ) ## ----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 ## ----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)) ## ----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 ## ----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 ## ----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) ## ----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 ## ----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) ## ----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_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 )) ## ----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_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") ## ----seed_example, eval=FALSE------------------------------------------------- # set.seed(123) # results <- FishPar(data = lf_data, resample = 1000) ## ----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 # )