1 Overview

FLicc is an user-friendly TMB implementation in the FLR and a merged generalization of both the multi-gear length-interval catch-curve framework developed in fishblicc (Medley 2025) and the length-based spawning potential ratio estimation in LBSRP (Hordyk et al. 2015. 2016).

FLicc retains the biological core of either of the original length-based equilibrium models, but adds a fast penalized-likelihood workflow in Template Model Builder (TMB), multi-year fitting, FLR integration, and a growing set of equilibrium and indicator tools. The original multi-gear formulation was designed to estimate mortality-at-length, selectivity, spawning potential ratio (SPR), yield-per-recruit from length compositions grouped by gear, and parallel computing utilities.

An overview with illustrations of the FLicc utilities is also provided here.

2 Installation

FLicc depends on FLCore, ggplotFL and TMB. The FLR packages are available from the FLR r-universe repository, while TMB is installed from CRAN. On Windows, compiling the TMB model requires Rtools.

2.1 Install required packages

install.packages("devtools")
install.packages("TMB", type = "source")
install.packages("ggplot2")

TMB is the estimation engine used by FLicc. Installing from source can help ensure compatibility when compiling and linking the local TMB model.

2.2 Windows: install Rtools if needed

On Windows, TMB requires Rtools. If it is not already installed, open the Rtools download page from R:

browseURL("https://cran.r-project.org/bin/windows/Rtools/")

After installing Rtools, restart R before installing or compiling FLicc.

2.3 Install FLR dependencies

install.packages(
  c("FLCore", "ggplotFL"),
  repos = c("https://flr.r-universe.dev", "https://cloud.r-project.org")
)

2.4 Install FLicc

devtools::install_github("henning-winker/FLicc")

2.5 Load the package

library(FLicc)

3 Basic functionality

3.1 Getting started

FLicc implements a length-based catch-curve assessment model for data-limited fisheries. The package is designed to work with the FLR data classes and provides tools to move from length-frequency observations to fitted selectivity, fishing mortality, spawning-potential ratio (SPR), length-based indicators and simple advice plots.

The core workflow is:

  1. prepare length-frequency data as FLQuants;
  2. specify life-history parameters;
  3. build an FLStockLen object with stocklen();
  4. fit the model using fiticc();
  5. inspect convergence, fit diagnostics, selectivity and SPR;
  6. derive equilibrium reference points and length-based advice indicators.

This vignette uses the built-in splendid alfonsino example and mirrors the main steps in test.R. It focuses on the single-fit workflow. Monte-Carlo ensemble fitting is not covered here.

library(FLicc)

4 Loading and inspecting length-frequency data

4.1 Built-in alfonsino data

The package contains an example length-frequency dataset for splendid alfonsino. The object lfd_alfonsino is stored as length-frequency data in FLR format.

data("alfonsino")

