Step 1. Format the data

cpue_data() turns the logbook (long format: one row per trip and species) into one row per operation with one column per species, merges species synonyms, and starts a processing log.

x <- cpue_data(sim$logbook, year = "year", species = "species", catch = "kg",
               effort = "hours", id = "trip", vessel = "vessel",
               lon = "lon", lat = "lat", covariates = "month",
               synonyms = c("Blackspot seabream" = "Pagellus bogaraveo"))
x <- subset(x, !is.na(lon), note = "positions reported")
x

## <cpue_data> 19218 operations, 10 species, years 2000 - 2024 
##   vessels: 40 
##   effort: hours 
##   position: lon / lat 
##   covariates: month 
##   main species: Scomber colias, Merluccius merluccius, Thunnus alalunga, Pagellus erythrinus, Pagellus bogaraveo ... 
## Processing log:
##                 step operations                             note
##  input (long format)      19574        46775 records, 10 species
##               subset      19218 positions reported (356 removed)

Guidelines

• Define the operation: the unit at which catch and effort are recorded (set, haul, day or trip). Here it is the trip, and about a fifth of the trips combine sets of two tactics; their catch composition is a mixture.
• Choose an effort measure that scales catch (here hours fished).
• Merge synonyms first (here blackspot seabream recorded under its common name in some trips); they change the species set.
• Every selection with subset() is written to the log (x$log).

Step 2. Explore the fleet

fl <- plot_fleet(x)
fl$records

fl$vessels

fl$vessel_years

Vessels enter and leave over the years. Each vessel’s catch profile shows how far the fleet is split into specialists:

groups <- list(demersal = sp_tab$species[sp_tab$group == "demersal"],
               pelagic = sp_tab$species[sp_tab$group == "pelagic"])
plot_vessel_profile(x, groups = groups)

head(vessel_summary(x), 8)

##   vessel records years years_qualified first last
## 1    V10     599    25              25  2000 2024
## 2    V16     599    25              25  2000 2024
## 3    V39     587    25              25  2000 2024
## 4    V38     585    25              25  2000 2024
## 5    V18     581    25              25  2000 2024
## 6    V01     578    25              25  2000 2024
## 7    V27     576    25              25  2000 2024
## 8    V05     575    25              25  2000 2024

Since the data are simulated, we can label the vessels with their true type:

type <- tapply(sim$trips$vessel_type, sprintf("V%02d", sim$trips$vessel),
               function(z) z[1])
plot_vessel_profile(x, groups = groups, by = type)

Guidelines

• Strong fleet turnover confounds vessel efficiency with year; random vessel effects help if vessels overlap in time.
• Specialists (vessels at either end of the profile plot) target the same species throughout; their targeting is absorbed by a vessel effect. Mixed vessels switch tactics with prices; that is the targeting that the catch composition has to capture.
• An indicator fleet (select_vessels()) keeps vessels that fished consistently, optionally with rules on the target species (min_years_target, min_share_target). Compare the index with and without the selection.

Step 3. Choose the species set

plot_species_curve(x, threshold = 0.95)

targets <- c("Pagellus bogaraveo", "Coryphaena hippurus")
x <- select_species(x, threshold = 0.95, keep = targets)
x$species_labels

##          Scomber.colias        Seriola.dumerili           Conger.conger 
##        "Scomber colias"      "Seriola dumerili"         "Conger conger" 
##        Thunnus.alalunga     Coryphaena.hippurus     Pagellus.erythrinus 
##      "Thunnus alalunga"   "Coryphaena hippurus"   "Pagellus erythrinus" 
##      Pagellus.bogaraveo   Merluccius.merluccius       Phycis.blennoides 
##    "Pagellus bogaraveo" "Merluccius merluccius"     "Phycis blennoides" 
##           Dentex.dentex 
##         "Dentex dentex"

# the species groups that are still in the data, for the selections below
demersal <- intersect(groups$demersal, x$species_labels)
pelagic <- intersect(groups$pelagic, x$species_labels)

Guidelines

