knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 6 )
This vignette simulates the original single data set from the MEE paper, fits it using the eight different methods and produces the resulting figures. Fitting and plotting functions can be used on a user's own data set.
library(sizeSpectra)
Generate random sample from bounded power-law distribution, where the resulting
x
is a vector of individual fish sizes (usually body masses).
n = 1000 # sample size b.known = -2 # known fixed value of b xmin.known = 1 # known fixed value of xmin xmax.known = 1000 # known fixed value of xmax set.seed(42) # To get the same observations for each run of code. x = rPLB(n, b = b.known, xmin = xmin.known, xmax = xmax.known) head(x)
Make the standard histogram with a break in the y-axis (by first fitting using the Llin method since that returns the required bin breaks, counts etc.).
num.bins = 8 # Suggested number of bins for standard histogram and Llin # method; Daan et al. used 8 bins. # hAAA is the h(istrogram) results for method AAA. # Llin method - plotting binned data on log-linear axes then fitting regression. hLlin = Llin.method(x, num.bins = num.bins) gap.barplot.cust(hLlin$counts, midpoints = hLlin$mids, breakpoints = hLlin$breaks, xlim = c(-10,max(hLlin$breaks)+10), col = rep("grey", length(hLlin$counts)), xaxs = "i", yaxs = "i" )
Calculate size spectra fits using all eight methods on the same data set:
eight.results = eightMethodsMEE(x, num.bins = num.bins)
Now plot Figure 2, where each panel corresponds to one of the methods. Most of the function code involves tailoring the plots to make good figures, and may need tweaking for another dataset.
eight.methods.plot(eight.results)
Recommended plots of (a) biomass size spectrum and (b) abundance size spectrum, fitted using MLE:
MLE.plots.recommend(x = x, b.MLE = eight.results$hMLE.list$b, confVals.MLE = eight.results$hMLE.list$confVals, hLBNbiom.list = eight.results$hLBNbiom.list)
Same as (h) above (Figure 2(h) in paper), but with a non-logarithmic y-axis, to help explain why the red curve on log-log axes in (h) does not pass through the maximum data point.
par(mai=c(0.8, 0.8, 0.2, 0.3)) MLE.plot(x, b = eight.results$hMLE.list$b, log="x")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.