Fisheries Management: Indicator Selection and ROC Analysis

source("~/Desktop/flr/FLCandy/R/saveTmp.R", echo=TRUE)

## 
## > sv <- function(..., list = character(), file = "~/Desktop/tmp/t.RData", 
## +     ascii = FALSE, version = NULL, envir = parent.frame(), compress = isT .... [TRUNCATED] 
## 
## > ld <- function(file = "~/Desktop/tmp/t.RData", envir = parent.frame(), 
## +     verbose = FALSE) {
## +     if (is.character(file)) {
## +         con <- gz .... [TRUNCATED]

Jump to Introduction

Back to Top

a

• Identify and prioritise Objectivesand trade-offs
• Selection of Hypotheses
• Simulate Data
• Conditioning the Operating Model based on data and knowledge
• Identifying candidate management strategies
• Running the Management Procedure as a feedback controller to simulate the long-term impact of management
• Selection of the Management Strategy that best meets the Objectives

Jump to Packages

setwd("~/Desktop/active/fao/Rmd/")

FLR packages

library(FLCore)
library(FLBRP)
library(FLasher)
library(FLife)
library(ggplotFL)

library(FLCandy)

Packages for extracting life history parameters from FishLife

library(rfishbase)
library(SPMpriors)
library(FishLife)

Utility packages, for data wrangling

library(plyr)
library(dplyr)
library(reshape)

Plotting packages

library(ggplot2)
theme_set(theme_bw())

library(ggpubr)
library(ggcorrplot)
library(GGally)

Installation

Libraries (i.e. packages) can be installed from cran or via gihtub using the remotes package for the FLR, FishBase and FishLife packages

install.packages("remotes")

library(remotes)

remotes::install_github("flr/FLCore")

There are compatability problems with FishBase, so install an earlier version

remotes::install_github("ropensci/rfishbase", ref = "d150f2e0f5")

remotes::install_github("james-thorson/FishLife")

Back to Top

Objectives

Evaluation

Compare performance against management objectives

• Safety: Low probability of stock collapse
• Status: Healthy stock most years
• Yield: High yields
• Variability: no big changes between years, seasons, …

Back to Top

Scenarios

Selection of Hypotheses

• Lake Tanganyika sprat (S. tanganicae) spawns several times a year, annual mortality is high and so it has a short life span.
• As juveniles grow, they move progressively closer to shore before moving offshore and start to appear in pelagic catches.
• Therefore, an operating model with a seasonal structure with multiple spawning episodes per year was developed.
• Each spawning episode produces a “subcohort” which can have different growth, maturity, and natural mortality.

Read in life history parameters from Fishbase using the FishLife package

pars = flmvn_traits(Genus = "Stolothrissa", Species = "tanganicae", Plot = FALSE)

     [,1]           [,2]         
[1,] "K"            "M"          
[2,] "Winfinity"    "Loo"        
[3,] "tmax"         "tm"         
[4,] "Lm"           "Temperature"
[5,] "ln_margsd"    "rho"        
[6,] "logitbound_h" "ln_r"       