• Rank species by occurrence (default) and check by weight.
• Always keep the species of interest, even if rare.
• Species characteristic of other tactics are valuable even if not abundant: they let the PCA separate the tactics.

Step 4. Nominal CPUE and catch composition

nom <- plot_nominal(x, targets)
nom$index

nom$positive

nom$distribution

plot_composition(x)

plot_cpue_map(x, res = 0.5, by_period = 8)

Expect the catch composition to shift towards pelagic species in 2008-2015, at the same time as the prices: a first sign of targeting.

Guidelines

• The share of positive trips and the distribution of positive CPUE guide the choice of error distribution (Tweedie or delta).
• Abrupt changes in composition at known changes in prices, regulations or gear point to targeting; gradual drifts may be abundance.

Step 6. Targeting covariates for the whole fleet

x <- add_targeting(x)
x$targeting$pca

## <cpue_pca> 19218 operations, 10 species, transform: fourth_root 
##  PC eigenvalue share cumulative
##   1      2.192 0.219      0.219
##   2      1.707 0.171      0.390
##   3      1.340 0.134      0.524
##   4      1.108 0.111      0.635
##   5      1.056 0.106      0.740
##   6      0.792 0.079      0.820
##   7      0.569 0.057      0.876
##   8      0.530 0.053      0.929
## Optimal Coordinates: 3 non-trivial PC(s); 4 clusters by the npc + 1 rule

plot_oc(x)

plot_biplot(x)

tactic_table(x)

##   FT operations share              dominant Scomber colias Seriola dumerili
## 1  1       9936 0.517        Scomber colias           38.3             11.6
## 2  2       3620 0.188      Thunnus alalunga            2.2              1.9
## 3  3       2849 0.148    Pagellus bogaraveo            1.6              1.6
## 4  4       2813 0.146 Merluccius merluccius            2.4              1.6
##   Conger conger Thunnus alalunga Coryphaena hippurus Pagellus erythrinus
## 1          10.8              0.7                 8.8                15.4
## 2           1.1             79.6                10.1                 2.2
## 3           3.2              0.7                 2.0                 9.0
## 4           2.4              1.0                 1.1                 2.6
##   Pagellus bogaraveo Merluccius merluccius Phycis blennoides Dentex dentex
## 1                0.5                   0.5               0.2          13.1
## 2                0.6                   0.7               0.4           1.2
## 3               75.5                   1.6               0.7           4.0
## 4                2.1                  73.2              10.9           2.7

pt <- plot_targeting(x, res = 0.5)
pt$time

pt$pcs_time

pt$map_tactic

Check the clusters against the true tactic of each trip (not possible with real data):

truth <- sim$trips[match(x$data$trip, sim$trips$trip), ]
round(prop.table(table(cluster = x$data$FT, true_tactic = truth$tactic), 1), 2)

##        true_tactic
## cluster Pelagic Reef Sediment
##       1    0.43 0.33     0.24
##       2    0.96 0.02     0.02
##       3    0.05 0.85     0.10
##       4    0.07 0.16     0.78

Across the whole fleet, the first axis is typically dominated by the contrast between pelagic and demersal trips (check the loadings). For a demersal species that contrast is mostly irrelevant: what matters is the targeting within the demersal trips. The worked examples therefore start by selecting the relevant tactics and recomputing the PCA.

Guidelines

• Name the tactics after their dominant species and check that they make sense to people who know the fishery.
• Number of PCs: Optimal Coordinates test, not AIC; the composition is not independent of the response, so AIC selects too many PCs.
• Continuous or discrete? PCs represent gradual changes and mixed trips; clusters assume distinct tactics.
• Run add_targeting() again after every selection; PCs and clusters depend on the operations included (select_tactics() does this).

