03: Indicators
FAO: Small Pelagics
L Kell
21 January, 2024
Jump to Empirical Harvest Control Rules
Introduction
This vignette provides examples of empirical indicators and how they can be used to provide catch advice as part of a Harvest Control Rule (HCR). The application of Receiver Operator Characteristics [ROC, @green1966signal] for selecting and indicators and tuning HCRs is also demonstrated. ROC analysis provides a methodological framework to ensure that indicators are effective in classification of stock status, but also sensitive to changes in trends, thereby improving decision-making in fisheries management.
For data-limited stocks empirical indicators, e.g. indices of abundance based on a survey or catch per unit effort (CPUE), are often used to provide management advice through the implementation of Harvest Control Rules (HCR).
Length-based indicators (LBIs) and reference points based on life history parameters have also been proposed for categorising stocks according to conservation and yield objectives [@10.1093/icesjms/fsac043]. For example,mean size, can be used to identifying growth over-fishing, as such stocks typically exhibit a lower proportion of older larger fish in the population, than for a stock being fished sustainably.
Empirical Harvest Control Rules
Emprirical indices can be used in simple HCRs, for example the “2-over-3” rule, where the last two values are compared to the preceding three, or by the application of a linear regression trend (Geromont and Butterworth 2015) .
X-Over-Y Rule
For data-limited stocks, empirical indicators such as survey indices or catch per unit effort (CPUE) can be used to set catches. For example, the “2-over-3” use the ratio of the average of the most recent 2 years to the average of the preceding 3 years, if the ratio is >1 then this implies the stock is increasing, and so an increase in catch is allowed.
The rule aims to keep stocks at their current level by multiplying recent catches by the trend in a biomass index:
\begin{equation}TAC_{t+1}=TAC_t−\frac{\sum\limits_{i=X}^{t-X+i} (I_i)/X}{\sum\limits_{i=Y}^{t-X-Y+i} (I_i)/Y}\end{equation}
where
\(TAC\_{+1}\) is the newly advised catch for year \(t+1\), \({TACy_t}\) is the previously advised catch, and \(I\) is an index of relative abundance.The rule can also be used for a reference period, rather than making \(Y\) a sliding window, as it can also be set it to a period of years (e.g. N years ago) when the stock was thought to be at a target level.
\begin{equation}TAC_{t+1}=TAC_t−\frac{\sum\limits_{i=X}^{t-X+i} (I_i)/X}{\sum\limits_{i=Y-N}^{t-X-Y+i} (I_i)/Y}\end{equation}
While these rules are straightforward and intuitive, guiding catch adjustments based on recent data trends, they do have limitations. Since without a specific target in mind, they are insufficient to effectively rebuild stock biomass when it is at low depletion levels.
Linear trend
The trend for the recent period (\(s\)) can be estimated from the slope of the linear regression of \(ln I_y\) against \(y^\prime\) for years \(y^\prime = y-p+1, y-p+2,...,y\) for an abundance index \(I\), where \(p\) is the number of years over which the slope is calculated.
For example the rule of Hillary, Ann Preece, and Davies (2013), which takes as input an index of abundance (\(I_{t}\)): the total allowable catch next year (\(TAC_{t+1}\)) is increased relative to the catch this year (\(TAC_t\)) if the trend is positive and decreased if the trend is negative.
\[\begin{equation} TAC_{t+1}=TAC_t\times \left\{\begin{array}{rcl} {1-k_1|\lambda|} & \mbox{for} & \lambda<0\\[0.35cm] {1+k_2\lambda} & \mbox{for} & \lambda\geq 0 \end{array}\right. \label{eqn:bft} \end{equation}\]
where \(\lambda\) is the slope in the regression of \(\ln I_t\) against year for the most recent \(n\) years and \(k_1\) and \(k_2\) are parameters. The rule actions asymmetry, as decreases in the index, do not necessarily result in the same relative change as an increase. The control parameters values (\(k_1\),\(k_2\)) are tuned, i.e. calibrated, to meet the management objectives.
If \(p\) is too small the trend estimates could fluctuate too much (tracking noise), but if too large, a HCR would not be able to react sufficiently rapidly to recent trends in resource abundance.
Composite rules
Rules can also combine different inducators, e.g. the \(rfb\) rule of Fischer, José, and Kell (2020), used by ICES for data poor stocks, combines length indicators and indices of abundance into a single rule
\[\begin{equation} A_y = C_{y-1}\times rfb \end{equation}\]
Where catch advice (\(A_{y}\)) based on the previous catch (\(C_{y-1}\)) multiplied by three components, \(r\) an index of biomass, \(f\) a LBI proxy for \(F\), and \(b\) which reduces the catch when the biomass index falls below a threshold.
The different components (\(r\), \(f\), \(b\)) are weighted reflecting their ability to track stock status, and to meet managemenet objectives.
Simulation
When selecting indicators or model-based trends for Harvest Control Rules (HCRs), it is crucial to demonstrate their ability to correctly classify stock status (as being overfished or experiencing overfishing) and detect trends. The selection and tuning process is vital for effective fisheries management and for achieving objectives. Therefore to assess indicator performance, we use a simulation-based approach.
For data-limited stocks lacking comprehensive assessments, we condition an operating model with life history information and employ an observation error model to simulate indicators, e.g. an index of abundance and a length-based indicator. We then compare this simulated indicators to the actual quantities.
To do this we
1 run an Operating Model without feedback, based on life history theory, to replicate fishery system dynamics.
2 generate empirical indices of abundance and exploitation level using an Observation Error Model.
3 use ROC analysis to evaluate how well the indicators classify stock status and detect trends.
This analysis helps identify the optimal combination of index and reference level for these tasks.
Packages
Load the required packages
FLR
library(FLCore)
library(FLBRP)
library(FLasher)
library(FLife)
library(ggplotFL)
library(FLCandy)FishLife
library(rfishbase)
library(SPMpriors)
library(FishLife)Plotting and Data wrangling
library(plyr)
library(dplyr)
library(reshape)library(ggplot2)
theme_set(theme_bw())
library(ggpubr)
library(ggcorrplot)
library(GGally)Operating Model Conditioning
Use life history parameters to condition a FLBRP
equilibrium object
par=lhPar(FLPar(
c(linf= 59.1, k =0.28, t0=-0.4,
a =0.01, b =3.2,
a50 =4.0, l50=43.25,
s =0.7),
units="NA"))parAn object of class "FLPar"
params
linf k t0 a b l50 a50 ato95 asym bg
59.10 0.28 -0.40 0.01 3.20 43.25 4.00 1.00 1.00 3.00
m1 m2 m3 s v sel1 sel2 sel3
0.55 -1.61 1.44 0.70 1000.00 5.00 1.00 5000.00
units: NA Create an FLBRP object
eq=lhEql(par)
plot(eq,refpts="msy")
Figure 1 Equilibrium dynamics as modelled by an
FLBRP object.
Set variation in \(F\), increasing and decreasing trends, to provide contrast
First set \(F=F_{MSY}\) for initial conditions
fbar(eq)=propagate(FLQuant(rep(1,50))%*%refpts(eq)["msy","harvest"],101)Set the \(F\) level to vary from 20% to 400% of \(F_{MSY}\) at the end of the time series, i.e. terminal \(F\)s.
Frange=c(seq(0.2,1,length.out=51),seq(1,4,length.out=51)[-1])
plot(Frange)
Figure 2 Terminal \(F\)s.
Fs=maply(Frange, function(x) seq(1,x,length.out=30))
Fs=t(Fs)
Ftrend=FLQuant(c(Fs),dimnames=list(year=21:50,iter=seq(101)))fbar(eq)=propagate(fbar(eq),101)
fbar(eq)[,ac(21:50)]=fbar(eq)[,ac(21:50)]%*%Ftrendplot(fbar(eq)%/%refpts(eq)["msy","harvest"],iter=c(17,67))+
ylab("F/Fmsy")+xlab("Year")
Figure 3 Trends in \(F\)s, CIs, with 2 example iterations.
Convert the equilibrium object to an FLStock to model
the time series dynamics
set.seed(234)
srDev =rlnoise(101,fbar(eq)%=%0,0.3,0)
om=as(eq,"FLStock")
om=propagate(om,dim(srDev)[6])Project so that the dynamics are correct, i.e. burn in.
om=ffwd(om,fbar=fbar(eq)[,-1],sr=eq,deviances=srDev)
om=window(om,start=11)[1] "maxfbar has been changed to accomodate new plusgroup"Summarise the dynamics, by creating FLQuants with
quantities of interest, i.e. divide the time series by \(MSY\) benchmarks.
metrics <- list(F =function(x, y) fbar( x)/ fmsy(y),
SSB =function(x, y) ssb( x)/ sbmsy(y),
Recruits=function(x, y) rec( x)/refpts(y)["msy","rec"],
Yield =function(x, y) catch(x)/refpts(y)["msy","yield"])
OM <- FLQuants(lapply(metrics, function(fn) fn(om, eq)))plot(OM,iter=1)+
geom_hline(yintercept=1, linetype=2)+
ylim(c(0, NA))+
theme_bw(16)+
theme(legend.position="none")
Figure 4 Operating Model.
Indicators
Empirical indices of abundance and exploitation level are simulated using an Observation Error Model.
Index of abundance
idx=ssb(om)%*%rlnorm(dim(ssb(om))[6],ssb(om)%=%0,0.3)plot(idx,iter=17)
Figure 5 An index of abundance with a worm, i.e. a single iteration.
Length-based indicator
lbi=FLQuants(indicators.len(om,
indicators='lmean',
params=FLPar(apply(par,1,median)),
metric='catch.n'))plot(lbi,iter=17)
Figure 6 An length-based indicator with a worm, i.e. a single iteration.
Classification Skill
Receiver Operating Characteristics (ROC) analysis is used to estimate the ability of indicators to assess status and trends. ROC curves, summarise the true positive rate (TPR) and the false positive rate (FPR). For example, in this case the TPR is the probability that a stock is overfished and the indicator correctly identifies it, while the FPR is when the fishing is sustainable but the indicator incorrectly assesses the stock as being overfished.
Index of Abundance
First get the true state, i.e. SSB
State =ssb(om)[,ac(21:40)]%/%refpts(eq)["msy","ssb"]Then the abundance index relative to a reference year, i.e. year 11 when \(F=F_{MSY}\)
Plotting the indicator against the true stock status, points in the top right-hand categorises classifications as true positives (TP), true negatives (TN), false positives (FP), or false negatives (FN). The plot is therefore a visual representation of a confusion matrix by
Figure 7. Pairwise scatter plot of index relative to reference verses \(SSB/B_{MSY}\)
ROC Curve
Receiver Operating Characteristics (ROC) can assess indicator performance in evaluating stock status. ROC curves sort indicator values, such as for overfishing (F/F_{MSY}>1) from high to low and for overfished (B/B_{MSY}<1) from low to high, comparing them to the known state from the Operating Model. Cumulative TPR and TNR are calculated for the ordered observed outcomes. TPR is then plotted against FPR for the ordered indicator. ROC curves represent power versus Type I Error, aiding in indicator selection and calibration by choosing the reference level with the best classification skill.
head(idxSkill) year iter ind label TP TN FP FN TPR FPR TSS
256 36 13 5.19 1 1 973 0 1046 0.000955 0.00000 0.000955
255 35 13 4.88 1 2 973 0 1045 0.001910 0.00000 0.001910
241 21 13 4.85 1 3 973 0 1044 0.002865 0.00000 0.002865
1630 30 82 4.78 0 3 972 1 1044 0.002865 0.00103 0.001838
250 30 13 4.46 1 4 972 1 1043 0.003820 0.00103 0.002793
36 36 2 4.15 1 5 972 1 1042 0.004776 0.00103 0.003748
Figure 8. ROC curve, black point is when index/reference =1, and blue is the point with the best skill score.
The ROC curve, a probability curve, uses the Area Under the Curve (AUC) as a key performance metric. For example, a coin toss generates a curve along the y=x line, resulting in an AUC of 0.5. The closer the AUC is to 1, the better the indicator ranks. The ROC curve also visually identifies the performance of a chosen indicator ratio (discriminant threshold).
The true skill statistic (TSS) to summarize classification skill, calculated as the sum of the true positive rate (TPR = TP / (TP + FP)) and the false positive rate (FPR = FP / (TP + FP)). A random guess results in a TSS of zero, perfect classifications score 1, random ones score 0, and predictions worse than random receive a negative score. Calibration, which maximizes the sum of TPR (sensitivity) and true negative rate (TNR, specificity), helps identify the optimal threshold. The second panel plots TSS against the reference level, negatively biased because increasing the reference level enhances TP and TSS.
Figure 9. Plot of TSS v index value.
Risks are asymmetric, meaning the risk of incorrectly indicating overfishing when the stock is sustainably exploited differs from the risk of failing to detect overfishing. Adjusting the threshold can alter sensitivity to false positives or negatives. Some indicators excel at detecting overfishing onset, while others perform better at identifying recovery. Beyond state assessment, ROC curves can evaluate trends, such as a metric based on the x-over-y rule of a linear trend from recent data points.
2-over-3 Rule
ggplot(model.frame(FLQuants(State =state,
Indicator=trnd)))+
geom_point(aes(State,Indicator))+
geom_hline(aes(yintercept=1),linetype=2)+
geom_vline(aes(xintercept=1),linetype=2)+
geom_text(aes(x=x,y=y,label=label),data=data.frame(x=c(0.1,0.1,2,2),y=c(0.1,2.75,2.75,0.1),label=c("TN","FP","TP","FN")),
size=6,col="blue")+
xlab("Operating Model Status")+ylab("Indicator/Reference")+
theme_bw()
Figure 10. Pairwise scatter plot of index relative to reference verses \(SSB/B_{MSY}\)
ggplot(idxSkill)+
geom_line(aes(FPR,TPR))+
geom_abline(aes(intercept=0,slope=1),linetype=2)+
geom_point(aes(FPR,TPR),data=idxSkill[order((idxSkill$ind-1)^2)[1],],size=5)+
geom_point(aes(FPR,TPR),data=idxSkill[max(idxSkill$TSS)==idxSkill$TSS,],size=5,col="blue")+
xlab("False Positive Rate (FPR)")+ylab("True Positive Rate (TPR)")
Figure 11. ROC curve, black point is when index/reference =1, and red is the point with the best skill score.
ggplot(idxSkill)+
geom_line(aes(ind,TSS))+
geom_vline(aes(xintercept=1),col="red")+
xlab("Reference Level")+ylab("True Skill Score (TSS)")
Figure 12. Plot of TSS v index value.
Reference Period
ggplot(model.frame(FLQuants(State =state,
Indicator=trnd)))+
geom_point(aes(State,Indicator))+
geom_hline(aes(yintercept=1),linetype=2)+
geom_vline(aes(xintercept=1),linetype=2)+
geom_text(aes(x=x,y=y,label=label),data=data.frame(x=c(0.1,0.1,2,2),y=c(0.1,2.75,2.75,0.1),label=c("TN","FP","TP","FN")),
size=6,col="blue")+
xlab("Operating Model Status")+ylab("Indicator/Reference")+
theme_bw()
Figure 13. Pairwise scatter plot of index relative to reference verses \(SSB/B_{MSY}\)
ggplot(idxSkill)+
geom_line(aes(FPR,TPR))+
geom_abline(aes(intercept=0,slope=1),linetype=2)+
geom_point(aes(FPR,TPR),data=idxSkill[order((idxSkill$ind-1)^2)[1],],size=5)+
geom_point(aes(FPR,TPR),data=idxSkill[max(idxSkill$TSS)==idxSkill$TSS,],size=5,col="blue")+
xlab("False Positive Rate (FPR)")+ylab("True Positive Rate (TPR)")
Figure 14. ROC curve, black point is when index/reference =1, and red is the point with the best skill score.
ggplot(idxSkill)+
geom_line(aes(ind,TSS))+
geom_vline(aes(xintercept=1),col="red")+
xlab("Reference Level")+ylab("True Skill Score (TSS)") 
Figure 15. Plot of TSS v index value.
Linear Regression
trnd<-function(object,n=5){
m =expand.grid(iter=dimnames(object)$iter,year=dimnames(object)$year[-seq(n-1)])
rtn=mdply(m, function(iter,year){
dat=as.data.frame(object[,rev(ac(an(ac(year))+1-seq(n))),,,,iter])
data.frame(data=lm(data~year,data=dat)$coefficients[2])})
as.FLQuant(rtn)}
state=trnd(ssb(om)[,ac(21:40)]/mean(ssb(om)[,ac(21:40)]))
trnd =trnd(idx[,ac(21:40)]/mean(idx[,ac(21:40)]))
idxSkill=FLCore::roc(state>0,trnd) More Information
- For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage 1.
Software Versions
- R version 4.2.1 (2022-06-23)
- FLCore: 2.6.19.9105
- FLPKG:
- Compiled: Sun Jan 21 12:38:52 2024
- Git Hash: 7bdba0c