Estimate catchabilities from combined data

Specify where to save outputs:

output_path <- "./fit_model_catchability"
make_folder(output_path, recursive = TRUE)

## directory contains: fit_cumulative_captures_barchart_capture_code.png, fit_cumulative_captures_point_capture_code.png, fit_empirical_captures_barchart_capture_code.png, fit_empirical_captures_point_capture_code.png, log10_catchabilities.png, log10_catchabilities_boxplot.png, log10_catchability_coef_fishery_group.png, pi_boxplot.png, prior_updates_N_BP.png, prior_updates_P_B.png, trace_plot.png, deaths.png, deaths_caption.tex, density_overlap.png, density_overlap_caption.tex, mdl_ini.rda, out.rda, sefra_data.rda, catchability_coefficients_fishery_group.tex, fit_captures_per_code_fishery_group.tex, fit_captures_per_species_group_fishery_group.tex, q_per_species_group_fishery_group.tex, q_per_species_group_fishery_group_alternative.tex

Create `sefraData` object

First, identify capture codes to be used in model:

mdl_codes <- sp_groups %>% filter(., ! code %in% codes_exclude)

Initialise the data object:

sefra_data <- sefraData(mdl_codes$code, identifier = "a_2024_SRA")

## capture codes input: including all contributing species

## constructed 'sefraData' object

species_names(sefra_data)

## Species

##    code                        common_name                    scientific_name
## 1   DIW                 Gibson's albatross      Diomedea antipodensis gibsoni
## 2   DQS               Antipodean albatross Diomedea antipodensis antipodensis
## 3   DIX                Wandering albatross                   Diomedea exulans
## 4   DBN                  Tristan albatross                 Diomedea dabbenena
## 5   DAM                Amsterdam albatross            Diomedea amsterdamensis
## 6   DIP           Southern royal albatross                Diomedea epomophora
## 7   DIQ           Northern royal albatross                  Diomedea sanfordi
## 8   DCR    Atlantic yellow-nosed albatross        Thalassarche chlororhynchos
## 9   TQH      Indian yellow-nosed albatross               Thalassarche carteri
## 10  DIM             Black-browed albatross           Thalassarche melanophris
## 11  TQW    Campbell black-browed albatross              Thalassarche impavida
## 12  DCU                      Shy albatross                 Thalassarche cauta
## 13  TWD New Zealand white-capped albatross          Thalassarche cauta steadi
## 14  DKS                 Salvin's albatross               Thalassarche salvini
## 15  DER           Chatham Island albatross               Thalassarche eremita
## 16  DIC              Grey-headed albatross           Thalassarche chrysostoma
## 17  DSB        Southern Buller's albatross       Thalassarche bulleri bulleri
## 18  DNB        Northern Buller's albatross        Thalassarche bulleri platei
## 19  PHU                    Sooty albatross                   Phoebetria fusca
## 20  PHE      Light-mantled sooty albatross              Phoebetria palpebrata
## 21  PCI                        Grey petrel                Procellaria cinerea
## 22  PRK                       Black petrel             Procellaria parkinsoni
## 23  PCW                    Westland petrel            Procellaria westlandica
## 24  PRO               White-chinned petrel         Procellaria aequinoctialis
## 25  PCN                  Spectacled petrel          Procellaria conspicillata
##           genus         family code_resolution id_species
## 1      Diomedea    Diomedeidae         species          1
## 2      Diomedea    Diomedeidae         species          2
## 3      Diomedea    Diomedeidae         species          3
## 4      Diomedea    Diomedeidae         species          4
## 5      Diomedea    Diomedeidae         species          5
## 6      Diomedea    Diomedeidae         species          6
## 7      Diomedea    Diomedeidae         species          7
## 8  Thalassarche    Diomedeidae         species          8
## 9  Thalassarche    Diomedeidae         species          9
## 10 Thalassarche    Diomedeidae         species         10
## 11 Thalassarche    Diomedeidae         species         11
## 12 Thalassarche    Diomedeidae         species         12
## 13 Thalassarche    Diomedeidae         species         13
## 14 Thalassarche    Diomedeidae         species         14
## 15 Thalassarche    Diomedeidae         species         15
## 16 Thalassarche    Diomedeidae         species         16
## 17 Thalassarche    Diomedeidae         species         17
## 18 Thalassarche    Diomedeidae         species         18
## 19   Phoebetria    Diomedeidae         species         19
## 20   Phoebetria    Diomedeidae         species         20
## 21  Procellaria Procellariidae         species         21
## 22  Procellaria Procellariidae         species         22
## 23  Procellaria Procellariidae         species         23
## 24  Procellaria Procellariidae         species         24
## 25  Procellaria Procellariidae         species         25

Assign biological inputs to the data object. Assignment functions are provided for each data type, and automatically select the required species:

n_breeding_pairs(sefra_data) <- N_BP       
adult_survival(sefra_data)   <- S_opt      
p_breeding(sefra_data)       <- P_B        
age_breeding(sefra_data)     <- A_curr     
p_nest(sefra_data)           <- p_nest
p_southern(sefra_data)       <- p_southern

We can view the input data using, for example:

N_BP %>% select(id_species, code, distribution, par1, par2)

##    id_species code distribution       par1      par2
## 1           1  DIW   log-normal    5400.00     0.050
## 2           2  DQS   log-normal    3255.00     0.050
## 3           3  DIX   log-normal   10072.00     0.050
## 4           4  DBN      weibull       9.25  1710.000
## 5           5  DAM   log-normal      60.00     0.100
## 6           6  DIP   log-normal    6347.00     0.140
## 7           7  DIQ   log-normal    4261.00     0.110
## 8           8  DCR   log-normal   26800.00     0.100
## 9           9  TQH   log-normal   33988.00     0.100
## 10         10  DIM   log-normal  681000.00     0.050
## 11         11  TQW   log-normal   24338.00     0.050
## 12         12  DCU   log-normal   15335.00     0.100
## 13         13  TWD   log-normal   71400.00     0.110
## 14         14  DKS   log-normal   58563.00     0.100
## 15         15  DER   log-normal    5294.00     0.010
## 16         16  DIC   log-normal   80633.00     0.050
## 17         17  DSB   log-normal   14320.00     0.050
## 18         18  DNB   log-normal   19354.00     0.050
## 19         19  PHU      weibull      23.20 13660.000
## 20         20  PHE   log-normal   20927.00     0.100
## 21         21  PCI   log-normal  105617.00     0.150
## 22         22  PRK   log-normal    5456.00     0.057
## 23         23  PCW   log-normal    6223.00     0.061
## 24         24  PRO   log-normal 1317300.00     0.100
## 25         25  PCN   log-normal   42000.00     0.096

n_breeding_pairs(sefra_data)

##     distribution       par1      par2
## DIW   log-normal    5400.00     0.050
## DQS   log-normal    3255.00     0.050
## DIX   log-normal   10072.00     0.050
## DBN      weibull       9.25  1710.000
## DAM   log-normal      60.00     0.100
## DIP   log-normal    6347.00     0.140
## DIQ   log-normal    4261.00     0.110
## DCR   log-normal   26800.00     0.100
## TQH   log-normal   33988.00     0.100
## DIM   log-normal  681000.00     0.050
## TQW   log-normal   24338.00     0.050
## DCU   log-normal   15335.00     0.100
## TWD   log-normal   71400.00     0.110
## DKS   log-normal   58563.00     0.100
## DER   log-normal    5294.00     0.010
## DIC   log-normal   80633.00     0.050
## DSB   log-normal   14320.00     0.050
## DNB   log-normal   19354.00     0.050
## PHU      weibull      23.20 13660.000
## PHE   log-normal   20927.00     0.100
## PCI   log-normal  105617.00     0.150
## PRK   log-normal    5456.00     0.057
## PCW   log-normal    6223.00     0.061
## PRO   log-normal 1317300.00     0.100
## PCN   log-normal   42000.00     0.096

Assign fishery groups:

fishery_groups(sefra_data) <- captures_o
fishery_groups(sefra_data)

## 4 fishery groups

##   fishery_group id_fishery_group
## 1     NZL (DOM)                1
## 2      NZL (JV)                2
## 3           JPN                3
## 4           TWN                4

View capture codes:

capture_codes(sefra_data)

## 21 capture codes:

## (empty captures data frame)

##    code id_code resolution id_resolution
## 1   DKS      14    species             1
## 2   DER      15    species             1
## 3   DIC      16    species             1
## 4   PHU      19    species             1
## 5   PHE      20    species             1
## 6   PCI      21    species             1
## 7   PCN      25    species             1
## 8   DRA      27    complex             2
## 9   DYN      28    complex             2
## 10  DST      29    complex             2
## 11  DBB      30    complex             2
## 12  DIB      31    complex             2
## 13  DWC      32    complex             2
## 14  PRZ      33    complex             2
## 15  DIZ      34      genus             3
## 16  THZ      35      genus             3
## 17  PHZ      36      genus             3
## 18  PTZ      37      genus             3
## 19  ALZ      38     family             4
## 20  PRX      39     family             4
## 21  BLZ      40      class             5

Assign species groups:

species_groups(sefra_data) <- overlap_o
species_groups(sefra_data)

## 5 species groups

##         species_group id_species_group
## 1 Wandering albatross                1
## 2     Royal albatross                2
## 3           Mollymawk                3
## 4     Sooty albatross                4
## 5       Medium petrel                5

