| Title: | Fit Bayesian Models to Bycatch Data |
|---|---|
| Description: | Fits models to bycatch data, and can expand to fleet-wide estimates, like in situations when fisheries are observed with less than 100% observer coverage. |
| Authors: | Eric Ward [aut, cre], Jason Jannot [aut], Trustees of Columbia University [cph] |
| Maintainer: | Eric Ward <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.7 |
| Built: | 2024-10-09 05:32:48 UTC |
| Source: | https://github.com/ericward-noaa/bycatch |
A DESCRIPTION OF THE PACKAGE
Stan Development Team (2023). RStan: the R interface to Stan. R package version 2.21.8. https://mc-stan.org
fit_bycatch is the primary function for fitting bycatch models to time series of takes and effort
fit_bycatch( formula, data, time = "year", effort = "effort", expansion_rate = NULL, family = c("poisson", "nbinom2", "poisson-hurdle", "nbinom2-hurdle", "lognormal", "gamma", "lognormal-hurdle", "gamma-hurdle", "normal", "normal-hurdle"), time_varying = FALSE, iter = 1000, chains = 3, control = list(adapt_delta = 0.9, max_treedepth = 20), ... )fit_bycatch( formula, data, time = "year", effort = "effort", expansion_rate = NULL, family = c("poisson", "nbinom2", "poisson-hurdle", "nbinom2-hurdle", "lognormal", "gamma", "lognormal-hurdle", "gamma-hurdle", "normal", "normal-hurdle"), time_varying = FALSE, iter = 1000, chains = 3, control = list(adapt_delta = 0.9, max_treedepth = 20), ... )
formula |
The model formula. |
data |
A data frame. |
time |
Named column of the 'data' data frame with the label for the time (e.g. year) variable |
effort |
Named column of the 'effort' variable in the data frame with the label for the fishing effort to be used in estimation of mean bycatch rate. This represents total observed effort |
expansion_rate |
The expansion rate to be used in generating distributions for unobserved sets. If NULL, defaults to 100% coverage (= 100) |
family |
Family for response distribution can be discrete ("poisson", "nbinom2", "poisson-hurdle","nbinom2-hurdle"), or continuous ("normal", "gamma","lognormal", "normal-hurdle", "gamma-hurdle", "lognormal-hurdle"). The default distribution is "poisson". The hurdle variants estimate the probability of zeros (theta) separately from the other models and use truncated distribution to model positive counts. All use a log link function. |
time_varying |
boolean TRUE/FALSE, whether to include time varying component (this is a random walk, analogous to making this a Dynamic linear model) |
iter |
the number of mcmc iterations, defaults to 1000 |
chains |
the number of mcmc chains, defaults to 3 |
control |
List to pass to |
... |
Any other arguments to pass to |
list of the data used to fit the model, the matrix of covariates, the expanded bycatch generated via the fit and simulations, and the fitted stan model
d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) loo::loo(fit$fitted_model)$estimates fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) plot_fitted(fit, xlab = "Year", ylab = "Fleet-level bycatch", include_points = TRUE ) # fit a negative binomial model, with more chains and control arguments fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # fit a time varying model fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = TRUE ) # include data for expansion to unobserved sets fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # fit a model with a lognormal distribution d$Takes <- rnorm(nrow(d), 5, 0.1) fit_ln <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "lognormal", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # add zeros and fit a delta-gamma distribution d$Takes <- rnorm(nrow(d), 5, 0.1) d$Takes[c(1, 5, 10)] <- 0 fit_ln <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "gamma-hurdle", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) )d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) loo::loo(fit$fitted_model)$estimates fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) plot_fitted(fit, xlab = "Year", ylab = "Fleet-level bycatch", include_points = TRUE ) # fit a negative binomial model, with more chains and control arguments fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # fit a time varying model fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = TRUE ) # include data for expansion to unobserved sets fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # fit a model with a lognormal distribution d$Takes <- rnorm(nrow(d), 5, 0.1) fit_ln <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "lognormal", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # add zeros and fit a delta-gamma distribution d$Takes <- rnorm(nrow(d), 5, 0.1) d$Takes[c(1, 5, 10)] <- 0 fit_ln <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "gamma-hurdle", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) )
get_expanded is a helper function to return a matrix of posterior predictive values for unobserved bycatch
get_expanded(fitted_model)get_expanded(fitted_model)
fitted_model |
Data and fitted model returned from estimation |
matrix (MCMC draws x time steps) of posterior predictive values for unobserved bycatch
d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", expansion_rate = "expansionRate", time_varying = FALSE ) expanded <- get_expanded(fit)d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", expansion_rate = "expansionRate", time_varying = FALSE ) expanded <- get_expanded(fit)
get_fitted returns df of observed bycatch estimates (lambda of Poisson), accounting for effort but not accounting for observer coverage
get_fitted(fitted_model, alpha = 0.05)get_fitted(fitted_model, alpha = 0.05)
fitted_model |
Data and fitted model returned from fit_bycatch(). If a hurdle model, then only then the plot returns the total bycatch rate (including zero and non-zero components). |
alpha |
The alpha level for the credible interval, defaults to 0.05 |
plot called from ggplot
d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) get_fitted(fit)d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) get_fitted(fit)
get_total is a helper function to return a matrix of total estimated bycatch
get_total(fitted_model)get_total(fitted_model)
fitted_model |
Data and fitted model returned from estimation |
matrix (MCMC draws x time steps) of posterior predictive values for unobserved bycatch
d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", expansion_rate = "expansionRate", time_varying = FALSE ) expanded <- get_total(fit)d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", expansion_rate = "expansionRate", time_varying = FALSE ) expanded <- get_total(fit)
plot_expanded is makes plots of the expanded bycatch estimates, accounting for observer coverage and effort
plot_expanded( fitted_model, xlab = "Time", ylab = "Events", show_total = TRUE, include_points = FALSE )plot_expanded( fitted_model, xlab = "Time", ylab = "Events", show_total = TRUE, include_points = FALSE )
fitted_model |
Data and fitted model returned from estimation |
xlab |
X-axis label for plot |
ylab |
Y-axis label for plot |
show_total |
Whether to show the total predicted bycatch (by default, this is TRUE) or just the expanded unobserved events (=FALSE) |
include_points |
whether or not to include raw bycatch events on plots, defaults to FALSE |
plot called from ggplot
d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", expansion_rate = "expansionRate", time_varying = FALSE ) plot_expanded( fitted_model = fit, xlab = "Year", ylab = "Fleet-level bycatch", include_points = TRUE )d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", expansion_rate = "expansionRate", time_varying = FALSE ) plot_expanded( fitted_model = fit, xlab = "Year", ylab = "Fleet-level bycatch", include_points = TRUE )
plot_fitted makes plots bycatch estimates (lambda of Poisson), accounting for effort but not accounting for observer coverage
plot_fitted( fitted_model, xlab = "Time", ylab = "Events", include_points = FALSE, alpha = 0.05 )plot_fitted( fitted_model, xlab = "Time", ylab = "Events", include_points = FALSE, alpha = 0.05 )
fitted_model |
Data and fitted model returned from fit_bycatch(). If a hurdle model, then only then the plot returns the total bycatch rate (including zero and non-zero components). |
xlab |
X-axis label for plot |
ylab |
Y-axis label for plot |
include_points |
whether or not to include raw bycatch events on plots, defaults to FALSE |
alpha |
The alpha level for the credible interval, defaults to 0.05 |
plot called from ggplot
d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) plot_fitted(fit, xlab = "Year", ylab = "Fleet-level bycatch", include_points = TRUE ) # fit a negative binomial model, with more chains and control arguments fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # fit a time varying model fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = TRUE ) # include data for expansion to unobserved sets fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) )d <- data.frame( "Year" = 2002:2014, "Takes" = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0), "expansionRate" = c(24, 22, 14, 32, 28, 25, 30, 7, 26, 21, 22, 23, 27), "Sets" = c(391, 340, 330, 660, 470, 500, 330, 287, 756, 673, 532, 351, 486) ) fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = FALSE ) plot_fitted(fit, xlab = "Year", ylab = "Fleet-level bycatch", include_points = TRUE ) # fit a negative binomial model, with more chains and control arguments fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) ) # fit a time varying model fit <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "poisson", time_varying = TRUE ) # include data for expansion to unobserved sets fit_nb <- fit_bycatch(Takes ~ 1, data = d, time = "Year", effort = "Sets", family = "nbinom2", expansion_rate = "expansionRate", time_varying = FALSE, iter = 2000, chains = 4, control = list(adapt_delta = 0.99, max_treedepth = 20) )