## How to write a custom Rel

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