In andrew-edwards/sizeSpectra: Fitting Size Spectra to Ecological Data Using Maximum Likelihood

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 6
)

library(sizeSpectra)

Simulate 10,000 data sets and fit them using eight methods, then plot histograms and confidence intervals of the estimates of the size-spectra exponent $b$. Reproduces MEE Figures 3 and 4 and Table 2. Also for the MLEfix method (Figures A.3 and A.4).

library(sizeSpectra)

Figure 3 -- Histograms of estimated exponent b for simulated data

The simulation of 10,000 data sets and fitting using the eight methods takes a while, and so results are saved within the package in eight.results.default. To re-run from scratch use the code in data-raw/simulate-data.R.

list2env(eight.results.default, envir=.GlobalEnv)    # Extract previously saved results


# Plotting
# brange = range( c( Llin.rep, LT.rep, LTplus1.rep, LBmiz.rep,
#    LBbiom.rep, LBNbiom.rep, MLE.rep))
# Comes out as -3.38, 0.233 for seed=42, num.reps = 10000
xrange = c(-3.5, 0.5)         # range of x-axis for histograms - actually to define bins
xbigticks = seq(-3.5, 0.5, 0.5) # NOW THE DEFAULT

xsmallticks = seq(xrange[1], xrange[2], by=0.1)  # ALSO now default
breakshist = seq(xrange[1], xrange[2], by=3/61)
         # bin ends at 0, 2 is centered. width of 4/81 from solving
         #  the mean of a bin equals 2, ((N+1)w + Nw)/2 = 2, and setting number
         #  of the bin that includes 2, N, =40 (which is how many you get with
         #  0.05, which looked good).

# figheight = 7 # 5.6
# figwidth = 5.7    # 5.7 inches for JAE
#postscript("fitting3rep.eps", height = figheight, width = figwidth,
#           horizontal=FALSE,  paper="special")
par(omi = c(0.12, 0.05, 0.12, 0.0))      # outer margins in inches
par(mfrow=c(4,2)) #7,1))

oldmai = par("mai")
par(mai=c(0.5, 0.5, 0.1, 0.3)) # Affects all figures if don't change agaiin
par(xaxs="i", yaxs="i")    # Have to define here for hist
par(mgp=c(2.0, 0.5, 0))    # puts axes labels closer I think
par(cex = 0.8)             # With no option all text comes out a bit small

vertCol = "red"            # Colour for vertical lines
vertThick = 1              # Thickness for vertical lines

xLim = c(-3.2, -1.5)
xLimLlin = xLim + 2
ylimA = c(0, 5500)
ylimLlin = c(0, 10000)
cexAxis = 0.9      # font size for axes labels to make the y ones fit okay

# Llin.rep has different breakhist so that 0 is a breakpoint. Same width as
#  others, just shifted.
hist(Llin.rep.df$slope, xlim=xLimLlin, breaks=breakshist - breakshist[31],
  xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimLlin)  #  ylim=c(0,1100),
axis(1, at=xbigticks, labels = xbigticks, mgp=c(1.7,0.7,0), cex.axis=cexAxis)  # big ticks
axis(1, at=xsmallticks, labels=rep("",length(xsmallticks)), tcl=-0.2)
axis(2, at=c(0, 5000, 10000),
         labels = c(0, 5000, 10000),
         mgp=c(1.7,0.7,0), cex.axis=cexAxis)  # big ticks labelled
axis(2, at=seq(0, 10000, 1000),
         labels = rep("", 11),
         mgp=c(1.7,0.7,0))  # big ticks unlabelled
abline(v=b.known, col=vertCol, lwd=vertThick)
inset = c(-0.08, -0.08)     # inset distance of legend
legend("topleft", "(a) Llin", bty="n", inset=inset)

# Think these three lines can be deleted
figlabpos = 0.93    # proportion of x and y axis lengths to put (a) in etc.
xlabpos = 0.75       # just play with a number, as now using pos=4 in text
text( xlabpos, figlabpos * 1100, "(a) Llin", pos=4)

