---
title: "Analysing real catch and effort data: a step-by-step guide"
author: "Henning Winker<br><span style=\"font-size: 0.8em;\">FAO Fisheries and Aquaculture Division, Assessment and Management Team (NFIFM)<br><a href=\"mailto:henning.winker@fao.org\">henning.winker@fao.org</a></span>"
output:
  rmarkdown::html_document:
    toc: true
    toc_float: true
    theme: readable
vignette: >
  %\VignetteIndexEntry{Analysing real catch and effort data: a step-by-step guide}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
                      fig.width = 8, fig.height = 4.5, dpi = 110,
                      eval = requireNamespace("cluster", quietly = TRUE))
library(cpueTargeting)
library(ggplot2)
```

This guide shows how to prepare fishery-dependent catch and effort data for
CPUE standardisation when fishers change what they target. It generalises the
protocol developed for blackspot seabream (GFCM) to any multispecies fishery,
from the raw logbook to the standardised index.

| Step | Question | Functions |
|---|---|---|
| 1. Format | Which column is what? One row per operation? | `cpue_data()`, `subset()` |
| 2. Fleet | Who fished when? Is reporting consistent? Which vessels are indicators? | `plot_fleet()`, `vessel_summary()`, `select_vessels()`, `plot_vessel_profile()` |
| 3. Species | Which species carry the targeting signal? | `species_curve()`, `plot_species_curve()`, `select_species()` |
| 4. Nominal CPUE | What do the raw data say? | `nominal_cpue()`, `plot_nominal()`, `plot_composition()`, `plot_cpue_map()` |
| 5. Subsetting (optional) | Which operations fished the target's habitat? | `stephens_maccall()`, `plot_stephens_maccall()`, `select_tactics()` (after Step 6) |
| 6. Targeting covariates | How did fishing tactics change? | `add_targeting()`, `plot_oc()`, `plot_biplot()`, `tactic_table()`, `plot_targeting()` |
| 7. Models | Base, Covariates, DPC, Cluster, spatial, vessel effects | `fit_cpue()` |
| 8. Diagnostics | Which terms matter, and do the models fit? Is the error distribution right? | `cpue_diagnostics()`, `plot_index()`, `plot_step()`, `plot_influence()`, `plot_residuals()`, `ppc_cpue()`, `plot_effects()`, `plot_vessel_effects()` |
| 9. Index | The index for the assessment | `cpue_index()`, `export_index()` |

Every step returns the same kind of object, and every step that removes
operations writes a line into a processing log, so the final data set can be
reported exactly.

## Example data

To show what each step should find, the example uses a logbook generated by
the bathymetric simulator (`simulate_bathymetric()`): four species live along
a depth gradient and 30 vessels shift effort between them as prices change
(years 1-10, 11-20, 21-30). In real data the truth is unknown; here it lets us
check that the tools recover the changes in targeting.

We turn the simulation into a typical logbook extract: one row per trip and
species, with catch in kg, hours fished, month, a coarse position, two rare
bycatch species, and one species recorded under two names.

```{r data}
sim <- simulate_bathymetric(seed = 123)
set.seed(1)
sets <- as_set_data(sim$Data)
sets$hours <- round(runif(nrow(sets), 4, 12), 1)
sets$month <- sample(1:12, nrow(sets), replace = TRUE)
names_sp <- c(Sp1 = "Eopsetta jordani", Sp2 = "Microstomus pacificus",
              Sp3 = "Sebastes crameri", Sp4 = "Anoplopoma fimbria")
long <- do.call(rbind, lapply(names(names_sp), function(s) {
  data.frame(trip = sets$SetID, year = sets$year + 1999, vessel = sets$vessel,
             hours = sets$hours, month = sets$month, lon = sets$Xobs, lat = sets$Yobs,
             species = names_sp[[s]], kg = sets[[s]] * sets$hours)
}))
# rare bycatch and a synonym
rare <- data.frame(trip = sample(sets$SetID, 400), species = "Raja spp.", kg = rexp(400))
rare2 <- data.frame(trip = sample(sets$SetID, 80), species = "Squalus suckleyi", kg = rexp(80))
extra <- rbind(rare, rare2)
extra <- cbind(long[match(extra$trip, long$trip), c("year", "vessel", "hours", "month",
                                                    "lon", "lat")], extra)