lfd_alfonsino
#> $ Trawl 
#> An object of class "FLQuant"
#> , , unit = unique, season = all, area = unique
#> 
#>     year
#> len  2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024
#>   10  0    0    0    0    0    0    0    0    0    0    0    0  
#>   11  1    0    1    1    0    2    0    0    0    0    0    0  
#>   12  2    3    1    5    4    2    5    8    2    2    3    2  
#>   13  7    8    7    7   18   11   10   11   10   11    7    9  
#>   14 17   13   19   28   24   21   32   20   18   29   23   34  
#>   15 36   32   25   48   37   50   31   41   37   45   29   60  
#>   16 53   38   32   82   58   42   58   59   41   58   49   61  
#>   17 59   37   48   68   62   43   60   52   33   61   67   67  
#>   18 66   45   55   51   67   57   52   45   50   56   71   71  
#>   19 41   48   51   51   54   62   54   47   52   59   52   55  
#>   20 57   42   35   32   45   41   64   39   45   54   48   36  
#>   21 54   36   38   24   30   41   32   37   57   37   44   40  
#>   22 28   33   39   27   21   38   33   30   29   26   35   24  
#>   23 22   36   33   20   24   31   15   34   32   19   20   17  
#>   24 21   28   32   11   19   19   19   20   22   12   18   13  
#>   25 20   19   22   15   13   10   10    6   17   11   10    7  
#>   26  7   22   20   16    8   10    5    9   18   11    8    2  
#>   27  5   20   14    5    6    8    7   10   14    3    7    2  
#>   28  6   17   10    2    5    2    3    7    4    4    4    1  
#>   29  1    6    7    3    2    4    2    6    8    3    5    3  
#>   30  0    6    5    1    3    3    3    2    2    1    1    0  
#>   31  1    5    5    2    0    3    1    4    1    0    2    0  
#>   32  0    4    2    0    0    3    1    2    4    1    0    0  
#>   33  0    1    2    1    1    0    3    0    1    0    0    0  
#>   34  0    1    1    0    0    0    1    0    3    1    0    0  
#>   35  0    1    0    1    1    0    2    1    2    0    0    0  
#>   36  0    0    0    0    0    1    0    1    1    0    0    0  
#>   37  0    2    0    1    0    0    0    1    0    0    0    0  
#>   38  0    1    0    1    1    0    0    1    0    0    0    0  
#>   39  0    0    0    0    1    0    0    0    0    0    0    0  
#>   40  0    0    0    1    0    0    0    2    0    0    0    0  
#>   41  0    0    0    0    0    0    0    1    1    0    0    0  
#>   42  0    0    0    0    0    0    0    1    0    0    0    0  
#>   43  0    0    0    0    0    0    1    2    0    0    1    0  
#>   44  0    0    0    0    0    0    0    1    0    0    0    0  
#>   45  0    0    0    0    0    0    0    0    0    0    0    0  
#>   46  0    0    0    0    0    0    0    2    0    0    0    0  
#>   47  0    0    0    0    0    0    0    2    0    0    0    0  
#>   48  0    0    0    0    0    0    0    0    0    0    0    0  
#>   49  0    0    0    0    0    0    0    0    0    0    0    0  
#>   50  0    0    0    0    0    0    0    0    0    0    0    0  
#>   51  0    0    0    0    0    0    0    0    0    0    0    0  
#>   52  0    0    0    0    0    0    0    0    0    0    0    0  
#>   53  0    0    0    0    0    0    0    0    0    0    0    0  
#>   54  0    0    0    0    0    0    0    0    0    0    0    0  
#>   55  0    0    0    0    0    0    0    0    0    0    0    0  
#>   56  0    0    0    0    0    0    0    0    0    0    0    0  
#>   57  0    0    0    0    0    0    0    0    0    0    0    0  
#>   58  0    0    0    0    0    0    0    0    0    0    0    0  
#>   59  0    0    0    0    0    0    0    0    0    0    0    0  
#>   60  0    0    0    0    0    0    0    0    0    0    0    0  
#>   61  0    0    0    0    0    0    0    0    0    0    0    0  
#>   62  0    0    0    0    0    0    0    0    0    0    0    0  
#>   63  0    0    0    0    0    0    0    0    0    0    0    0  
#>   64  0    0    0    0    0    0    0    0    0    0    0    0  
#>   65  0    0    0    0    0    0    0    0    0    0    0    0  
#>   66  0    0    0    0    0    0    0    0    0    0    0    0  
#>   67  0    0    0    0    0    0    0    0    0    0    0    0  
#> 
#> units:  cm 
#> 
#> $ Gillnet 
#> An object of class "FLQuant"
#> , , unit = unique, season = all, area = unique
#> 
#>     year
#> len  2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024
#>   10  0    0    0    0    0    0    0    0    0    0    0    0  
#>   11  0    0    0    0    0    0    0    0    0    0    0    0  
#>   12  0    0    0    0    0    0    0    0    0    0    0    0  
#>   13  0    0    0    0    0    0    0    0    0    0    0    0  
#>   14  0    0    0    0    0    0    0    0    0    0    0    0  
#>   15  0    0    0    0    0    0    0    0    0    0    0    0  
#>   16  0    0    0    0    0    0    1    0    0    0    0    0  
#>   17  1    0    0    0    0    0    2    0    0    0    0    1  
#>   18  0    0    0    0    0    0    0    0    0    0    0    0  
#>   19  1    1    1    1    1    0    2    1    1    0    0    3  
#>   20  5    1    1    1    1    2    4    0    1    3    2    5  
#>   21  6    6    3    6    4    2    6    4    4    5    2    9  
#>   22  9   11    8   10    5    4    9    5    6   11    6    8  
#>   23 20   16    6   17    5    5   10   21    8    8    7   14  
#>   24 22   20   20   16    7    8   15   14    9   15   15    9  
#>   25 32   39   35   18   18   14   13   33   21   23   18   18  
#>   26 42   38   35   42   20   24   28   25   22   29   19   15  
#>   27 42   45   45   35   21   30   22   23   19   24   37   29  
#>   28 28   50   50   29   31   25   35   27   28   33   34   35  
#>   29 34   37   29   40   34   25   33   37   27   35   46   36  
#>   30 33   41   36   36   37   46   24   27   32   34   24   24  
#>   31 21   37   31   43   33   41   29   38   32   31   33   35  
#>   32 21   24   40   28   38   36   30   25   39   34   27   40  
#>   33 27   26   19   33   57   36   37   26   20   28   19   28  
#>   34 22   14   17   17   31   38   32   23   33   27   15   15  
#>   35 16   18   17   17   26   39   26   18   32   34   24   29  
#>   36 16   12   12   19   24   17   27   21   29   25   20   18  
#>   37 22    6    4   22   19   15   22   28   25   16   14   13  
#>   38  9   20   16    8   12   23   23   25   28   15   20   25  
#>   39  9    8    7   13   11   13   12   13   18   13   11   13  
#>   40 13    6   12    5   15   16    6   15   19   17   15   15  
#>   41  7    2    8    5    7   12   13   10   11   12   13   16  
#>   42 10    2    8    4    6    6   12   10    8    8   19    5  
#>   43  4    5    3    3    5    6    4   10   10    4   13    7  
#>   44  6    2    5    4    6    5    2    6    5    8   11    6  
#>   45  5    4    5    5    4    3    5    2    6    1   12    8  
#>   46  2    1    4    6    7    2    4    2    4    3    7    6  
#>   47  2    2    9    6    3    2    5    0    2    1    4    1  
#>   48  5    0    4    2    3    1    2    7    0    0    4    4  
#>   49  4    0    3    1    4    2    1    3    0    0    2    3  
#>   50  1    1    1    2    1    1    1    1    1    1    2    2  
#>   51  1    2    4    3    1    1    2    0    0    1    0    1  
#>   52  1    2    0    0    0    0    0    0    0    0    2    1  
#>   53  0    0    0    1    2    0    0    0    0    1    1    0  
#>   54  1    0    1    1    0    0    1    0    0    0    1    2  
#>   55  0    1    0    0    1    0    0    0    0    0    1    1  
#>   56  0    0    0    1    0    0    0    0    0    0    0    0  
#>   57  0    0    0    0    0    0    0    0    0    0    0    0  
#>   58  0    0    0    0    0    0    0    0    0    0    0    0  
#>   59  0    0    1    0    0    0    0    0    0    0    0    0  
#>   60  0    0    0    0    0    0    0    0    0    0    0    0  
#>   61  0    0    0    0    0    0    0    0    0    0    0    0  
#>   62  0    0    0    0    0    0    0    0    0    0    0    0  
#>   63  0    0    0    0    0    0    0    0    0    0    0    0  
#>   64  0    0    0    0    0    0    0    0    0    0    0    0  
#>   65  0    0    0    0    0    0    0    0    0    0    0    0  
#>   66  0    0    0    0    0    0    0    0    0    0    0    0  
#>   67  0    0    0    0    0    0    0    0    0    0    0    0  
#> 
#> units:  cm

