In PaulRegular/MSP: Multispecies production model incorporating covariance
library(plotly)
library(viridis)
knitr::opts_chunk$set(echo = FALSE, fig.height = 7, fig.width = 8, cache = FALSE)
if (fit$call$K_groups == ~1) {
K_groups <- NULL
if (length(unique(fit$landings$species)) == 1) {
K_label <- "K"
} else {
K_label <- "All species"
}
} else {
K_groups <- as.formula(fit$call$K_groups)
K_label <- unique(fit$landings[, all.vars(K_groups)])
}
n_pe_covar <- length(all.vars(fit$call$pe_covariates))
n_K_covar <- length(all.vars(fit$call$K_covariates))
if (fit$call$species_cor != "none") {
n_spp_cor <- length(fit$par$logit_rho)
if (n_spp_cor > 1) {
sp_rho <- sp_nm_mat <- matrix(NA, nrow = nlevels(fit$pop$species), ncol = nlevels(fit$pop$species))
rownames(sp_rho) <- colnames(sp_rho) <- levels(fit$pop$species)
sp_rho[lower.tri(sp_rho)] <- inv_logit(fit$par$logit_rho, shift = TRUE)
sp_rho <- t(sp_rho)
sp_rho[lower.tri(sp_rho)] <- inv_logit(fit$par$logit_rho, shift = TRUE)
diag(sp_rho) <- 1
round(sp_rho, 2)
for (i in seq(nrow(sp_rho))) {
for (j in seq(ncol(sp_rho))) {
sp_nm_mat[i, j] <- paste(levels(fit$pop$species)[i], "-", levels(fit$pop$species)[j])
}
}
species_cor_label <- sp_nm_mat[lower.tri(sp_nm_mat)]
} else {
species_cor_label <- "All species"
}
} else {
n_spp_cor <- 0
species_cor_label <- "N/A"
}
n_time_cor <- as.numeric(fit$call$temporal_cor == "ar1")
temporal_cor_label <- ifelse(n_time_cor == 1, "All species", "N/A")
add_labels <- function(p, xlab = "", ylab = "", log_linear_button = TRUE) {
if (log_linear_button) {
## Modified from https://stackoverflow.com/a/60551263/4055438
buttons <- list(list(
active = 0,
type = "buttons",
buttons= list(
list(label = 'linear',
method = 'relayout',
args = list(list(yaxis = list(type = 'linear', title = ylab)))),
list(label = 'log',
method = 'relayout',
args = list(list(yaxis = list(type = 'log', title = ylab)))))))
} else {
buttons <- NULL
}
layout(p,
xaxis = list(title = xlab),
yaxis = list(title = ylab),
updatemenus = buttons)
}
Inputs
Row {.tabset}
Index
fit$index %>%
plot_ly() %>%
add_trace(x = ~year, y = ~index, color = ~survey,
colors = viridis::viridis(100), mode = "markers+lines",
type = "scatter") %>%
add_labels(xlab = "Year", ylab = "Index")
Landings
fit$landings %>%
plot_ly() %>%
add_trace(x = ~year, y = ~landings, color = ~species,
colors = viridis::viridis(100), mode = "markers+lines",
type = "scatter") %>%
add_labels(xlab = "Year", ylab = "Landings")
Covariates
r if (n_pe_covar == 0 && n_K_covar == 0) "N/A; no process error or K covariates supplied"
r if ((n_pe_covar + n_K_covar) > 1) "Plotting of more than one covariate not yet implemented"
if (n_pe_covar == 1) {
covar_name <- all.vars(fit$call$pe_covariates)
covariates <- unique(fit$landings[, c("year", covar_name), drop = FALSE])
}
if (n_K_covar == 1) {
covar_name <- all.vars(fit$call$K_covariates)
covariates <- unique(fit$landings[, c("year", covar_name), drop = FALSE])
}
names(covariates) <- c("year", "covariate")
covariates %>%
plot_ly(x = ~year, y = ~covariate) %>%
add_lines() %>%
layout(xaxis = list(title = "Year"),
yaxis = list(title = covar_name))
Residuals
Row {.tabset}
Residuals ~ predicted value
p <- fit$index %>%
plot_ly(color = ~species, colors = viridis::viridis(100)) %>%
layout(yaxis = list(title = "Standardized residuals"))
p %>% add_markers(x = ~log(pred), y = ~std_res)
~ year
p %>% add_markers(x = ~year, y = ~std_res)
~ survey
p %>% add_markers(x = ~survey, y = ~std_res)
~ species
p %>% add_markers(x = ~species, y = ~std_res)
Parameters
Row {.tabset}
Correlation
sd_rep <- fit$sd_rep
dsd <- sqrt(diag(sd_rep$cov.fixed))
cor_mat <- diag(1 / dsd) %*% sd_rep$cov.fixed %*% diag(1 / dsd)
rownames(cor_mat) <- paste0(seq(nrow(cor_mat)), ":", names(dsd))
colnames(cor_mat) <- paste0(seq(nrow(cor_mat)), ":", names(dsd))
# corrplot::corrplot.mixed(cor_mat, diag = "n", lower = "ellipse", upper = "number")
plot_ly(x = rownames(cor_mat), y = colnames(cor_mat), z = cor_mat) %>%
add_heatmap(colors = c("#B2182B", "white", "#2166AC")) %>%
colorbar(limits = c(-1.001, 1.001)) %>%
layout(yaxis = list(scaleanchor = "x", scaleratio = 1),
xaxis = list(constrain = "domain"))
log(K)
plot_post(prior_mean = fit$par$mean_log_K,
prior_sd = exp(fit$par$log_sd_log_K),
post_mean = fit$par$log_K,
post_sd = fit$se$log_K,
names = K_label,
xlab = "log(K)",
show_prior = fit$call$log_K_option$option == "prior")
log(r)
plot_post(prior_mean = fit$par$mean_log_r,
prior_sd = exp(fit$par$log_sd_log_r),
post_mean = fit$par$log_r,
post_sd = fit$se$log_r,
names = levels(fit$landings$species),
xlab = "log(r)",
show_prior = fit$call$log_r_option$option == "prior")
log(B0)
plot_post(prior_mean = fit$par$mean_log_B0,
prior_sd = exp(fit$par$log_sd_log_B0),
post_mean = fit$par$log_B0,
post_sd = fit$se$log_B0,
names = levels(fit$landings$species),
xlab = "log(B0)",
show_prior = fit$call$log_B0_option$option == "prior")
log(SDB)
plot_post(prior_mean = fit$par$mean_log_sd_B,
prior_sd = exp(fit$par$log_sd_log_sd_B),
post_mean = fit$par$log_sd_B,
post_sd = fit$se$log_sd_B,
names = levels(fit$landings$species),
xlab = "log(SD<sub>B</sub>)",
show_prior = fit$call$log_sd_B_option$option == "prior")
logit(rho)
r if (n_spp_cor == 0) "N/A; correlation of process errors across species was not estimated."
plot_post(prior_mean = fit$par$mean_logit_rho,
prior_sd = exp(fit$par$log_sd_logit_rho),
post_mean = fit$par$logit_rho,
post_sd = fit$se$logit_rho,
names = species_cor_label,
xlab = "logit(rho)", # trans_fun = function(x) inv_logit(x, shift = TRUE),
show_prior = fit$call$logit_rho_option$option == "prior")
logit(phi)
r if (n_time_cor == 0) "N/A; correlation of process errors across time was not estimated."
plot_post(prior_mean = fit$par$mean_logit_phi,
prior_sd = exp(fit$par$log_sd_logit_phi),
post_mean = fit$par$logit_phi,
post_sd = fit$se$logit_phi,
names = "logit(phi)",
xlab = "logit(phi)", # trans_fun = inv_logit,
show_prior = fit$call$logit_phi_option$option == "prior")
log(q)
plot_post(prior_mean = fit$par$mean_log_q,
prior_sd = exp(fit$par$log_sd_log_q),
post_mean = fit$par$log_q,
post_sd = fit$se$log_q,
names = levels(fit$index$survey),
xlab = "log(q)",
show_prior = fit$call$log_q_option$option == "prior")
log(SDI)
plot_post(prior_mean = fit$par$mean_log_sd_I,
prior_sd = exp(fit$par$log_sd_log_sd_I),
post_mean = fit$par$log_sd_I,
post_sd = fit$se$log_sd_I,
names = levels(fit$index$survey),
xlab = "log(SD<sub>I</sub>)",
show_prior = fit$call$log_sd_I_option$option == "prior")
betape
r if (n_pe_covar == 0) "N/A; no process error covariates supplied"
plot_post(prior_mean = fit$par$mean_pe_betas,
prior_sd = exp(fit$par$log_sd_pe_betas),
post_mean = fit$par$pe_betas,
post_sd = fit$se$pe_betas,
names = "pe_betas",
xlab = "pe_betas",
show_prior = fit$call$pe_betas_option$option == "prior")
betaK
r if (n_K_covar == 0) "N/A; no K covariates supplied"
plot_post(prior_mean = fit$par$mean_K_betas,
prior_sd = exp(fit$par$log_sd_K_betas),
post_mean = fit$par$K_betas,
post_sd = fit$se$K_betas,
names = "K_betas",
xlab = "K_betas",
show_prior = fit$call$K_betas_option$option == "prior")
Estimates (plot)
name_par <- c("log_K", "log_r", "log_B0", "log_sd_B",
"logit_rho", "logit_phi", "log_q", "log_sd_I",
"pe_betas", "K_betas")
par_fun <- function(nm, select = "par") {
if (grepl("beta", nm)) return(fit[[select]][[nm]])
if (grepl("log_", nm)) return(exp(fit[[select]][[nm]]))
if (nm == "logit_rho") return(inv_logit(fit[[select]][[nm]], shift = TRUE))
if (nm == "logit_phi") return(inv_logit(fit[[select]][[nm]], shift = FALSE))
}
named_par <- lapply(name_par, par_fun, select = "par")
named_par_cv <- lapply(name_par, par_fun, select = "se")
named_par_lwr <- lapply(name_par, par_fun, select = "par_lwr")
named_par_upr <- lapply(name_par, par_fun, select = "par_upr")
names(named_par) <- gsub("log_|logit_", "", name_par)
names(named_par$q) <- levels(fit$index$survey)
names(named_par$sd_I) <- levels(fit$index$survey)
names(named_par$K) <- K_label
names(named_par$r) <- levels(fit$index$species)
names(named_par$sd_B) <- levels(fit$index$species)
names(named_par$B0) <- levels(fit$index$species)
names(named_par$phi) <- temporal_cor_label
names(named_par$rho) <- species_cor_label
names(named_par$pe_betas) <- colnames(fit$tmb_dat$pe_covariates)
names(named_par$K_betas) <- colnames(fit$tmb_dat$K_covariates)
par_tab <- data.frame(parameter = rep(names(named_par), sapply(named_par, length)),
group = names(unlist(unname(named_par))),
estimate = unlist(unname(named_par)),
CV = unlist(named_par_cv),
lower = unlist(named_par_lwr),
upper = unlist(named_par_upr))
par_tab <- na.omit(par_tab) # will drop par lacking CV estimates (i.e. not estimated)
par_tab %>%
plot_ly(color = ~parameter, colors = viridis::viridis(100),
legendgroup = ~parameter, y = ~paste(parameter, "-", group)) %>%
add_segments(x = ~lower, xend = ~upper, yend = ~paste(parameter, "-", group),
showlegend = FALSE) %>%
add_markers(x = ~estimate) %>%
layout(xaxis = list(title = ""), yaxis = list(title = ""))
Estimates (table)
par_tab$parameter[duplicated(par_tab$parameter)] <- ""
knitr::kable(par_tab, digits = 3, row.names = FALSE)
Population trends
Row {.tabset}
Observed and predicted index
fit$index %>%
plot_ly(x = ~year, color = ~survey, colors = viridis::viridis(100),
legendgroup = ~survey) %>%
add_ribbons(ymin = ~pred_lwr, ymax = ~pred_upr, line = list(width = 0),
alpha = 0.2, showlegend = FALSE) %>%
add_lines(y = ~pred) %>%
add_markers(y = ~index, showlegend = FALSE) %>%
add_labels(xlab = "Year", ylab = "Index")
Process error
p <- fit$pop %>%
plot_ly(color = ~species, colors = viridis::viridis(100))
p %>% add_lines(x = ~year, y = ~pe) %>%
add_labels(xlab = "Year", ylab = "Process error")
r if (n_pe_covar >= 1) "> Includes covariate effects"
Standardized process error
p <- fit$pop %>%
plot_ly(color = ~species, colors = viridis::viridis(100))
p %>% add_lines(x = ~year, y = ~log_std_res_pe) %>%
layout(xaxis = list(title = "Year"),
yaxis = list(title = "Standardized process error"))
r if (n_pe_covar >= 1) "> Standardized process error = residual process error, in log space, not explained by the covariate effects divided by the standard deviation of the process." else "> Standardized process error = process error, in log space, divided by the standard deviation of the process."
Process error correlation
r if (n_spp_cor == 0) "N/A; correlation of process errors across species was not estimated."
r if (n_spp_cor == 1) "N/A; matrix of correlation was not estimated as rho was coupled across species."
plot_ly(x = rownames(sp_rho), y = rownames(sp_rho), z = ~sp_rho,
colors = c("#B2182B", "white", "#2166AC")) %>%
add_heatmap() %>% colorbar(limits = c(-1.001, 1.001), title = "ρ") %>%
layout(yaxis = list(scaleanchor = "x", scaleratio = 1),
xaxis = list(constrain = "domain"))
r if (n_spp_cor > 1 && n_pe_covar >= 1) "> Represents correlation in residual process error not explained by the covariate effects."
Biomass
tots <- fit$tot_pop
if (is.null(K_groups)) {
tots$group <- "K"
} else {
tots$group <- paste0("K", "-", tots[, all.vars(K_groups)])
}
fit$pop %>%
ungroup() %>%
plot_ly(x = ~year, color = ~species, colors = viridis::viridis(100),
legendgroup = ~species) %>%
add_ribbons(ymin = ~B_lwr, ymax = ~B_upr, line = list(width = 0),
alpha = 0.2, showlegend = FALSE) %>%
add_lines(y = ~B) %>%
add_lines(data = tots, x = ~year, y = ~K,
linetype = ~group, color = I("black"), inherit = FALSE) %>%
add_labels(xlab = "Year", ylab = "Biomass")
Total biomass
three_cols <- viridis::viridis(3)
fit$tot_pop %>%
plot_ly(x = ~year, frame = K_groups) %>%
add_ribbons(ymin = ~B_lwr, ymax = ~B_upr, line = list(width = 0),
alpha = 0.2, showlegend = FALSE, legendgroup = "B",
color = I(three_cols[2])) %>%
add_lines(y = ~B, name = "B", color = I(three_cols[2]),
legendgroup = "B") %>%
add_lines(y = ~K, name = "K", legendgroup = "K",
linetype = I(1), color = I(three_cols[1])) %>%
add_lines(y = ~K_lwr, legendgroup = "K", showlegend = FALSE,
linetype = I(3), color = I(three_cols[1]), size = I(1)) %>%
add_lines(y = ~K_upr, legendgroup = "K", showlegend = FALSE,
linetype = I(3), color = I(three_cols[1]), size = I(1)) %>%
add_lines(y = ~landings, name = "L", color = I(three_cols[3])) %>%
add_labels(xlab = "Year", ylab = "Biomass")
Harvest rate
fit$pop %>%
plot_ly(x = ~year, color = ~species, colors = viridis::viridis(100),
legendgroup = ~species) %>%
add_ribbons(ymin = ~F_lwr, ymax = ~F_upr, line = list(width = 0),
alpha = 0.2, showlegend = FALSE) %>%
add_lines(y = ~F) %>%
add_labels(xlab = "Year", ylab = "Harvest Rate")
PaulRegular/MSP documentation built on Dec. 16, 2024, 1:59 p.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.
library(plotly) library(viridis) knitr::opts_chunk$set(echo = FALSE, fig.height = 7, fig.width = 8, cache = FALSE)
if (fit$call$K_groups == ~1) { K_groups <- NULL if (length(unique(fit$landings$species)) == 1) { K_label <- "K" } else { K_label <- "All species" } } else { K_groups <- as.formula(fit$call$K_groups) K_label <- unique(fit$landings[, all.vars(K_groups)]) } n_pe_covar <- length(all.vars(fit$call$pe_covariates)) n_K_covar <- length(all.vars(fit$call$K_covariates)) if (fit$call$species_cor != "none") { n_spp_cor <- length(fit$par$logit_rho) if (n_spp_cor > 1) { sp_rho <- sp_nm_mat <- matrix(NA, nrow = nlevels(fit$pop$species), ncol = nlevels(fit$pop$species)) rownames(sp_rho) <- colnames(sp_rho) <- levels(fit$pop$species) sp_rho[lower.tri(sp_rho)] <- inv_logit(fit$par$logit_rho, shift = TRUE) sp_rho <- t(sp_rho) sp_rho[lower.tri(sp_rho)] <- inv_logit(fit$par$logit_rho, shift = TRUE) diag(sp_rho) <- 1 round(sp_rho, 2) for (i in seq(nrow(sp_rho))) { for (j in seq(ncol(sp_rho))) { sp_nm_mat[i, j] <- paste(levels(fit$pop$species)[i], "-", levels(fit$pop$species)[j]) } } species_cor_label <- sp_nm_mat[lower.tri(sp_nm_mat)] } else { species_cor_label <- "All species" } } else { n_spp_cor <- 0 species_cor_label <- "N/A" } n_time_cor <- as.numeric(fit$call$temporal_cor == "ar1") temporal_cor_label <- ifelse(n_time_cor == 1, "All species", "N/A") add_labels <- function(p, xlab = "", ylab = "", log_linear_button = TRUE) { if (log_linear_button) { ## Modified from https://stackoverflow.com/a/60551263/4055438 buttons <- list(list( active = 0, type = "buttons", buttons= list( list(label = 'linear', method = 'relayout', args = list(list(yaxis = list(type = 'linear', title = ylab)))), list(label = 'log', method = 'relayout', args = list(list(yaxis = list(type = 'log', title = ylab))))))) } else { buttons <- NULL } layout(p, xaxis = list(title = xlab), yaxis = list(title = ylab), updatemenus = buttons) }
Inputs
Row {.tabset}
Index
fit$index %>% plot_ly() %>% add_trace(x = ~year, y = ~index, color = ~survey, colors = viridis::viridis(100), mode = "markers+lines", type = "scatter") %>% add_labels(xlab = "Year", ylab = "Index")
Landings
fit$landings %>% plot_ly() %>% add_trace(x = ~year, y = ~landings, color = ~species, colors = viridis::viridis(100), mode = "markers+lines", type = "scatter") %>% add_labels(xlab = "Year", ylab = "Landings")
Covariates
r if (n_pe_covar == 0 && n_K_covar == 0) "N/A; no process error or K covariates supplied"
r if ((n_pe_covar + n_K_covar) > 1) "Plotting of more than one covariate not yet implemented"
if (n_pe_covar == 1) { covar_name <- all.vars(fit$call$pe_covariates) covariates <- unique(fit$landings[, c("year", covar_name), drop = FALSE]) } if (n_K_covar == 1) { covar_name <- all.vars(fit$call$K_covariates) covariates <- unique(fit$landings[, c("year", covar_name), drop = FALSE]) } names(covariates) <- c("year", "covariate") covariates %>% plot_ly(x = ~year, y = ~covariate) %>% add_lines() %>% layout(xaxis = list(title = "Year"), yaxis = list(title = covar_name))
Residuals
Row {.tabset}
Residuals ~ predicted value
p <- fit$index %>% plot_ly(color = ~species, colors = viridis::viridis(100)) %>% layout(yaxis = list(title = "Standardized residuals")) p %>% add_markers(x = ~log(pred), y = ~std_res)
~ year
p %>% add_markers(x = ~year, y = ~std_res)
~ survey
p %>% add_markers(x = ~survey, y = ~std_res)
~ species
p %>% add_markers(x = ~species, y = ~std_res)
Parameters
Row {.tabset}
Correlation
sd_rep <- fit$sd_rep dsd <- sqrt(diag(sd_rep$cov.fixed)) cor_mat <- diag(1 / dsd) %*% sd_rep$cov.fixed %*% diag(1 / dsd) rownames(cor_mat) <- paste0(seq(nrow(cor_mat)), ":", names(dsd)) colnames(cor_mat) <- paste0(seq(nrow(cor_mat)), ":", names(dsd)) # corrplot::corrplot.mixed(cor_mat, diag = "n", lower = "ellipse", upper = "number") plot_ly(x = rownames(cor_mat), y = colnames(cor_mat), z = cor_mat) %>% add_heatmap(colors = c("#B2182B", "white", "#2166AC")) %>% colorbar(limits = c(-1.001, 1.001)) %>% layout(yaxis = list(scaleanchor = "x", scaleratio = 1), xaxis = list(constrain = "domain"))
log(K)
plot_post(prior_mean = fit$par$mean_log_K, prior_sd = exp(fit$par$log_sd_log_K), post_mean = fit$par$log_K, post_sd = fit$se$log_K, names = K_label, xlab = "log(K)", show_prior = fit$call$log_K_option$option == "prior")
log(r)
plot_post(prior_mean = fit$par$mean_log_r, prior_sd = exp(fit$par$log_sd_log_r), post_mean = fit$par$log_r, post_sd = fit$se$log_r, names = levels(fit$landings$species), xlab = "log(r)", show_prior = fit$call$log_r_option$option == "prior")
log(B0)
plot_post(prior_mean = fit$par$mean_log_B0, prior_sd = exp(fit$par$log_sd_log_B0), post_mean = fit$par$log_B0, post_sd = fit$se$log_B0, names = levels(fit$landings$species), xlab = "log(B0)", show_prior = fit$call$log_B0_option$option == "prior")
log(SDB)
plot_post(prior_mean = fit$par$mean_log_sd_B, prior_sd = exp(fit$par$log_sd_log_sd_B), post_mean = fit$par$log_sd_B, post_sd = fit$se$log_sd_B, names = levels(fit$landings$species), xlab = "log(SD<sub>B</sub>)", show_prior = fit$call$log_sd_B_option$option == "prior")
logit(rho)
r if (n_spp_cor == 0) "N/A; correlation of process errors across species was not estimated."
plot_post(prior_mean = fit$par$mean_logit_rho, prior_sd = exp(fit$par$log_sd_logit_rho), post_mean = fit$par$logit_rho, post_sd = fit$se$logit_rho, names = species_cor_label, xlab = "logit(rho)", # trans_fun = function(x) inv_logit(x, shift = TRUE), show_prior = fit$call$logit_rho_option$option == "prior")
logit(phi)
r if (n_time_cor == 0) "N/A; correlation of process errors across time was not estimated."
plot_post(prior_mean = fit$par$mean_logit_phi, prior_sd = exp(fit$par$log_sd_logit_phi), post_mean = fit$par$logit_phi, post_sd = fit$se$logit_phi, names = "logit(phi)", xlab = "logit(phi)", # trans_fun = inv_logit, show_prior = fit$call$logit_phi_option$option == "prior")
log(q)
plot_post(prior_mean = fit$par$mean_log_q, prior_sd = exp(fit$par$log_sd_log_q), post_mean = fit$par$log_q, post_sd = fit$se$log_q, names = levels(fit$index$survey), xlab = "log(q)", show_prior = fit$call$log_q_option$option == "prior")
log(SDI)
plot_post(prior_mean = fit$par$mean_log_sd_I, prior_sd = exp(fit$par$log_sd_log_sd_I), post_mean = fit$par$log_sd_I, post_sd = fit$se$log_sd_I, names = levels(fit$index$survey), xlab = "log(SD<sub>I</sub>)", show_prior = fit$call$log_sd_I_option$option == "prior")
betape
r if (n_pe_covar == 0) "N/A; no process error covariates supplied"
plot_post(prior_mean = fit$par$mean_pe_betas, prior_sd = exp(fit$par$log_sd_pe_betas), post_mean = fit$par$pe_betas, post_sd = fit$se$pe_betas, names = "pe_betas", xlab = "pe_betas", show_prior = fit$call$pe_betas_option$option == "prior")
betaK
r if (n_K_covar == 0) "N/A; no K covariates supplied"
plot_post(prior_mean = fit$par$mean_K_betas, prior_sd = exp(fit$par$log_sd_K_betas), post_mean = fit$par$K_betas, post_sd = fit$se$K_betas, names = "K_betas", xlab = "K_betas", show_prior = fit$call$K_betas_option$option == "prior")
Estimates (plot)
name_par <- c("log_K", "log_r", "log_B0", "log_sd_B", "logit_rho", "logit_phi", "log_q", "log_sd_I", "pe_betas", "K_betas") par_fun <- function(nm, select = "par") { if (grepl("beta", nm)) return(fit[[select]][[nm]]) if (grepl("log_", nm)) return(exp(fit[[select]][[nm]])) if (nm == "logit_rho") return(inv_logit(fit[[select]][[nm]], shift = TRUE)) if (nm == "logit_phi") return(inv_logit(fit[[select]][[nm]], shift = FALSE)) } named_par <- lapply(name_par, par_fun, select = "par") named_par_cv <- lapply(name_par, par_fun, select = "se") named_par_lwr <- lapply(name_par, par_fun, select = "par_lwr") named_par_upr <- lapply(name_par, par_fun, select = "par_upr") names(named_par) <- gsub("log_|logit_", "", name_par) names(named_par$q) <- levels(fit$index$survey) names(named_par$sd_I) <- levels(fit$index$survey) names(named_par$K) <- K_label names(named_par$r) <- levels(fit$index$species) names(named_par$sd_B) <- levels(fit$index$species) names(named_par$B0) <- levels(fit$index$species) names(named_par$phi) <- temporal_cor_label names(named_par$rho) <- species_cor_label names(named_par$pe_betas) <- colnames(fit$tmb_dat$pe_covariates) names(named_par$K_betas) <- colnames(fit$tmb_dat$K_covariates) par_tab <- data.frame(parameter = rep(names(named_par), sapply(named_par, length)), group = names(unlist(unname(named_par))), estimate = unlist(unname(named_par)), CV = unlist(named_par_cv), lower = unlist(named_par_lwr), upper = unlist(named_par_upr)) par_tab <- na.omit(par_tab) # will drop par lacking CV estimates (i.e. not estimated) par_tab %>% plot_ly(color = ~parameter, colors = viridis::viridis(100), legendgroup = ~parameter, y = ~paste(parameter, "-", group)) %>% add_segments(x = ~lower, xend = ~upper, yend = ~paste(parameter, "-", group), showlegend = FALSE) %>% add_markers(x = ~estimate) %>% layout(xaxis = list(title = ""), yaxis = list(title = ""))
Estimates (table)
par_tab$parameter[duplicated(par_tab$parameter)] <- "" knitr::kable(par_tab, digits = 3, row.names = FALSE)
Population trends
Row {.tabset}
Observed and predicted index
fit$index %>% plot_ly(x = ~year, color = ~survey, colors = viridis::viridis(100), legendgroup = ~survey) %>% add_ribbons(ymin = ~pred_lwr, ymax = ~pred_upr, line = list(width = 0), alpha = 0.2, showlegend = FALSE) %>% add_lines(y = ~pred) %>% add_markers(y = ~index, showlegend = FALSE) %>% add_labels(xlab = "Year", ylab = "Index")
Process error
p <- fit$pop %>% plot_ly(color = ~species, colors = viridis::viridis(100)) p %>% add_lines(x = ~year, y = ~pe) %>% add_labels(xlab = "Year", ylab = "Process error")
r if (n_pe_covar >= 1) "> Includes covariate effects"
Standardized process error
p <- fit$pop %>% plot_ly(color = ~species, colors = viridis::viridis(100)) p %>% add_lines(x = ~year, y = ~log_std_res_pe) %>% layout(xaxis = list(title = "Year"), yaxis = list(title = "Standardized process error"))
r if (n_pe_covar >= 1) "> Standardized process error = residual process error, in log space, not explained by the covariate effects divided by the standard deviation of the process." else "> Standardized process error = process error, in log space, divided by the standard deviation of the process."
Process error correlation
r if (n_spp_cor == 0) "N/A; correlation of process errors across species was not estimated."
r if (n_spp_cor == 1) "N/A; matrix of correlation was not estimated as rho was coupled across species."
plot_ly(x = rownames(sp_rho), y = rownames(sp_rho), z = ~sp_rho, colors = c("#B2182B", "white", "#2166AC")) %>% add_heatmap() %>% colorbar(limits = c(-1.001, 1.001), title = "ρ") %>% layout(yaxis = list(scaleanchor = "x", scaleratio = 1), xaxis = list(constrain = "domain"))
r if (n_spp_cor > 1 && n_pe_covar >= 1) "> Represents correlation in residual process error not explained by the covariate effects."
Biomass
tots <- fit$tot_pop if (is.null(K_groups)) { tots$group <- "K" } else { tots$group <- paste0("K", "-", tots[, all.vars(K_groups)]) } fit$pop %>% ungroup() %>% plot_ly(x = ~year, color = ~species, colors = viridis::viridis(100), legendgroup = ~species) %>% add_ribbons(ymin = ~B_lwr, ymax = ~B_upr, line = list(width = 0), alpha = 0.2, showlegend = FALSE) %>% add_lines(y = ~B) %>% add_lines(data = tots, x = ~year, y = ~K, linetype = ~group, color = I("black"), inherit = FALSE) %>% add_labels(xlab = "Year", ylab = "Biomass")
Total biomass
three_cols <- viridis::viridis(3) fit$tot_pop %>% plot_ly(x = ~year, frame = K_groups) %>% add_ribbons(ymin = ~B_lwr, ymax = ~B_upr, line = list(width = 0), alpha = 0.2, showlegend = FALSE, legendgroup = "B", color = I(three_cols[2])) %>% add_lines(y = ~B, name = "B", color = I(three_cols[2]), legendgroup = "B") %>% add_lines(y = ~K, name = "K", legendgroup = "K", linetype = I(1), color = I(three_cols[1])) %>% add_lines(y = ~K_lwr, legendgroup = "K", showlegend = FALSE, linetype = I(3), color = I(three_cols[1]), size = I(1)) %>% add_lines(y = ~K_upr, legendgroup = "K", showlegend = FALSE, linetype = I(3), color = I(three_cols[1]), size = I(1)) %>% add_lines(y = ~landings, name = "L", color = I(three_cols[3])) %>% add_labels(xlab = "Year", ylab = "Biomass")
Harvest rate
fit$pop %>% plot_ly(x = ~year, color = ~species, colors = viridis::viridis(100), legendgroup = ~species) %>% add_ribbons(ymin = ~F_lwr, ymax = ~F_upr, line = list(width = 0), alpha = 0.2, showlegend = FALSE) %>% add_lines(y = ~F) %>% add_labels(xlab = "Year", ylab = "Harvest Rate")
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.