species_groups(sefra_data, print_species = TRUE)

## 5 species groups

## (with 25 species)

##    species       species_group id_species id_species_group
## 1      DIW Wandering albatross          1                1
## 2      DQS Wandering albatross          2                1
## 3      DIX Wandering albatross          3                1
## 4      DBN Wandering albatross          4                1
## 5      DAM Wandering albatross          5                1
## 6      DIP     Royal albatross          6                2
## 7      DIQ     Royal albatross          7                2
## 8      DCR           Mollymawk          8                3
## 9      TQH           Mollymawk          9                3
## 10     DIM           Mollymawk         10                3
## 11     TQW           Mollymawk         11                3
## 12     DCU           Mollymawk         12                3
## 13     TWD           Mollymawk         13                3
## 14     DKS           Mollymawk         14                3
## 15     DER           Mollymawk         15                3
## 16     DIC           Mollymawk         16                3
## 17     DSB           Mollymawk         17                3
## 18     DNB           Mollymawk         18                3
## 19     PHU     Sooty albatross         19                4
## 20     PHE     Sooty albatross         20                4
## 21     PCI       Medium petrel         21                5
## 22     PRK       Medium petrel         22                5
## 23     PCW       Medium petrel         23                5
## 24     PRO       Medium petrel         24                5
## 25     PCN       Medium petrel         25                5

Assign cryptic capture distributions:

cryptic_capture(sefra_data) <- cc_longline

## 4 fishery groups

## 5 species groups

cryptic_capture(sefra_data)

##    distribution fishery_group       species_group  par value
## 1    log-normal     NZL (DOM) Wandering albatross par1 1.420
## 2    log-normal      NZL (JV) Wandering albatross par1 1.420
## 3    log-normal           JPN Wandering albatross par1 1.420
## 4    log-normal           TWN Wandering albatross par1 1.420
## 5    log-normal     NZL (DOM)     Royal albatross par1 1.420
## 6    log-normal      NZL (JV)     Royal albatross par1 1.420
## 7    log-normal           JPN     Royal albatross par1 1.420
## 8    log-normal           TWN     Royal albatross par1 1.420
## 9    log-normal     NZL (DOM)           Mollymawk par1 1.420
## 10   log-normal      NZL (JV)           Mollymawk par1 1.420
## 11   log-normal           JPN           Mollymawk par1 1.420
## 12   log-normal           TWN           Mollymawk par1 1.420
## 13   log-normal     NZL (DOM)     Sooty albatross par1 1.420
## 14   log-normal      NZL (JV)     Sooty albatross par1 1.420
## 15   log-normal           JPN     Sooty albatross par1 1.420
## 16   log-normal           TWN     Sooty albatross par1 1.420
## 17   log-normal     NZL (DOM)       Medium petrel par1 1.420
## 18   log-normal      NZL (JV)       Medium petrel par1 1.420
## 19   log-normal           JPN       Medium petrel par1 1.420
## 20   log-normal           TWN       Medium petrel par1 1.420
## 21   log-normal     NZL (DOM) Wandering albatross par2 0.186
## 22   log-normal      NZL (JV) Wandering albatross par2 0.186
## 23   log-normal           JPN Wandering albatross par2 0.186
## 24   log-normal           TWN Wandering albatross par2 0.186
## 25   log-normal     NZL (DOM)     Royal albatross par2 0.186
## 26   log-normal      NZL (JV)     Royal albatross par2 0.186
## 27   log-normal           JPN     Royal albatross par2 0.186
## 28   log-normal           TWN     Royal albatross par2 0.186
## 29   log-normal     NZL (DOM)           Mollymawk par2 0.186
## 30   log-normal      NZL (JV)           Mollymawk par2 0.186
## 31   log-normal           JPN           Mollymawk par2 0.186
## 32   log-normal           TWN           Mollymawk par2 0.186
## 33   log-normal     NZL (DOM)     Sooty albatross par2 0.186
## 34   log-normal      NZL (JV)     Sooty albatross par2 0.186
## 35   log-normal           JPN     Sooty albatross par2 0.186
## 36   log-normal           TWN     Sooty albatross par2 0.186
## 37   log-normal     NZL (DOM)       Medium petrel par2 0.186
## 38   log-normal      NZL (JV)       Medium petrel par2 0.186
## 39   log-normal           JPN       Medium petrel par2 0.186
## 40   log-normal           TWN       Medium petrel par2 0.186

Observed capture and overlap data for estimation of catchabilities:

sefra_data <- sefra_data %>% data_prep(overlap_o, captures_o)

## Prepared observer captures and overlap data

Overlap data for full diagnostics outputs:

sefra_data <- sefra_data %>% data_prep(overlap_o)

## Prepared commercial (total) overlap data

Force taxonomic resolution of species and complex capture codes to be equivalent (i.e., so they share the same pi term):

sefra_data$n_resolutions <- 4
sefra_data$id_resolutions <- sefra_data$id_resolutions[!sefra_data$id_resolutions %in% 2L]
sefra_data$resolutions[sefra_data$resolutions %in% c("species", "complex")] <- "species"
sefra_data$code_resolution[sefra_data$code_resolution %in% 2L] <- 1L

Force all species to share the same pi vector:

sefra_data$n_pi <- 1L
sefra_data$id_pi <- rep(1, times = length(sefra_data$id_species))

save(sefra_data, file = file.path(output_path, "sefra_data.rda"))

Model run

Specify number of iterations:

n_iter <- 1E3

First, source and compile the model code. Because there are no commercial data included, this model will not generate a full suite of diagnostic outputs (only catchability and catches):

mdl <- sefra(sefra_data, likelihood = "cumulative")

## SEFRA-seabird model v2.4.0

Create initial values:

mdl_ini <- initial_values(sefra_data, mdl)

save(mdl_ini,  file = file.path(output_path, "mdl_ini.rda"))

Run model and check MCMC trace plots:

mdl_fit <- rstan::sampling(
  mdl, data = sefra_data, verbose = TRUE, init = function() mdl_ini, chains = 2, iter = n_iter,
  control = list(adapt_delta = 0.95, stepsize = 0.1, max_treedepth = 100), cores = 2
)

save(mdl_fit,  file = file.path(output_path, "mdl_fit.rda"))

Prelims for processing fitted model

Make folders and subfolders for outputs:

make_folder(file.path(output_path, "maps"), clean = TRUE)
make_folder(file.path(output_path, "figures"), clean = TRUE)
make_folder(file.path(output_path, "tables"), clean = TRUE)

Create and save sefraOutputs object that holds summarised model outputs:

out <- sefraOutputs(sefra_data, mdl_fit)

## constructed 'sefraOutputs' object

save(out, file = file.path(output_path, "out.rda"))

Model diagnostics

## Info useful for writing script
x <- rstan::extract(mdl_fit, pars = "lp__", permute = FALSE, include = FALSE)
x <- dimnames(x)[[3]]
x <- gsub("\\[.*$", "", x)
x <- unique(x)
quantity_names <- x

Discrete colour palette that is colour-blind friendly:

## Colour blind friendly palette by Okabe & Ito
okabe_cols <- c("#56B4E9", "#E69F00", "#009E73", "#F0E442", "#CC79A7", "#0072B2", "#D55E00", "#000000")
names(okabe_cols) <- c("light blue", "orange", "green", "yellow", "pink", "dark blue", "red", "black")

## function to get summary statistics (mean and 95% CI)
summary_fn <- function(x, var) {
    summarise(x, mean = mean({{ var }}), mq = median({{ var }}),
              lq = as.numeric(quantile({{ var }}, 0.025)), uq = as.numeric(quantile({{ var }}, 0.975)))
}

## summary function that gets quantiles to generate boxplot of relative mortalities
summary_fn_boxplot <- function(x, var) {
  summarise(x, lq = as.numeric(quantile({{ var }}, probs = 0.025)), lmq = as.numeric(quantile({{ var }}, probs = 0.125)),
             mq = as.numeric(quantile({{ var }}, probs = 0.5)),
             umq = as.numeric(quantile({{ var }}, probs = 0.875)), uq = as.numeric(quantile({{ var }}, probs = 0.975)))
}

Set colours for each species group:

cols_species_group <- okabe_cols[1:length(sefra_data[['species_groups']])]
names(cols_species_group) <- sefra_data[['species_groups']]

Specify order of groups for plots:

sp_group_levels <- species_groups(sefra_data) %>% filter(., !is.na(species_group)) %>% select(., id_species_group, species_group) %>% distinct(.) %>% with(., species_group)

## 5 species groups

fg_group_levels <- c("NZL (DOM)", "NZL (JV)", "JPN", "TWN")

Specify colours for taxonomic resolutions (e.g., used in plots of pi parameters):

taxo_resolution <- sp_groups %>% select(., taxonomic_resolution) %>% distinct(.) %>% with(., taxonomic_resolution)
taxo_resolution[taxo_resolution %in% "bird"] <- "class" 
cols_taxo_resolution <- okabe_cols[1:length(taxo_resolution)]
names(cols_taxo_resolution) <- taxo_resolution

Create variable for plots with taxonomic grouping (genus, family or class):

sp_groups <- sp_groups %>%
  mutate(taxo_group = case_when(
    taxonomic_resolution %in% c("species", "complex", "genus") ~ genus,
    taxonomic_resolution %in% "family" ~ family,
    .default = "Unspecified"
    ))

taxo_group_levels <- sp_groups %>% select(., taxo_group) %>% distinct(.) %>% with(., taxo_group)
  