A useful first step is to inspect the data as a long data.frame. This is especially important when moving between plain tabular input and FLQuants, because the model expects length, year and gear structure to be explicit.

lfd_long <- as.data.frame(lfd_alfonsino)
head(lfd_long)
#>   len year   unit season   area iter data qname
#> 1  10 2013 unique    all unique    1    0 Trawl
#> 2  11 2013 unique    all unique    1    1 Trawl
#> 3  12 2013 unique    all unique    1    2 Trawl
#> 4  13 2013 unique    all unique    1    7 Trawl
#> 5  14 2013 unique    all unique    1   17 Trawl
#> 6  15 2013 unique    all unique    1   36 Trawl

The long structure should contain the length bin, year, gear or unit, and observed frequency or count. The exact column names depend on the FLR object, but the key point is that each row represents one length-bin observation for a given year and gear.

unique(lfd_long$unit)
#> [1] unique
#> Levels: unique
range(lfd_long$len)
#> [1] 10 67
range(lfd_long$year)
#> [1] 2013 2024

4.2 From long input to wide input and back to FLQuants

A common input workflow starts from a long length-frequency table. lfd_long_to_wide() converts the long table into the wide format expected by FLQuantLen().

LFD.df <- lfd_long_to_wide(as.data.frame(lfd_alfonsino))

head(LFD.df)
#>   len  gear 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024
#> 1  10 Trawl    0    0    0    0    0    0    0    0    0    0    0    0
#> 2  11 Trawl    1    0    1    1    0    2    0    0    0    0    0    0
#> 3  12 Trawl    2    3    1    5    4    2    5    8    2    2    3    2
#> 4  13 Trawl    7    8    7    7   18   11   10   11   10   11    7    9
#> 5  14 Trawl   17   13   19   28   24   21   32   20   18   29   23   34
#> 6  15 Trawl   36   32   25   48   37   50   31   41   37   45   29   60
unique(LFD.df$gear)
#> [1] "Trawl"   "Gillnet"
unique(LFD.df$len)
#>  [1] 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
#> [26] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
#> [51] 60 61 62 63 64 65 66 67

The wide object is then converted to an FLQuants length object. Here midL = FALSE indicates that the length values represent lower length-bin boundaries rather than mid-points.

lfd <- FLQuantLen(LFD.df, unit = "cm", midL = FALSE)

The resulting object can again be displayed as a long table. This is a useful check that the FLQuants object contains the intended length, year and gear structure.

head(as.data.frame(lfd))
#>   len year   unit season   area iter data qname
#> 1  10 2013 unique    all unique    1    0 Trawl
#> 2  11 2013 unique    all unique    1    1 Trawl
#> 3  12 2013 unique    all unique    1    2 Trawl
#> 4  13 2013 unique    all unique    1    7 Trawl
#> 5  14 2013 unique    all unique    1   17 Trawl
#> 6  15 2013 unique    all unique    1   36 Trawl

4.3 Plotting observed length frequencies

Before fitting, always inspect the observed length distributions. plot_lfd() can show either relative frequencies or frequencies normalised to the maximum value.

plot_lfd(lfd, type = "relmax")

plot_lfd(lfd)

Large differences among gears are expected in multi-gear fisheries and are often informative. For example, one gear may mainly sample smaller fish, while another gear samples larger fish or the descending limb of a dome-shaped selectivity pattern.

5 Life-history input

5.1 Defining life-history parameters

FLicc requires life-history parameters for growth, natural mortality, maturity and weight-at-length. These are supplied as an FLPar object.

lhpar <- FLPar(
  linf = 55.7,
  k    = 0.08,
  M    = 0.162,
  L50  = 31.1859,
  a    = 0.004721956 / 1000,
  b    = 3.146168
)

lhpar
#> An object of class "FLPar"
#> params
#>     linf        k        M      L50        a        b 
#> 5.57e+01 8.00e-02 1.62e-01 3.12e+01 4.72e-06 3.15e+00 
#> units:  NA

The main quantities are:

  • linf: asymptotic length \(L_{\infty}\);
  • k: von Bertalanffy growth coefficient;
  • M: natural mortality;
  • L50: length at 50% maturity;
  • a and b: weight-length parameters in \(W(L) = aL^b\).

5.2 Biological equations

The weight-at-length relationship is

\[ W(L) = aL^b. \]

Maturity-at-length is represented by a logistic ogive,

\[ P_{mat}(L) = \frac{1}{1 + \exp[-\eta(L - L_{50})]}, \]

where \(L_{50}\) is the length at 50% maturity and \(\eta\) controls the slope of the maturity curve. In many applications the maturity curve is treated as fixed input rather than estimated from the length data.

Growth enters the length-based population model through the transition or growth formulation selected in settings. FLicc currently supports a gamma population model and a gtg growth-transition model. The practical role of these options is described below.

6 Building the FLStockLen object

6.1 Constant natural mortality

The stocklen() function combines length-frequency data and life-history parameters into an equilibrium length-structured stock object.

stklen <- stocklen(lfd, lhpar, m_model = "constant")

stklen
#> An object of class "FLStockLen"
#> 
#> Name: FLicc stocklen 
#> Description:  
#> Quant: len 
#> Dims:  len   year    unit    season  area    iter
#>  58  12  1   1   1   1   
#> 
#> Range:  min  max minyear maxyear minfbar maxfbar 
#>  10  67  2013    2024    10  67

The stock object stores biological schedules, length bins and derived quantities needed by the TMB model.