# relative indices of several fits against the true biomass
compare_indices <- function(fits, true) {
  idx <- do.call(rbind, lapply(names(fits), function(n) {
    i <- cpue_index(fits[[n]]$fit, fits[[n]]$model)
    data.frame(year = i$year, relative = i$relative, lower = i$lower,
               upper = i$upper, index = n)
  }))
  true$true <- true$true / mean(true$true)
  ggplot(idx, aes(x = year, y = relative, colour = index, fill = index)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.12, colour = NA) +
    geom_line(linewidth = 0.8) +
    geom_line(data = true, aes(x = year, y = true), inherit.aes = FALSE,
              colour = "black", linewidth = 1.1) +
    geom_hline(yintercept = 1, linetype = 3) +
    labs(x = "Year", y = "Relative CPUE", colour = NULL, fill = NULL,
         subtitle = "black: true biomass") + theme_bw()
}

Step 5. Remove the trips that could not catch seabream

select_tactics() keeps the clusters whose catch is at least 50% demersal species and recomputes the PCA and clusters on the remaining trips.

xd <- select_tactics(x, species = demersal, min_share = 0.5)
round(xd$targeting$previous$share, 2)    # demersal share of each previous cluster

##    1    2    3    4 
## 0.41 0.06 0.94 0.94

tail(xd$log, 3)

##             step operations
## 4  add targeting      19218
## 5 select tactics       5662
## 6  add targeting       5662
##                                                                                note
## 4                      3 PC(s) used (OC test: 3), 4 clusters (clara on composition)
## 5 clusters 3, 4 kept (clusters with >= 50% of catch from 6 species) (13556 removed)
## 6                      4 PC(s) used (OC test: 4), 5 clusters (clara on composition)

An alternative is the habitat-based selection of Stephens and MacCall (2004), which predicts the presence of seabream from the other species in each trip:

sm <- stephens_maccall(x, target = "Pagellus bogaraveo")
sm

## <stephens_maccall> target: Pagellus bogaraveo 
##   critical probability: 0.35
##   operations kept: 3208 of 19218 (17%)
##   target present in 63% of kept and 8% of dropped operations

pl <- plot_stephens_maccall(sm)
pl$coefficients

pl$kept

Guidelines

• select_tactics(): select on a group of species (here all demersal species), not on the target alone. Clusters come from the catch composition; keeping only clusters rich in the target would select trips by their catch of the target, and part of its abundance signal with them.
• The removed clusters should be clearly different (here almost no demersal catch), not borderline.
• Stephens and MacCall: check the share kept by year (pl$kept); a trend means the selection itself follows the targeting.
• Compare the index with and without the selection (below).

Step 6. Targeting covariates of the demersal trips

plot_oc(xd)

plot_biplot(xd)

tactic_table(xd)

##   FT operations share              dominant Scomber colias Seriola dumerili
## 1  1       1585 0.280 Merluccius merluccius            2.8              1.5
## 2  2       1581 0.279    Pagellus bogaraveo            2.2              1.3
## 3  3       1245 0.220    Pagellus bogaraveo            1.0              1.9
## 4  4        839 0.148 Merluccius merluccius            1.5              1.9
## 5  5        412 0.073     Phycis blennoides            0.0              1.1
##   Conger conger Thunnus alalunga Coryphaena hippurus Pagellus erythrinus
## 1           2.1              1.2                 0.9                 2.0
## 2           0.0              1.0                 2.9                 3.9
## 3           6.2              0.5                 1.1                14.0
## 4           2.9              0.6                 1.5                 3.5
## 5           2.6              0.0                 1.7                 7.4
##   Pagellus bogaraveo Merluccius merluccius Phycis blennoides Dentex dentex
## 1                2.4                  71.4              13.0           2.6
## 2               83.4                   1.6               0.7           3.0
## 3               68.8                   0.9               0.4           5.1
## 4                1.7                  83.5               0.1           2.8
## 5                3.7                   0.4              77.3           5.7

pd <- plot_targeting(xd, res = 0.5)
pd$time

pd$pcs_time

The new PCs describe the targeting within the demersal fishery, e.g. reef against sediment species and shallow against deep (check the loadings).

Step 7. Fit the models

Terms are added in a meaningful order: year only, then space, then the targeting covariates; all models have a random vessel effect.

