Nothing
R/RCM_rmd.R
In SAMtool: Stock Assessment Methods Toolkit
Defines functions
make_unique_names
rmd_RCM_likelihood_gradients
RCM_get_likelihoods
plot_composition_RCM
rmd_RCM_Hist_compare
rmd_RCM_retrospective
rmd_RCM_SR
rmd_RCM_SPR2
rmd_RCM_SPR
rmd_log_rec_dev
rmd_RCM_SSB_output
rmd_RCM_R_output
rmd_RCM_initD
rmd_RCM_index_output
rmd_RCM_fleet_output
rmd_RCM_sel
rmd_RCM_Find
rmd_RCM_Perr
rmd_RCM_D
rmd_RCM_R0
rmd_fit_comps
rmd_assess_resid2
rmd_assess_fit2
rmd_matplot
rmd_persp_plot
rmd_persp_plot <- function(x, y, z, xlab, ylab, zlab, phi, theta, expand, fig.cap, header = NULL) {
ans <- c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
paste0("persp(x = ", x, ", y = ", y, ", z = ", z, ", theta = ", theta, ", phi = ", phi, ", expand = ", expand, ", xlab = \"", xlab, "\",
ylab = \"", ylab, "\", zlab = \"", zlab, "\", ticktype = \"detailed\")"),
"```\n")
if (!is.null(header)) ans <- c(header, ans)
return(ans)
}
rmd_matplot <- function(x, y, col, xlab, ylab, legend.lab, type = "l", lty = 1, fig.cap, header = NULL) {
ans <- c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
paste0("yy <- ", y),
paste0("matplot(", x, ", ", y, ", type = \"", type, "\", lty = ", lty, ", col = ", col,
", ylim = c(0, 1.1 * max(yy, na.rm = TRUE)), xlab = \"", xlab, "\", ylab = \"", ylab, "\")"),
"abline(h = 0, col = \"grey\")",
paste0("if (ncol(yy) > 1) legend(\"topleft\", ", legend.lab, ", text.col = ", col, ")"),
"```\n")
if (!is.null(header)) ans <- c(header, ans)
return(ans)
}
# For RCM function
rmd_assess_fit2 <- function(year, obs, fit, fig.cap, label = fig.cap, match = FALSE) {
fig.cap2 <- paste0("Observed (black) and predicted (red) ", fig.cap, ".")
if (match) fig.cap2 <- paste(fig.cap2, "Predicted", fig.cap, "should match observed in this model.")
c(paste0("```{r, fig.cap = \"", fig.cap2, "\"}"),
paste0("plot_timeseries(", year, ", ", obs, ", ", fit, ", label = \"", label, "\")"),
"```\n")
}
# For RCM function
rmd_assess_resid2 <- function(year, obs, fit, fig.cap, label = fig.cap) {
fig.cap2 <- paste0("Index residuals (in log space) for ", fig.cap, ".")
c(paste0("```{r, fig.cap = \"", fig.cap2, "\"}"),
paste0("if (!all(is.na(", obs, "))) {"),
paste0(" istart <- which(!is.na(", obs, "))[1]"),
paste0(" istop <- which(!is.na(", obs, ")) %>% max()"),
paste0(" plot_residuals(", year, "[istart:istop], log(", obs, "[istart:istop]/", fit, "[istart:istop]), label = \"", label, "\")"),
"}",
"```\n")
}
rmd_fit_comps <- function(year, obs, fit, type = c("bubble_data", "annual", "bubble_residuals", "mean", "heat_residuals"),
ages = "NULL", CAL_bins = "NULL", N = "NULL", fig.cap,
bubble_adj = "1.5") {
type <- match.arg(type)
arg <- paste0("\"", type, "\", CAL_bins = ", CAL_bins, ", ages = ", ages)
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
paste0("ind_valid <- rowSums(", obs, ", na.rm = TRUE) > 0"),
paste0("if (any(ind_valid)) plot_composition(", year, "[ind_valid], ", obs, "[ind_valid, , drop = FALSE], ", fit, "[ind_valid, , drop = FALSE], plot_type = ", arg, ", N = ", N, "[ind_valid], bubble_adj = ", bubble_adj, ")"),
"```\n")
}
rmd_RCM_R0 <- function(fig.cap = "Histogram of R0 (unfished recruitment).") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"if (!is.null(OM@cpars$R0)) hist(OM@cpars$R0, main = \"\", xlab = expression(R[0]))",
"```\n")
}
rmd_RCM_D <- function(fig.cap = "Histogram of historical depletion.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"if (!is.null(OM@cpars$D)) hist(OM@cpars$D, main = \"\", xlab = \"Depletion\")",
"```\n")
}
rmd_RCM_Perr <- function(fig.cap = "Historical recruitment deviations among simulations.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"Perr <- OM@cpars$Perr_y[, (max_age+1):(max_age+nyears), drop = FALSE]",
"matplot(Year, t(Perr), type = \"l\", col = \"black\", xlab = \"Year\", ylab = \"Recruitment deviations\",",
" ylim = c(0, 1.1 * max(Perr)))",
"abline(h = 0, col = \"grey\")",
"```\n",
"",
"```{r, fig.cap = \"Future recruitment deviations (up to 5 simulations).\"}",
"Perr_future <- OM@cpars$Perr_y[, (max_age+nyears+1):(max_age+nyears+proyears)]",
"matplot(Year, t(Perr), type = \"l\", col = \"black\", xlab = \"Year\", ylab = \"Recruitment deviations\",",
" xlim = c(min(Year), max(Year) + proyears), ylim = c(0, 1.1 * max(c(Perr, Perr_future))))",
"matlines(max(Year) + 1:proyears, t(Perr_future[1:min(5, nrow(OM@cpars$Perr_y)), ]), type = \"l\")",
"abline(h = 0, col = \"grey\")",
"abline(v = max(Year), lty = 3)",
"```\n",
"",
"```{r, fig.cap = \"Annual mean and median of future recruitment deviations.\"}",
"matplot(Year, t(Perr), type = \"n\", xlab = \"Year\", ylab = \"Recruitment deviations\",",
" xlim = c(min(Year), max(Year) + proyears), ylim = c(0, 1.1 * max(c(Perr, Perr_future))))",
"abline(h = c(0, 1), col = \"grey\")",
"abline(v = max(Year), lty = 3)",
"matlines(Year, t(Perr), type = \"l\", col = \"black\")",
"lines(max(Year) + 1:proyears, apply(Perr_future, 2, mean), col = \"red\")",
"lines(max(Year) + 1:proyears, apply(Perr_future, 2, median), lty = 2)",
"legend(\"topleft\", c(\"Mean\", \"Median\"), col = c(\"red\", \"black\"), lty = c(1, 2))",
"```\n",
"```{r, fig.cap = \"Histogram of recruitment autocorrelation.\"}",
"if (!is.null(OM@cpars$AC)) hist(OM@cpars$AC, main = \"\", xlab = \"Recruitment Autocorrelation\")",
"```\n")
}
rmd_RCM_Find <- function(fig.cap = "Apical F from RCM model.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"matplot(Year, t(OM@cpars$Find), type = \"l\", col = \"black\", xlab = \"Year\", ylab = \"Apical F\")",
"abline(h = 0, col = \"grey\")",
"```\n")
}
rmd_RCM_sel <- function(fig.cap = "Operating model selectivity in the last historical year (all simulations).") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"vul <- sapply(report_list, getElement, \"vul_len\")",
"if (nsim == 1) V_plot <- matrix(OM@cpars$V[, , nyears], 1, byrow = TRUE) else V_plot <- OM@cpars$V[, , nyears]",
"matplot(age, t(V_plot), type = \"l\", col = \"black\", xlab = \"Age\", ylab = \"Selectivity\", ylim = c(0, 1.1))",
"abline(h = 0, col = \"grey\")",
"```\n")
}
rmd_RCM_fleet_output <- function(ff, f_name) {
ans <- c(paste("### ", f_name[ff], "\n"),
paste0("```{r, fig.cap = \"Length selectivity of ", f_name[ff], ".\"}"),
paste0("bl <- unique(RCMdata@sel_block[, ", ff, "])"),
"vul_bb <- list()",
"bl_col <- gplots::rich.colors(length(bl))",
"Year_legend <- character(length(bl))",
"for(bb in 1:length(bl)) {",
" vul_bb[[bb]] <- sapply(report_list, function(x) if (!is.null(x$vul_len)) x$vul_len[, bl[bb]] else NA)",
paste0(" Year_legend[bb] <- Year[RCMdata@sel_block[, ", ff, "] == bl[bb]] %>% range() %>% paste(collapse = \"-\")"),
"}",
"test <- vapply(vul_bb, function(x) all(!is.na(x)) && sum(x), logical(1))",
"if (any(test)) {",
paste0(" matplot(RCMdata@Misc$lbinmid, RCMdata@Misc$lbinmid, type = \"n\", xlab = \"Length\", ylim = c(0, 1), ylab = \"Selectivity of ", f_name[ff], "\")"),
" abline(h = 0, col = \"grey\")",
" for(bb in 1:length(bl)) {",
" if (test[bb]) matlines(RCMdata@Misc$lbinmid, vul_bb[[bb]], type = \"l\", col = bl_col[bb], lty = scenario$lty, lwd = scenario$lwd)",
" }",
" if (length(bl) > 1) legend(\"topright\", Year_legend, col = bl_col, lwd = 1)",
#"if (!is.null(scenario$names)) legend("topleft", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Age-based selectivity of ", f_name[ff], ".\"}"),
paste0("matplot(age, age, type = \"n\", xlab = \"Age\", ylim = c(0, 1), ylab = \"Selectivity of ", f_name[ff], "\")"),
"abline(h = 0, col = \"grey\")",
"",
"for(bb in 1:length(bl)) {",
paste0(" vul_bb_age <- do.call(rbind, lapply(report_list, function(x) x$vul[RCMdata@sel_block[, ", ff, "] == bl[bb], , ", ff, "])) %>% t()"),
" matlines(age, vul_bb_age, type = \"l\", col = bl_col[bb], lty = scenario$lty, lwd = scenario$lwd)",
"}",
"if (length(bl) > 1) legend(\"topleft\", Year_legend, col = bl_col, lwd = 1)",
"```\n",
"",
paste0("```{r, fig.cap = \"Fishing Mortality of ", f_name[ff], ".\"}"),
paste0("FM <- sapply(report_list, function(x) x$F[, ", ff, "])"),
paste0("matplot(Year, FM, type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, ylim = c(0, 1.1 * max(FM, na.rm = TRUE)), xlab = \"Year\", ylab = \"Fishing Mortality of ", f_name[ff], "\")"),
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) catch from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@Chist > 0, na.rm = TRUE)) {"),
paste0(" Cpred <- sapply(report_list, function(x) x$Cpred[, ", ff, "])"),
paste0(" Chist <- RCMdata@Chist[, ", ff, "]"),
" ylim <- c(0.9, 1.1) * range(c(Cpred, Chist), na.rm = TRUE)",
paste0(" plot(Year, Chist, type = \"o\", xlab = \"Year\", ylab = \"Catch of ", f_name[ff], "\", ylim = ylim)"),
paste0(" matlines(Year, Cpred, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)"),
"} else {",
paste0(" Cpred <- sapply(report_list, function(x) x$Cpred[, ", ff, "]/mean(x$Cpred[, ", ff, "]))"),
paste0(" matplot(Year, Cpred, col = scenario$col, type = \"l\", lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Predicted relative catch of ", f_name[ff], "\")"),
"}",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean ages from ", f_name[ff], ".\"}"),
paste0("MApred <- sapply(report_list, function(x) x$CAApred[, , ", ff, "] %*% age/x$CN[, ", ff, "])"),
paste0("MAobs <- (RCMdata@CAA[, , ", ff, "] %*% age)/rowSums(RCMdata@CAA[, , ", ff, "], na.rm = TRUE)"),
"ylim <- c(0.9, 1.1) * range(c(MApred, MAobs), na.rm = TRUE)",
"matplot(Year, MApred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean age\", ylim = ylim)",
paste0("if (any(RCMdata@CAA[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" lines(Year, MAobs, col = \"black\", typ = \"o\")"),
"}",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean lengths from ", f_name[ff], ".\"}"),
paste0("if (RCMdata@MS_type == \"length\" && any(RCMdata@MS[, ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" MLobs <- RCMdata@MS[, ", ff, "]"),
paste0(" MLpred <- sapply(report_list, function(x) x$MLpred[, ", ff, "])"),
paste0("} else if (any(RCMdata@CAL[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" MLobs <- (RCMdata@CAL[, , ", ff, "] %*% RCMdata@Misc$lbinmid)/rowSums(RCMdata@CAL[, , ", ff, "], na.rm = TRUE)"),
paste0(" MLpred <- sapply(report_list, function(x) (x$CALpred[, , ", ff, "] %*% RCMdata@Misc$lbinmid)/rowSums(x$CAL[, , ", ff, "], na.rm = TRUE))"),
"} else {",
" MLobs <- MLpred <- NA",
"}",
"if (!all(is.na(MLpred))) {",
" ylim <- c(0.9, 1.1) * range(c(MLpred, MLobs), na.rm = TRUE)",
" matplot(Year, MLpred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean length\", ylim = ylim)",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean weights from ", f_name[ff], ".\"}"),
paste0("if (RCMdata@MS_type == \"weight\" && any(RCMdata@MS[, ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" MWobs <- RCMdata@MS[, ", ff, "]"),
paste0("} else MWobs <- NA"),
"if (!all(is.na(MWobs))) {",
paste0(" MWpred <- sapply(report_list, function(x) x$MWpred[, ", ff, "])"),
" ylim <- c(0.9, 1.1) * range(c(MWpred, MWobs), na.rm = TRUE)",
" matplot(Year, MWpred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean weight\", ylim = ylim)",
" lines(Year, MWobs, col = \"black\", typ = \"o\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) age composition from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@CAA[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0("if (nsim == 1) CAA_plot <- array(x@CAA[, , , ", ff, "], c(1, nyears, max_age + 1)) else CAA_plot <- x@CAA[, , , ", ff, "]"),
paste0("plot_composition_RCM(Year, CAA_plot, RCMdata@CAA[, , ", ff, "], N = round(RCMdata@CAA_ESS[, ", ff, "], 1), ages = age, dat_col = scenario$col)"),
"}",
"```\n",
paste0("```{r, fig.cap = \"Predicted age composition from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@CAA[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0("plot_composition_RCM(Year, CAA_plot, ages = age, dat_col = scenario$col)"),
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) length composition from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@CAL[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0("if (nsim == 1) CAL_plot <- array(x@CAL[, , , ", ff, "], c(1, nyears, RCMdata@Misc$nlbin)) else CAL_plot <- x@CAL[, , , ", ff, "]"),
paste0("plot_composition_RCM(Year, fit = CAL_plot, dat = RCMdata@CAL[, , ", ff, "], N = round(RCMdata@CAL_ESS[, ", ff, "], 1), CAL_bins = length_bin, dat_col = scenario$col)"),
"}",
"```\n")
return(ans)
}
rmd_RCM_index_output <- function(sur, s_name) {
ans <- c(paste0("### ", s_name[sur], " \n"),
"",
paste0("```{r, fig.cap = \"Length selectivity of ", s_name[sur], ".\"}"),
"ivul_len <- sapply(report_list, function(x) if (!is.null(x$ivul_len)) x$ivul_len[, ", sur, "] else NA)",
"if (all(!is.na(ivul_len)) && sum(ivul_len, na.rm = TRUE)) {",
paste0(" matplot(RCMdata@Misc$lbinmid, ivul_len, type = \"l\", xlab = \"Length\", ylim = c(0, 1), ylab = \"Selectivity of ", s_name[sur], "\")"),
" abline(h = 0, col = \"grey\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Age-based selectivity of ", s_name[sur], " in last historical year.\"}"),
"if (!is.null(report_list[[1]]$ivul)) {",
paste0("ivul_ff_age <- sapply(report_list, function(x) x$ivul[nyears, , ", sur, "])"),
paste0("matplot(age, ivul_ff_age, type = \"l\", col = scenario$col2, xlab = \"Age\", ylim = c(0, 1), ylab = \"Selectivity of ", s_name[sur], "\")"),
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) index values for ", s_name[sur], ".\"}"),
paste0("Ipred <- sapply(report_list, function(x) x$Ipred[, ", sur, "])"),
paste0("matplot(Year, Ipred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, ylim = c(0, 1.1 * max(c(Ipred, RCMdata@Index[, ", sur, "]), na.rm = TRUE)), xlab = \"Year\", ylab = \"", s_name[sur], "\")"),
paste0("lines(Year, RCMdata@Index[, ", sur, "], col = \"black\", typ = \"o\")"),
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) index values for ", s_name[sur], ". Error bars indicate 95% confidence intervals for observed values.\"}"),
paste0("II <- RCMdata@Index[, ", sur, "]"),
"if (any(II > 0, na.rm = TRUE) && length(RCMdata@I_sd)) {",
" ind <- seq(min(which(!is.na(II))), max(which(!is.na(II))), 1)",
paste0(" err <- exp(log(II) + outer(RCMdata@I_sd[, ", sur, "], c(-1.96, 1.96)))"),
paste0(" matplot(Year[ind], Ipred[ind, , drop = FALSE], type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, ylim = c(0, 1.1 * max(c(Ipred[ind, ], II[ind], err[ind, ]), na.rm = TRUE)), xlab = \"Year\", ylab = \"", s_name[sur], "\")"),
" points(Year[ind], II[ind], lwd = 3, pch = 16)",
" arrows(Year[ind], y0 = err[ind, 1], y1 = err[ind, 2], length = 0, lwd = 1.5)",
" abline(h = 0, col = \"grey\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean ages from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAA) && any(RCMdata@IAA[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("MApred <- sapply(report_list, function(x) x$IAApred[, , ", sur, "] %*% age/rowSums(x$IAApred[, , ", sur, "]))"),
paste0("MAobs <- (RCMdata@IAA[, , ", sur, "] %*% age)/rowSums(RCMdata@IAA[, , ", sur, "], na.rm = TRUE)"),
"ylim <- c(0.9, 1.1) * range(c(MApred, MAobs), na.rm = TRUE)",
"matplot(Year, MApred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean age\", ylim = ylim)",
paste0("if (any(RCMdata@IAA[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0(" lines(Year, MAobs, col = \"black\", typ = \"o\")"),
"}",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean lengths from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAL) && any(RCMdata@IAL[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("MLpred <- sapply(report_list, function(x) x$IALpred[, , ", sur, "] %*% RCMdata@Misc$lbinmid/rowSums(x$IAApred[, , ", sur, "]))"),
paste0("MLobs <- (RCMdata@IAL[, , ", sur, "] %*% RCMdata@Misc$lbinmid)/rowSums(RCMdata@IAL[, , ", sur, "], na.rm = TRUE)"),
"ylim <- c(0.9, 1.1) * range(c(MLpred, MLobs), na.rm = TRUE)",
"matplot(Year, MLpred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean length\", ylim = ylim)",
paste0("if (any(RCMdata@IAL[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0(" lines(Year, MLobs, col = \"black\", typ = \"o\")"),
"}",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) age composition from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAA) && any(RCMdata@IAA[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("pred_IAA <- sapply(report_list, function(x) x$IAA[, , ", sur, "], simplify = \"array\") %>% aperm(perm = c(3, 1, 2))"),
paste0("plot_composition_RCM(Year, pred_IAA, RCMdata@IAA[, , ", sur, "], N = round(RCMdata@IAA_ESS[, ", sur, "], 1), ages = age, dat_col = scenario$col)"),
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) length composition from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAL) && any(RCMdata@IAL[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("pred_IAL <- sapply(report_list, function(x) x$IAL[, , ", sur, "], simplify = \"array\") %>% aperm(perm = c(3, 1, 2))"),
paste0("plot_composition_RCM(Year, pred_IAL, RCMdata@IAL[, , ", sur, "], N = round(RCMdata@IAL_ESS[, ", sur, "], 1), CAL_bins = length_bin, dat_col = scenario$col)"),
"}",
"```\n"
)
ans
}
rmd_RCM_initD <- function(fig.cap = "Histogram of initial depletion among all simulations.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"initD <- vapply(report_list, function(x) x$E[1]/x$E0_SR, numeric(1))",
"hist(initD, main = \"\", xlab = \"Initial depletion\")",
"```\n")
}
rmd_RCM_R_output <- function() {
c("```{r, fig.cap = \"Estimated recruitment among all simulations.\"}",
"R_out <- sapply(report_list, getElement, \"R\")",
"matplot(Yearplusone, R_out, ylim = c(0, 1.1 * max(R_out, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Recruitment\")",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_RCM_SSB_output <- function() {
c("```{r, fig.cap = \"Estimated spawning biomass among all simulations.\"}",
"E <- sapply(report_list, getElement, \"E\")",
"matplot(Yearplusone, E, ylim = c(0, 1.1 * max(E, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Spawning biomass\")",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Estimated spawning depletion among all simulations. Unfished spawning biomass is the value calculated from first year life history parameters.\"}",
"E_E0 <- sapply(report_list, function(x) x$E/x$E0_SR)",
"matplot(Yearplusone, E_E0, ylim = c(0, 1.1 * max(E_E0, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Spawning depletion\")",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Dynamic SSB0 among all simulations. Model is re-run assuming no historical catches.\"}",
"dyn_SSB0 <- sapply(report_list, function(x) x$dynamic_SSB0)",
"matplot(Yearplusone, dyn_SSB0, ylim = c(0, 1.1 * max(dyn_SSB0, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = expression(Dynamic~~SSB[0]))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_log_rec_dev <- function() {
c("```{r, fig.cap = \"Estimated recruitment deviations among all simulations.\"}",
"log_rec_dev2 <- sapply(report_list, getElement, \"log_rec_dev\")",
"matplot(Year, log_rec_dev2, type = \"n\", xlab = \"Year\", ylab = \"log-recruitment deviations\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, log_rec_dev2, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_RCM_SPR <- function() {
c("```{r, fig.cap = \"Equilibrium spawning potential ratio (SPR) calculated from the biological parameters and F-at-age in the corresponding year for all simulations.\"}",
"SPR_eq <- sapply(report_list, getElement, \"SPR_eq\")",
"matplot(Year, SPR_eq, type = \"n\", ylim = c(0, 1), xlab = \"Year\", ylab = \"Equilibrium SPR\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, SPR_eq, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Dynamic spawning potential ratio (SPR) calculated from the biological parameters and cumulative survival of the cohorts in the corresponding year for all simulations.\"}",
"SPR_dyn <- sapply(report_list, getElement, \"SPR_dyn\")",
"matplot(Year, SPR_dyn, type = \"n\", ylim = c(0, 1), xlab = \"Year\", ylab = \"Dynamic SPR\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, SPR_dyn, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_RCM_SPR2 <- function() {
c("```{r, fig.cap = \"Annual spawning potential ratio (SPR). Equilibrium SPR is calculated from the F-at-age in the corresponding year, while dynamic SPR is calculated from the cumulative survival of cohorts.\"}",
"plot(Year, report$SPR_eq, typ = \"o\", lty = 3, ylim = c(0, 1), xlab = \"Year\", ylab = \"Spawning potential ratio\")",
"abline(h = 0, col = \"grey\")",
"lines(Year, report$SPR_dyn, typ = \"o\", pch = 16)",
"legend(\"bottomleft\", c(\"Equilibrium SPR\", \"Dynamic SPR\"), lty = c(3, 1), pch = c(1, 16))",
"```\n")
}
rmd_RCM_SR <- function() {
c("```{r, fig.cap = \"Stock-recruit relationship and estimated recruitment. The red point indicates unfished values.\"}",
"if (!is.null(report$Rec_dev)) {",
" expectedR <- report$R/c(report$Rec_dev, 1)",
"} else if (OM@SRrel %in% c(1, 2)) {",
" expectedR <- R_pred(report$E, report[[\"h\"]], report[[\"R0\"]], report$E0_SR, switch(OM@SRrel, \"1\" = \"BH\", \"2\" = \"Ricker\"))",
"} else stop(\"Error in plotting recruitment\")",
"E_full <- seq(0, 1.1 * max(report$E, report$E0_SR), length.out = 100)",
"if (OM@SRrel %in% c(1, 2)) {",
" expectedR_full <- R_pred(E_full, report[[\"h\"]], report[[\"R0\"]], report$E0_SR, switch(OM@SRrel, \"1\" = \"BH\", \"2\" = \"Ricker\"))",
"} else {",
" expectedR_full <- R_pred(E_full, SR_type = \"Mesnil-Rochet\", Shinge = report$MRhinge, Rmax = report$MRRmax, gamma = report$MRgamma)",
"}",
"plot_SR(report$E, expectedR, report$R0, report$E0_SR, report$R)",
"lines(E_full, expectedR_full, lty = 2)",
"```\n",
"",
"```{r, fig.cap = \"Stock-recruit relationship with trajectory.\"}",
"plot_SR(report$E, expectedR, report$R0, report$E0_SR, report$R %>% structure(names = Yearplusone), trajectory = TRUE)",
"lines(E_full, expectedR_full, lty = 2)",
"```\n")
}
rmd_RCM_retrospective <- function(render_args) {
if (render_args$output_format == "html_document") {
x <- "summary(retro) %>% as.data.frame()"
} else {
x <- "summary(retro) %>% as.data.frame() %>% knitr::kable(format = \"markdown\")"
}
c("### Retrospective\n",
"```{r}",
x,
"plot(retro)",
"```\n")
}
rmd_RCM_Hist_compare <- function() {
c("### RCM-OM comparison {.tabset}\n",
"#### OM historical period\n\n",
"```{r, fig.cap = \"Apical F from the operating model.\"}",
"Hist_F <- apply(Hist@AtAge$F.Mortality, c(1, 3), max, na.rm = TRUE)",
"matplot(Year, t(Hist_F), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM Apical F\", ylim = c(0, 1.1 * max(Hist_F)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Spawning biomass (SSB) from the operating model.\"}",
"Hist_SSB <- apply(Hist@TSdata$SBiomass, 1:2, sum)",
"matplot(Year, t(Hist_SSB), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM SSB\", ylim = c(0, 1.1 * max(Hist_SSB)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Spawning biomass (SSB) relative to MSY from the operating model.\"}",
"matplot(Year, t(Hist_SSB/Hist@Ref$ReferencePoints$SSBMSY), col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, typ = \"l\", xlab = \"Year\", ylab = expression(OM~~SSB/SSB[MSY]), ylim = c(0, 1.1 * max(Hist_SSB/Hist@Ref$ReferencePoints$SSBMSY)))",
"abline(h = c(0, MSY_ref), col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Spawning biomass (SSB) relative to MSY in the most recent decade.\"}",
"if (length(Year) > 10) {",
" Yr_ind <- Year > max(Year) - 10",
" matplot(Year[Yr_ind], t(Hist_SSB[, Yr_ind, drop = FALSE]/Hist@Ref$ReferencePoints$SSBMSY), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = expression(OM~~SSB/SSB[MSY]), ylim = c(0, 1.1 * max(Hist_SSB[, Yr_ind, drop = FALSE]/Hist@Ref$ReferencePoints$SSBMSY)))",
" abline(h = c(0, MSY_ref), col = \"grey\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
"```{r, fig.cap = \"Spawning depletion from the operating model.\"}",
"matplot(Year, t(Hist_SSB/Hist@Ref$ReferencePoints$SSB0), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = expression(OM~~SSB/SSB[0]), ylim = c(0, 1.1 * max(Hist_SSB/Hist@Ref$ReferencePoints$SSB0)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Recruitment (age-0) from the operating model.\"}",
"Hist_R <- Hist@AtAge$Number[, 1, , ] %>% apply(c(1, 2), sum)",
"matplot(Year, t(Hist_R), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM Recruitment\", ylim = c(0, 1.1 * max(Hist_R)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Catch (total removals, including discards) from the operating model.\"}",
"Hist_C <- apply(Hist@TSdata$Removals, 1:2, sum)",
"matplot(Year, t(Hist_C), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM Catch\", ylim = c(0, 1.1 * max(Hist_C)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"#### Comparison\n\n",
"```{r, fig.cap = \"Apical F comparison between the OM and RCM.\"}",
"matplot(Year, t(Hist_F), typ = \"o\", pch = 16, col = \"red\", xlab = \"Year\", ylab = \"Apical F\", ylim = c(0, 1.1 * max(c(Hist_F, OM@cpars$Find))))",
"matlines(Year, t(OM@cpars$Find), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in apical F between the OM and RCM. Positive values indicate higher F in the OM.\"}",
"matplot(Year, t(Hist_F - OM@cpars$Find), typ = \"n\", xlab = \"Year\", ylab = \"Difference in apical F\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, t(Hist_F - OM@cpars$Find), col = \"black\")",
"```\n",
"",
"```{r, fig.cap = \"SSB comparison between the OM and RCM.\"}",
"matplot(Year, t(Hist_SSB), typ = \"o\", col = \"red\", pch = 16, xlab = \"Year\", ylab = \"SSB\",",
" ylim = c(0, 1.1 * max(c(Hist_SSB, x@SSB[sims, 1:OM@nyears]))))",
"matlines(Year, t(x@SSB[sims, 1:OM@nyears, drop = FALSE]), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in spawning biomass (SSB), relative to SSB0, between the OM and RCM, calculated as $(SSB^{OM}_y - SSB^{RCM}_y)/SSB^{OM}_0$. Positive values indicate higher SSB in the OM.\"}",
"matplot(Year, t((Hist_SSB - x@SSB[sims, 1:OM@nyears, drop = FALSE])/Hist@Ref$ReferencePoints$SSB0), typ = \"n\", xlab = \"Year\", ylab = \"Difference in relative SSB\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, t((Hist_SSB - x@SSB[sims, 1:OM@nyears, drop = FALSE])/Hist@Ref$ReferencePoints$SSB0), col = \"black\")",
"```\n",
"",
"```{r, fig.cap = \"Recruitment comparison between the OM and RCM.\"}",
"matplot(Year, t(Hist_R), typ = \"o\", col = \"red\", pch = 16, xlab = \"Year\", ylab = \"Recruitment\",",
" ylim = c(0, 1.1 * max(c(Hist_R, x@NAA[sims, 1:OM@nyears, 1]))))",
"matlines(Year, t(x@NAA[, 1:OM@nyears, 1][sims, , drop = FALSE]), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in recruitment (relative to R0) between the OM and RCM, calculated as $(R^{OM}_y - R^{RCM}_y)/R^{OM}_0$. Positive values indicate higher recruitment in the OM.\"}",
"matplot(Year, t(Hist_R/OM@cpars$R0 - x@NAA[, 1:OM@nyears, 1][sims, , drop = FALSE]/OM@cpars$R0), typ = \"n\", xlab = \"Year\", ylab = \"Difference in relative recruitment\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, t(Hist_R/OM@cpars$R0 - x@NAA[, 1:OM@nyears, 1][sims, , drop = FALSE]/OM@cpars$R0),",
" col = \"black\")",
"```\n",
"",
"```{r, fig.cap = \"Comparison of total removals between the OM and RCM.\"}",
"RCM_C <- sapply(x@Misc, function(xx) rowSums(xx$Cpred)) %>% t()",
"if (nrow(RCM_C) == 1) RCM_C <- replicate(x@OM@nsim, as.numeric(RCM_C)) %>% t()",
"matplot(Year, t(Hist_C), typ = \"o\", col = \"red\", pch = 16, xlab = \"Year\", ylab = \"Total removals\",",
" ylim = c(0, 1.1 * max(Hist_C, RCMdata@Chist)))",
"matlines(Year, t(RCM_C[sims, , drop = FALSE]), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in removals (relative to RCM predicted), calculated as $C^{OM}_y/C^{RCM}_y - 1$. Positive values indicate higher removals in the OM.\"}",
"Catch_difference <- t(Hist_C/RCM_C[sims, , drop = FALSE]) - 1",
"Catch_difference[is.infinite(Catch_difference)] <- 0",
"matplot(Year, Catch_difference, typ = \"n\", xlab = \"Year\", ylab = \"Difference in relative removals\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, Catch_difference, col = \"black\")",
"```\n")
}
plot_composition_RCM <- function(Year, fit, dat = NULL, CAL_bins = NULL, ages = NULL, annual_ylab = "Frequency",
annual_yscale = c("proportions", "raw"),
N = if (is.null(dat)) NULL else round(rowSums(dat)),
dat_linewidth = 2, dat_color = "black") {
old_par <- par(no.readonly = TRUE)
on.exit(par(old_par))
par(mfcol = c(4, 4), mar = rep(0, 4), oma = c(5.1, 5.1, 2.1, 2.1))
annual_yscale <- match.arg(annual_yscale)
if (is.null(CAL_bins)) data_type <- "age" else data_type <- "length"
if (data_type == 'length') {
data_val <- CAL_bins
data_lab <- "Length"
}
if (data_type == 'age') {
data_val <- if (is.null(ages)) 1:dim(fit)[3] else ages
data_lab <- "Age"
}
# Annual comps (fit vs. dat if available)
if (annual_yscale == "proportions") {
fit_plot <- apply(fit, 1:2, function(x) x/sum(x)) %>% aperm(c(2, 3, 1))
if (!is.null(dat)) {
dat_plot <- dat/rowSums(dat)
} else {
dat_plot <- dat
}
} else {
fit_plot <- fit
dat_plot <- dat
}
ylim <- c(0, 1.1 * max(fit_plot, dat_plot, na.rm = TRUE))
yaxp <- c(0, max(pretty(ylim, n = 4)), 4)
if (max(c(fit_plot, dat_plot), na.rm = TRUE) == 1) yaxp <- c(0, 1, 4)
las <- 1
for(i in 1:length(Year)) {
yaxt <- ifelse(i %% 16 %in% c(1:4), "s", "n") # TRUE = first column
xaxt <- ifelse(i < length(Year) & i %% 4 %in% c(1:3), "n", "s") # TRUE = first three rows
if (dim(fit_plot)[1] == 1) {
plot(data_val, fit_plot[, i, ], type = "n", ylim = ylim, yaxp = yaxp, xaxt = xaxt, yaxt = yaxt, las = las)
abline(h = 0, col = 'grey')
lines(data_val, fit_plot[, i, ], col = dat_color)
} else {
matplot(data_val, t(fit_plot[, i, ]), type = "n", ylim = ylim, yaxp = yaxp, xaxt = xaxt, yaxt = yaxt, las = las)
abline(h = 0, col = 'grey')
matlines(data_val, t(fit_plot[, i, ]), col = dat_color)
}
abline(h = 0, col = 'grey')
if (!is.null(dat)) lines(data_val, dat_plot[i, ], lwd = 1.5)
legend("topright", legend = c(Year[i], ifelse(is.null(N) | is.na(N[i]), "", paste("N =", N[i]))), bty = "n", xjust = 1)
if (i %% 16 == 1) {
mtext(data_lab, side = 1, line = 3, outer = TRUE)
mtext(annual_ylab, side = 2, line = 3.5, outer = TRUE)
}
}
invisible()
}
RCM_get_likelihoods <- function(x, LWT, f_name, s_name) {
if (inherits(x$nll_fleet, "array")) {
nll_fleet <- apply(x$nll_fleet, 2:3, sum) %>% t()
} else {
nll_fleet <- x$nll_fleet %>% t()
}
nll_fleet[is.na(nll_fleet)] <- 0
if (missing(f_name)) f_name <- paste("Fleet", 1:ncol(nll_fleet))
nll_fleet <- cbind(nll_fleet, rowSums(nll_fleet))
nll_fleet <- rbind(nll_fleet, colSums(nll_fleet))
colnames(nll_fleet) <- c(f_name, "Sum")
rownames(nll_fleet) <- c("Catch", "Equilibrium Catch", "CAA", "CAL", "Mean Size", "Sum")
lwt_fleet <- rbind(LWT$Chist, LWT$C_eq, LWT$CAA, LWT$CAL, LWT$MS) %>%
structure(dimnames = list(rownames(nll_fleet)[1:5], f_name))
if (inherits(x$nll_index, "array")) {
nll_index <- apply(x$nll_index, 2:3, sum) %>% t()
} else {
nll_index <- x$nll_index %>% t()
}
nll_index[is.na(nll_index)] <- 0
if (missing(s_name)) s_name <- paste("Index", 1:ncol(nll_index))
nll_index <- cbind(nll_index, rowSums(nll_index))
nll_index <- rbind(nll_index, colSums(nll_index))
colnames(nll_index) <- c(s_name, "Sum")
rownames(nll_index) <- c("Index", "IAA", "IAL", "Sum")
lwt_index <- rbind(LWT$Index, LWT$IAA, LWT$IAL) %>%
structure(dimnames = list(rownames(nll_index)[1:3], s_name))
tot <- c(
x$nll,
x$nll_log_rec_dev,
nll_fleet["Sum", "Sum"],
nll_index["Sum", "Sum"],
x$penalty,
x$prior
) %>%
matrix(ncol = 1) %>%
structure(dimnames = list(c("Total", "Recruitment Deviations", "Fleets", "Indices", "Penalty (High F)", "Priors"),
"Negative log-likelihood"))
res <- list(tot = tot,
nll_fleet = nll_fleet, lwt_fleet = lwt_fleet,
nll_index = nll_index, lwt_index = lwt_index) %>%
lapply(FUN = function(y) round(y, 2) %>% as.data.frame())
return(res)
}
rmd_RCM_likelihood_gradients <- function(f_name, s_name, do_index) {
header <- c("```{r}",
"obj <- x@mean_fit$obj",
"new_dat <- structure(obj$env$data, check.passed = NULL)",
"new_dat$nll_gr <- 1L",
"new_par <- as.list(SD, \"Estimate\") %>% structure(what = NULL)",
"",
"obj2 <- MakeADFun(data = new_dat, parameters = new_par, map = obj$env$map, random = obj$env$random, ADreport = TRUE,",
" DLL = obj$env$DLL, silent = obj$env$silent)",
"gr <- obj2$gr() %>% structure(dimnames = list(rownames = names(obj2$fn()), colnames = names(SD$par.fixed)))",
"par_names <- data.frame(par_type = colnames(gr), par = make_unique_names(names(SD$par.fixed)))",
"```\n\n")
if (!requireNamespace("ggplot2", quietly = TRUE) || !requireNamespace("reshape2", quietly = TRUE)) {
body <- c("#### NULL \n\nInstall the ggplot2 and reshape2 packages to plot gradients.\n\n")
} else {
fleet_lapply_fn <- function(ff) {
c(paste("####", f_name[ff], "\n"),
"```{r, fig.cap = \"Likelihood gradients (annual values by data type in columns) with respect to model parameters (rows).\"}",
paste0("gr_plot <- dplyr::filter(gr_fleet, Fleet == \"", f_name[ff], "\")"),
"if (nrow(gr_plot)) {",
" ggplot2::ggplot(gr_plot, ggplot2::aes(Year, Gradient, group = par, colour = par)) + ggplot2::facet_grid(par_type ~ data_type, scales = \"free_y\") +",
paste0(" ggplot2::geom_hline(yintercept = 0, linetype = 3) + ggplot2::geom_line() + ggplot2::theme_bw() + ggplot2::theme(legend.position = \"none\") + ggplot2::ggtitle(\"", f_name[ff], "\")"),
"}",
"```\n\n")
}
index_lapply_fn <- function(sur) {
c(paste("####", s_name[sur], "\n"),
"```{r, fig.cap = \"Likelihood gradients (annual values by data type in columns) with respect to model parameters (rows).\"}",
paste0("gr_plot <- dplyr::filter(gr_index, Index == \"", s_name[sur], "\")"),
"if (nrow(gr_plot)) {",
" ggplot2::ggplot(gr_plot, ggplot2::aes(Year, Gradient, group = par, colour = par)) + ggplot2::facet_grid(par_type ~ data_type, scales = \"free_y\") +",
paste0(" ggplot2::geom_hline(yintercept = 0, linetype = 3) + ggplot2::geom_line() + ggplot2::theme_bw() + ggplot2::theme(legend.position = \"none\") + ggplot2::ggtitle(\"", s_name[sur], "\")"),
"}",
"```\n\n")
}
f_plots <- lapply(1:length(f_name), fleet_lapply_fn) %>% unlist()
if (do_index) {
s_plots <- lapply(1:length(s_name), index_lapply_fn) %>% unlist()
} else {
s_plots <- NULL
}
body <- c("```{r}",
"gr_fleet <- gr[rownames(gr) == \"nll_fleet\", ] %>% array(dim = dim(report$nll_fleet) %>% c(length(SD$par.fixed))) %>%",
" structure(dimnames = list(Year = Year, Fleet = f_name, data_type = c(\"Catch\", \"Equilibrium Catch\", \"CAA\", \"CAL\", \"Mean size\"),",
" par = par_names$par)) %>%",
" reshape2::melt(value.name = \"Gradient\") %>% dplyr::left_join(par_names, by = \"par\") %>%",
" dplyr::group_by(data_type) %>% dplyr::filter(any(Gradient != 0))",
"",
"gr_index <- gr[rownames(gr) == \"nll_index\", ] %>% array(dim = dim(report$nll_index) %>% c(length(SD$par.fixed))) %>%",
" structure(dimnames = list(Year = Year, Index = s_name, data_type = c(\"Index\", \"CAA\", \"CAL\"),",
" par = par_names$par)) %>%",
" reshape2::melt(value.name = \"Gradient\") %>% dplyr::left_join(par_names, by = \"par\") %>%",
" dplyr::group_by(data_type) %>% dplyr::filter(any(Gradient != 0))",
"```\n\n", f_plots, s_plots)
}
if (requireNamespace("caret", quietly = TRUE)) {
jac <- c("#### Linear combos\n",
"",
"```{r}",
"gr_combo <- gr[rownames(gr) %in% c(\"nll_fleet\", \"nll_index\"), ] %>% caret::findLinearCombos()",
"if (is.null(gr_combo$remove)) {",
" print(\"Jacobian matrix is of full rank, according to caret::findLinearCombos().\")",
"} else {",
" combo_report <- data.frame(`Parameter Number` = gr_combo$remove, `Parameter Name` = par_names$par[gr_combo$remove])",
" print(\"Jacobian matrix is not of full rank, according to caret::findLinearCombos(). Reduce number of model parameters?\")",
" print(\"See table below:\")",
"}",
"```\n\n",
"```{r} \nif (!is.null(gr_combo$remove)) combo_report\n```\n\n",
"Output of caret::findLinearCombos():\n",
"```{r}\n gr_combo\n```\n\n")
} else {
jac <- c("#### Linear combos\n",
"",
"Install the caret package to evaluate if Jacobian matrix is of full rank.\n\n")
}
c(header, body, jac)
}
make_unique_names <- function(par_names) {
unique_names <- par_names %>% unique()
par_new <- lapply(unique_names, function(x) {
ind <- par_names == x
if (sum(ind) > 1) {
paste0(x, "_", 1:sum(ind))
} else x
})
do.call(c, par_new)
}
Try the SAMtool package in your browser
Any scripts or data that you put into this service are public.
SAMtool documentation built on Nov. 18, 2023, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions make_unique_names rmd_RCM_likelihood_gradients RCM_get_likelihoods plot_composition_RCM rmd_RCM_Hist_compare rmd_RCM_retrospective rmd_RCM_SR rmd_RCM_SPR2 rmd_RCM_SPR rmd_log_rec_dev rmd_RCM_SSB_output rmd_RCM_R_output rmd_RCM_initD rmd_RCM_index_output rmd_RCM_fleet_output rmd_RCM_sel rmd_RCM_Find rmd_RCM_Perr rmd_RCM_D rmd_RCM_R0 rmd_fit_comps rmd_assess_resid2 rmd_assess_fit2 rmd_matplot rmd_persp_plot
rmd_persp_plot <- function(x, y, z, xlab, ylab, zlab, phi, theta, expand, fig.cap, header = NULL) {
ans <- c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
paste0("persp(x = ", x, ", y = ", y, ", z = ", z, ", theta = ", theta, ", phi = ", phi, ", expand = ", expand, ", xlab = \"", xlab, "\",
ylab = \"", ylab, "\", zlab = \"", zlab, "\", ticktype = \"detailed\")"),
"```\n")
if (!is.null(header)) ans <- c(header, ans)
return(ans)
}
rmd_matplot <- function(x, y, col, xlab, ylab, legend.lab, type = "l", lty = 1, fig.cap, header = NULL) {
ans <- c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
paste0("yy <- ", y),
paste0("matplot(", x, ", ", y, ", type = \"", type, "\", lty = ", lty, ", col = ", col,
", ylim = c(0, 1.1 * max(yy, na.rm = TRUE)), xlab = \"", xlab, "\", ylab = \"", ylab, "\")"),
"abline(h = 0, col = \"grey\")",
paste0("if (ncol(yy) > 1) legend(\"topleft\", ", legend.lab, ", text.col = ", col, ")"),
"```\n")
if (!is.null(header)) ans <- c(header, ans)
return(ans)
}
# For RCM function
rmd_assess_fit2 <- function(year, obs, fit, fig.cap, label = fig.cap, match = FALSE) {
fig.cap2 <- paste0("Observed (black) and predicted (red) ", fig.cap, ".")
if (match) fig.cap2 <- paste(fig.cap2, "Predicted", fig.cap, "should match observed in this model.")
c(paste0("```{r, fig.cap = \"", fig.cap2, "\"}"),
paste0("plot_timeseries(", year, ", ", obs, ", ", fit, ", label = \"", label, "\")"),
"```\n")
}
# For RCM function
rmd_assess_resid2 <- function(year, obs, fit, fig.cap, label = fig.cap) {
fig.cap2 <- paste0("Index residuals (in log space) for ", fig.cap, ".")
c(paste0("```{r, fig.cap = \"", fig.cap2, "\"}"),
paste0("if (!all(is.na(", obs, "))) {"),
paste0(" istart <- which(!is.na(", obs, "))[1]"),
paste0(" istop <- which(!is.na(", obs, ")) %>% max()"),
paste0(" plot_residuals(", year, "[istart:istop], log(", obs, "[istart:istop]/", fit, "[istart:istop]), label = \"", label, "\")"),
"}",
"```\n")
}
rmd_fit_comps <- function(year, obs, fit, type = c("bubble_data", "annual", "bubble_residuals", "mean", "heat_residuals"),
ages = "NULL", CAL_bins = "NULL", N = "NULL", fig.cap,
bubble_adj = "1.5") {
type <- match.arg(type)
arg <- paste0("\"", type, "\", CAL_bins = ", CAL_bins, ", ages = ", ages)
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
paste0("ind_valid <- rowSums(", obs, ", na.rm = TRUE) > 0"),
paste0("if (any(ind_valid)) plot_composition(", year, "[ind_valid], ", obs, "[ind_valid, , drop = FALSE], ", fit, "[ind_valid, , drop = FALSE], plot_type = ", arg, ", N = ", N, "[ind_valid], bubble_adj = ", bubble_adj, ")"),
"```\n")
}
rmd_RCM_R0 <- function(fig.cap = "Histogram of R0 (unfished recruitment).") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"if (!is.null(OM@cpars$R0)) hist(OM@cpars$R0, main = \"\", xlab = expression(R[0]))",
"```\n")
}
rmd_RCM_D <- function(fig.cap = "Histogram of historical depletion.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"if (!is.null(OM@cpars$D)) hist(OM@cpars$D, main = \"\", xlab = \"Depletion\")",
"```\n")
}
rmd_RCM_Perr <- function(fig.cap = "Historical recruitment deviations among simulations.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"Perr <- OM@cpars$Perr_y[, (max_age+1):(max_age+nyears), drop = FALSE]",
"matplot(Year, t(Perr), type = \"l\", col = \"black\", xlab = \"Year\", ylab = \"Recruitment deviations\",",
" ylim = c(0, 1.1 * max(Perr)))",
"abline(h = 0, col = \"grey\")",
"```\n",
"",
"```{r, fig.cap = \"Future recruitment deviations (up to 5 simulations).\"}",
"Perr_future <- OM@cpars$Perr_y[, (max_age+nyears+1):(max_age+nyears+proyears)]",
"matplot(Year, t(Perr), type = \"l\", col = \"black\", xlab = \"Year\", ylab = \"Recruitment deviations\",",
" xlim = c(min(Year), max(Year) + proyears), ylim = c(0, 1.1 * max(c(Perr, Perr_future))))",
"matlines(max(Year) + 1:proyears, t(Perr_future[1:min(5, nrow(OM@cpars$Perr_y)), ]), type = \"l\")",
"abline(h = 0, col = \"grey\")",
"abline(v = max(Year), lty = 3)",
"```\n",
"",
"```{r, fig.cap = \"Annual mean and median of future recruitment deviations.\"}",
"matplot(Year, t(Perr), type = \"n\", xlab = \"Year\", ylab = \"Recruitment deviations\",",
" xlim = c(min(Year), max(Year) + proyears), ylim = c(0, 1.1 * max(c(Perr, Perr_future))))",
"abline(h = c(0, 1), col = \"grey\")",
"abline(v = max(Year), lty = 3)",
"matlines(Year, t(Perr), type = \"l\", col = \"black\")",
"lines(max(Year) + 1:proyears, apply(Perr_future, 2, mean), col = \"red\")",
"lines(max(Year) + 1:proyears, apply(Perr_future, 2, median), lty = 2)",
"legend(\"topleft\", c(\"Mean\", \"Median\"), col = c(\"red\", \"black\"), lty = c(1, 2))",
"```\n",
"```{r, fig.cap = \"Histogram of recruitment autocorrelation.\"}",
"if (!is.null(OM@cpars$AC)) hist(OM@cpars$AC, main = \"\", xlab = \"Recruitment Autocorrelation\")",
"```\n")
}
rmd_RCM_Find <- function(fig.cap = "Apical F from RCM model.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"matplot(Year, t(OM@cpars$Find), type = \"l\", col = \"black\", xlab = \"Year\", ylab = \"Apical F\")",
"abline(h = 0, col = \"grey\")",
"```\n")
}
rmd_RCM_sel <- function(fig.cap = "Operating model selectivity in the last historical year (all simulations).") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"vul <- sapply(report_list, getElement, \"vul_len\")",
"if (nsim == 1) V_plot <- matrix(OM@cpars$V[, , nyears], 1, byrow = TRUE) else V_plot <- OM@cpars$V[, , nyears]",
"matplot(age, t(V_plot), type = \"l\", col = \"black\", xlab = \"Age\", ylab = \"Selectivity\", ylim = c(0, 1.1))",
"abline(h = 0, col = \"grey\")",
"```\n")
}
rmd_RCM_fleet_output <- function(ff, f_name) {
ans <- c(paste("### ", f_name[ff], "\n"),
paste0("```{r, fig.cap = \"Length selectivity of ", f_name[ff], ".\"}"),
paste0("bl <- unique(RCMdata@sel_block[, ", ff, "])"),
"vul_bb <- list()",
"bl_col <- gplots::rich.colors(length(bl))",
"Year_legend <- character(length(bl))",
"for(bb in 1:length(bl)) {",
" vul_bb[[bb]] <- sapply(report_list, function(x) if (!is.null(x$vul_len)) x$vul_len[, bl[bb]] else NA)",
paste0(" Year_legend[bb] <- Year[RCMdata@sel_block[, ", ff, "] == bl[bb]] %>% range() %>% paste(collapse = \"-\")"),
"}",
"test <- vapply(vul_bb, function(x) all(!is.na(x)) && sum(x), logical(1))",
"if (any(test)) {",
paste0(" matplot(RCMdata@Misc$lbinmid, RCMdata@Misc$lbinmid, type = \"n\", xlab = \"Length\", ylim = c(0, 1), ylab = \"Selectivity of ", f_name[ff], "\")"),
" abline(h = 0, col = \"grey\")",
" for(bb in 1:length(bl)) {",
" if (test[bb]) matlines(RCMdata@Misc$lbinmid, vul_bb[[bb]], type = \"l\", col = bl_col[bb], lty = scenario$lty, lwd = scenario$lwd)",
" }",
" if (length(bl) > 1) legend(\"topright\", Year_legend, col = bl_col, lwd = 1)",
#"if (!is.null(scenario$names)) legend("topleft", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Age-based selectivity of ", f_name[ff], ".\"}"),
paste0("matplot(age, age, type = \"n\", xlab = \"Age\", ylim = c(0, 1), ylab = \"Selectivity of ", f_name[ff], "\")"),
"abline(h = 0, col = \"grey\")",
"",
"for(bb in 1:length(bl)) {",
paste0(" vul_bb_age <- do.call(rbind, lapply(report_list, function(x) x$vul[RCMdata@sel_block[, ", ff, "] == bl[bb], , ", ff, "])) %>% t()"),
" matlines(age, vul_bb_age, type = \"l\", col = bl_col[bb], lty = scenario$lty, lwd = scenario$lwd)",
"}",
"if (length(bl) > 1) legend(\"topleft\", Year_legend, col = bl_col, lwd = 1)",
"```\n",
"",
paste0("```{r, fig.cap = \"Fishing Mortality of ", f_name[ff], ".\"}"),
paste0("FM <- sapply(report_list, function(x) x$F[, ", ff, "])"),
paste0("matplot(Year, FM, type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, ylim = c(0, 1.1 * max(FM, na.rm = TRUE)), xlab = \"Year\", ylab = \"Fishing Mortality of ", f_name[ff], "\")"),
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) catch from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@Chist > 0, na.rm = TRUE)) {"),
paste0(" Cpred <- sapply(report_list, function(x) x$Cpred[, ", ff, "])"),
paste0(" Chist <- RCMdata@Chist[, ", ff, "]"),
" ylim <- c(0.9, 1.1) * range(c(Cpred, Chist), na.rm = TRUE)",
paste0(" plot(Year, Chist, type = \"o\", xlab = \"Year\", ylab = \"Catch of ", f_name[ff], "\", ylim = ylim)"),
paste0(" matlines(Year, Cpred, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)"),
"} else {",
paste0(" Cpred <- sapply(report_list, function(x) x$Cpred[, ", ff, "]/mean(x$Cpred[, ", ff, "]))"),
paste0(" matplot(Year, Cpred, col = scenario$col, type = \"l\", lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Predicted relative catch of ", f_name[ff], "\")"),
"}",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean ages from ", f_name[ff], ".\"}"),
paste0("MApred <- sapply(report_list, function(x) x$CAApred[, , ", ff, "] %*% age/x$CN[, ", ff, "])"),
paste0("MAobs <- (RCMdata@CAA[, , ", ff, "] %*% age)/rowSums(RCMdata@CAA[, , ", ff, "], na.rm = TRUE)"),
"ylim <- c(0.9, 1.1) * range(c(MApred, MAobs), na.rm = TRUE)",
"matplot(Year, MApred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean age\", ylim = ylim)",
paste0("if (any(RCMdata@CAA[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" lines(Year, MAobs, col = \"black\", typ = \"o\")"),
"}",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean lengths from ", f_name[ff], ".\"}"),
paste0("if (RCMdata@MS_type == \"length\" && any(RCMdata@MS[, ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" MLobs <- RCMdata@MS[, ", ff, "]"),
paste0(" MLpred <- sapply(report_list, function(x) x$MLpred[, ", ff, "])"),
paste0("} else if (any(RCMdata@CAL[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" MLobs <- (RCMdata@CAL[, , ", ff, "] %*% RCMdata@Misc$lbinmid)/rowSums(RCMdata@CAL[, , ", ff, "], na.rm = TRUE)"),
paste0(" MLpred <- sapply(report_list, function(x) (x$CALpred[, , ", ff, "] %*% RCMdata@Misc$lbinmid)/rowSums(x$CAL[, , ", ff, "], na.rm = TRUE))"),
"} else {",
" MLobs <- MLpred <- NA",
"}",
"if (!all(is.na(MLpred))) {",
" ylim <- c(0.9, 1.1) * range(c(MLpred, MLobs), na.rm = TRUE)",
" matplot(Year, MLpred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean length\", ylim = ylim)",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean weights from ", f_name[ff], ".\"}"),
paste0("if (RCMdata@MS_type == \"weight\" && any(RCMdata@MS[, ", ff, "] > 0, na.rm = TRUE)) {"),
paste0(" MWobs <- RCMdata@MS[, ", ff, "]"),
paste0("} else MWobs <- NA"),
"if (!all(is.na(MWobs))) {",
paste0(" MWpred <- sapply(report_list, function(x) x$MWpred[, ", ff, "])"),
" ylim <- c(0.9, 1.1) * range(c(MWpred, MWobs), na.rm = TRUE)",
" matplot(Year, MWpred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean weight\", ylim = ylim)",
" lines(Year, MWobs, col = \"black\", typ = \"o\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) age composition from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@CAA[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0("if (nsim == 1) CAA_plot <- array(x@CAA[, , , ", ff, "], c(1, nyears, max_age + 1)) else CAA_plot <- x@CAA[, , , ", ff, "]"),
paste0("plot_composition_RCM(Year, CAA_plot, RCMdata@CAA[, , ", ff, "], N = round(RCMdata@CAA_ESS[, ", ff, "], 1), ages = age, dat_col = scenario$col)"),
"}",
"```\n",
paste0("```{r, fig.cap = \"Predicted age composition from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@CAA[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0("plot_composition_RCM(Year, CAA_plot, ages = age, dat_col = scenario$col)"),
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) length composition from ", f_name[ff], ".\"}"),
paste0("if (any(RCMdata@CAL[, , ", ff, "] > 0, na.rm = TRUE)) {"),
paste0("if (nsim == 1) CAL_plot <- array(x@CAL[, , , ", ff, "], c(1, nyears, RCMdata@Misc$nlbin)) else CAL_plot <- x@CAL[, , , ", ff, "]"),
paste0("plot_composition_RCM(Year, fit = CAL_plot, dat = RCMdata@CAL[, , ", ff, "], N = round(RCMdata@CAL_ESS[, ", ff, "], 1), CAL_bins = length_bin, dat_col = scenario$col)"),
"}",
"```\n")
return(ans)
}
rmd_RCM_index_output <- function(sur, s_name) {
ans <- c(paste0("### ", s_name[sur], " \n"),
"",
paste0("```{r, fig.cap = \"Length selectivity of ", s_name[sur], ".\"}"),
"ivul_len <- sapply(report_list, function(x) if (!is.null(x$ivul_len)) x$ivul_len[, ", sur, "] else NA)",
"if (all(!is.na(ivul_len)) && sum(ivul_len, na.rm = TRUE)) {",
paste0(" matplot(RCMdata@Misc$lbinmid, ivul_len, type = \"l\", xlab = \"Length\", ylim = c(0, 1), ylab = \"Selectivity of ", s_name[sur], "\")"),
" abline(h = 0, col = \"grey\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Age-based selectivity of ", s_name[sur], " in last historical year.\"}"),
"if (!is.null(report_list[[1]]$ivul)) {",
paste0("ivul_ff_age <- sapply(report_list, function(x) x$ivul[nyears, , ", sur, "])"),
paste0("matplot(age, ivul_ff_age, type = \"l\", col = scenario$col2, xlab = \"Age\", ylim = c(0, 1), ylab = \"Selectivity of ", s_name[sur], "\")"),
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) index values for ", s_name[sur], ".\"}"),
paste0("Ipred <- sapply(report_list, function(x) x$Ipred[, ", sur, "])"),
paste0("matplot(Year, Ipred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, ylim = c(0, 1.1 * max(c(Ipred, RCMdata@Index[, ", sur, "]), na.rm = TRUE)), xlab = \"Year\", ylab = \"", s_name[sur], "\")"),
paste0("lines(Year, RCMdata@Index[, ", sur, "], col = \"black\", typ = \"o\")"),
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) index values for ", s_name[sur], ". Error bars indicate 95% confidence intervals for observed values.\"}"),
paste0("II <- RCMdata@Index[, ", sur, "]"),
"if (any(II > 0, na.rm = TRUE) && length(RCMdata@I_sd)) {",
" ind <- seq(min(which(!is.na(II))), max(which(!is.na(II))), 1)",
paste0(" err <- exp(log(II) + outer(RCMdata@I_sd[, ", sur, "], c(-1.96, 1.96)))"),
paste0(" matplot(Year[ind], Ipred[ind, , drop = FALSE], type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, ylim = c(0, 1.1 * max(c(Ipred[ind, ], II[ind], err[ind, ]), na.rm = TRUE)), xlab = \"Year\", ylab = \"", s_name[sur], "\")"),
" points(Year[ind], II[ind], lwd = 3, pch = 16)",
" arrows(Year[ind], y0 = err[ind, 1], y1 = err[ind, 2], length = 0, lwd = 1.5)",
" abline(h = 0, col = \"grey\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean ages from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAA) && any(RCMdata@IAA[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("MApred <- sapply(report_list, function(x) x$IAApred[, , ", sur, "] %*% age/rowSums(x$IAApred[, , ", sur, "]))"),
paste0("MAobs <- (RCMdata@IAA[, , ", sur, "] %*% age)/rowSums(RCMdata@IAA[, , ", sur, "], na.rm = TRUE)"),
"ylim <- c(0.9, 1.1) * range(c(MApred, MAobs), na.rm = TRUE)",
"matplot(Year, MApred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean age\", ylim = ylim)",
paste0("if (any(RCMdata@IAA[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0(" lines(Year, MAobs, col = \"black\", typ = \"o\")"),
"}",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) mean lengths from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAL) && any(RCMdata@IAL[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("MLpred <- sapply(report_list, function(x) x$IALpred[, , ", sur, "] %*% RCMdata@Misc$lbinmid/rowSums(x$IAApred[, , ", sur, "]))"),
paste0("MLobs <- (RCMdata@IAL[, , ", sur, "] %*% RCMdata@Misc$lbinmid)/rowSums(RCMdata@IAL[, , ", sur, "], na.rm = TRUE)"),
"ylim <- c(0.9, 1.1) * range(c(MLpred, MLobs), na.rm = TRUE)",
"matplot(Year, MLpred, type = \"l\", col = scenario$col, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Mean length\", ylim = ylim)",
paste0("if (any(RCMdata@IAL[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0(" lines(Year, MLobs, col = \"black\", typ = \"o\")"),
"}",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) age composition from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAA) && any(RCMdata@IAA[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("pred_IAA <- sapply(report_list, function(x) x$IAA[, , ", sur, "], simplify = \"array\") %>% aperm(perm = c(3, 1, 2))"),
paste0("plot_composition_RCM(Year, pred_IAA, RCMdata@IAA[, , ", sur, "], N = round(RCMdata@IAA_ESS[, ", sur, "], 1), ages = age, dat_col = scenario$col)"),
"}",
"```\n",
"",
paste0("```{r, fig.cap = \"Observed (black) and predicted (red) length composition from ", s_name[sur], ".\"}"),
paste0("if (length(RCMdata@IAL) && any(RCMdata@IAL[, , ", sur, "] > 0, na.rm = TRUE)) {"),
paste0("pred_IAL <- sapply(report_list, function(x) x$IAL[, , ", sur, "], simplify = \"array\") %>% aperm(perm = c(3, 1, 2))"),
paste0("plot_composition_RCM(Year, pred_IAL, RCMdata@IAL[, , ", sur, "], N = round(RCMdata@IAL_ESS[, ", sur, "], 1), CAL_bins = length_bin, dat_col = scenario$col)"),
"}",
"```\n"
)
ans
}
rmd_RCM_initD <- function(fig.cap = "Histogram of initial depletion among all simulations.") {
c(paste0("```{r, fig.cap = \"", fig.cap, "\"}"),
"initD <- vapply(report_list, function(x) x$E[1]/x$E0_SR, numeric(1))",
"hist(initD, main = \"\", xlab = \"Initial depletion\")",
"```\n")
}
rmd_RCM_R_output <- function() {
c("```{r, fig.cap = \"Estimated recruitment among all simulations.\"}",
"R_out <- sapply(report_list, getElement, \"R\")",
"matplot(Yearplusone, R_out, ylim = c(0, 1.1 * max(R_out, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Recruitment\")",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_RCM_SSB_output <- function() {
c("```{r, fig.cap = \"Estimated spawning biomass among all simulations.\"}",
"E <- sapply(report_list, getElement, \"E\")",
"matplot(Yearplusone, E, ylim = c(0, 1.1 * max(E, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Spawning biomass\")",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Estimated spawning depletion among all simulations. Unfished spawning biomass is the value calculated from first year life history parameters.\"}",
"E_E0 <- sapply(report_list, function(x) x$E/x$E0_SR)",
"matplot(Yearplusone, E_E0, ylim = c(0, 1.1 * max(E_E0, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"Spawning depletion\")",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Dynamic SSB0 among all simulations. Model is re-run assuming no historical catches.\"}",
"dyn_SSB0 <- sapply(report_list, function(x) x$dynamic_SSB0)",
"matplot(Yearplusone, dyn_SSB0, ylim = c(0, 1.1 * max(dyn_SSB0, na.rm = TRUE)), type = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = expression(Dynamic~~SSB[0]))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_log_rec_dev <- function() {
c("```{r, fig.cap = \"Estimated recruitment deviations among all simulations.\"}",
"log_rec_dev2 <- sapply(report_list, getElement, \"log_rec_dev\")",
"matplot(Year, log_rec_dev2, type = \"n\", xlab = \"Year\", ylab = \"log-recruitment deviations\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, log_rec_dev2, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_RCM_SPR <- function() {
c("```{r, fig.cap = \"Equilibrium spawning potential ratio (SPR) calculated from the biological parameters and F-at-age in the corresponding year for all simulations.\"}",
"SPR_eq <- sapply(report_list, getElement, \"SPR_eq\")",
"matplot(Year, SPR_eq, type = \"n\", ylim = c(0, 1), xlab = \"Year\", ylab = \"Equilibrium SPR\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, SPR_eq, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Dynamic spawning potential ratio (SPR) calculated from the biological parameters and cumulative survival of the cohorts in the corresponding year for all simulations.\"}",
"SPR_dyn <- sapply(report_list, getElement, \"SPR_dyn\")",
"matplot(Year, SPR_dyn, type = \"n\", ylim = c(0, 1), xlab = \"Year\", ylab = \"Dynamic SPR\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, SPR_dyn, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n")
}
rmd_RCM_SPR2 <- function() {
c("```{r, fig.cap = \"Annual spawning potential ratio (SPR). Equilibrium SPR is calculated from the F-at-age in the corresponding year, while dynamic SPR is calculated from the cumulative survival of cohorts.\"}",
"plot(Year, report$SPR_eq, typ = \"o\", lty = 3, ylim = c(0, 1), xlab = \"Year\", ylab = \"Spawning potential ratio\")",
"abline(h = 0, col = \"grey\")",
"lines(Year, report$SPR_dyn, typ = \"o\", pch = 16)",
"legend(\"bottomleft\", c(\"Equilibrium SPR\", \"Dynamic SPR\"), lty = c(3, 1), pch = c(1, 16))",
"```\n")
}
rmd_RCM_SR <- function() {
c("```{r, fig.cap = \"Stock-recruit relationship and estimated recruitment. The red point indicates unfished values.\"}",
"if (!is.null(report$Rec_dev)) {",
" expectedR <- report$R/c(report$Rec_dev, 1)",
"} else if (OM@SRrel %in% c(1, 2)) {",
" expectedR <- R_pred(report$E, report[[\"h\"]], report[[\"R0\"]], report$E0_SR, switch(OM@SRrel, \"1\" = \"BH\", \"2\" = \"Ricker\"))",
"} else stop(\"Error in plotting recruitment\")",
"E_full <- seq(0, 1.1 * max(report$E, report$E0_SR), length.out = 100)",
"if (OM@SRrel %in% c(1, 2)) {",
" expectedR_full <- R_pred(E_full, report[[\"h\"]], report[[\"R0\"]], report$E0_SR, switch(OM@SRrel, \"1\" = \"BH\", \"2\" = \"Ricker\"))",
"} else {",
" expectedR_full <- R_pred(E_full, SR_type = \"Mesnil-Rochet\", Shinge = report$MRhinge, Rmax = report$MRRmax, gamma = report$MRgamma)",
"}",
"plot_SR(report$E, expectedR, report$R0, report$E0_SR, report$R)",
"lines(E_full, expectedR_full, lty = 2)",
"```\n",
"",
"```{r, fig.cap = \"Stock-recruit relationship with trajectory.\"}",
"plot_SR(report$E, expectedR, report$R0, report$E0_SR, report$R %>% structure(names = Yearplusone), trajectory = TRUE)",
"lines(E_full, expectedR_full, lty = 2)",
"```\n")
}
rmd_RCM_retrospective <- function(render_args) {
if (render_args$output_format == "html_document") {
x <- "summary(retro) %>% as.data.frame()"
} else {
x <- "summary(retro) %>% as.data.frame() %>% knitr::kable(format = \"markdown\")"
}
c("### Retrospective\n",
"```{r}",
x,
"plot(retro)",
"```\n")
}
rmd_RCM_Hist_compare <- function() {
c("### RCM-OM comparison {.tabset}\n",
"#### OM historical period\n\n",
"```{r, fig.cap = \"Apical F from the operating model.\"}",
"Hist_F <- apply(Hist@AtAge$F.Mortality, c(1, 3), max, na.rm = TRUE)",
"matplot(Year, t(Hist_F), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM Apical F\", ylim = c(0, 1.1 * max(Hist_F)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Spawning biomass (SSB) from the operating model.\"}",
"Hist_SSB <- apply(Hist@TSdata$SBiomass, 1:2, sum)",
"matplot(Year, t(Hist_SSB), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM SSB\", ylim = c(0, 1.1 * max(Hist_SSB)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Spawning biomass (SSB) relative to MSY from the operating model.\"}",
"matplot(Year, t(Hist_SSB/Hist@Ref$ReferencePoints$SSBMSY), col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, typ = \"l\", xlab = \"Year\", ylab = expression(OM~~SSB/SSB[MSY]), ylim = c(0, 1.1 * max(Hist_SSB/Hist@Ref$ReferencePoints$SSBMSY)))",
"abline(h = c(0, MSY_ref), col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Spawning biomass (SSB) relative to MSY in the most recent decade.\"}",
"if (length(Year) > 10) {",
" Yr_ind <- Year > max(Year) - 10",
" matplot(Year[Yr_ind], t(Hist_SSB[, Yr_ind, drop = FALSE]/Hist@Ref$ReferencePoints$SSBMSY), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = expression(OM~~SSB/SSB[MSY]), ylim = c(0, 1.1 * max(Hist_SSB[, Yr_ind, drop = FALSE]/Hist@Ref$ReferencePoints$SSBMSY)))",
" abline(h = c(0, MSY_ref), col = \"grey\")",
" if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"}",
"```\n",
"",
"```{r, fig.cap = \"Spawning depletion from the operating model.\"}",
"matplot(Year, t(Hist_SSB/Hist@Ref$ReferencePoints$SSB0), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = expression(OM~~SSB/SSB[0]), ylim = c(0, 1.1 * max(Hist_SSB/Hist@Ref$ReferencePoints$SSB0)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Recruitment (age-0) from the operating model.\"}",
"Hist_R <- Hist@AtAge$Number[, 1, , ] %>% apply(c(1, 2), sum)",
"matplot(Year, t(Hist_R), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM Recruitment\", ylim = c(0, 1.1 * max(Hist_R)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"```{r, fig.cap = \"Catch (total removals, including discards) from the operating model.\"}",
"Hist_C <- apply(Hist@TSdata$Removals, 1:2, sum)",
"matplot(Year, t(Hist_C), typ = \"l\", col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd, xlab = \"Year\", ylab = \"OM Catch\", ylim = c(0, 1.1 * max(Hist_C)))",
"abline(h = 0, col = \"grey\")",
"if (!is.null(scenario$names)) legend(\"topleft\", scenario$names, col = scenario$col2, lty = scenario$lty, lwd = scenario$lwd)",
"```\n",
"",
"#### Comparison\n\n",
"```{r, fig.cap = \"Apical F comparison between the OM and RCM.\"}",
"matplot(Year, t(Hist_F), typ = \"o\", pch = 16, col = \"red\", xlab = \"Year\", ylab = \"Apical F\", ylim = c(0, 1.1 * max(c(Hist_F, OM@cpars$Find))))",
"matlines(Year, t(OM@cpars$Find), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in apical F between the OM and RCM. Positive values indicate higher F in the OM.\"}",
"matplot(Year, t(Hist_F - OM@cpars$Find), typ = \"n\", xlab = \"Year\", ylab = \"Difference in apical F\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, t(Hist_F - OM@cpars$Find), col = \"black\")",
"```\n",
"",
"```{r, fig.cap = \"SSB comparison between the OM and RCM.\"}",
"matplot(Year, t(Hist_SSB), typ = \"o\", col = \"red\", pch = 16, xlab = \"Year\", ylab = \"SSB\",",
" ylim = c(0, 1.1 * max(c(Hist_SSB, x@SSB[sims, 1:OM@nyears]))))",
"matlines(Year, t(x@SSB[sims, 1:OM@nyears, drop = FALSE]), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in spawning biomass (SSB), relative to SSB0, between the OM and RCM, calculated as $(SSB^{OM}_y - SSB^{RCM}_y)/SSB^{OM}_0$. Positive values indicate higher SSB in the OM.\"}",
"matplot(Year, t((Hist_SSB - x@SSB[sims, 1:OM@nyears, drop = FALSE])/Hist@Ref$ReferencePoints$SSB0), typ = \"n\", xlab = \"Year\", ylab = \"Difference in relative SSB\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, t((Hist_SSB - x@SSB[sims, 1:OM@nyears, drop = FALSE])/Hist@Ref$ReferencePoints$SSB0), col = \"black\")",
"```\n",
"",
"```{r, fig.cap = \"Recruitment comparison between the OM and RCM.\"}",
"matplot(Year, t(Hist_R), typ = \"o\", col = \"red\", pch = 16, xlab = \"Year\", ylab = \"Recruitment\",",
" ylim = c(0, 1.1 * max(c(Hist_R, x@NAA[sims, 1:OM@nyears, 1]))))",
"matlines(Year, t(x@NAA[, 1:OM@nyears, 1][sims, , drop = FALSE]), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in recruitment (relative to R0) between the OM and RCM, calculated as $(R^{OM}_y - R^{RCM}_y)/R^{OM}_0$. Positive values indicate higher recruitment in the OM.\"}",
"matplot(Year, t(Hist_R/OM@cpars$R0 - x@NAA[, 1:OM@nyears, 1][sims, , drop = FALSE]/OM@cpars$R0), typ = \"n\", xlab = \"Year\", ylab = \"Difference in relative recruitment\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, t(Hist_R/OM@cpars$R0 - x@NAA[, 1:OM@nyears, 1][sims, , drop = FALSE]/OM@cpars$R0),",
" col = \"black\")",
"```\n",
"",
"```{r, fig.cap = \"Comparison of total removals between the OM and RCM.\"}",
"RCM_C <- sapply(x@Misc, function(xx) rowSums(xx$Cpred)) %>% t()",
"if (nrow(RCM_C) == 1) RCM_C <- replicate(x@OM@nsim, as.numeric(RCM_C)) %>% t()",
"matplot(Year, t(Hist_C), typ = \"o\", col = \"red\", pch = 16, xlab = \"Year\", ylab = \"Total removals\",",
" ylim = c(0, 1.1 * max(Hist_C, RCMdata@Chist)))",
"matlines(Year, t(RCM_C[sims, , drop = FALSE]), col = \"black\")",
"abline(h = 0, col = \"grey\")",
"legend(\"topleft\", c(\"RCM\", \"OM\"), col = c(\"black\", \"red\"), pch = c(NA_integer_, 16), lwd = c(1, 1), bty = \"n\")",
"```\n",
"",
"```{r, fig.cap = \"Difference in removals (relative to RCM predicted), calculated as $C^{OM}_y/C^{RCM}_y - 1$. Positive values indicate higher removals in the OM.\"}",
"Catch_difference <- t(Hist_C/RCM_C[sims, , drop = FALSE]) - 1",
"Catch_difference[is.infinite(Catch_difference)] <- 0",
"matplot(Year, Catch_difference, typ = \"n\", xlab = \"Year\", ylab = \"Difference in relative removals\")",
"abline(h = 0, col = \"grey\")",
"matlines(Year, Catch_difference, col = \"black\")",
"```\n")
}
plot_composition_RCM <- function(Year, fit, dat = NULL, CAL_bins = NULL, ages = NULL, annual_ylab = "Frequency",
annual_yscale = c("proportions", "raw"),
N = if (is.null(dat)) NULL else round(rowSums(dat)),
dat_linewidth = 2, dat_color = "black") {
old_par <- par(no.readonly = TRUE)
on.exit(par(old_par))
par(mfcol = c(4, 4), mar = rep(0, 4), oma = c(5.1, 5.1, 2.1, 2.1))
annual_yscale <- match.arg(annual_yscale)
if (is.null(CAL_bins)) data_type <- "age" else data_type <- "length"
if (data_type == 'length') {
data_val <- CAL_bins
data_lab <- "Length"
}
if (data_type == 'age') {
data_val <- if (is.null(ages)) 1:dim(fit)[3] else ages
data_lab <- "Age"
}
# Annual comps (fit vs. dat if available)
if (annual_yscale == "proportions") {
fit_plot <- apply(fit, 1:2, function(x) x/sum(x)) %>% aperm(c(2, 3, 1))
if (!is.null(dat)) {
dat_plot <- dat/rowSums(dat)
} else {
dat_plot <- dat
}
} else {
fit_plot <- fit
dat_plot <- dat
}
ylim <- c(0, 1.1 * max(fit_plot, dat_plot, na.rm = TRUE))
yaxp <- c(0, max(pretty(ylim, n = 4)), 4)
if (max(c(fit_plot, dat_plot), na.rm = TRUE) == 1) yaxp <- c(0, 1, 4)
las <- 1
for(i in 1:length(Year)) {
yaxt <- ifelse(i %% 16 %in% c(1:4), "s", "n") # TRUE = first column
xaxt <- ifelse(i < length(Year) & i %% 4 %in% c(1:3), "n", "s") # TRUE = first three rows
if (dim(fit_plot)[1] == 1) {
plot(data_val, fit_plot[, i, ], type = "n", ylim = ylim, yaxp = yaxp, xaxt = xaxt, yaxt = yaxt, las = las)
abline(h = 0, col = 'grey')
lines(data_val, fit_plot[, i, ], col = dat_color)
} else {
matplot(data_val, t(fit_plot[, i, ]), type = "n", ylim = ylim, yaxp = yaxp, xaxt = xaxt, yaxt = yaxt, las = las)
abline(h = 0, col = 'grey')
matlines(data_val, t(fit_plot[, i, ]), col = dat_color)
}
abline(h = 0, col = 'grey')
if (!is.null(dat)) lines(data_val, dat_plot[i, ], lwd = 1.5)
legend("topright", legend = c(Year[i], ifelse(is.null(N) | is.na(N[i]), "", paste("N =", N[i]))), bty = "n", xjust = 1)
if (i %% 16 == 1) {
mtext(data_lab, side = 1, line = 3, outer = TRUE)
mtext(annual_ylab, side = 2, line = 3.5, outer = TRUE)
}
}
invisible()
}
RCM_get_likelihoods <- function(x, LWT, f_name, s_name) {
if (inherits(x$nll_fleet, "array")) {
nll_fleet <- apply(x$nll_fleet, 2:3, sum) %>% t()
} else {
nll_fleet <- x$nll_fleet %>% t()
}
nll_fleet[is.na(nll_fleet)] <- 0
if (missing(f_name)) f_name <- paste("Fleet", 1:ncol(nll_fleet))
nll_fleet <- cbind(nll_fleet, rowSums(nll_fleet))
nll_fleet <- rbind(nll_fleet, colSums(nll_fleet))
colnames(nll_fleet) <- c(f_name, "Sum")
rownames(nll_fleet) <- c("Catch", "Equilibrium Catch", "CAA", "CAL", "Mean Size", "Sum")
lwt_fleet <- rbind(LWT$Chist, LWT$C_eq, LWT$CAA, LWT$CAL, LWT$MS) %>%
structure(dimnames = list(rownames(nll_fleet)[1:5], f_name))
if (inherits(x$nll_index, "array")) {
nll_index <- apply(x$nll_index, 2:3, sum) %>% t()
} else {
nll_index <- x$nll_index %>% t()
}
nll_index[is.na(nll_index)] <- 0
if (missing(s_name)) s_name <- paste("Index", 1:ncol(nll_index))
nll_index <- cbind(nll_index, rowSums(nll_index))
nll_index <- rbind(nll_index, colSums(nll_index))
colnames(nll_index) <- c(s_name, "Sum")
rownames(nll_index) <- c("Index", "IAA", "IAL", "Sum")
lwt_index <- rbind(LWT$Index, LWT$IAA, LWT$IAL) %>%
structure(dimnames = list(rownames(nll_index)[1:3], s_name))
tot <- c(
x$nll,
x$nll_log_rec_dev,
nll_fleet["Sum", "Sum"],
nll_index["Sum", "Sum"],
x$penalty,
x$prior
) %>%
matrix(ncol = 1) %>%
structure(dimnames = list(c("Total", "Recruitment Deviations", "Fleets", "Indices", "Penalty (High F)", "Priors"),
"Negative log-likelihood"))
res <- list(tot = tot,
nll_fleet = nll_fleet, lwt_fleet = lwt_fleet,
nll_index = nll_index, lwt_index = lwt_index) %>%
lapply(FUN = function(y) round(y, 2) %>% as.data.frame())
return(res)
}
rmd_RCM_likelihood_gradients <- function(f_name, s_name, do_index) {
header <- c("```{r}",
"obj <- x@mean_fit$obj",
"new_dat <- structure(obj$env$data, check.passed = NULL)",
"new_dat$nll_gr <- 1L",
"new_par <- as.list(SD, \"Estimate\") %>% structure(what = NULL)",
"",
"obj2 <- MakeADFun(data = new_dat, parameters = new_par, map = obj$env$map, random = obj$env$random, ADreport = TRUE,",
" DLL = obj$env$DLL, silent = obj$env$silent)",
"gr <- obj2$gr() %>% structure(dimnames = list(rownames = names(obj2$fn()), colnames = names(SD$par.fixed)))",
"par_names <- data.frame(par_type = colnames(gr), par = make_unique_names(names(SD$par.fixed)))",
"```\n\n")
if (!requireNamespace("ggplot2", quietly = TRUE) || !requireNamespace("reshape2", quietly = TRUE)) {
body <- c("#### NULL \n\nInstall the ggplot2 and reshape2 packages to plot gradients.\n\n")
} else {
fleet_lapply_fn <- function(ff) {
c(paste("####", f_name[ff], "\n"),
"```{r, fig.cap = \"Likelihood gradients (annual values by data type in columns) with respect to model parameters (rows).\"}",
paste0("gr_plot <- dplyr::filter(gr_fleet, Fleet == \"", f_name[ff], "\")"),
"if (nrow(gr_plot)) {",
" ggplot2::ggplot(gr_plot, ggplot2::aes(Year, Gradient, group = par, colour = par)) + ggplot2::facet_grid(par_type ~ data_type, scales = \"free_y\") +",
paste0(" ggplot2::geom_hline(yintercept = 0, linetype = 3) + ggplot2::geom_line() + ggplot2::theme_bw() + ggplot2::theme(legend.position = \"none\") + ggplot2::ggtitle(\"", f_name[ff], "\")"),
"}",
"```\n\n")
}
index_lapply_fn <- function(sur) {
c(paste("####", s_name[sur], "\n"),
"```{r, fig.cap = \"Likelihood gradients (annual values by data type in columns) with respect to model parameters (rows).\"}",
paste0("gr_plot <- dplyr::filter(gr_index, Index == \"", s_name[sur], "\")"),
"if (nrow(gr_plot)) {",
" ggplot2::ggplot(gr_plot, ggplot2::aes(Year, Gradient, group = par, colour = par)) + ggplot2::facet_grid(par_type ~ data_type, scales = \"free_y\") +",
paste0(" ggplot2::geom_hline(yintercept = 0, linetype = 3) + ggplot2::geom_line() + ggplot2::theme_bw() + ggplot2::theme(legend.position = \"none\") + ggplot2::ggtitle(\"", s_name[sur], "\")"),
"}",
"```\n\n")
}
f_plots <- lapply(1:length(f_name), fleet_lapply_fn) %>% unlist()
if (do_index) {
s_plots <- lapply(1:length(s_name), index_lapply_fn) %>% unlist()
} else {
s_plots <- NULL
}
body <- c("```{r}",
"gr_fleet <- gr[rownames(gr) == \"nll_fleet\", ] %>% array(dim = dim(report$nll_fleet) %>% c(length(SD$par.fixed))) %>%",
" structure(dimnames = list(Year = Year, Fleet = f_name, data_type = c(\"Catch\", \"Equilibrium Catch\", \"CAA\", \"CAL\", \"Mean size\"),",
" par = par_names$par)) %>%",
" reshape2::melt(value.name = \"Gradient\") %>% dplyr::left_join(par_names, by = \"par\") %>%",
" dplyr::group_by(data_type) %>% dplyr::filter(any(Gradient != 0))",
"",
"gr_index <- gr[rownames(gr) == \"nll_index\", ] %>% array(dim = dim(report$nll_index) %>% c(length(SD$par.fixed))) %>%",
" structure(dimnames = list(Year = Year, Index = s_name, data_type = c(\"Index\", \"CAA\", \"CAL\"),",
" par = par_names$par)) %>%",
" reshape2::melt(value.name = \"Gradient\") %>% dplyr::left_join(par_names, by = \"par\") %>%",
" dplyr::group_by(data_type) %>% dplyr::filter(any(Gradient != 0))",
"```\n\n", f_plots, s_plots)
}
if (requireNamespace("caret", quietly = TRUE)) {
jac <- c("#### Linear combos\n",
"",
"```{r}",
"gr_combo <- gr[rownames(gr) %in% c(\"nll_fleet\", \"nll_index\"), ] %>% caret::findLinearCombos()",
"if (is.null(gr_combo$remove)) {",
" print(\"Jacobian matrix is of full rank, according to caret::findLinearCombos().\")",
"} else {",
" combo_report <- data.frame(`Parameter Number` = gr_combo$remove, `Parameter Name` = par_names$par[gr_combo$remove])",
" print(\"Jacobian matrix is not of full rank, according to caret::findLinearCombos(). Reduce number of model parameters?\")",
" print(\"See table below:\")",
"}",
"```\n\n",
"```{r} \nif (!is.null(gr_combo$remove)) combo_report\n```\n\n",
"Output of caret::findLinearCombos():\n",
"```{r}\n gr_combo\n```\n\n")
} else {
jac <- c("#### Linear combos\n",
"",
"Install the caret package to evaluate if Jacobian matrix is of full rank.\n\n")
}
c(header, body, jac)
}
make_unique_names <- function(par_names) {
unique_names <- par_names %>% unique()
par_new <- lapply(unique_names, function(x) {
ind <- par_names == x
if (sum(ind) > 1) {
paste0(x, "_", 1:sum(ind))
} else x
})
do.call(c, par_new)
}
Try the SAMtool package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.