plot_lw(stklen)

plot_mat(stklen)

plot_m(stklen)

6.2 Natural mortality options

The shape of natural mortality-at-length can be changed using m_model. The examples below compare the available options.

m_models <- c("constant", "inverse", "Lorenzen", "Gislason")

m_stks <- FLStocks(lapply(m_models, function(x) {
  stocklen(lfd, lhpar, m_model = x)
}))
names(m_stks) <- m_models

plot_m(m_stks)

The natural mortality schedule is an important structural assumption. A constant \(M\) assumes the same mortality rate for all lengths,

\[ M(L) = M. \]

Length-dependent options allow smaller fish to experience higher natural mortality. These options can be useful when the biology of the species or the assessment context suggests that mortality should decrease with size.

7 Model formulation

7.1 Gear-specific fishing mortality and selectivity

FLicc represents gear-specific fishing mortality as the product of an apical fishing mortality and a selectivity curve:

\[ F_g(L) = F_g S_g(L), \]

where \(g\) indexes gear, \(F_g\) is the gear-specific apical fishing mortality, and \(S_g(L)\) is selectivity-at-length scaled to a maximum of one.

Total fishing mortality-at-length is

\[ F(L) = \sum_g F_g S_g(L), \]

and total mortality-at-length is

\[ Z(L) = M(L) + F(L). \]

7.2 Selectivity functions

The selectivity function is specified separately for each gear with the sel_fun argument in fiticc().

7.2.1 Logistic selectivity

Logistic selectivity is appropriate for gears that approach full selection above a certain length:

\[ S(L) = \frac{1}{1 + \exp[-(L - L_{50}^{sel}) / \delta]}. \]

Here \(L_{50}^{sel}\) is the length at 50% selectivity and \(\delta\) controls the steepness of the curve.

7.2.2 Double-sided normal selectivity

Double-sided normal selectivity is useful for dome-shaped gears:

\[ S(L) = \begin{cases} \exp\left[-\frac{(L - \mu)^2}{2\sigma_L^2}\right], & L < \mu, \\ \exp\left[-\frac{(L - \mu)^2}{2\sigma_R^2}\right], & L \ge \mu, \end{cases} \]

where \(\mu\) is the modal selected length and \(\sigma_L\) and \(\sigma_R\) control the ascending and descending limbs.

7.3 Spawning-potential ratio

SPR is the central depletion and status quantity in this vignette. It is calculated as spawning output per recruit under fishing relative to the unfished spawning output per recruit:

\[ SPR = \frac{\sum_L N_F(L) W(L) P_{mat}(L)}{\sum_L N_0(L) W(L) P_{mat}(L)}. \]

Here \(N_F(L)\) is equilibrium numbers-at-length under fishing, and \(N_0(L)\) is numbers-at-length in the unfished population. Lower values of SPR indicate stronger fishing impact on reproductive potential.

8 Fitting the model

8.1 Basic multi-gear fit

The main fitting function is fiticc(). In this example there are two gears. The first is represented by double-sided normal selectivity and the second by logistic selectivity.

set.seed(123)

system.time({
  fit <- fiticc(
    lfd,
    stklen,
    sel_fun = c("dsnormal", "logistic"),
    catch_by_gear = c(0.7, 0.3),
    settings = list(
      prior_sigmaF = c(log(0.5), 0.3, 1),
      CVL = 0.1
    ),
    n_restart = 10
  )
})
#>    user  system elapsed 
#>    0.41    0.00    0.71

The main arguments are:

  • lfd: observed length-frequency data as FLQuants;
  • stklen: biological and length-structured stock object from stocklen();
  • sel_fun: one selectivity function per gear;
  • catch_by_gear: relative catch or catch contribution by gear;
  • settings: model and prior settings passed to the TMB model;
  • n_restart: number of optimisation restarts used to improve convergence.

The gear catch weights should sum to one or be interpretable as relative gear contributions. They determine how gear-specific selectivity and fishing mortality are combined into the fishery-level population effect.

8.2 Parameter and setting options

The settings list controls several important assumptions. Common settings include:

Setting Role Typical use
pop_model Population-at-length formulation "gamma" or "gtg"
obs_model Observation likelihood for length data "mn", "nb", or "dm"
CVL Coefficient of variation in length-at-age or growth-transition spread Fixed input or starting value
CVL.sd Prior standard deviation for CVL on transformed scale Penalise weakly or strongly
linf.sd Prior standard deviation for \(L_\infty\) on transformed scale Stabilise growth input
Mk.sd Prior standard deviation for \(M/k\) on transformed scale Stabilise mortality-growth ratio
prior_sigmaF Prior on interannual variation in fishing mortality Used in multi-year fits

For a basic single-fit example, CVL and the selectivity structure are often the most influential settings. For fishblicc-like applications, the gamma population model and negative-binomial observation model are commonly used. For LBSPR-like applications, the gtg population model and multinomial observation model provide a faster approximation.

9 Fit diagnostics and output

9.1 Log-likelihood and information criteria

The fitted likelihood can be extracted and used to compute information criteria.

ll <- LLflicc(fit)
ll[[1]]
#> [1] -35786.55
AIC(ll)
#> [1] 71633.1
BIC(ll)
#> [1] 71790.25

These values are most useful for comparing alternative models fitted to the same data, for example different selectivity assumptions or observation likelihoods.

9.2 Convergence diagnostics

Always inspect convergence before interpreting biological results.