fit_sbr <- fit_cpue(xd, species = "Pagellus bogaraveo",
                    models = c("Base", "Base+S", "Cluster+S", "DPC+S"),
                    vessel = TRUE, spatial_k = 40, engine = "bam")
fit_sbr

## <cpue_fit> Pagellus bogaraveo - tweedie with random vessel effects 
##   operations: 5662 ; years: 2000 - 2024 
##   index averaged over: FT (and positions in spatial models) 
## Models:
##   Base         y ~ fyear + s(vessel_re, bs = "re")
##   Base+S       y ~ fyear + s(pos_x, pos_y, k = 40) + s(vessel_re, bs = "re")
##   Cluster+S    y ~ fyear + s(pos_x, pos_y, k = 40) + FT + s(vessel_re, bs = "re")
##   DPC+S        y ~ fyear + s(pos_x, pos_y, k = 40) + s(PC1, k = 4) + s(PC2, k = 4) + s(PC3, k = 4) + s(PC4, k = 4) + s(vessel_re, bs = "re")

Step 8. Diagnostics

cpue_diagnostics(fit_sbr)

##       model component    n     AIC     BIC dev_explained  edf tweedie_p
## 1      Base      main 5662 23850.1 24218.9        0.0427 53.1     1.435
## 2    Base+S      main 5662 22713.8 23332.0        0.2251 90.2     1.392
## 3 Cluster+S      main 5662 18824.1 19423.9        0.6409 87.6     1.314
## 4     DPC+S      main 5662 17529.8 18195.4        0.7321 95.4     1.336
##   vessel_sd   dAIC
## 1     0.253 6320.4
## 2     0.254 5184.0
## 3     0.271 1294.4
## 4     0.278    0.0

plot_index(fit_sbr, true = true_sbr)

plot_step(fit_sbr)

plot_influence(fit_sbr, "DPC+S")

res <- plot_residuals(fit_sbr, "DPC+S", by = "FT")
res$qq

res$year

res$by

pp_sbr <- ppc_cpue(fit_sbr, models = c("Base", "DPC+S"))
pp_sbr

## Predictive check: 50 replicates, with parameter uncertainty, tweedie 
## 
##  model  stat observed median  lower  upper tail_prob
##   Base zeros    0.473  0.472  0.456  0.486      0.48
##   Base  mean    2.950  2.950  2.820  3.140      0.52
##   Base   q90   12.000 12.300 11.600 12.800      0.82
##   Base   q99   24.300 24.000 22.200 25.700      0.38
##  DPC+S zeros    0.473  0.477  0.470  0.487      0.80
##  DPC+S  mean    2.950  2.990  2.890  3.110      0.86
##  DPC+S   q90   12.000 12.800 12.400 13.400      1.00
##  DPC+S   q99   24.300 23.500 22.000 24.800      0.12

pp_sbr$plots$distribution

pp_sbr$plots$zeros

pp_sbr$plots$tails

eff <- plot_effects(fit_sbr, "DPC+S")
eff$numeric

eff$space

plot_vessel_effects(fit_sbr, "DPC+S")

fit_sbr_all <- fit_cpue(x, species = "Pagellus bogaraveo", models = "DPC+S",
                        vessel = TRUE, spatial_k = 40, engine = "bam")
x_sm <- add_targeting(subset(x, sm$selected, note = "Stephens & MacCall, seabream"))
fit_sbr_sm <- fit_cpue(x_sm, species = "Pagellus bogaraveo", models = "DPC+S",
                       vessel = TRUE, spatial_k = 40, engine = "bam")
compare_indices(list(`demersal tactics` = list(fit = fit_sbr, model = "DPC+S"),
                     `Stephens & MacCall` = list(fit = fit_sbr_sm, model = "DPC+S"),
                     `all trips` = list(fit = fit_sbr_all, model = "DPC+S")),
                true_sbr)

Step 9. The index

idx_sbr <- cpue_index(fit_sbr, "DPC+S")
knitr::kable(idx_sbr, digits = 3)

export_index(fit_sbr, "seabream_index_DPC.csv", model = "DPC+S")

Step 5. Keep the pelagic trips