cols_taxo_group <- okabe_cols[1:length(taxo_group_levels)]
names(cols_taxo_group) <- taxo_group_levels

Trace plots of model parameters (and estimated quantities)

Trace plots can be used to visualise convergence. For example, we can plot trace diagnostics for the biological parameters using:

trace_plot(mdl_fit, pars = c("p_breeding"), labels = paste0('P[', sefra_data$species, ']'))

trace_plot(mdl_fit, pars = c("age_breeding"), labels = paste0('A[', sefra_data$species, ']'))

trace_plot(mdl_fit, pars = c("adult_survival"), labels = paste0('S[', sefra_data$species, ']'))

trace_plot(mdl_fit, pars = c("n_breeding_pairs"), labels = paste0('N[', sefra_data$species, ']')) + facet_wrap(~parameter, scales = "free_y")

Trace plots of catchability parameters:

#trace_plot(mdl_fit, pars = c("alpha_intercept"))
#trace_plot(mdl_fit, pars = c("alpha_fishery_group"))
#trace_plot(mdl_fit, pars = c("alpha_species_group"))
trace_plot(mdl_fit, pars = c("q")) + facet_wrap(~parameter, scales = "free_y") + theme(axis.text.y = element_blank())

Trace plots of probability of live capture:

#trace_plot(mdl_fit, pars = c("beta_intercept"))
#trace_plot(mdl_fit, pars = c("beta_fishery_group"))
#trace_plot(mdl_fit, pars = c("beta_species_group"))
trace_plot(mdl_fit, pars = c("p_live")) + facet_wrap(~parameter, scales = "free_y") + theme(axis.text.y = element_blank())

Trace plots of probability of observed capture:

par_name <- "p_obs_logit"
#if("p_obs_cde" %in% quantity_names) par_name <- "p_obs_cde"
trace_plot(mdl_fit, pars = par_name) + facet_wrap(~parameter, scales = "free_y") + theme(axis.text.y = element_blank())

Updates to biological input distributions

To visualise prior updates to the biological input distributions, we use, for example:

gg <- plot_prior_update(mdl_fit, par = "n_breeding_pairs", labels = sefra_data$species)
gg

ggsave(gg, file = file.path(output_path, "figures", "prior_updates_N_BP.png"), width = 8, height = 12)
ggsave(gg, file = file.path(output_path, "figures", "prior_updates_N_BP.pdf"), width = 8, height = 12)

gg <- plot_prior_update(mdl_fit, par = "n_breeding_pairs", labels = sefra_data$species) + scale_y_log10()
gg

ggsave(gg, file = file.path(output_path, "figures", "prior_updates_N_BP_log10.png"), width = 8, height = 12)
ggsave(gg, file = file.path(output_path, "figures", "prior_updates_N_BP_log10.pdf"), width = 8, height = 12)

gg <- plot_prior_update(mdl_fit, par = "p_breeding", labels = sefra_data$species)
gg

ggsave(gg, file = file.path(output_path, "figures", "prior_updates_P_B.png"), width = 8, height = 12)
ggsave(gg, file = file.path(output_path, "figures", "prior_updates_P_B.pdf"), width = 8, height = 12)

Plots of rhat and Neff

Diagnostic plots also exist for Rhat and Neff. For example:

plot_rhat(mdl_fit, pars = c("alpha_fishery_group", "alpha_species_group", "beta_fishery_group", "beta_species_group"))

#plot_neff(mdl_fit, pars = c("alpha_fishery_group", "alpha_species_group", "beta_fishery_group", "beta_species_group"))

Trace plots of model parameters

Trace plot of summary statistics:

x <- rstan::extract(mdl_fit, pars = paste0("traces[", 1:4, "]"), permute = FALSE)

dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), chain = paste0("chain_", 1:dim(x)[2]),
                    parameter = c("Catchability", "Prob. live", "Vuln. breeders", "lambda"))
x <- reshape2::melt(x)

gg <- ggplot(x) +
  geom_line(aes(iteration, value, col = chain)) +
  scale_color_viridis_d(option = "D", begin = 0.2, end = 0.8) +
  facet_wrap(vars(parameter), ncol = 1, scales = "free_y") +
  theme_bw(base_size = 14) +
  theme(axis.text.y = element_blank(), axis.title.y = element_blank()) +
  guides(col = "none") + xlab("Posterior sample")
gg

ggsave(gg, file = file.path(output_path, "figures", "trace_plot.png"), width = 6, height = 10)
ggsave(gg, file = file.path(output_path, "figures", "trace_plot.pdf"), width = 6, height = 10)

Observed vs predicted captures

Perhaps the most useful diagnostic is provided from examination of fits to the capture data, and a specific function exists for this purpose:

captures(sefra_data, mdl_fit)

## $captures
## # A tibble: 84 × 7
##    fishery_group code  captures     hat   med low_ci upp_ci
##    <fct>         <chr>    <int>   <dbl> <dbl>  <dbl>  <dbl>
##  1 NZL (DOM)     DKS          1  12.2      12   6        20
##  2 NZL (DOM)     DER          0   0.537     0   0         2
##  3 NZL (DOM)     DIC          1   0.059     0   0         1
##  4 NZL (DOM)     PHU          0   0         0   0         0
##  5 NZL (DOM)     PHE          0   0         0   0         0
##  6 NZL (DOM)     PCI          4  17.4      17   9        27
##  7 NZL (DOM)     PCN          0   0         0   0         0
##  8 NZL (DOM)     DRA         11  10.4      10   3.98     20
##  9 NZL (DOM)     DYN          0   0         0   0         0
## 10 NZL (DOM)     DST        151 100.      100  78       124
## # ℹ 74 more rows
## 
## $captures_inclusive
## # A tibble: 84 × 7
##    fishery_group code  captures    hat   med low_ci upp_ci
##    <chr>         <chr>    <int>  <dbl> <dbl>  <dbl>  <dbl>
##  1 JPN           DKS          0  12.7     13    6      21 
##  2 JPN           DER          3   2.43     2    0       6 
##  3 JPN           DIC        840 735.     734  665     807 
##  4 JPN           PHU        134 146.     146  117     177 
##  5 JPN           PHE         95  74.3     74   55.0    96 
##  6 JPN           PCI        152 107.     106   81     135.
##  7 JPN           PCN         44  24.1     24   13      37 
##  8 JPN           DRA         21  21.1     21   10      35 
##  9 JPN           DYN        205 235.     235  198.    275 
## 10 JPN           DST        796 771.     770  705.    845 
## # ℹ 74 more rows

‘Cumulative’ captures for a given capture code include all captures from other capture codes recorded at finer taxonomic resolutions. For example, cumulative captures for ‘PTZ’ (Procellaria genus) includes captures identified as an unspecified Procellaria species, as well as all captures identified at a species or complex level within the Procellaria genus (e.g., black petrels).

‘Empirical’ captures for a given capture code only include captures identified as that specific capture code. For example, empirical captures for ‘PTZ’ (Procellaria genus) only includes captures with an assigned capture code of ‘PTZ’ - i.e., identified as an unspecified Procellaria species.

The model can be fitted to either cumulative or empirical captures (specified by the likelihood argument to the sefra function). Comparisons of captures can be informative for both cumulative and empirical captures, regardless of which type of captures the model was fitted to. However, model fits should only be assessed using the type of capture that the model was fitted to.

Comparisons of total empirical and cumulative captures (per capture code ), and total estimated captures (per species):

est_catch_spp <- posterior(mdl_fit, pars = c('E_MDL_captures'))
est_catch_code_emp <- posterior(mdl_fit, pars = c('E_EMP_captures'))
est_catch_code_csm <- posterior(mdl_fit, pars = c('E_CSM_captures'))
est_pi <- posterior(mdl_fit, pars = c('p_obs'))

## these should be equal:
sum(est_catch_spp)

## [1] 7681132

sum(est_catch_code_csm[, ,  nrow(mdl_codes)])

## [1] 7681132

sum(est_catch_code_emp)

## [1] 7681132

Observed vs predicted cumulative captures by capture group

Get data-frame with observed and predicted cumulative captures by code:

x <- captures(sefra_data, mdl_fit)[["captures_inclusive"]]
y <- select(filter(sp_groups, !(code %in% codes_exclude)), code, common_name, species_group, taxo_group)

x <- x %>% left_join(y)

## Joining with `by = join_by(code)`

x <- x %>% mutate(
  species_group = factor(species_group, levels = sp_group_levels),
  taxo_group = factor(taxo_group, levels = taxo_group_levels),
  fishery_group = factor(fishery_group, levels = fg_group_levels)
)

x <- x %>% mutate(., across(c(captures, hat, med, low_ci, upp_ci), function(a) a / sefra_data[["n_year_k"]]))

plt_dat_point <- x
plt_dat_bar   <- x %>% tidyr::pivot_longer(captures:hat) %>% mutate(low_ci = case_when(name == "captures" ~ NA, TRUE ~ low_ci), upp_ci = case_when(name == "captures" ~ NA, TRUE ~ upp_ci), med = case_when(name == "captures" ~ NA, TRUE ~ med))

Barchart of observed vs predicted cumulative captures by capture code and fishery group:

## Get colour palette
plt_cols <- scales::viridis_pal(alpha = 0.7, option = "D", begin = 0.1, end = 0.9)(2)
names(plt_cols) <- c("captures", "hat")

