# remotes::install_github("noaa-afsc/SparseNUTS")
if (dir.exists(file.path("..", "R")) && requireNamespace("pkgload", quietly = TRUE)) {
pkgload::load_all("..", quiet = TRUE)
} else {
library(decamod)
}
library(TMBhelper)
library(tidyverse)
library(kableExtra)
library(SparseNUTS)
theme_set(theme_bw())
options(mc.cores = parallel::detectCores() - 1) # Specify the number of coresThe surplus production model (SPM) is a state-space implementation of the Schaefer surplus production model (Schaefer 1954). The model is described in detail in the decamod package documentation. Here, we provide an example of how to fit the SPM to catch and CPUE data using TMB and the adnuts package for Bayesian inference.
where is the observed catch and is the surplus production defined as:
where is the carrying capacity, is the intrinsic growth rate, and is a shape parameter. When , the model reduces to the Schaefer model (Schaefer 1954). More general formulations, including the Pella-Tomlinson model, allow for asymmetric surplus production (Hilborn and Walters 1992; Quinn and Deriso 1999).
Load inputs
Define model
Define the data list.
catch <- read_csv(system.file("extdata", "CRA2_catch_data_input.csv", package = "decamod")) %>%
# catch <- read_csv("CRA2_catch_data_input.csv") %>%
# catch <- catch_CRA2 %>%
filter(Year > 1978) %>%
mutate(TotalCatch = Commercial + Recreational + Illegal + Customary) %>%
select(Year, Season, Catch = TotalCatch) %>%
arrange(Year, Season) %>%
mutate(Period = 1:n())
ys <- catch %>%
# mutate(Year = as.character(Year)) %>%
select(Year, Season, Period)
logbook_CRA2 <- read_csv(system.file("extdata", "CRA2_LB_CPUE_legal_2025.csv", package = "decamod"))
cpue <- bind_rows(
fsu_CRA2 %>% mutate(Index = "FSU"),
celr_CRA2 %>% mutate(Index = "CELR"),
logbook_CRA2 %>% mutate(Index = "Logbook")
) %>%
filter(Year > 1978) %>%
separate(Year, into = c("Year", "season"), sep = "_") %>%
mutate(Year = as.numeric(Year), Season = ifelse(season == "AW", 1, 2)) %>%
# mutate(Year = as.character(Year)) %>%
mutate(q = case_when(
Index == "FSU" & season == "AW" ~ 1,
Index == "FSU" & season == "SS" ~ 2,
Index == "CELR" & season == "AW" ~ 3,
Index == "CELR" & season == "SS" ~ 4,
Index == "Logbook" & season == "AW" ~ 5,
Index == "Logbook" & season == "SS" ~ 6
)) %>%
mutate(logSD = log(1 + SD / Mean)) %>%
left_join(ys, by = join_by(Year, Season)) %>%
select(Year, Season, Period, Median, logSD, q, Index)
data <- list(
n_time = nrow(catch),
catch_time = catch$Period,
catch_season = catch$Season,
catch_obs = catch$Catch,
cpue_time = cpue$Period,
cpue_season = cpue$Season,
cpue_q = cpue$q,
cpue_obs = cpue$Median,
cpue_sd = cpue$logSD
)
n_q <- length(unique(cpue$q))
n_time <- max(data$catch_time)Define the parameters list.
parameters <- list(
log_r = log(0.35),
log_z = log(1),
log_K = log(5000),
logit_D0 = qlogis(0.5), # logit ensures 0 < D0 < 1, qlogis accepts values up to 0.999999
log_q = log(rep(0.001, n_q)),
log_cpue_pow = log(1),
log_cpue_pro = log(rep(1e-6, n_q)),
cpue_creep = 0.01,
log_sigmap = log(0.001),
log_B = log(rep(2000, n_time - 1))
# epsilon = rep(0, n_time)
)
priors <- list()
priors[["log_r"]] <- list(type = "normal", par1 = log(0.35), par2 = 0.5, index = which("log_r" == names(parameters)))
priors[["log_z"]] <- list(type = "normal", par1 = 0, par2 = 0.1, index = which("log_z" == names(parameters)))
priors[["log_K"]] <- list(type = "normal", par1 = log(1000), par2 = 1, index = which("log_K" == names(parameters)))
priors[["logit_D0"]] <- list(type = "normal", par1 = qlogis(0.5), par2 = 0.5, index = which("logit_D0" == names(parameters)))
priors[["log_q"]] <- list(type = "normal", par1 = -5, par2 = 2, index = which("log_q" == names(parameters)))
priors[["log_cpue_pow"]] <- list(type = "normal", par1 = 0, par2 = 0.5, index = which("log_cpue_pow" == names(parameters)))
priors[["log_cpue_pro"]] <- list(type = "normal", par1 = 0, par2 = 0.5, index = which("log_cpue_pro" == names(parameters)))
priors[["log_sigmap"]] <- list(type = "normal", par1 = 0, par2 = 0.5, index = which("log_sigmap" == names(parameters)))
evaluate_priors(parameters, priors)
#> [1] -2401.345
# plot_priors(priors, bounds)
data$priors <- priorsCreate the AD object using the MakeADFun function from the TMB framework (Kristensen et al. 2016).
# obj <- MakeADFun(func = cmb(spm_model, data), parameters = parameters, map = map, random = "epsilon")
# obj <- MakeADFun(func = cmb(spm_model2, data), parameters = parameters, map = map, random = "log_B")
obj <- MakeADFun(func = cmb(spm_model, data), parameters = parameters, map = map, random = "log_B")
obj$simulate(par = obj$par)$cpue_obs
#> [1] 2.687577 3.027571 3.541143 3.491146 3.515212 3.543480 3.479302 4.215388
#> [9] 3.629795 4.250220 5.105758 3.899659 3.642400 4.429896 3.793459 4.491205
#> [17] 4.074774 5.132360 4.846928 5.514685 4.863845 3.280657 4.817914 4.637801
#> [25] 5.563388 4.547952 4.092199 5.134533 5.262241 4.821345 4.526684 5.286203
#> [33] 6.511813 5.854187 6.127017 3.908947 5.545008 5.765723 4.386044 4.126689
#> [41] 5.667079 4.812655 4.234969 5.124201 4.233965 5.123837 3.912046 4.837231
#> [49] 4.813287 4.707006 4.910167 4.730031 5.915419 4.278650 5.112183 7.274353
#> [57] 5.539699 4.376253 6.030540 5.828605 5.412883 7.934622 5.185701 5.619452
#> [65] 5.824710 6.490928 5.725824 5.267796 5.420655 6.488855 5.781846 5.209975
#> [73] 5.141072 6.208150 6.937789 5.901593 5.916523 6.278203 4.992709 6.614596
#> [81] 4.650385 5.890326 5.848492 4.550509 6.602852 5.780458 6.041071 4.964126
#> [89] 4.675593 5.286262 5.198266 5.890184 5.622123 5.570120 4.906593 4.792782
#> [97] 4.261422 5.062471 4.732164 3.983256 5.412246 5.781185 4.503419 5.047918
#> [105] 6.733603 4.715853 6.029188 5.481357 5.627525 4.160644 6.374323 4.866754
#> [113] 5.040196 4.995101 4.507915 4.752334 7.121900 4.939067 6.212892 3.807689
#> [121] 4.371059 5.101213 5.340937 5.150036 6.051449 5.111527 8.032227 5.881890
#> [129] 6.468808 7.614346 5.625994 5.525706 4.257260 5.429393 4.512276 6.133499
#> [137] 4.919313 8.362317 4.758794 5.480985 7.958114 7.846789 6.043545 4.891042List of parameters that are on:
Define parameter bounds.
Lwr <- rep(-Inf, length(obj$par))
Upr <- rep(Inf, length(obj$par))
Lwr[grep("log_r", names(obj$par))] <- -10
Upr[grep("log_r", names(obj$par))] <- 10
Lwr[grep("log_z", names(obj$par))] <- log(0)
Upr[grep("log_z", names(obj$par))] <- log(2)
Lwr[grep("log_K", names(obj$par))] <- log(300)
Upr[grep("log_K", names(obj$par))] <- log(300000)
Lwr[grep("logit_D0", names(obj$par))] <- qlogis(0.3)
Upr[grep("logit_D0", names(obj$par))] <- qlogis(0.7)
Lwr[grep("log_cpue_pow", names(obj$par))] <- log(0)
Upr[grep("log_cpue_pow", names(obj$par))] <- log(2)
Lwr[grep("log_q", names(obj$par))] <- -25
Upr[grep("log_q", names(obj$par))] <- 1
Lwr[grep("log_sigmap", names(obj$par))] <- log(0)
Upr[grep("log_sigmap", names(obj$par))] <- log(0.5)
Lwr[grep("log_cpue_pro", names(obj$par))] <- log(0)
Upr[grep("log_cpue_pro", names(obj$par))] <- log(2)
bounds <- data.frame(lower = Lwr, upper = Upr)
tibble(Parameter = names(obj$par), LB = Lwr, Initial = obj$par, UB = Upr)
#> # A tibble: 11 × 4
#> Parameter LB Initial UB
#> <chr> <dbl> <dbl> <dbl>
#> 1 log_r -10 -1.05 10
#> 2 log_z -Inf 0 0.693
#> 3 log_K 5.70 8.52 12.6
#> 4 logit_D0 -0.847 0 0.847
#> 5 log_q -25 -6.91 1
#> 6 log_q -25 -6.91 1
#> 7 log_q -25 -6.91 1
#> 8 log_q -25 -6.91 1
#> 9 log_q -25 -6.91 1
#> 10 log_q -25 -6.91 1
#> 11 log_sigmap -Inf -6.91 -0.693Optimise using the nlminb function:
control <- list(eval.max = 100000, iter.max = 100000)
opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr,
control = control, lower = bounds$lower, upper = bounds$upper)
opt <- nlminb(start = opt$par, objective = obj$fn, gradient = obj$gr,
control = control, lower = bounds$lower, upper = bounds$upper)
obj$par <- opt$par
obj$env$last.par <- opt$par
obj$env$last.par.best <- opt$par
report <- obj$report(obj$env$last.par.best)
# check_estimability(obj = obj)
tibble(par = names(obj$env$last.par.best), val = exp(obj$env$last.par.best))
#> # A tibble: 11 × 2
#> par val
#> <chr> <dbl>
#> 1 log_r 0.261
#> 2 log_z 0.988
#> 3 log_K 2449.
#> 4 logit_D0 1.73
#> 5 log_q 0.000725
#> 6 log_q 0.000853
#> 7 log_q 0.000839
#> 8 log_q 0.00114
#> 9 log_q 0.000827
#> 10 log_q 0.00109
#> 11 log_sigmap 0.126
# tibble(par = names(obj$par), val = obj$env$last.par.best, gr = obj$gr()) %>%
# arrange(-gr)
# check_estimability(obj = obj)
biomass_df <- data.frame(Year = catch$Year, Season = catch$Season,
Pred = report$B,
Expected = report$Bmu) %>%
mutate(Season = ifelse(Season == 1, "AW", "SS"))
set.seed(123)
sim_df <- NULL
for (i in 1:100) {
df <- data.frame(Year = catch$Year, Season = catch$Season,
Sim = c(NA, exp(obj$simulate(par = obj$env$last.par.best)$log_B)),
iter = i)
sim_df <- bind_rows(sim_df, df)
}
sim_df <- sim_df %>%
mutate(Season = ifelse(Season == 1, "AW", "SS"))
ggplot(data = biomass_df, aes(x = as.numeric(Year), y = Pred)) +
geom_point(aes(y = Expected, color = "Expected")) +
geom_line(data = sim_df, aes(y = Sim, color = "Simulated", group = iter), alpha = 0.5) +
geom_line(aes(color = "Predicted")) +
facet_grid(Season ~ .) +
scale_x_continuous(breaks = seq(1980, 2050, by = 5)) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
labs(x = "Fishing year", y = "Biomass (tonnes)", color = NULL)
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
load(system.file("extdata", "base_biomass_AW.rda", package = "decamod"))
biomass_AW <- biomass_df %>% filter(Season == "AW")
biomass_AW$Base <- base_biomass_AW
ggplot(data = biomass_AW, aes(x = as.numeric(Year), y = Pred)) +
geom_point(aes(y = Expected, color = "Expected")) +
geom_line(aes(color = "Predicted")) +
geom_line(aes(y = Base, color = "Base")) +
scale_x_continuous(breaks = seq(1980, 2050, by = 5)) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
labs(x = "Fishing year", y = "AW Biomass (tonnes)", color = NULL)
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
index_labels <- c("FSU", "FSU", "CELR", "CELR", "Logbook", "Logbook")
cpue_df <- data.frame(Period = data$cpue_time, Season = data$cpue_season,
Index = index_labels[data$cpue_q],
Obs = data$cpue_obs,
Pred = report$cpue_pred) %>%
left_join(ys, by = join_by(Period, Season)) %>%
mutate(Index = factor(Index, levels = c("FSU", "CELR", "Logbook"))) %>%
mutate(Season = ifelse(Season == 1, "AW", "SS"))
set.seed(456)
sim_df <- NULL
for (i in 1:100) {
df <- data.frame(Period = data$cpue_time, Season = data$cpue_season,
Index = index_labels[data$cpue_q],
Sim = obj$simulate(par = obj$env$last.par.best)$cpue_obs,
iter = i)
sim_df <- bind_rows(sim_df, df)
}
sim_df <- sim_df %>%
left_join(ys, by = join_by(Period, Season)) %>%
mutate(Index = factor(Index, levels = c("FSU", "CELR", "Logbook"))) %>%
mutate(Season = ifelse(Season == 1, "AW", "SS"))
ggplot(data = cpue_df, aes(x = as.numeric(Year), y = Obs)) +
geom_point(aes(color = "Observed")) +
geom_line(data = sim_df, aes(y = Sim, color = "Simulated", group = iter), alpha = 0.5) +
geom_line(aes(y = Pred, group = interaction(Index, Season), color = "Predicted"), linewidth = 1.2) +
facet_grid(Index ~ Season, scales = "free_y") +
scale_x_continuous(breaks = seq(1980, 2050, by = 5)) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
labs(x = "Fishing year", y = "CPUE", color = NULL)
#> `geom_line()`: Each group consists of only one observation.
#> ℹ Do you need to adjust the group aesthetic?
tibble(par = names(opt$par), value = exp(opt$par))
#> # A tibble: 11 × 2
#> par value
#> <chr> <dbl>
#> 1 log_r 0.261
#> 2 log_z 0.988
#> 3 log_K 2449.
#> 4 logit_D0 1.73
#> 5 log_q 0.000725
#> 6 log_q 0.000853
#> 7 log_q 0.000839
#> 8 log_q 0.00114
#> 9 log_q 0.000827
#> 10 log_q 0.00109
#> 11 log_sigmap 0.126
# obj$simulate(par = obj$env$last.par.best)$cpue_obs
# osa <- oneStepPredict(obj = obj, method = "oneStepGeneric", trace = FALSE)
# osa <- oneStepPredict(obj = obj, observation.name = "cpue_obs", method = "oneStepGeneric", trace = FALSE)
# qqnorm(osa$res); abline(0, 1)
# chk <- checkConsistency(obj = obj, hessian = TRUE, estimate = TRUE, n = 200)
# chk
# s <- summary(chk)
# s
# s$marginal$p.valueBayesian inference
Using Cole Monnahan’s adnuts package, we can run the MCMC sampler. The adnuts package is a wrapper around the NUTS sampler in Stan. The objective function is evaluated using automatic differentiation (Kristensen et al. 2016).
# mcmc <- sample_sparse_tmb(
# obj = obj, metric = "auto", iter = 4000, warmup = 3000, chains = 4, cores = 4,
# control = list(adapt_delta = 0.99), init = "last.par.best", globals = decamod_globals())
# mcmc1 <- sample_snuts(
# obj = obj, metric = "auto", num_samples = 1000, chains = 4, cores = 4,
# control = list(adapt_delta = 0.99), init = "last.par.best",
# # lower = Lwr, upper = Upr,
# globals = list(get_growth_increment = get_growth_increment, drobust = drobust,
# jac_exp = jac_exp, jac_invlogit = jac_invlogit))
# library(priorsense)
# lp <- get_posterior(object = obj, posterior = mcmc1, pars = c("loglikes", "priors"), type = "array")
# samps <- extract_samples(fit = mcmc1, inc_lp = FALSE, as.list = FALSE)
# samps$lprior <- as.numeric(lp$priors[,,1])
# samps$log_lik <- as.numeric(lp$loglikes[,,1])
# psd <- create_priorsense_data(x = samps)
# powerscale_sensitivity(psd)
# powerscale_plot_dens(psd, variable = "log_galpha")
# powerscale_plot_ecdf(psd)
# powerscale_plot_quantities(psd)
#
# x <- example_powerscale_model()
# psd <- create_priorsense_data(x$draws)
# powerscale_sensitivity(psd)
# powerscale_plot_dens(psd)
# powerscale_plot_ecdf(psd)
# powerscale_plot_quantities(psd)
# plot_sampler_params(fit = mcmc1, plot = TRUE)
# pairs_rtmb(fit = mcmc1, order = "mismatch", pars = 1:5)
# pairs_rtmb(fit = mcmc1, order = "slow", pars = 1:5)
# plot_marginals(fit = mcmc1)
# plot_uncertainties(fit = mcmc1, log = TRUE, plot = TRUE)Likelihood profiles
The likelihood-profile helper needs to be expanded so the user can list what to report, as the current version does not work for the growth model.
# lincomb <- numeric(length(names(obj$par)))
# lc <- lincomb; lc[1] <- 1
# # prof_galpha1 <- tmbprofile2(obj = obj, lincomb = lc)
# prof_galpha1 <- TMB::tmbprofile(obj = obj, lincomb = lc)
# lc <- lincomb; lc[2] <- 1
# prof_galpha2 <- TMB::tmbprofile(obj = obj, lincomb = lc)
# plot(prof_galpha1)
# plot(prof_galpha2)
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.
Quinn, Terrance J., and Richard B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press.
Schaefer, Milner B. 1954. “Some Aspects of the Dynamics of Populations Important to the Management of the Commercial Marine Fisheries.” Bulletin of the Inter-American Tropical Tuna Commission 1 (2): 27–56.