Assessment with FLIBM objects
1 Setup
The following packages are used in the tutorial:
# required packages
library(FLIBM)
library(FLa4a)
library(FLBRP)
library(spict)
library(TropFishR)
library(MLZ)2 Introduction
The examples below all start with a synthetic stock simulated with
FLIBM (stk1).
data(stk1) # load
stk <- stk1 # rename
plot(stk$stock.a@stock.n) # plot stock numbers
3 Age-based (a4a)
Collapse to yearly FLStock
As a first step, we simplify the object to yearly basis and trim the years to include those with some noticible fishing (and catches). Several
# collapse season dimension of age-based FLStock
stkTrue <- simplifySeason(stk)
# exclude zero age class & redefine plusgroup to collapse rare ages
stkTrue <- stkTrue[ac(1:range(stkTrue)["max"]),]
stkTrue <- setPlusGroup(stkTrue, plusgroup = 4, na.rm=TRUE)
# create time series of SSB, and compare to yearly estimates
tmp2 <- as.data.frame(apply((stkTrue@stock.n * stkTrue@stock.wt * stkTrue@mat), 2, sum, na.rm = TRUE)) #as.data.frame(ssb(stkTrue))
# spwn_wt <- weighted.mean(0:11/12, stk$rec$params$season_wt)
tmp2$date <- yeardec2date(tmp2$year)#+spwn_wt)
tmp3 <- as.data.frame(apply((stk$stock.a@stock.n * stk$stock.a@stock.wt * stk$stock.a@mat),
MARGIN = c(2:6), sum, na.rm=TRUE))
tmp3$date <- yeardec2date(tmp3$year + (an(tmp3$season)-1)/12)
ggplot(data = tmp3) + aes(x = date, y = data) + geom_line() +
geom_point(mapping = aes(x=date, y=data), data = tmp2)
# fill in missing slots
stkTrue@catch.n[] <- replace(stkTrue@catch.n, is.na(stkTrue@catch.n), 0)
stkTrue@landings.n <- stkTrue@catch.n
stkTrue@landings.wt <- stkTrue@catch.wt
stkTrue@discards.n[] <- 0
stkTrue@discards.wt[] <- 0
stock(stkTrue) <- computeStock(stkTrue)
landings(stkTrue) <- computeLandings(stkTrue)
catch(stkTrue) <- computeCatch(stkTrue)
discards(stkTrue) <- computeDiscards(stkTrue)
# plot yearly FLStock
plot(stkTrue)
Create survey indicies
# Create Index
set.seed(1)
idx <- FLIndex(index=stock.n(stkTrue)/1000, type="numbers")
idx@catch.n <- idx@index
idx@effort[] <- 1
range(idx)[c("startf","endf")] <- c(0,0.01)
# add observation error
index(idx) <- index(idx)*exp(rnorm(1, index(idx)-index(idx), 0.1))
plot(idx)
3.1 Assessment For All (a4a)
Adding smoothers to assessment help to reduce degrees of freedom. In the following example, spline terms have been used in the submodels for recruitment and yearly fishing mortality patterns. Age fishing mortality patterns are left free, allowing for individual F values for each age. Since there are only 3 ages in the object, this is does not represent a large number of coefficients, however stocks with more ages might benefit from an additional spline for catchability at age, as one would not expect large deviations between adjacent ages. In fact, here we see very little deviation in the age pattern.
fmod <- formula( ~ factor(age) + s(year, k=10) ) # fishing mortality submodel
srmod <- formula( ~ s(year, k=6)) # recruitment submodel
qmod <- list( ~ factor(age)) # survey catchability submodel
tmp <- sca(stkTrue, idx, srmodel = srmod, fmodel = fmod, qmodel = qmod)
stk_a4a <- stkTrue + tmp
# plot
plot(FLStocks(true=stkTrue, a4a=stk_a4a))
wireframe(harvest(stk_a4a), zlab="F")
3.2 Reference points
Stock-recruitment relationship
srbh <- fmle(as.FLSR(stk_a4a, model="bevholt"), method="L-BFGS-B",
lower=c(1e-6, 1e-6), upper=c(max(rec(stk_a4a)) * 3, Inf))
plot(srbh) 
Calculate reference points
# Calculate the reference points
brp <- brp(FLBRP(stk_a4a, srbh))
Fmsy <- c(refpts(brp)["msy","harvest"])
msy <- c(refpts(brp)["msy","yield"])
Bmsy <- c(refpts(brp)["msy","ssb"])
Bpa <- 0.5*Bmsy
Blim <- Bpa/1.4
plot(brp)
4 Biomass-based (SPiCT)
## Catches
Caa <- apply(stk$stock.a@catch.n * stk$stock.a@catch.wt,
1:2, sum, na.rm=TRUE)
Cyear <- as.data.frame(apply(Caa, 2, sum))
## Biomass indices
# Baa <- stk$stock.a@stock.n * stk$stock.a@stock.wt
# # sum exploited part (ages > 1)
# Bseas <- apply(Baa[2:7,],c(2,4), sum, na.rm=TRUE)
# Bseas <- as.data.frame(Bseas)
# Bseas$season <- as.numeric(as.character(Bseas$season))
# Bseas$yeardec <- Bseas$year +
# ((Bseas$season-1)/max(Bseas$season))
# Bseas <- Bseas[order(Bseas$yeardec),]
# Bseas$yeardec <- round(Bseas$yeardec, 3)
# scaled fishing mortality
Fsc <- apply(stk$stock.a@harvest, MARGIN = c(2:6),
FUN = function(x){x/max(x, na.rm = TRUE)})
tmp <- apply(stk$stock.a@harvest, MARGIN = c(2,4),
FUN = function(x){x/max(x, na.rm = TRUE)})
tmp <- apply(tmp, MARGIN = 1, sum, na.rm = TRUE)
tmp <- tmp/max(tmp)
tmp2 <- stk$stock.a@stock.n * stk$stock.a@stock.wt * c(tmp)
Bseas <- as.data.frame(apply(tmp2, c(2,4), sum, na.rm = T))
Bseas$season <- as.numeric(as.character(Bseas$season))
Bseas$yeardec <- Bseas$year +
((Bseas$season-1)/max(Bseas$season))
Bseas <- Bseas[order(Bseas$yeardec),]
Bseas$yeardec <- round(Bseas$yeardec, 3)
# start of year
idx1 <- subset(Bseas, season == 1)
# mid-year (only for second half of time-series)
idx2 <- subset(Bseas, season == 7 & year %in% 1990:1999)
## build input
inp <- list()
# catches
inp$timeC <- Cyear$year
inp$obsC <- Cyear$data
inp$dtc <- 1
# indices
inp$timeI <- list(idx1$yeardec, idx2$yeardec)
inp$obsI <- list(idx1$data, idx2$data)
inp$dteuler <- 1/12
# necessary to fix catchability, to indicate that index
# represents absolute biomass
inp$ini$logq <- log(c(1,1))
inp$phases$logq <- -1
# inp$dtpredc <- 1
# check for errors
inp <- check.inp(inp)
# fit SPiCT and plot summary
fit <- fit.spict(inp)
# plot summary
op <- par(mfcol=c(2,2), mar = c(3,3,2,4), mgp = c(2,0.5,0), ps=8)
plotspict.biomass(fit)
plotspict.f(fit)
plotspict.production(fit)
plotspict.fb(fit)
par(op)Compare biomass:
tmp <- as.data.frame(get.par("logB", fit, exp = TRUE))
tmp$yeardec <- round(as.numeric(rownames(tmp)),3)
tmp <- merge(x = Bseas, y = tmp, all.x=TRUE)
op <- par(mar = c(4,4,1,4))
plot(data ~ yeardec, tmp, t="l", ylim = c(0, max(tmp$data)))
grid()
text(max(tmp$yeardec), 0.95*max(tmp$data),
labels = paste("r.sq =",
round(cor(tmp$est, tmp$data)^2,4)),
pos = 2)
par(new = TRUE)
plot(est ~ yeardec, tmp, t="l", ylim = c(0, max(tmp$est)),
col=2, lty=2, axes=FALSE, xlab="", ylab="")
axis(4, col=2, col.axis=2); mtext("est", col=2, side=4,
line=3)
par(op)5 Length-based (LBI, MLZ, Thompson-Bell)
Fit growth parameters
lfq <- flquant2lfq(stk$stock.l@catch.n)
pal <- colorRampPalette(c("grey30",5,7,2), bias=2)
with(lfq, image(x=dates, y=midLengths, z=t(catch), col=pal(100)))
lfq2 <- TropFishR::lfqModify(lfq, bin_size = 2, years = 1998:1999)
lfq2 <- TropFishR::lfqRestructure(lfq2)
lfq2 <- ELEFAN_GA(lfq2, seasonalised = FALSE,
low_par = list(Linf = 60, K = 0.1, ta = 0),
up_par = list(Linf = 100, K = 2, ta = 1),
popSize = 80, maxiter = 100, pmutation = 0.2,
run = 20, monitor = FALSE,
MA = 5, parallel = TRUE, seed = 1)Genetic algorithm is running. This might take some time.
unlist(lfq2$par) Linf K t_anchor phiL
75.6063028 0.5155421 0.1982915 3.4693801 plot(lfq2, Fname = "rcounts")
5.1 MLZ
# make MLZ class
lfq3 <- TropFishR::lfqModify(lfq, aggregate = "year", bin_size = 2)
dim(lfq3$catch); length(lfq3$midLengths); length(lfq3$dates)
lfq3$par <- lfq2$par
plot(lfq3, Fname = "catch")
new.dataset <- new("MLZ_data",
Year = as.numeric(format(lfq3$dates, "%Y")),
Len_bins = lfq3$midLengths,
Len_matrix = t(lfq3$catch),
length.units = "cm")
new.dataset@Lc <- stk$harvest$params$L50
new.dataset <- calc_ML(new.dataset)
# convert t_anchor to t0
t_anchor2t0 <- function(par = NULL, Lrecr = NULL){
if(is.null(par) | is.null(Lrecr)){
stop("Error: must provide both 'par' and 'Lrecr'")
}
if(Lrecr > par$Linf){stop("Error: 'Lrecr' must be lower than 'par$Linf'")}
# add seasonalized pars if needed
seasonalized <- !is.null(par$C)
if(!seasonalized){par$ts <- 0; par$C <- 0}
t <- seq(0, 3, 0.01)
Lt <- VBGF(param = par, t = t)
trecr <- t[min(which((Lt - Lrecr)>0))]
# calc_trecr(par = par, Lrecr = Lrecr, Lt = tail(Lt,1), t = tail(t,1)) # same
t0 <- par$t_anchor - trecr
if(seasonalized){
ts <- (par$ts - trecr)%%1
} else {
ts <- 0
}
par_age <- par
par_age$t_anchor <- NULL
par_age$t0 <- t0
par_age$ts <- ts
return(par_age)
}
# new.dataset@MeanLength
# new.dataset@ss
new.dataset@vbLinf <- c(lfq3$par$Linf)
new.dataset@vbK <- lfq3$par$K[1]
new.dataset@vbt0 <- t_anchor2t0(par = lfq3$par, Lrecr = 0)$t0
plot(new.dataset)
model1 <- ML(new.dataset, ncp = 0, figure = FALSE)
model2 <- ML(new.dataset, ncp = 1, figure = FALSE)
model3 <- ML(new.dataset, ncp = 2, figure = FALSE)
compare_models(list(model1, model2, model3))
summary(model3)
# With Effort
yrs <- 1960:1999
steepness <- 0.75
FMmax <- 1
FMs <- FMmax / (1 + exp(-steepness * (yrs - 1985) ))
new.dataset@Effort <- FMs[which(yrs >= 1980)]
new.dataset@M <- 0.7
pred <- data.frame(Year = yrs[yrs >= 1980])
breaks <- c(-Inf,
model3@estimates[(model3@n.changepoint+2):
(model3@n.changepoint+2+model3@n.changepoint-1),1],
Inf)
Zs <- model3@estimates[seq(model3@n.changepoint+1),1]
pred$Z <- Zs[cut(pred$Year, breaks = breaks)]
plot(yrs, FMs+0.7, t = "o", ylim = c(0,max(c(pred$Z, FMs+0.7))))
lines(Z ~ Year, pred, t="o", col=2)
legend("topleft", legend = c("True", "MLZ est."), col=1:2, pch=1, lty = 1)
# Mean length with effort mortality estimator (not working yet)
# MLeffort(new.dataset, start = list(q = 0.1, M = 0.7), n_age = 6, estimate.M = FALSE)VPA - variable fishing mortality (F)
TO COME
** Yield-per-recruit (YPR) and estimation of Fmax, F01, and SPR**
6 References
L. T. Kell, I. Mosqueira, P. Grosjean, J-M. Fromentin, D. Garcia, R. Hillary, E. Jardim, S. Mardle, M. A. Pastoors, J. J. Poos, F. Scott, R. D. Scott; FLR: an open-source framework for the evaluation and development of management strategies. ICES J Mar Sci 2007; 64 (4): 640-646. doi: 10.1093/icesjms/fsm012.
7 More information
- You can submit bug reports, questions or suggestions on this tutorial at https://github.com/marchtaylor/FLIBM/issues.
- Alternatively, send a pull request to https://github.com/marchtaylor/FLIBM/.
7.1 Software Versions
- R version 4.4.2 (2024-10-31 ucrt)
- data.table: 1.16.2
- FLa4a: 1.8.3.9001
- FLAssess: 2.6.3
- FLBRP: 2.5.9.9025
- FLCore: 2.6.20.9102
- FLXSA: 2.6.6.9001
- FLIBM: 0.3.3
- ggplotFL: 2.7.0.9133
- spict: 1.3.8
- TropFishR: 1.6.5
- Compiled: 2025-Jan-16
7.2 License
This document is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license.