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