In haddonm/datalowSA: A Package to Faciliate Application of Data poor Stock Assessments
knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE)
options(knitr.kable.NA = "",
knitr.table.format = "pandoc")
options("show.signif.stars"=FALSE,"stringsAsFactors"=FALSE,
"max.print"=50000,"width"=240)
library(datalowSA)
suppressPackageStartupMessages(library(knitr))
Introduction
Which Stock Assessment?
Which stock assessment method to apply to fisheries for data-poor to data-moderate species will depend upon what fisheries and biological data are available but also, importantly, on what management objectives need to be met within the jurisdiction in question. It may be the case that the fishery for a particular species is of sufficient size and value to warrant on-going monitoring and management towards some defined goal for the stock. In such a case the assessment used should obviously be capable of generating some notion of the current state of the fishery and indicate what management actions may be required to eventually achieve the agreed management goals. But some fisheries may be so minor that trying to actively manage them would be inefficient. Nevertheless, to meet the requirements of the Status of key Australian Fish Stocks (SAFS) one still requires some form of defensible stock assessment capable of determining whether the current level of fishing is sustainable.
Age-Structured Production Model
Introduction to ASPM
The age-structured production model (ASPM or aspm) is literally a surplus production model which is based upon an age-structured model of production rather than an accumulated biomass model (see the vignette on spm).
There are some specific data requirements for fitting an age-structured production model to fishery data. Data from the fishery need to included, as a minimum, an accurate catch time-series plus an index of relative abundance for at least some of the years within the catch time-series. In addition, information (or defensible assumptions) is needed for the species concerned in relation to the description of its natural mortality, its growth, its maturation, and the selectivity of the fishery (maturity and selectivity could be knife-edge). If just the catches and CPUE data are available then one might try fitting a simple, aggregated biomass surplus production model. But if these biological data and information are available then an age-structured production model opens the way to on-going improvements with respect to the inclusion of occasional age-composition data or other observations that could be predicted by a suitable model and hence included in the model fitting process.
More details on age-structured production models can be found in Punt et al. (1995). The model equations are provided in the appendix. It would be helpful to the user to read the spm vignette as the theory in there also applies to age-structured production models.
A Typical Workflow
A typical workflow for using an age-structured production model might be something like:
- read in the available data and use checkdata to ensure it can be used with aspm.
- search for suitable initial parameter values using dynamics, aspmLL, and plotASPM, this will include deciding on the use of a two parameter model (no initial depletion) or a three parameter model accounting for an initial depletion level.
- given suitable initial parameters use fitASPM or more basically optim to fit the model to the available data.
- once successfully fitted it is best to plot the solution to determine visually how close the model fit is to the data using plotASPM. One approach to improving this is to include confidence intervals around the index of relative abundance (cpue). This is done first using getLNCI and then including the output of this into plotASPM or, for a closer look in plotceASPM. Eventually improved confidence intervals for the model outputs can be obtained using boostrap samples (see below).
- A better test of the robustness of the solution is to test how sensitive it is to the initial conditions. This can be done by randomly varying the initial parameters and determining whether the model fitting outcomes vary. Suitable example code is given in the vignette.
- After finally deciding on the overall optimum solution it would be sensible to use the optimum parameters to determine the implied production curve so that management statistics such as MSY, the target catch, and the limit and target reference points can be defined. This can be done using getProductionC, the C post-fix denoting this is a C++ function (used for speed), the results of which can be plotted and summarized using prodASPM.
- One can produce a classical phase plot of predicted biomass vs harvest rate to aid in the determination of the stock status, and for this we would use aspmphaseplot.
- Finally, to obtain a better notion of the uncertainty in the analysis, and we can do that using bootASPM, which facilitates the application of a bootstrap analysis.
An Example
The data requirements for the aspm are described above and for the next example we will use the dataspm built in data set. First we will examine an example using only two parameters (assuming the population begins in an unfished state) and then extend the model fitting to include the possibility of an initially depleted state. The two parameters being fitted are the average unfished recruitment level and the standard deviation of the errors around the CPUE data.
We load dataspm, which contains a dataframe of the catches and CPUE by year (fish) and other parameters used for the age-structured production model.
# library(datalowSA)
data(dataspm)
fish <- dataspm$fish
props <- dataspm$props # length-, weight-, maturity- and selectivity-at-age
While a simple, aggregated biomass surplus production model only requires the specification of the species name in the glb data object an aspm requires information on the growth, selectivity, weight-at-age, steepness of the stock-recruit relationship and natural mortality. The global parameters (glb), in addition to the spsname, need to contain the population ages (maxage and ages), natural mortality (M), von Bertalanffy growth parameters (Linf, K and t0), weight-at-age parameters (Waa and Wab), age at 50% maturity, delta (M50a and deltaM), age that 50% of the population are selected by the gear and delta (sela50 and deltaS) and the steepness of the stock-recruitment relationship. The number of years (nyrs) of over which catch and CPUE are available (including missing years) is calculated from the fish dataframe. A starting value for the log of the initial recruitment (R0) also needs to be provided, although this will be estimated along with the standard deviation of the errors around the CPUE data.
We inspect the global parameters specified in dataspm.
(glb <- dataspm$glb)
Just as with the spm model Fitting an aspm model entails first finding initial parameter estimates that lead to predicted cpue time-series that approximate the observed series. Thus, one might begin something like this and uses aspmLL to determine the negative log-likelihood (-verLL), then dynamics to calculate the dynamics of the aspm, and finally plotASPM to illustrate the quality of fit to help decide whether changes are required to the initial guesses.:
pars <- c(12.9,0.25)
aspmLL(pars,infish=fish,inglb=glb,inprops=props)
fishery <- dynamics(pars,infish=fish,inglb=glb,inprops = props)
plotASPM(fishery)
Figure 1. The outcome of a first guess at a two parameter version of the aspm. clearly the fit is very poor and the strong trends in the predicted cpue, the log-normal residuals, and the annual harvest rate bumping up against the built-in upper limit of 0.85 across numerous years, provide a strong indication that the initial guess at unfished average recruitment parameter (R0) is too small. Try increasing it slowly to see its effect on the model fit to the data.
Once reasonable starting values have been found for the parameters ($R_0$ the unfished average recruitment and $\hat\sigma_{I}$, the standard deviation associated with fitting the observed cpue) then an attempt at fitting the model to the data formally can be made using code something like this:
pars <- c(13.7,0.19)
ans <- fitASPM(pars,infish=fish,inglb=glb,inprops=props)
outoptim(ans) # a tidier way of printing the list output from optim
fishery <- dynamics(ans$par,infish=fish,inglb=glb,inprops = props)
Table 1. The output from the fitASPM function and the dynamics function.
kable(fishery,digits=c(0,3,3,3,3,3,4,4,3))
Note the predicted catches are identical to the observed catches. The catches are assumed to be known accurately and so as to ensure a close match between the predicted and observed catches aspmLL has a simple sum-of-squared deviations penalty built into it (try aspmLL, without brackets in the console). This is usually sufficient to force the solution to generate a close match once plausible parameter values are found. Note that with respect to reproduction only the average unfished recruitment is estimated. In its current form the aspm cannot take into accoutn strong and weak cohorts; this remains a very simple model of the dynamics.
To visualize the output of the model fitting we can plot some of the variables from the fishery information generated by the dynamics function.
plotASPM(fishery,CI=NA)
Figure 2. The outcome of fitting a two parameter age-structured production model to the available data. Included is a plot of the catch history, the predicted spawning biomass, the CPUE, the harvest rate, the log-normal residuals for the cpue, and the predicted depletion.
The two-parameter model manages to capture the main cpue trends but fails to capture some of the more obvious and more rapid consistent changes in cpue (Figure 2). The example fishery in question was known to have been fished prior to 1985 so in the year the data begin to be available the stock can be expected to be depleted to some extent. Thus, an alternative might be to fit the model using three parameters. The first approach, with only one parameter of real interest, required data from the beginning of the fishery. However, there are many fisheries for which data are only available after the fishery has been running for a number of years. In such cases it is necessary to estimate the level of depletion and its effect upon recruitment and thus requires two parameters of interest to the dynamics. The first is, as before, the unfished recruitment level, $R_0$, but then we use a parameter that defines the initial depletion at the start of the observations ($D_{init}$. If this parameter is present in pars a search is made for the constant harvest rate that when applied to the initial unfished stock leads to the predicted initial depletion level, only then does the model fitting proceed. Any initial depletion will influence the recruitment depending on the assumed steepness of the stock recruitment relationship, which is assumed to be a Beverton-Holt relationship. The dataspm data set is not particularly suited to a two-parameter model even though that arrangement was able to provide a result; these assessments should not be done automatically, it always takes some background knowledge to ensure that any methods or choices applied are valid.
As an alternative two-parameter example, the deep water fishery data data(fishdat) was a fishery with catch data from the very start of the fishery, which means it is better suited to using a simple two parameter model
data(fishdat)
fish <- fishdat$fish
glb <- fishdat$glb
props <- fishdat$props
pars <- c(14,0.3)
ans <- fitASPM(pars,infish=fish,inglb=glb,inprops = props)
str(ans)
The output contains the estimates of the optimum parameters, the final log-likelihood estimate, the number of iterations it needed to find the optimum, and some diagnostic information. The statement convergence: int 0 implies the solution appears valid, and no warning messages is also encouraging; but see later concerning how to test the robustness of such model fitting. If we put the fitted optimum parameters into the function dynamics we can see the time-series of the more important population variables implied by the model fit.
fishery <- dynamics(ans$par,infish=fish,inglb=glb,inprops = props)
print(fishery)
The use of fitASPM is basically shorthand for using bestL <- optim(pars, aspmLL, method="Nelder-Mead", infish=fish, inglb=glb, inprops=props, control=list(maxit = 1000, parscale = c(10,1))) twice in a row. Examine its code by using fitASPM without brackets in the R console window.
Note that the harvest rate in 1996 appears to have bumped up against the upper limit of 0.85 hard wired into the R-code (try typing dynamics into the R console without brackets to see the code). So that 1180 tonnes catch in 1996 was likely to be damaging. However, the predicted catch is only slightly less than the reported catch so this spike in harvest rate has some support in the data. Plotting up the fishery results enables the visualization of the impact of the fishery and he effect of cutting back catches.
ceCI <- getLNCI(fishery[,"PredCE"],ans$par[2])
plotASPM(fishery,CI=ceCI)
Figure 3. The outcome of fitting a two parameter age-structured production model to the deep-water fishery data. Only the CPUE with confidence intervals are plotted.
The modelling suggests that once catches were reduced to an average of about 246 t between 1997 - 2005 the stock began to very slowly recover, then, after catches were further reduced from 2007 onwards (due to a deepwater closure and cessation of targeting) the rate of recover was predicted to have increased. The model predicts that the stock breached above the 0.2$B_{0}$ limit reference point in about 2010. However, extrapolation beyond the data must always be treated with a great deal of caution. The uncertainty about the catch rates is relatively large, especially because of the records from 1992, 2005, and 2006 deviate so much from the rest of the trend. Without confirmation from other data the predicted recovery from 2007 onwards is purely dirven by the predictions from the fitted model. The predicted recovery should not be accepted on the basis of the model fit alone, there needs to be other data confirming such a recovery.
Whether the predicted recruitment from the model actually happened would require auxilliary information or data to corroborate the prediction. With Orange Roughy, because the fishery was so short lived and they only mature between 31 - 35 years of age, the biology suggested that the fishery is still receiving unfished recruitment levels, which can be expected to decline once the recruits produced from the depleted population are supposed to begin entering the fishery.
However, given the lowest point of the stock is predicted to have occurred in 1996 the expected minimum in recruitment should only occur in about 2026 - 2030.
One way of estimating confidence intervals around the cpue is to use the standard deviation estimates from the likelihood fitting of the CPUE (parameter 2) and set up log-normal error bars, but later we will consider bootstrap percentile confidence intervals as an alternative.
It is not surprising the bounds on the predicted CPUE become vary wide in the regions where there are no cpue data. But even where there are data the bounds are wide. A better estimation of the uncertainty is more likely to be generated by the bootstrap analysis. The variation expressed is primarily driven by the elevated values in 1992 and in 2005 and 2006. The period from 1989 - 2006 is assumed to related to when aggregations were not being targeted but occasionally, no doubt, smaller aggregations would have added heterogeneity to the cpue data. Certainly given this uncertainty it remains questionable whether one could validly claim that the likelihood of the limit reference point being passed was high based only on these data. It would be best to test the robustness of this result to the initial conditions by trialing the model fit using a large array of initial parameter values just to see whether the optimum result was repeatable and stable. We will consider this after introduicing the three parameter models.
A three parameter model
Returning to the dataspm data set it is possible to use a three parameter model to fit this data accepting that the observations began when the stock had already been fished and would be expected to be partially depleted. There are details that need attention if we are to assume such a model structure. These are mathematical models and they can exhibit mathematical behaviour, such as negative recruitment or initial depletions much greater than 1.0, unless such behaviour is controlled. To avoid initial depletions > 1.0 we have implemented a slightly different maximum likelihood function so we should use aspmPENLL instead of aspmLL. In this case, where the initial depletion is estimated to be a long way below 1.0 it might not matter, but such penalties can stabilize even such supposedly safe parameter sets.
data(dataspm)
fish <- dataspm$fish
glb <- dataspm$glb
props <- dataspm$props
pars <- c(14,0.19,0.6) # Fit 3 par__aspm__with penalty
# pars <- c(13.2794439,0.1731744,0.4933178) # for a second time through
scalepar <- magnitude(pars)
bestL <- optim(pars,aspmPENLL,method="Nelder-Mead",
infish=fish,inglb=glb,inprops=props,
control=list(maxit = 1000,parscale=scalepar))
outoptim(bestL)
fisheryPen <- dynamics(bestL$par,infish=fish,inglb=glb,inprops=props)
ceCI <- getLNCI(fisheryPen[,"PredCE"],bestL$par[2])
plotASPM(fisheryPen,CI=ceCI)
Figure 4. The outcome of fitting a three-parameter age-structured production model to the slope fishery data.
Once again the model fails to capture the more rapid changes in the predicted dynamics but does capture the general trends through time (Figure 4). Unlike the two-parameter model it predicts a final depletion close to 50% rather than 60% but this time suggests the starting depletion was about 65% rather than 100%. However, the -ve log-likelihood in this case is -9.95 rather than -7.58 as with the two-parameter model, indicating a slightly better fit (an improvement > 1.96 for each additional parameter suggests a better fit; Venzon and Moolgavkar, 1988).
Despite the improved fit to the data the confidence intervals around the CPUE remain very large. Thus, using the three parameter model one should test the robustness of this fit to the initial conditions to determine whether or not the outcome is stable following the model fitting or whether there is variation (uncertainty) and if so how much.
Testing the Robustness of the Model Fit
The sensitivity of the model fit to the initial parameter values is important to explore to gain greater confidence that the solution one finds is a global optimum rather than some local minima on the log-likelihood surface.
To test for robustness of the model fit we can use the original optimal model parameters or the original guesses, add variation to them, and re-fit the model. This process should enable an analysis of the stability of the modelling outcomes. If the optimum parameters are used then add more variation to enusre the parameter space is covered. The first parameter is $Log \left( R_0 \right)$ so to simplify the selection of random variations away from the original it helps to return that parameter to the linear scale and only when finished return it to the log-scale.
set.seed(12335) # to get repeatable results, normally you would not do this
data(fishdat)
fish <- fishdat$fish
glb <- fishdat$glb
props <- fishdat$props
pars <- c(14,0.3)
out <- robustASPM(pars,fish,glb,props,scaler=20,N=15,console=FALSE)
str(out)
print(round(out$results,4))
Starting with the deep water fishery data fishdat we find that 11 out of 15 generate one solution, which appears to be optimum, while the remaining four, which all began with highly unlikely first guesses (iLike, the initial likelihood was large) all gave implausible outcomes. It would be sensible to explore this lack of robustness further by using many more iterations. However, given the variation in the cpue data this is not a surprising result.
If we test the robustness of the model fit to the dataspm data set (a three parameter model) similar outcomes arise.
set.seed(12235) # to get repeatable results, normally you would not do this
data(dataspm)
fish <- dataspm$fish
glb <- dataspm$glb
props <- dataspm$props
pars <- c(14,0.2,0.6)
out <- robustASPM(pars,fish,glb,props,scaler=15,N=10,console=FALSE)
print(round(out$results,3))
print(round(out$range,3))
Here we find that four final negative log-likelihoods differ from the optimum, although in this case the differences are not too far from the optimum. Very slight differences in the parameters even with the optimum -veLL lead to small differences in the derived statistics such as MSY and $B_0$. Once again the variation in the cpue data is what leads to this instability. Whatever the case it is to be hoped that these examples illustrate that one should never accept the final result of fitting a model even if the diagnostics look acceptable (the plot, the -veLL value, and optim gives convergence = 0). Without testing the robustness it is possible that one is settling for only a local minima. This is one reason why it is usually a good idea to run a fitting routine twice, once from the initial parameter guesses, the second time from the solution of the first time.
When testing the robustness ideally one would run very many trials (at least 100 to allow for proportional attribution of variation), in which case it becomes a reasonable proposition to plot the results. The correlations between the parameters can also be calculated (they tend ot be very high).
cor(out$results[,c("LnR0","Depl","-veLL","MSY")]) # correlations between outputs
#plotprep(width=8,height=6)
intensity <- 2 # how many points overlapping = maximum colour
pairs(out$results[,c("LnR0","Depl","-veLL","MSY")],pch=16,
col=rgb(1,0,0,1/intensity),font=7,font.labels = 7)
Figure 5. The correlations between outputs from repeated trials starting from different initial parameter values. Usually one would use many more trials than the example of 10, then these plots might be more informative. Histograms of these values might also indicate the variation present.
The Production Curve and Statistics
Using two runs through the optim function each time the median of the different trials is very similar to the optimum model fit so we will use those values to determine the production curve predicted by the model We can then use that to estimate the biomass at the target (default = 0.48$B_{0}$) and at the limit reference point of 0.2$B_{0}$. In addition, by estimating the yield expected at those reference points and dividing that through by the biomass at those reference points we can calculate the target and limit harvest rate reference points. The contents of prod can be used to determine other statistics such as the sustainable yield over the range of the current predicted depletion levels.
data(dataspm)
fish <- dataspm$fish
glb <- dataspm$glb
props <- dataspm$props
pars <- c(13.75,0.189667,0.6) # Fit 3 par__aspm__with penalty
bestL <- optim(pars,aspmPENLL,method="Nelder-Mead",
infish=fish,inglb=glb,inprops=props,
control=list(maxit = 1000, parscale = c(10,1,0.1)))
# two times through
bestL <- optim(bestL$par,aspmPENLL,method="Nelder-Mead",
infish=fish,inglb=glb,inprops=props,
control=list(maxit = 1000, parscale = c(10,1,0.1)))
par <- bestL$par
print(par)
prod <- getProductionC(exp(par[1]),fish,glb,props,
Hrg=c(0.01,0.45,0.005),nyr=50)
head(round(prod,3),6)
tail(round(prod,3),6)
anspen <- prodASPM(prod,target=0.48,console=FALSE,plot=TRUE)
round(anspen,3)
Figure 6. Production curves for the optimum fitting three parameter age-structured production model fitted to the slope fishery data in dataspm. The target in this case is 0.48$B_{0}$ designated by the vertical green lines. The results contained within anspen are used as labels. In this case it is suggesting that $B_{MSY}$ is down at 0.243$B_{0}$ so in this case using a target of 0.48$B_{0}$ means that the harvest rate, and presumably effort, would be halved, the stock kept at a much higher presumably more resilient level, and the catch only reduced on average by about 18%.
A Phase Plot
The final part of age-structured production modelling would entail generating a phase plot of predicted biomass against the predicted harvest rates. The previous functions and analyses will provide all the information we require to feed into the function aspmphaseplot.
# plotprep(width=7,height=5.5)
fisheryPen <- dynamics(bestL$par,infish=fish,inglb=glb,inprops=props)
outs <- aspmphaseplot(fisheryPen,prod,anspen,Blim=0.2,fnt=7)
Figure 7. Phase plot of predicted biomass vs predicted harvest rate for the optimum fitting three parameter age-structured production model fitted to the slope fishery data in dataspm. The target in this case at 0.48$B_{0}$ is designated by the green lines, while the limit reference points are designated by the red lines.
The phase plot (Figure 7) suggests that the biomass is a little below the target but the fishing mortality is very close to its target. In addition the fishery appears relatively stable at present indicating it is not declining. In the SAFS system this fishery could defensibly be claimed to be sustainable although the uncertainty in the analysis would need to be noted explicitly.
Characterization of Uncertainty
When only fitting to CPUE it is possible to use many replicate bootstrap samples followed by reanalysis to generate a detailed characterization of uncertainty. The following example code illustrates the approach. First we need to obtain the optimum solution.
# library(datalowSA)
library(datalowSA)
data(dataspm)
fish <- dataspm$fish
glb <- dataspm$glb
props <- dataspm$props
pars <- c(13.5,0.18,0.5)
bestL <- fitASPM(pars,fish,glb,props,callfun=aspmPENLL)
fishery <- dynamics(bestL$par,fish,glb,props)
kable(fishery,digits=c(0,1,1,3,3,3,3,3,3))
Having run the model through optim twice inside fitASPM the optimum fit is used to characterize the dynamics using dynamics. The basis of the bootstrap sampling is that the log-normal residuals (CPUE/PredCE) are randomly sampled with replacement with each such bootstrap then being multiplied by the optimum model's predicted CPUE. If, for example, we take the original residuals and multiply them by the original predicted CPUE we would re-generate the original observed CPUE. All we are doing in the bootstrap procedure is reordering the residuals by randomly resampling them with replacement. The 'with replacement' bit implies that some values may be omitted and others may be repeated more than once.
Such bootstrap samples are generated within bootASPM. This function generates replicate numbers of optimal fitting parameters in param, estimates of unfished biomass in B0, and finally a matrix of five time-series of Spawning Biomass, Fully selected harvest rate, each bootstrap CPUE series, the optimum predicted CPUE, and the depletion level through time. Here we are only running 100 replicates so as to speed the process, but in a real analysis one might use at least 1000 replicates
reps <- 20
starttime <- Sys.time()
answer <- bootASPM(fish,glb,props,bestL$par,iter=reps)
Sys.time() - starttime
str(answer,max.level=1)
Once the bootstraps are completed there are multiple ways of displaying such information. Initially one can generate classical percentile confidence intervals from the bootstrap replicates (Haddon, 2011).
yrs <- fishery[,"Year"]
nyrs <- length(yrs)
par(mfrow=c(2,2),mai=c(0.45,0.45,0.05,0.05))
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
label <- names(answer$result[1,1,])
label <- label[-3] # remove CPUE
numvar <- length(label)
bootvar <- answer$result[,nyrs,label[1]]
for (i in 1:numvar) { # i=3
bootvar <- answer$result[,nyrs,label[i]]
quantCI <- quantile(bootvar,probs=c(0.05,0.5,0.95),na.rm=TRUE)
hist(bootvar,breaks=30,main="",xlab=label[i],col="red")
abline(v=quantCI,col=c(4,4,4),lwd=c(1,2,1))
}
Figure 7. Histograms of the final years' spawning biomass, fully selected harvest rates, predicted CPUE, and the stock depletion level. Of course 20 replicates is completely inadequate but each bootstrap replicate can take a significant time (note the time taken in the example). One thing that can be noted is the asymmetrical percentile confidence bounds.
With only 20 replicates no conclusions can be drawn but the plots still illustrate the principle behind the bootstraps. The percentile confidence intervals can illustrate the uncertainty in the assessments and the potential risk of falling below limit reference points.
pickvar <- "Deplete"
bootvar <- answer$result[,,pickvar]
yrs <- as.numeric(colnames(bootvar))
nyrs <- length(yrs)
quantCI <- t(apply(bootvar,2,quants))
kable(quantCI,digits=c(3,3,3,3,3,3))
ymax <- getmaxy(bootvar)
par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05))
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
plot(yrs,bootvar[1,],type="n",lwd=1,col=0,ylim=c(0,ymax),
panel.first = grid(),xlab="",ylab=pickvar)
for (i in 1:reps) lines(yrs,bootvar[i,],lwd=1,col="grey")
lines(yrs,quantCI[,"50%"],lwd=2,col="red")
arrows(x0=yrs,y0=quantCI[,"5%"],y1=quantCI[,"95%"],
col=2,lwd=1,length=0.035,angle=90,code=3)
Figure 8. The bootstrapped trajectories of stock depletion of the dataspm data set. Note that 20 replicates are far too few to provide sensible or valid percentile confidence intervals.
The output from the bootASPM function includes the bootstrap optimum parameters. These can be used along with the fish, glb, and props objects from the data set used to generate productivity curves and determine target catches, MSY, and other fishery outputs for each set of parameters. This means that percentile confidence intervals can be generated for such assessment outputs.
Management Advice with aspm
If working with a species that requires on-going management then it is necessary to produce advice with respect to acceptable catches that will lead to a sustainable fishery or whatever other management goal is in place for the fishery. To generate such advice formal harvest strategies are required to allow the outputs from the assessment to be converted into a recommended biological catch. This may then be modified by fishery managers taking into account potential rates of change within a fishery or social or economic drivers of management decisions. It was possible to put forward suggestions for new harvest strategies using the catch-MSY method because none were available previously and that put forward was only a suggestion for a possible consideration. Putting forward a proposed harvest control rule for the spm approach without consultation with jurisdictional fisheries managers could produce suggestions incompatible with a particular jurisdictions objectives. There are harvest control rules that can be used once limit and target reference points are agreed upon and these can be utilized where considered appropriate. The SAFS process, however, does not currently require a target reference point even though most harvest control rules do require one.
References
Dick, E.J. and A.D. MacCall (2011) Depletion-based stock reduction analysis: a catch-based method for determining sustainable yields for data-poor fish stocks. Fisheries Research 110(2): 331-341
Haddon, M. (2014) Tier 4 analyses in the SESSF, including deep water species. Data from 1986 – 2012. Pp 352 – 461 in Tuck, G.N. (ed) (2014) Stock Assessment for the Southern and Eastern Scalefish and Shark Fishery 2013. Part 2. Australian Fisheries Management Authority and CSIRO Marine and Atmospheric Research, Hobart. 313p.
Haddon, M., Klaer, N., Wayte, S., and G. Tuck (2015) Options for Tier 5 approaches in the SESSF and identification of when data support for harvest strategies are inappro-priate. CSIRO. FRDC Final Report 2013/200. Hobart. 115p.
Kimura, D.K. and J.V. Tagart (1982) Stock Reduction Analysis, another solution to the catch equations. Canadian Journal of Fisheries and Aquatic Sciences 39: 1467 - 1472.
Kimura, D.K., Balsiger, J.W., and Ito, D.H. 1984. Generalized stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 41: 1325–1333.
Little, L.R., Wayte, S.E., Tuck, G.N., Smith, A.D.M., Klaer, N., Haddon, M., Punt, A.E., Thomson, R., Day, J. and M. Fuller (2011) Development and evaluation of a cpue-based harvest control rule for the southern and eastern scalefish and shark fishery of Australia. ICES Journal of Marine Science 68(8): 1699-1705.
Martell, S. and R. Froese (2013) A simple method for estimating MSY from catch and resilience. Fish and Fisheries 14: 504-514
Venzon, D.J. and S.H. Moolgavkar (1988) A method for computing profile-likelihood-based confidence intervals. Applied Statistics, 37: 87-94.
Punt, A.E., Butterworth, D.S. and A.J. Penney (1995) Stock assessment and risk analysis for the South Atlantic population of albacore Thunnus alalunga using an age-structured production model South African Journal of Marine Science 16: 287-310. http://dx.doi.org/10.2989/025776195784156476
R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. see also https://cran.r-project.org/
RStudio (2016) www.rstudio.com
Schaefer, M.B. (1954) Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Bulletin, Inter-American Tropical Tuna Commission, 1: 25-56.
Schaefer, M.B. (1957) A study of the dynamics of the fishery for yellowfin tuna in the Eastern Tropical Pacific Ocean. Bulletin, Inter-American Tropical Tuna Commission, 2: 247-285
Walters, C.J., Martell, S.J.D. and J. Korman (2006) A stochastic approach to stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 63: 212 - 223.
Appendix: Age-Structured Production Model Equations
Initiation of an Age-Structured Model
At equilibrium, in an un-exploited population, the age-structure is assumed to be the result of natural mortality acting alone upon constant average unfished levels of recruitment. The equilibrium result would be a stable age distribution determined by those constant average recruitments and natural mortality. At the start of a time series, let us say in year 1, this is defined as:
$$
N_{a,1} =
\begin{cases}
N_{0,1} = R_0 & a = 0
\
N_{a-1,1}e^{-M} & 1 \leq a < a_{x}
\
N_{a_{x}-1,1} e^{-M}/(1-e^{-M}) & a = a_{x}
\end{cases}
$$
where $N_{a,1}$ is the numbers of age $a$, in year 1, $a_{x}$ is the maximum age modelled (the plus-group), and $M$ is the instantaneous rate of natural mortality. In a pre-exploitation population there is no fishing mortality and the final component the above equation (where $a$ = $a_{x}$), is referred to as the plus group because it is the series which combines ages $a_x$ and all older ages that are not modelled explicitly. This requires the inclusion of the $(1 – e^{-M})$ divisor to force the equation to be the sum of an exponential series. The $N_{0,1}$ is the constant unfished recruitment level, $R_0$. Sometimes this also has an $e^{-M}$ term, depending on the timing of spawning. If the natural mortality term is included then the estimated $R_0$ value will be somewhat higher than if it is omitted (by $1/e^{-M}$), so it is usually simpler to omit it. This stable age distribution can also be obtained by first calculating the numbers-at-age for a recruitment of 1, or the numbers-at-age per recruit, and then multiplying that vectors of numbers by $R_0$., which is how it is implemented in datalowSA::dynamics
Biological Characteristics
Length-at-age of fish is defined by the von Bertalanffy growth function:
$$ L_a = L_{\infty}{(1-e^{-k(a-t_0)})} $$
where $L_a$ is the mean length at age $a$, $L_{\infty}$ is the
asymptotic average maximum length, $k$ is the grow rate coefficient, and $t_0$
is the length at age zero.
The mass-at-age relationship is defined as:
$$ w_{a} = W_{aa}L^{W_{ab}} $$
where $w_{a}$ is the mass at age $a$, and $W_{aa}$ and $W_{ab}$ are the coefficients that define the power relationship between length and mass.
Spawning Stock Recruitment Relationship
The biomass $A_0$ can be defined as the mature stock biomass that would develop given a constant recruitment level of one (i.e. $N_{0,1} = 1$ in the above equation). Thus, at a biomass of $A_0$, distributed across a stable age distribution, the resulting average recruitment level would be $R_0 = 1$. $A_0$ acts as a scaling factor in the recruitment equations by providing the link between $R_0$ and $B_0$
$$ A_0 = \sum_{a=1}^{a_x} n_{a,1}m_aw_a $$
where $m_i$ is the proportion mature at age $a$, $n_{a,1}$ is the virgin number of animals per recruit of age $a$ in year 1, and $w_a$ is the weight of an animal of age $a$. The average unfished recruitment level, $R_0$, is directly related to the virgin mature, or recruited, biomass, $B_0$
$$ R_0 = B_0/A_0 $$
By determining $A_0$, from a constant recruitment level of one, the recruitment levels from realistic $B_0$ levels can be obtained by applying the above equation. Once $R_0$ has been determined the unfished number at age distribution can be obtained by substituting $R_0$ into the first equation. The spawning stock – recruitment relationship can be described by the deterministic form of the Beverton – Holt relationship:
$$ R_{y+1} = \frac{aB_y^{Sp}}{b+B_y^{Sp}} $$
where $B_y^{Sp}$ is the mature, or spawning biomass in the in year $y$.
A re-parameterization of the Beverton-Holt parameters in terms of steepness, $h$, and $B_0$ is to specify $a$ and $b$ such that:
$$ a=\frac{4hR_0}{5h-1} \qquad \text{and} \qquad b=\frac{B_0(1-h)}{5h-1}$$
Using this re-parameterization the the number of recruits produced in year $y$ from the spawning biomass in year $y-1$ is:
$$ N_{0,y} = \frac{4hR_0B_{y-1}^{Sp}}{(1-h)B_0+(5h-1)B_{y-1}^{Sp}}. $$
Stock dynamics
To describe the dynamics subsequent to population initiation (i.e. the generation of $N_{a,y}$, the number at age $a$ in year $y$, for years other than 0), requires the inclusion of the stock recruitment relationship and the impact of fishing mortality. Not all age classes are necessarily fully selected, thus the fishing mortality term must be multiplied by the selectivity associated with the fishing gear for age $a$, $s_a$, described by a logistic curve:
$$ s_a = \frac{1}{\left( 1+e^{(\frac{a-a_{50}}{\delta})}\right)} $$
where $a_{50}$ is the age at which 50% of individuals are selected by the fishing gear, and $\delta$ is a parameter that determines the width or steepness of the selectivity ogive. Such logistic curves are also used to describe the development of maturity within he population but in such a case the $a_{50}$ refers to the age at 50% maturity.
A term is also needed for the recruitment in each year (stock-recruit relationship above), and this is assumed to be a function of the spawning biomass of the stock at the end of the previous year $y$, $B_y^{Sp}$.
The spawning biomass for a year $y$ is:
$$ B_y^{Sp} = \sum_{a=0}^{a_x} w_a m_a N_{a,y}$$
If this is applied to the unfished stable age distribution this would provide an estimate of the unfished spawning biomass-per-recruit. When using difference equations (rather than continuous differential equations) the dynamics of the fishery, in terms of the order in which growth, natural, and fishing mortality occur, are important when defining how the numbers at age change. If the transition of numbers at age in year $y$ into numbers at age in year $y+1$ is made in a number of steps this simplifies the calculation of internally consistent estimates of exploitable biomass, catch rates, and harvest rates. If it is assumed that the dynamics of a population entails that fish first grow from year $y-1$ to year $y$, then undergo half of natural mortality before they are fished and only then undergo the final half of natural mortality this would imply two steps to define the transition from one year to the next. The first step entails recruitment, growth from each age class to the next, and the application of the effect of half of natural mortality:
$$
N_{a,y^*} =
\begin{cases}
N_{0,y} & a = 0
\
N_{a-1,y-1}e^{-M/2} & 1 \leq a < a_{x} -1
\
\left(N_{a_{x}-1,y-1}+N_{a_{x},y-1}\right) e^{-M/2} & a = a_{x}
\end{cases}
$$
where $N_{0,y}$ is defined by the stock - recruit relationship, ages 1 to $a_x$-1 are modelled by adding 1.0 to the previous year's ages 0 to $a_x$ – 2 and imposing the survivorship from half the natural mortality, and the plus group ($a_x$) is modelled by adding 1.0 to the previous year's age $a_x$ - 1 and adding those to the numbers in the previous year's age $a_x$ and then applying the survivorship from half the natural mortality. The above equation thus leads to the mid-year exploitable biomass (mid-year being the reason for the $e^{-M/2}$) in year $y$ being defined as:
$$ B_y^{E} = \sum_{a=0}^{a_x} w_a s_a N_{a,y^*}$$
The dynamics within any year are completed by the application of the survivorship following fishing mortality across all ages (expressed as an annual harvest rate), followed by the survivorship following the remainder of natural mortality. Natural mortality is not applied directly to the new recruits until they grow into the next year:
$$
N_{a,y} =
\begin{cases}
N_{0,y^} & a = 0
\
N_{a,y^}\left(1-s_a\hat{H_y}\right) e^{-M/2} & 1 \leq a \leq a_{x}
\end{cases}
$$
In the above equation, the $N_{a,y}$ refer the numbers in age $a$ at the end of year $y$ (i.e. after all the dynamics have occurred). The predicted harvest rate, \hat{H_y}., given an observed or recommended catch level in year $y$, $C_y$, is estimated as
$$ \hat{H_y} = \frac{C_y} {B_y^E}$$
where $B_y^E$ is defined above. The catch at age, in numbers, is therefore defined by:
$$ C_{a,y}^N = N_{a,y^*} s_a \hat{H_y}$$
and the total catch by mass is the sum of the separate catches at age multiplied by their respective average weights for all ages:
$$ C_y = \sum_{a=0}^{a_x} w_a C_{a,y}^N$$
Predicted catch rates also derive from the exploitable biomass and the average catchability coefficient, $q$:
$$ I_y = qB_y^E. $$
Likelihoods
Maximum likelihood methods, as the name dictates entail maximizing the likelihood of the available data given the model and a proposed set of parameters. Very often the likelihoods involved when fitting models are very small numbers. To avoid rounding errors (even when using 64 bit computers) it is standard to use log-likelihoods rather than likelihoods (in that way the log-likelihoods can be individually added together rather than multiply the individual likelihoods). Additionally, rather than maximizing a log-likelihood, minimization often best matches our intuitions about model fitting and this involves minimizing the negative log-likelihood. The full log-normal negative log likelihood for the aspm is similar to that used for the spm but with a few parameter changes and it estimates the ${\hat\sigma_{I}}$ directly rather than using a closed form:
$$-veLL\left( data|{R}{0},\hat{\sigma{I}} \right)=-\sum\limits_{t}Ln\left[{\frac{1}{{I}{t}\sqrt{2\pi \hat{\sigma{I} }}}{{e}^{\frac{-\left( Ln{{I}{t}}-Ln{{\hat{I}}{t}} \right)}{2{{{\hat{\sigma_{I} }}}^{2}}}}}} \right]$$
haddonm/datalowSA documentation built on Nov. 5, 2023, 6:40 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE) options(knitr.kable.NA = "", knitr.table.format = "pandoc") options("show.signif.stars"=FALSE,"stringsAsFactors"=FALSE, "max.print"=50000,"width"=240) library(datalowSA) suppressPackageStartupMessages(library(knitr))
Introduction
Which Stock Assessment?
Which stock assessment method to apply to fisheries for data-poor to data-moderate species will depend upon what fisheries and biological data are available but also, importantly, on what management objectives need to be met within the jurisdiction in question. It may be the case that the fishery for a particular species is of sufficient size and value to warrant on-going monitoring and management towards some defined goal for the stock. In such a case the assessment used should obviously be capable of generating some notion of the current state of the fishery and indicate what management actions may be required to eventually achieve the agreed management goals. But some fisheries may be so minor that trying to actively manage them would be inefficient. Nevertheless, to meet the requirements of the Status of key Australian Fish Stocks (SAFS) one still requires some form of defensible stock assessment capable of determining whether the current level of fishing is sustainable.
Age-Structured Production Model
Introduction to ASPM
The age-structured production model (ASPM or aspm) is literally a surplus production model which is based upon an age-structured model of production rather than an accumulated biomass model (see the vignette on spm).
There are some specific data requirements for fitting an age-structured production model to fishery data. Data from the fishery need to included, as a minimum, an accurate catch time-series plus an index of relative abundance for at least some of the years within the catch time-series. In addition, information (or defensible assumptions) is needed for the species concerned in relation to the description of its natural mortality, its growth, its maturation, and the selectivity of the fishery (maturity and selectivity could be knife-edge). If just the catches and CPUE data are available then one might try fitting a simple, aggregated biomass surplus production model. But if these biological data and information are available then an age-structured production model opens the way to on-going improvements with respect to the inclusion of occasional age-composition data or other observations that could be predicted by a suitable model and hence included in the model fitting process.
More details on age-structured production models can be found in Punt et al. (1995). The model equations are provided in the appendix. It would be helpful to the user to read the spm vignette as the theory in there also applies to age-structured production models.
A Typical Workflow
A typical workflow for using an age-structured production model might be something like:
- read in the available data and use checkdata to ensure it can be used with aspm.
- search for suitable initial parameter values using dynamics, aspmLL, and plotASPM, this will include deciding on the use of a two parameter model (no initial depletion) or a three parameter model accounting for an initial depletion level.
- given suitable initial parameters use fitASPM or more basically optim to fit the model to the available data.
- once successfully fitted it is best to plot the solution to determine visually how close the model fit is to the data using plotASPM. One approach to improving this is to include confidence intervals around the index of relative abundance (cpue). This is done first using getLNCI and then including the output of this into plotASPM or, for a closer look in plotceASPM. Eventually improved confidence intervals for the model outputs can be obtained using boostrap samples (see below).
- A better test of the robustness of the solution is to test how sensitive it is to the initial conditions. This can be done by randomly varying the initial parameters and determining whether the model fitting outcomes vary. Suitable example code is given in the vignette.
- After finally deciding on the overall optimum solution it would be sensible to use the optimum parameters to determine the implied production curve so that management statistics such as MSY, the target catch, and the limit and target reference points can be defined. This can be done using getProductionC, the C post-fix denoting this is a C++ function (used for speed), the results of which can be plotted and summarized using prodASPM.
- One can produce a classical phase plot of predicted biomass vs harvest rate to aid in the determination of the stock status, and for this we would use aspmphaseplot.
- Finally, to obtain a better notion of the uncertainty in the analysis, and we can do that using bootASPM, which facilitates the application of a bootstrap analysis.
An Example
The data requirements for the aspm are described above and for the next example we will use the dataspm built in data set. First we will examine an example using only two parameters (assuming the population begins in an unfished state) and then extend the model fitting to include the possibility of an initially depleted state. The two parameters being fitted are the average unfished recruitment level and the standard deviation of the errors around the CPUE data.
We load dataspm, which contains a dataframe of the catches and CPUE by year (fish) and other parameters used for the age-structured production model.
# library(datalowSA) data(dataspm) fish <- dataspm$fish props <- dataspm$props # length-, weight-, maturity- and selectivity-at-age
While a simple, aggregated biomass surplus production model only requires the specification of the species name in the glb data object an aspm requires information on the growth, selectivity, weight-at-age, steepness of the stock-recruit relationship and natural mortality. The global parameters (glb), in addition to the spsname, need to contain the population ages (maxage and ages), natural mortality (M), von Bertalanffy growth parameters (Linf, K and t0), weight-at-age parameters (Waa and Wab), age at 50% maturity, delta (M50a and deltaM), age that 50% of the population are selected by the gear and delta (sela50 and deltaS) and the steepness of the stock-recruitment relationship. The number of years (nyrs) of over which catch and CPUE are available (including missing years) is calculated from the fish dataframe. A starting value for the log of the initial recruitment (R0) also needs to be provided, although this will be estimated along with the standard deviation of the errors around the CPUE data.
We inspect the global parameters specified in dataspm.
(glb <- dataspm$glb)
Just as with the spm model Fitting an aspm model entails first finding initial parameter estimates that lead to predicted cpue time-series that approximate the observed series. Thus, one might begin something like this and uses aspmLL to determine the negative log-likelihood (-verLL), then dynamics to calculate the dynamics of the aspm, and finally plotASPM to illustrate the quality of fit to help decide whether changes are required to the initial guesses.:
pars <- c(12.9,0.25) aspmLL(pars,infish=fish,inglb=glb,inprops=props) fishery <- dynamics(pars,infish=fish,inglb=glb,inprops = props) plotASPM(fishery)
Figure 1. The outcome of a first guess at a two parameter version of the aspm. clearly the fit is very poor and the strong trends in the predicted cpue, the log-normal residuals, and the annual harvest rate bumping up against the built-in upper limit of 0.85 across numerous years, provide a strong indication that the initial guess at unfished average recruitment parameter (R0) is too small. Try increasing it slowly to see its effect on the model fit to the data.
Once reasonable starting values have been found for the parameters ($R_0$ the unfished average recruitment and $\hat\sigma_{I}$, the standard deviation associated with fitting the observed cpue) then an attempt at fitting the model to the data formally can be made using code something like this:
pars <- c(13.7,0.19) ans <- fitASPM(pars,infish=fish,inglb=glb,inprops=props) outoptim(ans) # a tidier way of printing the list output from optim fishery <- dynamics(ans$par,infish=fish,inglb=glb,inprops = props)
Table 1. The output from the fitASPM function and the dynamics function.
kable(fishery,digits=c(0,3,3,3,3,3,4,4,3))
Note the predicted catches are identical to the observed catches. The catches are assumed to be known accurately and so as to ensure a close match between the predicted and observed catches aspmLL has a simple sum-of-squared deviations penalty built into it (try aspmLL, without brackets in the console). This is usually sufficient to force the solution to generate a close match once plausible parameter values are found. Note that with respect to reproduction only the average unfished recruitment is estimated. In its current form the aspm cannot take into accoutn strong and weak cohorts; this remains a very simple model of the dynamics.
To visualize the output of the model fitting we can plot some of the variables from the fishery information generated by the dynamics function.
plotASPM(fishery,CI=NA)
Figure 2. The outcome of fitting a two parameter age-structured production model to the available data. Included is a plot of the catch history, the predicted spawning biomass, the CPUE, the harvest rate, the log-normal residuals for the cpue, and the predicted depletion.
The two-parameter model manages to capture the main cpue trends but fails to capture some of the more obvious and more rapid consistent changes in cpue (Figure 2). The example fishery in question was known to have been fished prior to 1985 so in the year the data begin to be available the stock can be expected to be depleted to some extent. Thus, an alternative might be to fit the model using three parameters. The first approach, with only one parameter of real interest, required data from the beginning of the fishery. However, there are many fisheries for which data are only available after the fishery has been running for a number of years. In such cases it is necessary to estimate the level of depletion and its effect upon recruitment and thus requires two parameters of interest to the dynamics. The first is, as before, the unfished recruitment level, $R_0$, but then we use a parameter that defines the initial depletion at the start of the observations ($D_{init}$. If this parameter is present in pars a search is made for the constant harvest rate that when applied to the initial unfished stock leads to the predicted initial depletion level, only then does the model fitting proceed. Any initial depletion will influence the recruitment depending on the assumed steepness of the stock recruitment relationship, which is assumed to be a Beverton-Holt relationship. The dataspm data set is not particularly suited to a two-parameter model even though that arrangement was able to provide a result; these assessments should not be done automatically, it always takes some background knowledge to ensure that any methods or choices applied are valid.
As an alternative two-parameter example, the deep water fishery data data(fishdat) was a fishery with catch data from the very start of the fishery, which means it is better suited to using a simple two parameter model
data(fishdat) fish <- fishdat$fish glb <- fishdat$glb props <- fishdat$props pars <- c(14,0.3) ans <- fitASPM(pars,infish=fish,inglb=glb,inprops = props) str(ans)
The output contains the estimates of the optimum parameters, the final log-likelihood estimate, the number of iterations it needed to find the optimum, and some diagnostic information. The statement convergence: int 0 implies the solution appears valid, and no warning messages is also encouraging; but see later concerning how to test the robustness of such model fitting. If we put the fitted optimum parameters into the function dynamics we can see the time-series of the more important population variables implied by the model fit.
fishery <- dynamics(ans$par,infish=fish,inglb=glb,inprops = props) print(fishery)
The use of fitASPM is basically shorthand for using bestL <- optim(pars, aspmLL, method="Nelder-Mead", infish=fish, inglb=glb, inprops=props, control=list(maxit = 1000, parscale = c(10,1))) twice in a row. Examine its code by using fitASPM without brackets in the R console window.
Note that the harvest rate in 1996 appears to have bumped up against the upper limit of 0.85 hard wired into the R-code (try typing dynamics into the R console without brackets to see the code). So that 1180 tonnes catch in 1996 was likely to be damaging. However, the predicted catch is only slightly less than the reported catch so this spike in harvest rate has some support in the data. Plotting up the fishery results enables the visualization of the impact of the fishery and he effect of cutting back catches.
ceCI <- getLNCI(fishery[,"PredCE"],ans$par[2]) plotASPM(fishery,CI=ceCI)
Figure 3. The outcome of fitting a two parameter age-structured production model to the deep-water fishery data. Only the CPUE with confidence intervals are plotted.
The modelling suggests that once catches were reduced to an average of about 246 t between 1997 - 2005 the stock began to very slowly recover, then, after catches were further reduced from 2007 onwards (due to a deepwater closure and cessation of targeting) the rate of recover was predicted to have increased. The model predicts that the stock breached above the 0.2$B_{0}$ limit reference point in about 2010. However, extrapolation beyond the data must always be treated with a great deal of caution. The uncertainty about the catch rates is relatively large, especially because of the records from 1992, 2005, and 2006 deviate so much from the rest of the trend. Without confirmation from other data the predicted recovery from 2007 onwards is purely dirven by the predictions from the fitted model. The predicted recovery should not be accepted on the basis of the model fit alone, there needs to be other data confirming such a recovery. Whether the predicted recruitment from the model actually happened would require auxilliary information or data to corroborate the prediction. With Orange Roughy, because the fishery was so short lived and they only mature between 31 - 35 years of age, the biology suggested that the fishery is still receiving unfished recruitment levels, which can be expected to decline once the recruits produced from the depleted population are supposed to begin entering the fishery. However, given the lowest point of the stock is predicted to have occurred in 1996 the expected minimum in recruitment should only occur in about 2026 - 2030.
One way of estimating confidence intervals around the cpue is to use the standard deviation estimates from the likelihood fitting of the CPUE (parameter 2) and set up log-normal error bars, but later we will consider bootstrap percentile confidence intervals as an alternative.
It is not surprising the bounds on the predicted CPUE become vary wide in the regions where there are no cpue data. But even where there are data the bounds are wide. A better estimation of the uncertainty is more likely to be generated by the bootstrap analysis. The variation expressed is primarily driven by the elevated values in 1992 and in 2005 and 2006. The period from 1989 - 2006 is assumed to related to when aggregations were not being targeted but occasionally, no doubt, smaller aggregations would have added heterogeneity to the cpue data. Certainly given this uncertainty it remains questionable whether one could validly claim that the likelihood of the limit reference point being passed was high based only on these data. It would be best to test the robustness of this result to the initial conditions by trialing the model fit using a large array of initial parameter values just to see whether the optimum result was repeatable and stable. We will consider this after introduicing the three parameter models.
A three parameter model
Returning to the dataspm data set it is possible to use a three parameter model to fit this data accepting that the observations began when the stock had already been fished and would be expected to be partially depleted. There are details that need attention if we are to assume such a model structure. These are mathematical models and they can exhibit mathematical behaviour, such as negative recruitment or initial depletions much greater than 1.0, unless such behaviour is controlled. To avoid initial depletions > 1.0 we have implemented a slightly different maximum likelihood function so we should use aspmPENLL instead of aspmLL. In this case, where the initial depletion is estimated to be a long way below 1.0 it might not matter, but such penalties can stabilize even such supposedly safe parameter sets.
data(dataspm) fish <- dataspm$fish glb <- dataspm$glb props <- dataspm$props pars <- c(14,0.19,0.6) # Fit 3 par__aspm__with penalty # pars <- c(13.2794439,0.1731744,0.4933178) # for a second time through scalepar <- magnitude(pars) bestL <- optim(pars,aspmPENLL,method="Nelder-Mead", infish=fish,inglb=glb,inprops=props, control=list(maxit = 1000,parscale=scalepar)) outoptim(bestL) fisheryPen <- dynamics(bestL$par,infish=fish,inglb=glb,inprops=props) ceCI <- getLNCI(fisheryPen[,"PredCE"],bestL$par[2]) plotASPM(fisheryPen,CI=ceCI)
Figure 4. The outcome of fitting a three-parameter age-structured production model to the slope fishery data.
Once again the model fails to capture the more rapid changes in the predicted dynamics but does capture the general trends through time (Figure 4). Unlike the two-parameter model it predicts a final depletion close to 50% rather than 60% but this time suggests the starting depletion was about 65% rather than 100%. However, the -ve log-likelihood in this case is -9.95 rather than -7.58 as with the two-parameter model, indicating a slightly better fit (an improvement > 1.96 for each additional parameter suggests a better fit; Venzon and Moolgavkar, 1988).
Despite the improved fit to the data the confidence intervals around the CPUE remain very large. Thus, using the three parameter model one should test the robustness of this fit to the initial conditions to determine whether or not the outcome is stable following the model fitting or whether there is variation (uncertainty) and if so how much.
Testing the Robustness of the Model Fit
The sensitivity of the model fit to the initial parameter values is important to explore to gain greater confidence that the solution one finds is a global optimum rather than some local minima on the log-likelihood surface.
To test for robustness of the model fit we can use the original optimal model parameters or the original guesses, add variation to them, and re-fit the model. This process should enable an analysis of the stability of the modelling outcomes. If the optimum parameters are used then add more variation to enusre the parameter space is covered. The first parameter is $Log \left( R_0 \right)$ so to simplify the selection of random variations away from the original it helps to return that parameter to the linear scale and only when finished return it to the log-scale.
set.seed(12335) # to get repeatable results, normally you would not do this data(fishdat) fish <- fishdat$fish glb <- fishdat$glb props <- fishdat$props pars <- c(14,0.3) out <- robustASPM(pars,fish,glb,props,scaler=20,N=15,console=FALSE) str(out) print(round(out$results,4))
Starting with the deep water fishery data fishdat we find that 11 out of 15 generate one solution, which appears to be optimum, while the remaining four, which all began with highly unlikely first guesses (iLike, the initial likelihood was large) all gave implausible outcomes. It would be sensible to explore this lack of robustness further by using many more iterations. However, given the variation in the cpue data this is not a surprising result.
If we test the robustness of the model fit to the dataspm data set (a three parameter model) similar outcomes arise.
set.seed(12235) # to get repeatable results, normally you would not do this data(dataspm) fish <- dataspm$fish glb <- dataspm$glb props <- dataspm$props pars <- c(14,0.2,0.6) out <- robustASPM(pars,fish,glb,props,scaler=15,N=10,console=FALSE) print(round(out$results,3)) print(round(out$range,3))
Here we find that four final negative log-likelihoods differ from the optimum, although in this case the differences are not too far from the optimum. Very slight differences in the parameters even with the optimum -veLL lead to small differences in the derived statistics such as MSY and $B_0$. Once again the variation in the cpue data is what leads to this instability. Whatever the case it is to be hoped that these examples illustrate that one should never accept the final result of fitting a model even if the diagnostics look acceptable (the plot, the -veLL value, and optim gives convergence = 0). Without testing the robustness it is possible that one is settling for only a local minima. This is one reason why it is usually a good idea to run a fitting routine twice, once from the initial parameter guesses, the second time from the solution of the first time.
When testing the robustness ideally one would run very many trials (at least 100 to allow for proportional attribution of variation), in which case it becomes a reasonable proposition to plot the results. The correlations between the parameters can also be calculated (they tend ot be very high).
cor(out$results[,c("LnR0","Depl","-veLL","MSY")]) # correlations between outputs #plotprep(width=8,height=6) intensity <- 2 # how many points overlapping = maximum colour pairs(out$results[,c("LnR0","Depl","-veLL","MSY")],pch=16, col=rgb(1,0,0,1/intensity),font=7,font.labels = 7)
Figure 5. The correlations between outputs from repeated trials starting from different initial parameter values. Usually one would use many more trials than the example of 10, then these plots might be more informative. Histograms of these values might also indicate the variation present.
The Production Curve and Statistics
Using two runs through the optim function each time the median of the different trials is very similar to the optimum model fit so we will use those values to determine the production curve predicted by the model We can then use that to estimate the biomass at the target (default = 0.48$B_{0}$) and at the limit reference point of 0.2$B_{0}$. In addition, by estimating the yield expected at those reference points and dividing that through by the biomass at those reference points we can calculate the target and limit harvest rate reference points. The contents of prod can be used to determine other statistics such as the sustainable yield over the range of the current predicted depletion levels.
data(dataspm) fish <- dataspm$fish glb <- dataspm$glb props <- dataspm$props pars <- c(13.75,0.189667,0.6) # Fit 3 par__aspm__with penalty bestL <- optim(pars,aspmPENLL,method="Nelder-Mead", infish=fish,inglb=glb,inprops=props, control=list(maxit = 1000, parscale = c(10,1,0.1))) # two times through bestL <- optim(bestL$par,aspmPENLL,method="Nelder-Mead", infish=fish,inglb=glb,inprops=props, control=list(maxit = 1000, parscale = c(10,1,0.1))) par <- bestL$par print(par) prod <- getProductionC(exp(par[1]),fish,glb,props, Hrg=c(0.01,0.45,0.005),nyr=50) head(round(prod,3),6) tail(round(prod,3),6) anspen <- prodASPM(prod,target=0.48,console=FALSE,plot=TRUE) round(anspen,3)
Figure 6. Production curves for the optimum fitting three parameter age-structured production model fitted to the slope fishery data in dataspm. The target in this case is 0.48$B_{0}$ designated by the vertical green lines. The results contained within anspen are used as labels. In this case it is suggesting that $B_{MSY}$ is down at 0.243$B_{0}$ so in this case using a target of 0.48$B_{0}$ means that the harvest rate, and presumably effort, would be halved, the stock kept at a much higher presumably more resilient level, and the catch only reduced on average by about 18%.
A Phase Plot
The final part of age-structured production modelling would entail generating a phase plot of predicted biomass against the predicted harvest rates. The previous functions and analyses will provide all the information we require to feed into the function aspmphaseplot.
# plotprep(width=7,height=5.5) fisheryPen <- dynamics(bestL$par,infish=fish,inglb=glb,inprops=props) outs <- aspmphaseplot(fisheryPen,prod,anspen,Blim=0.2,fnt=7)
Figure 7. Phase plot of predicted biomass vs predicted harvest rate for the optimum fitting three parameter age-structured production model fitted to the slope fishery data in dataspm. The target in this case at 0.48$B_{0}$ is designated by the green lines, while the limit reference points are designated by the red lines.
The phase plot (Figure 7) suggests that the biomass is a little below the target but the fishing mortality is very close to its target. In addition the fishery appears relatively stable at present indicating it is not declining. In the SAFS system this fishery could defensibly be claimed to be sustainable although the uncertainty in the analysis would need to be noted explicitly.
Characterization of Uncertainty
When only fitting to CPUE it is possible to use many replicate bootstrap samples followed by reanalysis to generate a detailed characterization of uncertainty. The following example code illustrates the approach. First we need to obtain the optimum solution.
# library(datalowSA) library(datalowSA) data(dataspm) fish <- dataspm$fish glb <- dataspm$glb props <- dataspm$props pars <- c(13.5,0.18,0.5) bestL <- fitASPM(pars,fish,glb,props,callfun=aspmPENLL) fishery <- dynamics(bestL$par,fish,glb,props) kable(fishery,digits=c(0,1,1,3,3,3,3,3,3))
Having run the model through optim twice inside fitASPM the optimum fit is used to characterize the dynamics using dynamics. The basis of the bootstrap sampling is that the log-normal residuals (CPUE/PredCE) are randomly sampled with replacement with each such bootstrap then being multiplied by the optimum model's predicted CPUE. If, for example, we take the original residuals and multiply them by the original predicted CPUE we would re-generate the original observed CPUE. All we are doing in the bootstrap procedure is reordering the residuals by randomly resampling them with replacement. The 'with replacement' bit implies that some values may be omitted and others may be repeated more than once.
Such bootstrap samples are generated within bootASPM. This function generates replicate numbers of optimal fitting parameters in param, estimates of unfished biomass in B0, and finally a matrix of five time-series of Spawning Biomass, Fully selected harvest rate, each bootstrap CPUE series, the optimum predicted CPUE, and the depletion level through time. Here we are only running 100 replicates so as to speed the process, but in a real analysis one might use at least 1000 replicates
reps <- 20 starttime <- Sys.time() answer <- bootASPM(fish,glb,props,bestL$par,iter=reps) Sys.time() - starttime str(answer,max.level=1)
Once the bootstraps are completed there are multiple ways of displaying such information. Initially one can generate classical percentile confidence intervals from the bootstrap replicates (Haddon, 2011).
yrs <- fishery[,"Year"] nyrs <- length(yrs) par(mfrow=c(2,2),mai=c(0.45,0.45,0.05,0.05)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) label <- names(answer$result[1,1,]) label <- label[-3] # remove CPUE numvar <- length(label) bootvar <- answer$result[,nyrs,label[1]] for (i in 1:numvar) { # i=3 bootvar <- answer$result[,nyrs,label[i]] quantCI <- quantile(bootvar,probs=c(0.05,0.5,0.95),na.rm=TRUE) hist(bootvar,breaks=30,main="",xlab=label[i],col="red") abline(v=quantCI,col=c(4,4,4),lwd=c(1,2,1)) }
Figure 7. Histograms of the final years' spawning biomass, fully selected harvest rates, predicted CPUE, and the stock depletion level. Of course 20 replicates is completely inadequate but each bootstrap replicate can take a significant time (note the time taken in the example). One thing that can be noted is the asymmetrical percentile confidence bounds.
With only 20 replicates no conclusions can be drawn but the plots still illustrate the principle behind the bootstraps. The percentile confidence intervals can illustrate the uncertainty in the assessments and the potential risk of falling below limit reference points.
pickvar <- "Deplete" bootvar <- answer$result[,,pickvar] yrs <- as.numeric(colnames(bootvar)) nyrs <- length(yrs) quantCI <- t(apply(bootvar,2,quants)) kable(quantCI,digits=c(3,3,3,3,3,3))
ymax <- getmaxy(bootvar) par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) plot(yrs,bootvar[1,],type="n",lwd=1,col=0,ylim=c(0,ymax), panel.first = grid(),xlab="",ylab=pickvar) for (i in 1:reps) lines(yrs,bootvar[i,],lwd=1,col="grey") lines(yrs,quantCI[,"50%"],lwd=2,col="red") arrows(x0=yrs,y0=quantCI[,"5%"],y1=quantCI[,"95%"], col=2,lwd=1,length=0.035,angle=90,code=3)
Figure 8. The bootstrapped trajectories of stock depletion of the dataspm data set. Note that 20 replicates are far too few to provide sensible or valid percentile confidence intervals.
The output from the bootASPM function includes the bootstrap optimum parameters. These can be used along with the fish, glb, and props objects from the data set used to generate productivity curves and determine target catches, MSY, and other fishery outputs for each set of parameters. This means that percentile confidence intervals can be generated for such assessment outputs.
Management Advice with aspm
If working with a species that requires on-going management then it is necessary to produce advice with respect to acceptable catches that will lead to a sustainable fishery or whatever other management goal is in place for the fishery. To generate such advice formal harvest strategies are required to allow the outputs from the assessment to be converted into a recommended biological catch. This may then be modified by fishery managers taking into account potential rates of change within a fishery or social or economic drivers of management decisions. It was possible to put forward suggestions for new harvest strategies using the catch-MSY method because none were available previously and that put forward was only a suggestion for a possible consideration. Putting forward a proposed harvest control rule for the spm approach without consultation with jurisdictional fisheries managers could produce suggestions incompatible with a particular jurisdictions objectives. There are harvest control rules that can be used once limit and target reference points are agreed upon and these can be utilized where considered appropriate. The SAFS process, however, does not currently require a target reference point even though most harvest control rules do require one.
References
Dick, E.J. and A.D. MacCall (2011) Depletion-based stock reduction analysis: a catch-based method for determining sustainable yields for data-poor fish stocks. Fisheries Research 110(2): 331-341
Haddon, M. (2014) Tier 4 analyses in the SESSF, including deep water species. Data from 1986 – 2012. Pp 352 – 461 in Tuck, G.N. (ed) (2014) Stock Assessment for the Southern and Eastern Scalefish and Shark Fishery 2013. Part 2. Australian Fisheries Management Authority and CSIRO Marine and Atmospheric Research, Hobart. 313p.
Haddon, M., Klaer, N., Wayte, S., and G. Tuck (2015) Options for Tier 5 approaches in the SESSF and identification of when data support for harvest strategies are inappro-priate. CSIRO. FRDC Final Report 2013/200. Hobart. 115p.
Kimura, D.K. and J.V. Tagart (1982) Stock Reduction Analysis, another solution to the catch equations. Canadian Journal of Fisheries and Aquatic Sciences 39: 1467 - 1472.
Kimura, D.K., Balsiger, J.W., and Ito, D.H. 1984. Generalized stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 41: 1325–1333.
Little, L.R., Wayte, S.E., Tuck, G.N., Smith, A.D.M., Klaer, N., Haddon, M., Punt, A.E., Thomson, R., Day, J. and M. Fuller (2011) Development and evaluation of a cpue-based harvest control rule for the southern and eastern scalefish and shark fishery of Australia. ICES Journal of Marine Science 68(8): 1699-1705.
Martell, S. and R. Froese (2013) A simple method for estimating MSY from catch and resilience. Fish and Fisheries 14: 504-514
Venzon, D.J. and S.H. Moolgavkar (1988) A method for computing profile-likelihood-based confidence intervals. Applied Statistics, 37: 87-94.
Punt, A.E., Butterworth, D.S. and A.J. Penney (1995) Stock assessment and risk analysis for the South Atlantic population of albacore Thunnus alalunga using an age-structured production model South African Journal of Marine Science 16: 287-310. http://dx.doi.org/10.2989/025776195784156476
R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. see also https://cran.r-project.org/
RStudio (2016) www.rstudio.com
Schaefer, M.B. (1954) Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Bulletin, Inter-American Tropical Tuna Commission, 1: 25-56.
Schaefer, M.B. (1957) A study of the dynamics of the fishery for yellowfin tuna in the Eastern Tropical Pacific Ocean. Bulletin, Inter-American Tropical Tuna Commission, 2: 247-285
Walters, C.J., Martell, S.J.D. and J. Korman (2006) A stochastic approach to stock reduction analysis. Canadian Journal of Fisheries and Aquatic Sciences 63: 212 - 223.
Appendix: Age-Structured Production Model Equations
Initiation of an Age-Structured Model
At equilibrium, in an un-exploited population, the age-structure is assumed to be the result of natural mortality acting alone upon constant average unfished levels of recruitment. The equilibrium result would be a stable age distribution determined by those constant average recruitments and natural mortality. At the start of a time series, let us say in year 1, this is defined as:
$$ N_{a,1} = \begin{cases} N_{0,1} = R_0 & a = 0 \ N_{a-1,1}e^{-M} & 1 \leq a < a_{x} \ N_{a_{x}-1,1} e^{-M}/(1-e^{-M}) & a = a_{x} \end{cases} $$
where $N_{a,1}$ is the numbers of age $a$, in year 1, $a_{x}$ is the maximum age modelled (the plus-group), and $M$ is the instantaneous rate of natural mortality. In a pre-exploitation population there is no fishing mortality and the final component the above equation (where $a$ = $a_{x}$), is referred to as the plus group because it is the series which combines ages $a_x$ and all older ages that are not modelled explicitly. This requires the inclusion of the $(1 – e^{-M})$ divisor to force the equation to be the sum of an exponential series. The $N_{0,1}$ is the constant unfished recruitment level, $R_0$. Sometimes this also has an $e^{-M}$ term, depending on the timing of spawning. If the natural mortality term is included then the estimated $R_0$ value will be somewhat higher than if it is omitted (by $1/e^{-M}$), so it is usually simpler to omit it. This stable age distribution can also be obtained by first calculating the numbers-at-age for a recruitment of 1, or the numbers-at-age per recruit, and then multiplying that vectors of numbers by $R_0$., which is how it is implemented in datalowSA::dynamics
Biological Characteristics
Length-at-age of fish is defined by the von Bertalanffy growth function:
$$ L_a = L_{\infty}{(1-e^{-k(a-t_0)})} $$
where $L_a$ is the mean length at age $a$, $L_{\infty}$ is the asymptotic average maximum length, $k$ is the grow rate coefficient, and $t_0$ is the length at age zero.
The mass-at-age relationship is defined as:
$$ w_{a} = W_{aa}L^{W_{ab}} $$
where $w_{a}$ is the mass at age $a$, and $W_{aa}$ and $W_{ab}$ are the coefficients that define the power relationship between length and mass.
Spawning Stock Recruitment Relationship
The biomass $A_0$ can be defined as the mature stock biomass that would develop given a constant recruitment level of one (i.e. $N_{0,1} = 1$ in the above equation). Thus, at a biomass of $A_0$, distributed across a stable age distribution, the resulting average recruitment level would be $R_0 = 1$. $A_0$ acts as a scaling factor in the recruitment equations by providing the link between $R_0$ and $B_0$
$$ A_0 = \sum_{a=1}^{a_x} n_{a,1}m_aw_a $$
where $m_i$ is the proportion mature at age $a$, $n_{a,1}$ is the virgin number of animals per recruit of age $a$ in year 1, and $w_a$ is the weight of an animal of age $a$. The average unfished recruitment level, $R_0$, is directly related to the virgin mature, or recruited, biomass, $B_0$
$$ R_0 = B_0/A_0 $$
By determining $A_0$, from a constant recruitment level of one, the recruitment levels from realistic $B_0$ levels can be obtained by applying the above equation. Once $R_0$ has been determined the unfished number at age distribution can be obtained by substituting $R_0$ into the first equation. The spawning stock – recruitment relationship can be described by the deterministic form of the Beverton – Holt relationship:
$$ R_{y+1} = \frac{aB_y^{Sp}}{b+B_y^{Sp}} $$ where $B_y^{Sp}$ is the mature, or spawning biomass in the in year $y$.
A re-parameterization of the Beverton-Holt parameters in terms of steepness, $h$, and $B_0$ is to specify $a$ and $b$ such that:
$$ a=\frac{4hR_0}{5h-1} \qquad \text{and} \qquad b=\frac{B_0(1-h)}{5h-1}$$
Using this re-parameterization the the number of recruits produced in year $y$ from the spawning biomass in year $y-1$ is:
$$ N_{0,y} = \frac{4hR_0B_{y-1}^{Sp}}{(1-h)B_0+(5h-1)B_{y-1}^{Sp}}. $$
Stock dynamics
To describe the dynamics subsequent to population initiation (i.e. the generation of $N_{a,y}$, the number at age $a$ in year $y$, for years other than 0), requires the inclusion of the stock recruitment relationship and the impact of fishing mortality. Not all age classes are necessarily fully selected, thus the fishing mortality term must be multiplied by the selectivity associated with the fishing gear for age $a$, $s_a$, described by a logistic curve:
$$ s_a = \frac{1}{\left( 1+e^{(\frac{a-a_{50}}{\delta})}\right)} $$
where $a_{50}$ is the age at which 50% of individuals are selected by the fishing gear, and $\delta$ is a parameter that determines the width or steepness of the selectivity ogive. Such logistic curves are also used to describe the development of maturity within he population but in such a case the $a_{50}$ refers to the age at 50% maturity.
A term is also needed for the recruitment in each year (stock-recruit relationship above), and this is assumed to be a function of the spawning biomass of the stock at the end of the previous year $y$, $B_y^{Sp}$.
The spawning biomass for a year $y$ is:
$$ B_y^{Sp} = \sum_{a=0}^{a_x} w_a m_a N_{a,y}$$
If this is applied to the unfished stable age distribution this would provide an estimate of the unfished spawning biomass-per-recruit. When using difference equations (rather than continuous differential equations) the dynamics of the fishery, in terms of the order in which growth, natural, and fishing mortality occur, are important when defining how the numbers at age change. If the transition of numbers at age in year $y$ into numbers at age in year $y+1$ is made in a number of steps this simplifies the calculation of internally consistent estimates of exploitable biomass, catch rates, and harvest rates. If it is assumed that the dynamics of a population entails that fish first grow from year $y-1$ to year $y$, then undergo half of natural mortality before they are fished and only then undergo the final half of natural mortality this would imply two steps to define the transition from one year to the next. The first step entails recruitment, growth from each age class to the next, and the application of the effect of half of natural mortality:
$$ N_{a,y^*} = \begin{cases} N_{0,y} & a = 0 \ N_{a-1,y-1}e^{-M/2} & 1 \leq a < a_{x} -1 \ \left(N_{a_{x}-1,y-1}+N_{a_{x},y-1}\right) e^{-M/2} & a = a_{x} \end{cases} $$
where $N_{0,y}$ is defined by the stock - recruit relationship, ages 1 to $a_x$-1 are modelled by adding 1.0 to the previous year's ages 0 to $a_x$ – 2 and imposing the survivorship from half the natural mortality, and the plus group ($a_x$) is modelled by adding 1.0 to the previous year's age $a_x$ - 1 and adding those to the numbers in the previous year's age $a_x$ and then applying the survivorship from half the natural mortality. The above equation thus leads to the mid-year exploitable biomass (mid-year being the reason for the $e^{-M/2}$) in year $y$ being defined as:
$$ B_y^{E} = \sum_{a=0}^{a_x} w_a s_a N_{a,y^*}$$
The dynamics within any year are completed by the application of the survivorship following fishing mortality across all ages (expressed as an annual harvest rate), followed by the survivorship following the remainder of natural mortality. Natural mortality is not applied directly to the new recruits until they grow into the next year:
$$ N_{a,y} = \begin{cases} N_{0,y^} & a = 0 \ N_{a,y^}\left(1-s_a\hat{H_y}\right) e^{-M/2} & 1 \leq a \leq a_{x} \end{cases} $$ In the above equation, the $N_{a,y}$ refer the numbers in age $a$ at the end of year $y$ (i.e. after all the dynamics have occurred). The predicted harvest rate, \hat{H_y}., given an observed or recommended catch level in year $y$, $C_y$, is estimated as
$$ \hat{H_y} = \frac{C_y} {B_y^E}$$
where $B_y^E$ is defined above. The catch at age, in numbers, is therefore defined by:
$$ C_{a,y}^N = N_{a,y^*} s_a \hat{H_y}$$
and the total catch by mass is the sum of the separate catches at age multiplied by their respective average weights for all ages:
$$ C_y = \sum_{a=0}^{a_x} w_a C_{a,y}^N$$
Predicted catch rates also derive from the exploitable biomass and the average catchability coefficient, $q$:
$$ I_y = qB_y^E. $$
Likelihoods
Maximum likelihood methods, as the name dictates entail maximizing the likelihood of the available data given the model and a proposed set of parameters. Very often the likelihoods involved when fitting models are very small numbers. To avoid rounding errors (even when using 64 bit computers) it is standard to use log-likelihoods rather than likelihoods (in that way the log-likelihoods can be individually added together rather than multiply the individual likelihoods). Additionally, rather than maximizing a log-likelihood, minimization often best matches our intuitions about model fitting and this involves minimizing the negative log-likelihood. The full log-normal negative log likelihood for the aspm is similar to that used for the spm but with a few parameter changes and it estimates the ${\hat\sigma_{I}}$ directly rather than using a closed form:
$$-veLL\left( data|{R}{0},\hat{\sigma{I}} \right)=-\sum\limits_{t}Ln\left[{\frac{1}{{I}{t}\sqrt{2\pi \hat{\sigma{I} }}}{{e}^{\frac{-\left( Ln{{I}{t}}-Ln{{\hat{I}}{t}} \right)}{2{{{\hat{\sigma_{I} }}}^{2}}}}}} \right]$$
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.