hist(LT.rep.df$slope, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(b) LT", bty="n", inset=inset)

hist(LTplus1.rep.df$slope, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(c) LTplus1", bty="n", inset=inset)

hist(LBmiz.rep.df$slope - 1, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(d) LBmiz", bty="n", inset=inset)

hist(LBbiom.rep.df$slope - 2, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(e) LBbiom", bty="n", inset=inset)

hist(LBNbiom.rep.df$slope - 1, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(f) LBNbiom", bty="n", inset=inset)


hist(LCD.rep.df$slope - 1, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(g) LCD", bty="n", inset=inset)

hist(MLE.rep.df$b, xlim=xLim, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = ylimA)  #  ylim=c(0,1100),
histAxes()
legend("topleft", "(h) MLE", bty="n", inset=inset)

# text(-2.5, 0, expression(paste("Estimate of ", italic(b))))
mtext(expression(paste("Estimate of ", italic(b))), side=1, outer=TRUE, line=-1)
# mtext("hello", side=1, outer=TRUE, line=-1)

Figure A.3 for just the MLE and MLEfix methods:

par(omi = c(0.12, 0.05, 0.12, 0.0))      # outer margins in inches
par(mfrow=c(2,1))

par(mai=c(0.5, 0.5, 0.1, 0.3))
par(xaxs="i", yaxs="i")    # Have to define here for hist
par(mgp=c(2.0, 0.5, 0))    # puts axes labels closer I think
par(cex = 0.8)             # With no option all text comes out a bit small

xLimFix = c(-2.5, -1.5)
yLimFix = c(0, 6000)
hist(MLE.rep.df$b, xlim=xLimFix, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = yLimFix)  #  ylim=c(0,1100),
histAxes()
axis(2, at=6000, labels = 6000, mgp=c(1.7,0.7,0), cex.axis=cexAxis)
                                        # big tick label
legend("topleft", "(a) MLE", bty="n", inset=0.5*inset)

hist(MLEfix.rep.df$b, xlim=xLimFix, breaks=breakshist, xlab="", ylab="Frequency",
  main="", axes=FALSE, ylim = yLimFix)  #  ylim=c(0,1100),
histAxes()
axis(2, at=6000, labels = 6000, mgp=c(1.7,0.7,0), cex.axis=cexAxis)
                                        # big tick label
legend("topleft", "(b) MLEfix", bty="n", inset=0.5*inset)

mtext(expression(paste("Estimate of ", italic(b))), side=1, outer=TRUE, line=-1)

Table 2, summarising the results for each method:

|Method | Slope represents | 5\% quantile |Median |Mean |95\% quantile |Percentage below true| |:----- | :--------------- | -----------: |-----: |---: |-----------: |--------------------:| |Llin | - | r qqtab(Llin.rep.df$slope, dig=2, true=b.known, quants=c(0.05, 0.95)) | |LT | b | r qqtab(LT.rep.df$slope, true = b.known, quants=c(0.05, 0.95)) | |LTplus1| b | r qqtab(LTplus1.rep.df$slope, true = b.known, quants=c(0.05, 0.95)) | |LBmiz | b+1 | r qqtab(LBmiz.rep.df$slope - 1, true = b.known, quants=c(0.05, 0.95)) | |LBbiom | b+2 | r qqtab(LBbiom.rep.df$slope - 2, true = b.known, quants=c(0.05, 0.95)) | |LBNbiom| b+1 | r qqtab(LBNbiom.rep.df$slope - 1, true = b.known, quants=c(0.05, 0.95)) | |LCD | b+1 | r qqtab(LCD.rep.df$slope - 1, true = b.known, quants=c(0.05, 0.95)) | |MLE | b | r qqtab(MLE.rep.df$b, true = b.known, quants=c(0.05, 0.95)) |

Figure 4 -- Confidence intervals of estimated exponent b for simulated data

Note that MEPS_reproduce_2.html uses a single function for a similar plot; that could be modified for this plot if desired.

par(omi = c(0.14, 0, 0.1, 0.15))      # outer margins in inches
par(mfrow=c(4,2)) #7,1))

oldmai = par("mai")    #   0.6732 0.5412 0.5412 0.2772  inches I think,
                       #    may be indpt of fig size
par(mai=c(0.3, 0.5, 0.08, 0))  # Affects all four figures if don't change agaiin
par(xaxs="i", yaxs="i")    # Have to define here for hist
par(mgp=c(2.0, 0.5, 0))    # puts axes labels closer I think
par(cex = 0.8)             # With no option all text comes out a bit small
vertThick = 1              # Thickness for vertical lines

# Each of these plots a panel for one method. Define xLim if the default
#  (integer-based calculation) is not suitable
Llin.rep.conf.sort = confPlot(Llin.rep.df[ ,-1], legName="(a) Llin",
    xLim = c(-0.25, 0.25)) #, colourCode=FALSE)

LT.rep.conf.sort = confPlot(LT.rep.df[ ,-1], legName="(b) LT", yLab="", yLabels=FALSE)

LTplus1.rep.conf.sort = confPlot(LTplus1.rep.df[ ,-1], legName="(c) LTplus1")

xLimCom = c(-2.65, -1.5)     # common width of axes for the remainder

# range(LBmiz.rep.conf) =  -1.5894281 -0.5268428. Was -1.6, -0.4
LBmiz.rep.conf.sort = confPlot(LBmiz.rep.df[ ,-1] - 1, legName="(d) LBmiz",
    xLim = xLimCom, yLab="", yLabels=FALSE)

# range(LBbiom.rep.conf) = -0.6106510  0.4370617
LBbiom.rep.conf.sort = confPlot(LBbiom.rep.df[ ,-1] - 2, legName="(e) LBbiom",
      xLim = xLimCom)

# range(LBNbiom.rep.conf) = -1.6106510 -0.5629383
LBNbiom.rep.conf.sort = confPlot(LBNbiom.rep.df[ ,-1] - 1, legName="(f) LBNbiom",
     yLab="", yLabels=FALSE, xLim = xLimCom)

# range(LCD.rep.conf) = -1.1744680 -0.8713301   was c(-1.2, -0.8)
LCD.rep.conf.sort = confPlot(LCD.rep.df[ ,-1] - 1, legName="(g) LCD",
    xLim = xLimCom, vertFirst = TRUE)

# c(-2.2, -1.8) - looks good, but better to have consistent
MLE.rep.conf.sort = confPlot(MLE.rep.df[ ,-1], legName="(h) MLE",
    xLim = xLimCom, yLab="", yLabels=FALSE)
mtext(expression(paste("Estimate of ", italic(b)), sep=""),
      side=1, outer=TRUE, line=-0.2, cex=0.8)

And Figure A.4 for the MLE and MLEfix methods:

# Confidence intervals for MLE and MLEfix methods.
#postscript("fitting3confMLEfix.eps", height = 0.8*figheight,
#           width = 0.8*figwidth,
#           horizontal=FALSE,  paper="special")
par(omi = c(0.12, 0.05, 0.12, 0.0))      # outer margins in inches
par(mfrow=c(2,1)) #7,1))

oldmai = par("mai")    #   0.6732 0.5412 0.5412 0.2772  inches I think,
                       #    may be indpt of fig size
par(mai=c(0.5, 0.5, 0.1, 0.3))  # Affects all four figures if don't change agaiin
par(xaxs="i", yaxs="i")    # Have to define here for hist
par(mgp=c(2.0, 0.5, 0))    # puts axes labels closer I think
par(cex = 0.8)             # With no option all text comes out a bit small

# xrangeMLEfix = c(-2.225, -1.775)  # Range for these two figures
# breakshistMLEfix = seq(xrangeMLEfix[1], xrangeMLEfix[2], by=3/61)

MLE.rep.conf.sort = confPlot(MLE.rep.df[ ,-1], legName="(a) MLE",
    xLim = xLimCom, insetVal = c(-0.04, -0.03), insetVal2 = c(-0.04, 0.05))


MLEfix.rep.conf.sort = confPlot(MLEfix.rep.df[ ,-1], legName="(b) MLEfix",
    xLim = xLimCom, insetVal = c(-0.04, -0.03), insetVal2 = c(-0.04, 0.05))

mtext(expression(paste("Estimate of ", italic(b)), sep=""),
      side=1, outer=TRUE, line=-0.2, cex=0.8)

Range of widths of confidence intervals for MLE method is r round(range(MLE.rep.df$confMax - MLE.rep.df$confMin), dig=3).

Range of widths of confidence intervals for MLEfix method is r range(MLEfix.rep.df$confMax - MLEfix.rep.df$confMin).

Figure 5 for the MLEbin is somewhat superceded by the more general results in the MEPS paper.

andrew-edwards/sizeSpectra documentation built on June 28, 2023, 7:09 p.m.