knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width=12, fig.height=8)
library(hookCompetition)
library(gfiphc)
library(sf)
library(inlabru)
library(ggpolypath)
library(tidyverse)
colsc <- function(...) {
  ggplot2::scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11,"RdYlBu")),
                       limits = range(..., na.rm=TRUE))
}

How to use the hookCompetition package

This markdown document will outline the steps needed to use the hookCompetition package.

Step 1 - Setting parameters for generating the Rmarkdown document

# Set the number of cores to use to generate the vignette
n_cores <- 1
# Set the number of samples to use for bootstrapping for Monte Carlo approximations to
# posterior quantities. Should be set much higher for real applications!
n_samp <- 100

Step 2 - Loading the data.

The function read_Data_hookcomp() reads in longline catch count data collected by the IPHC. The arguments to the function are:

data <- hookCompetition::read_Data_hookcomp(
  species_vec = c("yelloweye rockfish",
            "arrowtooth flounder",
            "lingcod",
            "north pacific spiny dogfish",
            "pacific cod",
            "pacific halibut",
            "redbanded rockfish",
            "sablefish",
            "walleye pollock"),
  simple_species_vec = c("yelloweye",
            "arrowtooth",
            "lingcod",
            "dogfish",
            "cod",
            "halibut",
            "redbanded",
            "sablefish",
            "walleye"),
  at_PBS = F,
  data_wd = '../',
  min_dist = 50,
  min_dist_to_boundary = 0,
  years_all_vec = c(1995,1996,2003:2012, 2014:2019),
  years_20_vec = c(1997:2002, 2013, 2020, 2021),
  survey_boundaries=NULL
)

The function read_Data_hookcomp() returns a list containing two entries: - Reduced_data_sp is a sf points containing all the set-level catch counts of the species specified (variable name N_it_yelloweye). - survey_boundaries are the survey boundaries projected into the same CRS as the Reduced_data_sp object.

To put the objects in the list into the current working environment, run the following:

list2env(data, globalenv())
rm(data)

Note: to use the package on custom data, ensure the following variable names are used inside the sf points object containing the catch counts:

Step 3 - Fit the indices without the spatio-temporal random effects.

We derive bootstrap indices for a set of test species using the function bootstrap_index_fun_run. We derive an overall index and a subregion index (with subregions defined by the survey boundaries).

The function bootstrap_index_fun_run takes as arguments:

bootstrap_ind <-
bootstrap_index_fun_run(species = 'yelloweye',
                          data = Reduced_data_sp,
                          survey_boundaries = survey_boundaries,
                          R=n_samp,
                          return=T,
                          ncpus=n_cores,
                          type='perc',
                          plot=T
                          )

The function bootstrap_index_fun_run returns a list containing 3 objects:

We derive hook-competition-adjusted bootstrap indices using the same function bootstrap_index_fun_run(), but adding the argument ICR_adjust=TRUE. We derive an overall index and a subregion index (with subregions defined by the survey boundaries).

bootstrap_ind_compfactor <-
bootstrap_index_fun_run(species = 'yelloweye',
                          data = Reduced_data_sp,
                          ICR_adjust=T,
                          survey_boundaries = survey_boundaries,
                          R=n_samp,
                          return=T,
                          ncpus=n_cores,
                          type='perc',
                          plot=T
                          )

Finally, we compute Poisson-lognormal model-based indices of relative abundance. These include the CPUE-based (unadjusted), hook-competition-adjusted (ICR-based), and censored likelihood-based indices.

To compute all three indices we use the function censored_index_fun. This function takes many arguments, many of which are optional. The non-optional arguments are:

The optional arguments are:

We now fit all three indices

CPUE_ind <-
censored_index_fun(
  data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=1.1, nthreads=n_cores, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = c(0.228, -2.346, 2.407, -0.586), n_knots=8, seed=0, verbose=F, n_trajectories=10,
  preserve_inter_regional_differences = F
)

ICR_ind <-
censored_index_fun(
  data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=T, cprop=1.1, nthreads=n_cores, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = NULL, n_knots=8, seed=0, verbose=F, n_trajectories=10,
  preserve_inter_regional_differences = F
)