flicc_convergence(fit)
#>                              check              value passed
#> 2           Gradient is acceptable           0.008151   TRUE
#> 3     Hessian is positive definite               TRUE   TRUE
#> 4        Standard errors available 30 / 30 parameters   TRUE
#> 5 Hessian eigenvalues are positive              8.606   TRUE
#>                                            interpretation
#> 2                 The fit is close to a stationary point.
#> 3                    The model curvature is well behaved.
#> 4      Uncertainty could be estimated for all parameters.
#> 5 No flat or negative curvature directions were detected.
fit$restart_log
#>   restart objective max_gradient improved
#> 1       0  35786.55  0.484589999     TRUE
#> 2       1  35786.55  0.095093100     TRUE
#> 3       2  35786.55  0.052737235     TRUE
#> 4       3  35786.55  0.008151018     TRUE
#> 5       4  35786.55  0.009442039    FALSE

The convergence summary checks optimiser status, maximum gradient, Hessian behaviour and standard-error availability. Large gradients or non-positive-definite Hessians indicate that the model may not be locally well identified or that the selected model structure is not well supported by the data.

9.3 Report object and derived quantities

The TMB report contains estimated and derived model quantities.

summary(fit$report)
#>               Length Class    Mode     
#> N             696    FLQuant  numeric  
#> w_at_len      696    FLQuant  numeric  
#> sel_gear        2    FLQuants list     
#> sel           696    FLQuant  numeric  
#> catch_by_gear   2    FLQuants list     
#> predlen_gear    2    FLQuants list     
#> predlen       696    FLQuant  numeric  
#> obslen          2    FLQuants list     
#> Mk            696    FLQuant  numeric  
#> Fk              2    FLQuants list     
#> F             696    FLQuant  numeric  
#> Fap            12    FLQuant  numeric  
#> M             696    FLQuant  numeric  
#> mat           696    FLQuant  numeric  
#> spr            12    FLQuant  numeric  
#> lhpar           9    FLPar    numeric  
#> selpars         2    FLPars   list     
#> pars            4    FLPar    numeric  
#> Lmid           58    -none-   numeric  
#> logLik          1    -none-   numeric  
#> pop_model       1    -none-   character
#> obs_model       1    -none-   character
#> FLReport        1    -none-   logical
fcur_flicc(fit)
#> [1] 0.1282293
sprcur_flicc(fit)
#> [1] 0.252691

fcur_flicc() extracts current fishing mortality, while sprcur_flicc() extracts current SPR.

10 Visual diagnostics

10.1 SPR trajectory or fitted SPR summary

plot_spr(fit)

10.2 Observed and predicted length compositions

plot_len(fit)

plot_len(fit, by_gear = TRUE, year = 2020:2024)

The overall plot is useful for checking whether the model reproduces the aggregate size structure. The by-gear plot is more diagnostic because it shows whether individual gear selectivity assumptions are plausible.

10.3 Selectivity

plot_sel(fit)

selpars_flicc(fit)
#> An object of class "FLPars"
#> [[1]]
#> An object of class "FLPar"
#> params
#>  mode   lsd   rsd 
#> 17.34  1.90  7.04 
#> units:  cm 
#> 
#> [[2]]
#> An object of class "FLPar"
#> params
#> SL50 SL95 
#> 26.5 31.3 
#> units:  cm 
#> 
#> Slot "names":
#> [1] "Trawl"   "Gillnet"
#> 
#> Slot "desc":
#> character(0)
#> 
#> Slot "lock":
#> [1] FALSE
do.call(c, lapply(selpars_flicc(fit), function(x) x))
#>    Trawl1    Trawl2    Trawl3  Gillnet1  Gillnet2 
#> 17.343277  1.902751  7.041593 26.487143 31.320355

The estimated selectivity curves should be checked against gear knowledge. Dome-shaped selectivity can be statistically attractive but biologically implausible if used without supporting evidence. Conversely, multiple gears that sample different length ranges can make dome-shaped and logistic patterns more identifiable than a single aggregated length composition.

11 Creating fitted FLStockLen output

The fitted model can be converted back to an FLStockLen object.

stkl <- flicc_stklen(fit)
plot_sel(stkl)

This is useful for downstream FLR workflows and for visualising fishery-weighted selectivity under the fitted model.

12 Equilibrium dynamics and reference points

12.1 Per-recruit equilibrium curves

Equilibrium curves are produced with eqstklen(). The steepness argument s is used when translating per-recruit quantities into MSY-proxy reference points.

eqstk <- eqstklen(fit, s = 0.75)
eqstk@refpts
#> An object of class "FLPar"
#> params
#>   Fcur   Fspr   Ycur SSBcur SPRcur   SPR0 
#>  0.128  0.081  0.937  0.185  0.253  3.323 
#> units:  NA
plot_eqcurves(eqstk)

The equilibrium curves show how SPR, relative spawning biomass and relative yield change with fishing mortality. They help interpret the fitted fishing mortality in relation to common SPR proxies.

12.2 Per-recruit biological reference points

prbrp <- prbrp_flicc(fit)

fmsy_flicc(prbrp, s = 0.75)
#> [1] 0.0972
fspr_flicc(fit, spr = 40)
#> [1] 0.08098899
fspr_flicc(fit, spr = 35)
#> [1] 0.09423604
round(spr_flicc(fit, F = fspr_flicc(fit, spr = 40)), 3)
#> [1] 0.4

fspr_flicc() finds the fishing mortality associated with a target SPR level. For example, spr = 40 corresponds to \(F_{SPR40\%}\). These quantities are often more defensible in data-limited settings than fully estimated MSY reference points.

13 Length-based indicators and advice plots

13.1 Length-based indicator concept

Length-based indicators in FLicc compare the observed proportion of large fish with the proportion expected under an equilibrium reference condition. The main idea is to translate a fitted SPR reference point into a length-structure benchmark and then evaluate whether the observed length composition contains more or fewer large fish than expected under that benchmark.

This follows the same logic as many age- or size-structure indicators: a threshold is first defined, then the observed proportion above that threshold is compared with a reference proportion. In FLicc, the threshold is commonly expressed as a fraction of a reference length, for example thresh = 0.75. Conceptually, this defines a large-fish threshold