long <- rbind(long, extra[, names(long)])
syn <- long$species == "Anoplopoma fimbria" & runif(nrow(long)) < 0.1
long$species[syn] <- "Sablefish"
logbook <- long[long$kg > 0, ]
head(logbook)
```

## Step 1. Format the data

`cpue_data()` takes the logbook in long format (one row per operation and
species, as above) or wide format (one row per operation, one column per
species) and returns one row per operation. Tell it which column plays which
role.

```{r format}
x <- cpue_data(logbook, year = "year", species = "species", catch = "kg",
               effort = "hours", id = "trip", vessel = "vessel",
               lon = "lon", lat = "lat", covariates = "month",
               synonyms = c("Sablefish" = "Anoplopoma fimbria"))
x
```

> **Guidelines**
>
> * Define the **operation** carefully: the unit at which effort and catch
>   composition are recorded (set, haul, fishing day or trip). Catch
>   composition reflects targeting best at the finest unit available; trips
>   that span several fishing grounds or gears mix tactics.
> * Choose an **effort** measure that scales catch (hours fished, hooks,
>   days at sea). Without effort, the catch per operation is used.
> * Merge **synonyms** and drop non-informative categories (bait, discards,
>   "unidentified") before anything else; they change the species set.
> * Use `subset()` for any further rule (e.g. single-day trips, a gear,
>   a period). Each call is recorded in the log.

```{r subset, eval = FALSE}
x <- subset(x, month %in% 1:12, note = "all months (example)")
```

## Step 2. Explore the fleet

Standardisation corrects for the non-random allocation of effort. First look
at how operations are spread over years and vessels.

```{r fleet, fig.height = 6}
fl <- plot_fleet(x, by = "month")
fl$records
fl$vessels
fl$vessel_years
```

> **Guidelines**
>
> * Strong **fleet turnover** (new vessels replacing old ones) confounds
>   vessel efficiency with year; random vessel effects help, but only if
>   vessels overlap in time.
> * An **indicator fleet** of vessels that fished consistently can make
>   reporting more homogeneous. `select_vessels(x, min_years, min_records)`
>   keeps vessels with at least `min_records` operations in at least
>   `min_years` years. Compare indices with and without the selection.
> * Rules on the **target species** keep vessels that fish its habitat
>   regularly: `target`, `min_years_target` (years with catches of the
>   target), `min_share_target` (share of the vessel's operations catching
>   it) and `min_catch_target`. `vessel_summary()` lists these quantities
>   for every vessel and `plot_vessel_profile()` shows each vessel's catch
>   split by species or species group, which separates specialists from
>   mixed vessels. Rules on the target select vessels by their catches of
>   it: prefer rules on years active and on the share of operations, and
>   always compare with the index of all vessels.

```{r vessels}
x_core <- select_vessels(x, min_years = 10, min_records = 5)
x_core$log
vs <- vessel_summary(x, target = "Anoplopoma fimbria")
head(vs)
x_target <- select_vessels(x, min_years = 10, min_records = 5,
                           target = "Anoplopoma fimbria", min_share_target = 0.3)
