Introduction

The inter-stock relationship functionality in MOM@Rel was developed with the intention to implement stochastic relationships fitted from R statistical models. For example, a linear regression of biomass of one species informs natural mortality in a prey species. The MSE simulation then generates M values for the prey species by sampling the linear regression model.

We may not want stochasticity in the inter-stock relationship, or we want to incorporate complex non-linear relationships that cannot be easily generated from a fitted R model (linear regression, GLM, etc.).

This article provides a template for users to develop a custom inter-stock relationship.

Create a function

First, we write a function whose arguments include the independent variables (see the multi-operating model article for the list of available options).

Here, we write a simple function where the total biomass determines the natural mortality of a stock that will be assigned later.

f <- function(B) {
  M <- ifelse(B > 1000, 0.2, 0.1)
  return(M)
}

Create a list

The function needs to be in a named list with auxiliary information. The required information includes the function and a character vector (labeled terms) that specify the relationship between the associated stocks. The first element is the dependent variable and subsequent elements are the independent variables.

A third entry in the list fitted.value is a placeholder for the predicted value during the simulation.

myRel <- list(
  f = f,
  terms = c("M_1", "B_2"),
  fitted.value = NA
)

Several additional items in the list are possible:

• If M is the dependent variable, the value can be age-specific. For example, specify M for age 3 with myRel$age <- 3 (remember that the age classes span zero to maxage).
• If recruitment deviates Perr_y is the dependent variable, then recruitment deviation in year y+1 is calculated from the independent variables in either the current year y with myRel$lag <- "current" (default), or within year y+1 with myRel$lag <- "next".

Class

Use the S3 system to give a class to the list.

class(myRel) <- "customRel"

Write generics

Both a predict and simulate generic need to be written for your custom Rel object:

The predict generic is a function that returns the value of the dependent variable from 2 arguments:

• object is the object of class customRel
• newdata will be a data frame containing the independent variable from the simulation

predict.customRel <- function(object, newdata) {
  if (missing(newdata)) {
    val <- NA
  } else {
    val <- object$f(newdata$`B_1`)
  }
  return(val)
}

newdata <- data.frame(B_1 = 100)
predict(myRel, newdata)

## [1] 0.1

multiMSE internally updates the list with the predicted value each time it is called (every year and simulation):

myRel$fitted.value <- predict(myRel, newdata)

The simulate generic is a function that will generate a value that will be used in the population dynamics. In this step, we can sample the value from a distribution or simply return the predicted value. Here we do the latter:

simulate.customRel <- function(object, stdev = 0) {
  val <- object$fitted.value
  val_sim <- val * exp(rlnorm(1, sd = stdev))
  return(val_sim)
}

These functions and objects should be in the global environment or in an R package (and exported) for multiMSE.