\[ L^* = q L_{ref}, \]

where \(q\) is the user-specified threshold fraction and \(L_{ref}\) is the relevant reference length scale used internally by the fitted length model. For a given year \(y\) and gear \(g\), the observed proportion of fish above the threshold is

\[ P^{obs}_{g,y}(L \geq L^*) = \frac{\sum_{i: L_i \geq L^*} y_{g,y,i}}{\sum_i y_{g,y,i}}, \]

where \(y_{g,y,i}\) is the observed count or effective count in length bin \(i\). The corresponding expected proportion under the fitted model is

\[ P^{fit}_{g,y}(L \geq L^*) = \sum_{i: L_i \geq L^*} \hat{\pi}_{g,y,i}, \]

where \(\hat{\pi}_{g,y,i}\) is the fitted expected length-composition probability. A simple proportional indicator can then be written as

\[ LBI_{g,y} = \frac{P^{obs}_{g,y}(L \geq L^*)}{P^{ref}_{g,y}(L \geq L^*)}, \]

where \(P^{ref}\) may be defined from the fitted or SPR-reference length composition, depending on the diagnostic being used. Values below 1 indicate that the observed sample contains fewer large fish than expected under the reference structure, while values above 1 indicate relatively more large fish.

This indicator is useful because it converts the full fitted length model into an interpretable diagnostic of size structure. It is not a replacement for the fitted SPR estimate, but it provides a complementary view of whether the observed length composition is consistent with the expected proportion of larger individuals.

LBIidx <- LBIspr(fit, thresh = 0.75)
LBIidx
#> An object of class "FLIndices"
#> [[1]]
#> An object of class "FLIndexBiomass"
#> 
#> Name:  
#> Description:  
#> Distribution :   
#> Quant: age 
#> Dims:  age   year    unit    season  area    iter
#>  1   12  1   1   1   1   
#> 
#> Range:  min  max pgroup  minyear maxyear startf  endf 
#>  NA  NA  NA  2013    2024    0.4 0.6 
#> 
#> 
#> [[2]]
#> An object of class "FLIndexBiomass"
#> 
#> Name:  
#> Description:  
#> Distribution :   
#> Quant: age 
#> Dims:  age   year    unit    season  area    iter
#>  1   12  1   1   1   1   
#> 
#> Range:  min  max pgroup  minyear maxyear startf  endf 
#>  NA  NA  NA  2013    2024    0.4 0.6 
#> 
#> 
#> Slot "names":
#> [1] "Trawl"   "Gillnet"
#> 
#> Slot "desc":
#> character(0)
#> 
#> Slot "lock":
#> [1] FALSE

plot_LBIspr() summarises the indicator through time. This is useful when the fitted model spans multiple years and the objective is to track whether the large-fish component is increasing, stable or declining.

plot_LBIspr(fit, thresh = 0.75)

plot_LBIp() provides a more direct visualisation of the proportions used by the indicator. For example, the following plot illustrates the observed and reference proportions for the most recent year in the example dataset.

plot_LBIp(fit, year = 2024)

In practice, plot_LBIp() is often the easiest diagnostic to explain to non-technical users because it shows the proportion of fish above the length threshold directly. The time-series indicator from plot_LBIspr() is more useful for comparing years or fleets.

13.2 Status and advice plots

plot_spr(fit)

plot_lbfao(fit)

A simplified FLStockR object can be created for length-based advice plots.

stk <- flicc2FLStockR(fit)
stk@refpts
#> An object of class "FLPar"
#> params
#>  Fspr  Bspr   Bpa  Blim 
#> 0.081 0.400 0.200 0.100 
#> units:  NA
plot_LBAdvice(stk)

The same object can be expressed relative to an MSY proxy.

stkr <- flicc2FLStockR(fit, rel = TRUE)
plot_LBAdvice(stkr, panel = 1) + ylim(0, 1.5)

14 Fitting each year separately

FLicc can also fit each year independently using by_year = TRUE. This is closer to an LBSPR-style workflow, where each annual length composition produces its own estimate.

fit.y <- fiticc(
  lfd,
  stklen,
  sel_fun = c("dsnormal", "logistic"),
  catch_by_gear = c(0.7, 0.3),
  settings = list(
    prior_sigmaF = c(log(0.5), 0.3, 1)
  ),
  by_year = TRUE
)
stky <- flicc_stklen(fit.y)
plot_sel(stky)

plot_spr(list(all.yr = fit, each.y = fit.y))

The joint multi-year fit shares information across years and can stabilise selectivity and fishing mortality estimates. The by-year fit is more flexible but can be noisier, especially when sample sizes are low or gear composition changes among years.

15 Handling missing gear-year length data

Real multi-gear length-frequency datasets are often unbalanced. Some gears may not have been sampled in all years, or a year may be missing for one gear but available for another. FLicc permits these missing gear-year combinations to be represented explicitly as NA values in the input FLQuants object. This is preferable to replacing missing samples with zeros, because zeros can be true observations within an observed length composition, whereas a missing gear-year sample means that no likelihood contribution should be made for that gear and year.

The example below creates a copy of the length-frequency data and removes selected gear-year combinations by setting them to NA. The model is then fitted in the same way as before.

lfd_na <- lfd
lfd_na$Gillnet[, ac(c(2020, 2022))][] <- NA
lfd_na$Trawl[, ac(c(2020:2021))][] <- NA

fit.na <- fiticc(
  lfd_na,
  stklen,
  sel_fun = c("dsnormal", "logistic"),
  catch_by_gear = c(0.7, 0.3),
  settings = list(
    prior_sigmaF = c(log(0.5), 0.3, 1),
    CVL = 0.1
  )
)