## reorder capture codes for plot using ordering in sp_groups
plt_dat_bar <- plt_dat_bar %>% mutate(., code = factor(code, levels = sp_groups$code))

gg <- plt_dat_bar %>% ggplot(.) + 
  geom_col(aes(x = code, y = value, col = name, fill = name), position = "dodge") +
  #geom_errorbar(aes(x = code, ymin = low_ci, ymax = upp_ci, col = name), position = "dodge", alpha = 1) +
  scale_fill_manual("", values = plt_cols, aesthetics = c("fill", "color"), labels = c("Observed", "Posterior prediction")) +
  xlab("Capture code") + ylab("Average annual captures") + 
  facet_wrap(vars(fishery_group), ncol = 1, scales = "free_y") +
  theme_bw() + theme(legend.position = "bottom")
gg

ggsave(gg, file = file.path(output_path, "figures", "fit_cumulative_captures_barchart_capture_code.png"), width = 8)

## Saving 8 x 7 in image

ggsave(gg, file = file.path(output_path, "figures", "fit_cumulative_captures_barchart_capture_code.pdf"), width = 8)

## Saving 8 x 7 in image

Plots of observed vs. predicted cumulative captures by species group and fishery group:

gg <- plt_dat_point %>%
  #filter(., obs > 0 & mq > 0) %>%
  ggplot() +
  geom_pointrange(aes(x = captures, y = med, ymin = low_ci, ymax = upp_ci, colour = taxo_group), alpha = 0.5, size = 1.0) +
  scale_colour_manual("", values = cols_taxo_group) +
  geom_abline(slope = 1) +
  facet_wrap(fishery_group ~ .) +
  theme_bw(base_size = 18) +
  theme(legend.position = "bottom", legend.title = element_text(size = 14), legend.text = element_text(size = 12)) +
  labs(x = "Observed number of annual captures", y = "Predicted number of annual observed captures") +
  theme(aspect.ratio = 1) + scale_y_log10() + scale_x_log10()
gg

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

## Warning in scale_x_log10(): log-10 transformation introduced infinite values.

ggsave(gg, file = file.path(output_path, "figures", "fit_cumulative_captures_point_capture_code.png"), width = 10)

## Saving 10 x 10 in image

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

## Warning in scale_x_log10(): log-10 transformation introduced infinite values.

ggsave(gg, file = file.path(output_path, "figures", "fit_cumulative_captures_point_capture_code.pdf"), width = 10)

## Saving 10 x 10 in image

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

## Warning in scale_x_log10(): log-10 transformation introduced infinite values.

Observed vs predicted empirical captures by capture group

Get data-frame with observed and predicted empirical captures by code:

x <- captures(sefra_data, mdl_fit)[["captures"]]
y <- select(filter(sp_groups, !(code %in% codes_exclude)), code, common_name, species_group, taxo_group)

x <- x %>% left_join(y)

## Joining with `by = join_by(code)`

x <- x %>% mutate(
  species_group = factor(species_group, levels = sp_group_levels),
  taxo_group = factor(taxo_group, levels = taxo_group_levels),
  fishery_group = factor(fishery_group, levels = fg_group_levels)
)

x <- x %>% mutate(., across(c(captures, hat, med, low_ci, upp_ci), function(a) a / sefra_data[["n_year_k"]]))

plt_dat_point <- x
plt_dat_bar   <- x %>% tidyr::pivot_longer(captures:hat) %>% mutate(low_ci = case_when(name == "captures" ~ NA, TRUE ~ low_ci), upp_ci = case_when(name == "captures" ~ NA, TRUE ~ upp_ci), med = case_when(name == "captures" ~ NA, TRUE ~ med))

Barchart of observed vs predicted empirical captures by capture code and fishery group:

## Get colour palette
plt_cols <- scales::viridis_pal(alpha = 0.7, option = "D", begin = 0.1, end = 0.9)(2)
names(plt_cols) <- c("captures", "hat")

## reorder capture codes for plot using ordering in sp_groups
plt_dat_bar <- plt_dat_bar %>% mutate(., code = factor(code, levels = sp_groups$code))

gg <- plt_dat_bar %>% ggplot(.) + 
  geom_col(aes(x = code, y = value, col = name, fill = name), position = "dodge") +
  #geom_errorbar(aes(x = code, ymin = low_ci, ymax = upp_ci, col = name), position = "dodge", alpha = 1) +
  scale_fill_manual("", values = plt_cols, aesthetics = c("fill", "color"), labels = c("Observed", "Posterior prediction")) +
  xlab("Capture code") + ylab("Average annual captures") + 
  facet_wrap(vars(fishery_group), ncol = 1, scales = "free_y") +
  theme_bw() + theme(legend.position = "bottom")
gg

ggsave(gg, file = file.path(output_path, "figures", "fit_empirical_captures_barchart_capture_code.png"), width = 8)

## Saving 8 x 7 in image

ggsave(gg, file = file.path(output_path, "figures", "fit_empirical_captures_barchart_capture_code.pdf"), width = 8)

## Saving 8 x 7 in image

Plots of observed vs. predicted empirical captures by species group and fishery group:

gg <- plt_dat_point %>%
  #filter(., obs > 0 & mq > 0) %>%
  ggplot() +
  geom_pointrange(aes(x = captures, y = med, ymin = low_ci, ymax = upp_ci, colour = taxo_group), alpha = 0.5, size = 1.0) +
  scale_colour_manual("", values = cols_taxo_group) +
  geom_abline(slope = 1) +
  facet_wrap(fishery_group ~ .) +
  theme_bw(base_size = 18) +
  theme(legend.position = "bottom", legend.title = element_text(size = 14), legend.text = element_text(size = 12)) +
  labs(x = "Observed number of annual captures", y = "Predicted number of annual observed captures") +
  theme(aspect.ratio = 1) + scale_y_log10() + scale_x_log10()
gg

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

## Warning in scale_x_log10(): log-10 transformation introduced infinite values.

ggsave(gg, file = file.path(output_path, "figures", "fit_empirical_captures_point_capture_code.png"), width = 10)

## Saving 10 x 10 in image

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

## Warning in scale_x_log10(): log-10 transformation introduced infinite values.

ggsave(gg, file = file.path(output_path, "figures", "fit_empirical_captures_point_capture_code.pdf"), width = 10)

## Saving 10 x 10 in image

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.

## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.

## Warning in scale_x_log10(): log-10 transformation introduced infinite values.

Tables of estimated quantities

Estimated captures

Estimated captures by species group:

x <- out[["captures_inclusive"]]

## Add id_code and species group information
x <- x %>% mutate(., id_code = as.numeric(code))
x <- x %>% mutate(., id_species_group = sefra_data[["species_to_species_group"]][id_code])
x <- x %>% left_join(., species_groups(sefra_data))

## 5 species groups

## Joining with `by = join_by(id_species_group)`

## Prepare observed captures
x <- x %>% #filter(., code %in% sefra_data$species) %>%
  filter(., name %in% c("emp", "sim")) %>%
  group_by(., fishery_group, species_group, name, iter) %>%
  summarise(., value = sum(value)) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group', 'species_group', 'name'.
## You can override using the `.groups` argument.

x <- x %>% group_by(., fishery_group, species_group, name) %>%
  summarise(., mean = mean(value), mq = median(value),
            lq = as.numeric(quantile(value, probs = 0.025)),
            uq = as.numeric(quantile(value, probs = 0.975))) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group', 'species_group'. You can
## override using the `.groups` argument.

## reorder species groups for plot
x$species_group <- factor(x$species_group, levels = sefra_data[['species_groups']])

## Restructure to have observed captures as a variable
tab_dat <- x %>%
  filter(., name %in% "emp") %>%
  select(., fishery_group, species_group, obs = mean)

tab_dat <- x %>% filter(., !name %in% "emp") %>%
  left_join(tab_dat, ., by = c("fishery_group", "species_group"))

## Combine estimate and 95% CI
tab_dat <- tab_dat %>% mutate(., est = paste0(round(mean, 1), " (", round(lq, 1), "-", round(uq, 1), ")"))

tab_dat <- tab_dat %>% tidyr::pivot_wider(id_cols = species_group, values_from = c('obs', 'est'), names_from = fishery_group) %>%
  arrange(., species_group)

## Reorder variables
tab_dat <- tab_dat %>% select('species_group', contains('DOM'), contains('JV'), contains('JPN'), contains('TWN'))

tab <- kable(tab_dat, format = 'latex',
             col.names = c('', rep(c('Observed', 'Estimated'), times = sefra_data$n_fishery_groups)),
             row.names = FALSE,
             caption = "Fit to (cumulative) captures by species and fisheries group.",
             label = 'fit_captures_by_species_group_fishery_group',
             align = c('l', rep('r', times = ncol(tab_dat) - 1)),
             linesep = "", booktabs = TRUE) %>%
    add_header_above(c("Species group", "NZL (DOM)" = 2, "NZL (JV)" = 2, "JPN" = 2, "TWN" = 2), align = c('l', 'c', 'c', 'c', 'c')) %>%
    kable_styling(font_size = 10, latex_options = "hold_position") %>%
  sub("\\\\cmidrule\\(l\\{3pt\\}r\\{3pt\\}\\)\\{1-1\\}", "", .) %>%
  sub("\\\\toprule", "", .) %>%
  sub("\\\\bottomrule", "", .)

