This is a working document, detailing the various data formats and assessment model configurations explored in finding the best model for the newly combined Plaice stock, across the Sub Divisions 21 - 32.
Throughout the document, the code used to undertake each step precedes the output and can be accessed by clicking on the drop-down arrows to the right of the page.
Data compilation is done externally, and this document deals only with the finalised input files. All models are build using the State-Space Assessment Model (SAM), and in this version of the document, all model input files and configurations are downloaded directly from the publicly available versions on stockassessment.org.
The file downloads and runs a series of assessments. Starting from a basic adaptation of the previous ple.27.3a.21-23 stock-annex, the document works through various methods of incorporating temporal change in the stock’s biological characteristics, before exploring the use of SAM’s “fleets” to differentiate between management areas (SD21 = North Sea; SDs22-32 = Baltic Sea) for the purpose of providing independent TAC advice directly from the model forecasts.
Along the way, we interrogate the model fits using various diagnostic tools, such as residual plots, leave-one-out analyses, and retrospective plots with mohn’s rho statistics. Some explanation is given to the benefits and limitations of each model, until a conclusion about which of the model structures investigated in this early evaluation will be taken forward for further model tuning and investigation with discard survivals and recreational data.
To run the scripts embedded in this document (if you would like to run the assessment yourself) you require the following packages and functions:
require(data.table)
require(reshape2)
require(knitr)
require(scales)
require(bookdown)
require(ggplot2)
require(ggthemes)
require(stringr)
require(lattice)
require(plyr)
require(plotly)
require(icesAdvice)
require("stockassessment") # available from https://github.com/fishfollower/SAM
require(parallel)
require(sf)
myprop <- function(df, nomo = "nomo", NoAtLngt = "CANoAtLngt") {
t <- sum(df$nomo)/sum(df$CANoAtLngt)
}
catchFrac <- function(x, nm, w, frac) {
F <- getF(x)
Z <- F + nm
N <- getN(x)
C <- F/Z * (1 - exp(-Z)) * N
return(sum(frac * w * C))
}
We also need to set a few variables manually.
The first is the DataYear, which is the year immediately preceeding the assessment year or the year for which the most recent data is available.
The second is the catch we (ICES) advised as TAC for the last assessment LastIcesAdvice, in tonnes. NOTE this value is taken from the advice sheet, to compare changes in advice, not changes in advice relative to actual Allowable Catches.
The third is whether this run is exploratory, or to be used to generate the final advice. This option will determine the number of iterations used to determine the short-term forecasts.
DataYear <- 2023
LastIcesAdvice <- 21735 # Advised catch for areas 21-32 (from stock splitting table in advice)
isFinal <- FALSE
# isFinal <- TRUE
This current document is generated as an exploratory run, where forecast simulations (not yet implemented) are run with 5000 iterations, except for the Fmsy, which is always run with 50,000 iterations. The absolute values in this output may differ from those in a final report and from final assessments. This is due to the fact that the final model for advice will run on stockassessment.org, where rounding rules, random number generators and other minor deviations will cause results to differ in minutia.
This assessment utilises the basic settings that were applied to the previous ple.21-23 stock, which already ran a SAM assessment. Some components are updated, such as the new indices produced from an updated index model for the new stock boundaries.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit = fitfromweb("ple.27.21-32_WKBPLAICE_2024_BASE")
basefit <- fitfromweb("ple.27.21-23_WGBFAS_2024_SPALY_v1")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. Last year’s assessment results are plotted alongside the current assessment, for comparison.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum <- as.data.frame(summary(fit))
asum$Year <- as.integer(row.names(summary(fit)))
asum <- cbind(asum, catchtable(fit, obs.show = TRUE))
colnames(asum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
# === Create dataframe of last year's fit for plotting ====
basum <- as.data.frame(summary(basefit))
basum$Year <- as.integer(row.names(summary(basefit)))
basum <- cbind(basum, catchtable(basefit, obs.show = T))
colnames(basum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 2.1: Base model estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment (configured similarly to the previous ple.27.21-23 stock assessment), and green is the previous ple.27.21-23 stock assessment.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 2.2: Base model annual fishing mortality estimates for ple.27.21.23 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment (configured similarly to the previous ple.27.21-23 stock assessment), and green is the previous ple.27.21-23 stock assessment.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum, mapping = aes(x = Year,
y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 2.3: Base model annual recruitment estimates for ple.27.21.23 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment (configured similarly to the previous ple.27.21-23 stock assessment), and green is the previous ple.27.21-23 stock assessment.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum, mapping = aes(x = Year,
y = CatchEst), colour = ebpal[1]) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum,
mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11], fill = ebpal[11]) +
ylab("Catch (tonnes)") + theme_clean())
Figure 2.4: Base model annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21.23 and point annual observations (points). Purple is the current assessment (configured similarly to the previous ple.27.21-23 stock assessment), and green is the previous ple.27.21-23 stock assessment.
### Data Prep
F <- as.data.frame(faytable(fit))
F$Year <- rownames(F)
rownames(F) <- NULL
FLng <- melt(data = F, id.vars = "Year", variable.name = "Age", value.name = "F_atAge")
FLng$Year <- as.numeric(as.character(FLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = FLng, mapping = aes(x = Year, y = F_atAge, colour = Age)) +
scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(FLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 2.5: Partial F at age; base model
sdplot(fit, marg = c(5, 4, 1, 1))
Figure 2.6: Standard Devations by Fleet; base model
First we can look at how the model fits the observations per age in the catches, as well as in the quarter one survey index and the quarter three/four survey index.
## Extract data from model object
logObs <- split(fit$data$logobs, ceiling(seq_along(fit$data$logobs)/fit$data$maxAgePerFleet[1]))
logPred <- split(fit$rep$predObs, ceiling(seq_along(fit$rep$predObs)/fit$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs[[i]])
logObs_df <- rbind(logObs_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred[[i]])
logPred_df <- rbind(logPred_df, temp_df)
}
## Plot
logPred_df$age <- as.character(logPred_df$age)
logObs_df$age <- as.character(logObs_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_df[logObs_df$fleet == 1, ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_df[logPred_df$fleet == 1, ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 2.7: Model fit to catch data by age (values on the model link function scale).
# fitplot(fit, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_df[logObs_df$fleet == 2, ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_df[logPred_df$fleet == 2, ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 2.8: Model fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_df[logObs_df$fleet == 3, ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_df[logPred_df$fleet == 3, ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 2.9: Model fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit, fleets=3)
Now we can also calculate the residuals of our observations from our model. When we plot these by age over time, we can see that there is an excess of over-estimation in the numbers of fish in recent years. This trend may be seen in one fleet, where seasonal selectivity or changes in fishing patterns may mean that one fleet catches fewer fish than expected, however, we can see this across all three fleets and for almost all ages. This indicates that there is probably some un-accounted for mortality, meaning the model estimates that there should be more fish available to the different fleets, than we observe being caught.
resid_base <- residuals(fit)
resid_base_df <- data.frame(year = resid_base$year, fleet = resid_base$fleet, age = resid_base$age,
observation = resid_base$observation, mean = resid_base$mean, residual = resid_base$residual)
resid_base_df$fleetName <- ifelse(resid_base_df$fleet == 1, attributes(resid_base)$fleetNames[1],
ifelse(resid_base_df$fleet == 2, attributes(resid_base)$fleetNames[2], ifelse(resid_base_df$fleet ==
3, attributes(resid_base)$fleetNames[3], NA)))
resid_base_df$fleetAltName <- ifelse(resid_base_df$fleet == 1, "Residual Fishing",
ifelse(resid_base_df$fleet == 2, "Q1 Surveys", ifelse(resid_base_df$fleet ==
3, "Q3/4 Surveys", NA)))
ggplotly(ggplot(resid_base_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 2.10: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual).
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. Q1 surveys), this is especially important.
if (!all(fit$conf$obsCorStruct == "ID")) {
corplot(fit)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 2.11: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom).
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
In the case of this model, we can see that both surveys contribute to increased estimations of recruitment. This is likely due to the much smaller mesh size in the survey gears meaning a much better catchability of the smaller sizes and hence, age 1 fish.
LO <- leaveout(fit)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(summary(LO$`w.o. Q34IBTS+BITS+CodSD21-25`))),
SSB = q3mat$SSB, Fbar = q3mat$`Fbar(3-5)`, R_age1 = q3mat$`R(age 1)`, CatchEst = catchtable(LO$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = nrow(q3mat)))
asum$series <- rep("full", times = nrow(asum))
asumi <- asum[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum <- rbind(woq1, woq34, asumi)
# ssbplot(LO)
lossb <- ggplotly(ggplot() + geom_line(data = losum, mapping = aes(x = Year, y = SSB,
colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year, ymin = SSBlow,
ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof <- ggplotly(ggplot() + geom_line(data = losum, mapping = aes(x = Year, y = Fbar,
colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year, ymin = Flow,
ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec <- ggplotly(ggplot() + geom_line(data = losum, mapping = aes(x = Year, y = R_age1,
colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca <- ggplotly(ggplot() + geom_line(data = losum, mapping = aes(x = Year, y = CatchEst,
colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum, mapping = aes(x = Year,
y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") +
scale_colour_manual(values = c("black", ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb, lof, lorec, loca, nrows = 2, shareX = TRUE, titleY = TRUE),
showlegend = FALSE)
Figure 2.12: Leave-one-out re-fits (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how a model performs in different data situations.
RETRO <- retro(fit, year = 5)
rho_base <- mohn(RETRO)
## Make RETROs in better plotting format
ret_df <- asum
for (i in 1:length(RETRO)) {
tsum <- as.data.frame(summary(RETRO[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO[[i]])))
tsum <- cbind(tsum, catchtable(RETRO[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
ret_df <- rbind(ret_df, tsum)
}
ret_df$series <- factor(ret_df$series, levels = c("full", "1", "2", "3", "4", "5"))
In this case the “peels” of each subsequent model re-fit remain within the main model’s bounds of uncertainty, and the Mohn’s rho statistics are acceptable.
# ssbplot(RETRO)
retssb <- layout(ggplotly(ggplot() + geom_line(data = ret_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + annotate(geom = "text",
y = max(asum$SSBhigh) * 0.85, x = ((max(asum$Year) - min(asum$Year)) * 0.2) +
min(asum$Year), label = paste0("Mohn's Rho = ", round(rho_base[2], digits = 3))) +
scale_colour_manual(values = c("black", ebpal[(length(ebpal) - 4):length(ebpal)])) +
theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf <- layout(ggplotly(ggplot() + geom_line(data = ret_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + annotate(geom = "text",
y = max(asum$Fhigh) * 0.85, x = ((max(asum$Year) - min(asum$Year)) * 0.8) + min(asum$Year),
label = paste0("Mohn's Rho = ", round(rho_base[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec <- layout(ggplotly(ggplot() + geom_line(data = ret_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + annotate(geom = "text",
y = max(asum$Rhigh) * 0.85, x = ((max(asum$Year) - min(asum$Year)) * 0.2) + min(asum$Year),
label = paste0("Mohn's Rho = ", round(rho_base[1], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca <- ggplotly(ggplot() + geom_line(data = ret_df, mapping = aes(x = Year, y = CatchEst,
colour = series)) + geom_ribbon(data = asum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum, mapping = aes(x = Year,
y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") +
scale_colour_manual(values = c("black", ebpal[(length(ebpal) - 4):length(ebpal)])) +
theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95)))
ret_fp <- style(subplot(retssb, retf, retrec, retca, nrows = 2, shareX = FALSE, shareY = FALSE,
titleY = TRUE), showlegend = FALSE, traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 2.13: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the base model, configured similarly to the old Ple.27.21-32 stock.
This direct translation of the Ple.27.21-23 stock annex to the new stock seems to work reasonably well. However, there remain a few issues with this:
Annually observed biological data, such as age, length, weight and maturity, can vary greatly due to observation error and small sample sizes. Therefore using each year’s observations as “truth” to fit the model to will introduce a lot of noise that cannot be accounted for. Often, fixed weights at age and maturity ogives are used for all years in an assessment model, as was the case in the previous Ple.27.21-23 stock and maintained in the above “baseline” assessment.
To overcome the excess variation from noisey annual data, while still accounting for temporal change in things like weight at age and maturity at age, there are multiple options. One can try to smooth the data, trusting that the interannually variability averages out when comparing neighbouring years, then put these smoothed data into the model fit as “truths”. Alternately one can allow the model to fit to the observations as observations, with associated errors. The former is often simpler and easier to explain to other experts and any stakeholders interested in such details, while the latter incorporates more of the uncertaintity derived from the observations into the model estimations themselves. However, with all of the extra parameters that need estimating, when including observation error in the model, there is a risk that the model will not converge.
These are trade-offs both in the operational and scientific robustness sense. In this section, we investigate ways of incorporating temporal change that are operationally realistic / acceptable, and as robust as possible.
This assessment builds on the previous one by providing the Assessment model with the annual observations of stock biological data (maturity ogives and stock weights at age) instead of a single time-series average. The assessment model then uses these as observations in fitting the model.
The first attempt, which included in-model smoothing of both maturity ogives and stock weights at age, resulted in model non-convergence. Coupling variance between older ages for these two variables did not remedy the non-convergence. Therefore, because the largest and most monotonic change has been seen in the stock weights at age, we continued with a model that uses five year sliding window means to create the input data for maturity ogives, while using the in-model smoothing over the raw annual estimates for stock weight at age
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_biopar = fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioPar_sw")
basefit <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BASE")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_biopar <- as.data.frame(summary(fit_biopar))
asum_biopar$Year <- as.integer(row.names(summary(fit_biopar)))
asum_biopar <- cbind(asum_biopar, catchtable(fit_biopar, obs.show = TRUE))
colnames(asum_biopar) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
# === Create dataframe of last year's fit for plotting ====
basum <- as.data.frame(summary(basefit))
basum$Year <- as.integer(row.names(summary(basefit)))
basum <- cbind(basum, catchtable(basefit, obs.show = T))
colnames(basum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_biopar, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 3.1: Biopar model estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the baseline assessment, described first in this document.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_biopar, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.2: Biopar model annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the baseline assessment, described first in this document.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_biopar,
mapping = aes(x = Year, y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_biopar,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) +
theme_clean())
Figure 3.3: Biopar model annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the baseline assessment, described first in this document.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_biopar,
mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) + geom_ribbon(data = asum_biopar,
mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1],
alpha = 0.3) + geom_point(data = asum_biopar, mapping = aes(x = Year, y = CatchObs),
shape = 3, colour = ebpal[11], fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 3.4: Biopar model annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the baseline assessment, described first in this document.
### Data Prep
F_biopar <- as.data.frame(faytable(fit_biopar))
F_biopar$Year <- rownames(F_biopar)
rownames(F_biopar) <- NULL
F_bioparLng <- melt(data = F_biopar, id.vars = "Year", variable.name = "Age", value.name = "F_atAge")
F_bioparLng$Year <- as.numeric(as.character(F_bioparLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_bioparLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_bioparLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 3.5: Partial F at age; Biopar model
sdplot(fit_biopar, marg = c(5, 4, 1, 1))
Figure 3.6: Standard Devations by Fleet; biopar model
First we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_biopar <- split(fit_biopar$data$logobs, ceiling(seq_along(fit_biopar$data$logobs)/fit_biopar$data$maxAgePerFleet[1]))
logPred_biopar <- split(fit_biopar$rep$predObs, ceiling(seq_along(fit_biopar$rep$predObs)/fit_biopar$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_biopar_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_biopar)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_biopar[[i]])
logObs_biopar_df <- rbind(logObs_biopar_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_biopar_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_biopar)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_biopar[[i]])
logPred_biopar_df <- rbind(logPred_biopar_df, temp_df)
}
## Plot
logPred_biopar_df$age <- as.character(logPred_biopar_df$age)
logObs_biopar_df$age <- as.character(logObs_biopar_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_biopar_df[logObs_biopar_df$fleet == 1, ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_biopar_df[logPred_biopar_df$fleet == 1, ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.7: Biopar model fit to catch data by age (values on the model link function scale).
# fitplot(fit_biopar, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_biopar_df[logObs_biopar_df$fleet == 2, ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_biopar_df[logPred_biopar_df$fleet == 2, ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.8: Biopar model fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_biopar, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_biopar_df[logObs_biopar_df$fleet == 3, ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_biopar_df[logPred_biopar_df$fleet == 3, ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.9: Biopar model fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_biopar, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_biopar <- residuals(fit_biopar)
resid_biopar_df <- data.frame(year = resid_biopar$year, fleet = resid_biopar$fleet,
age = resid_biopar$age, observation = resid_biopar$observation, mean = resid_biopar$mean,
residual = resid_biopar$residual)
resid_biopar_df$fleetName <- ifelse(resid_biopar_df$fleet == 1, attributes(resid_biopar)$fleetNames[1],
ifelse(resid_biopar_df$fleet == 2, attributes(resid_biopar)$fleetNames[2], ifelse(resid_biopar_df$fleet ==
3, attributes(resid_biopar)$fleetNames[3], NA)))
resid_biopar_df$fleetAltName <- ifelse(resid_biopar_df$fleet == 1, "Residual Fishing",
ifelse(resid_biopar_df$fleet == 2, "Q1 Surveys", ifelse(resid_biopar_df$fleet ==
3, "Q3/4 Surveys", NA)))
Again we see strong residual patterns indicating that this model over-estimates what all fleets will catch/sample in the latest years. This remains probably due to unaccounted for natural mortality.
ggplotly(ggplot(resid_biopar_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 3.10: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual) from the biopar model.
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. Q1 surveys), this is especially important.
if (!all(fit_biopar$conf$obsCorStruct == "ID")) {
corplot(fit_biopar)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 3.11: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the biopar model.
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
LO_biopar <- leaveout(fit_biopar)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_biopar$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_biopar$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_biopar$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_biopar$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(summary(LO_biopar$`w.o. Q34IBTS+BITS+CodSD21-25`))),
SSB = q3mat$SSB, Fbar = q3mat$`Fbar(3-5)`, R_age1 = q3mat$`R(age 1)`, CatchEst = catchtable(LO_biopar$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = nrow(q3mat)))
asum_biopar$series <- rep("full", times = nrow(asum_biopar))
asumi_biopar <- asum_biopar[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_biopar <- rbind(woq1, woq34, asumi_biopar)
# ssbplot(LO)
lossb_biopar <- ggplotly(ggplot() + geom_line(data = losum_biopar, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_biopar <- ggplotly(ggplot() + geom_line(data = losum_biopar, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_biopar <- ggplotly(ggplot() + geom_line(data = losum_biopar, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_biopar <- ggplotly(ggplot() + geom_line(data = losum_biopar, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum_biopar,
mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) +
ylab("Catch (tonnes)") + scale_colour_manual(values = c("black", ebpal[8:9])) +
theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95)))
layout(subplot(lossb_biopar, lof_biopar, lorec_biopar, loca_biopar, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 3.12: Leave-one-out re-fits for the biopar model (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
RETRO_biopar <- retro(fit_biopar, year = 5)
rho_biopar <- mohn(RETRO_biopar)
## Make RETROs in better plotting format
ret_biopar_df <- asum_biopar
ret_biopar_df$series <- rep("full", nrow(ret_biopar_df))
for (i in 1:length(RETRO_biopar)) {
tsum <- as.data.frame(summary(RETRO_biopar[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_biopar[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_biopar[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
ret_biopar_df <- rbind(ret_biopar_df, tsum)
}
ret_biopar_df$series <- factor(ret_biopar_df$series, levels = c("full", "1", "2",
"3", "4", "5"))
In this model run, we see that the majority of peels remain within the confidence bounds of the main model, indicating that they don’t significantly differ when the model is refit to less data, however, the first peel on F is very close to coming out of these bounds. The Mohn’s rho values are OK, but F is looking a little high.
# ssbplot(RETRO)
retssb_biopar <- layout(ggplotly(ggplot() + geom_line(data = ret_biopar_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + annotate(geom = "text",
y = max(asum_biopar$SSBhigh) * 0.85, x = ((max(asum_biopar$Year) - min(asum_biopar$Year)) *
0.2) + min(asum_biopar$Year), label = paste0("Mohn's Rho = ", round(rho_biopar[2],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_biopar <- layout(ggplotly(ggplot() + geom_line(data = ret_biopar_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + annotate(geom = "text",
y = max(asum_biopar$Fhigh) * 0.85, x = ((max(asum_biopar$Year) - min(asum_biopar$Year)) *
0.8) + min(asum_biopar$Year), label = paste0("Mohn's Rho = ", round(rho_biopar[3],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_biopar <- layout(ggplotly(ggplot() + geom_line(data = ret_biopar_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + annotate(geom = "text",
y = max(asum_biopar$Rhigh) * 0.85, x = ((max(asum_biopar$Year) - min(asum_biopar$Year)) *
0.2) + min(asum_biopar$Year), label = paste0("Mohn's Rho = ", round(rho_biopar[1],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_biopar <- ggplotly(ggplot() + geom_line(data = ret_biopar_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_biopar, mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum_biopar,
mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) +
ylab("Catch (tonnes)") + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_biopar, retf_biopar, retrec_biopar, retca_biopar,
nrows = 2, shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE,
traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 3.13: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising internal smoothing during model fitting for stock weight at age and maturity.
The use of the in-model biological parameter estimation seems to work reasonably well. However, there remain a few issues with this:
This assessment uses a simpler approach to smoothing the noisy annual observation data for stock-weights and maturity ogives. That is a sliding window mean, which is potentially easier to explain and defend to different audiences. In preparing the input files, the sliding window mean simply averages the preceding x-1 years’ observations with each subsequent year to reduce the effect of annual variability. In this case x = 3, for a three year sliding window.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_3ysw = fitfromweb("ple.27.21-32_WKBPLAICE_2024_3ysw")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_3ysw <- as.data.frame(summary(fit_3ysw))
asum_3ysw$Year <- as.integer(row.names(summary(fit_3ysw)))
ct_3ysw <- as.data.frame(catchtable(fit_3ysw, obs.show = TRUE))
ct_3ysw$Year <- as.integer(rownames(ct_3ysw))
rownames(ct_3ysw) <- NULL
ct_3ysw <- rbind(ct_3ysw, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_3ysw <- merge(x = asum_3ysw, y = ct_3ysw, by = "Year")
colnames(asum_3ysw) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, as in the biopar model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_3ysw, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_3ysw, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 3.14: Three year sliding window model estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the baseline assessment, described first in this document.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_3ysw, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_3ysw, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.15: Three year sliding window model annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the baseline assessment, described first in this document.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_3ysw, mapping = aes(x = Year,
y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_3ysw, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.16: Three year sliding window model annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the baseline assessment, described first in this document.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_3ysw[asum_3ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) +
geom_ribbon(data = asum_3ysw[asum_3ysw$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_3ysw[asum_3ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 3.17: Three year sliding window model annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the baseline assessment, described first in this document.
### Data Prep
F_3ysw <- as.data.frame(faytable(fit_3ysw))
F_3ysw$Year <- rownames(F_3ysw)
rownames(F_3ysw) <- NULL
F_3yswLng <- melt(data = F_3ysw, id.vars = "Year", variable.name = "Age", value.name = "F_atAge")
F_3yswLng$Year <- as.numeric(as.character(F_3yswLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_3yswLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_3yswLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 3.18: Partial F at age; Three year sliding window model
sdplot(fit_3ysw, marg = c(5, 4, 1, 1))
Figure 3.19: Standard Devations by Fleet; base model
First we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_3ysw <- split(fit_3ysw$data$logobs, ceiling(seq_along(fit_3ysw$data$logobs)/fit_3ysw$data$maxAgePerFleet[1]))
logPred_3ysw <- split(fit_3ysw$rep$predObs, ceiling(seq_along(fit_3ysw$rep$predObs)/fit_3ysw$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_3ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_3ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_3ysw[[i]])
logObs_3ysw_df <- rbind(logObs_3ysw_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_3ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_3ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_3ysw[[i]])
logPred_3ysw_df <- rbind(logPred_3ysw_df, temp_df)
}
## Plot
logPred_3ysw_df$age <- as.character(logPred_3ysw_df$age)
logObs_3ysw_df$age <- as.character(logObs_3ysw_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_3ysw_df[logObs_3ysw_df$fleet == 1 & logObs_3ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_3ysw_df[logPred_3ysw_df$fleet == 1 & logPred_3ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.20: Three year sliding window model fit to catch data by age (values on the model link function scale).
# fitplot(fit_3ysw, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_3ysw_df[logObs_3ysw_df$fleet == 2 & logObs_3ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_3ysw_df[logPred_3ysw_df$fleet == 2 & logPred_3ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.21: Three year sliding window model fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_3ysw, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_3ysw_df[logObs_3ysw_df$fleet == 3 & logObs_3ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_3ysw_df[logPred_3ysw_df$fleet == 3 & logPred_3ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.22: Three year sliding window model fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_3ysw, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_3ysw <- residuals(fit_3ysw)
resid_3ysw_df <- data.frame(year = resid_3ysw$year, fleet = resid_3ysw$fleet, age = resid_3ysw$age,
observation = resid_3ysw$observation, mean = resid_3ysw$mean, residual = resid_3ysw$residual)
resid_3ysw_df$fleetName <- ifelse(resid_3ysw_df$fleet == 1, attributes(resid_3ysw)$fleetNames[1],
ifelse(resid_3ysw_df$fleet == 2, attributes(resid_3ysw)$fleetNames[2], ifelse(resid_3ysw_df$fleet ==
3, attributes(resid_3ysw)$fleetNames[3], NA)))
resid_3ysw_df$fleetAltName <- ifelse(resid_3ysw_df$fleet == 1, "Residual Fishing",
ifelse(resid_3ysw_df$fleet == 2, "Q1 Surveys", ifelse(resid_3ysw_df$fleet ==
3, "Q3/4 Surveys", NA)))
In this model, the residual patterns in the later years remains strong, indicating that we still need to account for that missing mortality. . .
ggplotly(ggplot(resid_3ysw_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 3.23: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the three year sliding window model.
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. No fleets match this criteria in this model configuration.), this is especially important.
if (!all(fit_3ysw$conf$obsCorStruct == "ID")) {
corplot(fit_3ysw)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
## [1] "No correlation structure configured for residuals across age."
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
LO_3ysw <- leaveout(fit_3ysw)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_3ysw$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_3ysw$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_3ysw$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_3ysw$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(q3mat[(2002:DataYear) - 2001, ])),
SSB = q3mat[(2002:DataYear) - 2001, "SSB"], Fbar = q3mat[(2002:DataYear) - 2001,
"Fbar(3-5)"], R_age1 = q3mat[(2002:DataYear) - 2001, "R(age 1)"], CatchEst = catchtable(LO_3ysw$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_3ysw$series <- rep("full", times = nrow(asum_3ysw))
asumi_3ysw <- asum_3ysw[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_3ysw <- rbind(woq1, woq34, asumi_3ysw)
# ssbplot(LO)
lossb_3ysw <- ggplotly(ggplot() + geom_line(data = losum_3ysw, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_3ysw, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_3ysw <- ggplotly(ggplot() + geom_line(data = losum_3ysw, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_3ysw, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_3ysw <- ggplotly(ggplot() + geom_line(data = losum_3ysw, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_3ysw, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_3ysw <- ggplotly(ggplot() + geom_line(data = losum_3ysw, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_3ysw[asum_3ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_3ysw[asum_3ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb_3ysw, lof_3ysw, lorec_3ysw, loca_3ysw, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 3.24: Leave-one-out re-fits for the 3ysw model (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
In this case, our model fit included the within-assessement year Q1 survey index as a tuning fleet, therefore, we are doing the peels and calculating the Mohn’s rho on data truncated by one year, so that we aren’t judging the model fit with the year in which it doesn’t have catch data nor the second index.
RETRO_3ysw <- retro(fit_3ysw, year = 5)
rho_3ysw <- mohn(RETRO_3ysw, lag = 1)
## Make RETROs in better plotting format
ret_3ysw_df <- asum_3ysw[asum_3ysw$Year != max(asum_3ysw$Year), ]
for (i in 1:length(RETRO_3ysw)) {
tsum <- as.data.frame(summary(RETRO_3ysw[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_3ysw[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_3ysw[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_3ysw_df <- rbind(ret_3ysw_df, tsum)
}
ret_3ysw_df$series <- factor(ret_3ysw_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_3ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_3ysw_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_3ysw[asum_3ysw$Year != max(asum_3ysw$Year),
], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) +
annotate(geom = "text", y = max(ret_3ysw_df$SSBhigh) * 0.85, x = ((max(ret_3ysw_df$Year) -
min(ret_3ysw_df$Year)) * 0.2) + min(ret_3ysw_df$Year), label = paste0("Mohn's Rho = ",
round(rho_3ysw[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_3ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_3ysw_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_3ysw[asum_3ysw$Year !=
max(asum_3ysw$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_3ysw_df$Fhigh) * 0.85, x = ((max(ret_3ysw_df$Year) -
min(ret_3ysw_df$Year)) * 0.8) + min(ret_3ysw_df$Year), label = paste0("Mohn's Rho = ",
round(rho_3ysw[3], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_3ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_3ysw_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_3ysw[asum_3ysw$Year !=
max(asum_3ysw$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_3ysw_df$Rhigh) * 0.85, x = ((max(ret_3ysw_df$Year) -
min(ret_3ysw_df$Year)) * 0.2) + min(ret_3ysw_df$Year), label = paste0("Mohn's Rho = ",
round(rho_3ysw[1], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_3ysw <- ggplotly(ggplot() + geom_line(data = ret_3ysw_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_3ysw[asum_3ysw$Year !=
max(asum_3ysw$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_3ysw[asum_3ysw$Year !=
max(asum_3ysw$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_3ysw, retf_3ysw, retrec_3ysw, retca_3ysw, nrows = 2,
shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 *
4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0, yaxis = list(automargin = TRUE), xaxis = list(automargin = TRUE),
yaxis2 = list(automargin = TRUE), xaxis2 = list(automargin = TRUE), yaxis3 = list(automargin = TRUE),
xaxis3 = list(automargin = TRUE), yaxis4 = list(automargin = TRUE), xaxis4 = list(automargin = TRUE)))
Figure 3.25: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising three year sliding window means for stock weight at age and maturity.
The use of a three year sliding window mean on the biological data (explicitly, the stock weights at age and the maturity ogive) decreases the estimates of SSB in later years, probably due to the observed reduction in weights at age. However, this relatively short window does increases the variability in model estimates of SSB, F and Recruitment. Furthermore, this model does not account for different catches in different management areas, and therefore, advice generated from short-term forecasts will have to be manually split based on basic ratios.
This assessment uses a simpler approach to smoothing the noisy annual observation data for stock-weights and maturity ogives. That is a sliding window mean, which is potentially easier to explain and defend to different audiences. In preparing the input files, the sliding window mean simply averages the preceding x-1 years’ observations with each subsequent year to reduce the effect of annual variability. In this case x = 5, for a five year sliding window.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_5ysw = fitfromweb("ple.27.21-32_WKBPLAICE_2024_5ysw")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_5ysw <- as.data.frame(summary(fit_5ysw))
asum_5ysw$Year <- as.integer(row.names(summary(fit_5ysw)))
ct_5ysw <- as.data.frame(catchtable(fit_5ysw, obs.show = TRUE))
ct_5ysw$Year <- as.integer(rownames(ct_5ysw))
rownames(ct_5ysw) <- NULL
ct_5ysw <- rbind(ct_5ysw, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_5ysw <- merge(x = asum_5ysw, y = ct_5ysw, by = "Year")
colnames(asum_5ysw) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, as in the biopar model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5ysw, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_5ysw, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 3.26: Five year sliding window model estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the baseline assessment, described first in this document.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5ysw, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_5ysw, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.27: Five year sliding window model annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the baseline assessment, described first in this document.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5ysw, mapping = aes(x = Year,
y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_5ysw, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.28: Five year sliding window model annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the baseline assessment, described first in this document.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5ysw[asum_5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) +
geom_ribbon(data = asum_5ysw[asum_5ysw$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_5ysw[asum_5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 3.29: Five year sliding window model annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the baseline assessment, described first in this document.
### Data Prep
F_5ysw <- as.data.frame(faytable(fit_5ysw))
F_5ysw$Year <- rownames(F_5ysw)
rownames(F_5ysw) <- NULL
F_5yswLng <- melt(data = F_5ysw, id.vars = "Year", variable.name = "Age", value.name = "F_atAge")
F_5yswLng$Year <- as.numeric(as.character(F_5yswLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_5yswLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_5yswLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 3.30: Partial F at age; Five year sliding window model
sdplot(fit_5ysw, marg = c(5, 4, 1, 1))
Figure 3.31: Standard Devations by Fleet; base model
First we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_5ysw <- split(fit_5ysw$data$logobs, ceiling(seq_along(fit_5ysw$data$logobs)/fit_5ysw$data$maxAgePerFleet[1]))
logPred_5ysw <- split(fit_5ysw$rep$predObs, ceiling(seq_along(fit_5ysw$rep$predObs)/fit_5ysw$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_5ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_5ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_5ysw[[i]])
logObs_5ysw_df <- rbind(logObs_5ysw_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_5ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_5ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_5ysw[[i]])
logPred_5ysw_df <- rbind(logPred_5ysw_df, temp_df)
}
## Plot
logPred_5ysw_df$age <- as.character(logPred_5ysw_df$age)
logObs_5ysw_df$age <- as.character(logObs_5ysw_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_5ysw_df[logObs_5ysw_df$fleet == 1 & logObs_5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_5ysw_df[logPred_5ysw_df$fleet == 1 & logPred_5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.32: Five year sliding window model fit to catch data by age (values on the model link function scale).
# fitplot(fit_5ysw, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_5ysw_df[logObs_5ysw_df$fleet == 2 & logObs_5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_5ysw_df[logPred_5ysw_df$fleet == 2 & logPred_5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.33: Five year sliding window model fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_5ysw, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_5ysw_df[logObs_5ysw_df$fleet == 3 & logObs_5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_5ysw_df[logPred_5ysw_df$fleet == 3 & logPred_5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.34: Five year sliding window model fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_5ysw, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_5ysw <- residuals(fit_5ysw)
resid_5ysw_df <- data.frame(year = resid_5ysw$year, fleet = resid_5ysw$fleet, age = resid_5ysw$age,
observation = resid_5ysw$observation, mean = resid_5ysw$mean, residual = resid_5ysw$residual)
resid_5ysw_df$fleetName <- ifelse(resid_5ysw_df$fleet == 1, attributes(resid_5ysw)$fleetNames[1],
ifelse(resid_5ysw_df$fleet == 2, attributes(resid_5ysw)$fleetNames[2], ifelse(resid_5ysw_df$fleet ==
3, attributes(resid_5ysw)$fleetNames[3], NA)))
resid_5ysw_df$fleetAltName <- ifelse(resid_5ysw_df$fleet == 1, "Residual Fishing",
ifelse(resid_5ysw_df$fleet == 2, "Q1 Surveys", ifelse(resid_5ysw_df$fleet ==
3, "Q3/4 Surveys", NA)))
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. Q3/4 surveys), this is especially important.
if (!all(fit_5ysw$conf$obsCorStruct == "ID")) {
corplot(fit_5ysw)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 3.35: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the 5ysw model.
ggplotly(ggplot(resid_5ysw_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 3.36: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the five year sliding window model.
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
LO_5ysw <- leaveout(fit_5ysw)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_5ysw$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_5ysw$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_5ysw$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_5ysw$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(q3mat[(2002:DataYear) - 2001, ])),
SSB = q3mat[(2002:DataYear) - 2001, "SSB"], Fbar = q3mat[(2002:DataYear) - 2001,
"Fbar(3-5)"], R_age1 = q3mat[(2002:DataYear) - 2001, "R(age 1)"], CatchEst = catchtable(LO_5ysw$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_5ysw$series <- rep("full", times = nrow(asum_5ysw))
asumi_5ysw <- asum_5ysw[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_5ysw <- rbind(woq1, woq34, asumi_5ysw)
# ssbplot(LO)
lossb_5ysw <- ggplotly(ggplot() + geom_line(data = losum_5ysw, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_5ysw, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_5ysw <- ggplotly(ggplot() + geom_line(data = losum_5ysw, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_5ysw, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_5ysw <- ggplotly(ggplot() + geom_line(data = losum_5ysw, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_5ysw, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_5ysw <- ggplotly(ggplot() + geom_line(data = losum_5ysw, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_5ysw[asum_5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_5ysw[asum_5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb_5ysw, lof_5ysw, lorec_5ysw, loca_5ysw, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 3.37: Leave-one-out re-fits for the 5ysw model (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
In this case, our model fit included the within-assessement year Q1 survey index as a tuning fleet, therefore, we are doing the peels and calculating the Mohn’s rho on data truncated by one year, so that we aren’t judging the model fit with the year in which it doesn’t have catch data nor the second index.
RETRO_5ysw <- retro(fit_5ysw, year = 5)
rho_5ysw <- mohn(RETRO_5ysw, lag = 1)
## Make RETROs in better plotting format
ret_5ysw_df <- asum_5ysw[asum_5ysw$Year != max(asum_5ysw$Year), ]
for (i in 1:length(RETRO_5ysw)) {
tsum <- as.data.frame(summary(RETRO_5ysw[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_5ysw[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_5ysw[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_5ysw_df <- rbind(ret_5ysw_df, tsum)
}
ret_5ysw_df$series <- factor(ret_5ysw_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_5ysw_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_5ysw[asum_5ysw$Year != max(asum_5ysw$Year),
], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) +
annotate(geom = "text", y = max(ret_5ysw_df$SSBhigh) * 0.85, x = ((max(ret_5ysw_df$Year) -
min(ret_5ysw_df$Year)) * 0.2) + min(ret_5ysw_df$Year), label = paste0("Mohn's Rho = ",
round(rho_5ysw[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_5ysw_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_5ysw[asum_5ysw$Year !=
max(asum_5ysw$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_5ysw_df$Fhigh) * 0.85, x = ((max(ret_5ysw_df$Year) -
min(ret_5ysw_df$Year)) * 0.8) + min(ret_5ysw_df$Year), label = paste0("Mohn's Rho = ",
round(rho_5ysw[3], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_5ysw_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_5ysw[asum_5ysw$Year !=
max(asum_5ysw$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_5ysw_df$Rhigh) * 0.85, x = ((max(ret_5ysw_df$Year) -
min(ret_5ysw_df$Year)) * 0.2) + min(ret_5ysw_df$Year), label = paste0("Mohn's Rho = ",
round(rho_5ysw[1], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_5ysw <- ggplotly(ggplot() + geom_line(data = ret_5ysw_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_5ysw[asum_5ysw$Year !=
max(asum_5ysw$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_5ysw[asum_5ysw$Year !=
max(asum_5ysw$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_5ysw, retf_5ysw, retrec_5ysw, retca_5ysw, nrows = 2,
shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 *
4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 3.38: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising five year sliding window means for stock weight at age and maturity.
The use of a five year sliding window mean on the biological data (explicitly, the stock weights at age and the maturity ogive) decreases the estimates of SSB in later years but not by as much as the three year sliding window mean, nor the biopar model. Again, this decrease is probably due to the observed reduction in weights at age, but the smaller magnitude is due to the larger window mean, which responds slower to the changes. However, this longer window means that the variability in model estimates of SSB, F and Recruitment, remains much lower. Furthermore, this model does not account for different catches in different management areas, and therefore, advice generated from short-term forecasts will have to be manually split based on basic ratios.
This assessment uses a simpler approach to smoothing the noisy annual observation data for stock-weights and maturity ogives. That is a sliding window mean, which is potentially easier to explain and defend to different audiences. In preparing the input files, the sliding window mean simply averages the preceding x-1 years’ observations with each subsequent year to reduce the effect of annual variability. In this case x = 5, for a five year sliding window.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_5yswSV = fitfromweb("ple.27.21-32_WKBPLAICE_2024_5ysw_survCor")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_5yswSV <- as.data.frame(summary(fit_5yswSV))
asum_5yswSV$Year <- as.integer(row.names(summary(fit_5yswSV)))
ct_5yswSV <- as.data.frame(catchtable(fit_5yswSV, obs.show = TRUE))
ct_5yswSV$Year <- as.integer(rownames(ct_5yswSV))
rownames(ct_5yswSV) <- NULL
ct_5yswSV <- rbind(ct_5yswSV, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_5yswSV <- merge(x = asum_5yswSV, y = ct_5yswSV, by = "Year")
colnames(asum_5yswSV) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, as in the biopar model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5yswSV, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 3.39: Five year sliding window model (survey variation) estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the baseline assessment, described first in this document.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5yswSV, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.40: Five year sliding window model (survey variation) annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the baseline assessment, described first in this document.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5yswSV,
mapping = aes(x = Year, y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_5yswSV,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) +
theme_clean())
Figure 3.41: Five year sliding window model (survey variation) annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the baseline assessment, described first in this document.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_5yswSV[asum_5yswSV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) +
geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_5yswSV[asum_5yswSV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 3.42: Five year sliding window model (survey variation) annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the baseline assessment, described first in this document.
### Data Prep
F_5yswSV <- as.data.frame(faytable(fit_5yswSV))
F_5yswSV$Year <- rownames(F_5yswSV)
rownames(F_5yswSV) <- NULL
F_5yswSVLng <- melt(data = F_5yswSV, id.vars = "Year", variable.name = "Age", value.name = "F_atAge")
F_5yswSVLng$Year <- as.numeric(as.character(F_5yswSVLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_5yswSVLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_5yswSVLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 3.43: Partial F at age; Five year sliding window model (survey variation)
sdplot(fit_5yswSV, marg = c(5, 4, 1, 1))
Figure 3.44: Standard Devations by Fleet.
First we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_5yswSV <- split(fit_5yswSV$data$logobs, ceiling(seq_along(fit_5yswSV$data$logobs)/fit_5yswSV$data$maxAgePerFleet[1]))
logPred_5yswSV <- split(fit_5yswSV$rep$predObs, ceiling(seq_along(fit_5yswSV$rep$predObs)/fit_5yswSV$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_5yswSV_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_5yswSV)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_5yswSV[[i]])
logObs_5yswSV_df <- rbind(logObs_5yswSV_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_5yswSV_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_5yswSV)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_5yswSV[[i]])
logPred_5yswSV_df <- rbind(logPred_5yswSV_df, temp_df)
}
## Plot
logPred_5yswSV_df$age <- as.character(logPred_5yswSV_df$age)
logObs_5yswSV_df$age <- as.character(logObs_5yswSV_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_5yswSV_df[logObs_5yswSV_df$fleet == 1 & logObs_5yswSV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_5yswSV_df[logPred_5yswSV_df$fleet == 1 & logPred_5yswSV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.45: Five year sliding window model (survey variation) fit to catch data by age (values on the model link function scale).
# fitplot(fit_5yswSV, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_5yswSV_df[logObs_5yswSV_df$fleet == 2 & logObs_5yswSV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_5yswSV_df[logPred_5yswSV_df$fleet == 2 & logPred_5yswSV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.46: Five year sliding window model (survey variation) fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_5yswSV, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_5yswSV_df[logObs_5yswSV_df$fleet == 3 & logObs_5yswSV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_5yswSV_df[logPred_5yswSV_df$fleet == 3 & logPred_5yswSV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.47: Five year sliding window model (survey variation) fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_5yswSV, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_5yswSV <- residuals(fit_5yswSV)
resid_5yswSV_df <- data.frame(year = resid_5yswSV$year, fleet = resid_5yswSV$fleet,
age = resid_5yswSV$age, observation = resid_5yswSV$observation, mean = resid_5yswSV$mean,
residual = resid_5yswSV$residual)
resid_5yswSV_df$fleetName <- ifelse(resid_5yswSV_df$fleet == 1, attributes(resid_5yswSV)$fleetNames[1],
ifelse(resid_5yswSV_df$fleet == 2, attributes(resid_5yswSV)$fleetNames[2], ifelse(resid_5yswSV_df$fleet ==
3, attributes(resid_5yswSV)$fleetNames[3], NA)))
resid_5yswSV_df$fleetAltName <- ifelse(resid_5yswSV_df$fleet == 1, "Residual Fishing",
ifelse(resid_5yswSV_df$fleet == 2, "Q1 Surveys", ifelse(resid_5yswSV_df$fleet ==
3, "Q3/4 Surveys", NA)))
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. Q3/4 surveys), this is especially important.
if (!all(fit_5yswSV$conf$obsCorStruct == "ID")) {
corplot(fit_5yswSV)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 3.48: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the 5yswSV model.
ggplotly(ggplot(resid_5yswSV_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 3.49: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the five year sliding window model (with survey variation).
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
LO_5yswSV <- leaveout(fit_5yswSV)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_5yswSV$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_5yswSV$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_5yswSV$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_5yswSV$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(q3mat[(2002:DataYear) - 2001, ])),
SSB = q3mat[(2002:DataYear) - 2001, "SSB"], Fbar = q3mat[(2002:DataYear) - 2001,
"Fbar(3-5)"], R_age1 = q3mat[(2002:DataYear) - 2001, "R(age 1)"], CatchEst = catchtable(LO_5yswSV$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_5yswSV$series <- rep("full", times = nrow(asum_5yswSV))
asumi_5yswSV <- asum_5yswSV[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_5yswSV <- rbind(woq1, woq34, asumi_5yswSV)
# ssbplot(LO)
lossb_5yswSV <- ggplotly(ggplot() + geom_line(data = losum_5yswSV, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_5yswSV <- ggplotly(ggplot() + geom_line(data = losum_5yswSV, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_5yswSV <- ggplotly(ggplot() + geom_line(data = losum_5yswSV, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_5yswSV <- ggplotly(ggplot() + geom_line(data = losum_5yswSV, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_5yswSV[asum_5yswSV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb_5yswSV, lof_5yswSV, lorec_5yswSV, loca_5yswSV, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 3.50: Leave-one-out re-fits for the 5yswSV model (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
In this case, our model fit included the within-assessement year Q1 survey index as a tuning fleet, therefore, we are doing the peels and calculating the Mohn’s rho on data truncated by one year, so that we aren’t judging the model fit with the year in which it doesn’t have catch data nor the second index.
RETRO_5yswSV <- retro(fit_5yswSV, year = 5)
rho_5yswSV <- mohn(RETRO_5yswSV, lag = 1)
## Make RETROs in better plotting format
ret_5yswSV_df <- asum_5yswSV[asum_5yswSV$Year != max(asum_5yswSV$Year), ]
for (i in 1:length(RETRO_5yswSV)) {
tsum <- as.data.frame(summary(RETRO_5yswSV[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_5yswSV[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_5yswSV[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_5yswSV_df <- rbind(ret_5yswSV_df, tsum)
}
ret_5yswSV_df$series <- factor(ret_5yswSV_df$series, levels = c("full", "1", "2",
"3", "4", "5"))
# ssbplot(RETRO)
retssb_5yswSV <- layout(ggplotly(ggplot() + geom_line(data = ret_5yswSV_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year !=
max(asum_5yswSV$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_5yswSV_df$SSBhigh) *
0.85, x = ((max(ret_5yswSV_df$Year) - min(ret_5yswSV_df$Year)) * 0.2) + min(ret_5yswSV_df$Year),
label = paste0("Mohn's Rho = ", round(rho_5yswSV[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_5yswSV <- layout(ggplotly(ggplot() + geom_line(data = ret_5yswSV_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year !=
max(asum_5yswSV$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_5yswSV_df$Fhigh) *
0.85, x = ((max(ret_5yswSV_df$Year) - min(ret_5yswSV_df$Year)) * 0.8) + min(ret_5yswSV_df$Year),
label = paste0("Mohn's Rho = ", round(rho_5yswSV[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_5yswSV <- layout(ggplotly(ggplot() + geom_line(data = ret_5yswSV_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year !=
max(asum_5yswSV$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_5yswSV_df$Rhigh) *
0.85, x = ((max(ret_5yswSV_df$Year) - min(ret_5yswSV_df$Year)) * 0.2) + min(ret_5yswSV_df$Year),
label = paste0("Mohn's Rho = ", round(rho_5yswSV[1], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_5yswSV <- ggplotly(ggplot() + geom_line(data = ret_5yswSV_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year !=
max(asum_5yswSV$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_5yswSV[asum_5yswSV$Year !=
max(asum_5yswSV$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_5yswSV, retf_5yswSV, retrec_5yswSV, retca_5yswSV,
nrows = 2, shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE,
traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 3.51: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising five year sliding window means for stock weight at age and maturity and survey indices’ coefficients of variaton.
As noted in the previous section: The use of a five year sliding window mean on the biological data (explicitly, the stock weights at age and the maturity ogive) decreases the estimates of SSB in later years but not by as much as the three year sliding window mean, nor the biopar model. Again, this decrease is probably due to the observed reduction in weights at age, but the smaller magnitude is due to the larger window mean, which responds slower to the changes. However, this longer window means that the variability in model estimates of SSB, F and Recruitment, remains much lower. Furthermore, this model does not account for different catches in different management areas, and therefore, advice generated from short-term forecasts will have to be manually split based on basic ratios.
Weighting the inclusion of the survey indices’ by their coefficients of variation, does not greatly impact the estimates (e.g. SSB in 2023 changes by 0%), nor the uncertainty (e.g. the uncertainty range around this same figure changes by 0%).
This assessment combines two approaches to utilise the internal smoothing that SAM can do for biological data, while not over-parameterising the model. That is we utilise the biiopar smoothing for the more monotonic and faster changing stock-weights-at-age data, while we retain the simpler sliding window mean for the maturity ogives.
In preparing the input files, the sliding window mean simply averages the preceding x-1 years’ observations with each subsequent year to reduce the effect of annual variability. In this case x = 5, for a five year sliding window.
The biopar option is configured to treat all ages independently.
In this iteration, we also extend the timeseries to include the within-assessment year Q1 survey data.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_bio_5y_SV = fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_survCor")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
Because we are using the BIOPAR option in SAM to smooth the input of annually varying stock weights at age, we can first investigate the result of this part of the fit.
## Extract SW Fits from the model object
sw_fits <- as.data.frame(exp(fit_bio_5y_SV$pl$logSW)[1:length(2002:(DataYear + 1)),
])
colnames(sw_fits) <- as.character(c(1:7))
sw_fits$year <- factor(as.character(c(2002:(DataYear + 1))), levels = as.character(c((2002:(DataYear +
1)))))
sw_fitsL <- reshape(data = sw_fits, varying = as.character(c(1:7)), v.names = "swFits",
timevar = "Age", times = as.character(c(1:7)), idvar = "year", direction = "long")
## Extract SDs from the fit
sw_SDs <- as.data.frame(exp(fit_bio_5y_SV$plsd$logSW)[1:length(2002:(DataYear + 1)),
])
colnames(sw_SDs) <- as.character(c(1:7))
sw_SDs$year <- factor(as.character(c(2002:(DataYear + 1))), levels = as.character(c((2002:(DataYear +
1)))))
sw_SDsL <- reshape(data = sw_SDs, varying = as.character(c(1:7)), v.names = "swSDs",
timevar = "Age", times = as.character(c(1:7)), idvar = "year", direction = "long")
## Extract SW observations from the model object
sw_obs <- as.data.frame(fit_bio_5y_SV$data$stockMeanWeight)[1:length(2002:(DataYear +
1)), ]
colnames(sw_obs) <- as.character(c(1:7))
sw_obs$year <- factor(as.character(c(2002:(DataYear + 1))), levels = as.character(c((2002:(DataYear +
1)))))
sw_obsL <- reshape(data = sw_obs, varying = as.character(c(1:7)), v.names = "swObs",
timevar = "Age", times = as.character(c(1:7)), idvar = "year", direction = "long")
sw_mod <- merge(sw_fitsL, sw_SDsL, by = c("year", "Age"), all = TRUE)
sw_mod <- merge(sw_mod, sw_obsL, by = c("year", "Age"), all = TRUE)
sw_mod$year <- as.numeric(as.character(sw_mod$year))
ggplotly(ggplot(data = sw_mod) + geom_line(mapping = aes(x = year, y = swFits, colour = Age)) +
geom_point(mapping = aes(x = year, y = swObs, colour = Age), shape = 21) + facet_wrap(. ~
Age) + theme_few() + theme(axis.text.x = element_text(angle = 45, hjust = 1,
vjust = 1)) + scale_x_continuous(breaks = function(x) pretty(x, n = 5)) + scale_color_manual(values = c(ebpal)) +
scale_fill_manual(values = c(ebpal)) + guides(colour = "none", fill = "none"))
Figure 3.52: Fits to the stock weights at age, from the model with biopar on stock weights at age and sliding window on maturity. Annual observations are points, model fits are solid lines.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_bio_5y_SV <- as.data.frame(summary(fit_bio_5y_SV))
asum_bio_5y_SV$Year <- as.integer(row.names(summary(fit_bio_5y_SV)))
ct_bio_5y_SV <- as.data.frame(catchtable(fit_bio_5y_SV, obs.show = TRUE))
ct_bio_5y_SV$Year <- as.integer(rownames(ct_bio_5y_SV))
rownames(ct_bio_5y_SV) <- NULL
ct_bio_5y_SV <- rbind(ct_bio_5y_SV, data.frame(Estimate = NA, Low = NA, High = NA,
sop.catch = NA, Year = as.integer(2024)))
asum_bio_5y_SV <- merge(x = asum_bio_5y_SV, y = ct_bio_5y_SV, by = "Year")
colnames(asum_bio_5y_SV) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, as in the biopar model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 3.53: Model with Biopar on SWA and five year SWM on MO (survey variation) estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the baseline assessment, described first in this document.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 3.54: Model with Biopar on SWA and five year SWM on MO (survey variation) annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the baseline assessment, described first in this document.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV,
mapping = aes(x = Year, y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) +
theme_clean())
Figure 3.55: Model with Biopar on SWA and five year SWM on MO (survey variation) annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the baseline assessment, described first in this document.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) +
geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 3.56: Model with Biopar on SWA and five year SWM on MO (survey variation) annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the baseline assessment, described first in this document.
In the below plots, we’ve swapped out the baseline comparison for the model that uses the 5 year sliding window method for both weights at age and maturity ogives.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = asum_5yswSV, mapping = aes(x = Year, y = SSB),
colour = ebpal[2]) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV,
mapping = aes(x = Year, y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) +
theme_clean()
ggplotly(pssb)
Figure 3.57: Biopar and five year sliding window model (survey variation) estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the assessment using five year sliding windows for both stock weights at age and maturity ogives.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = asum_5yswSV, mapping = aes(x = Year, y = Fbar),
colour = ebpal[2]) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV,
mapping = aes(x = Year, y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) +
theme_clean())
Figure 3.58: Biopar and five year sliding window model (survey variation) annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the assessment using five year sliding windows for both stock weights at age and maturity ogives.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = asum_5yswSV, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = asum_5yswSV, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV,
mapping = aes(x = Year, y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) +
theme_clean())
Figure 3.59: Biopar and five year sliding window model (survey variation) annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the assessment using five year sliding windows for both stock weights at age and maturity ogives.
ggplotly(ggplot() + geom_line(data = asum_5yswSV[asum_5yswSV$Year != (DataYear +
1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[2]) + geom_ribbon(data = asum_5yswSV[asum_5yswSV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) +
geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 3.60: Biopar and five year sliding window model (survey variation) annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the assessment using five year sliding windows for both stock weights at age and maturity ogives.
### Data Prep
F_bio_5y_SV <- as.data.frame(faytable(fit_bio_5y_SV))
F_bio_5y_SV$Year <- rownames(F_bio_5y_SV)
rownames(F_bio_5y_SV) <- NULL
F_bio_5y_SVLng <- melt(data = F_bio_5y_SV, id.vars = "Year", variable.name = "Age",
value.name = "F_atAge")
F_bio_5y_SVLng$Year <- as.numeric(as.character(F_bio_5y_SVLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_bio_5y_SVLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_bio_5y_SVLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 3.61: Partial F at age; Model with Biopar on SWA and five year SWM on MO (survey variation)
sdplot(fit_bio_5y_SV, marg = c(5, 4, 1, 1))
Figure 3.62: Standard Devations by Fleet; base model
First we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_bio_5y_SV <- split(fit_bio_5y_SV$data$logobs, ceiling(seq_along(fit_bio_5y_SV$data$logobs)/fit_bio_5y_SV$data$maxAgePerFleet[1]))
logPred_bio_5y_SV <- split(fit_bio_5y_SV$rep$predObs, ceiling(seq_along(fit_bio_5y_SV$rep$predObs)/fit_bio_5y_SV$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_bio_5y_SV_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_bio_5y_SV)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_bio_5y_SV[[i]])
logObs_bio_5y_SV_df <- rbind(logObs_bio_5y_SV_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_bio_5y_SV_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_bio_5y_SV)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_bio_5y_SV[[i]])
logPred_bio_5y_SV_df <- rbind(logPred_bio_5y_SV_df, temp_df)
}
## Plot
logPred_bio_5y_SV_df$age <- as.character(logPred_bio_5y_SV_df$age)
logObs_bio_5y_SV_df$age <- as.character(logObs_bio_5y_SV_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_bio_5y_SV_df[logObs_bio_5y_SV_df$fleet == 1 & logObs_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_bio_5y_SV_df[logPred_bio_5y_SV_df$fleet == 1 & logPred_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.63: Model with Biopar on SWA and five year SWM on MO (survey variation) fit to catch data by age (values on the model link function scale).
# fitplot(fit_bio_5y_SV, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_bio_5y_SV_df[logObs_bio_5y_SV_df$fleet == 2 & logObs_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_bio_5y_SV_df[logPred_bio_5y_SV_df$fleet == 2 & logPred_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.64: Model with Biopar on SWA and five year SWM on MO (survey variation) fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_bio_5y_SV, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_bio_5y_SV_df[logObs_bio_5y_SV_df$fleet == 3 & logObs_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_bio_5y_SV_df[logPred_bio_5y_SV_df$fleet == 3 & logPred_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.65: Model with Biopar on SWA and five year SWM on MO (survey variation) fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_bio_5y_SV, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_bio_5y_SV <- residuals(fit_bio_5y_SV)
resid_bio_5y_SV_df <- data.frame(year = resid_bio_5y_SV$year, fleet = resid_bio_5y_SV$fleet,
age = resid_bio_5y_SV$age, observation = resid_bio_5y_SV$observation, mean = resid_bio_5y_SV$mean,
residual = resid_bio_5y_SV$residual)
resid_bio_5y_SV_df$fleetName <- ifelse(resid_bio_5y_SV_df$fleet == 1, attributes(resid_bio_5y_SV)$fleetNames[1],
ifelse(resid_bio_5y_SV_df$fleet == 2, attributes(resid_bio_5y_SV)$fleetNames[2],
ifelse(resid_bio_5y_SV_df$fleet == 3, attributes(resid_bio_5y_SV)$fleetNames[3],
NA)))
resid_bio_5y_SV_df$fleetAltName <- ifelse(resid_bio_5y_SV_df$fleet == 1, "Residual Fishing",
ifelse(resid_bio_5y_SV_df$fleet == 2, "Q1 Surveys", ifelse(resid_bio_5y_SV_df$fleet ==
3, "Q3/4 Surveys", NA)))
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. Q3/4 surveys), this is especially important.
if (!all(fit_bio_5y_SV$conf$obsCorStruct == "ID")) {
corplot(fit_bio_5y_SV)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 3.66: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the bio_5y_SV model.
ggplotly(ggplot(resid_bio_5y_SV_df) + geom_point(mapping = aes(x = year, y = age,
size = abs(residual), colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") +
scale_colour_manual(values = c(`TRUE` = ebpal[8], `FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1,
9), breaks = c(1:9)) + guides(colour = FALSE) + theme_few())
Figure 3.67: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with Biopar on SWA and five year SWM on MO (with survey variation).
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
LO_bio_5y_SV <- leaveout(fit_bio_5y_SV)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_bio_5y_SV$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_bio_5y_SV$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_bio_5y_SV$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_bio_5y_SV$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(q3mat[(2002:DataYear) - 2001, ])),
SSB = q3mat[(2002:DataYear) - 2001, "SSB"], Fbar = q3mat[(2002:DataYear) - 2001,
"Fbar(3-5)"], R_age1 = q3mat[(2002:DataYear) - 2001, "R(age 1)"], CatchEst = catchtable(LO_bio_5y_SV$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_bio_5y_SV$series <- rep("full", times = nrow(asum_bio_5y_SV))
asumi_bio_5y_SV <- asum_bio_5y_SV[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst",
"series")]
# =====
losum_bio_5y_SV <- rbind(woq1, woq34, asumi_bio_5y_SV)
# ssbplot(LO)
lossb_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb_bio_5y_SV, lof_bio_5y_SV, lorec_bio_5y_SV, loca_bio_5y_SV, nrows = 2,
shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 3.68: Leave-one-out re-fits for the Model with Biopar on SWA and five year SWM on MO (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
In this case, our model fit included the within-assessement year Q1 survey index as a tuning fleet, therefore, we are doing the peels and calculating the Mohn’s rho on data truncated by one year, so that we aren’t judging the model fit with the year in which it doesn’t have catch data nor the second index.
RETRO_bio_5y_SV <- retro(fit_bio_5y_SV, year = 5)
rho_bio_5y_SV <- mohn(RETRO_bio_5y_SV, lag = 1)
## Make RETROs in better plotting format
ret_bio_5y_SV_df <- asum_bio_5y_SV[asum_bio_5y_SV$Year != max(asum_bio_5y_SV$Year),
]
for (i in 1:length(RETRO_bio_5y_SV)) {
tsum <- as.data.frame(summary(RETRO_bio_5y_SV[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_bio_5y_SV[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_bio_5y_SV[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_bio_5y_SV_df <- rbind(ret_bio_5y_SV_df, tsum)
}
ret_bio_5y_SV_df$series <- factor(ret_bio_5y_SV_df$series, levels = c("full", "1",
"2", "3", "4", "5"))
# ssbplot(RETRO)
retssb_bio_5y_SV <- layout(ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_bio_5y_SV_df$SSBhigh) *
0.85, x = ((max(ret_bio_5y_SV_df$Year) - min(ret_bio_5y_SV_df$Year)) * 0.2) +
min(ret_bio_5y_SV_df$Year), label = paste0("Mohn's Rho = ", round(rho_bio_5y_SV[2],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_bio_5y_SV <- layout(ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_bio_5y_SV_df$Fhigh) *
0.85, x = ((max(ret_bio_5y_SV_df$Year) - min(ret_bio_5y_SV_df$Year)) * 0.8) +
min(ret_bio_5y_SV_df$Year), label = paste0("Mohn's Rho = ", round(rho_bio_5y_SV[3],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_bio_5y_SV <- layout(ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_bio_5y_SV_df$Rhigh) *
0.85, x = ((max(ret_bio_5y_SV_df$Year) - min(ret_bio_5y_SV_df$Year)) * 0.2) +
min(ret_bio_5y_SV_df$Year), label = paste0("Mohn's Rho = ", round(rho_bio_5y_SV[1],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3,
colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_bio_5y_SV, retf_bio_5y_SV, retrec_bio_5y_SV, retca_bio_5y_SV,
nrows = 2, shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE,
traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 3.69: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising Biopar for SWA and five year sliding window means for maturity and survey indices’ coefficients of variaton.
The combination of the biopar smoothing on the stock weights at age and the sliding window mean on the maturity ogives, results in more variability in SSB estimates, than the five year sliding window applied to both. The extension of the time series to include the Q1 survey from the in-assessment-year, leads to a higher SSB estimate in later years.
Large negative residuals in the most recent years, across most (sometimes all) ages and all fleets, are persistent across all options of model structure and model tuning. This suggests that there may be unexplained mortality in our model and that we may need to adjust the fixed mortality values that we’ve inherited from the previous two independent plaice stocks.
In this model, we are utilising the ability to operate with different fleets in SAM, to differentiate between the fishing that occurs in different management areas that divide the biological stock.
In the case of plaice in the Baltic, the Kattegat (SD21) is included in the same biological stock as teh rest of the Baltic Sea. Therefore, in this section we will split catch data by management area and treat them as independent fleets.
The model displayed here is only one iteration. We have not invested time in this option, as we saw some issues coming from: - The low sample sizes for the NS management area as a fleet. - Operationalising a model of this structure with the resources we have. - we do not have project funding to develop this assessment. - Tuning this model appropriately will take much more time and specific people to help find the right combinations of configurations across the various fleets and the degree to which to divide biological data between the fleets. - All of the follow-on work (reference points, short-term forecasts) would likely need a lot of tweaking and extra work to get them running.
In this example, we haven’t extended the time series to include the in-assessment-year Q1 survey data. We have included the survey indices coefficients of variation.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_NSBS = fitfromweb("ple.27.21-32_WKBPLAICE_2024_FleetNSBS_new2")
# =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_NSBS <- as.data.frame(summary(fit_NSBS))
asum_NSBS$Year <- as.integer(row.names(summary(fit_NSBS)))
ct_NSBS <- as.data.frame(catchtable(fit_NSBS, obs.show = TRUE))
ct_NSBS$Year <- as.integer(rownames(ct_NSBS))
rownames(ct_NSBS) <- NULL
ct_NSBS <- rbind(ct_NSBS, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_NSBS <- merge(x = asum_NSBS, y = ct_NSBS, by = "Year")
colnames(asum_NSBS) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, as in the biopar model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = SSB), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_NSBS, mapping = aes(x = Year,
y = SSB), colour = ebpal[1]) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 4.1: Model with management areas as fleets estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes). Purple is the current assessment, and green is the baseline assessment, described first in this document.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) +
geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_NSBS, mapping = aes(x = Year,
y = Fbar), colour = ebpal[1]) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 4.2: Model with management areas as fleets annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the current assessment, and green is the baseline assessment, described first in this document.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = R_age1),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Rlow,
ymax = Rhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_NSBS, mapping = aes(x = Year,
y = R_age1), colour = ebpal[1]) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 4.3: Model with management areas as fleets annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the current assessment, and green is the baseline assessment, described first in this document.
ggplotly(ggplot() + geom_line(data = basum, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[2]) + geom_ribbon(data = basum, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_NSBS[asum_NSBS$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) +
geom_ribbon(data = asum_NSBS[asum_NSBS$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_NSBS[asum_NSBS$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 4.4: Model with management areas as fleets annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the current assessment, and green is the baseline assessment, described first in this document.
### Data Prep
F_NSBS <- as.data.frame(faytable(fit_NSBS))
F_NSBS$Year <- rownames(F_NSBS)
rownames(F_NSBS) <- NULL
F_NSBSLng <- melt(data = F_NSBS, id.vars = "Year", variable.name = "Age", value.name = "F_atAge")
F_NSBSLng$Year <- as.numeric(as.character(F_NSBSLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_NSBSLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_NSBSLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 4.5: Partial F at age; Model with management areas as fleets
sdplot(fit_NSBS, marg = c(5, 4, 1, 1))
Figure 4.6: Standard Devations by Fleet; base model
First we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_NSBS <- split(fit_NSBS$data$logobs, ceiling(seq_along(fit_NSBS$data$logobs)/fit_NSBS$data$maxAgePerFleet[1]))
logPred_NSBS <- split(fit_NSBS$rep$predObs, ceiling(seq_along(fit_NSBS$rep$predObs)/fit_NSBS$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_NSBS_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_NSBS)) {
# Calculate the group and year
year <- ((i - 1) %/% 4) + 1
fleet <- ((i - 1) %% 4) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_NSBS[[i]])
logObs_NSBS_df <- rbind(logObs_NSBS_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_NSBS_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_NSBS)) {
# Calculate the group and year
year <- ((i - 1) %/% 4) + 1
fleet <- ((i - 1) %% 4) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_NSBS[[i]])
logPred_NSBS_df <- rbind(logPred_NSBS_df, temp_df)
}
## Plot
logPred_NSBS_df$age <- as.character(logPred_NSBS_df$age)
logObs_NSBS_df$age <- as.character(logObs_NSBS_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_NSBS_df[logObs_NSBS_df$fleet == 1 & logObs_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_NSBS_df[logPred_NSBS_df$fleet == 1 & logPred_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 4.7: Model with management areas as fleets fit to catch data from the North Sea management area (SD21) by age (values on the model link function scale).
# fitplot(fit_NSBS, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_NSBS_df[logObs_NSBS_df$fleet == 2 & logObs_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_NSBS_df[logPred_NSBS_df$fleet == 2 & logPred_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 4.8: Model with management areas as fleets fit to catch data from the Baltic Sea management area (SDs 21-32) by age (values on the model link function scale).
# fitplot(fit_NSBS, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_NSBS_df[logObs_NSBS_df$fleet == 3 & logObs_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_NSBS_df[logPred_NSBS_df$fleet == 3 & logPred_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 4.9: Model with management areas as fleets fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_NSBS, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_NSBS_df[logObs_NSBS_df$fleet == 4 & logObs_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_NSBS_df[logPred_NSBS_df$fleet == 4 & logPred_NSBS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 4.10: Model with management areas as fleets fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_NSBS, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_NSBS <- residuals(fit_NSBS)
resid_NSBS_df <- data.frame(year = resid_NSBS$year, fleet = resid_NSBS$fleet, age = resid_NSBS$age,
observation = resid_NSBS$observation, mean = resid_NSBS$mean, residual = resid_NSBS$residual)
resid_NSBS_df$fleetName <- ifelse(resid_NSBS_df$fleet == 1, attributes(resid_NSBS)$fleetNames[1],
ifelse(resid_NSBS_df$fleet == 2, attributes(resid_NSBS)$fleetNames[2], ifelse(resid_NSBS_df$fleet ==
3, attributes(resid_NSBS)$fleetNames[3], ifelse(resid_NSBS_df$fleet == 4,
attributes(resid_NSBS)$fleetNames[4], NA))))
resid_NSBS_df$fleetAltName <- ifelse(resid_NSBS_df$fleet == 1, "North Sea Fishing",
ifelse(resid_NSBS_df$fleet == 2, "Baltic Sea Fishing", ifelse(resid_NSBS_df$fleet ==
3, "Q1 Surveys", ifelse(resid_NSBS_df$fleet == 4, "Q3/4 Surveys", NA))))
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. No fleets match this criteria in this model configuration.), this is especially important.
if (!all(fit_NSBS$conf$obsCorStruct == "ID")) {
corplot(fit_NSBS)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
## [1] "No correlation structure configured for residuals across age."
ggplotly(ggplot(resid_NSBS_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 4.11: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the model with management areas as fleets.
A leave-one-out analysis is a form of sensitivity analysis, showing the impact the data from each tuning fleet has on the estimation of the key variables being estimated; namely SSB, F and recruitment.
First we must run the leave-one-out analysis which refits the model in two iterations, removing one survey at a time. Then we can plot each of these new model fits over the full model to see the impact the removal of each has.
LO_NSBS <- leaveout(fit_NSBS)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_NSBS$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_NSBS$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_NSBS$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_NSBS$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(q3mat[(2002:DataYear) - 2001, ])),
SSB = q3mat[(2002:DataYear) - 2001, "SSB"], Fbar = q3mat[(2002:DataYear) - 2001,
"Fbar(3-5)"], R_age1 = q3mat[(2002:DataYear) - 2001, "R(age 1)"], CatchEst = catchtable(LO_NSBS$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat))))
asum_NSBS$series <- rep("full", times = nrow(asum_NSBS))
asumi_NSBS <- asum_NSBS[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_NSBS <- rbind(woq1, woq34, asumi_NSBS)
# ssbplot(LO)
lossb_NSBS <- ggplotly(ggplot() + geom_line(data = losum_NSBS, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_NSBS <- ggplotly(ggplot() + geom_line(data = losum_NSBS, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_NSBS <- ggplotly(ggplot() + geom_line(data = losum_NSBS, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_NSBS <- ggplotly(ggplot() + geom_line(data = losum_NSBS, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_NSBS[asum_NSBS$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_NSBS[asum_NSBS$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb_NSBS, lof_NSBS, lorec_NSBS, loca_NSBS, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 4.12: Leave-one-out re-fits for the NSBS model (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
In this case, our model fit included the within-assessement year Q1 survey index as a tuning fleet, therefore, we are doing the peels and calculating the Mohn’s rho on data truncated by one year, so that we aren’t judging the model fit with the year in which it doesn’t have catch data nor the second index.
RETRO_NSBS <- retro(fit_NSBS, year = 5)
rho_NSBS <- mohn(RETRO_NSBS)
## Make RETROs in better plotting format
ret_NSBS_df <- asum_NSBS
for (i in 1:length(RETRO_NSBS)) {
tsum <- as.data.frame(summary(RETRO_NSBS[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_NSBS[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_NSBS[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
ret_NSBS_df <- rbind(ret_NSBS_df, tsum)
}
ret_NSBS_df$series <- factor(ret_NSBS_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_NSBS <- layout(ggplotly(ggplot() + geom_line(data = ret_NSBS_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + annotate(geom = "text",
y = max(ret_NSBS_df$SSBhigh) * 0.85, x = ((max(ret_NSBS_df$Year) - min(ret_NSBS_df$Year)) *
0.2) + min(ret_NSBS_df$Year), label = paste0("Mohn's Rho = ", round(rho_NSBS[2],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_NSBS <- layout(ggplotly(ggplot() + geom_line(data = ret_NSBS_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + annotate(geom = "text",
y = max(ret_NSBS_df$Fhigh) * 0.85, x = ((max(ret_NSBS_df$Year) - min(ret_NSBS_df$Year)) *
0.8) + min(ret_NSBS_df$Year), label = paste0("Mohn's Rho = ", round(rho_NSBS[3],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_NSBS <- layout(ggplotly(ggplot() + geom_line(data = ret_NSBS_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + annotate(geom = "text",
y = max(ret_NSBS_df$Rhigh) * 0.85, x = ((max(ret_NSBS_df$Year) - min(ret_NSBS_df$Year)) *
0.2) + min(ret_NSBS_df$Year), label = paste0("Mohn's Rho = ", round(rho_NSBS[1],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_NSBS <- ggplotly(ggplot() + geom_line(data = ret_NSBS_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_NSBS, mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(asum_NSBS,
mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) +
ylab("Catch (tonnes)") + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_NSBS, retf_NSBS, retrec_NSBS, retca_NSBS, nrows = 2,
shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 *
4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 4.13: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising management areas as fleets.
Based on discussion in the working group, the model that utilises SAM’s BioPar option for fitting to observed stock weights at age and a five year sliding window for the maturity ogives was selected to continue investigations.
However, like most other model runs, this model retains a lot of negative residuals in recent years, across all fleets, indicating that catches were below what was estimated for all fleets. This was concerning, as there appeared to be no basis in the data for the model estimations in the latest years. Based on this we decided to re-visit the inherited natural mortalities from the Ple.27.21-23 stock. These were time invariant and set to 0.2 for age 1 and 0.1 for all other ages.
In this section we change only the input values for the natural mortality and we do so based on three hypotheses: 1. That natural mortality is higher for younger ages and that it follows a general model elaborated by Gislason et al., 2010. a. this model is based on life-history parameters and size. b. this model is generally applicable across the whole time series. 2. That the shape of the decline in natural mortality with age from Gislason et al.’s model is correct, but that for the time invariant version of these natural mortalities, the absolute values may be different, due to un-calculated weighting in the numbers used to calculate the time-invariant mean. 3. That a Gislason model paramaterised by life-history traits collated across the whole time-series can be applied to annual length-at-age data to accurately produce time and age varying mortalities. a. these annually varying mortalities may be noisy due to noise in annual length samples, therefore a five-year sliding window mean is applied to them before input to SAM.
Already with the new mortalities we see reduced directionality in the residuals for all fleets in recent years. Fits to observed catches also look much better.
fit_nmG <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_SV_nmG")
## Extract data from model object
logObs_nmG <- split(fit_nmG$data$logobs, ceiling(seq_along(fit_nmG$data$logobs)/fit_nmG$data$maxAgePerFleet[1]))
logPred_nmG <- split(fit_nmG$rep$predObs, ceiling(seq_along(fit_nmG$rep$predObs)/fit_nmG$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmG[[i]])
logObs_nmG_df <- rbind(logObs_nmG_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmG[[i]])
logPred_nmG_df <- rbind(logPred_nmG_df, temp_df)
}
## Plot
logPred_nmG_df$age <- as.character(logPred_nmG_df$age)
logObs_nmG_df$age <- as.character(logObs_nmG_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_df[logObs_nmG_df$fleet == 1 & logObs_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_df[logPred_nmG_df$fleet == 1 & logPred_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.1: Model with time invariant Gislason natural mortalities fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_df[logObs_nmG_df$fleet == 2 & logObs_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_df[logPred_nmG_df$fleet == 2 & logPred_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.2: Model with time invariant Gislason natural mortalities fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_df[logObs_nmG_df$fleet == 3 & logObs_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_df[logPred_nmG_df$fleet == 3 & logPred_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.3: Model with time invariant Gislason natural mortalities fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG, fleets=3)
resid_nmG <- residuals(fit_nmG)
resid_nmG_df <- data.frame(year = resid_nmG$year, fleet = resid_nmG$fleet, age = resid_nmG$age,
observation = resid_nmG$observation, mean = resid_nmG$mean, residual = resid_nmG$residual)
resid_nmG_df$fleetName <- ifelse(resid_nmG_df$fleet == 1, attributes(resid_nmG)$fleetNames[1],
ifelse(resid_nmG_df$fleet == 2, attributes(resid_nmG)$fleetNames[2], ifelse(resid_nmG_df$fleet ==
3, attributes(resid_nmG)$fleetNames[3], NA)))
resid_nmG_df$fleetAltName <- ifelse(resid_nmG_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG_df$fleet == 3,
"Q3/4 Surveys", NA)))
if (!all(fit_nmG$conf$obsCorStruct == "ID")) {
corplot(fit_nmG)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 5.4: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the model with time invariant Gislason natural mortalities.
ggplotly(ggplot(resid_nmG_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 5.5: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the model with time invariant Gislason natural mortalities.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_nmG <- as.data.frame(summary(fit_nmG))
asum_nmG$Year <- as.integer(row.names(summary(fit_nmG)))
ct_nmG <- as.data.frame(catchtable(fit_nmG, obs.show = TRUE))
ct_nmG$Year <- as.integer(rownames(ct_nmG))
rownames(ct_nmG) <- NULL
ct_nmG <- rbind(ct_nmG, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG <- merge(x = asum_nmG, y = ct_nmG, by = "Year")
colnames(asum_nmG) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG$series <- rep("full", times = nrow(asum_nmG))
# =====
RETRO_nmG <- retro(fit_nmG, year = 5)
rho_nmG <- mohn(RETRO_nmG, lag = 1)
## Make RETROs in better plotting format
ret_nmG_df <- asum_nmG[asum_nmG$Year != max(asum_nmG$Year), ]
for (i in 1:length(RETRO_nmG)) {
tsum <- as.data.frame(summary(RETRO_nmG[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmG_df <- rbind(ret_nmG_df, tsum)
}
ret_nmG_df$series <- factor(ret_nmG_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_nmG <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year != max(asum_nmG$Year),
], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) +
annotate(geom = "text", y = max(ret_nmG_df$SSBhigh) * 0.85, x = ((max(ret_nmG_df$Year) -
min(ret_nmG_df$Year)) * 0.2) + min(ret_nmG_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmG[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmG <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year != max(asum_nmG$Year),
], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) +
annotate(geom = "text", y = max(ret_nmG_df$Fhigh) * 0.85, x = ((max(ret_nmG_df$Year) -
min(ret_nmG_df$Year)) * 0.8) + min(ret_nmG_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmG[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmG <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year !=
max(asum_nmG$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG_df$Rhigh) * 0.85, x = ((max(ret_nmG_df$Year) -
min(ret_nmG_df$Year)) * 0.2) + min(ret_nmG_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmG[1], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmG <- ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year !=
max(asum_nmG$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG[asum_nmG$Year != max(asum_nmG$Year),
], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) +
ylab("Catch (tonnes)") + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmG, retf_nmG, retrec_nmG, retca_nmG, nrows = 2, shareX = FALSE,
shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 5.6: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising time invariant, Gislason natural mortalities.
# === Scale natural mortality by a range of multipliers ====
profs = list()
nm = fit_nmG$data$natMor
mults = seq(0.1, 2, by = 0.1)
## Re-fit to scaled mortalities
for (mult in mults) {
fit_nmG$data$natMor = nm * mult
profs[[as.character(mult)]] = stockassessment:::refit(fit_nmG, silent = TRUE)
}
# ====
# === Aggregate AIC values and plot results ====
ll = sapply(profs, AIC)
plot(mults, ll)
Figure 5.7: Likelihood profile of scaled time invariant, Gislason natural mortalities. Y-axis = AIC, x-axis = scaling factor on nm.
# ====
We can now investigate if averaging the Gislason natural mortalities over the time series creates some scaling issues by ignoring weightings of numbers of observations contributing to each annual estimate. To partially address this we can investigate scaling these time invariant values up and down, to find where the best fit of this “shape” of natural mortalities across ages is.
By scaling the nm values up and down across a range of multipliers from 0.1 to 2, in increments of 0.1, we create 20 different runs from which we can compare AIC values. Here we find that a scaling value of 1.2 is responsible for the optimum model fit.
Therefore, below we investigate the model fit and residuals for a model with the Gislason, time invariant mortalities at age that have been scaled by a factor of 1.2 across all ages.
fit_nmGS <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG_S")
## Extract data from model object
logObs_nmGS <- split(fit_nmGS$data$logobs, ceiling(seq_along(fit_nmGS$data$logobs)/fit_nmGS$data$maxAgePerFleet[1]))
logPred_nmGS <- split(fit_nmGS$rep$predObs, ceiling(seq_along(fit_nmGS$rep$predObs)/fit_nmGS$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmGS_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmGS)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmGS[[i]])
logObs_nmGS_df <- rbind(logObs_nmGS_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmGS_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmGS)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmGS[[i]])
logPred_nmGS_df <- rbind(logPred_nmGS_df, temp_df)
}
## Plot
logPred_nmGS_df$age <- as.character(logPred_nmGS_df$age)
logObs_nmGS_df$age <- as.character(logObs_nmGS_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmGS_df[logObs_nmGS_df$fleet == 1 & logObs_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGS_df[logPred_nmGS_df$fleet == 1 & logPred_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.8: Model with SCALED time invariant Gislason natural mortalities fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmGS, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmGS_df[logObs_nmGS_df$fleet == 2 & logObs_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGS_df[logPred_nmGS_df$fleet == 2 & logPred_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.9: Model with SCALED time invariant Gislason natural mortalities fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmGS, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmGS_df[logObs_nmGS_df$fleet == 3 & logObs_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGS_df[logPred_nmGS_df$fleet == 3 & logPred_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.10: Model with SCALED time invariant Gislason natural mortalities fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmGS, fleets=3)
resid_nmGS <- residuals(fit_nmGS)
resid_nmGS_df <- data.frame(year = resid_nmGS$year, fleet = resid_nmGS$fleet, age = resid_nmGS$age,
observation = resid_nmGS$observation, mean = resid_nmGS$mean, residual = resid_nmGS$residual)
resid_nmGS_df$fleetName <- ifelse(resid_nmGS_df$fleet == 1, attributes(resid_nmGS)$fleetNames[1],
ifelse(resid_nmGS_df$fleet == 2, attributes(resid_nmGS)$fleetNames[2], ifelse(resid_nmGS_df$fleet ==
3, attributes(resid_nmGS)$fleetNames[3], NA)))
resid_nmGS_df$fleetAltName <- ifelse(resid_nmGS_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmGS_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmGS_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmGS$conf$obsCorStruct == "ID")) {
corplot(fit_nmGS)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 5.11: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with SCALED time invariant Gislason natural mortalities.
ggplotly(ggplot(resid_nmGS_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 5.12: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with SCALED time invariant Gislason natural mortalities.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_nmGS <- as.data.frame(summary(fit_nmGS))
asum_nmGS$Year <- as.integer(row.names(summary(fit_nmGS)))
ct_nmGS <- as.data.frame(catchtable(fit_nmGS, obs.show = TRUE))
ct_nmGS$Year <- as.integer(rownames(ct_nmGS))
rownames(ct_nmGS) <- NULL
ct_nmGS <- rbind(ct_nmGS, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmGS <- merge(x = asum_nmGS, y = ct_nmGS, by = "Year")
colnames(asum_nmGS) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmGS$series <- rep("full", times = nrow(asum_nmGS))
# =====
RETRO_nmGS <- retro(fit_nmGS, year = 5)
rho_nmGS <- mohn(RETRO_nmGS, lag = 1)
## Make RETROs in better plotting format
ret_nmGS_df <- asum_nmGS[asum_nmGS$Year != max(asum_nmGS$Year), ]
for (i in 1:length(RETRO_nmGS)) {
tsum <- as.data.frame(summary(RETRO_nmGS[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmGS[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmGS[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmGS_df <- rbind(ret_nmGS_df, tsum)
}
ret_nmGS_df$series <- factor(ret_nmGS_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_nmGS <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year != max(asum_nmGS$Year),
], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) +
annotate(geom = "text", y = max(ret_nmGS_df$SSBhigh) * 0.85, x = ((max(ret_nmGS_df$Year) -
min(ret_nmGS_df$Year)) * 0.2) + min(ret_nmGS_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmGS[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmGS <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_nmGS_df$Fhigh) * 0.85, x = ((max(ret_nmGS_df$Year) -
min(ret_nmGS_df$Year)) * 0.8) + min(ret_nmGS_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmGS[3], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmGS <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_nmGS_df$Rhigh) * 0.85, x = ((max(ret_nmGS_df$Year) -
min(ret_nmGS_df$Year)) * 0.2) + min(ret_nmGS_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmGS[1], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmGS <- ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmGS, retf_nmGS, retrec_nmGS, retca_nmGS, nrows = 2,
shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 *
4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 5.13: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising SCALED time invariant, Gislason natural mortalities.
An alternate way to dealing with the issue of taking the unweighted means over the time series is to not use the full time-series means. While the life-history parameters \(L_{\infty}\) and \(K\) are calculated from all ages observations (q1 surveys) for the whole time series, they can be applied in the Gislason model to specific years, by utilising the lengths at age for each year. Due to the sampling rates these data will be inherently noisey, just as is the case for the stock weights at age and the maturity ogives, hence we need some method of smoothing out this noise. In this example we utilise a five-year sliding window mean on the annual estimates of natural mortality at age, to allow mortality at age to vary in time without introducing too much interannual variation.
fit_nmG5ysw <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG5ysw")
## Extract data from model object
logObs_nmG5ysw <- split(fit_nmG5ysw$data$logobs, ceiling(seq_along(fit_nmG5ysw$data$logobs)/fit_nmG5ysw$data$maxAgePerFleet[1]))
logPred_nmG5ysw <- split(fit_nmG5ysw$rep$predObs, ceiling(seq_along(fit_nmG5ysw$rep$predObs)/fit_nmG5ysw$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmG5ysw[[i]])
logObs_nmG5ysw_df <- rbind(logObs_nmG5ysw_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmG5ysw[[i]])
logPred_nmG5ysw_df <- rbind(logPred_nmG5ysw_df, temp_df)
}
## Plot
logPred_nmG5ysw_df$age <- as.character(logPred_nmG5ysw_df$age)
logObs_nmG5ysw_df$age <- as.character(logObs_nmG5ysw_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5ysw_df[logObs_nmG5ysw_df$fleet == 1 & logObs_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5ysw_df[logPred_nmG5ysw_df$fleet == 1 & logPred_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.14: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5ysw, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5ysw_df[logObs_nmG5ysw_df$fleet == 2 & logObs_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5ysw_df[logPred_nmG5ysw_df$fleet == 2 & logPred_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.15: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5ysw, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5ysw_df[logObs_nmG5ysw_df$fleet == 3 & logObs_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5ysw_df[logPred_nmG5ysw_df$fleet == 3 & logPred_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 5.16: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5ysw, fleets=3)
resid_nmG5ysw <- residuals(fit_nmG5ysw)
resid_nmG5ysw_df <- data.frame(year = resid_nmG5ysw$year, fleet = resid_nmG5ysw$fleet,
age = resid_nmG5ysw$age, observation = resid_nmG5ysw$observation, mean = resid_nmG5ysw$mean,
residual = resid_nmG5ysw$residual)
resid_nmG5ysw_df$fleetName <- ifelse(resid_nmG5ysw_df$fleet == 1, attributes(resid_nmG5ysw)$fleetNames[1],
ifelse(resid_nmG5ysw_df$fleet == 2, attributes(resid_nmG5ysw)$fleetNames[2],
ifelse(resid_nmG5ysw_df$fleet == 3, attributes(resid_nmG5ysw)$fleetNames[3],
NA)))
resid_nmG5ysw_df$fleetAltName <- ifelse(resid_nmG5ysw_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG5ysw_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG5ysw_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5ysw$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5ysw)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 5.17: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities.
ggplotly(ggplot(resid_nmG5ysw_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 5.18: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_nmG5ysw <- as.data.frame(summary(fit_nmG5ysw))
asum_nmG5ysw$Year <- as.integer(row.names(summary(fit_nmG5ysw)))
ct_nmG5ysw <- as.data.frame(catchtable(fit_nmG5ysw, obs.show = TRUE))
ct_nmG5ysw$Year <- as.integer(rownames(ct_nmG5ysw))
rownames(ct_nmG5ysw) <- NULL
ct_nmG5ysw <- rbind(ct_nmG5ysw, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG5ysw <- merge(x = asum_nmG5ysw, y = ct_nmG5ysw, by = "Year")
colnames(asum_nmG5ysw) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5ysw$series <- rep("full", times = nrow(asum_nmG5ysw))
# =====
RETRO_nmG5ysw <- retro(fit_nmG5ysw, year = 5)
rho_nmG5ysw <- mohn(RETRO_nmG5ysw, lag = 1)
## Make RETROs in better plotting format
ret_nmG5ysw_df <- asum_nmG5ysw[asum_nmG5ysw$Year != max(asum_nmG5ysw$Year), ]
for (i in 1:length(RETRO_nmG5ysw)) {
tsum <- as.data.frame(summary(RETRO_nmG5ysw[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5ysw[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5ysw[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmG5ysw_df <- rbind(ret_nmG5ysw_df, tsum)
}
ret_nmG5ysw_df$series <- factor(ret_nmG5ysw_df$series, levels = c("full", "1", "2",
"3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5ysw_df$SSBhigh) *
0.85, x = ((max(ret_nmG5ysw_df$Year) - min(ret_nmG5ysw_df$Year)) * 0.2) + min(ret_nmG5ysw_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5ysw[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmG5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5ysw_df$Fhigh) *
0.85, x = ((max(ret_nmG5ysw_df$Year) - min(ret_nmG5ysw_df$Year)) * 0.8) + min(ret_nmG5ysw_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5ysw[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmG5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5ysw_df$Rhigh) *
0.85, x = ((max(ret_nmG5ysw_df$Year) - min(ret_nmG5ysw_df$Year)) * 0.2) + min(ret_nmG5ysw_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5ysw[1], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmG5ysw <- ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3,
colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmG5ysw, retf_nmG5ysw, retrec_nmG5ysw, retca_nmG5ysw,
nrows = 2, shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE,
traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 5.19: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising time variant (5 year sliding window mean), Gislason natural mortalities.
After interrogating the residuals and the model fits to observations, we can also compare the AICs of these models, as we’ve only modified the input data. We’ll also aggregate all of the Mohn’s rho values for easier comparisons.
nmAIC <- data.frame(`Type of M` = c("Time invariant, Gislason", "Optimally scaled, time invariant, Gislason",
"Time varying (5ysw), Gislason"), AIC = c(AIC(fit_nmG), AIC(fit_nmGS), AIC(fit_nmG5ysw)),
No.Params = c(length(coef(fit_nmG)), length(coef(fit_nmG)), length(coef(fit_nmG))),
Rho.SSB = c(rho_nmG[2], rho_nmGS[2], rho_nmG5ysw[2]), Rho.F = c(rho_nmG[3], rho_nmGS[3],
rho_nmG5ysw[3]), Rho.Rec = c(rho_nmG[1], rho_nmGS[1], rho_nmG5ysw[1]))
kable(nmAIC, caption = "Measures of model fit under various hypotheses of natural mortality.")
Type.of.M | AIC | No.Params | Rho.SSB | Rho.F | Rho.Rec |
---|---|---|---|---|---|
Time invariant, Gislason | 266.3056 | 36 | 0.0961543 | -0.0177057 | 0.2170079 |
Optimally scaled, time invariant, Gislason | 262.1331 | 36 | 0.0778651 | 0.0072884 | 0.1811353 |
Time varying (5ysw), Gislason | 262.8897 | 36 | 0.0625855 | -0.0177211 | 0.1195901 |
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = SSB,
series = "Historic nm"), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh, series = "Historic nm"),
fill = ebpal[1], alpha = 0.3) + geom_line(data = asum_nmG, mapping = aes(x = Year,
y = SSB, series = "Time-invariant Gislason nm"), colour = ebpal[2]) + geom_ribbon(data = asum_nmG,
mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh, series = "Time-invariant Gislason nm"),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_nmGS, mapping = aes(x = Year,
y = SSB, series = "Scaled, time-invariant Gislason nm"), colour = ebpal[3]) +
geom_ribbon(data = asum_nmGS, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh,
series = "Scaled, time-invariant Gislason nm"), fill = ebpal[3], alpha = 0.3) +
geom_line(data = asum_nmG5ysw, mapping = aes(x = Year, y = SSB, series = "Time varying (5ysw) Gislason nm"),
colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh, series = "Time varying (5ysw) Gislason nm"), fill = ebpal[4],
alpha = 0.3) + theme_clean()
ggplotly(pssb, tooltip = c("series", "x", "y"))
Figure 5.20: Estimated spawning stock biomass for ple.27.21-32 and 95% confidence intervals (tonnes). Purple is the model with the historic mortalities, green is with time invariant Gislason mortalities, orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = Fbar),
colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + geom_line(data = asum_nmG,
mapping = aes(x = Year, y = Fbar), colour = ebpal[2]) + geom_ribbon(data = asum_nmG,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = ebpal[2], alpha = 0.3) +
geom_line(data = asum_nmGS, mapping = aes(x = Year, y = Fbar), colour = ebpal[3]) +
geom_ribbon(data = asum_nmGS, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5ysw, mapping = aes(x = Year,
y = Fbar), colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[4], alpha = 0.3) + theme_clean())
Figure 5.21: Annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the model with the historic mortalities, green is with time invariant Gislason mortalities, orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = R_age1),
colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) + geom_line(data = asum_nmG,
mapping = aes(x = Year, y = R_age1), colour = ebpal[2]) + geom_ribbon(data = asum_nmG,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = ebpal[2], alpha = 0.3) +
geom_line(data = asum_nmGS, mapping = aes(x = Year, y = R_age1), colour = ebpal[3]) +
geom_ribbon(data = asum_nmGS, mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5ysw, mapping = aes(x = Year,
y = R_age1), colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[4], alpha = 0.3) + theme_clean())
Figure 5.22: Five year sliding window model annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the model with the historic mortalities, green is with time invariant Gislason mortalities, orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities.
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh), fill = ebpal[1], alpha = 0.3) + geom_line(data = asum_nmG,
mapping = aes(x = Year, y = CatchEst), colour = ebpal[2]) + geom_ribbon(data = asum_nmG,
mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh), fill = ebpal[2],
alpha = 0.3) + geom_line(data = asum_nmGS, mapping = aes(x = Year, y = CatchEst),
colour = ebpal[3]) + geom_ribbon(data = asum_nmGS, mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5ysw,
mapping = aes(x = Year, y = CatchEst), colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw,
mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh), fill = ebpal[4],
alpha = 0.3) + geom_point(data = asum_5ysw[asum_5ysw$Year != (DataYear + 1),
], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11], fill = ebpal[11]) +
ylab("Catch (tonnes)") + theme_clean())
Figure 5.23: Annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the model with the historic mortalities, green is with time invariant Gislason mortalities, orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities.
Based on these initial model runs we (during the benchmark meeting) have made the following conclusions:
The next working document will investigate the various options for modifying natural mortality, by interrogating model diagnostics further.
Subsequently, we will investigate incorporating discard survivals and perhaps German recreational harvests.