Input data (sardine and hake stocks and environmental variables)

Description

Assessment summary table for the hake stock in ICES subareas 4, 6, and 7, and ICES divisions 3.a, 8.a, 8.b, and 8.d (northern hake stock), as well as for the sardine stocks: northern stock (ICES divisions 8.a, 8.b, and 8.d) and southern stock (ICES divisions 8.c and 9.a).

The data corresponds to the assessment conducted in 2021 and was derived using the icesSAG package. Environmental data is also reported, including the Atlantic Multi-decadal Oscillation (AMO) and the North Atlantic Oscillation (NAO) indices.

Usage

data(knobi_dataset)

Format

A list where items correspond to the assessment summary tables for the ICES northern hake stock, the ICES southern and northern sardine stocks, and the environmental data frame.

Details

The stock's data can be obtained using the following R code:

# Hake northern stock data:  
# hake_n <- icesSAG::getSAG(stock = "hke.27.3a46-8abd", year = 2021)
# hake_n <- hake_n[-nrow(hake_n),]  

# Sardine southern stock data:  
# sardine_s <- icesSAG::getSAG(stock = "pil.27.8c9a", year = 2021)[1:43,]  

# Sardine northern stock data:  
# sardine_n <- icesSAG::getSAG(stock = "pil.27.8abd", year = 2021)[1:21,]  

The AMO data, i.e. the Atlantic Multi-decadal Oscillation index, and the NAO data, i.e. the North Atlantic Oscillation index — both large-scale climate variability patterns that affect the entire North Atlantic Ocean — were downloaded from the UCAR website: https://climatedataguide.ucar.edu/.

Source

getSAG function of the icesSAG package and the UCAR website: https://climatedataguide.ucar.edu/.

KBPM environmental analysis

Description

Analyze and model the relationships between surplus production (SP) and environmental covariable(s) to test whether productivity changes in response to environmental fluctuations. The analysis is conducted in three steps:

(1) Correlation Analysis: Assess the correlation between the standardized environmental variable(s), at different delays (lags), and the KBPM residuals using Pearson's correlation or autoregressive models. See details.

(2) Variable Selection: Determine which lagged environmental variable(s) will be included in the environmental KBPM models.

(3) Environmental Fitting: Fit the KBPM model incorporating environmental effects as both additive and multiplicative effects (see details).

Arguments

Details

It is important to mention that the environmental variable(s), in a first step, are standardized, in order to make their scale and magnitude comparable. To do this, each variable is subtracted from its mean and divided by its standard deviation.

The additive environmental model adds the following term on the right-hand side of Eq. (1) or Eq. (2) described in knobi_fit: c X_{t - lag} B_t, being X_{t - lag} the environmental variable at time t - lag and B_t the biomass or SSB at time t.

The multiplicative environmental model multiplies the right-hand side of Eq. (1) or Eq. (2) by \exp(c X_{t - lag}).

In the case of these models, the estimated biological reference points (BRPs) correspond to a value of the scaled environmental variable equal to the mean of the time series, i.e., X_t = 0, which cancels out the effect of the parameter c. The estimates of the remaining parameters included in Eq. (1) or Eq. (2), and therefore for the BRPs as well, will differ from the base model due to the inclusion of environmental effects. For more details, such as the calculation of BRPs as a function of the environmental variable, see the vignettes.

If ar_cor = TRUE, the correlation analysis between the knobi_fit residuals and the environmental variable(s) is conducted as follows:

First, an AR model is fitted to the KBPM base residuals:

r_t = \sum_{i = 1}^{\rho} \beta_i r_{t - i} + \epsilon_t

where r_t is the KBPM base residual for year t, and \rho is the AR model order, estimated as the maximum lag at which the absolute value of the residuals' partial autocorrelation exceeds qnorm(0.975) / \sqrt{N_r}, with N_r being the length of the residuals series.

AR models are then fitted to the residuals incorporating each lagged environmental variable X_{t - lag} as an explanatory covariate:

r_t = \sum_{i = 1}^{\rho} \beta_i r_{t - i} + X_{t - lag} + \epsilon_t

for lag = 0, 1, \ldots, nlag. Then, we have one autoregressive model for each lag.

The lagged environmental variable whose AR model has the lowest AIC is selected for inclusion in the environmental KBPM fit. If none improve the AIC compared to the base AR model, no environmental covariate is added.

It is important to highlight that the results include the analysis of the AR model only with the base model residuals in order to determine the need for coupling environmental information, considering that it would not be necessary if this model shows a lower AIC, even reducing the number of parameters to fit.

Value

A list containing the results of the three-step environmental analysis:

• add: estimates of the additive model parameters.

• mult: estimates of the multiplicative model parameters.

• BRPs: reference points (BRPs) estimates for each model (see details).

• df: data frame with the information used in the fit.

• selected_var: environmental variable(s) used in the fit.

• selected_lag: data frame providing the time lag of the environmental variable(s) in the KBPM fit.

• lag_cor: correlation between the environmental variable(s) and the KBPM residuals for each time lag.

• env_aic: if ar_cor = TRUE, AIC values of each autoregressive model (see details).

• scaled_var: standardized environmental variable(s) used in the fit.

• plots3D: list with 3D plot objects (if plot3d = TRUE).

• residuals: Pearson residuals from each model fit (base, additive, and multiplicative).

• performance_metrics: array of performance and accuracy measures for each model, including those from the error_table output of knobi_fit, plus an F-test comparing environmental models to the base model.

The function also returns plots in the graphics window and saves them (if plot_out = TRUE) to the specified directory. The first plot shows the correlation analysis between the environmental variable(s) and the base residuals. The second shows fitted SP values from all models. If multicovar = FALSE and plot3d = TRUE, 3D plots are also returned. If multicovar = TRUE, a plot with Pearson correlations among environmental variables is included.

Author(s)

• Anxo Paz

• Marta Cousido-Rocha

• Santiago Cerviño López

• M. Grazia Pennino

Examples

# First, run the example of knobi_fit function

# Then, provide environmental data series

Env <- knobi_dataset$Env

# The environmental data series must start in the first year of the KBPM fit data
# minus the provided nlag or lag arguments
years <- knobi_results$df$Year # See knobi_fit example to obtain the knobi_results object
ind <- which(Env[,1]==years[1])
ind1 <- which(Env[,1]==years[length(years)])
nlag <- 5
Env <- Env[(ind-nlag):ind1,]

# Now we create the environmental list
data <- list(env=data.frame(AMO=Env$AMO,NAO=Env$NAO),
           years=Env$years)
control <- list(nlag=nlag)

knobi_environmental <- knobi_env(knobi_results,data,control)
knobi_environmental
knobi_environmental$BRPs  # use the '$' to access to all the fit information

KBPM fit

Description

This function, that is the main function of the knobi package, fits a type of surplus production models named known-biomass production models (KBPM) (MacCall, 2002). The surplus production curve is fitted using the observed catch time series and the biomass or SSB (Spawning Stock Biomass) estimates derived from the fit of another stock assessment model.

Arguments

Details

The KBPMs implemented in the current package are explained below.

Schaefer model (1954) (Eq. 1):

SP_{t} = r \overline{B}_{t} (1-(\overline{B}_{t}/K))

where SP_{t} is the surplus production, \overline{B}_{t} is the average biomass or SSB (mean of two consecutive years), r is the population growth rate parameter, and K is the carrying capacity. The subscript t denotes the time (years).

Pella and Tomlinson model (1969) (Eq. 2):

SP_{t} = (r/p) \overline{B}_{t} (1-(\overline{B}_{t}/K)^{p})

where SP_{t} is the surplus production, \overline{B}_{t} is the average biomass or SSB, r is the population growth rate parameter, K is the carrying capacity and p is the asymmetry parameter. The subscript t denotes the time (years).

