knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

1. Introduction

The catchfunction package (which we refer to as the ABC To TAC and Commercial Harvest, aka ATTACH, model: R package rename forthcoming) was created for the Alaska Climate Integrated Modeling Project (ACLIM) by Amanda Faig (University of Washington; School of Aquatic Fisheries and Sciences) and Alan Haynie (NOAA; NMFS). This function, in a nutshell, takes Bering Sea (BS) acceptable biological catch (ABC) as input and uses a series of regression estimates to predict total allowable catch (TAC) and from that the commercial harvest in the Bering Sea, based on ABC, TAC, and catch data from 1992 to 2017.

If you have yet to install the R package (or if you want to upate your package, e.g. to see if a newer version has been released) run the following code:

install.packages("devtools")
library("devtools")
install_github("amandafaig/catchfunction")
library("catchfunction")

Otherwise, simply load the library:

library(catchfunction)

2. Predicting TAC

We created ATTACH as a two-step model. In the first step, TAC is predicted from ABC. The user passes, as inputs, the Bering Sea ABCs into the model for as many species as are defined in their biological models; up to the 22 species under the BSAI 2 million ton ecosystem cap ("the ecosystem cap"). To see the full list of species included in the cap, see the list of arguments in the help file, or in the pdf manual.

?catch_function

Since the ecosystem cap is for the entire BSAI, but the ABC input is only for the Bering Sea, the ATTACH model first calculates BSAI ABCs from the BS ABC inputs. (The ABC input is only for the Bering Sea is this is the scope of the ACLIM project.) We calculate BSAI ABCs assuming that the Aleutian Island ABC rises and falls (relative to the historical mean) in proportion to how the (user defined) Bering Sea ABC compares to it's historical mean. For ABCs left undefined by the user, we assume those species ABCs are at their historical means.

The entire set of BSAI ABCs is then passed into the first step of the model: estimating TAC. Each year the North Pacific Fishery Management Council ("the Council") sets the TAC of individual stocks based on the ABC estimates for the individual stocks. The ecosystem cap mandates that the Council must ensure the sum of these TACs does not exceed 2 million metric tons. The stock assessment ABCs and the TACs set by the Council are published annually in the Alaska Region's groundfish harvest specifications. We used data from 1992 to 2017 to create the current version of the model (version 1).

TACs for each stock were estimated statistically using a log-linear model. For $j = 1, 2, ... n$ stocks, the general model for stock $i$ took the form:

$$ \ln(TAC_{i,t}) = \alpha_{i} + \beta_{i}\ln(ABC_{i,t}) + \sum_{j \neq i}^n \beta_{j} ABC_{j,t} + \sum_{k=1}^m \beta_{k}I_{k,t} + \varepsilon_{i,t} $$

where $\alpha_{i}$ is the stock-specific intercept for species $i$, $\beta_i$ is the elasticity of the TAC of species $i$ with respect to its own ABC, and $\beta_j$ relates the ABC of some species $j \neq i$ to the TAC of species $i$. The effect ($\beta_k$) of $k =1, 2, ... m$ events or policy changes (e.g., changes in management, area closures, or implementation of catch share programs) on TAC was also estimated where $I_{k,t}$ is an indicator variable for event $k$ in year $t$, and $\varepsilon_{i,t}$ denotes the residual error for the prediction in year $t$. How the errors handled is discussed more in section 4. To see each individual equation see Haynie et al. (in prep).

ATTACH uses the predicted coefficients to predict TAC from ABCs. The events/policies are assumed to reflect the last year of the dataset. So, for example, the Amendment 80 indicator variable is set to 1 in predictions, while the Steller Sea Lion closure of 2011 to 2014 is set to 0 (and the Steller Sea Lion limited reopening is set to 1).

When the model is passed a set of historical ABCs, the predicted TACs add up to less than 2 million metric tons, since necessarily these combinations led to a net TAC at or below the ecosystem cap. (Even this is not guaranteed, however, due to prediction error.) When the set of ABCs input into ATTACH is not a historical set, it is possible the predicted TACs based on the regression estimates alone could together exceed the ecosystem cap. To ensure ATTACH does not violate the ecosystem cap, we check the sum of the TACs and, if they exceed 2 million metric tons, we decrease all TACs proportionally, except for that of BS and AI Sablefish, BSAI Shortraker rockfish, and BSAI Northern rockfish.

3. Predicting Catch