The by-gear length-composition plot shows the observed bars only where data are available. Missing gear-year combinations have no grey observed bars, but the fitted blue prediction is still shown because the model can still predict the expected length composition from the estimated population and selectivity model.

plot_len(fit.na, by_gear = TRUE, year = 2019:2024)

The SPR trajectory from the fit with missing gear-year samples can be compared with the original fit. Similar trajectories indicate that the missing samples are not being treated as zero catches or zero counts. Instead, the available length-composition data continue to inform the multi-year model.

plot_spr(list(fit = fit, NAs = fit.na))

This behaviour is useful for applied assessments because it allows an explicitly rectangular year-by-gear input object while preserving the distinction between true zero counts and missing length samples.

16 Population and observation model options

FLicc separates two modelling choices that are sometimes bundled together in length-based assessment software:

  1. the population model, which determines the expected numbers- or proportions-at-length under growth, natural mortality, fishing mortality and selectivity; and
  2. the observation model, which links the observed length-frequency data to the expected length composition.

This separation is useful because the same biological and selectivity assumptions can be fitted with different likelihoods, and the same observation model can be used with different approximations to the equilibrium population-at-length calculation.

16.1 Population model: expected length composition

For a given length bin \(i\) and gear \(g\), the fitted length composition is driven by the expected catch-at-length. In simplified notation,

\[ C_{ig} = \frac{F_{ig}}{Z_i}\left(N_i - N_{i+1}\right), \]

where \(N_i\) is the expected number of fish entering length interval \(i\), \(F_{ig}\) is fishing mortality for gear \(g\) in that interval, and

\[ Z_i = M_i + \sum_g F_{ig} \]

is total mortality-at-length. Expected proportions for the length-composition likelihood are then obtained by normalising expected catch-at-length:

\[ \pi_{ig} = \frac{C_{ig}}{\sum_i C_{ig}}. \]

FLicc currently supports two population model options through settings$pop_model.

16.1.1 GTG population model (pop_model = "gtg")

The growth-type-group, or GTG, formulation is the default fast option. It approximates growth variability by representing the population as a mixture of growth trajectories or length-transition groups. Within each group, survival through length intervals can be written recursively as

\[ N_{i+1,h} = N_{i,h}\exp(-Z_{i,h}\Delta t_{i,h}), \]

where \(h\) indexes the growth group and \(\Delta t_{i,h}\) is the expected time spent in length interval \(i\). The modelled population-at-length is then obtained by summing over growth groups,

\[ N_i = \sum_h \omega_h N_{i,h}, \]

with growth-group weights \(\omega_h\).

The GTG formulation is close in spirit to LBSPR-style length-based per-recruit models. It is generally fast, stable and recommended as the default option for routine FLicc analyses, especially with a multinomial observation model.

settings = list(pop_model = "gtg", obs_model = "mn")

16.1.2 Gamma population model (pop_model = "gamma")

The gamma formulation follows the original fishblicc-style calculation more closely by integrating over variability in asymptotic length. In simplified form,

\[ N_i = \int f(L_\infty) S_i(L_\infty)\,dL_\infty, \]

where \(f(L_\infty)\) describes variation in asymptotic length and \(S_i(L_\infty)\) is survival to length interval \(i\) conditional on a given growth trajectory. In practice this integral is approximated numerically.

This approach is more directly aligned with the fishblicc formulation and may be useful when reproducing fishblicc examples or when sensitivity to the population-at-length approximation is being explored. It is typically slower than the GTG option.

settings = list(pop_model = "gamma", obs_model = "nb")

16.2 Observation model: likelihood for length data

Let \(y_{ig}\) be the observed count in length bin \(i\) for gear \(g\), \(n_g = \sum_i y_{ig}\) the total effective sample size for that gear, and \(\pi_{ig}\) the expected model proportion.

16.2.1 Multinomial observation model (obs_model = "mn")

The multinomial likelihood is the default and is consistent with many length-composition applications:

\[ \mathbf{y}_g \sim \text{Multinomial}(n_g, \boldsymbol{\pi}_g). \]

The negative log-likelihood contribution is proportional to

\[ -\log L_g = -\sum_i y_{ig}\log(\pi_{ig}), \]

ignoring constants. This option is fast and usually a good first choice when effective sample sizes have been chosen carefully.

16.2.2 Negative-binomial observation model (obs_model = "nb")

The negative-binomial likelihood treats length-bin counts as overdispersed counts around expected values,

\[ y_{ig} \sim \text{NB}(\mu_{ig}, \phi), \qquad \mu_{ig} = n_g\pi_{ig}, \]

where \(\phi\) is an overdispersion parameter. One common variance parameterisation is

\[ \text{Var}(y_{ig}) = \mu_{ig} + \frac{\mu_{ig}^2}{\phi}. \]

This option is useful for fishblicc-style comparisons and for data where nominal sample sizes may overstate the true information content. It can be slower and, in some applications, less numerically stable than the multinomial option.

16.2.3 Dirichlet-multinomial observation model (obs_model = "dm")

The Dirichlet-multinomial likelihood allows overdispersion in the composition vector rather than modelling bins independently:

\[ \mathbf{y}_g \sim \text{DirichletMultinomial}(n_g, \boldsymbol{\pi}_g, \theta), \]

where \(\theta\) controls extra-multinomial variation. Conceptually, the expected proportions are treated as uncertain composition probabilities before counts are drawn. This can be useful when length-composition data are clustered, correlated or otherwise more variable than expected under a multinomial likelihood.

16.3 Effective sample size

Length-composition likelihoods are sensitive to the assumed effective sample size. FLicc allows sample sizes to be adjusted by gear using lfdess().

lfd.ess <- lfdess(lfd, ess.g = c(Trawl = 250, Gillnet = 350))
plot_lfd(lfd.ess)