It should be noted that, in order to associate the average biomass or SSB with the catch series, the SSB time series must be one year longer than the catch series. If both time series have the same length, the last year of catch data will be excluded, and a warning will be displayed in the R console.

The recruitment values and the F_input estimates are included to get an overview of the stock status, but are not needed for the KBPM fitting. Similarly, the input reference points are only used for comparison purposes.

KBPMs have also proven their usefulness for the multispecies management objectives. The KBPM approach can be applied to analyze the joint dynamics of the targeted fish species within a community, using the aggregated biomass and catch data, defining in this way a simple data-limited ecosystem model to assess the ecosystem status (see example).

Value

The results of the KBPM fit include the following:

• params: model parameters estimates.

• df: the data used for the model.

• BRPs: biological reference points estimates:

  • K: carrying capacity.

  • MSY: maximum sustainable yield (MSY).

  • B_MSY: biomass at MSY.

  • F_MSY: fishing mortality at MSY.

  • MSYoverK: ratio of MSY and K.

• residuals: Pearson's residuals calculated as (observations-estimates)/sd(observations).

• performance_metrics: data frame of accuracy and model performance measures:

  • SER: standard error of the regression, calculated as the root of the rate between the sum of residual squares and the degrees of freedom of the regression.

  • R-squared: coefficient of determination.

  • adj-R-squared: adjusted coefficient of determination.

  • AIC: Akaike information criterion.

  • RMSE: root mean squared error (observed vs. estimated values).

  • MAPE: mean absolute percentage error (observed vs. estimated values).

• input: the input list updated with the annual average biomass ($input$Average_Biomass), the surplus production ($input$SP) and the F estimates derived from the fit ($input$F_output).

• control: the updated control settings of the fit.

• optimx: list of some results provided by optimx:

  • value: value of the objective function in the minimization process.

  • convergence: integer code. ‘0’ indicates successful completion in the optimization.

  • message: character string giving any additional information returned by the optimizer, or NULL.

The plots are displayed in the plot window and also saved (if plot_out=TRUE) in the provided directory or in the current one (if plot_dir is not provided). The following input quantities are plotted: fishing mortality time series, SSB or biomass, surplus production and catch time series. Plots of catch over fishing mortality, fishing mortality over SSB (or biomass), and catch over SSB (or biomass) time series with a smooth line from a "loess" regression are also available. Plot of input-output time series of fishing mortality are also provided and includes horizontal lines for both, input and output, fishing mortalities at MSY (one line if input F_MSY is NULL). Fishing mortality relative to F_MSY is also plotted including an horizontal line at one (note that values greater than 1 indicate an F greater than F_MSY). The analogous SSB (or biomass) plots are also reported. On the other hand, the fitted surplus production curve is plotted twice with the SSB (or biomass) and SP observations (first) and with the catch and SP observations (second). Finally, a plot with the KBPM residuals is reported.

Author(s)

• Anxo Paz

• Marta Cousido-Rocha

• Santiago Cerviño López

• Maria Grazia Pennino

References

Schaefer, M.B. (1954). Some Aspects of the Dynamics of Populations Important to the Management of the Commercial Marine Fisheries. Bulletin of the Inter-American Tropical Tuna Commission. 1:26-56.

Pella, J.J., Tomlinson, P.K. (1969). A generalized stock-production model. Bulletin of the Inter-American Tropical Tuna Commission. 13:421–58.

MacCall, A. (2002). Use of Known-Biomass Production Models to Determine Productivity of West Coast Groundfish Stocks. North American Journal of Fisheries Management, 22, 272-279.

Examples

library(knobi)

data(knobi_dataset)
hake_n <- knobi_dataset$hake_n


# Then, create the data object.

data<-list()
data$SSB<-hake_n$SSB # We take the SSB in our data.
data$Catch<-hake_n$catches # We take the catch in our data.
data$F_input<-hake_n$F # We take the F in our data.
# Reference points estimates from ICES stock assessment model:
# ICES. 2021. Working Group for the Bay of Biscay and the Iberian Waters Ecoregion
# (WGBIE). ICES Scientific Reports. 3:48.1101 pp.
data$RP<-list(F_MSY=0.259, B_MSY=207398, MSY=75052, K=NA)
# In this case, B_MSY corresponds to SSB_MSY, since control$method<-"SSB"
# (see control list below).
data$years<-hake_n$Year    # Years corresponding to the catch values
# (can be different than the years corresponding to SSB or biomass).