# The following model does not converge
# censored_ind <-
# censored_index_fun(
#   data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=0.95, nthreads=n_cores, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = NULL, n_knots=8, seed=0, verbose=F, n_trajectories=10
#     )

# Use upper bound
# Model now converges
censored_ind <-
censored_index_fun(
  data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=0.95, nthreads=n_cores, keep=F, use_upper_bound=TRUE, upper_bound_quantile=0.95, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = c(0.228, -2.346, 2.407, -0.586), n_knots=8, seed=0, verbose=F, n_trajectories=10,
  preserve_inter_regional_differences = F
    )
# We don't run the models in the vignette as INLA crashes GitHub Actions
# Instead, we load pre-made files here
list2env(hookCompetition::vignettes_precompiled_data, globalenv())

ggplot2::ggplot(CPUE_inddf, ggplot2::aes(x=.data$year ,y=.data$mean, ymin=.data$q0.025, ymax=.data$q0.975)) +
        ggplot2::geom_point() + ggplot2::geom_errorbar() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~.data$region) +
        ggplot2::ggtitle(paste0('Overdispersed ','Poisson Index ','yelloweye'))

ggplot2::ggplot(ICR_inddf, ggplot2::aes(x=.data$year ,y=.data$mean, ymin=.data$q0.025, ymax=.data$q0.975)) +
        ggplot2::geom_point() + ggplot2::geom_errorbar() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~.data$region) +
        ggplot2::ggtitle(paste0('Overdispersed ','ICR-Adjusted' ,' Poisson Index ','yelloweye'))

ggplot2::ggplot(censored_inddf, ggplot2::aes(x=.data$year ,y=.data$mean, ymin=.data$q0.025, ymax=.data$q0.975)) +
        ggplot2::geom_point() + ggplot2::geom_errorbar() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~.data$region) +
        ggplot2::ggtitle(paste0('Overdispersed ','Censored' ,' Poisson Index ','yelloweye'), subtitle = paste0('censorship proportion ',0.95, ', censorship from data in upper ',max(0,(1-0.95)*100), '% of values removed'))

ggplot2::ggplot(CPUE_trajectory, ggplot2::aes(x=.data$year ,y=.data$mean, colour=.data$MC_ind, group=.data$MC_ind)) +
          ggplot2::geom_line() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~.data$region) +
          ggplot2::ggtitle(paste0('Overdispersed ',' Poisson Index ','yelloweye')) +
          viridis::scale_color_viridis()

ggplot2::ggplot(ICR_trajectory, ggplot2::aes(x=.data$year ,y=.data$mean, colour=.data$MC_ind, group=.data$MC_ind)) +
          ggplot2::geom_line() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~.data$region) +
          ggplot2::ggtitle(paste0('Overdispersed ','ICR-Adjusted', ' Poisson Index ','yelloweye')) +
          viridis::scale_color_viridis()

ggplot2::ggplot(censored_trajectory, ggplot2::aes(x=.data$year ,y=.data$mean, colour=.data$MC_ind, group=.data$MC_ind)) +
          ggplot2::geom_line() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~.data$region) +
          ggplot2::ggtitle(paste0('Overdispersed ','Censored' ,' Poisson Index ','yelloweye'),subtitle = paste0('censorship proportion ',0.95, ', censorship from data in upper ',max(0,(1-0.95)*100), '% of values removed')) +
          viridis::scale_color_viridis()

The function censored_index_fun() returns a list containing up to 8 entries:

Note that the censored Poisson indices can also be fit using the sdmTMB package using the function censored_index_fun_sdmTMB. The arguments to the function are similar, with the exclusion of priors and n_knots due to penalized splines or first-order random walk structures being placed on the temporal effects. If the relative abundance trend is assumed to be smooth - perhaps for modelling a long-living species with small home range, then set the argument smoothing_spline=TRUE. If set to FALSE, then a random walk temporal structure is assumed. This allows for larger changes to occur each year (e.g. Pacific Hake summer abundance which changes dramatically from year-year due to summer migration).

Note that the lack of priors penalizing the likelihood on the IID random effects means that we must remove 1996 and 1997 from consideration within the censored Poisson model to make it converge.

