Fit a "fisheries" GP model

Fits an alternate version of the GP model with one additonal parameter (catchability), designed for use in fisheries applications.

fitGP_fish(
  data,
  y,
  m,
  h,
  z = NULL,
  pop = NULL,
  time = NULL,
  scaling = c("global", "local", "none"),
  initpars = NULL,
  modeprior = 1,
  fixedpars = NULL,
  rhofixed = NULL,
  rhomatrix = NULL,
  augdata = NULL,
  linprior = "none",
  predictmethod = NULL,
  newdata = NULL,
  xname = "escapement",
  ytrans = c("none", "log", "gr1", "gr2"),
  bfixed = NULL,
  bshared = TRUE
)

Arguments

data

A data frame, required.

y

The response variable (required). Typically CPUE.

m

Lags of the response variable (required). Typically lags of CPUE.

h

Lags of the variable to be multipled by b and subtracted from m (required). Typically harvest. The dimensions of m and h must match.

z

Other predictor variables, e.g. covariates, that are unmodified (optional).

pop

Identifies separate populations (optional, if not supplied, defaults to 1 population). Population values can be either numeric, character, or factor.

time

The time index (recommended). If not supplied, defaults to a numeric index.

scaling

How the variables should be scaled (see Details). Scaling can be "global" (default), "local" (only applicable if more than 1 pop), or "none". Scaling is applied to y and each x. For a mix of scalings, do it manually beforehand and set scaling to "none". All outputs are backtransformed to original scale.

initpars

Starting values for hyperparameters (see Details) in the order c(phi,ve,sigma2,rho), in constrained (non-transformed) space. Should be a numeric vector with length ncol(x)+3. Defaults used if not supplied: c(rep(0.1,ncol(x)), 0.001, 1-0.001, 0.5)

modeprior

This value is used by the phi prior and sets the expected number of modes over the unit interval. Defaults to 1.

fixedpars

Fixes values of the hyperparameters phi, ve, and sigma2 (if desired). Should be a numeric vector with length ncol(x)+2 in the order c(phi,ve,sigma2) in constrained (non-transformed) space. Hyperparameters to estimate should be input as NA. The value of rho is fixed separately using argument rhofixed.

rhofixed

Fixes the rho parameter, if desired (see Details).

rhomatrix

Symmetrical square matrix of fixed pairwise rho values to use, with 1's on the diagonal. The rows and columns should be named with the population identifier. The output of getrhomatrix can be used here. All populations in the dataset must appear in rhomatrix. Populations in rhomatrix that are not in the dataset are allowed and will not be used.

augdata

A data frame with augmentation data (see Details).

linprior

Fit GP model to the residuals of a linear relationship between y and the first variable of x (prior to scaling) and backtransforms as appropriate. If y is growth rate, acts like a Ricker prior. Defaults to "none". Option "local" fits a separate regression for each population, "global" fits a single regression.

predictmethod

Using the training data, loo does leave-one-out prediction, lto does leave-timepoint-out prediction, and sequential does sequential (leave-future-out) prediction.

newdata

Data frame containing the same columns supplied in the original model.

xname

What the composite variable m-bh should be called. Defaults to "escapement".

ytrans

Tranformation to apply to y before fitting (m remains untransformed). Either "none" (default), "log", "gr1", or "gr2". R2 is calculated on y in its original units.

bfixed

Fixes b and bypasses optimization. If there are multiple pops, can be a scalar (shared across pops), or a named vector with different b's for each pop.

bshared

If there are multiple pops, should they share the same b (TRUE, default) or have different values of b (FALSE)? Note that optimizing multiple b values can be slow. Ignored if bfixed is supplied. If bshared=FALSE and either pop is not supplied or there is only one pop, reverts back to a single b estimate.

Value

A list (class GP and GPpred) with the same elements as fitGP and with additonal element b, the names of m, h, z stored under inputs, and the composite variable (escapement) included in the insampresults table.

Details

This fits a GP model of the form $$y=f(m-bh,z)$$ where \(y\) is catch per unit effort (CPUE), \(m\) are lags of CPUE, \(h\) are lags of harvest (in numbers or biomass), \(b\) is catchability (scalar), and \(z\) are optional covariates. CPUE is assumed proportional to total biomass (or numbers), with proportionality constant \(b\). The composite variable \(m-bh\) is the biomass or number of individuals remaining after harvesting (escapement).

Parameter \(b\) is found using optimize applied to the posterior likelihood. Alternatively, a fixed value for \(b\) can be provided under bfixed.

If fitting a hierarchical model, the default behavior is to estimate a single value of\(b\) shared by all pops (bshared=TRUE). You can fix a single shared value of \(b\) by providing a single value under bfixed. Alternatively, you can estimate different values of \(b\) for each population by setting bshared=FALSE. This will use the Nelder-Mead method of optim (and will be quite a bit slower than the single parameter optimization). You can fix different values of \(b\) for each population by supplying a named vector under bfixed.

Parameter ytrans applies a transformation to \(y\) before fitting the model. Inputs y and m should be in untransformed CPUE. ytrans="none" (the default) will apply no tranformation, ytrans="log" with compute \(log(y)\), ytrans="gr1" will compute \(log(y_t/m_{t-1})\), and ytrans="gr2" will compute \(log(y_t/(m_{t-1}-bh_{t-1})\).

Using this method requires the use of data with pre-generated lags (option A1 in fitGP). For more details on fitting a fisheries model and an example see the vignette. For more elaboration on the inputs, (e.g. scaling, augdata) see fitGP.