Wrangling the data
We first wrangle the data into the format required for using the MLEbins method. Similar wrangling would need to be done for users’ data to get the required inputs for the likelihood calculations.
For demonstration purposes we will fit the size spectrum for just one species group (Cephalopoda) and for one strata, the fishing grounds (fg). So we first filter to obtain just those values:
dat <- dplyr::filter(mediterranean_data,
strata == "fg",
group == "Cephalopoda")
dat
# A tibble: 157 × 5
strata group species length number
<fct> <fct> <fct> <dbl> <dbl>
1 fg Cephalopoda Abralia veranyi 42 58.8
2 fg Cephalopoda Todarodes sagittatus 294 7.35
3 fg Cephalopoda Todarodes sagittatus 185 7.35
4 fg Cephalopoda Todarodes sagittatus 269 7.35
5 fg Cephalopoda Todaropsis eblanae 116 7.35
6 fg Cephalopoda Todaropsis eblanae 126 7.35
7 fg Cephalopoda Todaropsis eblanae 87 7.35
8 fg Cephalopoda Todaropsis eblanae 130 7.35
9 fg Cephalopoda Todaropsis eblanae 94 7.35
10 fg Cephalopoda Todaropsis eblanae 120 7.35
# ℹ 147 more rows
summary(dat)
strata group species length
fg :157 Actinopterygii: 0 Illex coindetii :53 Min. : 17.00
ntr : 0 Cephalopoda :157 Todaropsis eblanae :38 1st Qu.: 46.00
baseline: 0 Chondrichthyes: 0 Abralia veranyi :37 Median : 95.00
Crustacea : 0 Rossia macrosoma :11 Mean : 90.62
Todarodes sagittatus: 9 3rd Qu.:115.00
Sepia orbignyana : 7 Max. :335.00
(Other) : 2
number
Min. : 7.151
1st Qu.: 7.204
Median : 7.353
Mean :12.505
3rd Qu.: 7.435
Max. :61.345
Thus, dat contains data for just one strata
and one group, but multiple species.
The lengths are in mm, with integer values measured to nearest mm using standard lab metric tape. We convert these length values to minima and maxima of bins, with the measured lengths representing the midpoints. This keeps track of the uncertainty in the ‘true’ length of individuals and is needed for the MLEbins method, which also allows for other species groups having different resolution.
This conversion is done using calc_bin_breaks(), which
takes a tibble that contains a column of measurements and appends
minimum and maximum values for each bin (see
?calc_bin_breaks):
dat_with_breaks <- calc_bin_breaks(dat,
bin_width = 1) %>%
dplyr::rename(bin_count = number)
dat_with_breaks
# A tibble: 157 × 7
strata group species length bin_count length_bin_min length_bin_max
<fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
1 fg Cephalopoda Abralia veranyi 42 58.8 41.5 42.5
2 fg Cephalopoda Todarodes sagi… 294 7.35 294. 294.
3 fg Cephalopoda Todarodes sagi… 185 7.35 184. 186.
4 fg Cephalopoda Todarodes sagi… 269 7.35 268. 270.
5 fg Cephalopoda Todaropsis ebl… 116 7.35 116. 116.
6 fg Cephalopoda Todaropsis ebl… 126 7.35 126. 126.
7 fg Cephalopoda Todaropsis ebl… 87 7.35 86.5 87.5
8 fg Cephalopoda Todaropsis ebl… 130 7.35 130. 130.
9 fg Cephalopoda Todaropsis ebl… 94 7.35 93.5 94.5
10 fg Cephalopoda Todaropsis ebl… 120 7.35 120. 120.
# ℹ 147 more rowsWe then convert data of counts of species in length bins to counts in
body-mass bins calculated using the species-specific length-weight
coefficients. This is done using
length_bins_to_body_mass_bins()
dat_joined <-
length_bins_to_body_mass_bins(dat_with_breaks,
mediterranean_length_weight_coefficients,
length_data_unit = "mm")
dat_joined
# A tibble: 157 × 11
strata group species length bin_count length_bin_min length_bin_max alpha beta
<fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fg Cepha… Abrali… 42 58.8 4.15 4.25 0.689 1.42
2 fg Cepha… Todaro… 294 7.35 29.4 29.4 0.0102 3.31
3 fg Cepha… Todaro… 185 7.35 18.4 18.6 0.0102 3.31
4 fg Cepha… Todaro… 269 7.35 26.8 27.0 0.0102 3.31
5 fg Cepha… Todaro… 116 7.35 11.6 11.6 0.0996 2.79
6 fg Cepha… Todaro… 126 7.35 12.6 12.6 0.0996 2.79
7 fg Cepha… Todaro… 87 7.35 8.65 8.75 0.0996 2.79
8 fg Cepha… Todaro… 130 7.35 13.0 13.0 0.0996 2.79
9 fg Cepha… Todaro… 94 7.35 9.35 9.45 0.0996 2.79
10 fg Cepha… Todaro… 120 7.35 12.0 12.0 0.0996 2.79
# ℹ 147 more rows
# ℹ 2 more variables: weight_bin_min <dbl>, weight_bin_max <dbl>We then retain just the columns needed to fit using the MLEbins
method for input into determine_xmin_and_fit_mlebins() and
also sum the counts for the same combinations of species and weight
bins, which is done in mediterranean_for_mlebins(). That
function should work for similarly structured data (requires
group and strata, but the function can be
easily adapted if needed). Since we repeated the same types of analysis
for multiple species groups and strata combinations, it made sense to
write a function. As an extra step here we also are going to look at
only small cephalopods below 20.9 g; see page 8 of Appendix B of [Quevedo
et al. (2026)] for details.
data_cephsmall_fg <- mediterranean_for_mlebins(dat_joined) %>%
dplyr::filter(bin_min < 20.9) %>%
dplyr::arrange(bin_min)
data_cephsmall_fg
# A tibble: 29 × 4
species bin_min bin_max bin_count
<fct> <dbl> <dbl> <dbl>
1 Abralia veranyi 2.46 2.60 61.3
2 Abralia veranyi 2.60 2.75 64.6
3 Abralia veranyi 2.89 3.05 7.35
4 Illex coindetii 3.14 3.35 43.9
5 Abralia veranyi 3.35 3.51 14.8
6 Abralia veranyi 3.51 3.67 64.8
7 Abralia veranyi 3.67 3.83 79.4
8 Illex coindetii 3.79 4.03 43.2
9 Abralia veranyi 3.83 3.99 58.4
10 Abralia veranyi 3.99 4.16 29.4
# ℹ 19 more rowsNote the overlapping bins in rows 7 and 8, and also in rows 8 and 9, because the bins correspond to different species which have different length-weight coefficients. Hence the need for the MLEbins method.
We are now ready to fit the data. The
determine_xmin_and_fit_mlebins() function calculates a
histogram to find the peak from which to start the bounded power-law fit
(Quevedo
et al., 2026). It then fits the size spectrum using
fit_size_spectrum_mlebins(); if you wish to specify
x_min then you would just use
fit_size_spectrum_mlebins(). This can take a few minutes,
so I have saved res_cephsmall_fg as a data object in
sizeSpectraFit (this next chunk is not actually run, don’t tell anyone…;
so do make sure to take out eval = FALSE if you use this
code as a template):
The results are
res_cephsmall_fg
$mlebins_fit
$b_mle
[1] -2.997532
$b_conf
[1] -3.162692 -2.836292
$data
# A tibble: 25 × 7
species bin_min bin_max bin_count count_gte_bin_min low_count high_count
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Illex coindetii 3.14 3.35 43.9 952. 908. 952.
2 Abralia veranyi 3.35 3.51 14.8 908. 893. 908.
3 Abralia veranyi 3.51 3.67 64.8 893. 829. 893.
4 Abralia veranyi 3.67 3.83 79.4 829. 706. 829.
5 Illex coindetii 3.79 4.03 43.2 749. 618. 829.
6 Abralia veranyi 3.83 3.99 58.4 706. 648. 749.
7 Abralia veranyi 3.99 4.16 29.4 648. 618. 691.
8 Abralia veranyi 4.16 4.33 75.9 618. 542. 618.
9 Abralia veranyi 4.33 4.49 7.44 542. 535. 542.
10 Abralia veranyi 4.49 4.67 14.7 535. 520. 535.
# ℹ 15 more rows
$x_min
[1] 3.139746
$x_max
[1] 18.04885
$n
[1] 951.977
$method
[1] "MLEbins"
attr(,"class")
[1] "size_spectrum_mlebins" "list"
$h
$breaks
[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
$mids
[1] 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5
[17] 16.5 17.5
$counts
[1] 0.000 0.000 133.257 333.841 289.859 175.430 58.278 0.000 65.590 7.312
[11] 0.000 0.000 0.000 0.000 7.204 0.000 7.312 7.151
$xname
[1] "Total counts in each bin"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
attr(,"class")
[1] "determine_xmin_and_fit_mlebins" "list"This is a list of two objects (each is a list):
$mlebins_fit contains the details of the maximum likelihood
fitting similar to the other methods, and \(h\) is an object of class
histogram that contains details of the histogram used to
determine \(x_{min}\). (Use
remove_hist(res_cephsmall_fg) to remove the histogram
object and retain the class; this also works on lists of results from
multiple applications of
determine_xmin_and_fit_mlebins()).
The object has the class determine_xmin_and_fit_mlebins
such that our customised
plot.determine_xmin_and_fit_mlebins() automatically gets
used:
The top panel is a histogram of total counts of minimum body sizes,
using 1-g body-size bins, used to determine \(x_{min}\). The left-most red bin is the
modal bin, and \(x_{x_min}\) for this
group-strata combination is the minimum value within that 1-g bin.
Middle and bottom panels (as suggested in Figure 7 of [2]) are the fit of the bounded power-law distribution using the MLEbins method with linear and logarithmic y-axes, respectively. For each body-mass bin, the green horizontal line show the range of body sizes, with its value on the y-axis corresponding to the total number of individuals in bins whose minima are \(\geq\) the bin’s minimum; the green lines help to distinguish each bin when bins are overlapping. The vertical span of each grey rectangle shows the possible range of number of individuals with body masses \(\geq\) body mass of individuals in that bin (horizontal span matches the green line). Red curves are fits of the MLE, with dashed lines showing 95% confidence intervals.
These figure are admittedly a little complex, but retain the full uncertainty arising from binned data and conversions through species-specific length-weight relationships. As can be seen, some of the resulting rectangles are quite large, and simply using a single point for each would suggest higher accuracy than is warranted.
Do see the help files for relevant functions, and feel free to email me (or make an issue) with any questions.