We now fit all three indices within sdmTMB

CPUE_ind_sdmTMB <-
censored_index_fun_sdmTMB(
  data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=1.1, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, seed=0, verbose=F, n_trajectories=10,
  preserve_inter_regional_differences = F, time_effect = 'spline'
)

ICR_ind_sdmTMB <-
censored_index_fun_sdmTMB(
  data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=T, cprop=1.1, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, seed=0, verbose=F, n_trajectories=10,
  preserve_inter_regional_differences = F, time_effect = 'spline'
)

censored_ind_sdmTMB <-
censored_index_fun_sdmTMB(
  data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=0.95, keep=F, use_upper_bound=TRUE, upper_bound_quantile=0.95, plot=T, allyears=F, station_effects=T, seed=0, verbose=F, n_trajectories=10,
  preserve_inter_regional_differences = F, time_effect = 'spline'
    )

Bonus

The package also allows for spatio-temporal indices to be fit

Step 5 load the R objects required for fitting spatio-temporal models. These include: a delauney triangulation mesh for computing the GMRF approximation to the SPDE (cite Lindgren) within r-INLA (cite Rue); an INLA spde object with user-specified pc priors specified on the range and standard deviation parameters. Note that the spde object is a barrier model (cite Bakka) which accounts for land as a barrier.

The function make_spatial_objects() achieves this. It takes required arguments:

It takes optional arguments:

spatial_obj <- make_spatial_objects(
  data=Reduced_data_sp,
  survey_boundaries = survey_boundaries,
  mesh_coast = NULL, hires_coast = NULL,
  prior_range = 20, prior_range_prob = 0.01,
  prior_sigma = 2, prior_sigma_prob = 0.01
)

# load it into the global envirionment
list2env(spatial_obj, globalenv())
rm(spatial_obj)

ggplot2::ggplot() + inlabru::gg(mesh) + inlabru::gg(as_Spatial(Reduced_data_sp))

The function make_spatial_objects() returns a list containing X entries:

Step 6 Fit the spatio-temporal indices

The following code would be run to make spatio-temporal indices of yelloweye abundance. We do not run the code to save time whilst generating this vignette. Instead, we load pre-compiled objects. Note, we removed map predictions from all years except 2000 to save space.

The function spatiotemp_censored_index_fun() is used to fit the spatio-temporal indices. It takes required arguments (the new ones are described):

The optional arguments are:

Note that these models can take a while to run (5-30 minutes). We advise you to set verbose=T so that you can assure yourself that something useful is happening behind the scenes :)

spatio_temporal_CPUE_ind <-
  spatiotemp_censored_index_fun(data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=1.1, nthreads=n_cores, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = c(0.385, -2.083, 2.711, -0.861, 4.000, 0.200, 3.400, 0.300, 2.938), n_knots=8, seed=0, mesh=mesh, spde=spde_mod, pixels=predict_pixels, covs=NULL, spatiotemporal=T, verbose = F, n_trajectories = 10,
  preserve_inter_regional_differences = F)


spatio_temporal_ICR_ind <-
  spatiotemp_censored_index_fun(data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=T, cprop=1.1, nthreads=n_cores, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = NULL, n_knots=8, seed=0, mesh=mesh, spde=spde_mod, pixels=predict_pixels, covs=NULL, spatiotemporal=T, verbose = F, n_trajectories = 10,
  preserve_inter_regional_differences = F)

spatio_temporal_censored_ind <-
  spatiotemp_censored_index_fun(data=Reduced_data_sp, survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=0.95, nthreads=n_cores, keep=F, use_upper_bound=TRUE, upper_bound_quantile=0.95, plot=T, allyears=F, station_effects=T, prior_event=HT.prior(), prior_station=HT.prior(), init_vals = NULL, n_knots=8, seed=0, mesh=mesh, spde=spde_mod, pixels=predict_pixels, covs=NULL, spatiotemporal=T, verbose = F, n_trajectories = 10,
  preserve_inter_regional_differences = F)

The function spatiotemp_censored_index_fun() returns a list containing up to 8 entries:

Now we plot maps of yelloweye density from the three spatio-temporal models in the year 2000 using the function spatio_temporal_plot_fun() which takes arguments:

