This working document is the second step in evaluating various assessment models for the newly combined Plaice stock, across the Sub Divisions 21 - 32. In the first document ( Ple.27.3a.21-32_Assessments ), Some of the more fundamental decisions about input data were evaluated and decided upon. In this document, we focus in on an issue that is new to the stocks in this area, namely natural mortality. Because of the rapidly increasing density derived from consecutive years of high recruitment and low fishing pressure, we see decreasing stock weights at age, lengths at age and we anticipate increasing mortality across ages. This is especially true in the younger ages, where recent record breaking cohorts continue to top each other. In this document we investigate the different ways that these new estimations of mortality can be introduced to the model and the effect that this has on the model fit and predictive power.
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 the DataYear, which is the year immediately preceeding the assessment year or the year for which the most recent data is available.
DataYear <- 2023
LastIcesAdvice <- 21735 # Advised catch for areas 21-32 (from stock splitting table in advice)
isFinal <- FALSE
# isFinal <- TRUE
In this working document, the baseline model was selected from the immediately preceding working document, titled: Ple.27.21-32_WKBPLAICE_Assessments.Rmd.
This assessment combines two approaches to utilise the internal smoothing that SAM can do for biological data, while not over-parameterising the model. That is we utilise the biiopar smoothing for the more monotonic and faster changing stock-weights-at-age data, while we retain the simpler sliding window mean for the maturity ogives.
In preparing the input files, the sliding window mean simply averages the preceding x-1 years’ observations with each subsequent year to reduce the effect of annual variability. In this case x = 5, for a five year sliding window.
The biopar option is configured to treat all ages independently.
In this and all future iterations, we extended the time series to include the within-assessment year Q1 survey data.
To run the model, we must first download the input data and configuration files from their public repository. The data in these files has been presented in another working document and in plenary, but in these downloads it has been reformatted to fit with SAM, commonly known as the Lowestoft, CEFAS, or xsa format.
# === Read in data and assign appropriate names ====
fit_bio_5y_SV = fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_survCor")
# olco <- fit_bio_5y_SV$conf$keyStockWeightObsVar
# fit_bio_5y_SV$conf$keyStockWeightObsVar <- earlierKSOV
# fit_bio_5y_SV$conf$keyStockWeightObsVar <- olco
# olco <- fit_bio_5y_SV$conf$keyVarObs fit_bio_5y_SV$conf$keyVarObs <-
# earlierKVO fit_bio_5y_SV$conf$keyVarObs <- olco
# fit_bio_5y_SV <- stockassessment:::refit(fit_bio_5y_SV) =====
This section interrogates the model fit and generates tables and figures that describe the fit. In the background, these are saved to disk but also printed here for visibility in the rendered document.
Because we are using the BIOPAR option in SAM to smooth the input of annually varying stock weights at age, we can first investigate the result of this part of the fit.
## Extract SW Fits from the model object
sw_fits <- as.data.frame(exp(fit_bio_5y_SV$pl$logSW)[1:length(2002:(DataYear + 1)),
])
colnames(sw_fits) <- as.character(c(1:7))
sw_fits$year <- factor(as.character(c(2002:(DataYear + 1))), levels = as.character(c((2002:(DataYear +
1)))))
sw_fitsL <- reshape(data = sw_fits, varying = as.character(c(1:7)), v.names = "swFits",
timevar = "Age", times = as.character(c(1:7)), idvar = "year", direction = "long")
## Extract SDs from the fit
sw_SDs <- as.data.frame(exp(fit_bio_5y_SV$plsd$logSW)[1:length(2002:(DataYear + 1)),
])
colnames(sw_SDs) <- as.character(c(1:7))
sw_SDs$year <- factor(as.character(c(2002:(DataYear + 1))), levels = as.character(c((2002:(DataYear +
1)))))
sw_SDsL <- reshape(data = sw_SDs, varying = as.character(c(1:7)), v.names = "swSDs",
timevar = "Age", times = as.character(c(1:7)), idvar = "year", direction = "long")
## Extract SW observations from the model object
sw_obs <- as.data.frame(fit_bio_5y_SV$data$stockMeanWeight)[1:length(2002:(DataYear +
1)), ]
colnames(sw_obs) <- as.character(c(1:7))
sw_obs$year <- factor(as.character(c(2002:(DataYear + 1))), levels = as.character(c((2002:(DataYear +
1)))))
sw_obsL <- reshape(data = sw_obs, varying = as.character(c(1:7)), v.names = "swObs",
timevar = "Age", times = as.character(c(1:7)), idvar = "year", direction = "long")
sw_mod <- merge(sw_fitsL, sw_SDsL, by = c("year", "Age"), all = TRUE)
sw_mod <- merge(sw_mod, sw_obsL, by = c("year", "Age"), all = TRUE)
sw_mod$year <- as.numeric(as.character(sw_mod$year))
ggplotly(ggplot(data = sw_mod) + geom_line(mapping = aes(x = year, y = swFits, colour = Age)) +
geom_point(mapping = aes(x = year, y = swObs, colour = Age), shape = 21) + facet_wrap(. ~
Age, scales = "free_y") + theme_few() + theme(axis.text.x = element_text(angle = 45,
hjust = 1, vjust = 1)) + scale_x_continuous(breaks = function(x) pretty(x, n = 5)) +
scale_color_manual(values = c(ebpal)) + scale_fill_manual(values = c(ebpal)) +
guides(colour = "none", fill = "none"))
Figure 2.1: Fits to the stock weights at age, from the model with biopar on stock weights at age and sliding window on maturity. Annual observations are points, model fits are solid lines. Note the varying y axes.
The baseline assessment described first in this document, is now used as a comparison for the changes made in this version of the model.
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_bio_5y_SV <- as.data.frame(summary(fit_bio_5y_SV))
asum_bio_5y_SV$Year <- as.integer(row.names(summary(fit_bio_5y_SV)))
ct_bio_5y_SV <- as.data.frame(catchtable(fit_bio_5y_SV, obs.show = TRUE))
ct_bio_5y_SV$Year <- as.integer(rownames(ct_bio_5y_SV))
rownames(ct_bio_5y_SV) <- NULL
ct_bio_5y_SV <- rbind(ct_bio_5y_SV, data.frame(Estimate = NA, Low = NA, High = NA,
sop.catch = NA, Year = as.integer(2024)))
asum_bio_5y_SV <- merge(x = asum_bio_5y_SV, y = ct_bio_5y_SV, by = "Year")
colnames(asum_bio_5y_SV) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow",
"SSBhigh", "Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
We can see in this model, as in the biopar model, that the changes in stock weights at age are probably reflected in the lower estimates of SSB, albeit with a higher interannual variability than is the case for the baseline model, which uses a fixed stock weight at age for all years.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x =
# c('MSY-Btrigger', 'Bpa', 'Blim'),levels = c('MSY-Btrigger', 'Bpa', 'Blim')))
pssb <- ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = SSB),
colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = ebpal[1], alpha = 0.3) + theme_clean()
ggplotly(pssb)
Figure 2.2: Model with Biopar on SWA and five year SWM on MO (survey variation) estimated spawning stock biomass for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (tonnes)
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = Fbar),
colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 2.3: Model with Biopar on SWA and five year SWM on MO (survey variation) annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = R_age1),
colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = ebpal[1], alpha = 0.3) + theme_clean())
Figure 2.4: Model with Biopar on SWA and five year SWM on MO (survey variation) annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers).
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV[asum_bio_5y_SV$Year != (DataYear +
1), ], mapping = aes(x = Year, y = CatchEst), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = ebpal[1], alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[11],
fill = ebpal[11]) + ylab("Catch (tonnes)") + theme_clean())
Figure 2.5: Model with Biopar on SWA and five year SWM on MO (survey variation) annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points).
### Data Prep
F_bio_5y_SV <- as.data.frame(faytable(fit_bio_5y_SV))
F_bio_5y_SV$Year <- rownames(F_bio_5y_SV)
rownames(F_bio_5y_SV) <- NULL
F_bio_5y_SVLng <- melt(data = F_bio_5y_SV, id.vars = "Year", variable.name = "Age",
value.name = "F_atAge")
F_bio_5y_SVLng$Year <- as.numeric(as.character(F_bio_5y_SVLng$Year))
### Plot
ggplotly(ggplot() + geom_line(data = F_bio_5y_SVLng, mapping = aes(x = Year, y = F_atAge,
colour = Age)) + scale_color_manual(values = ebpal, ) + scale_x_continuous(breaks = unique(F_bio_5y_SVLng$Year)) +
theme_few() + theme(axis.text.x = element_text(angle = 45)))
Figure 2.6: Partial F at age; Model with Biopar on SWA and five year SWM on MO (survey variation)
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: - The final model gradient: 135 - That there is a positive definite hessian: TRUE
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.
We can investigate the variance around the different catch/stock components.
sdplot(fit_bio_5y_SV, marg = c(5, 4, 1, 1))
Figure 2.7: Standard Devations by Fleet; base model
Next we can look at how the model fits the observations per age in the different fleets.
## Extract data from model object
logObs_bio_5y_SV <- split(fit_bio_5y_SV$data$logobs, ceiling(seq_along(fit_bio_5y_SV$data$logobs)/fit_bio_5y_SV$data$maxAgePerFleet[1]))
logPred_bio_5y_SV <- split(fit_bio_5y_SV$rep$predObs, ceiling(seq_along(fit_bio_5y_SV$rep$predObs)/fit_bio_5y_SV$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_bio_5y_SV_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_bio_5y_SV)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_bio_5y_SV[[i]])
logObs_bio_5y_SV_df <- rbind(logObs_bio_5y_SV_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_bio_5y_SV_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_bio_5y_SV)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_bio_5y_SV[[i]])
logPred_bio_5y_SV_df <- rbind(logPred_bio_5y_SV_df, temp_df)
}
## Plot
logPred_bio_5y_SV_df$age <- as.character(logPred_bio_5y_SV_df$age)
logObs_bio_5y_SV_df$age <- as.character(logObs_bio_5y_SV_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_bio_5y_SV_df[logObs_bio_5y_SV_df$fleet == 1 & logObs_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_bio_5y_SV_df[logPred_bio_5y_SV_df$fleet == 1 & logPred_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 2.8: Model with Biopar on SWA and five year SWM on MO (survey variation) fit to catch data by age (values on the model link function scale).
# fitplot(fit_bio_5y_SV, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_bio_5y_SV_df[logObs_bio_5y_SV_df$fleet == 2 & logObs_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_bio_5y_SV_df[logPred_bio_5y_SV_df$fleet == 2 & logPred_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 2.9: Model with Biopar on SWA and five year SWM on MO (survey variation) fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_bio_5y_SV, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_bio_5y_SV_df[logObs_bio_5y_SV_df$fleet == 3 & logObs_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_bio_5y_SV_df[logPred_bio_5y_SV_df$fleet == 3 & logPred_bio_5y_SV_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 2.10: Model with Biopar on SWA and five year SWM on MO (survey variation) fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_bio_5y_SV, fleets=3)
Now we can also calculate the residuals of our observations from our model.
resid_bio_5y_SV <- residuals(fit_bio_5y_SV)
resid_bio_5y_SV_df <- data.frame(year = resid_bio_5y_SV$year, fleet = resid_bio_5y_SV$fleet,
age = resid_bio_5y_SV$age, observation = resid_bio_5y_SV$observation, mean = resid_bio_5y_SV$mean,
residual = resid_bio_5y_SV$residual)
resid_bio_5y_SV_df$fleetName <- ifelse(resid_bio_5y_SV_df$fleet == 1, attributes(resid_bio_5y_SV)$fleetNames[1],
ifelse(resid_bio_5y_SV_df$fleet == 2, attributes(resid_bio_5y_SV)$fleetNames[2],
ifelse(resid_bio_5y_SV_df$fleet == 3, attributes(resid_bio_5y_SV)$fleetNames[3],
NA)))
resid_bio_5y_SV_df$fleetAltName <- ifelse(resid_bio_5y_SV_df$fleet == 1, "Residual Fishing",
ifelse(resid_bio_5y_SV_df$fleet == 2, "Q1 Surveys", ifelse(resid_bio_5y_SV_df$fleet ==
3, "Q3/4 Surveys", NA)))
We can also look for patterns of correlation in the residuals, indicating unexplained trends. Ideally, there will be no large patterns, especially those weighted asymmetrically in the figures. For fleets where configuration allows for a correlation structure between ages, from year to year (i.e. Q3/4 surveys), this is especially important.
if (!all(fit_bio_5y_SV$conf$obsCorStruct == "ID")) {
corplot(fit_bio_5y_SV)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 2.11: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the bio_5y_SV model.
ggplotly(ggplot(resid_bio_5y_SV_df) + geom_point(mapping = aes(x = year, y = age,
size = abs(residual), colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") +
scale_colour_manual(values = c(`TRUE` = ebpal[8], `FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1,
9), breaks = c(1:9)) + guides(colour = FALSE) + theme_few())
Figure 2.12: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with Biopar on SWA and five year SWM on MO (with survey variation).
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_bio_5y_SV <- jit(fit = fit_bio_5y_SV)
mt <- as.data.frame(modeltable(jit_bio_5y_SV))
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 |
---|---|---|---|
-127.608 | 36 | 327.216 | M1 |
-127.608 | 36 | 327.216 | M2 |
-127.608 | 36 | 327.216 | M3 |
-127.608 | 36 | 327.216 | M4 |
-127.608 | 36 | 327.216 | M5 |
-127.608 | 36 | 327.216 | M6 |
-127.608 | 36 | 327.216 | M7 |
-127.608 | 36 | 327.216 | M8 |
-127.608 | 36 | 327.216 | M9 |
-127.608 | 36 | 327.216 | M10 |
-127.608 | 36 | 327.216 | 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_bio_5y_SV <- leaveout(fit_bio_5y_SV)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_bio_5y_SV$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_bio_5y_SV$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_bio_5y_SV$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_bio_5y_SV$`w.o. Q1IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q1", times = nrow(q1mat)))
woq34 <- data.frame(Year = as.integer(row.names(q3mat[(2002:DataYear) - 2001, ])),
SSB = q3mat[(2002:DataYear) - 2001, "SSB"], Fbar = q3mat[(2002:DataYear) - 2001,
"Fbar(3-5)"], R_age1 = q3mat[(2002:DataYear) - 2001, "R(age 1)"], CatchEst = catchtable(LO_bio_5y_SV$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_bio_5y_SV$series <- rep("full", times = nrow(asum_bio_5y_SV))
asumi_bio_5y_SV <- asum_bio_5y_SV[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst",
"series")]
# =====
losum_bio_5y_SV <- rbind(woq1, woq34, asumi_bio_5y_SV)
# ssbplot(LO)
lossb_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Fplot(LO)
lof_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# recplot(LO)
lorec_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_bio_5y_SV, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh), fill = "black", alpha = 0.3) + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
# Catch plot (LO)
loca_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = losum_bio_5y_SV, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[8:9])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
layout(subplot(lossb_bio_5y_SV, lof_bio_5y_SV, lorec_bio_5y_SV, loca_bio_5y_SV, nrows = 2,
shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 2.13: Leave-one-out re-fits for the Model with Biopar on SWA and five year SWM on MO (without Q1 = blue, without Q3/4 = purple), overlain with full model estimates (black line and grey ribbon) of SSB (top left), F (top right), Recruitment (bottom left), and catch (bottom right; observations as yellow +).
Another way to interrogate the model’s ability to fit to data is to change the input data. By dropping years sequentially from the latest year, backwards, and refitting the model to these data, we can visualise how this “tuned” model performs in different data situations.
In this case, our model fit included the within-assessement year Q1 survey index as a tuning fleet, therefore, we are doing the peels and calculating the Mohn’s rho on data truncated by one year, so that we aren’t judging the model fit with the year in which it doesn’t have catch data nor the second index.
RETRO_bio_5y_SV <- retro(fit_bio_5y_SV, year = 5)
rho_bio_5y_SV <- mohn(RETRO_bio_5y_SV, lag = 1)
## Make RETROs in better plotting format
ret_bio_5y_SV_df <- asum_bio_5y_SV[asum_bio_5y_SV$Year != max(asum_bio_5y_SV$Year),
]
ret_bio_5y_SV_df$series <- rep("full", nrow(ret_bio_5y_SV_df))
for (i in 1:length(RETRO_bio_5y_SV)) {
tsum <- as.data.frame(summary(RETRO_bio_5y_SV[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_bio_5y_SV[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_bio_5y_SV[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_bio_5y_SV_df <- rbind(ret_bio_5y_SV_df, tsum)
}
ret_bio_5y_SV_df$series <- factor(ret_bio_5y_SV_df$series, levels = c("full", "1",
"2", "3", "4", "5"))
# ssbplot(RETRO)
retssb_bio_5y_SV <- layout(ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df,
mapping = aes(x = Year, y = SSB, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_bio_5y_SV_df$SSBhigh) *
0.85, x = ((max(ret_bio_5y_SV_df$Year) - min(ret_bio_5y_SV_df$Year)) * 0.2) +
min(ret_bio_5y_SV_df$Year), label = paste0("Mohn's Rho = ", round(rho_bio_5y_SV[2],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_bio_5y_SV <- layout(ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_bio_5y_SV_df$Fhigh) *
0.85, x = ((max(ret_bio_5y_SV_df$Year) - min(ret_bio_5y_SV_df$Year)) * 0.8) +
min(ret_bio_5y_SV_df$Year), label = paste0("Mohn's Rho = ", round(rho_bio_5y_SV[3],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_bio_5y_SV <- layout(ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df,
mapping = aes(x = Year, y = R_age1, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_bio_5y_SV_df$Rhigh) *
0.85, x = ((max(ret_bio_5y_SV_df$Year) - min(ret_bio_5y_SV_df$Year)) * 0.2) +
min(ret_bio_5y_SV_df$Year), label = paste0("Mohn's Rho = ", round(rho_bio_5y_SV[1],
digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_bio_5y_SV <- ggplotly(ggplot() + geom_line(data = ret_bio_5y_SV_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_bio_5y_SV[asum_bio_5y_SV$Year !=
max(asum_bio_5y_SV$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3,
colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_bio_5y_SV, retf_bio_5y_SV, retrec_bio_5y_SV, retca_bio_5y_SV,
nrows = 2, shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE,
traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 2.14: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising Biopar for SWA and five year sliding window means for maturity and survey indices’ coefficients of variaton.
The combination of the biopar smoothing on the stock weights at age and the sliding window mean on the maturity ogives, results in more variability in SSB estimates, than the five year sliding window applied to both. The extension of the time series to include the Q1 survey from the in-assessment-year, leads to a higher SSB estimate in later years.
Large negative residuals in the most recent years, across most (sometimes all) ages and all fleets, are persistent across all options of model structure and model tuning. This suggests that there may be unexplained mortality in our model and that we may need to adjust the fixed natural mortality values that we’ve inherited from the previous two independent plaice stocks.
Based on discussion in the working group, the model that utilises SAM’s BioPar option for fitting to observed stock weights at age and a five year sliding window for the maturity ogives was selected to continue investigations.
However, like most other model runs, this model retains a lot of negative residuals in recent years, across all fleets, indicating that catches were below what was estimated for all fleets. This was concerning, as there appeared to be no basis in the data for the model estimations in the latest years. Based on this we decided to re-visit the inherited natural mortalities from the Ple.27.21-23 and Ple.27.24-32 stocks. These were time invariant and set to 0.2 for age 1, and 0.1 for all other ages.
In this section we change only the input values for the natural mortality and we do so based on three hypotheses: 1. That natural mortality is higher for younger ages and that it follows a general model elaborated by Gislason et al., 2010. a. this model is based on life-history parameters and size. b. this model is generally applicable across the whole time series. 2. That the shape of the decline in natural mortality with age from Gislason et al.’s model is correct, but that for the time invariant version of these natural mortalities, the absolute values may be different, due to un-calculated weighting in the numbers used to calculate the time-invariant mean. 3. That a Gislason model paramaterised by life-history traits collated across the whole time-series can be applied to annual length-at-age data to accurately produce time and age varying mortalities. a. these annually varying mortalities may be noisy due to noise in annual length samples, therefore a five-year sliding window mean is applied to them before input to SAM.
fit_nmG <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_SV_nmG")
# fit_nmG <- stockassessment:::refit(fit_nmG)
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_nmG <- as.data.frame(summary(fit_nmG))
asum_nmG$Year <- as.integer(row.names(summary(fit_nmG)))
ct_nmG <- as.data.frame(catchtable(fit_nmG, obs.show = TRUE))
ct_nmG$Year <- as.integer(rownames(ct_nmG))
rownames(ct_nmG) <- NULL
ct_nmG <- rbind(ct_nmG, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG <- merge(x = asum_nmG, y = ct_nmG, by = "Year")
colnames(asum_nmG) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
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: - The final model gradient: 131 - That there is a positive definite hessian: TRUE
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.
Already with the new mortalities we see reduced directionality in the residuals for all fleets in recent years. Fits to observed catches also look much better.
## Extract data from model object
logObs_nmG <- split(fit_nmG$data$logobs, ceiling(seq_along(fit_nmG$data$logobs)/fit_nmG$data$maxAgePerFleet[1]))
logPred_nmG <- split(fit_nmG$rep$predObs, ceiling(seq_along(fit_nmG$rep$predObs)/fit_nmG$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmG[[i]])
logObs_nmG_df <- rbind(logObs_nmG_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmG[[i]])
logPred_nmG_df <- rbind(logPred_nmG_df, temp_df)
}
## Plot
logPred_nmG_df$age <- as.character(logPred_nmG_df$age)
logObs_nmG_df$age <- as.character(logObs_nmG_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_df[logObs_nmG_df$fleet == 1 & logObs_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_df[logPred_nmG_df$fleet == 1 & logPred_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.1: Model with time invariant Gislason natural mortalities fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_df[logObs_nmG_df$fleet == 2 & logObs_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_df[logPred_nmG_df$fleet == 2 & logPred_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.2: Model with time invariant Gislason natural mortalities fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_df[logObs_nmG_df$fleet == 3 & logObs_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_df[logPred_nmG_df$fleet == 3 & logPred_nmG_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.3: Model with time invariant Gislason natural mortalities fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG, fleets=3)
sdplot(fit_nmG, marg = c(5, 4, 1, 1))
Figure 3.4: Standard Devations by Fleet; Gislason, time-invariant mortality
resid_nmG <- residuals(fit_nmG)
resid_nmG_df <- data.frame(year = resid_nmG$year, fleet = resid_nmG$fleet, age = resid_nmG$age,
observation = resid_nmG$observation, mean = resid_nmG$mean, residual = resid_nmG$residual)
resid_nmG_df$fleetName <- ifelse(resid_nmG_df$fleet == 1, attributes(resid_nmG)$fleetNames[1],
ifelse(resid_nmG_df$fleet == 2, attributes(resid_nmG)$fleetNames[2], ifelse(resid_nmG_df$fleet ==
3, attributes(resid_nmG)$fleetNames[3], NA)))
resid_nmG_df$fleetAltName <- ifelse(resid_nmG_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG_df$fleet == 3,
"Q3/4 Surveys", NA)))
if (!all(fit_nmG$conf$obsCorStruct == "ID")) {
corplot(fit_nmG)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 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 invariant Gislason natural mortalities.
ggplotly(ggplot(resid_nmG_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 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 time invariant Gislason natural mortalities.
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_nmG <- jit(fit = fit_nmG)
mt <- as.data.frame(modeltable(jit_nmG))
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 |
---|---|---|---|
-97.153 | 36 | 266.306 | M1 |
-97.153 | 36 | 266.306 | M2 |
-97.153 | 36 | 266.306 | M3 |
-97.153 | 36 | 266.306 | M4 |
-97.153 | 36 | 266.306 | M5 |
-97.153 | 36 | 266.306 | M6 |
-97.153 | 36 | 266.306 | M7 |
-97.153 | 36 | 266.306 | M8 |
-97.153 | 36 | 266.306 | M9 |
-97.153 | 36 | 266.306 | M10 |
-97.153 | 36 | 266.306 | 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_nmG <- leaveout(fit_nmG)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG$`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_nmG$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG$series <- rep("full", times = nrow(asum_nmG))
asumi_nmG <- asum_nmG[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_nmG <- rbind(woq1, woq34, asumi_nmG)
# ssbplot(LO)
lossb_nmG <- ggplotly(ggplot() + geom_line(data = losum_nmG, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG, 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_nmG <- ggplotly(ggplot() + geom_line(data = losum_nmG, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG, 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_nmG <- ggplotly(ggplot() + geom_line(data = losum_nmG, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG, 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_nmG <- ggplotly(ggplot() + geom_line(data = losum_nmG, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG[asum_nmG$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_nmG, lof_nmG, lorec_nmG, loca_nmG, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 3.7: Leave-one-out re-fits for the Time invariant, Gislason natural mortalities (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_nmG <- as.data.frame(summary(fit_nmG))
asum_nmG$Year <- as.integer(row.names(summary(fit_nmG)))
ct_nmG <- as.data.frame(catchtable(fit_nmG, obs.show = TRUE))
ct_nmG$Year <- as.integer(rownames(ct_nmG))
rownames(ct_nmG) <- NULL
ct_nmG <- rbind(ct_nmG, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG <- merge(x = asum_nmG, y = ct_nmG, by = "Year")
colnames(asum_nmG) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG$series <- rep("full", times = nrow(asum_nmG))
# =====
RETRO_nmG <- retro(fit_nmG, year = 5)
rho_nmG <- mohn(RETRO_nmG, lag = 1)
## Make RETROs in better plotting format
ret_nmG_df <- asum_nmG[asum_nmG$Year != max(asum_nmG$Year), ]
for (i in 1:length(RETRO_nmG)) {
tsum <- as.data.frame(summary(RETRO_nmG[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmG_df <- rbind(ret_nmG_df, tsum)
}
ret_nmG_df$series <- factor(ret_nmG_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_nmG <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year != max(asum_nmG$Year),
], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) +
annotate(geom = "text", y = max(ret_nmG_df$SSBhigh) * 0.85, x = ((max(ret_nmG_df$Year) -
min(ret_nmG_df$Year)) * 0.2) + min(ret_nmG_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmG[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmG <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year != max(asum_nmG$Year),
], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = "black", alpha = 0.3) +
annotate(geom = "text", y = max(ret_nmG_df$Fhigh) * 0.85, x = ((max(ret_nmG_df$Year) -
min(ret_nmG_df$Year)) * 0.8) + min(ret_nmG_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmG[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmG <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year !=
max(asum_nmG$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG_df$Rhigh) * 0.85, x = ((max(ret_nmG_df$Year) -
min(ret_nmG_df$Year)) * 0.2) + min(ret_nmG_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmG[1], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmG <- ggplotly(ggplot() + geom_line(data = ret_nmG_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG[asum_nmG$Year !=
max(asum_nmG$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG[asum_nmG$Year != max(asum_nmG$Year),
], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5], fill = ebpal[5]) +
ylab("Catch (tonnes)") + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmG, retf_nmG, retrec_nmG, retca_nmG, nrows = 2, shareX = FALSE,
shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 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 utilising time invariant, Gislason natural mortalities.
# === Scale natural mortality by a range of multipliers ====
profs = list()
nm = fit_nmG$data$natMor
mults = seq(0.1, 2, by = 0.1)
## Re-fit to scaled mortalities
for (mult in mults) {
fit_nmG$data$natMor = nm * mult
profs[[as.character(mult)]] = stockassessment:::refit(fit_nmG, silent = TRUE)
}
# ====
# === Aggregate AIC values and plot results ====
ll = sapply(profs, AIC)
plot(mults, ll)
Figure 3.9: Likelihood profile of scaled time invariant, Gislason natural mortalities. Y-axis = AIC, x-axis = scaling factor on nm.
# ====
We can now investigate if averaging the Gislason natural mortalities over the time series creates some scaling issues by ignoring weightings of numbers of observations contributing to each annual estimate. To partially address this we can investigate scaling these time invariant values up and down, to find where the best fit of this “shape” of natural mortalities across ages is.
By scaling the nm values up and down across a range of multipliers from 0.1 to 2, in increments of 0.1, we create 20 different runs from which we can compare AIC values. Here we find that a scaling value of 1.2 is responsible for the optimum model fit.
Therefore, below we investigate the model fit and residuals for a model with the Gislason, time invariant mortalities at age that have been scaled by a factor of 1.2 across all ages.
fit_nmGS <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG_S")
# === Create dataframe of current year's fit for plotting ==== IC_agg_temp <-
# aggregate(CATON~Year, ic_clean[ic_clean$CatchCategory %in% c('Landings',
# 'Discards'), ], FUN = 'sum')
asum_nmGS <- as.data.frame(summary(fit_nmGS))
asum_nmGS$Year <- as.integer(row.names(summary(fit_nmGS)))
ct_nmGS <- as.data.frame(catchtable(fit_nmGS, obs.show = TRUE))
ct_nmGS$Year <- as.integer(rownames(ct_nmGS))
rownames(ct_nmGS) <- NULL
ct_nmGS <- rbind(ct_nmGS, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmGS <- merge(x = asum_nmGS, y = ct_nmGS, by = "Year")
colnames(asum_nmGS) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
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: - The final model gradient: 141 - That there is a positive definite hessian: TRUE
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_nmGS <- split(fit_nmGS$data$logobs, ceiling(seq_along(fit_nmGS$data$logobs)/fit_nmGS$data$maxAgePerFleet[1]))
logPred_nmGS <- split(fit_nmGS$rep$predObs, ceiling(seq_along(fit_nmGS$rep$predObs)/fit_nmGS$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmGS_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmGS)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmGS[[i]])
logObs_nmGS_df <- rbind(logObs_nmGS_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmGS_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmGS)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmGS[[i]])
logPred_nmGS_df <- rbind(logPred_nmGS_df, temp_df)
}
## Plot
logPred_nmGS_df$age <- as.character(logPred_nmGS_df$age)
logObs_nmGS_df$age <- as.character(logObs_nmGS_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmGS_df[logObs_nmGS_df$fleet == 1 & logObs_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGS_df[logPred_nmGS_df$fleet == 1 & logPred_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.10: Model with SCALED time invariant Gislason natural mortalities fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmGS, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmGS_df[logObs_nmGS_df$fleet == 2 & logObs_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGS_df[logPred_nmGS_df$fleet == 2 & logPred_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.11: Model with SCALED time invariant Gislason natural mortalities fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmGS, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmGS_df[logObs_nmGS_df$fleet == 3 & logObs_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGS_df[logPred_nmGS_df$fleet == 3 & logPred_nmGS_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.12: Model with SCALED time invariant Gislason natural mortalities fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmGS, fleets=3)
sdplot(fit_nmGS, marg = c(5, 4, 1, 1))
Figure 3.13: Standard Devations by Fleet; Gislason, scaled, time-invariant mortality
resid_nmGS <- residuals(fit_nmGS)
resid_nmGS_df <- data.frame(year = resid_nmGS$year, fleet = resid_nmGS$fleet, age = resid_nmGS$age,
observation = resid_nmGS$observation, mean = resid_nmGS$mean, residual = resid_nmGS$residual)
resid_nmGS_df$fleetName <- ifelse(resid_nmGS_df$fleet == 1, attributes(resid_nmGS)$fleetNames[1],
ifelse(resid_nmGS_df$fleet == 2, attributes(resid_nmGS)$fleetNames[2], ifelse(resid_nmGS_df$fleet ==
3, attributes(resid_nmGS)$fleetNames[3], NA)))
resid_nmGS_df$fleetAltName <- ifelse(resid_nmGS_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmGS_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmGS_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmGS$conf$obsCorStruct == "ID")) {
corplot(fit_nmGS)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 3.14: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with SCALED time invariant Gislason natural mortalities.
ggplotly(ggplot(resid_nmGS_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 3.15: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with SCALED time invariant Gislason natural mortalities.
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_nmGS <- jit(fit = fit_nmGS)
mt <- as.data.frame(modeltable(jit_nmGS))
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.067 | 36 | 262.133 | M1 |
-95.067 | 36 | 262.133 | M2 |
-95.067 | 36 | 262.133 | M3 |
-95.067 | 36 | 262.133 | M4 |
-95.067 | 36 | 262.133 | M5 |
-95.067 | 36 | 262.133 | M6 |
-95.067 | 36 | 262.133 | M7 |
-95.067 | 36 | 262.133 | M8 |
-95.067 | 36 | 262.133 | M9 |
-95.067 | 36 | 262.133 | M10 |
-95.067 | 36 | 262.133 | 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_nmGS <- leaveout(fit_nmGS)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmGS$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmGS$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmGS$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmGS$`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_nmGS$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmGS$series <- rep("full", times = nrow(asum_nmGS))
asumi_nmGS <- asum_nmGS[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_nmGS <- rbind(woq1, woq34, asumi_nmGS)
# ssbplot(LO)
lossb_nmGS <- ggplotly(ggplot() + geom_line(data = losum_nmGS, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmGS, 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_nmGS <- ggplotly(ggplot() + geom_line(data = losum_nmGS, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmGS, 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_nmGS <- ggplotly(ggplot() + geom_line(data = losum_nmGS, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmGS, 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_nmGS <- ggplotly(ggplot() + geom_line(data = losum_nmGS, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmGS[asum_nmGS$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_nmGS, lof_nmGS, lorec_nmGS, loca_nmGS, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 3.16: Leave-one-out re-fits for the Time invariant, scaled, Gislason natural mortalities (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_nmGS <- as.data.frame(summary(fit_nmGS))
asum_nmGS$Year <- as.integer(row.names(summary(fit_nmGS)))
ct_nmGS <- as.data.frame(catchtable(fit_nmGS, obs.show = TRUE))
ct_nmGS$Year <- as.integer(rownames(ct_nmGS))
rownames(ct_nmGS) <- NULL
ct_nmGS <- rbind(ct_nmGS, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmGS <- merge(x = asum_nmGS, y = ct_nmGS, by = "Year")
colnames(asum_nmGS) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmGS$series <- rep("full", times = nrow(asum_nmGS))
# =====
RETRO_nmGS <- retro(fit_nmGS, year = 5)
rho_nmGS <- mohn(RETRO_nmGS, lag = 1)
## Make RETROs in better plotting format
ret_nmGS_df <- asum_nmGS[asum_nmGS$Year != max(asum_nmGS$Year), ]
for (i in 1:length(RETRO_nmGS)) {
tsum <- as.data.frame(summary(RETRO_nmGS[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmGS[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmGS[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmGS_df <- rbind(ret_nmGS_df, tsum)
}
ret_nmGS_df$series <- factor(ret_nmGS_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_nmGS <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year != max(asum_nmGS$Year),
], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh), fill = "black", alpha = 0.2) +
annotate(geom = "text", y = max(ret_nmGS_df$SSBhigh) * 0.85, x = ((max(ret_nmGS_df$Year) -
min(ret_nmGS_df$Year)) * 0.2) + min(ret_nmGS_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmGS[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmGS <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_nmGS_df$Fhigh) * 0.85, x = ((max(ret_nmGS_df$Year) -
min(ret_nmGS_df$Year)) * 0.8) + min(ret_nmGS_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmGS[3], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmGS <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh), fill = "black",
alpha = 0.3) + annotate(geom = "text", y = max(ret_nmGS_df$Rhigh) * 0.85, x = ((max(ret_nmGS_df$Year) -
min(ret_nmGS_df$Year)) * 0.2) + min(ret_nmGS_df$Year), label = paste0("Mohn's Rho = ",
round(rho_nmGS[1], digits = 3))) + scale_colour_manual(values = c("black", ebpal[(length(ebpal) -
4):length(ebpal)])) + theme_clean() + theme(legend.position = "none", axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmGS <- ggplotly(ggplot() + geom_line(data = ret_nmGS_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmGS[asum_nmGS$Year !=
max(asum_nmGS$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = ebpal[5],
fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmGS, retf_nmGS, retrec_nmGS, retca_nmGS, nrows = 2,
shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE, traces = c(8:((7 *
4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 3.17: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising SCALED time invariant, Gislason natural mortalities.
An alternate way to dealing with the issue of taking the unweighted means over the time series is to not use the full time-series means. While the life-history parameters \(L_{\infty}\) and \(K\) are calculated from all ages observations (q1 surveys) for the whole time series, they can be applied in the Gislason model to specific years, by utilising the lengths at age for each year. Due to the sampling rates these data will be inherently noisey, just as is the case for the stock weights at age and the maturity ogives, hence we need some method of smoothing out this noise. In this example we utilise a five-year sliding window mean on the annual estimates of natural mortality at age, to allow mortality at age to vary in time without introducing too much interannual variation.
fit_nmG5ysw <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG5ysw")
# === Create dataframe of current year's fit for plotting ====
asum_nmG5ysw <- as.data.frame(summary(fit_nmG5ysw))
asum_nmG5ysw$Year <- as.integer(row.names(summary(fit_nmG5ysw)))
ct_nmG5ysw <- as.data.frame(catchtable(fit_nmG5ysw, obs.show = TRUE))
ct_nmG5ysw$Year <- as.integer(rownames(ct_nmG5ysw))
rownames(ct_nmG5ysw) <- NULL
ct_nmG5ysw <- rbind(ct_nmG5ysw, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG5ysw <- merge(x = asum_nmG5ysw, y = ct_nmG5ysw, by = "Year")
colnames(asum_nmG5ysw) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
# =====
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: - The final model gradient: 137 - That there is a positive definite hessian: TRUE
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_nmG5ysw <- split(fit_nmG5ysw$data$logobs, ceiling(seq_along(fit_nmG5ysw$data$logobs)/fit_nmG5ysw$data$maxAgePerFleet[1]))
logPred_nmG5ysw <- split(fit_nmG5ysw$rep$predObs, ceiling(seq_along(fit_nmG5ysw$rep$predObs)/fit_nmG5ysw$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG5ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG5ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmG5ysw[[i]])
logObs_nmG5ysw_df <- rbind(logObs_nmG5ysw_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG5ysw_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG5ysw)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmG5ysw[[i]])
logPred_nmG5ysw_df <- rbind(logPred_nmG5ysw_df, temp_df)
}
## Plot
logPred_nmG5ysw_df$age <- as.character(logPred_nmG5ysw_df$age)
logObs_nmG5ysw_df$age <- as.character(logObs_nmG5ysw_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5ysw_df[logObs_nmG5ysw_df$fleet == 1 & logObs_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5ysw_df[logPred_nmG5ysw_df$fleet == 1 & logPred_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.18: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG5ysw, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5ysw_df[logObs_nmG5ysw_df$fleet == 2 & logObs_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5ysw_df[logPred_nmG5ysw_df$fleet == 2 & logPred_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.19: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5ysw, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG5ysw_df[logObs_nmG5ysw_df$fleet == 3 & logObs_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG5ysw_df[logPred_nmG5ysw_df$fleet == 3 & logPred_nmG5ysw_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logPred,
colour = age)) +
facet_wrap(facets = "age") + #, scales = "free_y") +
scale_color_manual(values = ebpal, guide = FALSE) +
guides(colour=FALSE, shape=FALSE) +
theme_few())
Figure 3.20: Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG5ysw, fleets=3)
sdplot(fit_nmG5ysw, marg = c(5, 4, 1, 1))
Figure 3.21: Standard Devations by Fleet; Gislason, time-variant mortality
resid_nmG5ysw <- residuals(fit_nmG5ysw)
resid_nmG5ysw_df <- data.frame(year = resid_nmG5ysw$year, fleet = resid_nmG5ysw$fleet,
age = resid_nmG5ysw$age, observation = resid_nmG5ysw$observation, mean = resid_nmG5ysw$mean,
residual = resid_nmG5ysw$residual)
resid_nmG5ysw_df$fleetName <- ifelse(resid_nmG5ysw_df$fleet == 1, attributes(resid_nmG5ysw)$fleetNames[1],
ifelse(resid_nmG5ysw_df$fleet == 2, attributes(resid_nmG5ysw)$fleetNames[2],
ifelse(resid_nmG5ysw_df$fleet == 3, attributes(resid_nmG5ysw)$fleetNames[3],
NA)))
resid_nmG5ysw_df$fleetAltName <- ifelse(resid_nmG5ysw_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG5ysw_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG5ysw_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG5ysw$conf$obsCorStruct == "ID")) {
corplot(fit_nmG5ysw)
# setcap('Estimated correlations', 'Estimates correlations between age
# groups for each fleet') stampit(fit)
} else {
print("No correlation structure configured for residuals across age.")
}
Figure 3.22: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities.
ggplotly(ggplot(resid_nmG5ysw_df) + geom_point(mapping = aes(x = year, y = age, size = abs(residual),
colour = residual >= 0), alpha = 0.7) + facet_grid(rows = "fleetAltName") + scale_colour_manual(values = c(`TRUE` = ebpal[8],
`FALSE` = ebpal[9])) + scale_y_continuous(limits = c(-1, 9), breaks = c(1:9)) +
guides(colour = FALSE) + theme_few())
Figure 3.23: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with time varying (five year sliiding window mean), age varying, Gislason natural mortalities.
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_nmG5ysw <- jit(fit = fit_nmG5ysw)
mt <- as.data.frame(modeltable(jit_nmG5ysw))
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.445 | 36 | 262.89 | M1 |
-95.445 | 36 | 262.89 | M2 |
-95.445 | 36 | 262.89 | M3 |
-95.445 | 36 | 262.89 | M4 |
-95.445 | 36 | 262.89 | M5 |
-95.445 | 36 | 262.89 | M6 |
-95.445 | 36 | 262.89 | M7 |
-95.445 | 36 | 262.89 | M8 |
-95.445 | 36 | 262.89 | M9 |
-95.445 | 36 | 262.89 | M10 |
-95.445 | 36 | 262.89 | 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_nmG5ysw <- leaveout(fit_nmG5ysw)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG5ysw$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG5ysw$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG5ysw$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG5ysw$`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_nmG5ysw$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG5ysw$series <- rep("full", times = nrow(asum_nmG5ysw))
asumi_nmG5ysw <- asum_nmG5ysw[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_nmG5ysw <- rbind(woq1, woq34, asumi_nmG5ysw)
# ssbplot(LO)
lossb_nmG5ysw <- ggplotly(ggplot() + geom_line(data = losum_nmG5ysw, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5ysw, 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_nmG5ysw <- ggplotly(ggplot() + geom_line(data = losum_nmG5ysw, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5ysw, 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_nmG5ysw <- ggplotly(ggplot() + geom_line(data = losum_nmG5ysw, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5ysw, 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_nmG5ysw <- ggplotly(ggplot() + geom_line(data = losum_nmG5ysw, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5ysw[asum_nmG5ysw$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_nmG5ysw, lof_nmG5ysw, lorec_nmG5ysw, loca_nmG5ysw, nrows = 2,
shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 3.24: Leave-one-out re-fits for the Time varying with five year sliding window, Gislason natural mortalities (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_nmG5ysw <- as.data.frame(summary(fit_nmG5ysw))
asum_nmG5ysw$Year <- as.integer(row.names(summary(fit_nmG5ysw)))
ct_nmG5ysw <- as.data.frame(catchtable(fit_nmG5ysw, obs.show = TRUE))
ct_nmG5ysw$Year <- as.integer(rownames(ct_nmG5ysw))
rownames(ct_nmG5ysw) <- NULL
ct_nmG5ysw <- rbind(ct_nmG5ysw, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG5ysw <- merge(x = asum_nmG5ysw, y = ct_nmG5ysw, by = "Year")
colnames(asum_nmG5ysw) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG5ysw$series <- rep("full", times = nrow(asum_nmG5ysw))
# =====
RETRO_nmG5ysw <- retro(fit_nmG5ysw, year = 5)
rho_nmG5ysw <- mohn(RETRO_nmG5ysw, lag = 1)
## Make RETROs in better plotting format
ret_nmG5ysw_df <- asum_nmG5ysw[asum_nmG5ysw$Year != max(asum_nmG5ysw$Year), ]
for (i in 1:length(RETRO_nmG5ysw)) {
tsum <- as.data.frame(summary(RETRO_nmG5ysw[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG5ysw[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG5ysw[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmG5ysw_df <- rbind(ret_nmG5ysw_df, tsum)
}
ret_nmG5ysw_df$series <- factor(ret_nmG5ysw_df$series, levels = c("full", "1", "2",
"3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG5ysw_df$SSBhigh) *
0.85, x = ((max(ret_nmG5ysw_df$Year) - min(ret_nmG5ysw_df$Year)) * 0.2) + min(ret_nmG5ysw_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5ysw[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmG5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5ysw_df$Fhigh) *
0.85, x = ((max(ret_nmG5ysw_df$Year) - min(ret_nmG5ysw_df$Year)) * 0.8) + min(ret_nmG5ysw_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5ysw[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmG5ysw <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG5ysw_df$Rhigh) *
0.85, x = ((max(ret_nmG5ysw_df$Year) - min(ret_nmG5ysw_df$Year)) * 0.2) + min(ret_nmG5ysw_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG5ysw[1], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmG5ysw <- ggplotly(ggplot() + geom_line(data = ret_nmG5ysw_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
max(asum_nmG5ysw$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3,
colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmG5ysw, retf_nmG5ysw, retrec_nmG5ysw, retca_nmG5ysw,
nrows = 2, shareX = FALSE, shareY = FALSE, titleY = TRUE), showlegend = FALSE,
traces = c(8:((7 * 4) + 2)))
layout(ret_fp, legend = list(orientation = "h", x = 0.5, xanchor = "center", y = -0.05,
yanchor = "top", borderwidth = 0))
Figure 3.25: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising time variant (5 year sliding window mean), Gislason natural mortalities.
fit_nmGFb <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_SV_nmGFb")
# === Create dataframe of current year's fit for plotting ====
asum_nmGFb <- as.data.frame(summary(fit_nmGFb))
asum_nmGFb$Year <- as.integer(row.names(summary(fit_nmGFb)))
ct_nmGFb <- as.data.frame(catchtable(fit_nmGFb, obs.show = TRUE))
ct_nmGFb$Year <- as.integer(rownames(ct_nmGFb))
rownames(ct_nmGFb) <- NULL
ct_nmGFb <- rbind(ct_nmGFb, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmGFb <- merge(x = asum_nmGFb, y = ct_nmGFb, by = "Year")
colnames(asum_nmGFb) <- 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: - The final model gradient: 145 - That there is a positive definite hessian: TRUE
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_nmGFb <- split(fit_nmGFb$data$logobs, ceiling(seq_along(fit_nmGFb$data$logobs)/fit_nmGFb$data$maxAgePerFleet[1]))
logPred_nmGFb <- split(fit_nmGFb$rep$predObs, ceiling(seq_along(fit_nmGFb$rep$predObs)/fit_nmGFb$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmGFb_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmGFb)) {
# 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_nmGFb[[i]])
logObs_nmGFb_df <- rbind(logObs_nmGFb_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmGFb_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmGFb)) {
# 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_nmGFb[[i]])
logPred_nmGFb_df <- rbind(logPred_nmGFb_df, temp_df)
}
## Plot
logPred_nmGFb_df$age <- as.character(logPred_nmGFb_df$age)
logObs_nmGFb_df$age <- as.character(logObs_nmGFb_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmGFb_df[logObs_nmGFb_df$fleet == 1 & logObs_nmGFb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGFb_df[logPred_nmGFb_df$fleet == 1 & logPred_nmGFb_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 time invariant Gislason natural mortalities (FishBase) fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmGFb, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmGFb_df[logObs_nmGFb_df$fleet == 2 & logObs_nmGFb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGFb_df[logPred_nmGFb_df$fleet == 2 & logPred_nmGFb_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 time invariant Gislason natural mortalities (FishBase) fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmGFb, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmGFb_df[logObs_nmGFb_df$fleet == 3 & logObs_nmGFb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmGFb_df[logPred_nmGFb_df$fleet == 3 & logPred_nmGFb_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 time invariant Gislason natural mortalities (FishBase) fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmGFb, fleets=3)
sdplot(fit_nmGFb, marg = c(5, 4, 1, 1))
Figure 4.4: Standard Devations by Fleet; Gislason, time-invariant mortality (FishBase)
resid_nmGFb <- residuals(fit_nmGFb)
resid_nmGFb_df <- data.frame(year = resid_nmGFb$year, fleet = resid_nmGFb$fleet,
age = resid_nmGFb$age, observation = resid_nmGFb$observation, mean = resid_nmGFb$mean,
residual = resid_nmGFb$residual)
resid_nmGFb_df$fleetName <- ifelse(resid_nmGFb_df$fleet == 1, attributes(resid_nmGFb)$fleetNames[1],
ifelse(resid_nmGFb_df$fleet == 2, attributes(resid_nmGFb)$fleetNames[2], ifelse(resid_nmGFb_df$fleet ==
3, attributes(resid_nmGFb)$fleetNames[3], NA)))
resid_nmGFb_df$fleetAltName <- ifelse(resid_nmGFb_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmGFb_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmGFb_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmGFb$conf$obsCorStruct == "ID")) {
corplot(fit_nmGFb)
# 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 time invariant Gislason natural mortalities (FishBase).
ggplotly(ggplot(resid_nmGFb_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 time invariant 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_nmGFb <- jit(fit = fit_nmGFb)
mt <- as.data.frame(modeltable(jit_nmGFb))
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.82 | 36 | 265.64 | M1 |
-96.82 | 36 | 265.64 | M2 |
-96.82 | 36 | 265.64 | M3 |
-96.82 | 36 | 265.64 | M4 |
-96.82 | 36 | 265.64 | M5 |
-96.82 | 36 | 265.64 | M6 |
-96.82 | 36 | 265.64 | M7 |
-96.82 | 36 | 265.64 | M8 |
-96.82 | 36 | 265.64 | M9 |
-96.82 | 36 | 265.64 | M10 |
-96.82 | 36 | 265.64 | 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_nmGFb <- leaveout(fit_nmGFb)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmGFb$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmGFb$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmGFb$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmGFb$`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_nmGFb$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmGFb$series <- rep("full", times = nrow(asum_nmGFb))
asumi_nmGFb <- asum_nmGFb[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_nmGFb <- rbind(woq1, woq34, asumi_nmGFb)
# ssbplot(LO)
lossb_nmGFb <- ggplotly(ggplot() + geom_line(data = losum_nmGFb, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmGFb, 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_nmGFb <- ggplotly(ggplot() + geom_line(data = losum_nmGFb, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmGFb, 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_nmGFb <- ggplotly(ggplot() + geom_line(data = losum_nmGFb, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmGFb, 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_nmGFb <- ggplotly(ggplot() + geom_line(data = losum_nmGFb, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmGFb[asum_nmGFb$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmGFb[asum_nmGFb$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_nmGFb, lof_nmGFb, lorec_nmGFb, loca_nmGFb, nrows = 2, shareX = TRUE,
titleY = TRUE), showlegend = FALSE)
Figure 4.7: Leave-one-out re-fits for the Time invariant, 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_nmGFb <- as.data.frame(summary(fit_nmGFb))
asum_nmGFb$Year <- as.integer(row.names(summary(fit_nmGFb)))
ct_nmGFb <- as.data.frame(catchtable(fit_nmGFb, obs.show = TRUE))
ct_nmGFb$Year <- as.integer(rownames(ct_nmGFb))
rownames(ct_nmGFb) <- NULL
ct_nmGFb <- rbind(ct_nmGFb, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmGFb <- merge(x = asum_nmGFb, y = ct_nmGFb, by = "Year")
colnames(asum_nmGFb) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmGFb$series <- rep("full", times = nrow(asum_nmGFb))
# =====
RETRO_nmGFb <- retro(fit_nmGFb, year = 5)
rho_nmGFb <- mohn(RETRO_nmGFb, lag = 1)
## Make RETROs in better plotting format
ret_nmGFb_df <- asum_nmGFb[asum_nmGFb$Year != max(asum_nmGFb$Year), ]
for (i in 1:length(RETRO_nmGFb)) {
tsum <- as.data.frame(summary(RETRO_nmGFb[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmGFb[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmGFb[[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_nmGFb_df <- rbind(ret_nmGFb_df, tsum)
}
ret_nmGFb_df$series <- factor(ret_nmGFb_df$series, levels = c("full", "1", "2", "3",
"4", "5"))
# ssbplot(RETRO)
retssb_nmGFb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGFb_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmGFb[asum_nmGFb$Year !=
max(asum_nmGFb$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmGFb_df$SSBhigh) *
0.85, x = ((max(ret_nmGFb_df$Year) - min(ret_nmGFb_df$Year)) * 0.2) + min(ret_nmGFb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmGFb[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_nmGFb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGFb_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmGFb[asum_nmGFb$Year !=
max(asum_nmGFb$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmGFb_df$Fhigh) *
0.85, x = ((max(ret_nmGFb_df$Year) - min(ret_nmGFb_df$Year)) * 0.8) + min(ret_nmGFb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmGFb[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_nmGFb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmGFb_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmGFb[asum_nmGFb$Year !=
max(asum_nmGFb$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmGFb_df$Rhigh) *
0.85, x = ((max(ret_nmGFb_df$Year) - min(ret_nmGFb_df$Year)) * 0.2) + min(ret_nmGFb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmGFb[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_nmGFb <- ggplotly(ggplot() + geom_line(data = ret_nmGFb_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmGFb[asum_nmGFb$Year !=
max(asum_nmGFb$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmGFb[asum_nmGFb$Year !=
max(asum_nmGFb$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_nmGFb, retf_nmGFb, retrec_nmGFb, retca_nmGFb, 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 utilising time invariant, Gislason natural mortalities (FishBase).
# === Scale natural mortality by a range of multipliers ====
fit_nmGFb <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_SV_nmGFb")
profs = list()
nm = fit_nmGFb$data$natMor
mults = seq(0.1, 2, by = 0.1)
## Re-fit to scaled mortalities
for (mult in mults) {
fit_nmGFb$data$natMor = nm * mult
profs[[as.character(mult)]] = stockassessment:::refit(fit_nmGFb, silent = TRUE)
}
# ====
# === Aggregate AIC values and plot results ====
ll = sapply(profs, AIC)
plot(mults, ll)
Figure 4.9: Likelihood profile of scaled time invariant, Gislason natural mortalities (FishBase). Y-axis = AIC, x-axis = scaling factor on nm.
# ====
We can now investigate if averaging the Gislason natural mortalities over the time series creates some scaling issues by ignoring weightings of numbers of observations contributing to each annual estimate. To partially address this we can investigate scaling these time invariant values up and down, to find where the best fit of this “shape” of natural mortalities across ages is.
By scaling the nm values up and down across a range of multipliers from 0.1 to 2, in increments of 0.1, we create 20 different runs from which we can compare AIC values. Here we find that a scaling value of 1.2 is responsible for the optimum model fit.
Therefore, below we investigate the model fit and residuals for a model with the Gislason, time invariant mortalities at age that have been scaled by a factor of 1.2 across all ages.
fit_nmG_SFb <- fitfromweb("ple.27.21-32_WKBPLAICE_2024_BioParsw_nmG_SFb")
# === Create dataframe of current year's fit for plotting ====
asum_nmG_SFb <- as.data.frame(summary(fit_nmG_SFb))
asum_nmG_SFb$Year <- as.integer(row.names(summary(fit_nmG_SFb)))
ct_nmG_SFb <- as.data.frame(catchtable(fit_nmG_SFb, obs.show = TRUE))
ct_nmG_SFb$Year <- as.integer(rownames(ct_nmG_SFb))
rownames(ct_nmG_SFb) <- NULL
ct_nmG_SFb <- rbind(ct_nmG_SFb, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG_SFb <- merge(x = asum_nmG_SFb, y = ct_nmG_SFb, by = "Year")
colnames(asum_nmG_SFb) <- 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: - The final model gradient: 132 - That there is a positive definite hessian: TRUE
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_nmG_SFb <- split(fit_nmG_SFb$data$logobs, ceiling(seq_along(fit_nmG_SFb$data$logobs)/fit_nmG_SFb$data$maxAgePerFleet[1]))
logPred_nmG_SFb <- split(fit_nmG_SFb$rep$predObs, ceiling(seq_along(fit_nmG_SFb$rep$predObs)/fit_nmG_SFb$data$maxAgePerFleet[1]))
## Transform data to useable dataframes
### Initialize an empty data frame for observations
logObs_nmG_SFb_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logObs = numeric())
### Populate the data frame for observations
for (i in seq_along(logObs_nmG_SFb)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for observations
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logObs = logObs_nmG_SFb[[i]])
logObs_nmG_SFb_df <- rbind(logObs_nmG_SFb_df, temp_df)
}
### Initialize an empty data frame for predictions
logPred_nmG_SFb_df <- data.frame(age = integer(),
fleet = integer(),
year = integer(),
logPred = numeric())
### Populate the data frame for predictions
for (i in seq_along(logPred_nmG_SFb)) {
# Calculate the group and year
year <- ((i - 1) %/% 3) + 1
fleet <- ((i - 1) %% 3) + 1
### Create a temporary data frame and append to the main data frame for predictions
temp_df <- data.frame(age = 1:7,
fleet = fleet,
year = year+2001,
logPred = logPred_nmG_SFb[[i]])
logPred_nmG_SFb_df <- rbind(logPred_nmG_SFb_df, temp_df)
}
## Plot
logPred_nmG_SFb_df$age <- as.character(logPred_nmG_SFb_df$age)
logObs_nmG_SFb_df$age <- as.character(logObs_nmG_SFb_df$age)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_SFb_df[logObs_nmG_SFb_df$fleet == 1 & logObs_nmG_SFb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_SFb_df[logPred_nmG_SFb_df$fleet == 1 & logPred_nmG_SFb_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 SCALED time invariant Gislason natural mortalities (FishBase) fit to catch data by age (values on the model link function scale).
# fitplot(fit_nmG_SFb, fleets=1)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_SFb_df[logObs_nmG_SFb_df$fleet == 2 & logObs_nmG_SFb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_SFb_df[logPred_nmG_SFb_df$fleet == 2 & logPred_nmG_SFb_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 SCALED time invariant Gislason natural mortalities (FishBase) fit to Q1 survey data by age (values on the model link function scale).
# fitplot(fit_nmG_SFb, fleets=2)
ggplotly(ggplot()+
geom_point(data = logObs_nmG_SFb_df[logObs_nmG_SFb_df$fleet == 3 & logObs_nmG_SFb_df$year != (DataYear+1), ],
mapping = aes(x = year,
y = logObs,
colour = age),
shape = 21) +
geom_line(data = logPred_nmG_SFb_df[logPred_nmG_SFb_df$fleet == 3 & logPred_nmG_SFb_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.12: Model with SCALED time invariant Gislason natural mortalities (FishBase) fit to Q3/4 survey data by age (values on the model link function scale).
# fitplot(fit_nmG_SFb, fleets=3)
sdplot(fit_nmG_SFb, marg = c(5, 4, 1, 1))
Figure 4.13: Standard Devations by Fleet; Gislason, scaled, time-invariant mortality (FishBase)
resid_nmG_SFb <- residuals(fit_nmG_SFb)
resid_nmG_SFb_df <- data.frame(year = resid_nmG_SFb$year, fleet = resid_nmG_SFb$fleet,
age = resid_nmG_SFb$age, observation = resid_nmG_SFb$observation, mean = resid_nmG_SFb$mean,
residual = resid_nmG_SFb$residual)
resid_nmG_SFb_df$fleetName <- ifelse(resid_nmG_SFb_df$fleet == 1, attributes(resid_nmG_SFb)$fleetNames[1],
ifelse(resid_nmG_SFb_df$fleet == 2, attributes(resid_nmG_SFb)$fleetNames[2],
ifelse(resid_nmG_SFb_df$fleet == 3, attributes(resid_nmG_SFb)$fleetNames[3],
NA)))
resid_nmG_SFb_df$fleetAltName <- ifelse(resid_nmG_SFb_df$fleet == 1, "Residual Fishing",
ifelse(resid_nmG_SFb_df$fleet == 2, "Q1 Surveys", ifelse(resid_nmG_SFb_df$fleet ==
3, "Q3/4 Surveys", NA)))
if (!all(fit_nmG_SFb$conf$obsCorStruct == "ID")) {
corplot(fit_nmG_SFb)
# 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.14: Estimated correlations in residual variation between ages for each of the fishing fleet (top) and the two surveys (middle & bottom), from the Model with SCALED time invariant Gislason natural mortalities (FishBase).
ggplotly(ggplot(resid_nmG_SFb_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.15: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the Model with SCALED time invariant Gislason natural mortalities (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_nmG_SFb <- jit(fit = fit_nmG_SFb)
mt <- as.data.frame(modeltable(jit_nmG_SFb))
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.066 | 36 | 262.132 | M1 |
-95.066 | 36 | 262.132 | M2 |
-95.066 | 36 | 262.132 | M3 |
-95.066 | 36 | 262.132 | M4 |
-95.066 | 36 | 262.132 | M5 |
-95.066 | 36 | 262.132 | M6 |
-95.066 | 36 | 262.132 | M7 |
-95.066 | 36 | 262.132 | M8 |
-95.066 | 36 | 262.132 | M9 |
-95.066 | 36 | 262.132 | M10 |
-95.066 | 36 | 262.132 | 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_nmG_SFb <- leaveout(fit_nmG_SFb)
# === Get data from sam objects and generate useable dataframes ====
q1mat <- as.data.frame(summary(LO_nmG_SFb$`w.o. Q1IBTS+BITS+CodSD21-25`))
q3mat <- as.data.frame(summary(LO_nmG_SFb$`w.o. Q34IBTS+BITS+CodSD21-25`))
woq1 <- data.frame(Year = as.integer(row.names(summary(LO_nmG_SFb$`w.o. Q1IBTS+BITS+CodSD21-25`))),
SSB = q1mat$SSB, Fbar = q1mat$`Fbar(3-5)`, R_age1 = q1mat$`R(age 1)`, CatchEst = catchtable(LO_nmG_SFb$`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_nmG_SFb$`w.o. Q34IBTS+BITS+CodSD21-25`)[,
1], series = rep("wo_Q34", times = (nrow(q3mat) - 1)))
asum_nmG_SFb$series <- rep("full", times = nrow(asum_nmG_SFb))
asumi_nmG_SFb <- asum_nmG_SFb[, c("Year", "SSB", "Fbar", "R_age1", "CatchEst", "series")]
# =====
losum_nmG_SFb <- rbind(woq1, woq34, asumi_nmG_SFb)
# ssbplot(LO)
lossb_nmG_SFb <- ggplotly(ggplot() + geom_line(data = losum_nmG_SFb, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG_SFb, 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_nmG_SFb <- ggplotly(ggplot() + geom_line(data = losum_nmG_SFb, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG_SFb, 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_nmG_SFb <- ggplotly(ggplot() + geom_line(data = losum_nmG_SFb, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG_SFb, 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_nmG_SFb <- ggplotly(ggplot() + geom_line(data = losum_nmG_SFb, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
(DataYear + 1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG_SFb[asum_nmG_SFb$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_nmG_SFb, lof_nmG_SFb, lorec_nmG_SFb, loca_nmG_SFb, nrows = 2,
shareX = TRUE, titleY = TRUE), showlegend = FALSE)
Figure 4.16: Leave-one-out re-fits for the Time invariant, scaled, 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_nmG_SFb <- as.data.frame(summary(fit_nmG_SFb))
asum_nmG_SFb$Year <- as.integer(row.names(summary(fit_nmG_SFb)))
ct_nmG_SFb <- as.data.frame(catchtable(fit_nmG_SFb, obs.show = TRUE))
ct_nmG_SFb$Year <- as.integer(rownames(ct_nmG_SFb))
rownames(ct_nmG_SFb) <- NULL
ct_nmG_SFb <- rbind(ct_nmG_SFb, data.frame(Estimate = NA, Low = NA, High = NA, sop.catch = NA,
Year = as.integer(2024)))
asum_nmG_SFb <- merge(x = asum_nmG_SFb, y = ct_nmG_SFb, by = "Year")
colnames(asum_nmG_SFb) <- c("Year", "R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh",
"Fbar", "Flow", "Fhigh", "CatchEst", "Catchlow", "Catchhigh", "CatchObs")
asum_nmG_SFb$series <- rep("full", times = nrow(asum_nmG_SFb))
# =====
RETRO_nmG_SFb <- retro(fit_nmG_SFb, year = 5)
rho_nmG_SFb <- mohn(RETRO_nmG_SFb, lag = 1)
## Make RETROs in better plotting format
ret_nmG_SFb_df <- asum_nmG_SFb[asum_nmG_SFb$Year != max(asum_nmG_SFb$Year), ]
for (i in 1:length(RETRO_nmG_SFb)) {
tsum <- as.data.frame(summary(RETRO_nmG_SFb[[i]]))
tsum$Year <- as.integer(row.names(summary(RETRO_nmG_SFb[[i]])))
tsum <- cbind(tsum, catchtable(RETRO_nmG_SFb[[i]], obs.show = TRUE))
tsum$series <- as.character(rep(i, nrow(tsum)))
colnames(tsum) <- c("R_age1", "Rlow", "Rhigh", "SSB", "SSBlow", "SSBhigh", "Fbar",
"Flow", "Fhigh", "Year", "CatchEst", "Catchlow", "Catchhigh", "CatchObs",
"series")
tsum <- tsum[tsum$Year != max(tsum$Year), ]
ret_nmG_SFb_df <- rbind(ret_nmG_SFb_df, tsum)
}
ret_nmG_SFb_df$series <- factor(ret_nmG_SFb_df$series, levels = c("full", "1", "2",
"3", "4", "5"))
# ssbplot(RETRO)
retssb_nmG_SFb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_SFb_df, mapping = aes(x = Year,
y = SSB, colour = series)) + geom_ribbon(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
max(asum_nmG_SFb$Year), ], mapping = aes(x = Year, ymin = SSBlow, ymax = SSBhigh),
fill = "black", alpha = 0.2) + annotate(geom = "text", y = max(ret_nmG_SFb_df$SSBhigh) *
0.85, x = ((max(ret_nmG_SFb_df$Year) - min(ret_nmG_SFb_df$Year)) * 0.2) + min(ret_nmG_SFb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG_SFb[2], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Fplot(RETRO)
retf_nmG_SFb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_SFb_df, mapping = aes(x = Year,
y = Fbar, colour = series)) + geom_ribbon(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
max(asum_nmG_SFb$Year), ], mapping = aes(x = Year, ymin = Flow, ymax = Fhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG_SFb_df$Fhigh) *
0.85, x = ((max(ret_nmG_SFb_df$Year) - min(ret_nmG_SFb_df$Year)) * 0.8) + min(ret_nmG_SFb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG_SFb[3], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# recplot(RETRO)
retrec_nmG_SFb <- layout(ggplotly(ggplot() + geom_line(data = ret_nmG_SFb_df, mapping = aes(x = Year,
y = R_age1, colour = series)) + geom_ribbon(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
max(asum_nmG_SFb$Year), ], mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh),
fill = "black", alpha = 0.3) + annotate(geom = "text", y = max(ret_nmG_SFb_df$Rhigh) *
0.85, x = ((max(ret_nmG_SFb_df$Year) - min(ret_nmG_SFb_df$Year)) * 0.2) + min(ret_nmG_SFb_df$Year),
label = paste0("Mohn's Rho = ", round(rho_nmG_SFb[1], digits = 3))) + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(legend.position = "none",
axis.title.y.left = element_text(vjust = -0.05, hjust = 0.95))), showlegend = FALSE)
# Catch plot (RETRO)
retca_nmG_SFb <- ggplotly(ggplot() + geom_line(data = ret_nmG_SFb_df, mapping = aes(x = Year,
y = CatchEst, colour = series)) + geom_ribbon(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
max(asum_nmG_SFb$Year), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh),
fill = "black", alpha = 0.3) + geom_point(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
max(asum_nmG_SFb$Year), ], mapping = aes(x = Year, y = CatchObs), shape = 3,
colour = ebpal[5], fill = ebpal[5]) + ylab("Catch (tonnes)") + scale_colour_manual(values = c("black",
ebpal[(length(ebpal) - 4):length(ebpal)])) + theme_clean() + theme(axis.title.y.left = element_text(vjust = -0.05,
hjust = 0.95)))
ret_fp <- style(subplot(retssb_nmG_SFb, retf_nmG_SFb, retrec_nmG_SFb, retca_nmG_SFb,
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.17: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising SCALED time invariant, Gislason natural mortalities (FishBase).
An alternate way to dealing with the issue of taking the unweighted means over the time series is to not use the full time-series means. While the life-history parameters \(L_{\infty}\) and \(K\) are calculated from all ages observations (q1 surveys) for the whole time series, they can be applied in the Gislason model to specific years, by utilising the lengths at age for each year. Due to the sampling rates these data will be inherently noisey, just as is the case for the stock weights at age and the maturity ogives, hence we need some method of smoothing out this noise. In this example we utilise a five-year sliding window mean on the annual estimates of natural mortality at age, to allow mortality at age to vary in time without introducing too much interannual variation.
fit_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: - The final model gradient: 136 - That there is a positive definite hessian: TRUE
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 4.18: 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 4.19: 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 4.20: 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 4.21: 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 4.22: 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 4.23: One observation ahead residuals for the three fleets (red/pink = observation lower than model estimate, blue = observation higher than model estimate, size = magnitude of residual), from the 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 4.24: 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 4.25: Retrospective analyses for SSB (top-left), F for ages 3-5 (top-right), recruitment (bottom-left), and catch (bottom-right), from the model utilising time variant (5 year sliding window mean), Gislason natural mortalities (FishBase).
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("Historic M", "Time invariant, Gislason", "Optimally scaled, time invariant, Gislason",
"Time varying (5ysw), Gislason", "Time invariant, Gislason & FishBase", "Optimally scaled, time invariant, Gislason & FishBase",
"Time varying (5ysw), Gislason & FishBase"), AIC = c(AIC(fit_bio_5y_SV), AIC(fit_nmG),
AIC(fit_nmGS), AIC(fit_nmG5ysw), AIC(fit_nmGFb), AIC(fit_nmG_SFb), AIC(fit_nmG5Fb)),
No.Params = c(length(coef(fit_bio_5y_SV)), length(coef(fit_nmG)), length(coef(fit_nmGS)),
length(coef(fit_nmG5ysw)), length(coef(fit_nmGFb)), length(coef(fit_nmG_SFb)),
length(coef(fit_nmG5Fb))), Rho.SSB = c(rho_bio_5y_SV[2], rho_nmG[2], rho_nmGS[2],
rho_nmG5ysw[2], rho_nmGFb[2], rho_nmG_SFb[2], rho_nmG5Fb[2]), Rho.F = c(rho_bio_5y_SV[3],
rho_nmG[3], rho_nmGS[3], rho_nmG5ysw[3], rho_nmGFb[3], rho_nmG_SFb[3], rho_nmG5Fb[3]),
Rho.Rec = c(rho_bio_5y_SV[1], rho_nmG[1], rho_nmGS[1], rho_nmG5ysw[1], rho_nmGFb[1],
rho_nmG_SFb[1], rho_nmG5Fb[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 |
---|---|---|---|---|---|
Historic M | 327.2155 | 36 | 0.1201919 | -0.0259916 | 0.3158176 |
Time invariant, Gislason | 266.3056 | 36 | 0.0961543 | -0.0177057 | 0.2170079 |
Optimally scaled, time invariant, Gislason | 262.1331 | 36 | 0.0778651 | 0.0072884 | 0.1811353 |
Time varying (5ysw), Gislason | 262.8897 | 36 | 0.0625855 | -0.0177211 | 0.1195901 |
Time invariant, Gislason & FishBase | 265.6402 | 36 | 0.0945449 | -0.0156215 | 0.2137940 |
Optimally scaled, time invariant, Gislason & FishBase | 262.1320 | 36 | 0.0759293 | 0.0101150 | 0.1773563 |
Time varying (5ysw), Gislason & FishBase | 262.3204 | 36 | 0.0600907 | -0.0158922 | 0.1149631 |
Now that we’ve considered the model fits, diagnostics and tuning, we might as well have a look at how the models’ estimates compare.
# ssblines <- data.frame(yintercept=c(4730, 4370, 3635), Lines = factor(x = c("MSY-Btrigger", "Bpa", "Blim"),levels = c("MSY-Btrigger", "Bpa", "Blim")))
pssb <- ggplot() +
geom_line(data = asum_bio_5y_SV,
mapping = aes(x=Year,
y=SSB,
series="Historic nm"),
colour = ebpal[1]) +
geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh,
series="Historic nm"),
fill = ebpal[1],
alpha = 0.3) +
geom_line(data = asum_nmG,
mapping = aes(x=Year,
y=SSB,
series="Time-invariant Gislason nm"),
colour = ebpal[2]) +
geom_ribbon(data = asum_nmG,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh,
series="Time-invariant Gislason nm"),
fill = ebpal[2],
alpha = 0.3) +
geom_line(data = asum_nmGS,
mapping = aes(x=Year,
y=SSB,
series="Scaled, time-invariant Gislason nm"),
colour = ebpal[3]) +
geom_ribbon(data = asum_nmGS,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh,
series="Scaled, time-invariant Gislason nm"),
fill = ebpal[3],
alpha = 0.3) +
geom_line(data = asum_nmG5ysw,
mapping = aes(x=Year,
y=SSB,
series="Time varying (5ysw) Gislason nm"),
colour = ebpal[4]) +
geom_ribbon(data = asum_nmG5ysw,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh,
series="Time varying (5ysw) Gislason nm"),
fill = ebpal[4],
alpha = 0.3) +
geom_line(data = asum_nmGFb,
mapping = aes(x=Year,
y=SSB,
series="Time invariant, Gislason & FishBase"),
colour = ebpal[5]) +
geom_ribbon(data = asum_nmGFb,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh,
series="Time invariant, Gislason & FishBase"),
fill = ebpal[5],
alpha = 0.3) +
geom_line(data = asum_nmG_SFb,
mapping = aes(x=Year,
y=SSB,
series="Time invariant, scaled Gislason & FishBase"),
colour = ebpal[6]) +
geom_ribbon(data = asum_nmG_SFb,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh,
series="Time invariant, scaled Gislason & FishBase"),
fill = ebpal[6],
alpha = 0.3) +
geom_line(data = asum_nmG5Fb,
mapping = aes(x=Year,
y=SSB,
series="Time varying (5ysw) Gislason & FishBase"),
colour = ebpal[7]) +
geom_ribbon(data = asum_nmG5Fb,
mapping = aes(x=Year,
ymin=SSBlow,
ymax=SSBhigh),
series="Time varying (5ysw) Gislason & FishBase",
fill = ebpal[7],
alpha = 0.3) +
theme_clean()
ggplotly(pssb, tooltip = c("series", "x", "y"))
Figure 4.26: Estimated spawning stock biomass for ple.27.21-32 and 95% confidence intervals (tonnes). Purple is the model with the historic mortalities, dark green is with time invariant Gislason mortalities, dark orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities. Yellow, light-green, and light-orange, are in the same order as the previous three, except all mortalities are calculated with life-history data from FishBase, instead of survey observations.
Fishing pressure over time seems to respond to the changes in the estimations of SSB, while input data for catches (obviously) remain the same betwen the two models.
# flines <- data.frame(yintercept=c(0.31, 0.81, 1.00), Lines = factor(x =
# c('FMSY', 'Fpa', 'Flim'), levels = c('FMSY', 'Fpa', 'Flim')))
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = Fbar,
series = "Historic nm"), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh, series = "Historic nm"), fill = ebpal[1],
alpha = 0.3) + geom_line(data = asum_nmG, mapping = aes(x = Year, y = Fbar, series = "Time-invariant Gislason nm"),
colour = ebpal[2]) + geom_ribbon(data = asum_nmG, mapping = aes(x = Year, ymin = Flow,
ymax = Fhigh, series = "Time-invariant Gislason nm"), fill = ebpal[2], alpha = 0.3) +
geom_line(data = asum_nmGS, mapping = aes(x = Year, y = Fbar, series = "Scaled, time-invariant Gislason nm"),
colour = ebpal[3]) + geom_ribbon(data = asum_nmGS, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh, series = "Scaled, time-invariant Gislason nm"), fill = ebpal[3],
alpha = 0.3) + geom_line(data = asum_nmG5ysw, mapping = aes(x = Year, y = Fbar,
series = "Time varying (5ysw) Gislason nm"), colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw,
mapping = aes(x = Year, ymin = Flow, ymax = Fhigh, series = "Time varying (5ysw) Gislason nm"),
fill = ebpal[4], alpha = 0.3) + geom_line(data = asum_nmGFb, mapping = aes(x = Year,
y = Fbar, series = "Time invariant, Gislason & FishBase"), colour = ebpal[5]) +
geom_ribbon(data = asum_nmGFb, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh,
series = "Time invariant, Gislason & FishBase"), fill = ebpal[5], alpha = 0.3) +
geom_line(data = asum_nmG_SFb, mapping = aes(x = Year, y = Fbar, series = "Time invariant, scaled Gislason & FishBase"),
colour = ebpal[6]) + geom_ribbon(data = asum_nmG_SFb, mapping = aes(x = Year,
ymin = Flow, ymax = Fhigh, series = "Time invariant, scaled Gislason & FishBase"),
fill = ebpal[6], alpha = 0.3) + geom_line(data = asum_nmG5Fb, mapping = aes(x = Year,
y = Fbar, series = "Time varying (5ysw) Gislason & FishBase"), colour = ebpal[7]) +
geom_ribbon(data = asum_nmG5Fb, mapping = aes(x = Year, ymin = Flow, ymax = Fhigh,
series = "Time varying (5ysw) Gislason & FishBase"), fill = ebpal[7], alpha = 0.3) +
theme_clean())
Figure 4.27: Annual fishing mortality estimates for ple.27.21-32 ages 3-5 and point wise 95% confidence intervals are shown by line and shaded area. Purple is the model with the historic mortalities, dark green is with time invariant Gislason mortalities, dark orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities. Yellow, light-green, and light-orange, are in the same order as the previous three, except all mortalities are calculated with life-history data from FishBase, instead of survey observations.
Estimations of recruitment in the most recent years fall, as SSB also falls in this version of the model, again probably due to the reduction in stock weights at age.
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV, mapping = aes(x = Year, y = R_age1,
series = "Historic nm"), colour = ebpal[1]) + geom_ribbon(data = asum_bio_5y_SV,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "Historic nm"), fill = ebpal[1],
alpha = 0.3) + geom_line(data = asum_nmG, mapping = aes(x = Year, y = R_age1,
series = "Time-invariant Gislason nm"), colour = ebpal[2]) + geom_ribbon(data = asum_nmG,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "Time-invariant Gislason nm"),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_nmGS, mapping = aes(x = Year,
y = R_age1, series = "Scaled, time-invariant Gislason nm"), colour = ebpal[3]) +
geom_ribbon(data = asum_nmGS, mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh,
series = "Scaled, time-invariant Gislason nm"), fill = ebpal[3], alpha = 0.3) +
geom_line(data = asum_nmG5ysw, mapping = aes(x = Year, y = R_age1, series = "Time varying (5ysw) Gislason nm"),
colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw, mapping = aes(x = Year,
ymin = Rlow, ymax = Rhigh, series = "Time varying (5ysw) Gislason nm"), fill = ebpal[4],
alpha = 0.3) + geom_line(data = asum_nmGFb, mapping = aes(x = Year, y = R_age1,
series = "Time invariant, Gislason & FishBase"), colour = ebpal[5]) + geom_ribbon(data = asum_nmGFb,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "Time invariant, Gislason & FishBase"),
fill = ebpal[5], alpha = 0.3) + geom_line(data = asum_nmG_SFb, mapping = aes(x = Year,
y = R_age1, series = "Time invariant, scaled Gislason & FishBase"), colour = ebpal[6]) +
geom_ribbon(data = asum_nmG_SFb, mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh,
series = "Time invariant, scaled Gislason & FishBase"), fill = ebpal[6],
alpha = 0.3) + geom_line(data = asum_nmG5Fb, mapping = aes(x = Year, y = R_age1,
series = "Time varying (5ysw) Gislason & FishBase"), colour = ebpal[7]) + geom_ribbon(data = asum_nmG5Fb,
mapping = aes(x = Year, ymin = Rlow, ymax = Rhigh, series = "Time varying (5ysw) Gislason & FishBase"),
fill = ebpal[7], alpha = 0.3) + theme_clean())
Figure 4.28: Five year sliding window model annual recruitment estimates for ple.27.21-32 and point wise 95% confidence intervals are shown by line and shaded area (numbers). Purple is the model with the historic mortalities, dark green is with time invariant Gislason mortalities, dark orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities. Yellow, light-green, and light-orange, are in the same order as the previous three, except all mortalities are calculated with life-history data from FishBase, instead of survey observations.
ggplotly(ggplot() + geom_line(data = asum_bio_5y_SV[asum_bio_5y_SV$Year != (DataYear +
1), ], mapping = aes(x = Year, y = CatchEst, series = "Historic nm"), colour = ebpal[1]) +
geom_ribbon(data = asum_bio_5y_SV[asum_bio_5y_SV$Year != (DataYear + 1), ], mapping = aes(x = Year,
ymin = Catchlow, ymax = Catchhigh, series = "Historic nm"), fill = ebpal[1],
alpha = 0.3) + geom_line(data = asum_nmG[asum_nmG$Year != (DataYear + 1),
], mapping = aes(x = Year, y = CatchEst, series = "Time-invariant Gislason nm"),
colour = ebpal[2]) + geom_ribbon(data = asum_nmG[asum_nmG$Year != (DataYear +
1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh, series = "Time-invariant Gislason nm"),
fill = ebpal[2], alpha = 0.3) + geom_line(data = asum_nmGS[asum_nmGS$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "Scaled, time-invariant Gislason nm"),
colour = ebpal[3]) + geom_ribbon(data = asum_nmGS[asum_nmGS$Year != (DataYear +
1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh, series = "Scaled, time-invariant Gislason nm"),
fill = ebpal[3], alpha = 0.3) + geom_line(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "Time varying (5ysw) Gislason nm"),
colour = ebpal[4]) + geom_ribbon(data = asum_nmG5ysw[asum_nmG5ysw$Year != (DataYear +
1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh, series = "Time varying (5ysw) Gislason nm"),
fill = ebpal[4], alpha = 0.3) + geom_point(data = asum_nmG5ysw[asum_nmG5ysw$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchObs), shape = 3, colour = "black",
fill = "black") + geom_line(data = asum_nmGFb[asum_nmGFb$Year != (DataYear +
1), ], mapping = aes(x = Year, y = CatchEst, series = "Time invariant, Gislason & FishBase"),
colour = ebpal[5]) + geom_ribbon(data = asum_nmGFb[asum_nmGFb$Year != (DataYear +
1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh, series = "Time invariant, Gislason & FishBase"),
fill = ebpal[5], alpha = 0.3) + geom_line(data = asum_nmG_SFb[asum_nmG_SFb$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "Time invariant, scaled Gislason & FishBase"),
colour = ebpal[6]) + geom_ribbon(data = asum_nmG_SFb[asum_nmG_SFb$Year != (DataYear +
1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh, series = "Time invariant, scaled Gislason & FishBase"),
fill = ebpal[6], alpha = 0.3) + geom_line(data = asum_nmG5Fb[asum_nmG5Fb$Year !=
(DataYear + 1), ], mapping = aes(x = Year, y = CatchEst, series = "Time varying (5ysw) Gislason & FishBase"),
colour = ebpal[7]) + geom_ribbon(data = asum_nmG5Fb[asum_nmG5Fb$Year != (DataYear +
1), ], mapping = aes(x = Year, ymin = Catchlow, ymax = Catchhigh, series = "Time varying (5ysw) Gislason & FishBase"),
fill = ebpal[7], alpha = 0.3) + ylab("Catch (tonnes)") + theme_clean())
Figure 4.29: Annual catch estimates and 95% confidence intervals (line and shaded area, respectively) for ple.27.21-32 and point annual observations (points). Purple is the model with the historic mortalities, dark green is with time invariant Gislason mortalities, dark orange is with scaled-time-invariant Gislason mortalities, and blue is with time varying Gislason mortalities. Yellow, light-green, and light-orange, are in the same order as the previous three, except all mortalities are calculated with life-history data from FishBase, instead of survey observations.
From the diagnostics of these plots we can see that the inclusion of the new natural mortality estimates greatly improves the model fits. Furthermore, we see that having age and time varying natural mortality provides the best fit overall. This is likely due to the fact that the Gislason mortalities are based on both life-history traits and size, and therefore, as the size at age has been decreasing over time, the time varying nm better matches this change than does a fixed, average nm based on average size. The option of time-varying nm is also internally coherent, as we have included time-varying stock weights-at-age as the best representation of stock development and so we should reflect this development in our estimations of natural mortality as well.
Hence, as discussed in the workshop plenary, there are three lines of reasoning for utilising the time & age varying Gislason et al (2010) natural mortality estimates:
In the next working document we will investigate the utility of setting maturity at age-1 to zero, based on basic biological knowledge of the species, and the fact that there may be some age reading errors meaning that older fish are being aged to 1. Furthermore, we will investigate a range of ways of including discard survival in the assessment process before finally attempting to include recently available German recreational fisheries data.