The output from the TAC prediction step is then passed to another sub-model, the catch prediction step. Catch estimates are based on TAC data (from the aforementioned groundfish harvest specifications) and catch data from the Catch Accounting System. As in the first step, we use data from 1992 to 2017.

Catches for each stock were estimated statistically using a log-linear model. For $j = 1, 2, ... n$ stocks, the general model for stock $i$ took the form:

$$ \ln(catch_{i,t}) = \alpha_{i} + \beta_{i}\ln(TAC_{i,t}) + \sum_{j \neq i}^n \beta_{j} TAC_{j,t} + \sum_{k=1}^m \beta_{k}I_{k,t} + \varepsilon_{i,t} $$

where $\alpha_{i}$ is the stock-specific intercept for species $i$, $\beta_i$ is the elasticity of the catch of species $i$ with respect to its own TAC, and $\beta_j$ relates the TAC of some species $j \neq i$ to the catch of species $i$. The effect ($\beta_k$) of $k =1, 2, ... m$ events or policy changes (e.g., changes in management, area closures, or implementation of catch share programs) on TAC was also estimated where $I_{k,t}$ is an indicator variable for event $k$ in year $t$, and $\varepsilon_{i,t}$ denotes the residual error for the prediction in year $t$. How the errors handled is discussed more in section 4. To see each individual equation see Haynie et al. (in prep).

The catch estimate for a given species can exceed it's own TAC. This is because the TAC measure we use is the TAC set at the beginning of the season, and in-season management can adjust TAC to an extent. In ATTACH we check that catch does not exceed the BSAI wide ABC and that the ecosystem cap, but otherwise allow predicted catch to exceed TAC.

4. The Ensemble

ATTACH is an ensemble of three models that include different explanatory variables that fit data better for different species. We chose the ensemble rather than any individual model because it better captures possible environmental and policy uncertainty and is therefore more likely to be robust to ABC combinations and individual ABC levels outside of historical bounds. The three models in the ensemble differ in the error structures in both the TAC and catch estimation equations. In all of the models, the errors in the TAC estimation stage are linked via Seemingly Unrelated Regressions (SUR; Zellner 1962), a common econometric modeling technique. The format of the SUR here is a set of linear regression equations that are valid independently for each species; exogenous shocks which affect one species are then assumed to affect all the included species, implying the error terms are correlated across species. This error covariance is especially strong in the TAC stage, where the 2 million ton cap ensures the errors of all species sum to zero. In the catch estimation stage, model 1 assumes that each log-linear regression is independent, model 2 has two groups of SUR-linked regressions (representing species that are typically caught concurrently), and model 3 includes a group of species that are caught concurrently.

The ensemble averages the estimated catch of the three models equally before returning the estimated Bering Sea catch to the user. Only the estimates for the species whose ABCs were specified by the user are returned.

5. Performance

As a rule, it is easier to predict TAC and catch in directed fisheries (e.g. pollock, Pacific cod, yellowfin sole) than in bycatch fisheries. In particular in the catch prediction stage, there is a large degree of randomness in some of the bycatch fisheries. Despite this, in all fisheries using the ATTACH estimate provides a better estimate of catch than the standard assumption made in biological modeling, which is to set catch equal to ABC. In some fisheries outside of the BSAI this may not be a terrible assumption but in the BSAI, for most species managed under the ecosystem cap, this would result in the chronic overestimation of harvest.

The following figure illustrates how the ATTACH predictions compare to setting catch equal to ABC in selected species, created originally for Reum et al. (in prep).

Comparison of observed catches in the Bering Sea-Aluetian Islands against predicted catches from the ATTACH fishery submodel (orange circles) and Allowable Biological Catch (ABC) values from stock assessments (blue symbols).  ATTACH accounts for policy considerations in the setting of catch limits and provides a more accurate representation of realized catches relative to information on ABC alone.  For each species, Pearson’s correlation coefficient (R) and the proportional bias (pb) are provided. The dashed black line corresponds to the 1:1 relationship and linear trend lines (orange and blue) are provided for each relationship. Data and predictions span the period 1992 to 2017.

For each stock, performance of the models was evaluated using leave-one-out cross-validation (LOO-CV). We estimated the coefficients of each model using all but one year in the data-set, and then used the estimators to predict the TAC and catch of the omitted year. We followed this process for each year and saved the predicted catches of that year. We calculated the difference between the predicted and actual catch for each year (1992-2017) and used this to create a variety of metrics: a simple sum of differences; a sum of percent differences; sums of squared differences, weighted by value, ABC, TAC, or catch; and the sum of squared percent differences. We calculated these same metrics for in-sample predictions as well (that is, predictions were made for each year based on the complete sample of years).