This does not change the observed proportions, but changes the weight assigned to each gear’s length composition in the likelihood. This is often more appropriate than using raw sample sizes when fish are measured in clustered samples, trips, hauls or markets.

16.4 Practical combinations

The two most useful starting configurations are:

Purpose pop_model obs_model Comment
Default FLicc / LBSPR-style workflow "gtg" "mn" Fast, stable, good first model
fishblicc-style replication "gamma" "nb" Slower, closer to original fishblicc setup
Overdispersed compositions "gtg" or "gamma" "dm" Useful sensitivity run

16.5 Comparing fishblicc-like and LBSPR-like settings

The following code compares a fishblicc-like gamma + nb model with a faster gtg + mn model after setting gear-specific effective sample sizes.

system.time({
  fit.gamma.nb <- fiticc(
    lfd.ess,
    stklen,
    sel_fun = c("dsnormal", "logistic"),
    catch_by_gear = c(0.7, 0.3),
    settings = list(pop_model = "gamma", obs_model = "nb")
  )
})
#>    user  system elapsed 
#>    4.86    0.97    8.64

system.time({
  fit.gtg.mn <- fiticc(
    lfd.ess,
    stklen,
    sel_fun = c("dsnormal", "logistic"),
    catch_by_gear = c(0.7, 0.3),
    settings = list(pop_model = "gtg", obs_model = "mn")
  )
})
#>    user  system elapsed 
#>    0.45    0.00    0.66

plot_spr(list(
  gamma.nb = fit.gamma.nb,
  gtg.mn   = fit.gtg.mn
))

Differences between the two fits indicate sensitivity to the population-at-length approximation and observation-error model. In a routine analysis, such comparisons are useful diagnostics rather than automatic model-selection rules.

17 Reproducing the fishblicc Bombay duck example

FLicc includes a fishblicc-style Bombay duck example with three gears. This example is useful for checking consistency with fishblicc assumptions.

data("fishblicc_example")

lhpar_fishblicc <- rbind(
  lhpar_fishblicc,
  FLPar(M = lhpar_fishblicc["k"] * lhpar_fishblicc["Mk"])
)

stklen_fishblicc <- stocklen(
  lfd_fishblicc,
  lhpar_fishblicc,
  m_model = "inverse",
  reflen = 21.9
)

The fit below uses penalties on \(L_\infty\), \(M/k\) and CVL to approximate the original fishblicc setup.

fit_fishblicc <- fiticc(
  lfd_fishblicc,
  stklen_fishblicc,
  sel_fun = c("dsnormal", "dsnormal", "dsnormal"),
  catch_by_gear = c(0.1802070, 0.2101353, 0.6096577),
  settings = list(
    pop_model = "gamma",
    obs_model = "nb",
    CVL = 0.14,
    CVL.sd = 0.1,
    linf.sd = 2 / 40,
    Mk.sd = 0.1
  )
)
flicc_convergence(fit_fishblicc)
#>                              check              value passed
#> 2           Gradient is acceptable           0.002302   TRUE
#> 3     Hessian is positive definite               TRUE   TRUE
#> 4        Standard errors available 16 / 16 parameters   TRUE
#> 5 Hessian eigenvalues are positive              6.945   TRUE
#>                                            interpretation
#> 2                 The fit is close to a stationary point.
#> 3                    The model curvature is well behaved.
#> 4      Uncertainty could be estimated for all parameters.
#> 5 No flat or negative curvature directions were detected.
sprcur_flicc(fit_fishblicc)
#> [1] 0.3408298
plot_sel(fit_fishblicc)

plot_len(fit_fishblicc, by_gear = TRUE)

plot_len(fit_fishblicc)

eqstk_fishblicc <- eqstklen(fit_fishblicc, s = 0.7)
eqstk_fishblicc@refpts
#> An object of class "FLPar"
#> params
#>   Fcur   Fspr   Ycur SSBcur SPRcur   SPR0 
#>  1.098  0.925  0.945  0.262  0.341  1.517 
#> units:  NA
plot_eqcurves(eqstk_fishblicc)

18 Practical guidance

18.1 Choosing selectivity functions

Use gear knowledge wherever possible. Logistic selectivity is appropriate when larger fish remain fully selected. Double-sided normal selectivity is appropriate when the gear or fishing process is expected to be dome-shaped. Avoid adding dome-shaped selectivity only to improve fit unless the pattern is supported by the fishery.

18.2 Choosing natural mortality-at-length

Start with m_model = "constant" for transparent baseline analyses. Explore length-dependent alternatives when there is biological justification or when sensitivity to juvenile mortality is important. Report the choice clearly because \(M(L)\) strongly affects the inferred survival curve and SPR.

18.3 Choosing observation models

Use obs_model = "mn" for a simple multinomial length-composition likelihood, "nb" where overdispersion in counts is expected, and "dm" where Dirichlet-multinomial overdispersion is preferred. The negative-binomial and Dirichlet-multinomial options are often more robust when sample sizes are nominal rather than truly effective.

18.4 Priors and penalties

Priors or penalties can stabilise weakly identified biological parameters. In FLicc these are commonly provided through standard-deviation settings such as linf.sd, Mk.sd and CVL.sd. Smaller standard deviations imply stronger penalties around the input values. Strong priors should only be used when the life-history information is reliable.

18.5 Interpreting SPR

SPR is a relative reproductive potential metric, not a direct estimate of absolute biomass. In length-only or length-dominated assessments, SPR is often more robust than absolute biomass because it depends on the equilibrium length structure implied by growth, mortality, maturity and selectivity. However, SPR remains sensitive to life-history assumptions and selectivity.

19 Summary

This vignette demonstrated a complete basic FLicc workflow:

The next step after this basic workflow is to propagate uncertainty using Monte-Carlo ensembles of life-history and model assumptions. That ensemble workflow is intentionally left for a separate vignette.