# Now we define the control.

control<-list(
  pella = "TRUE")   # Logical. TRUE for Pella-Tomlinson model.
                    #          FALSE for Schaefer model.

# Finally, we can fit the model

knobi_results<-knobi_fit(data,control,plot_out=TRUE,plot_filename="results")
# Note that a warning is shown for the reduction of the SSB vector 
# so that the length is the same as the catch length

knobi_results
knobi_results$BRPs  # use the '$' to access to all the fit information


## Fitting multispecific KBPM

# Below, a multistock approximation aggregating the
# northern and southern stocks of sardine is performed.

# Firstly, read southern stock data
sardine1 <- knobi_dataset$sardine_s

# Secondly, read northern stock data
sardine2 <- knobi_dataset$sardine_n

# Extract common years of data in both stocks
index <- which(sardine1$Year %in% sardine2$Year)
sardine1 <- sardine1[index,]

# Create a data.frame where the SSB and the catch are
# the sum of such data in the two stocks
years<-sardine1$Year
sardine <- data.frame(years=years,SSB=sardine1$SSB+sardine2$SSB,
                      catch=sardine1$catches+sardine2$catches)

# Once the total SSB and catch are available
# we follow previous KBPM illustration
data<-list()
data$SSB<-sardine$SSB
data$Catch<-sardine$catch
data$years<-sardine$years

knobi_results2<-knobi_fit(data)
# In this case, in addition, a series of warnings are shown due to 
# the SP being less than 0 for some year.
knobi_results2

KBPM projections

Description

This function forecasts population and fishery dynamics under a defined set of management targets. It projects the time series of biomass (or SSB) and the surplus production for a set of future catch or fishing mortality scenarios.

Arguments

Value

A list containing the projection results.

• df: data frame containing the projected time series. The rows correspond to the projection years, while the columns represent stock quantities: biomass or SSB, surplus production (SP), catch (C), and fishing mortality (F). Three additional columns specify the catch or fishing mortality scenario (Sc), the model used, and the environmental scenario (EnvSc).

• plots: list containing the plots with the projections. Each plot is a panel plot that includes four subplots representing biomass or SSB, surplus production, fishing mortality, and catches for each catch or fishing mortality scenario and for each environmental scenario (if they is provided).

The resulting plots are displayed in the plot window and are also saved (if plot_out="TRUE") in the provided directory or in the same directory as link{knobi_fit}.

Author(s)

• Anxo Paz

• Marta Cousido-Rocha

• Santiago Cerviño López

• M. Grazia Pennino

Examples

### Projecting through catch with no environmental information

# Then, create the data frame containing the selected catch for the projected
# years. In this illustration, within each scenario, the catch values are
# constant through the projected years. Three scenarios are considered:
# (i) catch value equal to the last historical catch multiplied by 1,
# (ii) last historical catch multiplied by 1.2 and
# (iii) last historical catch multiplied by 0.8.

catch<-rep(knobi_results$df$C[length(knobi_results$df$C)],5)

C<-data.frame(catch=catch,
              catch08=0.8*catch,
              catch12=1.2*catch)

# Then, knobi_proj function can be applied

knobi_proj(knobi_results, c=C)


### With environmental information

# In this case, in addition to the previous example, the 'knobi_env' example
# has to be run at first, where AMO variable was selected in the fit

# We include the future values of the environmental variable(s) in a data
# frame containing the environmental covariable values for the projected
# years. Three scenarios are considered:
# (i) Constant AMO equal to last year's AMO,
# (ii) Constant AMO equal to last year's AMO with a 50% increment
# (iii) Constant AMO equal to last year's AMO with a 50% decrease