# plot the posterior mean yelloweye density on a map in year 2000
spatio_temporal_plot_fun(
  spatio_temporal_CPUE_plotdf, year=2000, variable='mean',
  hires_COAST=Coast_hires, plot_figure=T, return_figure=F
)

# plot the posterior mean yelloweye density on a map in year 2000
spatio_temporal_plot_fun(
  spatio_temporal_ICR_plotdf, year=2000, variable='mean',
  hires_COAST=Coast_hires, plot_figure=T, return_figure=F
)

# plot the posterior mean yelloweye density on a map in year 2000
spatio_temporal_plot_fun(
  spatio_temporal_censored_plotdf, year=2000, variable='mean',
  hires_COAST=Coast_hires, plot_figure=T, return_figure=F
)
# We don't run the models in the vignette as INLA crashes GitHub Actions
# Instead, we load pre-made files here
ggplot2::ggplot(spatio_temporal_CPUE_inddf, ggplot2::aes(x=year ,y=mean, ymin=q0.025, ymax=q0.975)) +
        ggplot2::geom_point() + ggplot2::geom_errorbar() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~region) +
        ggplot2::ggtitle(paste0('Spatiotemporal Overdispersed Poisson Index ','Yelloweye'))

ggplot2::ggplot(spatio_temporal_ICR_inddf, ggplot2::aes(x=year ,y=mean, ymin=q0.025, ymax=q0.975)) +
        ggplot2::geom_point() + ggplot2::geom_errorbar() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~region) +
                ggplot2::ggtitle(paste0('Spatiotemporal Overdispersed ICR-based Poisson Index ','Yelloweye'))

ggplot2::ggplot(spatio_temporal_censored_inddf, ggplot2::aes(x=year ,y=mean, ymin=q0.025, ymax=q0.975)) +
        ggplot2::geom_point() + ggplot2::geom_errorbar() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~region) +
        ggplot2::ggtitle(paste0('Spatiotemporal Overdispersed Censored Poisson Index yelloweye'), subtitle = 'censorship proportion 0.9, censorship from data in upper 10% of values removed')

ggplot2::ggplot(spatio_temporal_CPUE_trajectory, ggplot2::aes(x=year ,y=mean, colour=MC_ind, group=MC_ind)) +
          ggplot2::geom_line() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~region) +
          ggplot2::ggtitle(paste0('Overdispersed',' Poisson Index ','yelloweye')) +
          viridis::scale_color_viridis()

ggplot2::ggplot(spatio_temporal_ICR_trajectory, ggplot2::aes(x=year ,y=mean, colour=MC_ind, group=MC_ind)) +
          ggplot2::geom_line() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~region) +
          ggplot2::ggtitle(paste0('Overdispersed ICR-based',' Poisson Index ','yelloweye')) +
          viridis::scale_color_viridis()

ggplot2::ggplot(spatio_temporal_censored_trajectory, ggplot2::aes(x=year ,y=mean, colour=MC_ind, group=MC_ind)) +
          ggplot2::geom_line() + ggplot2::ylab('Catch rate index') + ggplot2::facet_grid(~region) +
          ggplot2::ggtitle(paste0('Overdispersed censored',' Poisson Index ','yelloweye'),
                  subtitle = 'Censorship defined as 95% or more baits removed \nCensorship removed from upper 10% of values each year') +          viridis::scale_color_viridis()
# Plot the pre-made files
ggplot2::ggplot() + inlabru::gg(spatio_temporal_CPUE_plotdf[which(spatio_temporal_CPUE_plotdf$year==2000), c('mean')]) + inlabru::gg(sf::as_Spatial(Coast_hires)) + colsc(as.matrix(spatio_temporal_CPUE_plotdf[which(spatio_temporal_CPUE_plotdf$year==2000), c('mean')])) + ggplot2::ggtitle('CPUE-based yelloweye abundance in year 2000')