table(attr(x_target, "vessels")$selected)
plot_vessel_profile(x)
```

## Step 3. Choose the species set

The catch composition carries the targeting signal, but rare species add
noise and many zero proportions. The cumulative species curve ranks species
by how often they are recorded; the species needed to reach a cumulative share
(95% by default) form the species set.

```{r species}
species_curve(x)
plot_species_curve(x, threshold = 0.95)
```

```{r select-species}
x <- select_species(x, threshold = 0.95)
x$species_labels
```

> **Guidelines**
>
> * Use **occurrence** (default) as the main criterion and check
>   `by = "weight"`: a species caught rarely but in large amounts can be a
>   target of a distinct tactic.
> * Always **keep the species of interest** (`keep = "..."`), even if rare.
> * Species that are characteristic of other tactics are valuable even when
>   not abundant; they are what lets the PCA separate tactics.
> * `pool_other = TRUE` keeps the dropped species as one `Other` column. Try
>   it if many operations would be lost.
> * Operations that caught none of the selected species are removed (their
>   composition is undefined). The log reports how many.

## Step 4. Nominal CPUE and catch composition

```{r nominal, fig.height = 5}
nom <- plot_nominal(x)
nom$index
nom$positive
nom$distribution
```

Changes in the catch composition over time are the first sign of changes in
targeting (or of diverging trends in abundance).

```{r composition}
plot_composition(x)
plot_cpue_map(x, res = 5, by_period = 10)
```

In this simulation prices change in 2010 and 2020 while the true biomass
trends change smoothly, so shifts in composition at those years reflect
targeting.

> **Guidelines**
>
> * Look at the **proportion of positive operations** and the **distribution
>   of positive CPUE**; they guide the choice between a Tweedie and a delta
>   model.
> * Abrupt changes in composition that coincide with regulations, markets
>   or gear changes point to targeting, not abundance.

## Step 5. Habitat-based subsetting (optional)

Stephens and MacCall (2004) keep only operations likely to have fished the
habitat of the target species, predicted from the presence of the other
species with a binomial GLM. It is useful when many operations could not have
caught the species at all.

Here we use the least frequently recorded of the selected species as the
target.

```{r sm}
target <- tail(species_curve(x)$label, 1)
target
sm <- stephens_maccall(x, target = target)
sm
pl <- plot_stephens_maccall(sm)
pl$threshold
pl$coefficients
pl$kept
```

```{r sm-subset}
x_sm <- subset(x, sm$selected, note = paste("Stephens & MacCall,", target))
```

> **Guidelines**
>
> * Subsetting **discards data** and the subset is defined by the catch
>   composition, so it can change over time with targeting. Check
>   `pl$kept`: a trend in the share of operations kept means the subset
>   itself follows the targeting.
> * It removes operations that could not catch the species, but it does not
>   correct for different catch rates **within** the habitat; combine it
>   with the targeting covariates of Step 6 when targeting varies within the
>   subset.
> * The method needs a target that is absent from a fair share of
>   operations.

## Step 6. Targeting covariates from the catch composition

`add_targeting()` converts the catch of each operation into proportions,
applies a fourth-root transformation, runs a PCA, chooses the number of
non-trivial PCs with the Optimal Coordinates test, and clusters the operations
into fishing tactics (`n_pc + 1` clusters, Winker et al. 2014). It adds the PC
scores (`PC1`, `PC2`, ...) and the tactic (`FT`) to the data.

```{r targeting}
x <- add_targeting(x)
x$targeting$pca
```

### How many PCs?

```{r oc}
plot_oc(x)
```

A PC is retained if its eigenvalue is at least 1 (Kaiser-Guttman) and at
least the value predicted by the line through the following eigenvalues
(Optimal Coordinates). Model selection criteria such as AIC are not suitable
here: the composition is not independent of the response, so they select too
many PCs.

### What do the PCs and clusters represent?

```{r biplot, fig.height = 6}
plot_biplot(x)
tactic_table(x)
```

```{r targeting-plots, fig.height = 5}
pt <- plot_targeting(x, by = "month", res = 5)
pt$loadings
pt$composition
pt$time
pt$pcs_time
pt$map_tactic
```

In this simulation, the shares of tactics should change at the price
changes of 2010 and 2020: this is the targeting that the nominal CPUE cannot
separate from abundance.

> **Guidelines**
>
> * **Name the tactics** after their dominant species (`tactic_table()`) and
>   check that they make sense to people who know the fishery.
> * **Continuous or discrete?** PCs (DPC) represent gradual changes and
>   mixtures of tactics; clusters assume distinct tactics. Gear-defined
>   tactics favour clusters; continuous habitat gradients favour PCs.
> * **Diverging abundance trends.** When the relative abundance of species
>   changes strongly, the catch composition also carries an abundance
>   signal. PCs can then absorb part of the true trend of the species whose
>   relative abundance changes most, and more PCs absorb more. Coarser
>   clusters are less affected. Compare indices with 1, 2 and `n_pc` PCs and
>   with clusters, and be cautious when they differ.
> * **Robustness.** Compare clustering methods (`method = "kmeans"`,
>   `"hclust"`), the input (`on = "pcs"`), and the species set from Step 3.
> * Run `add_targeting()` **after** the last subsetting step; PCs and
>   clusters depend on the operations included.

### Removing irrelevant tactics (optional)

When the fleet uses tactics that cannot catch the species of interest (e.g.
pelagic trolling for a demersal species), these trips add zeros that carry no
information on its abundance, and the PCA spends its first axes on the
contrast with them. `select_tactics()` removes clusters, by number or by
species content, and recomputes the PCA and clusters on the remaining
operations:

```{r select-tactics, eval = FALSE}
tactic_table(x)                       # identify the irrelevant cluster(s)
x2 <- select_tactics(x, drop = 3)     # by number, or by species content:
x2 <- select_tactics(x, species = c("Sp A", "Sp B"), min_share = 0.5)
tactic_table(x2)
```

> * Select on a **group of species**, not on the target alone: the clusters
>   come from the catch composition, so selecting clusters rich in the
>   target would select operations by their catch of the target.
> * Removed clusters should be clearly different, not borderline; check the
>   share removed by year and compare the index with and without the
>   selection.

## Step 7. Fit the standardisation models

`fit_cpue()` fits a set of GAMs (`mgcv`) to the CPUE of one species and
derives an annual index from each. Every model contains the year effect and
differs in how it accounts for targeting:

| Model | Terms besides year |
|---|---|
| `Base` | none |
| `Covariates` | `terms` you give (e.g. month, port, gear) |
| `DPC` | `terms` + smooths of the first `n_pc` PC scores |
| `Cluster` | `terms` + fishing tactic (`FT`) as a factor |
| `...+S` | adds a smoother of the reported positions |

`vessel = TRUE` adds a random vessel intercept to every model. The default
family is Tweedie, whose power parameter is estimated together with the
model (also with random vessel effects); delta-lognormal and delta-gamma
models are alternatives.

Here we standardise the deep-water species, whose price rises over the three
decades so that it is targeted more and more.

```{r fit}
fit <- fit_cpue(x, species = "Anoplopoma fimbria",
                models = c("Base", "Base+S", "Cluster+S", "DPC", "DPC+S"),
                vessel = TRUE, spatial_k = 30, engine = "bam")
