Measurements of individual body sizes – MLE method
The simplest data consists of just a vector of measurements, such as body masses or body lengths. We fit the size spectrum by calculating the maximum likelihood estimate (the MLE method) of the size-spectrum exponent, along with its 95% confidence intervals.
As an example we use the simulated data from a bounded power-law
distribution as used for Figures 1 and 2 of our MEE
paper. The values are included in sizeSpectraFit as the data object
sim_vec (which stands for simulated vector), and we assume
they are individual body masses measured in grammes. Note that the fit
will be great because the data are sampled from the expected power-law
distribution.
First, verify that there are 1000 values and print the first 20 values to take a look:
length(sim_vec)
[1] 1000
sim_vec[1:20]
[1] 11.613220 15.658847 1.400273 5.869136 2.786321 2.077175 3.785753 1.155444
[9] 2.909813 3.382489 1.842584 3.551050 15.091509 1.342594 1.858151 16.413491
[17] 43.952570 1.132977 1.903030 2.271552The main fitting function is fit_size_spectrum(), which
will automatically use the MLE method here (via
fit_size_spectrum.numeric()) because the data are a
numeric vector of individual values:
res <- fit_size_spectrum(sim_vec)
res
$b_mle
[1] -2.029697
$b_conf
[1] -2.097417 -1.964317
$x
[1] 11.613220 15.658847 1.400273 5.869136 2.786321 2.077175 3.785753 1.155444
[9] 2.909813 3.382489
$x_min
[1] 1.000239
$x_max
[1] 398.7787
$n
[1] 1000
$method
[1] "MLE"The resulting res is a list that contains:
b_mle: the MLE of the exponent \(b\),b_conf: the 95% confidence interval of \(b\),x: the original data (the printing here is set up to only show the first 10 values),x_minandx_max: the MLEs of \(x_{min}\) and \(x_{max}\), which are simply the minimum and maximum of the data,method: tells us that the MLE method was used.
The list res also has class
size_spectrum_numeric which helps with our plotting and
printing functions.
To plot the PLB fit and the original data as an Individual Size Distribution (ISD):
as
recommended in Figure 6b of our MEE
paper. Each individual body mass is shown as a point (note the log
scales), the solid red line is the maximum likelihood estimate (MLE) of
the size-spectrum exponent \(b\), and
the dashed lines are the fits from using the low and high values of the
95% confidence interval of \(b\). The
MLE for \(b\) and the sample size \(n\) are also explicitly given. The default
x-axis label assumes the body masses are in g, but this can easily be
changed.
The plotting is very customisable as described in the help function.
Since the res object has:
the plot(res) command uses our custom function
plot.size_spectrum_numeric(); this is an R feature. Thus,
the help function needed is
?plot.size_spectrum_numeric().
The style argument gives other styles of plotting, for
example:
gives the same plot but with a linear y-axis.
To show both the above plots together, with panels automatically labelled as (a) and (b):
which is similar to Figure 7 of our MEPS
paper but for nonbinned data.
The fourth option is only suitable when the data represent body masses. It shows the normalised biomass:
The black horizontal lines indicate the normalised biomass within each bin on the y-axis, and the span on the x-axis shows the range of the bin. We use bars, rather than the usual points (e.g. Figure 6(f) of our MEE paper) to emphasise that the actual values of body size in each bin span a range; the points somewhat imply that the midpoint of the bin is more relevant.
Successive bins double in size, such that they appear of equal width on the log scale. The normalised biomass within each one is calculated as the sum of the biomass in the bin divided by the bin width; there is no uncertainty in this calculation here because we know the individual body sizes (this is not so if the data are already binned, see later). The horizontal red lines show the estimated normalised biomass in each bin based on the MLE of the exponent \(b\), and the pink boxes span the normalised biomass estimated for all values within the 95% confidence interval of \(b\). The straight lines are also shown for the MLE of \(b\) (solid line) and the 95% confidence interval estimates (dashed lines). Note that the straight lines fit the PLB distribution, whose definition includes the \(x_{min}\) and \(x_{max}\) bounds (see MEE paper equation 1) so the lines are bounded, but the top bin here goes beyond \(x_{max}\), even though there are no actual body masses there (the bin happens to extend out further).
The fit here is clearly very good, but remember this is expected as the data are simulated from a known PLB distribution and are not noisy, unlike real data.
To show the biomass plot and the ISD plot, which is an improved version of the suggested plot in Figure 6 of our MEE paper:
Determine \(x_{min}\) when the data are not simply declining
A simple histogram of your data may not look like a declining
function; in fisheries examples this can be due to selectivity, with the
fishing gear not catching all the smaller organisms. If you have a
similar type of justification for this happening and still want to fit a
declining PLB distribution to a vector of values, we have the built in
function determine_xmin_and_fit() to determine a value of
\(x_{min}\) to use and then fit using
the MLE method. This is the nonbinned version of the approach described
in Quevedo
et al. (2026).
Let’s make up a mock dataset, with some uniformly distributed values between 0.25 and 10, and then a PLB distribution for values \(\geq 8\):
set.seed(42)
mock_data <- c(runif(100, 0.25, 2),
rPLB(1000, -2, x_min = 2, x_max = 100))
hist(mock_data, breaks = 40)
So
even a rudimentary histogram shows a unimodal distribution. To fit
this:
We can plot and then print the results:
The red bars in the histogram show the values that are used to fit the PLB distribution, with \(x_{min}\) being the lowest value in the lowest red bin.
res_mock
$mle_fit
$b_mle
[1] -2.056786
$b_conf
[1] -2.135446 -1.980196
$x
[1] 5.177597 2.540697 2.538831 3.231885 26.180265 35.308150 7.274266 7.106840
[9] 4.210935 2.004465
$x_min
[1] 2.000794
$x_max
[1] 97.04711
$n
[1] 1000
$method
[1] "MLE"
$h
$breaks
[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[27] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
[53] 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
[79] 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
$counts
[1] 38 62 366 172 105 54 43 31 26 25 19 12 12 7 7 8 7 9 7 9
[21] 8 6 4 4 4 3 3 6 2 4 1 2 1 0 2 5 1 2 1 0
[41] 1 1 0 0 0 0 0 0 1 3 3 0 0 2 1 0 0 0 0 0
[61] 1 0 1 0 0 0 0 0 1 0 0 0 1 0 2 0 0 0 0 0
[81] 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 1
$density
[1] 0.0345454545 0.0563636364 0.3327272727 0.1563636364 0.0954545455 0.0490909091
[7] 0.0390909091 0.0281818182 0.0236363636 0.0227272727 0.0172727273 0.0109090909
[13] 0.0109090909 0.0063636364 0.0063636364 0.0072727273 0.0063636364 0.0081818182
[19] 0.0063636364 0.0081818182 0.0072727273 0.0054545455 0.0036363636 0.0036363636
[25] 0.0036363636 0.0027272727 0.0027272727 0.0054545455 0.0018181818 0.0036363636
[31] 0.0009090909 0.0018181818 0.0009090909 0.0000000000 0.0018181818 0.0045454545
[37] 0.0009090909 0.0018181818 0.0009090909 0.0000000000 0.0009090909 0.0009090909
[43] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[49] 0.0009090909 0.0027272727 0.0027272727 0.0000000000 0.0000000000 0.0018181818
[55] 0.0009090909 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[61] 0.0009090909 0.0000000000 0.0009090909 0.0000000000 0.0000000000 0.0000000000
[67] 0.0000000000 0.0000000000 0.0009090909 0.0000000000 0.0000000000 0.0000000000
[73] 0.0009090909 0.0000000000 0.0018181818 0.0000000000 0.0000000000 0.0000000000
[79] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0009090909
[85] 0.0000000000 0.0009090909 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[91] 0.0000000000 0.0000000000 0.0000000000 0.0009090909 0.0000000000 0.0000000000
[97] 0.0000000000 0.0009090909
$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 18.5 19.5 20.5 21.5 22.5 23.5 24.5 25.5 26.5 27.5 28.5 29.5 30.5 31.5
[33] 32.5 33.5 34.5 35.5 36.5 37.5 38.5 39.5 40.5 41.5 42.5 43.5 44.5 45.5 46.5 47.5
[49] 48.5 49.5 50.5 51.5 52.5 53.5 54.5 55.5 56.5 57.5 58.5 59.5 60.5 61.5 62.5 63.5
[65] 64.5 65.5 66.5 67.5 68.5 69.5 70.5 71.5 72.5 73.5 74.5 75.5 76.5 77.5 78.5 79.5
[81] 80.5 81.5 82.5 83.5 84.5 85.5 86.5 87.5 88.5 89.5 90.5 91.5 92.5 93.5 94.5 95.5
[97] 96.5 97.5
$xname
[1] "Body length (x), mm"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
attr(,"class")
[1] "determine_xmin_and_fit" "list"The res_mock object is list of two objects:
$mle_fit: contains the MLE results as in the simpler example above$h: the histogram results as used for calculating \(x_{min}\) and for plotting.
The object is of class determine_xmin_and_fit, and so
the plotting function plot.determine_xmin_and_fit() is
automatically used by plot(res_mock) to give the above
plot.
Note that the true value of \(b = -2\) falls within the confidence interval.
gives
the same plot but with a linear y-axis.
which is similar to, but simpler than, Figure 7 of our 
and
is interpreted the same way as the earlier plot for the MLE method (but
here the bins are already predefined, because that is how the data are
provided).
The y-axes happen to have the same numbers (though they mean different
things) for this simulated data set, though this will not generally be
the case.