ggplot2::ggplot() + inlabru::gg(spatio_temporal_ICR_plotdf[which(spatio_temporal_ICR_plotdf$year==2000), c('mean')]) + inlabru::gg(sf::as_Spatial(Coast_hires)) + colsc(as.matrix(spatio_temporal_ICR_plotdf[which(spatio_temporal_ICR_plotdf$year==2000), c('mean')])) + ggplot2::ggtitle('ICR-based yelloweye abundance in year 2000')

ggplot2::ggplot() + inlabru::gg(spatio_temporal_censored_plotdf[which(spatio_temporal_censored_plotdf$year==2000), c('mean')]) + inlabru::gg(sf::as_Spatial(Coast_hires)) + colsc(as.matrix(spatio_temporal_censored_plotdf[which(spatio_temporal_censored_plotdf$year==2000), c('mean')])) + ggplot2::ggtitle('Censored Poisson yelloweye abundance in year 2000',subtitle = 'Censorship defined as 95% or more baits removed \nCensorship removed from upper 10% of values each year')

Note that we can also fit spatiotemporal models in sdmTMB using spatio_temporal_censored_ind_sdmTMB(). It takes the same arguments with the exception of priors.

spatio_temporal_CPUE_ind_sdmTMB <-
  spatiotemp_censored_index_fun_sdmTMB(data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=1.1, keep=T, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, seed=0, mesh=mesh, pixels=predict_pixels, covs=NULL, spatiotemporal='AR1', verbose = T, n_trajectories = 10, preserve_inter_regional_differences = F, grf_priors = sdmTMBpriors(
  matern_s = pc_matern(range_gt = 20, sigma_lt = 2),
  matern_st = pc_matern(range_gt = 20, sigma_lt = 2)),
  time_effect = 'unstructured')


spatio_temporal_ICR_ind_sdmTMB <-
  spatiotemp_censored_index_fun_sdmTMB(data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=T, cprop=1.1, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1.1, plot=T, allyears=F, station_effects=T, seed=0, mesh=mesh, pixels=predict_pixels, covs=NULL, spatiotemporal='AR1', verbose = T, n_trajectories = 10,
  preserve_inter_regional_differences = F, grf_priors = sdmTMBpriors(
  matern_s = pc_matern(range_gt = 20, sigma_lt = 2),
  matern_st = pc_matern(range_gt = 20, sigma_lt = 2)),
  prev_fit = spatio_temporal_CPUE_ind_sdmTMB$mod,
  time_effect = 'unstructured')

spatio_temporal_censored_ind_sdmTMB <-
  spatiotemp_censored_index_fun_sdmTMB(data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=0.95, keep=F, use_upper_bound=TRUE, upper_bound_quantile=0.95, plot=T, allyears=F, station_effects=T, seed=0, mesh=mesh, pixels=predict_pixels, covs=NULL, spatiotemporal='AR1', verbose = T, n_trajectories = 10,
  preserve_inter_regional_differences = F, grf_priors = sdmTMBpriors(
  matern_s = pc_matern(range_gt = 20, sigma_lt = 2),
  matern_st = pc_matern(range_gt = 20, sigma_lt = 2)),
  prev_fit = spatio_temporal_CPUE_ind_sdmTMB$mod,
  time_effect = 'unstructured')

# Remove the upper bound and upper quantile
spatio_temporal_censored_ind_sdmTMB2 <-
    spatiotemp_censored_index_fun_sdmTMB(data=Reduced_data_sp %>% filter(year>=1998), survey_boundaries=survey_boundaries, species='yelloweye', M=n_samp, return=T, ICR_adjust=F, cprop=0.95, keep=F, use_upper_bound=FALSE, upper_bound_quantile=1, plot=T, allyears=F, station_effects=T, seed=0, mesh=mesh, pixels=predict_pixels, covs=NULL, spatiotemporal='AR1', verbose = T, n_trajectories = 10,
                                         preserve_inter_regional_differences = F,
                                         grf_priors = sdmTMBpriors(
                                             matern_s = pc_matern(range_gt = 20, sigma_lt = 2),
                                             matern_st = pc_matern(range_gt = 20, sigma_lt = 2)),
                                         prev_fit = spatio_temporal_CPUE_ind_sdmTMB$mod,
  time_effect = 'unstructured')


pbs-assess/hookCompetition documentation built on April 27, 2023, 12:47 p.m.