This working document is the third step in evaluating various assessment models for the newly combined Plaice stock, across the Sub Divisions 21 - 32. In the first assessment working document ( Ple.27.3a.21-32_Assessments ), Some of the more fundamental decisions about input data were evaluated and decided upon. In the second assessment working document ( Ple.27.3a.21-32_WKBPLAICE_NewMortality.RMD ), we focused in on an issue that is new to the stocks in this area, namely time-varying natural mortality. From this work we selected a model that uses Gislason et al (2010) natural mortalities calculated from life-history traits that were found from literature (via fishbase.se), and the varying length at age in the stock. Now, in this document, we investigate three different variations of input data to the selected model:
In scenario 1, we investigate the effect of setting maturity at age 1 to zero. In the maturity ogives calculated from the surveys, we see that ~25% of age1 fish are found to be mature. This conflicts with fundamental knowledge of European plaice biology, even in the specific context of this stock. Furhtermore, there are some issues around aging fish with some systematic errors reported for counting the first year’s ring. Hence some individuals that are sampled as being mature, may be incorrectly aged as age 1, when they are in fact age 2. Because of the large influxes of recruitment, having a significant proportion of recruits as mature leads to a volatile inflation of the SSB and a perception of the stock as larger than it really is. By manually setting maturity at age 1 to zero, we over-write any errors introduced during sampling with expert knowledge that age 1 fish are extremely unlikely to be mature, which reduces this over-estimation and interannual volatility of SSB.
In scenario 2, we run a type of sensitivity analysis but base the range of values in the sensitivity profile on the best available knowledge for discard survival. Discard survival is difficult to estimate for this species and this stock for multiple reasons. Survival varies according to a range of factors, only some of which we can account for and for only a subset of which do we have literature based observations of survival. For example, Survival after discards has shown to be correlated with:
In our data compilation process we can account for some of these to greater or lesser extents:
However, there are many more combinations of these factors than what is addressed in studies reported in primary or grey literature. Therefore, we’ve used the studies at face value, where they match some combinations of these factors, then used estimations for the gaps. The estimations are based on copying temporal trends from different fleet segments, where seasonal variation has been reported, and applying them to those segments where survival rates have only been reported for one season. Furthermore, because fleet segments in this stock are aggregated into “active” and “passive”, we take averages (unweighted) across different metiers where these are reported.
For more details on the evidence used and how it is applied, see another of the working documents from this workshop.
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.
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 catch data is available.
DataYear <- 2023
LastIcesAdvice <- 21735 # Advised catch for areas 21-32 (from stock splitting table in advice)
isFinal <- FALSE
# isFinal <- TRUE
The baseline model we’ve selected for this working document, comes from the investigations into the different versions of natural mortality we investigated in the previous working document: Ple.27.3a.21-32_WKBPLAICE_NewMortality.RMD.
This baseline model is called ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG5Fb and can be found on stockassessment.org.
fit_nmG5Fb <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG5Fb")
# === Create dataframe of current year's fit for plotting ====
asum_nmG5Fb <- as.data.frame(summary(fit_nmG5Fb))
asum_nmG5Fb$Year <- as.integer(row.names(summary(fit_nmG5Fb)))
ct_nmG5Fb <- as.data.frame(catchtable(fit_nmG5Fb, obs.show = TRUE))
ct_nmG5Fb$Year <- as.integer(rownames(ct_nmG5Fb))
rownames(ct_nmG5Fb) <- NULL
ct_nmG5Fb <- rbind(ct_nmG5Fb, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG5Fb <- merge(x = asum_nmG5Fb, y = ct_nmG5Fb, by = "Year")
colnames(asum_nmG5Fb) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
When the model is run we can easily see from warnings/errors if there is a convergence issue. However, we can also confirm this, explicitly:
Furthermore, SAM does not utilise “bounds” when fitting the model, and therefore, as standard check of whether model parameters are approaching their bounds is irrelevant.
## Extract data from model object
logObs_nmG5Fb <- split(fit_nmG5Fb$data$logobs, ceiling(seq_along(fit_nmG5Fb$data$logobs)/fit_nmG5Fb$data$maxAgePerFleet[1]))
logPred_nmG5Fb <- split(fit_nmG5Fb$rep$predObs, ceiling(seq_along(fit_nmG5Fb$rep$predObs)/fit_nmG5Fb$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5Fb_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5Fb)) {
# 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_nmG5Fb[[i]])
logObs_nmG5Fb_df <- rbind(logObs_nmG5Fb_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5Fb_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5Fb)) {
# 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_nmG5Fb[[i]])
logPred_nmG5Fb_df <- rbind(logPred_nmG5Fb_df, temp_df)
}
## Plot
logPred_nmG5Fb_df$age <- as.character(logPred_nmG5Fb_df$age)
logObs_nmG5Fb_df$age <- as.character(logObs_nmG5Fb_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_df[logObs_nmG5Fb_df$fleet == 1 & logObs_nmG5Fb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_df[logPred_nmG5Fb_df$fleet == 1 & logPred_nmG5Fb_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 2.1: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities (FishBase) fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_df[logObs_nmG5Fb_df$fleet == 2 & logObs_nmG5Fb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_df[logPred_nmG5Fb_df$fleet == 2 & logPred_nmG5Fb_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 2.2: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities (FishBase) fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_df[logObs_nmG5Fb_df$fleet == 3 & logObs_nmG5Fb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_df[logPred_nmG5Fb_df$fleet == 3 & logPred_nmG5Fb_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 2.3: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities (FishBase) fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb, fleets=3)
sdplot(fit_nmG5Fb, marg = c(5, 4, 1, 1))
Figure 2.4: Standard Devations by Fleet; Gislason, time-varying mortality (FishBase)
resid_nmG5Fb <- residuals(fit_nmG5Fb)
resid_nmG5Fb_df <- data.frame(year = resid_nmG5Fb$year, fleet = resid_nmG5Fb$fleet,
age = resid_nmG5Fb$age, observation = resid_nmG5Fb$observation, mean = resid_nmG5Fb$mean,
residual = resid_nmG5Fb$residual)
resid_nmG5Fb_df$fleetName <- ifelse(resid_nmG5Fb_df$fleet == 1, attributes(resid_nmG5Fb)$fleetNames[1],
ifelse(resid_nmG5Fb_df$fleet == 2, attributes(resid_nmG5Fb)$fleetNames[2], ifelse(resid_nmG5Fb_df$fleet ==
3, attributes(resid_nmG5Fb)$fleetNames[3], NA)))
resid_nmG5Fb_df$fleetAltName <- ifelse(resid_nmG5Fb_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG5Fb_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG5Fb_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5Fb$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5Fb)
# 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.5: 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 (FishBase).
ggplotly(ggplot(resid_nmG5Fb_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.6: 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 (FishBase).
We can also test to see if the model is converging on some local minimum (i.e. it’s fitting to some solution close to initialising values that represents noise and not the global solution, that is the proper solution). To do this, we add random noise to the initial parameter values to see if the model will converge on a different minimum in it’s data-space. The easiest way to investigate this is to see if the model fits vary alot depending on the starting values:
jit_nmG5Fb <- jit(fit = fit_nmG5Fb)
mt <- as.data.frame(modeltable(jit_nmG5Fb))
mt$model <- rownames(mt)
rownames(mt) <- NULL
kable(x = mt, digits = 3, caption = "Measures of fit for a series of model refits with jitter applied to the input parameters.")
log(L) | #par | AIC | model |
---|---|---|---|
-95.16 | 36 | 262.32 | M1 |
-95.16 | 36 | 262.32 | M2 |
-95.16 | 36 | 262.32 | M3 |
-95.16 | 36 | 262.32 | M4 |
-95.16 | 36 | 262.32 | M5 |
-95.16 | 36 | 262.32 | M6 |
-95.16 | 36 | 262.32 | M7 |
-95.16 | 36 | 262.32 | M8 |
-95.16 | 36 | 262.32 | M9 |
-95.16 | 36 | 262.32 | M10 |
-95.16 | 36 | 262.32 | M11 |
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_nmG5Fb <- leaveout(fit_nmG5Fb)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG5Fb$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG5Fb$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG5Fb$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG5Fb$`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_nmG5Fb$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG5Fb$series <- rep("full", times = nrow(asum_nmG5Fb))
asumi_nmG5Fb <- asum_nmG5Fb[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_nmG5Fb <- rbind(woq1, woq34, asumi_nmG5Fb)
# ssbplot(LO)
lossb_nmG5Fb <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb, 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_nmG5Fb <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb, 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_nmG5Fb <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb, 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_nmG5Fb <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb[asum_nmG5Fb$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_nmG5Fb, lof_nmG5Fb, lorec_nmG5Fb, loca_nmG5Fb, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 2.7: Leave-one-out re-fits for the Time varying, Gislason natural mortalities calculated from life-history off of fishbase (without Q1 survey = blue, without Q3/4 survey = 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 +).
# === 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_nmG5Fb <- as.data.frame(summary(fit_nmG5Fb))
asum_nmG5Fb$Year <- as.integer(row.names(summary(fit_nmG5Fb)))
ct_nmG5Fb <- as.data.frame(catchtable(fit_nmG5Fb, obs.show = TRUE))
ct_nmG5Fb$Year <- as.integer(rownames(ct_nmG5Fb))
rownames(ct_nmG5Fb) <- NULL
ct_nmG5Fb <- rbind(ct_nmG5Fb, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG5Fb <- merge(x = asum_nmG5Fb, y = ct_nmG5Fb, by = "Year")
colnames(asum_nmG5Fb) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5Fb$series <- rep("full", times = nrow(asum_nmG5Fb))
# =====
RETRO_nmG5Fb <- retro(fit_nmG5Fb, year = 5)
rho_nmG5Fb <- mohn(RETRO_nmG5Fb, lag = 1)
## Make RETROs in better plotting format
ret_nmG5Fb_df <- asum_nmG5Fb[asum_nmG5Fb$Year != max(asum_nmG5Fb$Year), ]
for (i in 1:length(RETRO_nmG5Fb)) {
tsum <- as.data.frame(summary(RETRO_nmG5Fb[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5Fb[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5Fb[[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_nmG5Fb_df <- rbind(ret_nmG5Fb_df, tsum)
}
ret_nmG5Fb_df$series <- factor(ret_nmG5Fb_df$series, levels = c("full", "1", "2",
"3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5Fb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
max(asum_nmG5Fb$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5Fb_df$SSBhigh) *
0.85, x = ((max(ret_nmG5Fb_df$Year) - min(ret_nmG5Fb_df$Year)) * 0.2) + min(ret_nmG5Fb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5Fb[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_nmG5Fb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
max(asum_nmG5Fb$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_df$Fhigh) *
0.85, x = ((max(ret_nmG5Fb_df$Year) - min(ret_nmG5Fb_df$Year)) * 0.8) + min(ret_nmG5Fb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5Fb[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_nmG5Fb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
max(asum_nmG5Fb$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_df$Rhigh) *
0.85, x = ((max(ret_nmG5Fb_df$Year) - min(ret_nmG5Fb_df$Year)) * 0.2) + min(ret_nmG5Fb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5Fb[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_nmG5Fb <- ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
max(asum_nmG5Fb$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
max(asum_nmG5Fb$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_nmG5Fb, retf_nmG5Fb, retrec_nmG5Fb, retca_nmG5Fb,
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.8: 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 (FishBase).
This model is called ple.27.21-32_WKBPLAICE_2024_BP_nmG5Fb_mo1 and can be found on stockassessment.org.
fit_nmG5Fb_mo1 <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BP_nmG5Fb_mo1")
# === Create dataframe of current year's fit for plotting ====
asum_nmG5Fb_mo1 <- as.data.frame(summary(fit_nmG5Fb_mo1))
asum_nmG5Fb_mo1$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1)))
ct_nmG5Fb_mo1 <- as.data.frame(catchtable(fit_nmG5Fb_mo1, obs.show = TRUE))
ct_nmG5Fb_mo1$Year <- as.integer(rownames(ct_nmG5Fb_mo1))
rownames(ct_nmG5Fb_mo1) <- NULL
ct_nmG5Fb_mo1 <- rbind(ct_nmG5Fb_mo1, data.frame(Estimate = NA, Low = NA, High = NA,
sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1 <- merge(x = asum_nmG5Fb_mo1, y = ct_nmG5Fb_mo1, by = "Year")
colnames(asum_nmG5Fb_mo1) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
When the model is run we can easily see from warnings/errors if there is a convergence issue. However, we can also confirm this, explicitly:
Furthermore, SAM does not utilise “bounds” when fitting the model, and therefore, as standard check of whether model parameters are approaching their bounds is irrelevant.
## Extract data from model object
logObs_nmG5Fb_mo1 <- split(fit_nmG5Fb_mo1$data$logobs, ceiling(seq_along(fit_nmG5Fb_mo1$data$logobs)/fit_nmG5Fb_mo1$data$maxAgePerFleet[1]))
logPred_nmG5Fb_mo1 <- split(fit_nmG5Fb_mo1$rep$predObs, ceiling(seq_along(fit_nmG5Fb_mo1$rep$predObs)/fit_nmG5Fb_mo1$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5Fb_mo1_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5Fb_mo1)) {
# 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_nmG5Fb_mo1[[i]])
logObs_nmG5Fb_mo1_df <- rbind(logObs_nmG5Fb_mo1_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5Fb_mo1_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5Fb_mo1)) {
# 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_nmG5Fb_mo1[[i]])
logPred_nmG5Fb_mo1_df <- rbind(logPred_nmG5Fb_mo1_df, temp_df)
}
## Plot
logPred_nmG5Fb_mo1_df$age <- as.character(logPred_nmG5Fb_mo1_df$age)
logObs_nmG5Fb_mo1_df$age <- as.character(logObs_nmG5Fb_mo1_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_df[logObs_nmG5Fb_mo1_df$fleet == 1 & logObs_nmG5Fb_mo1_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_df[logPred_nmG5Fb_mo1_df$fleet == 1 & logPred_nmG5Fb_mo1_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.1: Model with Age1 maturity set to zero, fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_df[logObs_nmG5Fb_mo1_df$fleet == 2 & logObs_nmG5Fb_mo1_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_df[logPred_nmG5Fb_mo1_df$fleet == 2 & logPred_nmG5Fb_mo1_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.2: Model with Age1 maturity set to zero, fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_df[logObs_nmG5Fb_mo1_df$fleet == 3 & logObs_nmG5Fb_mo1_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_df[logPred_nmG5Fb_mo1_df$fleet == 3 & logPred_nmG5Fb_mo1_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.3: Model with Age1 maturity set to zero, fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1, fleets=3)
sdplot(fit_nmG5Fb_mo1, marg = c(5, 4, 1, 1))
Figure 3.4: Standard Devations by Fleet; Age1 maturity set to zero.
resid_nmG5Fb_mo1 <- residuals(fit_nmG5Fb_mo1)
resid_nmG5Fb_mo1_df <- data.frame(year = resid_nmG5Fb_mo1$year, fleet = resid_nmG5Fb_mo1$fleet,
age = resid_nmG5Fb_mo1$age, observation = resid_nmG5Fb_mo1$observation, mean = resid_nmG5Fb_mo1$mean,
residual = resid_nmG5Fb_mo1$residual)
resid_nmG5Fb_mo1_df$fleetName <- ifelse(resid_nmG5Fb_mo1_df$fleet == 1, attributes(resid_nmG5Fb_mo1)$fleetNames[1],
ifelse(resid_nmG5Fb_mo1_df$fleet == 2, attributes(resid_nmG5Fb_mo1)$fleetNames[2],
ifelse(resid_nmG5Fb_mo1_df$fleet == 3, attributes(resid_nmG5Fb_mo1)$fleetNames[3],
NA)))
resid_nmG5Fb_mo1_df$fleetAltName <- ifelse(resid_nmG5Fb_mo1_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG5Fb_mo1_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG5Fb_mo1_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5Fb_mo1$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5Fb_mo1)
# 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.5: 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 (FishBase).
ggplotly(ggplot(resid_nmG5Fb_mo1_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.6: 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 age1 maturities set to zero.
We can also test to see if the model is converging on some local minimum (i.e. it’s fitting to some solution close to initialising values that represents noise and not the global solution, that is the proper solution). To do this, we add random noise to the initial parameter values to see if the model will converge on a different minimum in it’s data-space. The easiest way to investigate this is to see if the model fits vary alot depending on the starting values:
jit_nmG5Fb_mo1 <- jit(fit = fit_nmG5Fb_mo1)
mt <- as.data.frame(modeltable(jit_nmG5Fb_mo1))
mt$model <- rownames(mt)
rownames(mt) <- NULL
kable(x = mt, digits = 3, caption = "Measures of fit for a series of model refits with jitter applied to the input parameters.")
log(L) | #par | AIC | model |
---|---|---|---|
-95.16 | 36 | 262.32 | M1 |
-95.16 | 36 | 262.32 | M2 |
-95.16 | 36 | 262.32 | M3 |
-95.16 | 36 | 262.32 | M4 |
-95.16 | 36 | 262.32 | M5 |
-95.16 | 36 | 262.32 | M6 |
-95.16 | 36 | 262.32 | M7 |
-95.16 | 36 | 262.32 | M8 |
-95.16 | 36 | 262.32 | M9 |
-95.16 | 36 | 262.32 | M10 |
-95.16 | 36 | 262.32 | M11 |
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_nmG5Fb_mo1 <- leaveout(fit_nmG5Fb_mo1)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG5Fb_mo1$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG5Fb_mo1$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG5Fb_mo1$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG5Fb_mo1$`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_nmG5Fb_mo1$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG5Fb_mo1$series <- rep("full", times = nrow(asum_nmG5Fb_mo1))
asumi_nmG5Fb_mo1 <- asum_nmG5Fb_mo1[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst",
"series")]
# =====
losum_nmG5Fb_mo1 <- rbind(woq1, woq34, asumi_nmG5Fb_mo1)
# ssbplot(LO)
lossb_nmG5Fb_mo1 <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1, 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_nmG5Fb_mo1 <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1, 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_nmG5Fb_mo1 <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1, 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_nmG5Fb_mo1 <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$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_nmG5Fb_mo1, lof_nmG5Fb_mo1, lorec_nmG5Fb_mo1, loca_nmG5Fb_mo1,
nrows = 2, shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 3.7: Leave-one-out re-fits for the model with age1 maturities set to zero (without Q1 survey = blue, without Q3/4 survey = 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 +).
# === 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_nmG5Fb_mo1 <- as.data.frame(summary(fit_nmG5Fb_mo1))
asum_nmG5Fb_mo1$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1)))
ct_nmG5Fb_mo1 <- as.data.frame(catchtable(fit_nmG5Fb_mo1, obs.show = TRUE))
ct_nmG5Fb_mo1$Year <- as.integer(rownames(ct_nmG5Fb_mo1))
rownames(ct_nmG5Fb_mo1) <- NULL
ct_nmG5Fb_mo1 <- rbind(ct_nmG5Fb_mo1, data.frame(Estimate = NA, Low = NA, High = NA,
sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1 <- merge(x = asum_nmG5Fb_mo1, y = ct_nmG5Fb_mo1, by = "Year")
colnames(asum_nmG5Fb_mo1) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5Fb_mo1$series <- rep("full", times = nrow(asum_nmG5Fb_mo1))
# =====
RETRO_nmG5Fb_mo1 <- retro(fit_nmG5Fb_mo1, year = 5)
rho_nmG5Fb_mo1 <- mohn(RETRO_nmG5Fb_mo1, lag = 1)
## Make RETROs in better plotting format
ret_nmG5Fb_mo1_df <- asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year != max(asum_nmG5Fb_mo1$Year),
]
for (i in 1:length(RETRO_nmG5Fb_mo1)) {
tsum <- as.data.frame(summary(RETRO_nmG5Fb_mo1[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5Fb_mo1[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5Fb_mo1[[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_nmG5Fb_mo1_df <- rbind(ret_nmG5Fb_mo1_df, tsum)
}
ret_nmG5Fb_mo1_df$series <- factor(ret_nmG5Fb_mo1_df$series, levels = c("full", "1",
"2", "3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5Fb_mo1 <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_df,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
max(asum_nmG5Fb_mo1$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_df$SSBhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_df$Year) - min(ret_nmG5Fb_mo1_df$Year)) * 0.2) +
min(ret_nmG5Fb_mo1_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1[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_nmG5Fb_mo1 <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_df,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
max(asum_nmG5Fb_mo1$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_df$Fhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_df$Year) - min(ret_nmG5Fb_mo1_df$Year)) * 0.8) +
min(ret_nmG5Fb_mo1_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1[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_nmG5Fb_mo1 <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_df,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
max(asum_nmG5Fb_mo1$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_df$Rhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_df$Year) - min(ret_nmG5Fb_mo1_df$Year)) * 0.2) +
min(ret_nmG5Fb_mo1_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1[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_nmG5Fb_mo1 <- ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
max(asum_nmG5Fb_mo1$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
max(asum_nmG5Fb_mo1$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_nmG5Fb_mo1, retf_nmG5Fb_mo1, retrec_nmG5Fb_mo1, retca_nmG5Fb_mo1,
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.8: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model with age1 maturities set to zero.
ss_nmG5Fb_mo1 <- simstudy(fit = fit_nmG5Fb_mo1, nsim = 10)
This model is called ple.27.21-32_WKBPLAICE_2024_nmG5Fb_mo1_dslow and can be found on stockassessment.org.
fit_nmG5Fb_mo1_dslow <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_nmG5Fb_mo1_dslow")
# === Create dataframe of current year's fit for plotting ====
asum_nmG5Fb_mo1_dslow <- as.data.frame(summary(fit_nmG5Fb_mo1_dslow))
asum_nmG5Fb_mo1_dslow$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1_dslow)))
ct_nmG5Fb_mo1_dslow <- as.data.frame(catchtable(fit_nmG5Fb_mo1_dslow, obs.show = TRUE))
ct_nmG5Fb_mo1_dslow$Year <- as.integer(rownames(ct_nmG5Fb_mo1_dslow))
rownames(ct_nmG5Fb_mo1_dslow) <- NULL
ct_nmG5Fb_mo1_dslow <- rbind(ct_nmG5Fb_mo1_dslow, data.frame(Estimate = NA, Low = NA,
High = NA, sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1_dslow <- merge(x = asum_nmG5Fb_mo1_dslow, y = ct_nmG5Fb_mo1_dslow,
by = "Year")
colnames(asum_nmG5Fb_mo1_dslow) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
When the model is run we can easily see from warnings/errors if there is a convergence issue. However, we can also confirm this, explicitly:
Furthermore, SAM does not utilise “bounds” when fitting the model, and therefore, as standard check of whether model parameters are approaching their bounds is irrelevant.
## Extract data from model object
logObs_nmG5Fb_mo1_dslow <- split(fit_nmG5Fb_mo1_dslow$data$logobs, ceiling(seq_along(fit_nmG5Fb_mo1_dslow$data$logobs)/fit_nmG5Fb_mo1_dslow$data$maxAgePerFleet[1]))
logPred_nmG5Fb_mo1_dslow <- split(fit_nmG5Fb_mo1_dslow$rep$predObs, ceiling(seq_along(fit_nmG5Fb_mo1_dslow$rep$predObs)/fit_nmG5Fb_mo1_dslow$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5Fb_mo1_dslow_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5Fb_mo1_dslow)) {
# 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_nmG5Fb_mo1_dslow[[i]])
logObs_nmG5Fb_mo1_dslow_df <- rbind(logObs_nmG5Fb_mo1_dslow_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5Fb_mo1_dslow_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5Fb_mo1_dslow)) {
# 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_nmG5Fb_mo1_dslow[[i]])
logPred_nmG5Fb_mo1_dslow_df <- rbind(logPred_nmG5Fb_mo1_dslow_df, temp_df)
}
## Plot
logPred_nmG5Fb_mo1_dslow_df$age <- as.character(logPred_nmG5Fb_mo1_dslow_df$age)
logObs_nmG5Fb_mo1_dslow_df$age <- as.character(logObs_nmG5Fb_mo1_dslow_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dslow_df[logObs_nmG5Fb_mo1_dslow_df$fleet == 1 & logObs_nmG5Fb_mo1_dslow_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dslow_df[logPred_nmG5Fb_mo1_dslow_df$fleet == 1 & logPred_nmG5Fb_mo1_dslow_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.1: Model with Low discard survival estimates, fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dslow, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dslow_df[logObs_nmG5Fb_mo1_dslow_df$fleet == 2 & logObs_nmG5Fb_mo1_dslow_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dslow_df[logPred_nmG5Fb_mo1_dslow_df$fleet == 2 & logPred_nmG5Fb_mo1_dslow_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.2: Model with Low discard survival estimates, fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dslow, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dslow_df[logObs_nmG5Fb_mo1_dslow_df$fleet == 3 & logObs_nmG5Fb_mo1_dslow_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dslow_df[logPred_nmG5Fb_mo1_dslow_df$fleet == 3 & logPred_nmG5Fb_mo1_dslow_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.3: Model with low discard survival estimates, fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dslow, fleets=3)
sdplot(fit_nmG5Fb_mo1_dslow, marg = c(5, 4, 1, 1))
Figure 4.4: Standard Devations by Fleet; low discard survival estimates.
resid_nmG5Fb_mo1_dslow <- residuals(fit_nmG5Fb_mo1_dslow)
resid_nmG5Fb_mo1_dslow_df <- data.frame(year = resid_nmG5Fb_mo1_dslow$year, fleet = resid_nmG5Fb_mo1_dslow$fleet,
age = resid_nmG5Fb_mo1_dslow$age, observation = resid_nmG5Fb_mo1_dslow$observation,
mean = resid_nmG5Fb_mo1_dslow$mean, residual = resid_nmG5Fb_mo1_dslow$residual)
resid_nmG5Fb_mo1_dslow_df$fleetName <- ifelse(resid_nmG5Fb_mo1_dslow_df$fleet ==
1, attributes(resid_nmG5Fb_mo1_dslow)$fleetNames[1], ifelse(resid_nmG5Fb_mo1_dslow_df$fleet ==
2, attributes(resid_nmG5Fb_mo1_dslow)$fleetNames[2], ifelse(resid_nmG5Fb_mo1_dslow_df$fleet ==
3, attributes(resid_nmG5Fb_mo1_dslow)$fleetNames[3], NA)))
resid_nmG5Fb_mo1_dslow_df$fleetAltName <- ifelse(resid_nmG5Fb_mo1_dslow_df$fleet ==
1, "Residual Fishing", ifelse(resid_nmG5Fb_mo1_dslow_df$fleet == 2, "Q1 Surveys",
ifelse(resid_nmG5Fb_mo1_dslow_df$fleet == 3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5Fb_mo1_dslow$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5Fb_mo1_dslow)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 4.5: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with low discard survival estimates.
ggplotly(ggplot(resid_nmG5Fb_mo1_dslow_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.6: 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 low discard survival estimates.
We can also test to see if the model is converging on some local minimum (i.e. it’s fitting to some solution close to initialising values that represents noise and not the global solution, that is the proper solution). To do this, we add random noise to the initial parameter values to see if the model will converge on a different minimum in it’s data-space. The easiest way to investigate this is to see if the model fits vary alot depending on the starting values:
jit_nmG5Fb_mo1_dslow <- jit(fit = fit_nmG5Fb_mo1_dslow)
mt <- as.data.frame(modeltable(jit_nmG5Fb_mo1_dslow))
mt$model <- rownames(mt)
rownames(mt) <- NULL
kable(x = mt, digits = 3, caption = "Measures of fit for a series of model refits with jitter applied to the input parameters.")
log(L) | #par | AIC | model |
---|---|---|---|
-94.497 | 36 | 260.994 | M1 |
-94.497 | 36 | 260.994 | M2 |
-94.497 | 36 | 260.994 | M3 |
-94.497 | 36 | 260.994 | M4 |
-94.497 | 36 | 260.994 | M5 |
-94.497 | 36 | 260.994 | M6 |
-94.497 | 36 | 260.994 | M7 |
-94.497 | 36 | 260.994 | M8 |
-94.497 | 36 | 260.994 | M9 |
-94.497 | 36 | 260.994 | M10 |
-94.497 | 36 | 260.994 | M11 |
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_nmG5Fb_mo1_dslow <- leaveout(fit_nmG5Fb_mo1_dslow)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG5Fb_mo1_dslow$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG5Fb_mo1_dslow$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG5Fb_mo1_dslow$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG5Fb_mo1_dslow$`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_nmG5Fb_mo1_dslow$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG5Fb_mo1_dslow$series <- rep("full", times = nrow(asum_nmG5Fb_mo1_dslow))
asumi_nmG5Fb_mo1_dslow <- asum_nmG5Fb_mo1_dslow[, c("Year", "SSB", "Fbar", "R_age1",
"CatchEst", "series")]
# =====
losum_nmG5Fb_mo1_dslow <- rbind(woq1, woq34, asumi_nmG5Fb_mo1_dslow)
# ssbplot(LO)
lossb_nmG5Fb_mo1_dslow <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow,
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_nmG5Fb_mo1_dslow <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow,
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_nmG5Fb_mo1_dslow <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow,
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_nmG5Fb_mo1_dslow <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$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_nmG5Fb_mo1_dslow, lof_nmG5Fb_mo1_dslow, lorec_nmG5Fb_mo1_dslow,
loca_nmG5Fb_mo1_dslow, nrows = 2, shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 4.7: Leave-one-out re-fits for the model with age1 maturities set to zero (without Q1 survey = blue, without Q3/4 survey = 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 +).
# === 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_nmG5Fb_mo1_dslow <- as.data.frame(summary(fit_nmG5Fb_mo1_dslow))
asum_nmG5Fb_mo1_dslow$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1_dslow)))
ct_nmG5Fb_mo1_dslow <- as.data.frame(catchtable(fit_nmG5Fb_mo1_dslow, obs.show = TRUE))
ct_nmG5Fb_mo1_dslow$Year <- as.integer(rownames(ct_nmG5Fb_mo1_dslow))
rownames(ct_nmG5Fb_mo1_dslow) <- NULL
ct_nmG5Fb_mo1_dslow <- rbind(ct_nmG5Fb_mo1_dslow, data.frame(Estimate = NA, Low = NA,
High = NA, sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1_dslow <- merge(x = asum_nmG5Fb_mo1_dslow, y = ct_nmG5Fb_mo1_dslow,
by = "Year")
colnames(asum_nmG5Fb_mo1_dslow) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5Fb_mo1_dslow$series <- rep("full", times = nrow(asum_nmG5Fb_mo1_dslow))
# =====
RETRO_nmG5Fb_mo1_dslow <- retro(fit_nmG5Fb_mo1_dslow, year = 5)
rho_nmG5Fb_mo1_dslow <- mohn(RETRO_nmG5Fb_mo1_dslow, lag = 1)
## Make RETROs in better plotting format
ret_nmG5Fb_mo1_dslow_df <- asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year != max(asum_nmG5Fb_mo1_dslow$Year),
]
for (i in 1:length(RETRO_nmG5Fb_mo1_dslow)) {
tsum <- as.data.frame(summary(RETRO_nmG5Fb_mo1_dslow[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5Fb_mo1_dslow[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5Fb_mo1_dslow[[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_nmG5Fb_mo1_dslow_df <- rbind(ret_nmG5Fb_mo1_dslow_df, tsum)
}
ret_nmG5Fb_mo1_dslow_df$series <- factor(ret_nmG5Fb_mo1_dslow_df$series, levels = c("full",
"1", "2", "3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5Fb_mo1_dslow <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dslow_df,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
max(asum_nmG5Fb_mo1_dslow$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dslow_df$SSBhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dslow_df$Year) - min(ret_nmG5Fb_mo1_dslow_df$Year)) *
0.2) + min(ret_nmG5Fb_mo1_dslow_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dslow[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_nmG5Fb_mo1_dslow <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dslow_df,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
max(asum_nmG5Fb_mo1_dslow$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dslow_df$Fhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dslow_df$Year) - min(ret_nmG5Fb_mo1_dslow_df$Year)) *
0.8) + min(ret_nmG5Fb_mo1_dslow_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dslow[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_nmG5Fb_mo1_dslow <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dslow_df,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
max(asum_nmG5Fb_mo1_dslow$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dslow_df$Rhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dslow_df$Year) - min(ret_nmG5Fb_mo1_dslow_df$Year)) *
0.2) + min(ret_nmG5Fb_mo1_dslow_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dslow[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_nmG5Fb_mo1_dslow <- ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dslow_df,
mapping = aes(x = Year, y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
max(asum_nmG5Fb_mo1_dslow$Year), ], mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
max(asum_nmG5Fb_mo1_dslow$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_nmG5Fb_mo1_dslow, retf_nmG5Fb_mo1_dslow, retrec_nmG5Fb_mo1_dslow,
retca_nmG5Fb_mo1_dslow, 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.8: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model with low discard survival estimates.
This model is called ple.27.21-32_WKBPLAICE_2024_nmG5Fb_mo1_dsmed and can be found on stockassessment.org.
fit_nmG5Fb_mo1_dsmed <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_nmG5Fb_mo1_dsmed")
# === Create dataframe of current year's fit for plotting ====
asum_nmG5Fb_mo1_dsmed <- as.data.frame(summary(fit_nmG5Fb_mo1_dsmed))
asum_nmG5Fb_mo1_dsmed$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1_dsmed)))
ct_nmG5Fb_mo1_dsmed <- as.data.frame(catchtable(fit_nmG5Fb_mo1_dsmed, obs.show = TRUE))
ct_nmG5Fb_mo1_dsmed$Year <- as.integer(rownames(ct_nmG5Fb_mo1_dsmed))
rownames(ct_nmG5Fb_mo1_dsmed) <- NULL
ct_nmG5Fb_mo1_dsmed <- rbind(ct_nmG5Fb_mo1_dsmed, data.frame(Estimate = NA, Low = NA,
High = NA, sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1_dsmed <- merge(x = asum_nmG5Fb_mo1_dsmed, y = ct_nmG5Fb_mo1_dsmed,
by = "Year")
colnames(asum_nmG5Fb_mo1_dsmed) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
When the model is run we can easily see from warnings/errors if there is a convergence issue. However, we can also confirm this, explicitly:
Furthermore, SAM does not utilise “bounds” when fitting the model, and therefore, as standard check of whether model parameters are approaching their bounds is irrelevant.
## Extract data from model object
logObs_nmG5Fb_mo1_dsmed <- split(fit_nmG5Fb_mo1_dsmed$data$logobs, ceiling(seq_along(fit_nmG5Fb_mo1_dsmed$data$logobs)/fit_nmG5Fb_mo1_dsmed$data$maxAgePerFleet[1]))
logPred_nmG5Fb_mo1_dsmed <- split(fit_nmG5Fb_mo1_dsmed$rep$predObs, ceiling(seq_along(fit_nmG5Fb_mo1_dsmed$rep$predObs)/fit_nmG5Fb_mo1_dsmed$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5Fb_mo1_dsmed_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5Fb_mo1_dsmed)) {
# 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_nmG5Fb_mo1_dsmed[[i]])
logObs_nmG5Fb_mo1_dsmed_df <- rbind(logObs_nmG5Fb_mo1_dsmed_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5Fb_mo1_dsmed_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5Fb_mo1_dsmed)) {
# 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_nmG5Fb_mo1_dsmed[[i]])
logPred_nmG5Fb_mo1_dsmed_df <- rbind(logPred_nmG5Fb_mo1_dsmed_df, temp_df)
}
## Plot
logPred_nmG5Fb_mo1_dsmed_df$age <- as.character(logPred_nmG5Fb_mo1_dsmed_df$age)
logObs_nmG5Fb_mo1_dsmed_df$age <- as.character(logObs_nmG5Fb_mo1_dsmed_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dsmed_df[logObs_nmG5Fb_mo1_dsmed_df$fleet == 1 & logObs_nmG5Fb_mo1_dsmed_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dsmed_df[logPred_nmG5Fb_mo1_dsmed_df$fleet == 1 & logPred_nmG5Fb_mo1_dsmed_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 medium discard survival estimates, fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dsmed, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dsmed_df[logObs_nmG5Fb_mo1_dsmed_df$fleet == 2 & logObs_nmG5Fb_mo1_dsmed_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dsmed_df[logPred_nmG5Fb_mo1_dsmed_df$fleet == 2 & logPred_nmG5Fb_mo1_dsmed_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 medium discard survival estimates, fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dsmed, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dsmed_df[logObs_nmG5Fb_mo1_dsmed_df$fleet == 3 & logObs_nmG5Fb_mo1_dsmed_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dsmed_df[logPred_nmG5Fb_mo1_dsmed_df$fleet == 3 & logPred_nmG5Fb_mo1_dsmed_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.11: Model with medium discard survival estimates, fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dsmed, fleets=3)
sdplot(fit_nmG5Fb_mo1_dsmed, marg = c(5, 4, 1, 1))
Figure 4.12: Standard Devations by Fleet; medium discard survival estimates.
resid_nmG5Fb_mo1_dsmed <- residuals(fit_nmG5Fb_mo1_dsmed)
resid_nmG5Fb_mo1_dsmed_df <- data.frame(year = resid_nmG5Fb_mo1_dsmed$year, fleet = resid_nmG5Fb_mo1_dsmed$fleet,
age = resid_nmG5Fb_mo1_dsmed$age, observation = resid_nmG5Fb_mo1_dsmed$observation,
mean = resid_nmG5Fb_mo1_dsmed$mean, residual = resid_nmG5Fb_mo1_dsmed$residual)
resid_nmG5Fb_mo1_dsmed_df$fleetName <- ifelse(resid_nmG5Fb_mo1_dsmed_df$fleet ==
1, attributes(resid_nmG5Fb_mo1_dsmed)$fleetNames[1], ifelse(resid_nmG5Fb_mo1_dsmed_df$fleet ==
2, attributes(resid_nmG5Fb_mo1_dsmed)$fleetNames[2], ifelse(resid_nmG5Fb_mo1_dsmed_df$fleet ==
3, attributes(resid_nmG5Fb_mo1_dsmed)$fleetNames[3], NA)))
resid_nmG5Fb_mo1_dsmed_df$fleetAltName <- ifelse(resid_nmG5Fb_mo1_dsmed_df$fleet ==
1, "Residual Fishing", ifelse(resid_nmG5Fb_mo1_dsmed_df$fleet == 2, "Q1 Surveys",
ifelse(resid_nmG5Fb_mo1_dsmed_df$fleet == 3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5Fb_mo1_dsmed$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5Fb_mo1_dsmed)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 4.13: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with medium discard mortality.
ggplotly(ggplot(resid_nmG5Fb_mo1_dsmed_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.14: 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 medium discard mortality.
We can also test to see if the model is converging on some local minimum (i.e. it’s fitting to some solution close to initialising values that represents noise and not the global solution, that is the proper solution). To do this, we add random noise to the initial parameter values to see if the model will converge on a different minimum in it’s data-space. The easiest way to investigate this is to see if the model fits vary alot depending on the starting values:
jit_nmG5Fb_mo1_dsmed <- jit(fit = fit_nmG5Fb_mo1_dsmed)
mt <- as.data.frame(modeltable(jit_nmG5Fb_mo1_dsmed))
mt$model <- rownames(mt)
rownames(mt) <- NULL
kable(x = mt, digits = 3, caption = "Measures of fit for a series of model refits with jitter applied to the input parameters.")
log(L) | #par | AIC | model |
---|---|---|---|
-95.389 | 36 | 262.779 | M1 |
-95.389 | 36 | 262.779 | M2 |
-95.389 | 36 | 262.779 | M3 |
-95.389 | 36 | 262.779 | M4 |
-95.389 | 36 | 262.779 | M5 |
-95.389 | 36 | 262.779 | M6 |
-95.389 | 36 | 262.779 | M7 |
-95.389 | 36 | 262.779 | M8 |
-95.389 | 36 | 262.779 | M9 |
-95.389 | 36 | 262.779 | M10 |
-95.389 | 36 | 262.779 | M11 |
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_nmG5Fb_mo1_dsmed <- leaveout(fit_nmG5Fb_mo1_dsmed)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG5Fb_mo1_dsmed$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG5Fb_mo1_dsmed$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG5Fb_mo1_dsmed$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG5Fb_mo1_dsmed$`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_nmG5Fb_mo1_dsmed$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG5Fb_mo1_dsmed$series <- rep("full", times = nrow(asum_nmG5Fb_mo1_dsmed))
asumi_nmG5Fb_mo1_dsmed <- asum_nmG5Fb_mo1_dsmed[, c("Year", "SSB", "Fbar", "R_age1",
"CatchEst", "series")]
# =====
losum_nmG5Fb_mo1_dsmed <- rbind(woq1, woq34, asumi_nmG5Fb_mo1_dsmed)
# ssbplot(LO)
lossb_nmG5Fb_mo1_dsmed <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dsmed,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed,
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_nmG5Fb_mo1_dsmed <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dsmed,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed,
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_nmG5Fb_mo1_dsmed <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dsmed,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed,
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_nmG5Fb_mo1_dsmed <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dsmed,
mapping = aes(x = Year, y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$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_nmG5Fb_mo1_dsmed, lof_nmG5Fb_mo1_dsmed, lorec_nmG5Fb_mo1_dsmed,
loca_nmG5Fb_mo1_dsmed, nrows = 2, shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 4.15: Leave-one-out re-fits for the model with medium discard survival (without Q1 survey = blue, without Q3/4 survey = 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 +).
# === 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_nmG5Fb_mo1_dsmed <- as.data.frame(summary(fit_nmG5Fb_mo1_dsmed))
asum_nmG5Fb_mo1_dsmed$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1_dsmed)))
ct_nmG5Fb_mo1_dsmed <- as.data.frame(catchtable(fit_nmG5Fb_mo1_dsmed, obs.show = TRUE))
ct_nmG5Fb_mo1_dsmed$Year <- as.integer(rownames(ct_nmG5Fb_mo1_dsmed))
rownames(ct_nmG5Fb_mo1_dsmed) <- NULL
ct_nmG5Fb_mo1_dsmed <- rbind(ct_nmG5Fb_mo1_dsmed, data.frame(Estimate = NA, Low = NA,
High = NA, sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1_dsmed <- merge(x = asum_nmG5Fb_mo1_dsmed, y = ct_nmG5Fb_mo1_dsmed,
by = "Year")
colnames(asum_nmG5Fb_mo1_dsmed) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5Fb_mo1_dsmed$series <- rep("full", times = nrow(asum_nmG5Fb_mo1_dsmed))
# =====
RETRO_nmG5Fb_mo1_dsmed <- retro(fit_nmG5Fb_mo1_dsmed, year = 5)
rho_nmG5Fb_mo1_dsmed <- mohn(RETRO_nmG5Fb_mo1_dsmed, lag = 1)
## Make RETROs in better plotting format
ret_nmG5Fb_mo1_dsmed_df <- asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year != max(asum_nmG5Fb_mo1_dsmed$Year),
]
for (i in 1:length(RETRO_nmG5Fb_mo1_dsmed)) {
tsum <- as.data.frame(summary(RETRO_nmG5Fb_mo1_dsmed[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5Fb_mo1_dsmed[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5Fb_mo1_dsmed[[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_nmG5Fb_mo1_dsmed_df <- rbind(ret_nmG5Fb_mo1_dsmed_df, tsum)
}
ret_nmG5Fb_mo1_dsmed_df$series <- factor(ret_nmG5Fb_mo1_dsmed_df$series, levels = c("full",
"1", "2", "3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5Fb_mo1_dsmed <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dsmed_df,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
max(asum_nmG5Fb_mo1_dsmed$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dsmed_df$SSBhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dsmed_df$Year) - min(ret_nmG5Fb_mo1_dsmed_df$Year)) *
0.2) + min(ret_nmG5Fb_mo1_dsmed_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dsmed[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_nmG5Fb_mo1_dsmed <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dsmed_df,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
max(asum_nmG5Fb_mo1_dsmed$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dsmed_df$Fhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dsmed_df$Year) - min(ret_nmG5Fb_mo1_dsmed_df$Year)) *
0.8) + min(ret_nmG5Fb_mo1_dsmed_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dsmed[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_nmG5Fb_mo1_dsmed <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dsmed_df,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
max(asum_nmG5Fb_mo1_dsmed$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dsmed_df$Rhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dsmed_df$Year) - min(ret_nmG5Fb_mo1_dsmed_df$Year)) *
0.2) + min(ret_nmG5Fb_mo1_dsmed_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dsmed[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_nmG5Fb_mo1_dsmed <- ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dsmed_df,
mapping = aes(x = Year, y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
max(asum_nmG5Fb_mo1_dsmed$Year), ], mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
max(asum_nmG5Fb_mo1_dsmed$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_nmG5Fb_mo1_dsmed, retf_nmG5Fb_mo1_dsmed, retrec_nmG5Fb_mo1_dsmed,
retca_nmG5Fb_mo1_dsmed, 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.16: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model with medium discard survival.
This model is called ple.27.21-32_WKBPLAICE_2024_nmG5Fb_mo1_dshig and can be found on stockassessment.org.
fit_nmG5Fb_mo1_dshig <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_nmG5Fb_mo1_dshig")
# === Create dataframe of current year's fit for plotting ====
asum_nmG5Fb_mo1_dshig <- as.data.frame(summary(fit_nmG5Fb_mo1_dshig))
asum_nmG5Fb_mo1_dshig$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1_dshig)))
ct_nmG5Fb_mo1_dshig <- as.data.frame(catchtable(fit_nmG5Fb_mo1_dshig, obs.show = TRUE))
ct_nmG5Fb_mo1_dshig$Year <- as.integer(rownames(ct_nmG5Fb_mo1_dshig))
rownames(ct_nmG5Fb_mo1_dshig) <- NULL
ct_nmG5Fb_mo1_dshig <- rbind(ct_nmG5Fb_mo1_dshig, data.frame(Estimate = NA, Low = NA,
High = NA, sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1_dshig <- merge(x = asum_nmG5Fb_mo1_dshig, y = ct_nmG5Fb_mo1_dshig,
by = "Year")
colnames(asum_nmG5Fb_mo1_dshig) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
When the model is run we can easily see from warnings/errors if there is a convergence issue. However, we can also confirm this, explicitly:
Furthermore, SAM does not utilise “bounds” when fitting the model, and therefore, as standard check of whether model parameters are approaching their bounds is irrelevant.
## Extract data from model object
logObs_nmG5Fb_mo1_dshig <- split(fit_nmG5Fb_mo1_dshig$data$logobs, ceiling(seq_along(fit_nmG5Fb_mo1_dshig$data$logobs)/fit_nmG5Fb_mo1_dshig$data$maxAgePerFleet[1]))
logPred_nmG5Fb_mo1_dshig <- split(fit_nmG5Fb_mo1_dshig$rep$predObs, ceiling(seq_along(fit_nmG5Fb_mo1_dshig$rep$predObs)/fit_nmG5Fb_mo1_dshig$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5Fb_mo1_dshig_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5Fb_mo1_dshig)) {
# 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_nmG5Fb_mo1_dshig[[i]])
logObs_nmG5Fb_mo1_dshig_df <- rbind(logObs_nmG5Fb_mo1_dshig_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5Fb_mo1_dshig_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5Fb_mo1_dshig)) {
# 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_nmG5Fb_mo1_dshig[[i]])
logPred_nmG5Fb_mo1_dshig_df <- rbind(logPred_nmG5Fb_mo1_dshig_df, temp_df)
}
## Plot
logPred_nmG5Fb_mo1_dshig_df$age <- as.character(logPred_nmG5Fb_mo1_dshig_df$age)
logObs_nmG5Fb_mo1_dshig_df$age <- as.character(logObs_nmG5Fb_mo1_dshig_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dshig_df[logObs_nmG5Fb_mo1_dshig_df$fleet == 1 & logObs_nmG5Fb_mo1_dshig_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dshig_df[logPred_nmG5Fb_mo1_dshig_df$fleet == 1 & logPred_nmG5Fb_mo1_dshig_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.17: Model with high discard survival estimates, fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dshig, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dshig_df[logObs_nmG5Fb_mo1_dshig_df$fleet == 2 & logObs_nmG5Fb_mo1_dshig_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dshig_df[logPred_nmG5Fb_mo1_dshig_df$fleet == 2 & logPred_nmG5Fb_mo1_dshig_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.18: Model with high discard survival estimates, fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dshig, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5Fb_mo1_dshig_df[logObs_nmG5Fb_mo1_dshig_df$fleet == 3 & logObs_nmG5Fb_mo1_dshig_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5Fb_mo1_dshig_df[logPred_nmG5Fb_mo1_dshig_df$fleet == 3 & logPred_nmG5Fb_mo1_dshig_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.19: Model with high discard survival estimates, fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5Fb_mo1_dshig, fleets=3)
sdplot(fit_nmG5Fb_mo1_dshig, marg = c(5, 4, 1, 1))
Figure 4.20: Standard Devations by Fleet; high discard survival estimates.
resid_nmG5Fb_mo1_dshig <- residuals(fit_nmG5Fb_mo1_dshig)
resid_nmG5Fb_mo1_dshig_df <- data.frame(year = resid_nmG5Fb_mo1_dshig$year, fleet = resid_nmG5Fb_mo1_dshig$fleet,
age = resid_nmG5Fb_mo1_dshig$age, observation = resid_nmG5Fb_mo1_dshig$observation,
mean = resid_nmG5Fb_mo1_dshig$mean, residual = resid_nmG5Fb_mo1_dshig$residual)
resid_nmG5Fb_mo1_dshig_df$fleetName <- ifelse(resid_nmG5Fb_mo1_dshig_df$fleet ==
1, attributes(resid_nmG5Fb_mo1_dshig)$fleetNames[1], ifelse(resid_nmG5Fb_mo1_dshig_df$fleet ==
2, attributes(resid_nmG5Fb_mo1_dshig)$fleetNames[2], ifelse(resid_nmG5Fb_mo1_dshig_df$fleet ==
3, attributes(resid_nmG5Fb_mo1_dshig)$fleetNames[3], NA)))
resid_nmG5Fb_mo1_dshig_df$fleetAltName <- ifelse(resid_nmG5Fb_mo1_dshig_df$fleet ==
1, "Residual Fishing", ifelse(resid_nmG5Fb_mo1_dshig_df$fleet == 2, "Q1 Surveys",
ifelse(resid_nmG5Fb_mo1_dshig_df$fleet == 3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5Fb_mo1_dshig$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5Fb_mo1_dshig)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 4.21: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with high discard mortality.
ggplotly(ggplot(resid_nmG5Fb_mo1_dshig_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.22: 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 high discard mortality.
We can also test to see if the model is converging on some local minimum (i.e. it’s fitting to some solution close to initialising values that represents noise and not the global solution, that is the proper solution). To do this, we add random noise to the initial parameter values to see if the model will converge on a different minimum in it’s data-space. The easiest way to investigate this is to see if the model fits vary alot depending on the starting values:
jit_nmG5Fb_mo1_dshig <- jit(fit = fit_nmG5Fb_mo1_dshig)
mt <- as.data.frame(modeltable(jit_nmG5Fb_mo1_dshig))
mt$model <- rownames(mt)
rownames(mt) <- NULL
kable(x = mt, digits = 3, caption = "Measures of fit for a series of model refits with jitter applied to the input parameters.")
log(L) | #par | AIC | model |
---|---|---|---|
-96.835 | 36 | 265.67 | M1 |
-96.835 | 36 | 265.67 | M2 |
-96.835 | 36 | 265.67 | M3 |
-96.835 | 36 | 265.67 | M4 |
-96.835 | 36 | 265.67 | M5 |
-96.835 | 36 | 265.67 | M6 |
-96.835 | 36 | 265.67 | M7 |
-96.835 | 36 | 265.67 | M8 |
-96.835 | 36 | 265.67 | M9 |
-96.835 | 36 | 265.67 | M10 |
-96.835 | 36 | 265.67 | M11 |
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_nmG5Fb_mo1_dshig <- leaveout(fit_nmG5Fb_mo1_dshig)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG5Fb_mo1_dshig$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG5Fb_mo1_dshig$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG5Fb_mo1_dshig$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG5Fb_mo1_dshig$`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_nmG5Fb_mo1_dshig$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG5Fb_mo1_dshig$series <- rep("full", times = nrow(asum_nmG5Fb_mo1_dshig))
asumi_nmG5Fb_mo1_dshig <- asum_nmG5Fb_mo1_dshig[, c("Year", "SSB", "Fbar", "R_age1",
"CatchEst", "series")]
# =====
losum_nmG5Fb_mo1_dshig <- rbind(woq1, woq34, asumi_nmG5Fb_mo1_dshig)
# ssbplot(LO)
lossb_nmG5Fb_mo1_dshig <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig,
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_nmG5Fb_mo1_dshig <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig,
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_nmG5Fb_mo1_dshig <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig,
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_nmG5Fb_mo1_dshig <- ggplotly(ggplot() + geom_line(data = losum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$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_nmG5Fb_mo1_dshig, lof_nmG5Fb_mo1_dshig, lorec_nmG5Fb_mo1_dshig,
loca_nmG5Fb_mo1_dshig, nrows = 2, shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 4.23: Leave-one-out re-fits for the model with high discard survival (without Q1 survey = blue, without Q3/4 survey = 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 +).
# === 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_nmG5Fb_mo1_dshig <- as.data.frame(summary(fit_nmG5Fb_mo1_dshig))
asum_nmG5Fb_mo1_dshig$Year <- as.integer(row.names(summary(fit_nmG5Fb_mo1_dshig)))
ct_nmG5Fb_mo1_dshig <- as.data.frame(catchtable(fit_nmG5Fb_mo1_dshig, obs.show = TRUE))
ct_nmG5Fb_mo1_dshig$Year <- as.integer(rownames(ct_nmG5Fb_mo1_dshig))
rownames(ct_nmG5Fb_mo1_dshig) <- NULL
ct_nmG5Fb_mo1_dshig <- rbind(ct_nmG5Fb_mo1_dshig, data.frame(Estimate = NA, Low = NA,
High = NA, sop.catch = NA, Year = as.integer(2024)))
asum_nmG5Fb_mo1_dshig <- merge(x = asum_nmG5Fb_mo1_dshig, y = ct_nmG5Fb_mo1_dshig,
by = "Year")
colnames(asum_nmG5Fb_mo1_dshig) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5Fb_mo1_dshig$series <- rep("full", times = nrow(asum_nmG5Fb_mo1_dshig))
# =====
RETRO_nmG5Fb_mo1_dshig <- retro(fit_nmG5Fb_mo1_dshig, year = 5)
rho_nmG5Fb_mo1_dshig <- mohn(RETRO_nmG5Fb_mo1_dshig, lag = 1)
## Make RETROs in better plotting format
ret_nmG5Fb_mo1_dshig_df <- asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year != max(asum_nmG5Fb_mo1_dshig$Year),
]
for (i in 1:length(RETRO_nmG5Fb_mo1_dshig)) {
tsum <- as.data.frame(summary(RETRO_nmG5Fb_mo1_dshig[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5Fb_mo1_dshig[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5Fb_mo1_dshig[[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_nmG5Fb_mo1_dshig_df <- rbind(ret_nmG5Fb_mo1_dshig_df, tsum)
}
ret_nmG5Fb_mo1_dshig_df$series <- factor(ret_nmG5Fb_mo1_dshig_df$series, levels = c("full",
"1", "2", "3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5Fb_mo1_dshig <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dshig_df,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
max(asum_nmG5Fb_mo1_dshig$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dshig_df$SSBhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dshig_df$Year) - min(ret_nmG5Fb_mo1_dshig_df$Year)) *
0.2) + min(ret_nmG5Fb_mo1_dshig_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dshig[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_nmG5Fb_mo1_dshig <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dshig_df,
mapping = aes(x = Year, y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
max(asum_nmG5Fb_mo1_dshig$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dshig_df$Fhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dshig_df$Year) - min(ret_nmG5Fb_mo1_dshig_df$Year)) *
0.8) + min(ret_nmG5Fb_mo1_dshig_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dshig[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_nmG5Fb_mo1_dshig <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dshig_df,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
max(asum_nmG5Fb_mo1_dshig$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5Fb_mo1_dshig_df$Rhigh) *
0.85, x = ((max(ret_nmG5Fb_mo1_dshig_df$Year) - min(ret_nmG5Fb_mo1_dshig_df$Year)) *
0.2) + min(ret_nmG5Fb_mo1_dshig_df$Year), label = paste0("Mohn's Rho = ", round(rho_nmG5Fb_mo1_dshig[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_nmG5Fb_mo1_dshig <- ggplotly(ggplot() + geom_line(data = ret_nmG5Fb_mo1_dshig_df,
mapping = aes(x = Year, y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
max(asum_nmG5Fb_mo1_dshig$Year), ], mapping = aes(x = Year, ymin = Catchlow,
ymax = Catchhigh), fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
max(asum_nmG5Fb_mo1_dshig$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_nmG5Fb_mo1_dshig, retf_nmG5Fb_mo1_dshig, retrec_nmG5Fb_mo1_dshig,
retca_nmG5Fb_mo1_dshig, 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.24: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model with high discard survival.
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("Base Model", "Age1 Maturity to Zero", "Low Discard Survival",
"Medium Discard Survival", "High Discard Survival"), AIC = c(AIC(fit_nmG5Fb),
AIC(fit_nmG5Fb_mo1), AIC(fit_nmG5Fb_mo1_dslow), AIC(fit_nmG5Fb_mo1_dsmed), AIC(fit_nmG5Fb_mo1_dshig)),
No.Params = c(length(coef(fit_nmG5Fb)), length(coef(fit_nmG5Fb_mo1)), length(coef(fit_nmG5Fb_mo1_dslow)),
length(coef(fit_nmG5Fb_mo1_dsmed)), length(coef(fit_nmG5Fb_mo1_dshig))),
Rho.SSB = c(rho_nmG5Fb[2], rho_nmG5Fb_mo1[2], rho_nmG5Fb_mo1_dslow[2], rho_nmG5Fb_mo1_dsmed[2],
rho_nmG5Fb_mo1_dshig[2]), Rho.F = c(rho_nmG5Fb[3], rho_nmG5Fb_mo1[3], rho_nmG5Fb_mo1_dslow[3],
rho_nmG5Fb_mo1_dsmed[3], rho_nmG5Fb_mo1_dshig[3]), Rho.Rec = c(rho_nmG5Fb[1],
rho_nmG5Fb_mo1[1], rho_nmG5Fb_mo1_dslow[1], rho_nmG5Fb_mo1_dsmed[1], rho_nmG5Fb_mo1_dshig[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 |
---|---|---|---|---|---|
Base Model | 262.3204 | 36 | 0.0600907 | -0.0158922 | 0.1149631 |
Age1 Maturity to Zero | 262.3204 | 36 | 0.0185523 | -0.0158922 | 0.1149631 |
Low Discard Survival | 260.9943 | 36 | 0.0084699 | -0.0130867 | 0.0944167 |
Medium Discard Survival | 262.7788 | 36 | 0.0068493 | -0.0108810 | 0.0906564 |
High Discard Survival | 265.6705 | 36 | 0.0051730 | -0.0074690 | 0.0866411 |
Now that we’ve considered the model fits, diagnostics and tuning, we might as well have a look at how the models’ estimates compare.
We can see from the graphs that attempts to account for discard survival by reducing the catches proportionally, have the result of scaling the estimates of SSB. In this case, in the latest data year (), we see a -15% change from the model not accounting for discard survival to the scenario with the lowest levels of survival. When we compare to the model with the highest levels of survival, we see an even greater change of -21%
# 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_nmG5Fb, mapping = aes(x = Year, y = SSB,
series = "Baseline"), colour = ebpal[1]) + geom_ribbon(data = asum_nmG5Fb, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1,
mapping = aes(x = Year, y = SSB, series = "Maturity at age1 set to zero"), colour = ebpal[2]) +
geom_ribbon(data = asum_nmG5Fb_mo1, mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dslow, mapping = aes(x = Year,
y = SSB, series = "Low Discard Survival"), colour = ebpal[3]) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = ebpal[3], alpha = 0.3) +
geom_line(data = asum_nmG5Fb_mo1_dsmed, mapping = aes(x = Year, y = SSB, series = "Medium Discard Survival"),
colour = ebpal[4]) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[4], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, y = SSB, series = "High Discard Survival"), colour = ebpal[13]) +
geom_ribbon(data = asum_nmG5Fb_mo1_dshig, mapping = aes(x = Year, ymin = SSBlow,
ymax = SSBhigh), fill = ebpal[13], alpha = 0.3) + theme_clean()
ggplotly(pssb, tooltip = c("series", "x", "y"))
Figure 5.1: Estimated spawning stock biomass for ple.27.21-32 and 95% confidence intervals (tonnes). Purple is the current working document’s baseline model, green is with maturity at age1 set to zero, orange is with Low discard Survival, blue is with medium discard survival, and red is with highdiscard survival.
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_nmG5Fb, mapping = aes(x = Year, y = Fbar,
series = "Baseline"), colour = ebpal[1]) + geom_ribbon(data = asum_nmG5Fb, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh, series = "Baseline"), fill = ebpal[1], alpha = 0.3) +
geom_line(data = asum_nmG5Fb_mo1, mapping = aes(x = Year, y = Fbar, series = "Maturity at age1 set to zero"),
colour = ebpal[2]) + geom_ribbon(data = asum_nmG5Fb_mo1, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh, series = "Maturity at age1 set to zero"), fill = ebpal[2],
alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dslow, mapping = aes(x = Year,
y = Fbar, series = "Low Discard Survival"), colour = ebpal[3]) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh, series = "Low Discard Survival"),
fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dsmed, mapping = aes(x = Year,
y = Fbar, series = "Medium Discard Survival"), colour = ebpal[4]) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh, series = "Medium Discard Survival"),
fill = ebpal[4], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dshig, mapping = aes(x = Year,
y = Fbar, series = "High Discard Survival"), colour = ebpal[13]) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh, series = "High Discard Survival"),
fill = ebpal[13], alpha = 0.3) + theme_clean())
Figure 5.2: 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 working document’s baseline model, green is with maturity at age1 set to zero, orange is with Low discard Survival, blue is with medium discard survival, and red is with highdiscard survival.
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_nmG5Fb, mapping = aes(x = Year, y = R_age1,
series = "Baseline"), colour = ebpal[1]) + geom_ribbon(data = asum_nmG5Fb, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh, series = "Baseline"), fill = ebpal[1], alpha = 0.3) +
geom_line(data = asum_nmG5Fb_mo1, mapping = aes(x = Year, y = R_age1, series = "Maturity at age1 set to zero"),
colour = ebpal[2]) + geom_ribbon(data = asum_nmG5Fb_mo1, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh, series = "Maturity at age1 set to zero"), fill = ebpal[2],
alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dslow, mapping = aes(x = Year,
y = R_age1, series = "Low Discard Survival"), colour = ebpal[3]) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "Low Discard Survival"),
fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dsmed, mapping = aes(x = Year,
y = R_age1, series = "Medium Discard Survival"), colour = ebpal[4]) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "Medium Discard Survival"),
fill = ebpal[4], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dshig, mapping = aes(x = Year,
y = R_age1, series = "High Discard Survival"), colour = ebpal[13]) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "High Discard Survival"),
fill = ebpal[13], alpha = 0.3) + theme_clean())
Figure 5.3: 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 working document’s baseline model, green is with maturity at age1 set to zero, orange is with Low discard Survival, blue is with medium discard survival, and red is with highdiscard survival.
ggplotly(ggplot() + geom_line(data = asum_nmG5Fb[asum_nmG5Fb$Year != (DataYear +
1), ], mapping = aes(x = Year, y = CatchEst, series = "Baseline"), colour = ebpal[1]) +
geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh, series = "Baseline"), fill = ebpal[1],
alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year != (DataYear +
1), ], mapping = aes(x = Year, y = CatchEst, series = "Maturity at age1 set to zero"),
colour = ebpal[2]) + geom_ribbon(data = asum_nmG5Fb_mo1[asum_nmG5Fb_mo1$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh,
series = "Maturity at age1 set to zero"), fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "Low Discard Survival"),
colour = ebpal[3]) + geom_ribbon(data = asum_nmG5Fb_mo1_dslow[asum_nmG5Fb_mo1_dslow$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh,
series = "Low Discard Survival"), fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "Medium Discard Survival"),
colour = ebpal[4]) + geom_ribbon(data = asum_nmG5Fb_mo1_dsmed[asum_nmG5Fb_mo1_dsmed$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh,
series = "Medium Discard Survival"), fill = ebpal[4], alpha = 0.3) + geom_line(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "High Discard Survival"),
colour = ebpal[13]) + geom_ribbon(data = asum_nmG5Fb_mo1_dshig[asum_nmG5Fb_mo1_dshig$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh,
series = "High Discard Survival"), fill = ebpal[13], alpha = 0.3) + ylab("Catch (tonnes)") +
theme_clean())
Figure 5.4: 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 working document’s baseline model, green is with maturity at age1 set to zero, orange is with Low discard Survival, blue is with medium discard survival, and red is with highdiscard survival.
The setting of age 1 proportion mature to zero, overwrites proportions determined from surveys with an assumption based on expert knowledge. This would not, ordinarily, be best practice. However, due to the difficulty in getting reliable age reading rings from the first year of life, one can assume that some of the fish reported as mature in age 1, are from the age 2 cohort. Furthermore, juvenile and pre-adult plaice are known to distribute themselves across depth gradients according to their size. In this area, plaice are known to spawn and settle as juveniles over extended periods crossing seasons, which introduces length, growth and survival differentials into the various within-year, “settlement cohorts”. This means that the surveys that operate in deeper waters are more likely to catch individuals from the larger end of the size spectrum for any given age cohort, further biasing the proportion of mature fish at any given age. All experts in the workshop agreed that having significant proportions of age 1 plaice mature was unrealistic and that accepting these numbers (which are likely a result of sampling bias and observation error) was probably over-inflating our estimation of SSB.
Once set to zero, the model fit utilising these data does not change (AICs = 262.32, 262.32, respectively), only the SSB estimates scale up and down. The SSB estimates fall by -23% when setting age1 maturity to zero, and will likely reduce the inter-annual variability of the SSB estimates because of the larger variances around the estimations of numbers at age1.
In this section of the document we utilised the upper, lower, and middle values of estimates of discard survival from the literature, where studies existed. For catch segments where there was no data, we borrowed or extrapolated known values and trends (respectively) to ensure all segments had some estimate of survival. This provided us with three scenarios upon which we could run a type of sensitivity analysis, to visualise the impact that accounting for discard survival might have on both the model fit and the perception of the stock.
Here we see that utilising the any of the new catch scenarios has not substantial impact on the model fit, with all AIC values being within a range of five from one another, and all retro-plots / Mohn’s rho statistics being acceptable.
However, there are substantial impacts on our perception of the stock size, where the low end of the published survival estimates results in substantial change in the SSB (-15%), and recruitment (-16%).
When we compare the effects of selecting discard survival rates from across the ranges of those reported in the literature, we can see that the variation in our model estimates of SSB are not as great, relative to not accounting for any survival. This can be seen in the percentage difference between the low and the high estimates of discard, relative to the estimate made for the model with no survival for SSB (5% between low and high survival, vs -15% betwen low survival to no survival).
However, this difference is much reduced when comparing fishing pressure (6% between low and high survival, vs -8% between low survival to no survival).
Based on the fact that there are likely biases and some errors in the collection of data for age 1 fish, and the fact that changing the maturity at age1 has no effect on the model fit, but down-scales the SSB estimates and likely reduces interannual variation in the estimates, then we recommend retaining this expert opinion based “over-write”.
As was stated at the beginning of the benchmark process, the information that exists for rates of survival in discarded plaice is weak, sparse and has proven to be massively dependent on factors that cannot be accounted for in an assessment process (i.e. time exposed to air, being the largest driver of survival).
In order to address a specific request from the European Commission to investigate incorporating discard survival into the assessment, we have borrowed information across various fleet strata and applied these ad-hoc across our catch dataset. The resultant model fits show that there is no substantial difference in the quality of the model fits (based on AIC) but that our perception of the stock is somewhat scaled, and scaled differently according to the indicator use (e.g. see SSB vs F discussion above).
Subsequent to this exercise, we realise that our status quo assumption is 0% survival, which does not seem to reflect reality, and that while studies on this topic are sparse, information does exist that we should not ignore. Our attempts to utilise this knowledge involve rudimentary ratio-based borrowing of information to cover the catch segments as best as possible, but whether this approach is robust enough to utilise in the production of advice remains up for debate with peers and reviewers.
Therefore, there remains a decision to be made, to select between the assessment without accounting for discard survival (but with age 1 maturity set to zero), or the medium discard survival scenario?