knobi_fit: 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.

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.

Arguments

data

A list containing the data.

Catch: catch time series observations.
Biomass: biomass time series estimated through a data-rich stock assessment model. If available, otherwise enter SSB in the next argument.
SSB: spawning stock biomass time series estimated through a data-rich stock assessment model. If available, otherwise introduce biomass in the previous argument.
years: optional. Years associated to the observed catch time series. By default, an increasing natural sequence from 1 to the total number of catch observations.
Stock: optional. Character string with the stock name for using in the plot subtitles.
Recruitment: optional. Recruitment time series estimated through a data-rich stock assessment model. See details.
F_input: optional. Fishing mortality time series estimated through a data-rich stock assessment model. See details.
RP: optional. Estimates of the following biological reference points derived from the fit of a data-rich stock assessment model (see details):
- MSY: Maximum Sustainable Yield.
- F_MSY: fishing mortality at MSY.
- B_MSY: biomass at MSY (or SSB at MSY depending on the value of the 'method' argument, see control list).
- K: carrying capacity.

control

A list containing the following control parameters.

pella: logical. If TRUE, the Pella-Tomlinson model is applied; if FALSE (default), the Schaefer model is used. See details.
start_r: optional. Starting value for the growth rate parameter r (intrinsic rate of natural increase). See details.
start_K: optional. Starting value for the carrying capacity parameter K (maximum population size where growth remains positive). See details.
start_p: optional. Starting value for the shape parameter p, used exclusively in the Pella-Tomlinson model. See details.
method: specifies whether the model is fitted using 'SSB' (spawning stock biomass) or 'Biomass'. This argument is only required if both time series, 'SSB' and 'Biomass', are provided.

plot_out

Logical. If TRUE, files are generated with plots for both the input time series and the fitting results. Defaults to FALSE, which means that the plots will be displayed only in the window.

plot_dir

optional. Directory for saving the plot files. Required when plot_out=TRUE. Current directory by default.

plot_filename

optional. Name of the folder that will contain the plot files. By default, "knobi_results". Required when plot_out=TRUE.

Author

Anxo Paz
Marta Cousido-Rocha
Santiago Cerviño López
Maria Grazia Pennino

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).

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

Run this code


# \donttest{

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

# }

Run the code above in your browser using DataLab