save_kable(tab, file = file.path(output_path, "tables", "fit_captures_per_species_group_fishery_group.tex"))
write.csv(tab_dat, file = file.path(output_path, "tables", "fit_captures_per_species_group_fishery_group.csv"))

Fit to (cumulative) captures by species and fisheries group.
Species group	NZL (DOM)		NZL (JV)		JPN		TWN
	Observed	Estimated	Observed	Estimated	Observed	Estimated	Observed	Estimated
Wandering albatross	2	12.7 (6-21)	0	1 (0-3)	1072	971 (894-1050)	92	175.7 (147-207)
Royal albatross	4	17.3 (9-27)	0	0.2 (0-2)	196	131.1 (104-162)	56	55 (40-72)
Mollymawk	826	748.2 (682-813)	152	130 (103-158)	8269	8420.2 (8174-8669)	1277	1182.8 (1090-1284)
Sooty albatross	417	428.5 (383-479)	76	81.5 (63-104)	5904	5867.2 (5667.9-6070.1)	909	886.2 (810-968)
Medium petrel	417	459.6 (408-518)	76	93.4 (71-117)	6127	6175.9 (5975-6387)	917	955.7 (876-1036)

Estimated captures by code:

x <- out[["captures_inclusive"]]

## Add id_code and species group information
x <- x %>% mutate(., id_code = as.numeric(code))
x <- x %>% mutate(., id_species_group = sefra_data[["species_to_species_group"]][id_code])
x <- x %>% left_join(., species_groups(sefra_data))

## 5 species groups

## Joining with `by = join_by(id_species_group)`

## Prepare observed captures
x <- x %>% #filter(., code %in% sefra_data$species) %>%
  filter(., name %in% c("emp", "sim")) %>%
  group_by(., fishery_group, species_group, code, name, iter) %>%
  summarise(., value = sum(value)) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group', 'species_group', 'code',
## 'name'. You can override using the `.groups` argument.

x <- x %>% group_by(., fishery_group, species_group, code, name) %>%
  summarise(., mean = mean(value), mq = median(value),
            lq = as.numeric(quantile(value, probs = 0.025)),
            uq = as.numeric(quantile(value, probs = 0.975))) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group', 'species_group', 'code'.
## You can override using the `.groups` argument.

## reorder species groups for plot
x$species_group <- factor(x$species_group, levels = sefra_data[['species_groups']])

## Restructure to have observed captures as a variable
tab_dat <- x %>%
  filter(., name %in% "emp") %>%
  select(., fishery_group, code, obs = mean)

tab_dat <- x %>% filter(., !name %in% "emp") %>%
  left_join(tab_dat, ., by = c("fishery_group", "code"))

## Combine estimate and 95% CI
tab_dat <- tab_dat %>% mutate(., est = paste0(round(mean, 1), " (", round(lq, 1), "-", round(uq, 1), ")"))

tab_dat <- tab_dat %>% tidyr::pivot_wider(id_cols = code, values_from = c('obs', 'est'), names_from = fishery_group) %>%
  arrange(., code)

## Add species names
tab_dat <- species_names(sefra_data) %>% select(., code, common_name) %>%
  left_join(tab_dat, ., by = "code")

## Species

# filter
tab_dat <- tab_dat %>% filter(!code %in% codes_exclude)

## Reorder variables
tab_dat <- tab_dat %>% select('code', contains('DOM'), contains('JV'), contains('JPN'), contains('TWN'))

tab <- kable(tab_dat, format = 'latex',
             col.names = c('', rep(c('Observed', 'Estimated'), times = sefra_data$n_fishery_groups)),
             row.names = FALSE,
             caption = "Fit to (cumulative) captures by capture code and fishery group.",
             label = 'fit_captures_by_code_fishery_group',
             align = c('l', rep('r', times = ncol(tab_dat) - 1)),
             linesep = "", booktabs = TRUE) %>%
    add_header_above(c("Code", "NZL (DOM)" = 2, "NZL (JV)" = 2, "JPN" = 2, "TWN" = 2), align = c('l', 'c', 'c', 'c', 'c')) %>%
    kable_styling(font_size = 10, latex_options = "hold_position") %>%
  sub("\\\\cmidrule\\(l\\{3pt\\}r\\{3pt\\}\\)\\{1-1\\}", "", .) %>%
  sub("\\\\toprule", "", .) %>%
  sub("\\\\bottomrule", "", .)

save_kable(tab, file = file.path(output_path, "tables", "fit_captures_per_code_fishery_group.tex"))
write.csv(tab_dat, file = file.path(output_path, "tables", "fit_captures_per_code_fishery_group.csv"))

Fit to (cumulative) captures by capture code and fisheries group.
Code	NZL (DOM)		NZL (JV)		JPN		TWN
	Observed	Estimated	Observed	Estimated	Observed	Estimated	Observed	Estimated
DKS	1	12.1 (5-20)	0	1 (0-3)	0	12.7 (6-21)	8	1.5 (0-4)
DER	0	0.5 (0-2)	0	0 (0-0)	3	2.4 (0-6)	0	0.1 (0-1)
DIC	1	0 (0-1)	0	0 (0-0)	840	735.2 (665-807)	17	110.1 (88-134)
PHU	0	0 (0-0)	0	0 (0-0)	134	146.3 (117-177)	61	49 (34-66)
PHE	0	0 (0-0)	0	0 (0-0)	95	74.3 (55-96)	6	15.1 (7-24)
PCI	4	17.3 (9-27)	0	0.2 (0-2)	152	107 (81-135)	3	17.3 (9-27)
PCN	0	0 (0-0)	0	0 (0-0)	44	24.1 (13-37)	53	37.7 (25-52)
DRA	11	10.2 (3-19)	0	0.3 (0-2)	21	21.1 (10-35)	4	3.6 (0-9)
DYN	0	0 (0-0)	0	0 (0-0)	205	234.9 (198-275)	148	78.4 (61-98)
DST	151	100.8 (82-124)	11	25 (15-36)	796	771.4 (705-845)	38	48 (34-62)
DBB	14	16.1 (9-25)	1	4.2 (1-9)	783	915.3 (838-1006)	113	45.5 (31-61)
DIB	125	133.6 (109-160)	62	26.7 (16-38)	789	730.7 (663-806)	4	35.4 (24-48)
DWC	27	24.6 (13-37)	0	0.4 (0-2)	342	341.8 (294-390)	54	52.4 (37-71)
PRZ	81	58.9 (43-77)	2	1.4 (0-5)	416	461.4 (412-513)	191	177 (145-208)
DIZ	40	37.6 (24-53)	0	0.9 (0-4)	389	390.8 (350-441)	58	61.7 (43-82)
THZ	292	283.9 (247-323)	74	69 (52-89)	3683	3674.5 (3519-3839.1)	328	352.8 (308-399)
PHZ	0	0 (0-0)	0	0 (0-0)	233	239.2 (204-280)	67	71.8 (52-93)
PTZ	85	82.4 (63-105)	2	2 (0-6)	612	639.3 (580-696)	272	256.4 (219-293)
ALZ	332	340.5 (299-386)	74	79.2 (61-102)	5127	5108.1 (4923-5297)	623	579.7 (520-642)
PRX	85	88 (68-110)	2	2.3 (0-6)	777	759.1 (693-829)	286	306.5 (266-349)
BLZ	417	459.6 (408-518)	76	93.4 (71-117)	6127	6175.9 (5975-6387)	917	955.7 (876-1036)

Estimated catchabilities

Estimated catchabilities (log-scale):

x <- rstan::extract(mdl_fit, pars = "q", permute = FALSE)
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), chain = paste0("chain_", 1:dim(x)[2]),
                    parameter = dimnames(x)[[3]])

x <- reshape2::melt(x)

x <- x %>% group_by(., parameter) %>%
  summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>%
  ungroup(.) %>%
  mutate(., across(c(mean, low_ci, high_ci), log10))

kable(x, digits = 2)

parameter	mean	low_ci	high_ci
q[1,1]	0.69	0.53	0.81
q[2,1]	-1.32	-2.18	-0.86
q[3,1]	1.18	1.11	1.25
q[4,1]	0.48	0.37	0.58
q[1,2]	0.52	0.22	0.74
q[2,2]	-1.12	-2.04	-0.59
q[3,2]	0.84	0.55	1.05
q[4,2]	0.37	-0.22	0.71
q[1,3]	0.34	0.27	0.40
q[2,3]	-0.67	-0.75	-0.60
q[3,3]	0.64	0.60	0.68
q[4,3]	-0.10	-0.16	-0.04
q[1,4]	0.77	-0.73	1.48
q[2,4]	-0.41	-1.81	0.34
q[3,4]	1.25	1.15	1.34
q[4,4]	0.71	0.59	0.83
q[1,5]	0.61	0.49	0.71
q[2,5]	-0.97	-1.40	-0.67
q[3,5]	0.50	0.41	0.59
q[4,5]	-0.18	-0.27	-0.09

Table of catchability estimates for reporting:

x <- rstan::extract(mdl_fit, pars = "q", permute = TRUE)[[1]]
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), fishery_group = sefra_data$fishery_groups,
                    species_group = sefra_data$species_groups)

x <- reshape2::melt(x) %>% group_by(fishery_group, species_group) %>% summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>% ungroup(.)

## `summarise()` has grouped output by 'fishery_group'. You can override using the
## `.groups` argument.

x <- x %>% tidyr::pivot_wider(values_from = c('mean', 'low_ci', 'high_ci'), names_from = fishery_group)

