if (dir.exists(file.path("..", "R")) && requireNamespace("pkgload", quietly = TRUE)) {
pkgload::load_all("..", quiet = TRUE)
} else {
library(decamod)
}
library(TMBhelper)
library(tidyverse)
library(kableExtra)
theme_set(theme_bw())Load inputs
Below the required R libraries are loaded. The decamod RTMB model is loaded along with several R functions using library(decamod) (Webber and McKenzie 2022, 2024). The kableExtra package is used for generating tables of inputs. The code theme_set(theme_bw()) alters the plot aesthetics. The model is implemented using automatic differentiation via RTMB (Kristensen et al. 2016).
Dimensions
Could have different midpoints by sex. Could have different breaks (i.e., not even). But would probably need to have the same number of length classes? In the future LFs could be estimated to 1mm and then aggregated however the user wants.
# # non equal bin widths - this is what we should aim for
# binwidth <- 2 # mm
# breaks <- c(30, 35, 40, seq(45, 92, binwidth))
# # non equal bin widths - this is to accommodate LFs for aggregation as
# # LFs are in 2mm bins
# binwidth <- 2 # mm
# breaks <- c(30, 34, 38, 42, seq(46, 92, binwidth))
# midpoint_l <- breaks[-length(breaks)] + diff(breaks) / 2
# breaks[length(breaks)] <- Inf
# bins <- cut(x = midpoint_l, breaks = breaks, right = FALSE)
breaks <- c(seq(30, 90, 2), Inf)
bins <- cut(x = seq(31, 91, 2), breaks = breaks, right = FALSE)
midpoint_l <- midcut(seq(31, 91, 2), 30, 92, 2)
bins <- cut(x = midpoint_l, breaks = breaks, right = FALSE)
tibble(bins, midpoint_l) %>% head()
#> # A tibble: 6 × 2
#> bins midpoint_l
#> <fct> <dbl>
#> 1 [30,32) 31
#> 2 [32,34) 33
#> 3 [34,36) 35
#> 4 [36,38) 37
#> 5 [38,40) 39
#> 6 [40,42) 41
tibble(bins, midpoint_l) %>% tail()
#> # A tibble: 6 × 2
#> bins midpoint_l
#> <fct> <dbl>
#> 1 [80,82) 81
#> 2 [82,84) 83
#> 3 [84,86) 85
#> 4 [86,88) 87
#> 5 [88,90) 89
#> 6 [90,Inf) 91
tibble(test = c(29, 30, 30.1, 31.9, 32, 35, 90, 91, 91.9, 120)) %>%
mutate(bin = cut(x = test, breaks = breaks, right = FALSE))
#> # A tibble: 10 × 2
#> test bin
#> <dbl> <fct>
#> 1 29 <NA>
#> 2 30 [30,32)
#> 3 30.1 [30,32)
#> 4 31.9 [30,32)
#> 5 32 [32,34)
#> 6 35 [34,36)
#> 7 90 [90,Inf)
#> 8 91 [90,Inf)
#> 9 91.9 [90,Inf)
#> 10 120 [90,Inf)
first_year <- 1945
last_year <- 2023
l_year <- first_year:last_year
l_season <- c("AW", "SS")
l_region <- "CRA 2"
l_sex <- c("Male", "Female")
l_fishery <- c("SL", "NSL")
l_length <- bins
l_cpue <- c("FSU", "CELR", "Logbook")
n_year <- length(l_year)
n_season <- length(l_season)
n_region <- length(l_region)
n_fishery <- length(l_fishery)
n_sex <- length(l_sex)
n_length <- length(midpoint_l)
dim_ytrsl <- c(n_year, n_season, n_region, n_sex, n_length)
dimnames_ytrsl <- list(year = l_year, season = l_season, region = l_region,
sex = l_sex, length = l_length)
dim_ytrfl <- c(n_year, n_season, n_region, n_fishery, n_length)
dimnames_ytrfl <- list(year = l_year, season = l_season, region = l_region,
fishery = l_fishery, length = l_length)
dim_ytrsfl <- c(n_year, n_season, n_region, n_sex, n_fishery, n_length)
dimnames_ytrsfl <- list(year = l_year, season = l_season, region = l_region,
sex = l_sex, fishery = l_fishery, length = l_length)Tags
Tag-recapture data are used to estimate growth increments following Schnute and Fournier (1980) and Fournier and Archibald (1982).
tags <- tags_CRA2 %>%
rename(gi_L1 = RelLngth, gi_L2 = RecLngth) %>%
mutate(gi_obs = gi_L2 - gi_L1, gi_liberty = DaysatLiberty / 365) %>%
mutate(gi_sex = as.integer(Sex)) %>%
mutate(gi_morph = gi_sex) %>%
filter(TagRelease == 0) %>%
# filter(RelYear %in% first_year:last_year) %>%
# filter(RecYear %in% first_year:last_year) %>%
filter(gi_liberty > 0.5, gi_obs > -10)
# select(gi_sex, gi_morph, gi_L1, gi_liberty, gi_obs)
which_growth_ytrs <- array(0, dim = c(n_year, n_season, n_region, n_sex))
which_growth_ytrs[,,,1] <- 1
which_growth_ytrs[,,,2] <- 2
col.names <- c("Sex", "Initial length (mm)", "Time at liberty (years)",
"Growth increment (mm)")
tags %>% select(gi_sex, gi_L1, gi_liberty, gi_obs) %>% head(6) %>%
kbl(digits = 1, col.names = col.names,
caption = "Example tag-recapture records used to estimate growth increments.") %>%
add_header_above(c("Covariates" = 3, "Data" = 1))Catch
catch <- catch_CRA2 %>%
mutate(SL = Commercial + Recreational, NSL = Illegal + Customary) %>%
filter(Year >= first_year)
catch_ytrf <- array(0, dim = c(n_year, n_season, n_region, n_fishery),
dimnames = list(year = l_year, season = l_season,
region = l_region, fishery = l_fishery))
for (y in 1:n_year) {
for (t in 1:n_season) {
if (any(catch$Season == t & catch$Year %in% l_year[y])) {
catch_ytrf[y,t,1,1] <- catch$SL[catch$Season == t & catch$Year %in% l_year[y]]
catch_ytrf[y,t,1,2] <- catch$NSL[catch$Season == t & catch$Year %in% l_year[y]]
}
}
}
kbl(tail(catch %>% select(-QMA)), digits = 1,
caption = "Recent catch inputs by fishing year, season, and fishery component.")| Year | Season | Commercial | Recreational | Illegal | Customary | SL | NSL |
|---|---|---|---|---|---|---|---|
| 2021 | 1 | 47.5 | 3.2 | 13.5 | 0.5 | 50.7 | 14.0 |
| 2021 | 2 | 32.5 | 28.5 | 9.3 | 4.5 | 61.0 | 13.8 |
| 2022 | 1 | 46.9 | 1.3 | 12.8 | 0.5 | 48.2 | 13.3 |
| 2022 | 2 | 33.1 | 11.6 | 9.0 | 4.5 | 44.7 | 13.5 |
| 2023 | 1 | 36.3 | 1.2 | 11.7 | 0.5 | 37.5 | 12.2 |
| 2023 | 2 | 43.7 | 10.4 | 14.0 | 4.5 | 54.1 | 18.5 |
Legal status and retention
Includes MLS, proportion of females in berry, legal status, and proportion retained:
mls_ytrsfl <- array(1, dim = dim_ytrsfl, dimnames = dimnames_ytrsfl)
for (t in 1:n_season) {
for (r in 1:n_region) {
for (y in which(l_year < 1988)) {
mls_ytrsfl[y, t, r, 1, 1,] <- ifelse(midpoint_l > 53, 1, 0)
mls_ytrsfl[y, t, r, 2, 1,] <- ifelse(midpoint_l > 58, 1, 0)
}
for (y in which(l_year %in% 1988:1991)) {
mls_ytrsfl[y, t, r, 1, 1,] <- ifelse(midpoint_l > 54, 1, 0)
mls_ytrsfl[y, t, r, 2, 1,] <- ifelse(midpoint_l > 58, 1, 0)
}
for (y in which(l_year > 1991)) {
mls_ytrsfl[y, t, r, 1, 1,] <- ifelse(midpoint_l > 54, 1, 0)
mls_ytrsfl[y, t, r, 2, 1,] <- ifelse(midpoint_l > 60, 1, 0)
}
}
}
berry_ytrfl <- array(1, dim = dim_ytrfl, dimnames = dimnames_ytrfl)
for (y in 1:n_year) {
for (t in 1:n_season) {
sub <- berry_CRA2 %>% filter(Season == t)
if (nrow(sub) > 0) {
berry_ytrfl[y,t,,1,] <- 1 - sub$mean_prop_berry
}
}
}
legal_ytrsfl <- mls_ytrsfl
for (y in 1:n_year) {
for (t in 1:n_season) {
for (r in 1:n_region) {
for (f in 1:n_fishery) {
legal_ytrsfl[y,t,r,2,f,] <- legal_ytrsfl[y,t,r,2,f,] * berry_ytrfl[y,t,r,f,]
}
}
}
}
retained_ytrsfl <- array(1, dim = dim_ytrsfl, dimnames = dimnames_ytrsfl)LFs
Length frequency data are weighted using the iterative re-weighting approach of Francis (2011). For now I have summed over IF and MF and rescaled to 1. I’ve also added code to sum over length bins if needed.
lf <- lf_CRA2 %>%
arrange(fyear, season, sex) %>%
mutate(sex = ifelse(sex == 1, 1, 2)) %>%
mutate(season = ifelse(season == "AW", 1, 2)) %>%
pivot_longer(cols = starts_with("tw_"), names_to = "length") %>%
mutate(length = parse_number(length)) %>%
mutate(bin = cut(x = length, breaks = breaks, right = FALSE)) %>%
mutate(value = value * weight_scaled) %>%
group_by(fyear, season, sex, bin) %>%
summarise(value = sum(value)) %>%
pivot_wider(names_from = bin, values_from = value) %>%
ungroup() %>%
mutate(weight_scaled = rowSums(select(., -c(fyear, season, sex)))) %>%
relocate(weight_scaled, .after = sex)
lf_obs <- as.matrix(lf %>% select(-c(fyear, season, sex, weight_scaled)))
lf_obs <- lf_obs / rowSums(lf_obs)Sex ratio
Just using proportion male.
CPUE
Catch per unit effort (CPUE) indices are standardised following Maunder and Punt (2004).
fsu <- fsu_CRA2 %>% mutate(q = 1, units = 0)
celr <- celr_CRA2 %>% mutate(q = 2, units = 0)
logbook <- logbook_CRA2 %>% mutate(q = 3, units = 1)
cpue <- bind_rows(fsu, celr, logbook) %>%
separate(Year, into = c("year", "season")) %>%
mutate(year = as.integer(year)) %>%
mutate(season = as.integer(ifelse(season == "AW", 1, 2))) %>%
mutate(logSD = log(1 + SD / Mean)) %>%
mutate(wt = case_when(q == 1 ~ 0.75, q == 2 ~ 1.35, q == 3 ~ 1.31)) %>%
mutate(pro = 0) # mutate(pro = 0.25)
cpue %>% select(year, season, Median, Mean, SD, logSD, q, units) %>%
tail(10) %>%
kbl(digits = 1,
caption = "Recent CPUE index observations and associated uncertainty inputs.")| year | season | Median | Mean | SD | logSD | q | units |
|---|---|---|---|---|---|---|---|
| 2019 | 1 | 1.1 | 1.1 | 0.2 | 0.2 | 3 | 1 |
| 2019 | 2 | 1.4 | 1.4 | 0.3 | 0.2 | 3 | 1 |
| 2020 | 1 | 1.3 | 1.3 | 0.2 | 0.2 | 3 | 1 |
| 2020 | 2 | 1.6 | 1.6 | 0.3 | 0.2 | 3 | 1 |
| 2021 | 1 | 1.4 | 1.4 | 0.2 | 0.2 | 3 | 1 |
| 2021 | 2 | 1.2 | 1.2 | 0.3 | 0.2 | 3 | 1 |
| 2022 | 1 | 1.2 | 1.2 | 0.2 | 0.2 | 3 | 1 |
| 2022 | 2 | 1.3 | 1.3 | 0.3 | 0.2 | 3 | 1 |
| 2023 | 1 | 1.3 | 1.3 | 0.2 | 0.2 | 3 | 1 |
| 2023 | 2 | 1.6 | 1.6 | 0.5 | 0.2 | 3 | 1 |
# plot(cpue$SD, cpue$logSD); abline(coef = c(0, 1), col = 2)
# plot(cpue$CV, cpue$logSD); abline(coef = c(0, 1), col = 2)
x <- seq(0.32, 0.43, length.out = 10)
wt <- 1.2
y <- sqrt(x^2 + 0.25^2) * 1.0 / wt
plot(x)
lines(y, col = 2) # effect of adding process error is to make sds more similar
Maturity
Maturity is set to zero for males.
maturity_ytrsl <- array(0, dim = dim_ytrsl, dimnames = dimnames_ytrsl)
for (y in 1:n_year) {
for (t in 1:n_season) {
for (r in 1:n_region) {
maturity_ytrsl[y,t,r,2,] <- logistic_ogive(x = midpoint_l,
x50 = 50.4482, x95 = 11.2724)
}
}
}
as.data.frame(ftable(maturity_ytrsl)) %>%
filter(year == last_year, season == "AW", sex == "Female") %>%
ggplot(aes(x = length, y = Freq)) +
geom_point() +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "Tail width bin (mm)", y = "Proportion mature") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
Weight
weight_a <- c(4.16e-06, 1.30e-05)
weight_b <- c(2.9354, 2.5452)
n_weight <- length(weight_a)
which_weight_ytrs <- array(0, dim = c(n_year, n_season, n_region, n_sex))
which_weight_ytrs[,,,1] <- 1
which_weight_ytrs[,,,2] <- 2
weight_il <- matrix(0, nrow = n_weight, ncol = n_length)
for (i in 1:n_weight) {
weight_il[i,] <- length_to_weight(midpoint_l, weight_a[i], weight_b[i])
}
weight_ytrsl <- get_weight(weight_il, which_weight_ytrs)
dimnames(weight_ytrsl) <- dimnames_ytrsl
plot(midpoint_l, weight_ytrsl[1,1,1,1,], type = "l", col = 4)
lines(midpoint_l, weight_ytrsl[1,1,1,2,], type = "l", col = 2)
Is M at midpoints same as M integrated over lengths try v small bins vs large bins.
midpoint2 <- seq(31, 91, 2)
M2 <- get_M_length(midpoint2, a = 60, Ma = 0.2, Mc = -1)
df2 <- tibble(midpoint2 = midpoint2, M = M2, N = 10000) %>%
mutate(NeM = N * exp(-M))
midpoint1 <- seq(31, 92, 1)
M1 <- get_M_length(midpoint1, a = 60, Ma = 0.2, Mc = -1)
df1 <- tibble(midpoint1 = midpoint1, M = M1, N = 5000) %>%
mutate(NeM = N * exp(-M)) %>%
mutate(id = cut(x = midpoint1, breaks = seq(30, 92, 2), labels = FALSE)) %>%
mutate(midpoint2 = midpoint2[id]) %>%
group_by(midpoint2) %>%
summarise(N = sum(N), NeM = sum(NeM))
sum(df1$N)
#> [1] 310000
sum(df2$N)
#> [1] 310000
sum(df1$NeM)
#> [1] 250700.3
sum(df2$NeM)
#> [1] 250191.4
plot(midpoint2, df1$NeM, type = "l")
lines(midpoint2, df2$NeM, col = 2)
Selectivity
Selectivity is modelled as a double-normal ogive (Quinn and Deriso 1999).
Define data list
Define the data list that can be passed to MakeADFun:
data <- list(
first_year = first_year,
n_year = n_year, n_season = n_season, n_region = n_region, n_sex = n_sex,
n_fishery = n_fishery,
midpoint_l = as.numeric(midpoint_l),
l_year = l_year, l_season = l_season, l_region = l_region, l_sex = l_sex,
l_fishery = l_fishery, l_length = l_length, l_cpue = l_cpue,
weight_ytrsl = weight_ytrsl, maturity_ytrsl = maturity_ytrsl,
# Selectivty
n_sel = 4, # M AW, F AW, M SS, F SS
sel_like = 0, # 0=off, 1=RE
# sel_type = c(1, 2, 1, 2), # 0=constant 1, 1=logistic, 2=logistic-capped
sel_type = c(2, 2, 2, 2), # 0=constant 1, 1=logistic, 2=logistic-capped
which_sel_ytrsf = which_sel_ytrsf,
# M
M_size = 60,
which_M_ytrs = array(1, dim = c(n_year, n_season, n_region, n_sex)),
# Recruitment
rec_like = 5, # 0=off, 1=RE, 2=RW, 3=AR1, 4=RE/RW, 5=RE/AR1
rsigma_bias = 0,
last_rdev_year = n_year,
# Growth
growth_like = 1, # 0=off, 1=normal, 2=robust-normal, 3=student-t, 4=Cauchy
growth_dd = 1, # density dependent growth: 0=off, 1=on
n_growth = 2,
n_tags = nrow(tags),
ga = 50, gb = 80,
gmin = 0.0001, # gmin = -10,
which_growth_ytrs = which_growth_ytrs,
# need which statements for each growth parameter
# Catch
catch_like = -1, # -1=Us, 0=Fs using NR, 1=normal, 2=lognormal
catch_obs = as.vector(catch_ytrf),
catch_cv = 0.06,
pen_wt = 1000,
match_old_u = 0,
catch_ytrf = catch_ytrf,
mls_ytrsfl = mls_ytrsfl,
legal_ytrsfl = legal_ytrsfl,
retained_ytrsfl = retained_ytrsfl,
handling_mortality_y = ifelse(l_year < 1990, 0.1, 0.05),
# CPUE
cpue_like = 2, # 0=off, 1=normal, 2=lognormal, 3=gamma
cpue_first_qcreep = c(1990, 1990, 1993) - first_year + 1,
cpue_obs = cpue$Median,
cpue_sd = cpue$logSD,
cpue_pro = cpue$pro,
cpue_wt = c(0.75, 0.7, 1)[cpue$q],
cpue_q = cpue$q,
cpue_units = cpue$units, # 0=numbers, 1=weight
cpue_fishery = rep(1, nrow(cpue)),
cpue_year = cpue$year - first_year + 1,
cpue_season = cpue$season, # AW, SS
cpue_region = rep(1, nrow(cpue)),
cpue_sex = rep(0, nrow(cpue)), # 0=all, 1=female, 2=male
# Sex-ratio
sexr_like = 2, # 0=off, 1=old, 2=logit-normal, 3=beta
sexr_obs = sexr$mean_scaled_1,
sexr_sd = sexr$weight_scaled,
sexr_wt = rep(800, nrow(sexr)),
sexr_year = sexr$fyear - first_year + 1,
sexr_season = sexr$season,
sexr_region = rep(1, nrow(sexr)),
sexr_fishery = rep(1, nrow(sexr)),
# LFs
lf_like = 1, # 1=multinomial, 2=Dirichlet
lf_obs = lf_obs,
lf_min = rep(8, nrow(lf)),
lf_max = rep(n_length - 4, nrow(lf)),
# lf_max = rep(n_length - 1, nrow(lf)),
lf_roll = rep(1, nrow(lf)),
lf_sd = lf$weight_scaled,
lf_wt = ifelse(lf$sex == 1, 1.204, 0.982),
lf_year = lf$fyear - first_year + 1,
lf_season = lf$season,
lf_region = rep(1, nrow(lf)),
lf_sex = lf$sex,
lf_fishery = rep(1, nrow(lf)),
# Switches
jacobian_on = 0, prior_on = 1, like_on = 1)
data <- c(data, as.list(tags))
names(data)
#> [1] "first_year" "n_year" "n_season"
#> [4] "n_region" "n_sex" "n_fishery"
#> [7] "midpoint_l" "l_year" "l_season"
#> [10] "l_region" "l_sex" "l_fishery"
#> [13] "l_length" "l_cpue" "weight_ytrsl"
#> [16] "maturity_ytrsl" "n_sel" "sel_like"
#> [19] "sel_type" "which_sel_ytrsf" "M_size"
#> [22] "which_M_ytrs" "rec_like" "rsigma_bias"
#> [25] "last_rdev_year" "growth_like" "growth_dd"
#> [28] "n_growth" "n_tags" "ga"
#> [31] "gb" "gmin" "which_growth_ytrs"
#> [34] "catch_like" "catch_obs" "catch_cv"
#> [37] "pen_wt" "match_old_u" "catch_ytrf"
#> [40] "mls_ytrsfl" "legal_ytrsfl" "retained_ytrsfl"
#> [43] "handling_mortality_y" "cpue_like" "cpue_first_qcreep"
#> [46] "cpue_obs" "cpue_sd" "cpue_pro"
#> [49] "cpue_wt" "cpue_q" "cpue_units"
#> [52] "cpue_fishery" "cpue_year" "cpue_season"
#> [55] "cpue_region" "cpue_sex" "sexr_like"
#> [58] "sexr_obs" "sexr_sd" "sexr_wt"
#> [61] "sexr_year" "sexr_season" "sexr_region"
#> [64] "sexr_fishery" "lf_like" "lf_obs"
#> [67] "lf_min" "lf_max" "lf_roll"
#> [70] "lf_sd" "lf_wt" "lf_year"
#> [73] "lf_season" "lf_region" "lf_sex"
#> [76] "lf_fishery" "jacobian_on" "prior_on"
#> [79] "like_on" "Event" "Sex"
#> [82] "RelMonth" "RelYear" "RecMonth"
#> [85] "RecYear" "DaysatLiberty" "gi_L1"
#> [88] "gi_L2" "TWConvertRel" "TWConvertRec"
#> [91] "TagRelease" "RelStatArea" "Condition"
#> [94] "Tagtype" "RecStatArea" "gi_obs"
#> [97] "gi_liberty" "gi_sex" "gi_morph"
# sqrt(data$cpue_sd^2 + data$cpue_pro^2) / data$cpue_wt # cpue_sigma
# hist(rlnorm(n = 10000, meanlog = log(2), sdlog = sqrt(data$cpue_sd^2 + data$cpue_pro^2) / data$cpue_wt))
# ii <- 11
# yy <- data$lf_obs[ii,4:31] / sum(data$lf_obs[ii,4:31])
# xx <- gtools::rdirichlet(n = 10000, alpha = yy * data$lf_sd[ii] * data$lf_wt[ii])
# plot(yy, type = "b", col = 2, ylim = c(0, 0.25))
# for (i in 1:100) lines(xx[i,] / sum(xx[i,]), col = 4)
# data$lf_obs <- round(data$lf_obs[,data$lf_min[1]:data$lf_max[1]] * data$lf_sd * data$lf_wt)
# data$lf_N <- rowSums(data$lf_obs)Model setup
Define the parameter list:
parameters <- list(
log_F0 = log(0), logit_U0 = qlogis(1e-06),
log_R0 = log(554439), log_Rsigma = log(0.25), Rrho = 0,
logit_h = qlogis(0.999), log_Rmu = log(32), log_Rsd = log(2),
log_Ma = log(c(0.187377, 0.187377)), Mc = c(0, 0), #Mc = -1,
log_q = log(c(0.00127899, 0.00242722, 1.27082e-06)),
log_pow = c(0, 0, 0),
qcreep = c(0.01, 0.01, 0.01),
log_galpha = log(c(6.89278, 4.98873)),
logit_gdiff = qlogis(c(0.351061, 0.229016)),
log_gshape = log(c(1.78938, 4.03173)),
log_gcv = log(c(0.428625, 0.731757)),
log_gobs = log(1.05115),
log_gdd = log(0),
log_sel50 = log(c(50.8998, 55.5424, 50.8998, 55.5424)),
log_sel95 = log(c(4.68424, 7.56075, 4.68424, 7.56075)),
log_seldome = log(c(1, 1, 1, 1)),
log_selmax = log(c(1, 0.5904295, 1, 0.724695)),
sel_sd50 = 0.1, sel_sd95 = 0.1, sel_sdmax = 0.1,
N0dev = rep(0, n_region * n_sex * n_length),
log_F = rep(log(0.02), length(data$catch_ytrf)),
Rdev_yr = array(c(
rep(0.5, 34),
0.443867, 0.13716, 0.309184, 0.374968, 0.162059, 0.0988453,
0.0256326, -0.1596, -0.138787, -0.144103, -0.137282, 0.220984, 0.573551,
0.766394, 0.57366, 0.188741, -0.205083, -0.452905, -0.412874, -0.403783,
-0.511264, -0.0934694, -0.0812145, -0.144924, 0.335367, -0.0463576,
-0.147839, -0.065563, -0.13226, -0.400957, 0.193345, -0.0966899, -0.379127,
-0.46455, -0.513768, -0.382686, 0.256375, -0.114664, -0.451873, -0.369532,
-0.362545, -0.789267, rep(-0.789267, 3)), dim = c(n_year, n_region))
)
priors <- list()
priors[["logit_U0"]] <- list(type = "normal", par1 = 0, par2 = 1.5, index = which("logit_U0" == names(parameters)))
priors[["log_Ma"]] <- list(type = "normal", par1 = log(0.12), par2 = 0.4, index = which("log_Ma" == names(parameters)))
priors[["log_galpha"]] <- list(type = "normal", par1 = log(2), par2 = 1.655705, index = which("log_galpha" == names(parameters))) # matches dnorm(2, 30)
priors[["logit_gdiff"]] <- list(type = "normal", par1 = -2, par2 = 1.5, index = which("logit_gdiff" == names(parameters)))
priors[["log_gshape"]] <- list(type = "normal", par1 = c(log(4.81), log(4.51)), par2 = c(0.2016748, 0.2142224), index = which("log_gshape" == names(parameters))) # matches dnorm(4.81, 1), dnorm(4.51, 1)
priors[["log_gcv"]] <- list(type = "normal", par1 = c(log(0.59), log(0.82)), par2 = c(0.439544, 0.3360953), index = which("log_gcv" == names(parameters))) # matches dnorm(0.59, 0.3), dnorm(0.82, 0.3)
priors[["log_gobs"]] <- list(type = "normal", par1 = log(1.48), par2 = 0.04990982, index = which("log_gobs" == names(parameters))) # matches dnorm(1.48, 0.074)
priors[["log_Rmu"]] <- list(type = "normal", par1 = log(32), par2 = 0.1, index = which("log_Rmu" == names(parameters)))
priors[["log_Rsd"]] <- list(type = "normal", par1 = log(2), par2 = 0.1, index = which("log_Rsd" == names(parameters)))
priors[["log_Rsigma"]] <- list(type = "normal", par1 = log(0.4), par2 = 0.3, index = which("log_Rsigma" == names(parameters)))
priors[["log_sel50"]] <- list(type = "normal", par1 = log(50), par2 = 0.5, index = which("log_sel50" == names(parameters))) # matches dnorm(50, 30)
priors[["log_sel95"]] <- list(type = "normal", par1 = log(10), par2 = 1.124519, index = which("log_sel95" == names(parameters))) # matches dnorm(10, 30)
priors[["log_selmax"]] <- list(type = "normal", par1 = log(1), par2 = 1, index = which("log_selmax" == names(parameters)))
data <- c(data, priors)Use TMB’s map option to turn parameters on/off:
map <- list()
map[["log_F0"]] <- factor(NA)
map[["logit_U0"]] <- factor(NA)
map[["logit_h"]] <- factor(rep(NA, length(parameters$logit_h)))
map[["N0dev"]] <- factor(rep(NA, length(parameters$N0dev)))
map[["log_Rsigma"]] <- factor(NA)
map[["Rrho"]] <- factor(NA)
map[["log_Rmu"]] <- factor(NA)
map[["log_Rsd"]] <- factor(NA)
map[["log_Ma"]] <- factor(c(1, 1))
map[["Mc"]] <- factor(rep(NA, length(parameters$Mc)))
map[["qcreep"]] <- factor(rep(NA, length(parameters$qcreep)))
# map[["log_pow"]] <- factor(c(NA, 1, 2))
map[["log_pow"]] <- factor(c(NA, NA, NA))
map[["log_gdd"]] <- factor(NA)
map[["log_sel50"]] <- factor(c(1, 2, 1, 2))
map[["log_sel95"]] <- factor(c(1, 2, 1, 2))
map[["log_selmax"]] <- factor(c(NA, 1, 2, 3))
map[["log_seldome"]] <- factor(c(NA, NA, NA, NA))
map[["sel_sd50"]] <- factor(NA)
map[["sel_sd95"]] <- factor(NA)
map[["sel_sdmax"]] <- factor(NA)
map[["log_F"]] <- factor(rep(NA, length(parameters$log_F)))Create the AD object using TMBs MakeADFun function from the TMB framework (Kristensen et al. 2016):
# check_data(data = data, parameters = parameters)
# obj <- MakeADFun(func = cmb(decamod, data), parameters = parameters, map = map, random = "Rdev_yr")
# data$rec_like <- 1 # 0=off, 1=RE, 2=RW, 3=AR1, 4=RE/RW, 5=RE/AR1
obj <- MakeADFun(func = cmb(decamod, data), parameters = parameters, map = map)Optimisation
List of parameters that are on:
The objective function value given the initial values:
obj$fn(obj$par)
#> [1] 12377.14Load the default parameter bounds:
Lwr <- rep(-Inf, length(obj$par))
Upr <- rep(Inf, length(obj$par))
Lwr[grep("log_R0", names(obj$par))] <- 1
Upr[grep("log_R0", names(obj$par))] <- 25
Lwr[grep("log_Ma", names(obj$par))] <- log(0.01)
Upr[grep("log_Ma", names(obj$par))] <- log(0.35)
Lwr[grep("log_q", names(obj$par))] <- -25
Upr[grep("log_q", names(obj$par))] <- 1
Lwr[grep("log_galpha", names(obj$par))] <- log(0)
Upr[grep("log_galpha", names(obj$par))] <- log(15)
Lwr[grep("logit_diff", names(obj$par))] <- qlogis(0.001)
Upr[grep("logit_diff", names(obj$par))] <- qlogis(0.99)
Lwr[grep("log_gshape", names(obj$par))] <- log(0.1)
Upr[grep("log_gshape", names(obj$par))] <- log(10)
Lwr[grep("log_gcv", names(obj$par))] <- log(0.01)
Upr[grep("log_gcv", names(obj$par))] <- log(2)
Lwr[grep("log_gobs", names(obj$par))] <- log(0.00001)
Upr[grep("log_gobs", names(obj$par))] <- log(5)
Lwr[grep("log_sel50", names(obj$par))] <- log(30)
Upr[grep("log_sel50", names(obj$par))] <- log(80)
Lwr[grep("log_sel95", names(obj$par))] <- log(1)
Upr[grep("log_sel95", names(obj$par))] <- log(50)
Lwr[grep("log_selmax", names(obj$par))] <- log(0.01)
Upr[grep("log_selmax", names(obj$par))] <- log(2)
Lwr[grep("Rdev_yr", names(obj$par))] <- -2.3
Upr[grep("Rdev_yr", names(obj$par))] <- 2.3
bounds <- data.frame(par = names(obj$par), lower = Lwr, value = obj$par, upper = Upr)Optimise using the nlminb function:
control <- list(eval.max = 100000, iter.max = 100000)
# opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, hessian = obj$he, control = control, lower = bounds$lower, upper = bounds$upper)
# opt <- nlminb(start = opt$par, objective = obj$fn, gradient = obj$gr, hessian = obj$he, control = control, lower = bounds$lower, upper = bounds$upper)
# check_estimability(obj = obj)
# save(data, parameters, bounds, map, obj, opt, file = "inst/extdata/opt.rda")
load(system.file("extdata", "opt.rda", package = "decamod"))
obj_orig <- obj
obj <- MakeADFun(func = cmb(decamod, data), parameters = parameters, map = map)
obj$par <- obj_orig$env$last.par.best
obj$env$last.par <- obj_orig$env$last.par.best
obj$env$last.par.best <- obj_orig$env$last.par.best
sim_check <- obj$simulate(par = obj$env$last.par.best)
sim_check$cpue_obs
#> [1] 1.1436197 1.4143439 1.0778170 1.6263713 1.1717140 1.4174619 1.1355674
#> [8] 1.0988176 0.7998008 1.1697327 0.8464932 1.3143685 1.2816506 0.8069272
#> [15] 0.6859810 0.9561318 0.8202065 0.8317777 0.5971988 0.8904222 0.7434740
#> [22] 0.7554992 1.0908253 1.2520284 1.3308942 0.8422874 0.8893723 1.0971338
#> [29] 0.8725407 1.3226972 1.2538763 1.5431672 1.4525257 1.5361200 2.3895805
#> [36] 1.4570919 1.8474460 1.6507856 2.1903790 1.0220731 1.5952363 1.2587090
#> [43] 1.1138957 1.1907391 0.9728322 0.7661218 0.7107993 0.6667983 0.9700323
#> [50] 0.9779460 0.9262793 0.6524625 0.9344182 0.9087093 0.6246443 0.7767931
#> [57] 1.0260233 1.2695448 0.7067377 0.7465646 0.9137000 0.6619899 0.8062588
#> [64] 0.5885949 1.2307322 0.5474687 0.6137540 0.7647464 0.8873836 0.6250065
#> [71] 0.6323406 0.5034142 0.5484842 0.6032764 0.6895966 0.6632288 0.9263637
#> [78] 1.2626022 1.6390232 1.4358115 0.7731432 1.1966847 1.2681892 1.4974291
#> [85] 1.2328537 1.5855560 1.2062674 2.4050799 1.6838910 1.3034910 0.9949625
#> [92] 1.0989364 1.0956776 1.5468076 1.1765856 1.6813467 0.6766573 0.9418232
#> [99] 0.7107431 0.7135743 0.6877821 0.8026488 0.8024451 0.9855816 0.4822100
#> [106] 1.1881128 0.6865255 0.9907488 0.8277532 1.0619371 0.8067730 0.9512680
#> [113] 0.7606435 0.9746290 1.0526323 0.9096840 0.9145889 1.3549585 0.5455945
#> [120] 0.9853760 0.6706854 0.6478939 0.6272543 0.6179945 0.4172797 0.7075730
#> [127] 0.5305094 0.7886473 0.7367754 0.8987394 0.6039221 0.8119592 1.2404362
#> [134] 0.9497365 1.2030458 1.3021654 1.2337470 1.4475957 1.3919887 2.1671320
#> [141] 1.4075028 2.2312320
sim_check$sexr_obs2
#> [1] 0.474401055 0.087457565 0.393509710 -0.165530715 0.546000209
#> [6] -0.009425154 0.282345309 0.101748119 0.293531497 -0.357501983
#> [11] -0.279755083 -0.334873857 -0.165246676 -0.496356967 -0.127422790
#> [16] -0.441823474 -0.003130806 -0.681917870 -0.425570484 -0.667469081
#> [21] 0.084243816 -0.259924072 0.113932394 -0.197644296 0.262788546
#> [26] -0.339261213 0.179793907 -0.101699325 0.153510932 -0.184047720
#> [31] 0.350354187 -0.200806254 0.260033524 -0.156533346 0.160754413
#> [36] -0.029474411 -0.080488908 -0.520912287 0.209093505 -0.381410818
#> [41] 0.255242852 -0.244519308 -0.222355390 -0.223425286 0.046831818
#> [46] -0.452434676 -0.085876599 -0.140466270 0.288777304 -0.070758859
#> [51] 0.453464805 -0.030812570 0.290304353 0.179895753 0.299678753
#> [56] 0.074203738 0.424729933 -0.104916891 0.350577174 -0.083020880
#> [61] 0.346212712 0.119644845
sim_check$lf_obs2
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
#> [1,] 1 1 2 3 9 13 8 5 6 7 5 6 5
#> [2,] 0 1 2 6 3 10 12 7 13 9 8 3 3
#> [3,] 2 1 6 8 19 11 9 12 6 5 3 5 4
#> [4,] 0 0 4 5 7 7 17 13 10 11 15 8 7
#> [5,] 0 2 5 12 9 18 16 16 12 5 7 3 3
#> [6,] 1 2 3 2 11 11 19 13 17 10 13 5 6
#> [7,] 1 1 8 9 14 13 16 11 13 4 7 6 1
#> [8,] 0 1 1 7 7 18 16 14 15 11 10 9 3
#> [9,] 0 2 3 3 12 20 11 8 10 9 8 2 3
#> [10,] 0 1 1 2 5 3 10 9 10 10 8 1 5
#> [11,] 0 0 3 7 8 13 8 11 7 7 8 0 3
#> [12,] 2 0 2 3 1 5 10 6 7 8 7 9 5
#> [13,] 0 2 1 9 9 11 10 7 15 14 8 2 2
#> [14,] 0 0 0 4 3 5 7 17 10 5 7 8 6
#> [15,] 1 2 2 2 6 15 13 4 5 8 6 6 4
#> [16,] 0 1 0 4 4 3 5 10 3 11 12 9 4
#> [17,] 1 0 3 8 4 10 6 12 11 7 3 3 5
#> [18,] 1 0 0 1 1 7 3 7 10 6 12 8 2
#> [19,] 0 0 1 5 4 6 3 6 7 6 4 6 6
#> [20,] 0 1 0 1 0 2 2 3 9 6 6 6 6
#> [21,] 0 0 2 4 4 10 10 8 6 6 8 6 2
#> [22,] 1 0 2 1 3 5 10 6 11 7 4 6 5
#> [23,] 1 0 0 6 4 6 7 8 6 5 6 7 10
#> [24,] 1 0 0 1 0 2 9 5 8 7 9 3 7
#> [25,] 0 0 1 7 6 9 9 8 9 5 12 4 4
#> [26,] 0 0 0 1 1 8 9 12 9 12 5 9 9
#> [27,] 0 0 2 7 8 3 8 7 5 4 4 1 6
#> [28,] 0 0 1 2 3 4 7 6 9 12 12 14 5
#> [29,] 0 0 2 5 7 10 7 11 12 4 4 3 4
#> [30,] 0 0 0 2 3 5 6 3 15 10 8 12 6
#> [31,] 0 0 2 8 14 6 10 11 6 3 3 4 2
#> [32,] 0 0 1 2 5 5 9 7 9 8 11 8 7
#> [33,] 0 1 3 3 11 11 8 6 5 4 4 7 1
#> [34,] 1 0 1 2 2 7 9 9 9 10 5 7 5
#> [35,] 1 1 2 7 10 5 11 7 3 3 7 3 2
#> [36,] 0 1 1 3 3 5 10 5 7 9 2 8 4
#> [37,] 0 2 1 6 10 11 9 8 8 9 2 0 4
#> [38,] 0 2 1 4 5 8 9 7 10 5 14 9 3
#> [39,] 0 1 4 2 6 4 10 8 10 4 9 5 4
#> [40,] 1 0 3 1 5 2 6 6 9 14 9 3 9
#> [41,] 0 1 5 10 10 10 11 5 3 6 5 3 0
#> [42,] 0 0 4 1 3 7 11 8 11 8 8 5 3
#> [43,] 3 1 1 5 10 11 8 9 9 6 6 3 3
#> [44,] 0 1 0 3 7 3 11 14 12 10 5 8 4
#> [45,] 0 2 2 5 9 13 12 9 8 4 2 1 3
#> [46,] 0 1 1 2 7 8 4 13 7 5 6 5 2
#> [47,] 0 2 1 8 12 14 5 4 7 7 8 2 3
#> [48,] 0 2 1 1 3 6 4 14 10 8 7 2 5
#> [49,] 0 0 3 4 10 18 10 10 7 6 3 4 4
#> [50,] 0 0 2 7 6 18 13 12 5 8 11 4 1
#> [51,] 2 1 2 6 11 16 13 8 11 3 4 6 3
#> [52,] 2 1 2 4 7 9 8 11 12 15 7 7 7
#> [53,] 1 2 1 9 13 13 15 10 7 7 4 3 3
#> [54,] 0 0 3 5 7 11 13 11 8 8 8 5 4
#> [55,] 1 2 2 6 11 17 14 5 11 9 3 1 1
#> [56,] 2 0 0 5 14 8 9 8 10 12 8 5 6
#> [57,] 0 1 3 5 7 14 11 9 7 4 7 2 2
#> [58,] 0 2 1 5 8 10 4 12 10 11 6 3 4
#> [59,] 2 0 1 8 7 7 3 7 9 9 11 4 3
#> [60,] 0 0 1 2 10 3 14 11 10 11 7 5 4
#> [61,] 0 0 1 7 11 10 11 11 8 2 2 1 2
#> [62,] 1 0 1 2 3 5 9 11 6 7 12 6 3
#> [63,] 0 0 3 8 12 8 7 11 6 4 5 3 1
#> [64,] 1 0 0 2 5 9 11 14 15 9 7 5 1
#> [65,] 1 1 2 7 13 15 9 8 3 9 2 1 1
#> [66,] 2 0 1 1 5 7 6 10 8 15 11 6 9
#> [67,] 0 0 3 9 9 15 9 11 5 3 2 6 4
#> [68,] 0 0 1 3 4 9 9 10 12 14 11 12 3
#> [69,] 0 1 1 2 8 18 9 12 7 5 3 4 3
#> [70,] 0 1 3 3 3 8 7 13 13 8 8 4 2
#> [71,] 0 1 6 7 12 9 7 6 10 9 6 2 1
#> [72,] 0 2 2 1 3 6 8 11 15 4 9 9 8
#> [73,] 0 2 3 4 13 12 5 7 9 6 4 2 1
#> [74,] 0 0 4 3 2 9 11 10 12 8 11 6 8
#> [75,] 0 0 3 7 11 17 9 8 5 6 4 4 0
#> [76,] 1 1 2 2 3 9 11 15 7 10 6 6 9
#> [77,] 0 3 5 7 14 17 8 18 5 7 3 0 2
#> [78,] 0 2 2 8 6 13 11 10 12 10 10 5 3
#> [79,] 0 2 1 7 18 17 12 14 9 6 5 6 2
#> [80,] 0 1 2 3 2 6 10 23 20 11 9 5 5
#> [81,] 0 2 1 7 13 20 16 8 6 7 4 2 0
#> [82,] 0 1 2 5 2 10 9 12 13 9 9 6 5
#> [83,] 0 1 6 5 19 9 16 8 9 10 4 2 1
#> [84,] 1 2 2 4 7 13 10 10 19 18 7 8 3
#> [85,] 1 2 2 6 8 12 14 13 10 11 3 3 2
#> [86,] 0 0 4 4 2 6 12 7 11 10 12 10 6
#> [87,] 4 0 4 11 15 17 8 14 8 5 4 1 0
#> [88,] 2 0 0 2 7 11 11 10 7 16 9 4 7
#> [89,] 1 0 1 5 8 8 9 11 10 3 4 5 1
#> [90,] 0 1 0 4 5 6 14 12 15 4 3 7 5
#> [91,] 0 3 2 3 11 5 12 9 5 10 4 1 2
#> [92,] 1 0 0 3 5 8 10 9 10 14 8 4 5
#> [93,] 0 2 2 5 7 11 11 4 11 2 3 3 1
#> [94,] 0 0 2 1 2 4 6 18 11 10 6 6 2
#> [95,] 0 3 0 10 7 9 8 4 9 8 3 3 2
#> [96,] 1 1 1 0 5 6 10 9 14 10 11 3 2
#> [97,] 0 2 6 6 7 7 8 12 3 7 6 1 2
#> [98,] 1 0 2 0 12 4 8 13 15 8 2 2 6
#> [99,] 1 1 2 7 10 10 14 8 7 3 2 2 1
#> [100,] 0 0 0 3 5 6 7 12 12 5 16 5 6
#> [101,] 0 1 4 7 10 12 6 12 5 3 5 7 0
#> [102,] 0 0 0 5 2 10 18 7 6 7 6 5 6
#> [103,] 0 1 3 5 13 12 9 6 10 3 3 1 0
#> [104,] 1 0 1 1 2 7 5 9 17 12 5 6 5
#> [105,] 0 0 1 4 7 12 12 4 5 6 3 2 2
#> [106,] 0 0 0 3 4 3 5 10 10 10 3 0 4
#> [107,] 0 0 1 4 5 8 12 8 4 6 3 5 3
#> [108,] 0 0 0 1 7 6 6 8 10 2 6 2 3
#> [109,] 1 0 0 5 11 6 10 6 10 4 4 5 4
#> [110,] 0 0 1 1 7 6 8 6 10 13 8 4 2
#> [111,] 0 1 1 7 7 12 8 8 10 6 5 6 2
#> [112,] 0 0 1 3 3 11 11 9 14 6 10 9 2
#> [113,] 0 0 3 4 8 8 5 6 9 5 6 3 3
#> [114,] 0 1 1 3 1 5 9 5 9 7 7 3 2
#> [115,] 0 0 1 1 5 6 7 9 11 12 4 4 2
#> [116,] 0 0 1 1 2 9 3 8 7 5 6 6 7
#> [117,] 0 0 1 1 6 4 5 3 2 3 4 5 4
#> [118,] 0 0 2 1 3 5 2 6 9 4 9 3 2
#> [119,] 0 0 0 4 2 9 6 3 2 7 4 6 5
#> [120,] 1 0 0 1 3 2 4 8 9 8 7 1 2
#> [121,] 0 2 3 4 4 5 5 5 5 12 4 3 3
#> [122,] 0 0 2 1 3 3 4 7 3 7 6 7 6
#> [123,] 0 0 3 0 5 8 3 8 4 4 5 4 3
#> [124,] 0 0 0 0 4 3 10 5 5 3 8 3 5
#> [,14] [,15] [,16] [,17] [,18] [,19] [,20]
#> [1,] 1 1 2 1 1 0 0
#> [2,] 3 1 3 0 1 1 2
#> [3,] 2 3 1 2 1 0 1
#> [4,] 2 4 2 2 1 0 1
#> [5,] 3 0 0 2 0 1 1
#> [6,] 2 4 4 2 2 0 1
#> [7,] 2 1 0 0 0 0 0
#> [8,] 3 2 1 1 0 1 1
#> [9,] 0 1 1 0 0 0 0
#> [10,] 6 4 2 2 1 1 1
#> [11,] 2 3 0 0 1 0 0
#> [12,] 4 1 2 1 0 1 2
#> [13,] 3 2 0 0 1 0 1
#> [14,] 6 2 1 2 0 1 0
#> [15,] 2 2 1 2 0 0 0
#> [16,] 4 4 2 0 0 0 0
#> [17,] 4 1 3 0 0 1 0
#> [18,] 6 1 2 1 0 0 0
#> [19,] 3 2 6 1 2 0 1
#> [20,] 6 0 0 1 1 0 0
#> [21,] 2 1 4 1 3 1 3
#> [22,] 4 1 3 1 0 1 1
#> [23,] 2 3 0 1 0 0 1
#> [24,] 9 5 1 1 0 1 0
#> [25,] 4 3 1 3 1 0 2
#> [26,] 8 3 5 3 3 0 1
#> [27,] 6 1 4 1 0 0 4
#> [28,] 3 2 3 0 1 0 1
#> [29,] 1 3 0 2 2 2 5
#> [30,] 5 6 8 3 0 1 1
#> [31,] 2 1 1 2 1 1 0
#> [32,] 3 2 0 4 3 2 0
#> [33,] 1 2 1 0 2 1 3
#> [34,] 5 4 3 3 1 0 0
#> [35,] 0 0 1 2 3 0 0
#> [36,] 5 9 3 3 2 0 2
#> [37,] 2 1 1 1 1 1 1
#> [38,] 9 3 0 1 1 0 0
#> [39,] 3 1 0 1 3 1 1
#> [40,] 3 5 7 3 3 0 1
#> [41,] 0 1 2 0 2 0 1
#> [42,] 5 4 4 3 1 1 1
#> [43,] 2 2 0 0 0 1 0
#> [44,] 4 1 2 3 2 0 1
#> [45,] 0 1 0 0 0 2 0
#> [46,] 4 1 5 2 2 1 1
#> [47,] 0 2 1 0 0 1 1
#> [48,] 6 3 1 2 1 0 2
#> [49,] 3 1 2 1 0 1 0
#> [50,] 3 2 2 1 0 2 1
#> [51,] 2 4 1 1 0 0 0
#> [52,] 3 2 5 2 0 0 3
#> [53,] 1 0 0 0 0 0 0
#> [54,] 3 1 5 1 0 0 0
#> [55,] 6 1 0 1 0 0 0
#> [56,] 3 1 2 1 0 2 0
#> [57,] 1 0 1 1 1 0 0
#> [58,] 2 1 0 1 1 0 0
#> [59,] 1 1 0 0 0 0 1
#> [60,] 2 0 0 0 0 1 0
#> [61,] 1 1 0 1 1 0 0
#> [62,] 5 2 3 1 0 0 0
#> [63,] 6 2 0 1 0 1 0
#> [64,] 4 1 1 0 0 0 0
#> [65,] 1 1 0 1 0 0 1
#> [66,] 6 1 0 0 0 0 0
#> [67,] 2 1 0 0 2 0 0
#> [68,] 4 0 3 1 0 1 0
#> [69,] 2 1 1 0 1 2 0
#> [70,] 2 2 4 1 0 0 1
#> [71,] 2 1 2 1 0 0 0
#> [72,] 3 3 1 2 0 1 2
#> [73,] 4 0 1 0 2 1 0
#> [74,] 0 4 1 0 0 0 0
#> [75,] 2 0 2 1 0 0 0
#> [76,] 4 3 2 3 0 0 0
#> [77,] 0 2 0 0 0 0 1
#> [78,] 2 1 0 0 1 0 0
#> [79,] 2 2 0 0 0 0 0
#> [80,] 2 8 1 1 0 0 0
#> [81,] 1 0 0 0 0 0 1
#> [82,] 5 7 1 4 0 0 0
#> [83,] 3 3 0 0 0 0 1
#> [84,] 2 2 4 0 1 0 0
#> [85,] 0 1 1 0 1 0 0
#> [86,] 2 2 1 2 2 0 0
#> [87,] 1 1 0 1 0 0 0
#> [88,] 2 8 1 2 3 1 0
#> [89,] 0 0 0 0 0 0 0
#> [90,] 2 2 0 0 0 0 0
#> [91,] 0 1 0 0 0 0 0
#> [92,] 2 2 2 0 0 0 0
#> [93,] 3 1 0 0 0 0 1
#> [94,] 2 1 0 2 0 1 0
#> [95,] 1 1 1 0 0 0 0
#> [96,] 1 0 1 0 0 0 1
#> [97,] 2 0 0 0 0 0 0
#> [98,] 4 1 1 0 0 0 0
#> [99,] 1 2 0 0 0 0 0
#> [100,] 4 0 2 0 0 0 0
#> [101,] 1 0 0 0 0 0 1
#> [102,] 4 1 1 0 0 0 0
#> [103,] 0 0 0 1 1 0 0
#> [104,] 3 0 1 0 0 0 0
#> [105,] 2 0 1 1 0 0 0
#> [106,] 0 0 2 0 1 0 0
#> [107,] 0 0 1 0 0 0 0
#> [108,] 2 0 2 1 0 0 0
#> [109,] 3 2 2 0 0 0 0
#> [110,] 5 5 1 0 0 0 0
#> [111,] 2 1 3 2 0 0 0
#> [112,] 2 1 3 1 0 0 0
#> [113,] 2 1 0 1 0 1 0
#> [114,] 3 1 1 1 2 0 0
#> [115,] 5 1 4 0 0 1 0
#> [116,] 6 2 2 1 0 0 0
#> [117,] 2 3 0 0 1 0 1
#> [118,] 0 2 2 0 1 0 0
#> [119,] 0 0 4 0 0 1 1
#> [120,] 5 4 1 1 1 0 0
#> [121,] 1 4 2 3 0 0 1
#> [122,] 4 6 1 1 0 2 0
#> [123,] 1 5 3 1 1 0 2
#> [124,] 6 4 4 1 0 1 0
plogis(sim_check$sexr_obs2)
#> [1] 0.6164249 0.5218505 0.5971273 0.4587116 0.6332071 0.4976437 0.5701211
#> [8] 0.5254151 0.5728605 0.4115644 0.4305138 0.4170552 0.4587821 0.3783972
#> [15] 0.4681873 0.3913066 0.4992173 0.3358334 0.3951846 0.3390638 0.5210485
#> [22] 0.4353824 0.5284523 0.4507491 0.5653217 0.4159889 0.5448278 0.4745971
#> [29] 0.5383025 0.4541175 0.5867035 0.4499664 0.5646445 0.4609464 0.5401023
#> [36] 0.4926319 0.4798886 0.3726389 0.5520838 0.4057867 0.5634665 0.4391729
#> [43] 0.4446391 0.4443749 0.5117058 0.3887821 0.4785440 0.4649411 0.5716968
#> [50] 0.4823177 0.6114627 0.4922975 0.5720706 0.5448530 0.5743640 0.5185424
#> [57] 0.6046145 0.4737948 0.5867575 0.4792567 0.5856989 0.5298756
plot_lfs(data, obj)
# oneStepPredict(obj, observation.name = "sexr_obs2", method = "oneStepGeneric", trace = FALSE)
res1 <- oneStepPredict(obj, observation.name = "lf_obs2", method = "cdf", discrete = TRUE, trace = FALSE)
# res <- oneStepPredict(obj = obj, observation.name = "lf_obs2", method = "oneStepGeneric", discrete = TRUE, discreteSupport = 0:55, trace = FALSE)
# res
# chk <- checkConsistency(obj, par = c("log_R0", "log_Ma", "log_q", "log_sel50", "log_sel95", "log_selmax", "Rdev_yr"))
# chk <- checkConsistency(obj, par = "log_Ma")
# chk <- checkConsistency(obj, estimate = TRUE)
# chk
# summary(chk)
exp(obj$par)[names(obj$par) %in% c("log_R0", "log_Ma", "log_pow")]
#> log_R0 log_Ma
#> 5.845468e+05 1.982332e-01
data.frame(par = names(obj$par), val = obj$par, gr = obj$gr()[1,]) %>%
arrange(-gr) %>%
head()
#> outer mgc: 0.0001895574
#> par val gr
#> 1 log_R0 13.2785921 1.895574e-04
#> 2 log_galpha 1.9010120 1.197014e-04
#> 3 log_galpha 1.6107058 5.351254e-05
#> 4 log_sel50 3.9288348 4.382421e-05
#> 5 log_gshape 0.9814145 1.193615e-05
#> 6 log_gcv -0.3052584 1.039047e-05
Calculate standard deviations of all model parameters, including non linear functions of random effects and parameters specified through the ADREPORT() macro from the user template. Check that the Hessian is positive definite. Try the analytical Hessian:
# Report <- sdreport(obj)
# print(Report$pdHess)
# he <- obj$he()
# he_inv <- solve(he)
# he_ch <- chol(he)
# ev <- eigen(he)
# range(ev$values)Model weights
Francis length-frequency iterative re-weighting (Francis 2011).
names(tags)
#> [1] "Event" "Sex" "RelMonth" "RelYear"
#> [5] "RecMonth" "RecYear" "DaysatLiberty" "gi_L1"
#> [9] "gi_L2" "TWConvertRel" "TWConvertRec" "TagRelease"
#> [13] "RelStatArea" "Condition" "Tagtype" "RecStatArea"
#> [17] "gi_obs" "gi_liberty" "gi_sex" "gi_morph"
plot_tag_residuals(data = data, obj = obj, type = "Pearson", xvar = "gi_L1")
#> # A tibble: 1 × 2
#> SDNR MAR
#> <dbl> <dbl>
#> 1 1.02 0.620
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
# plot_tag_residuals(data = data, obj = obj, type = "Pearson", xvar = "gi_liberty")
plot_tag_residuals(data = data, obj = obj, type = "Pearson", xvar = "RelStatArea")
#> # A tibble: 1 × 2
#> SDNR MAR
#> <dbl> <dbl>
#> 1 1.02 0.620
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_tag_residuals(data = data, obj = obj, type = "Pearson", xvar = "Condition")
#> # A tibble: 1 × 2
#> SDNR MAR
#> <dbl> <dbl>
#> 1 1.02 0.620
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_tag_residuals(data = data, obj = obj, type = "Pearson", xvar = "Tagtype")
#> # A tibble: 1 × 2
#> SDNR MAR
#> <dbl> <dbl>
#> 1 1.02 0.620
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_tag_residuals(data = data, obj = obj, type = "Pearson", xvar = "RelMonth")
#> # A tibble: 1 × 2
#> SDNR MAR
#> <dbl> <dbl>
#> 1 1.02 0.620
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
get_lf_Francis(data = data, obj = obj)
#> # A tibble: 2 × 7
#> # Groups: region, sex [2]
#> region sex sd n Francis_mult orig_wt new_wt
#> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
#> 1 1 1 1.21 62 0.696 1.20 0.838
#> 2 1 2 1.16 62 0.752 0.982 0.738
plot_cpue_residuals(data = data, obj = obj, type = "PIT")
#> # A tibble: 3 × 3
#> q SDNR MAR
#> <fct> <dbl> <dbl>
#> 1 q1 0.797 0.616
#> 2 q2 0.857 0.451
#> 3 q3 0.982 0.775
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_cpue_residuals(data = data, obj = obj, type = "OSA")
#> # A tibble: 3 × 3
#> q SDNR MAR
#> <fct> <dbl> <dbl>
#> 1 q1 0.797 0.616
#> 2 q2 0.857 0.451
#> 3 q3 0.982 0.775
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_sexr_residuals(data = data, obj = obj, type = "PIT")
#> SDNR MAR
#> 1 1.074625 0.5652161
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
# plot_sexr_residuals(data = data, obj = obj, type = "OSA")
Model fits
plot_growth_increment(data = data, parameters = parameters, obj = obj, tau = 1)
plot_growth_density(data = data, obj = obj)
nsim <- 5
plot_catch(data = data, obj = obj, nsim = nsim)
plot_catch_residuals(data = data, obj = obj)
plot_cpue(data = data, obj = obj, nsim = nsim)
plot_cpue_residuals(data = data, obj = obj, type = "PIT")
#> # A tibble: 3 × 3
#> q SDNR MAR
#> <fct> <dbl> <dbl>
#> 1 q1 0.797 0.616
#> 2 q2 0.857 0.451
#> 3 q3 0.982 0.775
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
Ideally, Probability Integral Transform (PIT) residuals should appear as a roughly uniform distribution when plotted on a histogram, meaning the points should be scattered randomly around a horizontal line with no clear patterns, indicating that the model is accurately capturing the data’s probability distribution across different values of the response variable.
When plotted on a histogram, the distribution of PIT residuals should resemble a uniform distribution between 0 and 1, with no significant clustering at any particular point.
A scatter plot of PIT residuals against fitted values should not show any discernible patterns, like funneling or fanning, which would suggest issues with heteroscedasticity (unequal variance). Large deviations from the uniform distribution, especially at the extreme ends, could indicate potential outliers or model misspecification.
plot_sexr(data = data, obj = obj, nsim = 15)
# plot_sexr_residuals(data = data, obj = obj, type = "Raw")
# plot_sexr_residuals(data = data, obj = obj, type = "old")
plot_sexr_residuals(data = data, obj = obj, type = "PIT")
#> SDNR MAR
#> 1 1.074625 0.5652161
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_lfs(data = data, obj = obj, nsim = nsim)
# plot_Fs(data = data, obj = obj)
plot_Us(data = data, obj = obj)
One step ahead residuals
Plot the one step ahead (OSA) residuals for the LFs using the R package compResidual (https://github.com/fishfollower/compResidual):
# TMB:::install.contrib("https://github.com/vtrijoulet/OSA_multivariate_dists/archive/main.zip")
# remotes::install_github("fishfollower/compResidual/compResidual", force=TRUE)
# library(compResidual)
# plot_lf_residuals(data = data, obj = obj, pick_one = 1)
# plot_lf_residuals(data = data, obj = obj, pick_one = 6)
# plot_lf_residuals(data = data, obj = obj, pick_one = 7)
# plot_lf_residuals(data = data, obj = obj, pick_one = 8)
osa_cpue <- oneStepPredict(obj = obj, observation.name = "cpue_obs")
#> [1] 142
#> [1] 141
#> [1] 140
#> [1] 139
#> [1] 138
#> [1] 137
#> [1] 136
#> [1] 135
#> [1] 134
#> [1] 133
#> [1] 132
#> [1] 131
#> [1] 130
#> [1] 129
#> [1] 128
#> [1] 127
#> [1] 126
#> [1] 125
#> [1] 124
#> [1] 123
#> [1] 122
#> [1] 121
#> [1] 120
#> [1] 119
#> [1] 118
#> [1] 117
#> [1] 116
#> [1] 115
#> [1] 114
#> [1] 113
#> [1] 112
#> [1] 111
#> [1] 110
#> [1] 109
#> [1] 108
#> [1] 107
#> [1] 106
#> [1] 105
#> [1] 104
#> [1] 103
#> [1] 102
#> [1] 101
#> [1] 100
#> [1] 99
#> [1] 98
#> [1] 97
#> [1] 96
#> [1] 95
#> [1] 94
#> [1] 93
#> [1] 92
#> [1] 91
#> [1] 90
#> [1] 89
#> [1] 88
#> [1] 87
#> [1] 86
#> [1] 85
#> [1] 84
#> [1] 83
#> [1] 82
#> [1] 81
#> [1] 80
#> [1] 79
#> [1] 78
#> [1] 77
#> [1] 76
#> [1] 75
#> [1] 74
#> [1] 73
#> [1] 72
#> [1] 71
#> [1] 70
#> [1] 69
#> [1] 68
#> [1] 67
#> [1] 66
#> [1] 65
#> [1] 64
#> [1] 63
#> [1] 62
#> [1] 61
#> [1] 60
#> [1] 59
#> [1] 58
#> [1] 57
#> [1] 56
#> [1] 55
#> [1] 54
#> [1] 53
#> [1] 52
#> [1] 51
#> [1] 50
#> [1] 49
#> [1] 48
#> [1] 47
#> [1] 46
#> [1] 45
#> [1] 44
#> [1] 43
#> [1] 42
#> [1] 41
#> [1] 40
#> [1] 39
#> [1] 38
#> [1] 37
#> [1] 36
#> [1] 35
#> [1] 34
#> [1] 33
#> [1] 32
#> [1] 31
#> [1] 30
#> [1] 29
#> [1] 28
#> [1] 27
#> [1] 26
#> [1] 25
#> [1] 24
#> [1] 23
#> [1] 22
#> [1] 21
#> [1] 20
#> [1] 19
#> [1] 18
#> [1] 17
#> [1] 16
#> [1] 15
#> [1] 14
#> [1] 13
#> [1] 12
#> [1] 11
#> [1] 10
#> [1] 9
#> [1] 8
#> [1] 7
#> [1] 6
#> [1] 5
#> [1] 4
#> [1] 3
#> [1] 2
#> [1] 1
dim(osa_cpue)
#> [1] 142 5
length(data$cpue_obs)
#> [1] 142
plot_cpue_residuals(data = data, obj = obj, type = "PIT")
#> # A tibble: 3 × 3
#> q SDNR MAR
#> <fct> <dbl> <dbl>
#> 1 q1 0.797 0.616
#> 2 q2 0.857 0.451
#> 3 q3 0.982 0.775
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
plot_cpue_residuals(data = data, obj = obj, type = "OSA")
#> # A tibble: 3 × 3
#> q SDNR MAR
#> <fct> <dbl> <dbl>
#> 1 q1 0.797 0.616
#> 2 q2 0.857 0.451
#> 3 q3 0.982 0.775
#> [1] "SDNR should be close to 1 and MAR should be close to 0.67"
Projections
# i1 <- which(c(first_year:last_year) == 1979)
# i1 <- which(c(first_year:last_year) == min(lf$fyear))
# i1 <- 1
# sim <- project_recruitment(object = obj, h = 5, y1 = NULL, y2 = NULL)
#
# library(forecast)
# i1 <- which(c(first_year:last_year) == 1979)
# i1 <- which(c(first_year:last_year) == min(lf$fyear))
# i1 <- 1
# Rdevs <- obj$report(obj$env$last.par.best)$Rdev_yr[i1:n_year,1]
# fit <- auto.arima(y = Rdevs, max.p = 0, max.d = 5, max.q = 5, approximation = FALSE, stepwise = FALSE, ic = "bic", trace = TRUE)
# fit <- auto.arima(y = Rdevs, max.p = 5, max.d = 0, max.q = 5, approximation = FALSE, stepwise = FALSE, ic = "bic", trace = TRUE)
# fit <- auto.arima(y = Rdevs, max.p = 5, max.d = 5, max.q = 0, approximation = FALSE, stepwise = FALSE, ic = "bic", trace = TRUE)
# fit <- auto.arima(y = Rdevs, approximation = FALSE, stepwise = FALSE, ic = "bic", trace = TRUE)
# checkresiduals(fit)
# autoplot(fit)
# autoplot(forecast(fit, h = 5))
# autoplot(simulate(object = fit, nsim = 5, future = TRUE, bootstrap = TRUE))
#
# fc <- forecast(object = fit, h = 5)
#
# df1 <- data.frame(x = 1:length(Rdevs), mean = Rdevs, group = 1)
# df2 <- data.frame(x = 80:84, mean = fc$mean, lower = fc$lower, upper = fc$upper, group = 2)
# df <- bind_rows(df1, df2)
#
# p <- ggplot(data = df, aes(x = x, y = mean, group = group)) +
# geom_ribbon(aes(ymin = lower.95., ymax = upper.95.), alpha = 0.5, color = NA, fill = "lightblue") +
# geom_ribbon(aes(ymin = lower.80., ymax = upper.80.), alpha = 0.5, color = NA, fill = "blue") +
# geom_line(color = "blue")
#
# for (i in 1:10) {
# sim <- simulate(object = fit, nsim = 5, future = TRUE, bootstrap = TRUE)
# dsim <- data.frame(x = 80:84, mean = sim, group = 3)
# p <- p + geom_line(data = dsim, size = 0.3)
# }
#
# pA likelihood profile (Hilborn and Walters 1992) for log_Ma can be done using:
# prof <- TMB::tmbprofile(obj = obj, name = "Rdev_yr[65,1]")
# lincomb <- as.numeric(names(obj$par) == "log_Ma")
# prof <- TMB::tmbprofile(obj = obj, lincomb = lincomb)
# plot(prof)
# confint(prof)
# df <- likelihood_profile(obj = obj, name = "log_Ma")
# lincomb <- rep(0, length(obj$par))
# names(obj$par)
# names(obj$par)[68]
# lincomb[68] <- 1
# prof_mle <- tmbprofile2(obj = obj, lincomb = lincomb)
# # prof_mcmc <- get_mcmc_profile(obj = obj, mcmc = mcmc1, name = "log_galpha[1]")
# names(obj$par)[names(obj$par) == "Rdev_yr"]
# prof_mcmc <- get_mcmc_profile(obj = obj, mcmc = mcmc1, name = "Rdev_yr[42]")
#
# plot_profile(x = prof_mle, lab = "log(Galpha)")
# ggsave(file = "like_prof1.png", width = 9)
# plot_profile(x = prof_mle, y = prof_mcmc, lab = "log(M)", rescale = FALSE)
# ggsave(file = "like_prof2.png", width = 9)
# plot_profile(x = prof_mle, y = prof_mcmc, lab = "log(M)")
# ggsave(file = "like_prof3.png", width = 9)
#
# prof_mle <- tmbprofile2(obj = obj, name = "log_Ma")
# prof_mcmc <- get_mcmc_profile(obj = obj, mcmc = mcmc1, name = "log_Ma")
#
# plot_profile(x = prof_mle, lab = "log(M)")
# ggsave(file = "like_prof1.png", width = 9)
# plot_profile(x = prof_mle, y = prof_mcmc, lab = "log(M)", rescale = FALSE)
# ggsave(file = "like_prof2.png", width = 9)
# plot_profile(x = prof_mle, y = prof_mcmc, lab = "log(M)")
# ggsave(file = "like_prof3.png", width = 9)
Fournier, David, and C. P. Archibald. 1982. “A General Theory for Analyzing Catch at Age Data.” Canadian Journal of Fisheries and Aquatic Sciences 39 (8): 1195–207. https://doi.org/10.1139/f82-157.
Francis, R. I. C. Chris. 2011. “Data Weighting in Statistical Fisheries Stock Assessment Models.” Canadian Journal of Fisheries and Aquatic Sciences 68 (6): 1124–38. https://doi.org/10.1139/f2011-025.
Hilborn, Ray, and Carl J. Walters. 1992. Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty. Chapman; Hall. https://doi.org/10.1007/978-1-4615-3598-0.
Kristensen, Kasper, Anders Nielsen, Casper W. Berg, Hans Skaug, and Bradley M. Bell. 2016. “TMB: Automatic Differentiation and Laplace Approximation.” Journal of Statistical Software 70 (5): 1–21. https://doi.org/10.18637/jss.v070.i05.
Maunder, Mark N., and André E. Punt. 2004. “Standardizing Catch and Effort Data: A Review of Recent Approaches.” Fisheries Research 70 (2–3): 141–59. https://doi.org/10.1016/j.fishres.2004.08.002.
Quinn, Terrance J., and Richard B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press.
Schnute, Jon, and David Fournier. 1980. “A New Approach to Estimating Populations and Stock Composition in Mark-Recapture Experiments.” Canadian Journal of Fisheries and Aquatic Sciences 37 (10): 1744–59. https://doi.org/10.1139/f80-220.
Webber, Darcy N., and Jeremy R. McKenzie. 2022. Rock Lobster (CRA 2) Stock Assessment for the 2021-22 Fishing Year. 2022/XX. New Zealand Fisheries Stock Assessment Report.
Webber, Darcy N., and Jeremy R. McKenzie. 2024. Rock Lobster (CRA 2) Stock Assessment for the 2023-24 Fishing Year. 2024/XX. New Zealand Fisheries Assessment Report.