fit
```

The index of each year is the expected CPUE averaged over a reference set of
covariate values: all observed positions, the fishing tactics (and months, if
in the model), the median PC scores, and the average vessel.

> **Guidelines**
>
> * Add terms **in a meaningful order** (e.g. Base, then space, then
>   targeting) so that the step plot shows what each term does.
> * For a seasonal effect use a cyclic smoother:
>   `terms = 's(month, bs = "cc", k = 8)', knots = list(month = c(0.5, 12.5))`.
> * Use `engine = "bam"` for large data sets (tens of thousands of
>   operations or more).
> * With many zeros and no clear separation between "could not catch" and
>   "did not catch", Tweedie is a natural default; compare with a delta
>   model (`family = "delta_lognormal"`) when the proportion of positive
>   catches changes strongly over time.

## Step 8. Diagnostics

### Model comparison

```{r diagnostics}
cpue_diagnostics(fit)
```

AIC and BIC are only comparable between models with the same data and
family, and they tend to favour many targeting terms because the catch
composition is not independent of the response. Use them together with the
following plots.

### Indices and steps

Since the data are simulated, the true biomass trend is known and can be
added to the plot.

```{r indices, fig.height = 5}
truth <- data.frame(year = 1:30 + 1999,
                    true = sim$Biomass$True[sim$Biomass$Sp == "Sp4"])