x <- x %>% select('species_group', contains('DOM'), contains('JV'), contains('JPN'), contains('TWN'))

tab <- kable(x, format = 'latex',
             digits = 2,
             col.names = c('', rep(c('mean', 'low', 'upp'), times = sefra_data$n_fishery_groups)),
             row.names = FALSE,
             caption = "Catchability",
             label = 'q_per_species_group_fishery_group',
             align = c('l', rep('r', times = ncol(x) - 1)),
             linesep = "", booktabs = TRUE) %>% kable_styling(font_size = 10, latex_options = "hold_position") 

if (grepl("combined", output_path)) {
    tab <- tab %>% add_header_above(c("Year", "NZL (DOM)" = 2, "NZL (JV)" = 2, "JPN" = 2, "TWN" = 2), align = c('l', rep('c', times = sefra_data$n_fishery_groups))) %>% sub("\\\\cmidrule\\(l\\{3pt\\}r\\{3pt\\}\\)\\{1-1\\}", "", .)
}
if (grepl("nzl", output_path)) {
    tab <- tab %>% add_header_above(c("Year", "NZL (DOM)" = 2, "NZL (JV)" = 2), align = c('l', rep('c', times = sefra_data$n_fishery_groups))) %>% sub("\\\\cmidrule\\(l\\{3pt\\}r\\{3pt\\}\\)\\{1-1\\}", "", .)
}
if (grepl("jpn", output_path)) {
    tab <- tab %>% add_header_above(c("Year", "JPN" = 2), align = c('l', rep('c', times = sefra_data$n_fishery_groups))) %>% sub("\\\\cmidrule\\(l\\{3pt\\}r\\{3pt\\}\\)\\{1-1\\}", "", .)
}
if (grepl("twn", output_path)) {
    tab <- tab %>% add_header_above(c("Year", "TWN" = 2), align = c('l', rep('c', times = sefra_data$n_fishery_groups))) %>% sub("\\\\cmidrule\\(l\\{3pt\\}r\\{3pt\\}\\)\\{1-1\\}", "", .)
}

tab <- tab  %>%
  sub("\\\\toprule", "", .) %>%
  sub("\\\\bottomrule", "", .)

save_kable(tab, file = file.path(output_path, "tables", "q_per_species_group_fishery_group.tex"))
write.csv(x,    file = file.path(output_path, "tables", "q_per_species_group_fishery_group.csv"))

Table of catchability estimates with CI in brackets:

x <- rstan::extract(mdl_fit, pars = "q", permute = TRUE)[[1]]
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), fishery_group = sefra_data$fishery_groups,
                    species_group = sefra_data$species_groups)

x <- reshape2::melt(x) %>% group_by(fishery_group, species_group) %>% summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>% ungroup(.)

## `summarise()` has grouped output by 'fishery_group'. You can override using the
## `.groups` argument.

x <- x %>% mutate(., value = paste0(round(mean, 2), " (", round(low_ci, 2), "-", round(high_ci, 2), ")"))

x <- x %>% tidyr::pivot_wider(id_cols = species_group, values_from = value, names_from = fishery_group)

tab <- kable(x, format = 'latex',
             col.names = stringr::str_to_sentence(colnames(x)),
             row.names = FALSE,
             caption = "Catchability",
             label = 'q_per_species_group_fishery_group',
             align = c('l', rep('r', times = ncol(x) - 1)),
             linesep = "", booktabs = TRUE) %>%
    kable_styling(font_size = 10, latex_options = "hold_position") %>%
  sub("\\\\toprule", "", .) %>%
  sub("\\\\bottomrule", "", .)

save_kable(tab, file = file.path(output_path, "tables", "q_per_species_group_fishery_group_alternative.tex"))
write.csv(x, file = file.path(output_path, "tables", "q_per_species_group_fishery_group_alternative.csv"))

Catchability
Species_group	Nzl (dom)	Nzl (jv)	Jpn	Twn
Wandering albatross	4.87 (3.41-6.51)	0.05 (0.01-0.14)	15.31 (12.74-17.92)	3.05 (2.32-3.79)
Royal albatross	3.3 (1.67-5.53)	0.08 (0.01-0.26)	6.88 (3.52-11.24)	2.32 (0.61-5.16)
Mollymawk	2.17 (1.87-2.49)	0.21 (0.18-0.25)	4.39 (3.96-4.83)	0.8 (0.7-0.91)
Sooty albatross	5.87 (0.18-30.27)	0.39 (0.02-2.2)	17.63 (13.97-22.1)	5.1 (3.92-6.73)
Medium petrel	4.07 (3.06-5.18)	0.11 (0.04-0.21)	3.18 (2.58-3.87)	0.66 (0.54-0.82)

Table of alpha_fishery_group coefficients (adjusted by alpha_intercept):

alpha_0 <- rstan::extract(mdl_fit, pars = "alpha_intercept", permute = TRUE)[[1]]
dimnames(alpha_0) <- list(iteration = 0.5 * n_iter + (1:dim(alpha_0)[1]), fishery_group = sefra_data$fishery_groups)

alpha_fg <- rstan::extract(mdl_fit, pars = "alpha_fishery_group", permute = TRUE)[[1]]
dimnames(alpha_fg) <- list(iteration = 0.5 * n_iter + (1:dim(alpha_fg)[1]), fishery_group = sefra_data$fishery_groups)

# take sum of catchability intercept and fishery group coefficients
# and (back)transform to nominal scale
alpha_fg_adj <- alpha_0 + alpha_fg
alpha_fg_adj <- exp(alpha_fg_adj)

alpha_fg_adj <- reshape2::melt(alpha_fg_adj)

tab_dat <- alpha_fg_adj %>% group_by(fishery_group) %>% summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>% ungroup(.) 

tab_dat <- tab_dat %>% mutate(., alpha_fg = paste0(round(mean, 2), " (", round(low_ci, 2), "-", round(high_ci, 2), ")"))

tab <- kable(tab_dat %>% select(., fishery_group, alpha_fg), format = 'latex',
             col.names = c("Fishery group", "alpha_fishery_group"),
             row.names = FALSE,
             caption = "Catchability coefficients per fishery group.",
             label = 'catchability_coefficients_fishery_group',
             align = c('l', rep('r', times = ncol(tab_dat) - 1)),
             linesep = "", booktabs = TRUE) %>%
    kable_styling(font_size = 10, latex_options = "hold_position") %>%
  sub("\\\\toprule", "", .) %>%
  sub("\\\\bottomrule", "", .)

save_kable(tab, file = file.path(output_path, "tables", "catchability_coefficients_fishery_group.tex"))
write.csv(tab_dat, file = file.path(output_path, "tables", "catchability_coefficients_fishery_group.csv"))

Catchability coefficients per fishery group.
Fishery group	alpha_fishery_group
NZL (DOM)	3.32 (1.86-5.35)
NZL (JV)	0.1 (0.04-0.2)
JPN	7.56 (6.54-8.61)
TWN	1.76 (1.37-2.15)

Catchability intercept:

x <- rstan::extract(mdl_fit, pars = "alpha_intercept", permute = FALSE)
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), chain = paste0("chain_", 1:dim(x)[2]),
                    parameter = dimnames(x)[[3]])

x <- reshape2::melt(x)

x <- x %>% group_by(., parameter) %>%
  summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>%
  ungroup(.)

kable(x, digits = 2)

parameter	mean	low_ci	high_ci
alpha_intercept[1]	0.34	0.08	0.61
alpha_intercept[2]	0.34	0.08	0.61
alpha_intercept[3]	0.34	0.08	0.61
alpha_intercept[4]	0.34	0.08	0.61

Fishery group catchability parameters:

x <- rstan::extract(mdl_fit, pars = "alpha_fishery_group", permute = FALSE)
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), chain = paste0("chain_", 1:dim(x)[2]),
                    parameter = dimnames(x)[[3]])

x <- reshape2::melt(x)

x <- x %>% group_by(., parameter) %>%
  summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>%
  ungroup(.)

kable(x, digits = 2)

parameter	mean	low_ci	high_ci
alpha_fishery_group[1]	0.82	0.36	1.27
alpha_fishery_group[2]	-2.72	-3.36	-2.14
alpha_fishery_group[3]	1.68	1.44	1.94
alpha_fishery_group[4]	0.22	-0.06	0.48

Species group catchability parameters

x <- rstan::extract(mdl_fit, pars = "alpha_species_group", permute = FALSE)
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), chain = paste0("chain_", 1:dim(x)[2]),
                    parameter = dimnames(x)[[3]])

x <- reshape2::melt(x)

x <- x %>% group_by(., parameter) %>%
  summarise(., mean = mean(value), low_ci = quantile(value, prob = 0.025), high_ci = quantile(value, prob = 0.975)) %>%
  ungroup(.)

kable(x, digits = 2)