xp <- select_tactics(x, species = pelagic, min_share = 0.5)
round(xp$targeting$previous$share, 2)

##    1    2    3    4 
## 0.59 0.94 0.06 0.06

tail(xp$log, 3)

##             step operations
## 4  add targeting      19218
## 5 select tactics      13556
## 6  add targeting      13556
##                                                                               note
## 4                     3 PC(s) used (OC test: 3), 4 clusters (clara on composition)
## 5 clusters 1, 2 kept (clusters with >= 50% of catch from 4 species) (5662 removed)
## 6                     2 PC(s) used (OC test: 2), 3 clusters (clara on composition)

Outside its habitat (demersal trips) a zero says nothing about dolphinfish abundance, and the PCA of the remaining trips can now describe the targeting among pelagic species.

Step 6. Targeting covariates of the pelagic trips

plot_oc(xp)

plot_biplot(xp)

tactic_table(xp)

##   FT operations share            dominant Scomber colias Seriola dumerili
## 1  1       6558 0.484      Scomber colias           48.3             14.3
## 2  2       3507 0.259    Thunnus alalunga            2.0              1.5
## 3  3       3491 0.258 Pagellus erythrinus           10.1              4.5
##   Conger conger Thunnus alalunga Coryphaena hippurus Pagellus erythrinus
## 1           8.4              1.1                11.9                 0.4
## 2           1.2             82.3                 9.3                 0.8
## 3          16.8              0.7                 1.3                56.4
##   Pagellus bogaraveo Merluccius merluccius Phycis blennoides Dentex dentex
## 1                0.6                   0.6               0.2          14.4
## 2                0.6                   0.7               0.4           1.2
## 3                0.3                   0.4               0.3           9.3

pp_t <- plot_targeting(xp, by = "month", res = 0.5)
pp_t$time

pp_t$by

## $pc

## 
## $tactic

The PCs should now separate offshore species (albacore, dolphinfish) from coastal ones (amberjack, mackerel); check the loadings.

Step 7. Fit the models

Month enters as a cyclic smoother in every model except Base; the index is averaged over the observed months.

fit_dol <- fit_cpue(xp, species = "Coryphaena hippurus",
                    models = c("Base", "Covariates", "Cluster+S", "DPC+S"),
                    terms = 's(month, bs = "cc", k = 6)',
                    knots = list(month = c(0.5, 12.5)),
                    vessel = TRUE, spatial_k = 40, engine = "bam")
fit_dol

## <cpue_fit> Coryphaena hippurus - tweedie with random vessel effects 
##   operations: 13556 ; years: 2000 - 2024 
##   index averaged over: FT, month (and positions in spatial models) 
## Models:
##   Base         y ~ fyear + s(vessel_re, bs = "re")
##   Covariates   y ~ fyear + s(month, bs = "cc", k = 6) + s(vessel_re, bs = "re")
##   Cluster+S    y ~ fyear + s(pos_x, pos_y, k = 40) + s(month, bs = "cc", k = 6) + FT + s(vessel_re, bs = "re")
##   DPC+S        y ~ fyear + s(pos_x, pos_y, k = 40) + s(month, bs = "cc", k = 6) + s(PC1, k = 4) + s(PC2, k = 4) + s(vessel_re, bs = "re")

Step 8. Diagnostics

cpue_diagnostics(fit_dol)

##        model component     n     AIC     BIC dev_explained  edf tweedie_p
## 1       Base      main 13556 30865.1 31348.0        0.0779 62.1     1.669
## 2 Covariates      main 13556 29756.3 30269.5        0.1919 66.1     1.644
## 3  Cluster+S      main 13556 25041.4 25801.8        0.5618 98.8     1.549
## 4      DPC+S      main 13556 19234.5 19953.3        0.8050 92.9     1.479
##   vessel_sd    dAIC
## 1     0.816 11630.5
## 2     0.791 10521.7
## 3     0.635  5806.9
## 4     0.261     0.0

plot_index(fit_dol, true = true_dol)

plot_step(fit_dol)

