library(FLCore)
library(FLXSA)
#> Error in library(FLXSA): there is no package called 'FLXSA'
library(FLasher)
#> Error in library(FLasher): there is no package called 'FLasher'
library(ggplotFL)
library(patchwork)
library(tidyverse)
# remotes::install_github("einarhjorleifsson/fishvice", dependencies = FALSE)
library(fishvice)little helper functions
# a retro functions
lh_retro <- function(end.year, fls, ind, cnt) {
window(fls, end = end.year) +
FLXSA(window(fls, end = end.year),
window(ind, end = end.year),
cnt)
}
# get the reference biomass next year
lh_1year_projection <- function(fls) {
# the mean recruitment model
sr <- as.FLSR(fls, model="geomean")
params(sr)['a',] <- exp(mean(log(rec(fls))))
# use same weights and Fay (gets muliplied by Flevel) as in the terminal year
proj <- stf(fls, nyears = 1, wts.nyears = 1)
last.year <- range(fls)["maxyear"][[1]]
fsq <- mean(fbar(fls)[,as.character(last.year)])
ctrl <-
fwdControl(data.frame(year = last.year + 1,
quant = "f",
value = fsq))
fwd(proj, control = ctrl, sr = sr)
}
# a compilation function (here only biomass, B4+)
lh_retro_B4p <- function(ryears, fls, ind, cnt) {
retro <- map(ryears, lh_retro, fls, ind, cnt)
# add one year projection to get the biomass in the ass year
# note the smb index in the ass year is NOT used in the fitting, need to check
if(projections) retro <- map(retro, lh_1year_projection)
names(retro) <- ryears
bio <-
map(retro, computeStock) %>%
map(as.data.frame) %>%
map(as_tibble) %>%
bind_rows(.id = "assyear") %>%
mutate(assyear = as.integer(assyear)) %>%
ungroup()
bio.last <-
bio %>%
filter(assyear == max(assyear)) %>%
select(year, bio = data)
bio <-
bio %>%
left_join(bio.last) %>%
mutate(r = log(data / bio))
# just so no need to alter code below
return(list(bio = bio))
}
# the plot
lh_ggplot <- function(res, last.point, bio.spaly) {
p1 <-
res %>%
ggplot() +
geom_line(data = bio.spaly,
aes(year, bio),
colour = "yellow",
lwd = 1) +
geom_line(aes(year, data, group = assyear)) +
geom_point(data = last.point, colour = "red",
aes(year, data)) +
facet_wrap(~ run) +
labs(x = NULL, y = "B4+")
p2 <-
res %>%
ggplot(aes(year, r, group = assyear)) +
geom_line() +
geom_point(data = last.point,
colour = "red") +
facet_wrap(~ run) +
labs(x = NULL, y = "contemporaneous\n vs current")
return(p1 + p2 + plot_layout(ncol = 1))
}
# doing retro statistic using the last ~5 years of the assessment
# is plain stupid, even though ices does so
lh_retro_performance <- function(d) {
left_join(d %>%
group_by(run) %>% summarise(bias98.21 = mean(r), cv98.21 = sd(r)),
d %>% filter(year %in% 1998:2017) %>%
group_by(run) %>% summarise(bias98.17 = mean(r), cv98.17 = sd(r)))
}get data
YRS <- 1985:2021 # it is only fair to use 1994 onwards, the time the
# estimates start, but here just trying to repeat the paper
fil <- "/u2/reikn/Tac/2022/01-cod/data-raw/Supplementary S3.xlsx"
shits <- readxl::excel_sheets(fil)
mc <-
expand_grid(year = YRS,
age = 1:14) %>%
left_join(bind_rows(readxl::read_excel(fil, sheet = shits[1]) %>% mutate(sur = "smb"),
readxl::read_excel(fil, sheet = shits[2]) %>% mutate(sur = "smh"))) %>%
spread(sur, Mc) %>%
mutate(mc = case_when(year < 1994 ~ 0.2,
is.na(smh) ~ smb,
TRUE ~ (smb + smh) / 2)) %>%
select(year, age, mc) %>%
arrange(year, age) %>%
group_by(year) %>%
# assume same mc in all year after last datayear
fill(mc, .direction = "downup") %>%
# 10+ same as 9
group_by(age) %>%
fill(mc) %>%
ungroup() %>%
filter(year %in% YRS)
rbx <- fishvice::mup_rbx("/u2/reikn/Tac/2022/01-cod/ass/mup/smx", scale = 1000)
bio.spaly <- rbx$rby %>% filter(year %in% min(YRS):2022) %>% select(year, bio)
# have not figured out how to set xsa up with survey data in year with no catch data
# besides of course the old-way, but can not be bothered with that now
rbx$rby <- rbx$rby %>% filter(year %in% YRS)
rbx$rbya <- rbx$rbya %>% filter(year %in% YRS)
# set maturity and stock weights such that SSB is actually B4+
rbx$rbya <-
rbx$rbya %>%
mutate(mat = replace_na(mat, 0),
mat = ifelse(age >= 4, 1, 0),
ssbW = sW, # 2022-05-28 from cW to sW as in paper
ssbW = replace_na(ssbW, 0))
rbx.mc <- rbx
# replace 0.2 with mc
rbx.mc$rbya <-
rbx.mc$rbya %>%
select(-m) %>%
left_join(mc %>% rename(m = mc))flr setup
# pf and pm set to 0, so "SSB" = B4+ is based on stock in numbers in start of
# year
# issue: age 14 should not be a plus group, for relative comparisons
# should though be ok
fls <- rbx %>% rbx_to_flstock(pf = 0, pm = 0, fage = c(5, 10))
# range(fls)["plusgroup"] <- NA
fls.mc <- rbx.mc %>% rbx_to_flstock(pf = 0, pm = 0, fage = c(5, 10))
# range(fls.mc)["plusgroup"] <- NA
ind.smb <-
rbx_to_flindex(rbx, "U1", time = c(3/12, 3.5/12)) %>% trim(year = YRS)
# range(ind.smb)["plusgroup"] <- NA
# ind.smh <-
# rbx_to_flindex(rbx, "U2", time = c(10/12, 10.5/12)) %>%
# trim(year = 1996:2021)
# only spring survey used in paper
ind <- FLIndices(ind.smb)xsa retrospectives
controls
# power or no power, that is the question
cntrp0 <- FLXSA.control(maxit = 70, rage = 0, qage = 14) # i would set qage to
#> Error in FLXSA.control(maxit = 70, rage = 0, qage = 14): could not find function "FLXSA.control"
cntrp5 <- FLXSA.control(maxit = 70, rage = 5, qage = 14) # age 9 or 10
#> Error in FLXSA.control(maxit = 70, rage = 5, qage = 14): could not find function "FLXSA.control"run retro
projections <- FALSE
ryears <- 1998:max(YRS)
retp0 <- lh_retro_B4p(ryears, fls, ind, cntrp0)
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
retp0.mc <- lh_retro_B4p(ryears, fls.mc, ind, cntrp0)
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
retp5 <- lh_retro_B4p(ryears, fls, ind, cntrp5)
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
retp5.mc <- lh_retro_B4p(ryears, fls.mc, ind, cntrp5)
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"results
clarification on runs:
~run, ~comment,
p002 No power, M = 0.20
p0mc No power, M = mc
p502 Power ages 1:5, M = 0.20
p5mc Power ages 1:5, M = mc
res <- bind_rows(retp0$bio %>% mutate(run = "p002"),
retp0.mc$bio %>% mutate(run = "p0mc"),
retp5$bio %>% mutate(run = "p502"),
retp5.mc$bio %>% mutate(run = "p5mc"))
#> Error in mutate(., run = "p002"): object 'retp0' not found
last.point <-
res %>%
filter(year == assyear + projections)
#> Error in filter(., year == assyear + projections): object 'res' not found
lh_ggplot(res, last.point, bio.spaly)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'ggplot': object 'res' not foundconclusions on the retrospectives
- The variable M (p0mc) show better retrospective performance than constant M (p002) if no power is assumed.
- Both variable and constant M show improved retrospectives if a power up to and including age 5 is assumed.
- When power is assumed, the retrospective performance of constant M (p502) outperformes the variable M (p5mc).
There may though be more under the hood in the variable M diagnostics that may warrant further investigation. Something beyond that done in the Mc article.
sanity check - the mc input
m(fls.mc) %>% trim(year = 1993:2021, age = 1:10) %>% round(2)
#> An object of class "FLQuant"
#> , , unit = unique, season = all, area = unique
#>
#> year
#> age 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
#> 1 0.20 0.06 0.08 0.06 0.08 0.02 0.08 0.11 0.05 0.13 0.19 0.13 0.31 0.14 0.19
#> 2 0.20 0.06 0.08 0.06 0.08 0.02 0.08 0.11 0.05 0.13 0.19 0.13 0.31 0.14 0.19
#> 3 0.20 0.08 0.09 0.12 0.15 0.09 0.11 0.15 0.12 0.15 0.27 0.22 0.25 0.22 0.20
#> 4 0.20 0.13 0.13 0.10 0.35 0.17 0.16 0.17 0.16 0.17 0.19 0.25 0.12 0.22 0.15
#> 5 0.20 0.18 0.14 0.09 0.24 0.26 0.17 0.17 0.07 0.17 0.25 0.23 0.28 0.26 0.18
#> 6 0.20 0.13 0.07 0.09 0.25 0.23 0.22 0.21 0.12 0.16 0.21 0.16 0.28 0.27 0.16
#> 7 0.20 0.14 0.02 0.11 0.17 0.20 0.14 0.15 0.16 0.16 0.24 0.19 0.23 0.31 0.20
#> 8 0.20 0.05 0.08 0.15 0.22 0.17 0.16 0.18 0.19 0.05 0.10 0.22 0.25 0.24 0.28
#> 9 0.20 0.19 0.02 0.20 0.27 0.29 0.29 0.22 0.10 0.40 0.18 0.37 0.18 0.37 0.37
#> 10 0.20 0.19 0.02 0.20 0.27 0.29 0.29 0.22 0.10 0.40 0.18 0.37 0.18 0.37 0.37
#> year
#> age 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
#> 1 0.17 0.27 0.38 0.18 0.11 0.29 0.15 0.16 0.42 0.32 0.22 0.26 0.17 0.19
#> 2 0.17 0.27 0.38 0.18 0.11 0.29 0.15 0.16 0.42 0.32 0.22 0.26 0.17 0.19
#> 3 0.21 0.26 0.42 0.24 0.20 0.24 0.25 0.26 0.27 0.29 0.42 0.15 0.20 0.19
#> 4 0.15 0.27 0.18 0.12 0.17 0.21 0.15 0.38 0.29 0.10 0.22 0.19 0.10 0.07
#> 5 0.18 0.27 0.27 0.09 0.19 0.15 0.25 0.23 0.28 0.15 0.17 0.20 0.18 0.07
#> 6 0.25 0.21 0.20 0.11 0.19 0.27 0.23 0.31 0.17 0.25 0.16 0.19 0.12 0.13
#> 7 0.21 0.19 0.18 0.15 0.19 0.23 0.27 0.39 0.27 0.16 0.16 0.20 0.11 0.18
#> 8 0.33 0.25 0.19 0.07 0.15 0.20 0.21 0.21 0.24 0.31 0.16 0.24 0.26 0.18
#> 9 0.23 0.47 0.25 0.16 0.09 0.17 0.18 0.20 0.22 0.24 0.15 0.29 0.19 0.13
#> 10 0.23 0.47 0.25 0.16 0.09 0.17 0.18 0.20 0.22 0.24 0.15 0.29 0.19 0.13
#>
#> units: NAresults of retro xsa with 1 year projections
… still not using the survey in the assessment year
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
#> Error in `map()`:
#> ℹ In index: 1.
#> Caused by error in `FLXSA()`:
#> ! could not find function "FLXSA"
#> Error in mutate(., run = "p002"): object 'retp0' not found
#> Error in filter(., year == assyear + projections): object 'res' not found
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'ggplot': object 'res' not found
#> Error in group_by(., run): object 'last.point' not foundchange log
2022-05-28:
* changed basis of reference biomass from cW to sW, reflecting what i think was
done in the paper. B4+ is actually based on catch weights. That basis
is of course wrong, but it is linked to the 0.20 multiplier in the HCR that
also includes a stabilizer.
Whatever the case, for some reasons XSA gives a much higher biomass when
using catch weights than muppet, or for that matter sam. Most likely related
to mishandling of the plusgroup. But these details are not an issue in this
doodle.
* added one year prediction, i.e. into the advise year, to reflect what was
actually done in the paper. note though that not using the survey index in
the advice year is still an issue (can not be done in XSA, i think).
* also simplified code, but that alone does not affecting previous versions
of this document
* added different periods for estimating retrospective performance. not yet
though repeating the 2 five year peels in the paper because that is not
kosher.
2022-05-27
* converted lh_function, now used purrr::map rather than lapply
* added window to the index as well, just for clarity
* neither affects the results nor the conclusion
* added conclusions and some minor additional text edits
2022-05-26
* first version - bb & co notified via email