parameter	mean	low_ci	high_ci
alpha_species_group[1,1]	0.40	-0.18	1.01
alpha_species_group[2,1]	-0.92	-2.30	0.26
alpha_species_group[3,1]	0.70	0.51	0.89
alpha_species_group[4,1]	0.55	0.22	0.86
alpha_species_group[1,2]	-0.01	-0.76	0.70
alpha_species_group[2,2]	-0.54	-2.08	0.77
alpha_species_group[3,2]	-0.13	-0.64	0.29
alpha_species_group[4,2]	0.15	-0.83	0.92
alpha_species_group[1,3]	-0.39	-0.92	0.14
alpha_species_group[2,3]	0.84	0.06	1.64
alpha_species_group[3,3]	-0.54	-0.70	-0.39
alpha_species_group[4,3]	-0.79	-1.01	-0.52
alpha_species_group[1,4]	-0.23	-2.33	1.69
alpha_species_group[2,4]	0.56	-1.31	2.63
alpha_species_group[3,4]	0.84	0.62	1.08
alpha_species_group[4,4]	1.06	0.77	1.38
alpha_species_group[1,5]	0.23	-0.31	0.81
alpha_species_group[2,5]	0.06	-0.79	0.93
alpha_species_group[3,5]	-0.87	-1.07	-0.67
alpha_species_group[4,5]	-0.97	-1.21	-0.70

Recover estimated catchabilities from model parameters:

## Get catchability parameters
qq <- rstan::extract(mdl_fit, pars = "q", permute = TRUE)[[1]]
dimnames(qq) <- list(iteration = 0.5 * n_iter + (1:dim(qq)[1]), fishery_group = sefra_data$fishery_groups,
                    species_group = sefra_data$species_groups)
qq <- reshape2::melt(qq) %>% rename(., q = value)

## beta_0 (intercept - same for all fishery groups)
b_0 <- rstan::extract(mdl_fit, pars = "alpha_intercept", permute = TRUE)[[1]]
dimnames(b_0) <- list(iteration = 0.5 * n_iter + (1:dim(b_0)[1]), fishery_group = sefra_data$fishery_groups)
b_0 <- reshape2::melt(b_0) %>% rename(., b_0 = value)

## beta_f
b_f <- rstan::extract(mdl_fit, pars = "alpha_fishery_group", permute = TRUE)[[1]]
dimnames(b_f) <- list(iteration = 0.5 * n_iter + (1:dim(b_f)[1]), fishery_group = sefra_data$fishery_groups)
b_f <- reshape2::melt(b_f) %>% rename(., b_f = value)

## beta_fg
b_fg <- rstan::extract(mdl_fit, pars = "alpha_species_group", permute = TRUE)[[1]]
dimnames(b_fg) <- list(iteration = 0.5 * n_iter + (1:dim(b_fg)[1]), fishery_group = sefra_data$fishery_groups,
                    species_group = sefra_data$species_groups)
b_fg <- reshape2::melt(b_fg) %>% rename(., b_fg = value)

## combine betas and recover q
qq_test <- b_fg
qq_test <- b_f %>% left_join(qq_test, ., by = c("iteration", "fishery_group"))
qq_test <- b_0 %>% left_join(qq_test, ., by = c("iteration", "fishery_group"))
qq_test <- qq_test %>% mutate(., q = exp(b_0 + b_f + b_fg))

## check recovered q is equal to that provided by the model
all.equal(qq$q, qq_test$q)

## [1] TRUE

Figures of estimated quantities

`pi` parameters

pi parameters estimate the probability of identifying captures of a specific taxa to different taxonomic resolutions.

Estimated pi parameters:

# transfrom from logit to probabilities
inv_logit <- binomial()$linkinv
plt_dat <- out[["p_obs_logit"]]
plt_dat <- plt_dat %>% mutate(., value = inv_logit(value))

# rescale to sum to 1
plt_dat <- plt_dat %>% group_by(., iter, fishery_group, pi) %>%
  mutate(., value = value / sum(value)) %>% ungroup(.)

plt_dat <- plt_dat %>%
  group_by(., fishery_group, pi, resolution) %>%
  summary_fn_boxplot(., value) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group', 'pi'. You can override
## using the `.groups` argument.

if(length(levels(plt_dat$pi)) == 1) {
  gg <- plt_dat %>%
    ggplot() +
    geom_boxplot(aes(x = fishery_group, ymin = lq, lower = lmq, middle = mq, upper = umq, ymax = uq, fill = resolution), stat = "identity", width = 0.75) +
    coord_cartesian(ylim = c(0, NA)) +
    scale_fill_manual("Taxonomic\nresolution", values = cols_taxo_resolution) +
    theme_bw(base_size = 12) +
    theme(legend.position = "right") +
    labs(x = "Fisheries group", y = "Probability of taxonomic resolutions for capture identifications")
  print(gg)
  ggsave(gg, file = file.path(output_path, "figures", "pi_boxplot.png"), width = 8, height = 6)
  ggsave(gg, file = file.path(output_path, "figures", "pi_boxplot.pdf"), width = 8, height = 6)
}

if(length(levels(plt_dat$pi)) > 1) {
  for(i in 1:length(fg_group_levels)) {
    gg <- plt_dat %>%
      filter(., fishery_group %in% fg_group_levels[i]) %>%
      ggplot() + geom_boxplot(aes(x = pi, ymin = lq, lower = lmq, middle = mq, upper = umq, ymax = uq, fill = resolution), stat = "identity", width = 0.75) +
      facet_wrap(vars(fishery_group)) +
      coord_cartesian(ylim = c(0, NA)) +
      scale_fill_manual("Taxonomic\nresolution", values = cols_taxo_resolution) +
      theme_bw(base_size = 12) +
      theme(legend.position = "right") +
      labs(x = "Pi group", y = "Probability of taxonomic resolutions")
    print(gg)
    ggsave(gg, file = file.path(output_path, "figures", paste0("pi_boxplot_fishery_", i, ".png")), width = 8, height = 6)
    ggsave(gg, file = file.path(output_path, "figures", paste0("pi_boxplot_fishery_", i, ".pdf")), width = 8, height = 6)
  }
}

Estimated catchabilities

Estimated catchabilities (log-scale):

x <- rstan::extract(mdl_fit, pars = "q", permute = TRUE)[[1]]
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), fishery_group = sefra_data$fishery_groups,
                    species_group = sefra_data$species_groups)

x <- reshape2::melt(x)

plt_dat <- x %>% group_by(., fishery_group, species_group) %>%
  summarise(., mean = mean(value), mq = median(value),
            lq = as.numeric(quantile(value, probs = 0.025)),
            uq = as.numeric(quantile(value, probs = 0.975))) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group'. You can override using the
## `.groups` argument.

gg <- plt_dat %>%
  ggplot(aes(x = species_group, y = mq, col = species_group)) +
  geom_pointrange(aes(ymin = lq, ymax = uq), alpha = 0.75, size = 0.75, linewidth = 0.75) +
  scale_colour_manual("", values = cols_species_group) +
  facet_wrap(fishery_group ~ ., ncol = 1) +
  scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "none") +
  labs(x = "Species group", y = "Estimated catchability (q)") +
  scale_y_log10()
gg

ggsave(gg, file = file.path(output_path, "figures", "log10_catchabilities.png"), width = 6, height = 10)
ggsave(gg, file = file.path(output_path, "figures", "log10_catchabilities.pdf"), width = 6, height = 10)

Boxplot of estimated catchabilities (natural-scale):

x <- rstan::extract(mdl_fit, pars = "q", permute = TRUE)[[1]]
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), fishery_group = sefra_data$fishery_groups,
                    species_group = sefra_data$species_groups)

x <- reshape2::melt(x)

plt_dat <- x %>% group_by(., fishery_group, species_group) %>%
  summary_fn_boxplot(., value) %>% ungroup()

## `summarise()` has grouped output by 'fishery_group'. You can override using the
## `.groups` argument.

gg <- plt_dat %>%
  ggplot() +
  geom_boxplot(aes(x = fishery_group, ymin = lq, lower = lmq, middle = mq, upper = umq, ymax = uq, fill = species_group), stat = "identity", width = 0.75) +
  scale_fill_manual("Species group", values = cols_species_group) +
  theme_bw(base_size = 12) +
  theme(legend.position = "right") +
  labs(x = "Fisheries group", y = "Estimated catchability (q)")
gg

ggsave(gg, file = file.path(output_path, "figures", "catchabilities_boxplot.png"), width = 8, height = 6)
ggsave(gg, file = file.path(output_path, "figures", "catchabilities_boxplot.pdf"), width = 8, height = 6)

Boxplot of estimated catchabilities (log-scale):

gg <- plt_dat %>%
  ggplot() +
  geom_boxplot(aes(x = fishery_group, ymin = lq, lower = lmq, middle = mq, upper = umq, ymax = uq, fill = species_group), stat = "identity", width = 0.75) +
  scale_fill_manual("Species group", values = cols_species_group) +
  theme_bw(base_size = 12) +
  theme(legend.position = "right") +
  labs(x = "Fisheries group", y = "Estimated catchability (q)") +
  scale_y_log10()
gg

ggsave(gg, file = file.path(output_path, "figures", "log10_catchabilities_boxplot.png"), width = 8, height = 6)
ggsave(gg, file = file.path(output_path, "figures", "log10_catchabilities_boxplot.pdf"), width = 8, height = 6)

Estimated catchabilities per fleet (log-scale):

x <- rstan::extract(mdl_fit, pars = "alpha_fishery_group", permute = TRUE)[[1]]
dimnames(x) <- list(iteration = 0.5 * n_iter + (1:dim(x)[1]), fishery_group = sefra_data$fishery_groups)

x <- reshape2::melt(x)

plt_dat <- x %>% group_by(., fishery_group) %>%
  summarise(., mean = mean(value), mq = median(value),
            lq = as.numeric(quantile(value, probs = 0.025)),
            uq = as.numeric(quantile(value, probs = 0.975))) %>% ungroup()