plot_influence(fit_dol, "DPC+S")

res_dol <- plot_residuals(fit_dol, "DPC+S", by = "month")
res_dol$qq

res_dol$year

res_dol$by

pp_dol_m <- ppc_cpue(fit_dol, models = c("Base", "DPC+S"), by = "month")
pp_dol_m$plots$zeros

pp_dol <- ppc_cpue(fit_dol, models = c("Base", "DPC+S"), by = "year")
pp_dol

## Predictive check: 50 replicates, with parameter uncertainty, tweedie 
## 
##  model  stat observed median  lower  upper tail_prob
##   Base zeros    0.715  0.715  0.710  0.723      0.54
##   Base  mean    0.954  0.967  0.907  1.030      0.62
##   Base   q90    9.610  9.190  8.660  9.730      0.06
##   Base   q99   25.800 24.700 21.500 27.800      0.20
##  DPC+S zeros    0.715  0.716  0.711  0.721      0.68
##  DPC+S  mean    0.954  0.988  0.940  1.020      0.90
##  DPC+S   q90    9.610 10.200  9.570 10.800      0.94
##  DPC+S   q99   25.800 25.600 23.800 27.700      0.46

pp_dol$plots$distribution

pp_dol$plots$zeros

pp_dol$plots$tails

fit_dol_dl <- fit_cpue(xp, species = "Coryphaena hippurus", models = "DPC+S",
                       family = "delta_lognormal",
                       terms = 's(month, bs = "cc", k = 6)',
                       knots = list(month = c(0.5, 12.5)),
                       vessel = TRUE, spatial_k = 40, engine = "bam")
pp_dl <- ppc_cpue(fit_dol_dl, by = "month")
rbind(data.frame(family = "Tweedie",
                 pp_dol_m$overall[pp_dol_m$overall$model == "DPC+S", ]),
      data.frame(family = "delta-lognormal", pp_dl$overall))

##            family model group  stat   observed     median      lower      upper
## 5         Tweedie DPC+S   All zeros  0.7148126  0.7162511  0.7110191  0.7211973
## 6         Tweedie DPC+S   All  mean  0.9539995  0.9877376  0.9401629  1.0223080
## 7         Tweedie DPC+S   All   q90  9.6052632 10.1706136  9.5709359 10.7872752
## 8         Tweedie DPC+S   All   q99 25.7524357 25.5919369 23.8012316 27.7413133
## 1 delta-lognormal DPC+S   All zeros  0.7148126  0.7152552  0.7085036  0.7207952
## 2 delta-lognormal DPC+S   All  mean  0.9539995  1.1408306  1.0670850  1.2713289
## 3 delta-lognormal DPC+S   All   q90  9.6052632 10.0403629  9.1416388 10.9862446
## 4 delta-lognormal DPC+S   All   q99 25.7524357 44.9223675 39.4275724 53.2216513
##   tail_prob
## 5      0.68
## 6      0.90
## 7      0.94
## 8      0.46
## 1      0.58
## 2      1.00
## 3      0.88
## 4      1.00

pp_dl$plots$distribution

pp_dl$plots$zeros

eff_dol <- plot_effects(fit_dol, "DPC+S")
eff_dol$numeric

plot_vessel_effects(fit_dol, "DPC+S")

fit_dol_all <- fit_cpue(x, species = "Coryphaena hippurus", models = "DPC+S",
                        terms = 's(month, bs = "cc", k = 6)',
                        knots = list(month = c(0.5, 12.5)),
                        vessel = TRUE, spatial_k = 40, engine = "bam")
compare_indices(list(`Tweedie, pelagic tactics` = list(fit = fit_dol, model = "DPC+S"),
                     `delta-lognormal, pelagic tactics` = list(fit = fit_dol_dl,
                                                              model = "DPC+S"),
                     `Tweedie, all trips` = list(fit = fit_dol_all, model = "DPC+S")),
                true_dol)

Step 9. The index

idx_dol <- cpue_index(fit_dol, "DPC+S")
knitr::kable(idx_dol, digits = 3)