plot_index(fit, true = truth)
```

```{r step, fig.height = 9}
plot_step(fit)
```

### Influence

The influence of a term in a year is the ratio between the expected CPUE
given that year's covariate values and given the values of all years
(Bentley et al. 2012). A term that changes the index has influence that
changes over time.

```{r influence}
plot_influence(fit, "DPC+S")
```

### Residuals, effects and vessels

For Tweedie models the residuals are randomised quantile residuals, which
should be standard normal if the model fits, despite the zeros.

```{r residuals, fig.height = 4}
res <- plot_residuals(fit, "DPC+S", by = "FT")
res$qq
res$fitted
res$year
```

### Predictive check of the error distribution

Residuals check each observation against its own predicted distribution.
A predictive check asks whether the fitted error distribution reproduces
the data **as a whole**: the share of zeros, the upper tail and the shape of
the positive catches. `ppc_cpue()` simulates replicate data sets from each
fitted model: for every operation the expected CPUE is computed (with the
coefficients drawn from their approximate posterior, so parameter
uncertainty is included) and a new value is drawn from the fitted Tweedie
(or delta) distribution. Each replicate is summarised like the observed
data, and the spread across replicates gives a 95% envelope. This is an
approximate posterior predictive check (Gelman et al. 1996); with
`param_uncertainty = FALSE` it is a plug-in (parametric bootstrap) check.

```{r ppc, fig.height = 5}
pp <- ppc_cpue(fit, models = c("Base", "DPC+S"), n_sim = 50)
pp
pp$plots$distribution
pp$plots$zeros
pp$plots$tails
```

How to read it:

* **Zeros** observed above the envelope: the model predicts too few zeros
  (e.g. a Poisson-type or lognormal-type error); below: too many.
* **Upper quantiles** observed below the envelope: the model's tail is too
  heavy (e.g. a negative binomial or a lognormal with large variance);
  above: too light.
* **Histogram of CPUE, zeros included** (`pp$plots$distribution`): observed
  bars (the zero bar leftmost) against the predictive distribution (points
  with 95% ranges). Systematic departures point to the wrong distribution or
  power parameter. CPUE spans orders of magnitude, so the histogram is on a
  fourth-root scale by default (zero stays next to the small catches);
  `scale = "log"` shows the shape of the positive catches best, with the
  zeros as a separate bar, and `scale = "raw"` the untransformed values.
* **By year.** Pooled over all operations, errors in opposite directions can
  cancel. The year effect fixes the mean CPUE of each year, not its zeros:
  a model without targeting covariates gives every operation in a year the
  same expected catch and gets the zeros wrong in years when the fleet
  fished elsewhere. A good targeting term brings the observed zeros inside
  the envelope. Other groupings: `by = "FT"`, `by = "month"`.
* `tail_prob` is the share of replicates at least as large as the observed
  value; values near 0 or 1 flag misfit.
* Information criteria cannot compare a Tweedie with a delta model; a
  predictive check can.

```{r effects, fig.height = 4}
eff <- plot_effects(fit, "DPC+S")
eff$numeric
eff$space
```

```{r vessel-effects, fig.height = 4}
plot_vessel_effects(fit, "DPC+S")
```

> **Guidelines**
>
> * Check the **influence** of the targeting terms: a strong trend in
>   influence means targeting changed and the term corrects for it.
> * Check that residuals show **no trend by year** and no pattern against
>   the targeting covariates or tactics.
> * Check that the error distribution reproduces the **zeros and the upper
>   tail by year** (`ppc_cpue()`).
> * Compare the **DPC and Cluster** indices. If they disagree, look at what
>   the PCs capture over time (Step 6) and whether the species' relative
>   abundances may have diverged.

## Step 9. The final index

```{r index}
cpue_index(fit, "DPC+S")
```

```{r export, eval = FALSE}
export_index(fit, "cpue_index_DPC.csv", model = "DPC+S")
```

`index` is in the units of the data (CPUE at the reference covariate
values), `relative` is scaled to its mean, and `se_log` is the standard
error of the log index, the usual input to stock assessment models.

## Processing log

The log records every step, how many operations it removed, and the
settings used.

```{r log}
x$log
```

## Reporting checklist

| Item | Report |
|---|---|
| Operation and effort | unit (set, trip, ...), effort measure |
| Data selection | all subsetting rules and the processing log |
| Species set | threshold, criterion (occurrence or weight), species kept |
| Fleet | indicator-fleet rules (activity and target species), number of vessels |
| Subsetting | Stephens & MacCall target, critical value, share kept over time; tactics removed (`select_tactics()`) |
| Targeting | transformation, eigenvalues, OC result, number of PCs and clusters, clustering method, tactic composition |
| Models | family, terms of each model, random effects, reference set of the index |
| Diagnostics | model table, step and influence plots, residuals, predictive check (zeros and tails by year), vessel effects |
| Index | table with nominal and standardised CPUE, log SE and CV |
| Sensitivity | indices with alternative species sets, numbers of PCs, clusters, families |

## References

Bentley, N., Kendrick, T.H., Starr, P.J. and Breen, P.A. (2012) Influence
plots and metrics: tools for better understanding fisheries
catch-per-unit-effort standardizations. ICES Journal of Marine Science 69:
84-88.

Dunn, P.K. and Smyth, G.K. (1996) Randomized quantile residuals. Journal of
Computational and Graphical Statistics 5: 236-244.

Gelman, A., Meng, X.-L. and Stern, H. (1996) Posterior predictive assessment
of model fitness via realized discrepancies. Statistica Sinica 6: 733-807.

Hoyle, S.D., Langley, A.D. and Campbell, R.A. (2014) Recommended approaches
for standardizing CPUE data from pelagic fisheries. WCPFC-SC10-2014/SA-IP-10.

Stephens, A. and MacCall, A. (2004) A multispecies approach to subsetting
logbook data for purposes of estimating CPUE. Fisheries Research 70: 299-310.

Winker, H., Kerwath, S.E. and Attwood, C.G. (2013) Comparison of two
approaches to standardize catch-per-unit-effort for targeting behaviour in a
multispecies hand-line fishery. Fisheries Research 139: 118-131.

Winker, H., Kerwath, S.E. and Attwood, C.G. (2014) Proof of concept for a
novel procedure to standardize multispecies catch and effort data. Fisheries
Research 155: 149-159.