gg <- plt_dat %>%
  ggplot(aes(x = fishery_group, y = mq, col = fishery_group)) +
  geom_pointrange(aes(ymin = lq, ymax = uq), alpha = 0.75, size = 0.75, linewidth = 0.75) +
  theme_bw(base_size = 12) +
  theme(legend.position = "none") +
  labs(x = "Fishery group", y = "Estimated catchability coefficient")
gg

ggsave(gg, file = file.path(output_path, "figures", "log10_catchability_coef_fishery_group.png"), width = 6, height = 10)
ggsave(gg, file = file.path(output_path, "figures", "log10_catchability_coef_fishery_group.pdf"), width = 6, height = 10)

Maps of density overlap, observed captures and deaths by species group

Maps of density overlap:

dfr <- out[["map_overlap"]]
dfr <- dfr %>% mutate(., cell = as.numeric(as.character(cell)))

## Sum overlap across species group (by fishery group)
dfr <- species_groups(sefra_data, print_species = TRUE) %>% select(., species, species_group) %>%
  left_join(dfr, .)

## 5 species groups

## (with 25 species)

## Joining with `by = join_by(species)`

dfr <- dfr %>% group_by(., fishery_group, species_group, cell) %>%
  summarise(., value = sum(value)) %>% data.frame(.)

## `summarise()` has grouped output by 'fishery_group', 'species_group'. You can
## override using the `.groups` argument.

## Factorise species groups for appropriate ordering
dfr <- dfr %>% mutate(., species_group = factor(species_group, levels = sefra_data$species_groups))

## Include spatial information
dfr <- dfr %>% left_join(., grid, by = c("cell" = "id_cell")) %>% st_as_sf()
dfr <- dfr %>% filter(., value > 0)

gg <- ggplot() + 
    geom_sf(data = dfr, aes(fill = as.numeric(value)), col = NA, size = NA) +
    facet_grid(species_group ~ fishery_group) +
    #geom_sf(data = grid, fill = NA, size = NA) + 
    geom_sf(data = southern_map, fill = "grey") + 
    coord_sf(xlim = southern_bbox[c(1, 3)], ylim = southern_bbox[c(2, 4)], expand = FALSE) +
    #guides(fill = "none") +
    labs(fill = "Density\nOverlap") +
    viridis::scale_fill_viridis(direction = -1, trans = "log10") +
    theme_bw() +
    theme(axis.title = element_blank(),
          legend.position = "bottom",
          legend.key.width = unit(1.5, "cm"),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 10))
gg

ggsave(gg, file = file.path(output_path, "maps", "density_overlap.png"), height = 12)

## Saving 8 x 12 in image

ggsave(gg, file = file.path(output_path, "maps", "density_overlap.pdf"), height = 12)

## Saving 8 x 12 in image

map_caption <- paste0("Map of estimated density overlap from observed effort, by species group and fishery group.")
write(paste0("\\caption{", map_caption, "}"), file = file.path(output_path, "maps", "density_overlap_caption.tex"))

Maps of predicted and observed captures. The model predicts captures per species. These captures are down-weighted according to the probability of observation, to allow comparison with the observations.

# map of the overlap and captures per species and fishery group
x <- out[["map_capture"]] %>% group_by(cell, fishery_group, species) %>% summarise(value = sum(value)) %>% ungroup() %>% mutate(id_cell = as.integer(as.character(cell)))
x <- x %>% left_join(grid, by = "id_cell") %>% filter(value > 0)
y <- captures_o %>% filter(id_species %in% 1:27) %>% group_by(id_cell, fishery_group, code) %>% summarise(n_captures = sum(n_captures)) %>% ungroup()
y <- y %>% filter(n_captures > 0)

dfr <- full_join(x, y, by = c("id_cell" = "id_cell", "species" = "code", "fishery_group" = "fishery_group")) %>% filter(!is.na(value)) %>% mutate(n_captures = replace_na(n_captures, 0))
dfr <- dfr %>% st_as_sf()

# maps per species
df1 <- split(dfr, ~species)

for (i in names(df1)) {

    dfr_captures_sim <- df1[[i]]
    dfr_captures_emp <- df1[[i]] %>% st_centroid() %>% filter(n_captures > 0)
    
    dfr_caption <- dfr_captures_emp %>% group_by(fishery_group) %>% summarise(n_captures = sum(n_captures)) %>% st_drop_geometry()
    dfr_caption <- purrr::map2(dfr_caption$fishery_group, dfr_caption$n_captures, function(x, y) paste0(x, " = ", y)) %>% unlist() %>% paste(collapse = ", ")
    
    map_caption <- paste0("Map of the sum of the observed and predicted captures per fishery group for ", i, ". Total captures per fishery group are: ", ifelse(dfr_caption == "", "0", dfr_caption), ".")
    map_label   <- paste0("captures_emp_vs_sim_per_species_", i)

    write(paste0("\\caption{", map_caption, "}"), file = file.path(output_path, "maps", paste0(map_label, "_caption.tex")))
    
    gg <- ggplot() + 
        geom_sf(data = dfr_captures_sim, aes(fill = value), col = NA) +
        geom_sf(data = dfr_captures_emp, aes(size = n_captures, alpha = n_captures), col = "tomato") + 
        facet_wrap(~fishery_group) + theme_sh() + scale_fill_viridis(direction = -1) +
        labs(fill = "Predicted\nCaptures", size = "Observed\nCaptures") + guides(alpha = "none")

    ggfile <- file.path(output_path, "maps", paste0(map_label, ".png"))
    ggsave(gg, file = ggfile)
}

Maps of deaths:

dfr <- out[["map_deaths"]]
dfr <- dfr %>% mutate(., cell = as.numeric(as.character(cell)))

## Sum deaths across species group (by fishery group)
dfr <- species_groups(sefra_data, print_species = TRUE)  %>% select(., species, species_group) %>%
  left_join(dfr, .)

## 5 species groups

## (with 25 species)

## Joining with `by = join_by(species)`

dfr <- dfr %>% group_by(., fishery_group, species_group, cell) %>%
  summarise(., value = sum(value)) %>% data.frame(.)

## `summarise()` has grouped output by 'fishery_group', 'species_group'. You can
## override using the `.groups` argument.

## Factorise species groups for appropriate ordering
dfr <- dfr %>% mutate(., species_group = factor(species_group, levels = sefra_data$species_groups))

## Include spatial information
dfr <- dfr %>% left_join(., grid, by = c("cell" = "id_cell")) %>% st_as_sf()
dfr <- dfr %>% filter(., value > 0)

gg <- ggplot() + 
    geom_sf(data = dfr, aes(fill = as.numeric(value)), col = NA, size = NA) +
    facet_grid(species_group ~ fishery_group) +
    #geom_sf(data = grid, fill = NA, size = NA) + 
    geom_sf(data = southern_map, fill = "grey") + 
    coord_sf(xlim = southern_bbox[c(1, 3)], ylim = southern_bbox[c(2, 4)], expand = FALSE) +
    #guides(fill = "none") +
    labs(fill = "Deaths") +
    viridis::scale_fill_viridis(direction = -1, trans = "log10", option = "inferno") +
    theme_bw() +
    theme(axis.title = element_blank(),
          legend.position = "bottom",
          legend.key.width = unit(1.5, "cm"),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 10))
gg

ggsave(gg, file = file.path(output_path, "maps", "deaths.png"), height = 12)

## Saving 8 x 12 in image

ggsave(gg, file = file.path(output_path, "maps", "deaths.pdf"), height = 12)

## Saving 8 x 12 in image

map_caption <- paste0("Map of estimated deaths resulting from observed fishing effort, by species group and fishery group.")
write(paste0("\\caption{", map_caption, "}"), file = file.path(output_path, "maps", "deaths_caption.tex"))

##            used  (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells  4219201 225.4   10923793   583.4   10923793   583.4
## Vcells 12945815  98.8 2772907496 21155.7 3853039619 29396.4

Fit SEFRA model to combined dataset from 2024 SRA

Charles Edwards and Tom Peatman

26 Mar 2025

Introduction

Prelims

Fisheries and structural data

Estimate catchabilities from combined data

Create `sefraData` object

Model run

Prelims for processing fitted model

Model diagnostics

Trace plots of model parameters (and estimated quantities)

Updates to biological input distributions

Plots of rhat and Neff

Trace plots of model parameters

Observed vs predicted captures

Observed vs predicted cumulative captures by capture group

Observed vs predicted empirical captures by capture group

Tables of estimated quantities

Estimated captures

Estimated catchabilities

Figures of estimated quantities

`pi` parameters

Estimated catchabilities

Maps of density overlap, observed captures and deaths by species group

Fit SEFRA model to combined dataset from 2024 SRA

Charles Edwards and Tom Peatman

26 Mar 2025

Introduction

Prelims

Fisheries and structural data

Estimate catchabilities from combined data

Create sefraData object

Model run

Prelims for processing fitted model

Model diagnostics

Trace plots of model parameters (and estimated quantities)

Updates to biological input distributions

Plots of rhat and Neff

Trace plots of model parameters

Observed vs predicted captures

Observed vs predicted cumulative captures by capture group

Observed vs predicted empirical captures by capture group

Tables of estimated quantities

Estimated captures

Estimated catchabilities

Figures of estimated quantities

pi parameters

Estimated catchabilities

Maps of density overlap, observed captures and deaths by species group

Create `sefraData` object

`pi` parameters