vignettes/background.Rmd
- background.RmdPrevious authors, including Gardner et al. (2008) and Martin et al. (2015) have assumed that rare even bycatch follows a Poisson process, where the observed bycatch events \(y_{t}\) are modeled according to some estimated bycatch rate \(\lambda\),
-\[p(y_{t}|\lambda) = e^{-\lambda}\frac{\lambda^{y_{t}}}{y_{t}!}\] where \(\lambda\) is the mean of the Poisson distribution. The mean rate parameter \(\lambda\) can be further decomposed into a per - event (e.g. set) parameter and the number of observed sets \(n_{t}\), \(\lambda = \theta * n_{t}\). In a GLM setting, the log-link function is often used to model effects of covariates, e.g. \(log(\lambda) = log(\theta) + b_{1}*x_{1} + b_{2}*x_{2}...\), and in this setting the known number of sets \(n_{t}\) are treated as an offset, \(log(\lambda) = log(\theta) + b_{1}*x_{1} + b_{2}*x_{2} + log(n_{t})\).
-Estimation of \(\lambda\) or \(\theta\) can be done in a maximum likelihood or Bayesian setting, and here we use Bayesian methods to estimate the posterior distribution of the parameters given observed bycatch data,
-\[p(\theta|\underline{\mathbf{y}}) \propto p(\theta)p(\underline{\mathbf{y}}|\theta)\]
-In addition to the Poisson model described here, we include the Negative Binomial model, which allows for the variance to be greater than the mean. We use the ‘nbinom2’ alternative parameterization of the Negative Binomial, where \(Var(Y) = u + \frac{u^2}{\tau}\), where \(\tau\) controls the degree of overdispersion. Two additional extensions to these models are the zero-inflated extensions. We adopt hurdle models to model extra zeros, so that \(p(y_{t}=0) = \phi\) and \(p(y_{t}>0) = 1-\phi\).
-As described above, any type of fixed effects can be included as a predictors of bycatch in the formula interface. Coefficients are estimated for all coefficients using the log-link, \(log(\lambda) = log(\theta) + b_{1}*x_{1} + b_{2}*x_{2} + log(n_{t})\), where \(x_{1}\) and \(x_{2}\) are covariates.
-In some cases, such as when covariates may be missing or not completely observed, it may be useful to include time effects as time - varying. We have implemented these as random effects, in link (log) space. We can modify the predictions of bycatch rate to be \(log(\lambda) = log(\theta) + b_{1}*x_{1} + b_{2}*x_{2} + log(n_{t}) + \omega_{t}\), where \(\omega_{t}\) is the latent time effect. We constrain the time effects to be a random walk, so that \(\omega_{t} \sim Normal(\omega_{t-1},\sigma_{\omega})\). Extensions of this model could include autoregressive coefficients to help aid in stationarity.
-A major goal of bycatch analyses with partially observed datasets is estimating the total unobserved bycatch. Because the total numbers of effort (sets, trips, etc) is known, and the number of observed units of effort is known, the numbers of unobserved units is also known, \(N_{t}-n_{t}\).
-Using the Bayesian framework, we can generate samples from the posterior predictive distribution of unobserved sets, using the same distribution that is assumed for the observed sets (Poisson, Negative Binomial, etc). This posterior predictive approach assumes that the distribution of unobserved bycatch is
-\[P(Y_{t}-y_{t} | \theta, N_{t}, n_{t}) = \int_{\theta} p(Y_{t}-y_{t} | \theta, N_{t} - n_{t})p(\theta | \underline{\mathbf{y}})d\theta\]
-The integration includes a product of two quantities. The first is the probability of individual bycatch values (e.g. 0, 1, 2, …) conditioned on the estimated parameters \(\theta\) and number of unobserved units \(N_{t}-n_{t}\), and the second represents the posterior distribution of the parameters given all data, \(p(\theta | \underline{\mathbf{y}})\). By generating large numbers of samples from the posterior predictive distribution, we can generate credible intervals and other summary statistics on the distribution.
-vignettes/expanding.Rmd
- expanding.Rmd
-#library(devtools)
-#devtools::install_github("ericward-noaa/bycatch")
-library(bycatch)
-set.seed(123)Previous authors, including Gardner et al. (2008) and Martin et al. (2015) have assumed that rare even bycatch follows a Poisson process, where the observed bycatch events \(y_{t}\) are modeled according to some estimated bycatch rate \(\lambda\),
-\[p(y_{t}|\lambda) = e^{-\lambda}\frac{\lambda^{y_{t}}}{y_{t}!}\] where \(\lambda\) is the mean of the Poisson distribution. The mean rate parameter \(\lambda\) can be further decomposed into a per - event (e.g. set) parameter and the number of observed sets \(n_{t}\), \(\lambda = \theta * n_{t}\). In a GLM setting, the log-link function is often used to model effects of covariates, e.g. \(log(\lambda) = log(\theta) + B_{1}*x_{1} + B_{2}*x_{2}...\), and in this setting the known number of sets \(n_{t}\) are treated as an offset, \(log(\lambda) = log(\theta) + B_{1}*x_{1} + B_{2}*x_{2} + log(n_{t})\).
-The reported effort (here, sets) represents the total events that are observed. Focusing on a single year of data, the expansion rate, \(p_{obs}\), represents the fraction (n / 100) of the sets that are observed. The total sets are also assumed known, \(N = e_{obs} + e_{unobs}\). So we are interested in estimating \(e_{unobs}\) given \(N\) and \(e_{obs}\), where \(N = e_{obs} / p_{obs}\) and \(e_{unobs} = e_{obs} * (1/p_{obs} - 1)\).
-
-# replace this with your own data frame
-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))We’ll start by fitting a model with constant bycatch rate,
-
-fit = fit_bycatch(Takes ~ 1, data=d, time="Year", effort="Sets", family="poisson",
- time_varying = FALSE)Using our example above, the observer coverage for that dataset was less than 100% and so our estimates need to be expanded to the fleetwide level. There are some important control arguments here that are left at defaults, but should maybe be changed for huge numbers of bycatch.
-expanded = expand(fit, coverage = d$expansionRate)And we can then plot these estimates. Like the previous function we can specify whether to include the raw points or not.
-
-plot_expanded(fitted_model=fit, expanded_estimates = expanded, xlab="Year", ylab = "Fleet-level bycatch", include_points = TRUE)We can also do things like summarize the expanded estimates in table form
-
-df = data.frame("time" = d[,"Year"],
- "mean" = apply(expanded, 2, mean),
- "median" = apply(expanded, 2, quantile, 0.5),
- "lower95" = apply(expanded, 2, quantile, 0.025),
- "upper95" = apply(expanded, 2, quantile, 0.975))
-
-write.table(df, "estimated_bycatch.csv", row.names=F, col.names=T, sep=",")Stan Development Team (2018). RStan: the R interface to Stan. R package version 2.18.2. http://mc-stan.org
+Stan Development Team (2023). RStan: the R interface to Stan. R package version 2.21.8. https://mc-stan.org