Extract the life history parameters required to parameterise the Operating Model, these are parameters for the individual growth (\(k\) and \(L_{\infty}\)), age at maturity (a_{50}\(, natural mortality (\)M\(), the steepness of the stock recruitment relationship (\)h$), and recruitment variability (sigR).

par = FLPar(c(cast(pars[[1]][, c("trait", "mu.sp")], ~trait)[, c("K", "Loo", "tm",
    "M", "h", "sigR")]))

Rename parameters, so they are compatible with the FLR package naming convention, and then save the expected values, and the covariances

dimnames(par)[[1]] = c("k", "linf", "a50", "m", "s", "sigmaR")
cov = pars[[2]][[1]][[1]][c("K", "Loo", "tm", "M", "h"), c("K", "Loo", "tm", "M",
    "h")]
dimnames(cov) = list(c("k", "linf", "a50", "m", "s"), c("k", "linf", "a50", "m",
    "s"))
S.tanganicae = list(param = par, cov = cov)

par = S.tanganicae[[1]]
par = rbind(lhPar(rbind(par[c("k", "linf", "a50", "s")], FLPar(t0 = -0.1))), par["m"])

par["a50"] = 1

par

An object of class "FLPar"
params
     linf         k        t0         a         b       l50       a50     ato95 
   9.9650    2.1581   -0.1000    0.0003    3.0000    7.8310    1.0000    1.0000 
     asym        bg        m1        m2        m3         s         v      sel1 
   1.0000    3.0000    0.5500   -1.6100    1.4400    0.9141 1000.0000    1.6141 
     sel2      sel3         m 
   1.0000 5000.0000    3.3364 
units:  NA 

Natural mortality rates (\(M\)) vary with body size and age, therefore, take the single value from FishBase and use the ‘generalized length-inverse mortality’ (GLIM, Lorenzen 2022) paradigm to model M-at-age, i.e.

\(log(M)= m_1 + m_2 \times log(L/L_{50}_{\infty}) + m_3 \times log(k)\)

mGlim <- function(object, par) {
    # M= exp(a + b ln(L/L∞) + c*lnK)
    exp(par["m1"] %+% (par["m2"] %*% log(object%/%par["linf"])) %+% (par["m3"] %*%
        log(par["k"])))
}

mGlimFn = function(object, par) {
    # lnM= a + b ln(L/L∞) + c*lnK
    len = wt2len(stock.wt(object), par)

    exp(par["m1"] %+% (par["m2"] %*% log(len%/%par["linf"])) %+% (par["m3"] %*% log(par["k"])))
}


par["m2"] = -1.3
par["m3"] = 1.08

Scale m1 so that \(M\) at \(l_{50}\) is equal to the ‘FishBase’ value

$ m_1 = log(M) - m_2log(L_{50}/L_{}) - m_3log(k)$

par["m1"] = log(par["m"]) %-% (par["m2"] %*% log(par["l50"]%/%par["linf"])) %-% (par["m3"] %*%
    log(par["k"]))

par

An object of class "FLPar"
params
     linf         k        t0         a         b       l50       a50     ato95 
   9.9650    2.1581   -0.1000    0.0003    3.0000    7.8310    1.0000    1.0000 
     asym        bg        m1        m2        m3         s         v      sel1 
   1.0000    3.0000    0.0608   -1.3000    1.0800    0.9141 1000.0000    1.6141 
     sel2      sel3         m 
   1.0000 5000.0000    3.3364 
units:  NA 

ages = FLQuant(0:10, dimnames = list(age = 0:10))
length = vonB(ages, par)

ggplot(as.data.frame(FLQuants(length = vonB(ages, par), weight = len2wt(length, par),
    mat = logistic(length, par), m = mGlim(length, par)))) + geom_line(aes(age, data)) +
    facet_wrap(~qname, scale = "free", ncol = 2)

Figure 1 Vectors of quantities-at-age.

ggplot(as.data.frame(FLQuants(mat = logistic(length, par), Selectivity = dnormal(ages,
    par)))) + geom_line(aes(age, data)) + facet_wrap(~qname, scale = "free", ncol = 2)

Figure 2 Selectivity and maturity-at-age.

Back to Top

Conditioning the Operating Model

sv(par)
ld()

par = lhPar(par[-19])
eq = lhEql(par, m = mGlimFn)
fbar(eq)[] = mean(m(eq)[4])
om = as(brp(eq), "FLStock")
om = fwd(om, f = fbar(om)[, -1], sr = eq)

trnd = FLQuant(c(t(maply(c(seq(5, 1, length.out = 51), seq(1, 0.2, length.out = 51)[-1]),
    function(x) seq(1, x, length.out = 21)))), dimnames = list(year = seq(21), iter = seq(101)))

kick = FLQuant(rep(c(trnd[, 21]), each = 21), dimnames = list(year = seq(21), iter = seq(101)))

oms = FLStocks(Constant = om)
oms[["Trend"]] = fwd(om, f = fbar(om)[, ac(51:71)] %*% trnd, sr = eq)
oms[["Kick"]] = fwd(om, f = fbar(om)[, ac(51:71)] %*% kick, sr = eq)
oms = FLStocks(llply(oms, window, start = 31, end = 71))

[1] "maxfbar has been changed to accomodate new plusgroup"
[1] "maxfbar has been changed to accomodate new plusgroup"
[1] "maxfbar has been changed to accomodate new plusgroup"

p1 = plot(window(oms, start = 51), metrics = list(invF = function(x) 1/fbar(x), Index = vb,
    F = fbar)) + theme_bw() + theme(legend.position = "bottom") + scale_fill_manual("Scenario",
    values = rainbow(3)) + scale_colour_manual("Scenario", values = rainbow(3)) +
    facet_grid(qname ~ stock, scale = "free_y")
p1

Figure 3 Operating Model, for the three trends in fishing mortality.

Back to Top

Generating Data

sv(par, oms)
ld()

ak = FLCore:::invALK(par, cv = 0.3, age = an(dimnames(oms[[1]])$age), bin = 0.1)
lfd = mydas:::lenSample(iter(catch.n(oms[["Trend"]]), 1), ak, nsample = 100)

dat = subset(transform(as.data.frame(lfd), lustrum = 5 * (year%/%5)), data > 0)

dt2 = ddply(dat, .(lustrum), with, {
    lmodeFn <- function(len, n, bins = 25) {

        dat = data.frame(bin = cut(len, breaks = bins), n = n)
        res = ddply(dat, .(bin), with, data.frame(freq = sum(n)))
        res = as.character(subset(res, freq == max(freq))[1, "bin"])
        unbin(res)$mid
    }

    unbin = function(x) {
        left = as.numeric(substr(x, 2, unlist(gregexpr(",", x)) - 1))
        right = as.numeric(substr(x, unlist(gregexpr(",", x)) + 1, nchar(x) - 1))
        mid = (left + right)/2

        data.frame(left = left, right = right, mid = mid, n = x)
    }

    lc = lmodeFn(len, data) * 0.5
    data.frame(lc = lc, lmean = weighted.mean(len, data * (len >= lc)))
})

gghistogram(dat, x = "len", weight = "data", bins = 30) + coord_flip() + scale_x_reverse() +
    facet_grid(. ~ lustrum, scale = "free") + xlab("Length (cm)") + geom_vline(aes(xintercept = lc),
    data = dt2, col = "red") + geom_vline(aes(xintercept = lmean), data = dt2, col = "blue") +
    scale_x_continuous(limits = c(0, 6)) + theme_bw(16) + theme(legend.position = "none",
    axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())

Figure 4 Length data

p2 = plot(window(oms, start = 51), metrics = list(LBI = function(x) {
    res = indicators.len(x[-1], indicators = c("lmean"), params = apply(par, 1, median),
        metric = "stock.n")[[1]]
    names(dimnames(res))[1] = "age"
    res
}, Index = function(x) rlnorm(dim(x)[6], log(vb(x)), 0.3), HR = function(x) apply(catch(x),
    c(2.6), sum)%/%vb(x))) + theme_bw() + theme(legend.position = "bottom") + facet_grid(qname ~
    ., scale = "free_y") + scale_fill_manual("Scenario", values = rainbow(3)) + scale_colour_manual("Scenario",
    values = rainbow(3)) + scale_y_continuous(limit = c(9, NA))
p2 + facet_grid(qname ~ stock, scale = "free_y")

Figure 5 Indices generated by the Observation Error Model, for the three trends in fishing mortality.

ts = ldply(oms, function(x) {
    rbind(cbind(quant = "hr", melt(apply(h(fbar(x)), c(2, 6), sum)[drop = T])), cbind(quant = "catch",
        melt(apply(catch(x), c(2, 6), sum)[drop = T])), cbind(quant = "biomass",
        melt(biomass(x)[drop = T])), cbind(quant = "ssb", melt(ssb(x)[drop = T])),
        cbind(quant = "vb", melt(vb(x)[, drop = T])), cbind(quant = "idx", melt((devH %*%
            biomass(x))[drop = T])), cbind(quant = "nidx", melt(quantSums(lenSamp(x,
            params = FLPar(apply(par, 1, median)), bin = 0.2, metric = "stock.n")[-seq(25)])[drop = T])),
        cbind(quant = "hr.hat", melt((catch(x)%/%(devH %*% biomass(x)))[drop = T])),
        cbind(quant = "lbi", melt(indicators.len(x, indicators = c("lbar"), params = FLPar(apply(par,
            1, median)), metric = "catch.n")[[1]][drop = T])))
})

ts = cast(ts, .id + year + iter ~ quant)
ts = transmute(merge(ts, subset(ts, year == 1990)[, -2], by = c(".id", "iter")),
    .id = .id, year = year, iter = iter, biomass = biomass.x/biomass.y, ssb = ssb.x/ssb.y,
    vb = vb.x/vb.y, catch = catch.x/catch.y, hr = hr.x/hr.y, hr.hat = hr.hat.x/hr.hat.y,
    idx = idx.x/idx.y, nidx = nidx.x/nidx.y, lbi = lbi.x/lbi.y)

dat = ddply(subset(ts, .id != "Constant" & year %in% 2001:2020), .(.id), transform,
    x = (ssb - mean(ssb))/var(ssb)^0.5, y = (nidx - mean(idx))/var(nidx)^0.5)
p1 = ggplot(dat) + geom_point(aes(x, y), size = 0.2) + geom_vline(aes(xintercept = 0),
    col = "blue") + geom_hline(aes(yintercept = 0), col = "blue") + facet_grid(. ~
    .id, scale = "free") + xlab("OM") + ylab("OEM") + theme_bw(20)

tss = ddply(dat, .(.id), with, data.frame(hr = mean((x > 0 & y > 0)) - mean(x < 0 &
    y > 0)))

dat = ddply(subset(ts, .id != "Constant" & year %in% 2001:2020), .(.id), transform,
    x = (hr - mean(hr))/var(hr)^0.5, y = (hr.hat - mean(hr.hat))/var(hr.hat)^0.5)
p2 = ggplot(dat) + geom_point(aes(x, y), size = 0.2) + geom_vline(aes(xintercept = 0),
    col = "blue") + geom_hline(aes(yintercept = 0), col = "blue") + facet_grid(. ~
    .id, scale = "free") + xlab("OM") + ylab("OEM") + theme_bw(20)

tss = ddply(dat, .(.id), with, data.frame(hr = mean((x > 0 & y > 0)) - mean(x < 0 &
    y > 0)))

dat = ddply(subset(ts, .id != "Constant" & year %in% 2001:2020), .(.id), transform,
    x = (hr - mean(hr))/var(hr)^0.5, y = (1/lbi - mean(1/lbi))/var(1/lbi)^0.5)
p3 = ggplot(dat) + geom_point(aes(x, y), size = 0.2) + geom_vline(aes(xintercept = 0),
    col = "blue") + geom_hline(aes(yintercept = 0), col = "blue") + facet_grid(. ~
    .id, scale = "free") + xlab("OM") + ylab("OEM") + theme_bw(20)

tss = ddply(dat, .(.id), with, data.frame(hr = mean((x > 0 & y > 0)) - mean(x < 0 &
    y > 0)))

ggarrange(p1, p2, p3, label.y = 1, widths = c(1, 1, 1), heights = c(6), nrow = 3,
    ncol = 1, labels = c("Index", "SBI", "Harvest Rate"), common.legend = TRUE, align = "hv") +
    theme_bw(12) + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
    axis.title.y = element_blank())

Figure 6 Comparison of Operating Model (OM) and Indicators (OEM).

Back to Top

Management Procedure

Back to Top

Running

Back to Top

More Information

Back to Top

References

Back to Top

Lorenzen, Kai. 2022. “Size-and Age-Dependent Natural Mortality in Fish Populations: Biology, Models, Implications, and a Generalized Length-Inverse Mortality Paradigm.” Fisheries Research 255: 106454.