FLIBM is a individual-based model (IBM) framework for the simulation
of fish or invertibrate populations, and should provide a flexible
structure that allows for modification given the population under study.
Presently, the IBM describes processes affecting life history in terms
of length (e.g. growth, mortality, maturity). This faciltates the
generation of length-frequency data (FLD), which is typically used in
many data-poor assessment approaches. Nevertheless, an individual’s age
is also followed in order to be able to record both length- and
age-based summary statistics in the form of FLStock
objects
(see FLR project: http://www.flr-project.org/). Following the FLR
framework allows for integration in other FLR-related recources.
The main sequence of processes in the IBM operating model are as follows:
Generally, the model should be run at a sufficiently high number of interannual time steps (default = 12, i.e. monthly) to produce realistic behaviour and summary indices.
The FLIBM imports the package FLCore
and recommends
ggplotFL
for use of summary plotting functions.
library(FLIBM)
library(FLCore)
library(ggplotFL)
library(data.table)
set.seed(42) # for reproducibility
The create.FLIBM
function is used to intialize an FLIBM
object. As with FLStock
and FLQuant
objects,
one defines 6 data dimensions:
Other arguments for create.FLIBM
are related to units of
numbers (n.units
) and weights (wt.units
).
stk <- create.FLIBM(
length = 0:85, age = 0:6,
year = ac(2000:2009),
season = ac(1:12),
n.units = "1e3", wt.units = "kg"
)
The basic structure of an FLIBM class is as follows:
summary(stk)
Length Class Mode
growth 2 -none- list
rec 5 -none- list
m 1 -none- list
harvest 2 -none- list
n.units 1 -none- character
wt.units 1 -none- character
make.inds 1 -none- function
inds 14 data.table list
stock.l 1 FLStock S4
stock.a 1 FLStock S4
length.a 840 FLQuant numeric
age.l 10320 FLQuant numeric
There are three main types of objects within the class:
FLStock
objects
obj$stock.l
and obj$stock.a
for length- and
age-based IBM summary statistics, respectively, which are recorded at
the start each time step.obj$growth
,
obj$rec
, obj$m
,
obj$harvest
).obj$inds
). It is a multi-level list class with
the dimensions “unit”, “area” and “iter”, following the framework of
FLStock objects
(e.g. obj$inds[[unit]][[area]][[iter]]
)The main objects types that are manipulated by the user are the controlling objects.
stk$rec$params['rmax'] <- 1e4
By default, FLStock objects are set up to have fbar ranges spanning
the full range of dimension 1 (age
and
length
).
range(stk$stock.a); range(stk$stock.l)
min max plusgroup minyear maxyear minfbar maxfbar
0 6 6 2000 2009 0 6
min max plusgroup minyear maxyear minfbar maxfbar
0 85 85 2000 2009 0 85
For plotting and further FLR-related functionality, it is advised to adjust these as appropriate, e.g.:
range(stk$stock.a)[c("minfbar", "maxfbar")] <- c(1,3)
range(stk$stock.a)
min max plusgroup minyear maxyear minfbar maxfbar
0 6 6 2000 2009 1 3
A spin-up of the ìnds`is also possible. At the moment this function uses the control information from the first simulation year.
stk <- spinup.FLIBM(obj = stk, nyearsslope = 10)
year = 1 | ssb = 2.072 | trend = 1
year = 2 | ssb = 123.339 | trend = 1
year = 3 | ssb = 851.973 | trend = 1
year = 4 | ssb = 1667.527 | trend = 1
year = 5 | ssb = 1927.63 | trend = 1
year = 6 | ssb = 1917.891 | trend = 1
year = 7 | ssb = 2061.884 | trend = 1
year = 8 | ssb = 2129.647 | trend = 1
year = 9 | ssb = 2086.101 | trend = 1
year = 10 | ssb = 2084.67 | trend = 1
year = 11 | ssb = 2036.07 | trend = 1
year = 12 | ssb = 2055.573 | trend = 1
year = 13 | ssb = 2038.006 | trend = 1
year = 14 | ssb = 2107.572 | trend = 0
Viewing the $inds
ggplot(stk$inds, aes(x = length)) +
theme(text = element_text(size=12)) +
geom_histogram(color = "black", fill = "white",
na.rm=TRUE)
stk <- adv.FLIBM(stk, years = ac(2000:2003), monitor = FALSE)
# plot ssb
df <- as.data.frame(stock.n(stk$stock.a)
* stock.wt(stk$stock.a) * mat(stk$stock.a))
df$date <- as.Date(paste(as.numeric(df$year),
as.numeric(as.character(df$season)), "01", sep = "-"),
format = "%Y-%m-%d")
df <- aggregate(x = df$data, by = list(date = df$date),
FUN = sum, na.rm=TRUE)
df$x[df$x == 0] <- NA
p <- ggplot(data = df, aes(x = date, y = x)) +
geom_line() +
scale_y_continuous(limits = c(0,NA)) +
labs(y = "SSB")
print(p)
# plot stock numbers
plot(stk$stock.a@stock.n)
Summary stats
plot(simplifySeason(stk))
# plot length-frequency
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)))
df <- data.frame(
age = as.numeric(dimnames(stk$stock.a@catch.n)$age),
n = c(apply(window(stk$stock.a@catch.n, start=2003, end=2003),
1, sum, na.rm=TRUE)))
plot(log(n) ~ age, df)
fit <- lm(log(n) ~ age, df, subset = age %in% 2:5)
points(log(n) ~ age, df, subset = age %in% 2:5, pch=20)
abline(fit, col=4)
usr <- par()$usr
text(0.9*usr[2], 0.9*usr[4],
labels = bquote(italic(Z)==.(sprintf("%.2f",-round(coef(fit)[2], 2)))),
pos = 2, col=4)
To come…
length at maturity
df <- as.data.frame(apply(stk$stock.l@mat[,"2003"],
1, mean, na.rm=TRUE))
names(df)[1] <- "length"
df <- na.omit(df)
df$length <- df$length+0.5 # bin mid-length
fit <- glm(data ~ length, data = df,
family = binomial(link = "logit"))
pred <- predict(fit, se.fit=TRUE)
df <- cbind(df, pred)
# Calculate confidence intervals
std <- qnorm(0.95 / 2 + 0.5)
df$ymin <- fit$family$linkinv(df$fit - std * df$se.fit)
df$ymax <- fit$family$linkinv(df$fit + std * df$se.fit)
df$fit <- fit$family$linkinv(df$fit) # Rescale to 0-1
plot(data ~ length, data = df, t="n", ylab="proportion mature")
points(data ~ length, data = df, col=adjustcolor(1, 0.5))
lines(fit ~ length, df, col=1)
lrPerc <- function(alpha, beta, p) (log(p/(1-p))-alpha)/beta
L50 <- lrPerc(alpha=coef(fit)[1], beta=coef(fit)[2], p=0.5)
lines(x=c(L50,L50,0), y=c(-100,0.5,0.5), lty=2, col=2)
text(x=L50, y=0.5, labels = bquote(L[mat]==.(round(L50,2))), pos=4, col=2 )
Mean length in catch:
df <- expand.grid(
length = as.numeric(dimnames(stk$stock.l@catch.n)[[1]])+0.5,
year = as.numeric(dimnames(stk$stock.l@catch.n)[[2]])
)
tmp <- apply(stk$stock.l@catch.n, 1:2, sum, na.rm=TRUE)
df$n = c(tmp)
df$lengthProd <- df$length * df$n
agg1 <- aggregate(cbind(lengthProd, n) ~ year, data = df,
FUN = sum, na.rm=TRUE)
agg1$meanLength <- agg1$lengthProd / agg1$n
plot(meanLength ~ year, agg1, t="b", ylim=range(df$length),
ylab = "mean length in catch (mm)")
abline(h = L50, lty=2, col=8)
Yearly fishing mortality:
df <- apply(stk$stock.a@harvest, 1:2, sum, na.rm=TRUE)
plot(df)
year <- as.numeric(dimnames(stk$stock.a@stock.n)[[2]])
steepness <- 1
FMmax <- 0.7*1.5
FMs <- FMmax / (1 + exp(-steepness * (year - 2004) ))
# spin-up with starting F
stk$harvest$params$FM <- FMs[1]
stk <- spinup.FLIBM(obj = stk, nyearsslope = 5, monitor = FALSE)
for(y in seq(year)){
stk$harvest$params$FM <- FMs[y]
stk <- adv.FLIBM(obj = stk, year = ac(year[y]), monitor = FALSE)
}
# plot stock numbers
plot(stk$stock.a@stock.n)
View biomass and spawning stock biomass development by year:
# B
dat1 <- as.data.frame(stock.n(stk$stock.a) * stock.wt(stk$stock.a),
date=TRUE)
dat1$age <- factor(dat1$age)
# SSB
dat2 <- as.data.frame(stock.n(stk$stock.a) * stock.wt(stk$stock.a) *
mat(stk$stock.a), date=TRUE)
dat2 <- aggregate(x = dat2$data, by = list(date = dat2$date),
FUN = sum, na.rm=TRUE)
dat2$age <- factor("0")
p <- ggplot2::ggplot(data = dat1,
aes(x=date, y=data, fill=age, group=age)) +
theme_set(theme_gray(base_size = 9)) +
geom_area(position = 'stack') +
scale_fill_brewer(palette="Spectral") +
geom_line(data = dat2, aes(x=date, y=x) ) +
ylab("Biomass (t)") +
theme(axis.text.x=element_text(angle = 90, hjust=1, vjust = 0.5))
print(p)
Mean length in catch:
Fishing mortality:
df <- apply(stk$stock.a@harvest, 1:2, sum, na.rm=TRUE)
plot(df)
More to come…
# pulsed recruitment
stk$rec$params$season_wt[] <- 0
stk$rec$params$season_wt[3:5] <- c(0.25, 1, 0.25)
# seasonal growth oscillation
stk$growth$params['C'] <- 0.75
# fishing mortality
stk$harvest$params['FM'] <- 0.7
# spinup and advance
stk <- adv.FLIBM(stk, years = ac(2000:2009), monitor = FALSE)
# plot stock numbers
plot(stk$stock.a@stock.n)
The function simplifySeason
can be used to collapse the
season dimension, resulting in the more typical yearly FLStock object.
The following example also removes the age zero group to define age = 1
as recruits. Without this change, recruits are not plotted since no
individuals of age = 0 are alive at the start of the year. This object
can be direcly plotted, providing a clear visualization of the main
stock attributed (recruitment, SSB, Catch, fishing mortality).
# summary stats
tmp <- simplifySeason(stk)
tmp <- tmp[ac(1:range(tmp)["max"]),] # remove age=0 group
plot(tmp)
LFQ
# plot length-frequency
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)))
# add log-normal noise to rec covar (one value per year)
stk$rec$covar[] <- rlnorm(n = dim(stk$rec$covar)[2], meanlog = 0, sdlog = 0.5)
stk <- adv.FLIBM(stk, years = ac(2000:2009), monitor = FALSE)
# plot stock numbers
plot(stk$stock.a@stock.n)
LFQ
# plot length-frequency
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)))
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.
This document is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license.