last_AMO <- Env$AMO[length(Env$AMO)]
env <- data.frame( AMOi=rep(last_AMO,5),
                   AMOii=rep(last_AMO*1.5,5),
                   AMOiii=rep(last_AMO*0.5,5))

# Based on the previous objects we can apply the projection function.

knobi_proj(knobi_results, knobi_environmental, c=C, env=env)


### Through fishing mortality without environmental information

# Alternatively, projections can be based on fishing mortality.
# The scenarios presented below have been created from the estimated F_msy of
# knobi_fit analysis.

fmsy<-knobi_results$BRPs['F_MSY']
ff<-rep(fmsy,8)
f<-data.frame(f=ff,f12=ff*1.2,f08=ff*0.8)

knobi_proj(knobi_results, f=f)


### Through fishing mortality with environmental information

knobi_proj(knobi_results, f=f[1:5,], env_results=env_results, env=env)


# In case of multicovar<-TRUE in knobi_env, a list is required in which
# each item is a data frame for each environmental scenario

env<-list(climate_1=data.frame(AMO=c(0.2,0.2,0.3,0.3,0.4),
                              NAO=c(0.2,0.2,0.3,0.3,0.4)),
          climate_2=data.frame(AMO=c(0.2,0.3,0.4,0.5,0.6),
                              NAO=c(0.2,0.2,0.3,0.3,0.4)))

knobi_proj(knobi_results, knobi_environmental2, c=C, env=env)

KBPM retrospective analysis

Description

This function performs a retrospective analysis that evaluates the robustness of the KBPM fit to the systematic deletion of recent data.

Arguments

Details

There are different options for defining retrospective fits:

(1) Usage of nR argument. This argument specifies the number of retrospective peels. By using this argument, it is implied that the retrospective fits will consist of systematically deleting the last year of data, up to the number of years specified by nR.

(2) Usage of yR argument. This argument specifies the final years of the catch time series for each of the retrospective models, providing greater flexibility in choosing the years from which to delete information. The number of retrospective fits will correspond to the length of the yR vector. Additionally, different starting years can be set using the yR0 argument.

If both arguments, nR and yR, are provided, the package will prioritize the use of yR.

As described in the knobi_env function details, in the case of the environmental models, both the estimated biological reference points and the plotted production curve correspond to a value of the scaled environmental variable equal to the mean of the time series, i.e. X_{t}=0, which cancels out the environmental effect in the equations defining both models. For more details, such as the calculation of BRPs as a function of the environmental variable, see vignettes.

Value

A list containing the results of the retrospective analysis, including parameter estimates and reference points for each model. The estimated surplus production curves from the retrospective analysis are also plotted, with a panel where each graph represents the curves for each model in case of considering environmental models. The plot is displayed in the plot window and saved (if plot_out=TRUE) in the specified directory or in the current directory.

Author(s)

• Anxo Paz

• Marta Cousido-Rocha

• Santiago Cerviño López

• M. Grazia Pennino

Examples

library(knobi)

# See knobi_fit example to obtain the knobi_results object
knobi_retrospectives<-knobi_retro(knobi_results,plot_out=T)  # default nR=5
knobi_retrospectives

knobi_retro(knobi_results,nR=3)
knobi_retro(knobi_results,yR=c(2010,2015))
knobi_retro(knobi_results,yR=c(2010,2015),yR0=c(1995,2000))

# See knobi_env example to obtain the env_results object
knobi_retro(knobi_results,env_results=knobi_environmental,yR=c(2010,2015),yR0=c(1995,2000))

Print a knobi object

Description

The default print method for knobi_fit, knobi_env and knobi_proj object

Arguments

Details

Prints out the formula and the parameters estimates of the base KBPM fit or the environmental KBPM fit, furthermore it also reports the KBPM projections.

Author(s)

• Anxo Paz

• Marta Cousido-Rocha

• Santiago Cerviño López

• M. Grazia Pennino

See Also

knobi_fit, knobi_env, knobi_proj