6. Scenarios

Whitefish Preference scenario (Scenario 2)

In the "whitefish preference"" scenario the Council allocates an extra 10% TAC to whitefish, decreasing flatfish TAC by up to 50% in order to do so and still comply with the ecosystem cap. (Whitefish are defined as Alaska pollock and Pacific cod; flatfish are defined as yellowfin sole, northern rock sole, Kamchatka flounder, Alaska plaice, Greenland turbot, Other flatfish, flathead sole, and arrowtooth flounder.)

We model this by first estimating TAC using the base model. We then define $x$ as minimum of: 10% whitefish TAC; 50% flatfish TAC; and the difference between whitefish ABC and TAC.

We increase pollock and Pacific cod TAC by raising both TACs equally until either the net increase in whitefish TAC is $x$, or one of the whitefish species' TAC is equal to it's respective ABC. If the latter occurs, we then increase only the species whose TAC is less than it's ABC until the net increase in whitefish TAC is equal to $x$.

We decrease flatfish TAC by $x$, spreading the decrease proportionally across the flatfish species.

We then pass these adjusted TACs to the catch portion of the model to estimate catch.

Flatfish Preference scenario (Scenario 3)

In the "flatfish preference" scenario the Council decreases whitefish TAC by up to 10%, and increases flatfish TAC by that amount. (Whitefish and flatfish defined as in scenario 2.)

We model this by first estimating TAC using the base model. We then define $y$ as the minimum of: 10% whitefish TAC; the difference between flatfish ABC and TAC; and the difference between whitefish TAC and the 'threshold'. (Below a threshold of approximately 1.5 million tons, we have found that the Council sets TAC equal to ABC in whitefish. For this scenario we adjust this threshold down by 10%, to 1.35 million tons.)

We increase all flatfish proportionally until either the net increase in flatfish TAC is equal to $y$, or one of the flatfish species' TAC is equal to it's ABC. If the latter occurs, we remove that species from the set and repeat the pattern until the net increase in TAC is equal to $y$.

We decrease whitefish TAC proportionally so that the net decrease in whitefish TAC is $y$.

We then pass the adjusted TACs to the catch portion of the model to estimate catch. Unlike in scenario 2, where we assume catch behaves as usual relative to TAC, in the flatfish preference scenario we assume the fleet is better able to catch flatfish TAC. We make this assumption because currently most flatfish species' catches are rarely limited by TAC. We assume catch is up to 30% higher for each flatfish species than it would normally be for any given TAC; still respecting the rule that catch for any given species cannot exceed it's respective ABC and net catch cannot exceed the ecosystem cap.

Leave one or more out of the ecosystem cap scenarios (Scenarios 5.X)

Scenarios 5.1 and 5.2 are variations of the answer to a question asked of us: what if we remove a species (or more) from the cap?

The answer to that question was broken out into parts because we felt the answer depended on what assumptions we made.

In scenario 5.1 we assume the ABC still influences the TAC of other species, even if the species is not managed under the cap. This means we assume the tons that would otherwise have been allocated to the now removed species are not reallocated to other species. i.e. the ecosystem cap would be effectively reduced from 2 million metric tons, in order to keep the ecosystem cap constraints on other species the same as before.

In scenario 5.2 we assume the cap remains 2 million metric tons, but the council only uses this extra leeway when ABCs are unusually high. Thus we assume that other species TAC calculations are as before, despite no cap in place. (Meaning that in most years the net TAC will be below 2 million metric tons, since the TAC calculation is based on Council decisions made with respect to the cap.) The main distinction between 5.1 and 5.2 comes into play in situations where the sum of the predicted TACs would have surpassed 2 million metric tons.

In both we apply the code described at the end of section 2, which would adjust TACs down to comply with the ecosystem cap, but in 5.1 we include what the excluded species' TAC would have been in that calculation, whereas in 5.2 we do not.

In all of the above we assume TAC = ABC for the removed species, and catch is some assumed fraction based on the multiplier provided by the user.

7. Examples

Example 1: Status Quo Catch

If you are interested in predicting what catch would be, assuming everything about the world stays as it was in 2017 (A80 and AFA alive and well, Steller Sea Lion closure partially reopened, etc...). You would simply use:

AP_BS_ABC = 2e6
ATF_BS_ABC = 2e5
YFS_BS_ABC = 2e5
catch_function(scenario = 1, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

In the above example we defined the ABCs before calling catch_function simply to illustrate that the numbers chosen were not significant. We will use these numbers for all subsequent examples for consistency. Scenario 1 is the status quo scenario. Scenarios 2+ are defined in the help file, and correspond to alternative scenarios like the whitefish and flatfish preference scenarios described in the previous section.

Example 2: Status Quo TAC

If for your purposes you only want the predicted TAC allocation for these species, you would simply use:

TAC_function(scenario = 1, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

Note, TAC_function is based on annual allocations and thus returns BSAI TAC for all species except Alaska Pollock, Greenland turbot, Pacific oeach perch, Sablefish and Other rockfish. This is because those are the species for which TAC is set on the BS level*.

*Pacific cod is also set on a BS level, but this is a more recent development, and so for now we still model TAC as being set on a BSAI level in Pacific cod.

Example 3: Whitefish Preference Catch (& TAC)

# Status Quo catch prediction (same as example 1)
catch_function(scenario = 1, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

# Whitefish preference catch prediction
catch_function(scenario = 2, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

As you can see above, pollock catch is now higher than it was in Example 1, and yellowfin sole catch is lower. It may seem strange that arrowtooth flounder catch is unchanged: arrowotooth flounder catch, historically, has not responded to changes in it's own TAC. As a result, while Arrowtooth flounder TAC is adjusted in the whitefish preference scenario as expected, catch is unchanged.

# Status quo TAC prediction (same as example 2)
 TAC_function(scenario = 1, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

# Whitefish preference TAC prediction
TAC_function(scenario = 2, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

Example 4: Series 5 scenarios

In the "5.x" Scenarios, we look at what would happen if a species (or more) left the ecosystem cap. In the following examples we assume arrowtooth and yellowfin are allowed to leave the cap.

We include a multiplier which assumes we catch 50% of the arrowtooth TAC/ABC (we assume any species not managed under the ecosystem cap has TAC set at ABC, as that is what happens in most fisheries not managed with a similar sort of cap) and 100% of the Yellowfin TAC/ABC is caught.

# Status Quo catch prediction (same as example 1)
catch_function(scenario = 1, 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

# Remove Arrowtooth and Yellowfin from the cap (in that TAC can equal ABC for them)
# assume 50% of Arrowtooth ABC/TAC is harvested, 100% of Yellowfin ABC/TAC.
# Assume all other species as affected by the cap as before.
catch_function(5.1, spptomult = c("Arrowtooth","Yellowfin"), multiplier = c(0.5,1), 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

# Remove Arrowtooth and Yellowfin from the cap; 
# assume 50% of Arrowtooth ABC/TAC is harvested, 100% of Yellowfin ABC/TAC.
catch_function(5.2, spptomult = c("Arrowtooth","Yellowfin"), multiplier = c(0.5,1), 
               Pollock = AP_BS_ABC, 
               Arrowtooth = ATF_BS_ABC, 
               Yellowfin = YFS_BS_ABC)

As you can see, in 5.1 Pollock's estimated catch is the same as in the status quo scenario: only when the predicted catch would have been lowered due to the 2 million ton catch does the new "room" under the ecosystem cap come into play. In 5.2 however, Pollock's estimated catch increases as the new "room" under the ecosystem cap is always utilized.

8. Future Steps

Manuscript for the work that laid the foundation for the ATTACH package

Alan Haynie, Amanda Faig, Kirstin Holsman, Anne Hollowed, Jonathan Reum, Steve Kasperski, and Mary Furuness. "Predicting future management allocations and catch under the Bering Sea and Aleutian Islands Ecosystem Cap." In Prep: contact alan.haynie@noaa.gov.

Publications using the ATTACH (aka catchfunction) package

  1. Jonathan Reum, Julia Blanchard, Kirstin Holsman, Kerim Aydin, Anne Babcock Hollowed, Albert J Hermann, Wei Cheng, Amanda Faig, Alan C Haynie, and Andre Punt. "Ensemble projections of future climate change impacts on the Eastern Bering Sea food web using a multispecies size spectrum model. Frontiers in Marine Science." Accepted: contact jonathan.reum@noaa.gov.

  2. Kirstin Holsman, Alan C Haynie, Anne Babcock Hollowed, Jonathan Reum, Kerim Aydin, Wei Cheng, Amanda Faig, Jim Ianelli, Kelly Kearney, Andre E Punt. "Ecosystem-based fisheries management forestalls climate-driven collapse." Submitted: contact kristin.holsman@noaa.gov



amandafaig/catchfunction documentation built on March 17, 2023, 11:52 p.m.