In this section we describe how to develop custom performance metric functions.

We’ll go through the `P50`

function in detail to see how it works:

`P50`

```
## function (MSEobj = NULL, Ref = 0.5, Yrs = NULL)
## {
## Yrs <- ChkYrs(Yrs, MSEobj)
## PMobj <- new("PMobj")
## PMobj@Name <- "Spawning Biomass relative to SBMSY"
## if (Ref != 1) {
## PMobj@Caption <- paste0("Prob. SB > ", Ref, " SBMSY (Years ",
## Yrs[1], " - ", Yrs[2], ")")
## }
## else {
## PMobj@Caption <- paste0("Prob. SB > SBMSY (Years ", Yrs[1],
## " - ", Yrs[2], ")")
## }
## PMobj@Ref <- Ref
## PMobj@Stat <- MSEobj@SB_SBMSY[, , Yrs[1]:Yrs[2]]
## PMobj@Prob <- calcProb(PMobj@Stat > PMobj@Ref, MSEobj)
## PMobj@Mean <- calcMean(PMobj@Prob)
## PMobj@MPs <- MSEobj@MPs
## PMobj
## }
## <bytecode: 0x000000001a531a20>
## <environment: namespace:MSEtool>
## attr(,"class")
## [1] "PM"
```

Functions of class `PM`

must have three arguments: `MSEobj`

, `Ref`

, and `Yrs`

:

`args(P50)`

```
## function (MSEobj = NULL, Ref = 0.5, Yrs = NULL)
## NULL
```

- The first argument
`MSEobj`

is obvious, an object of class`MSE`

to calculate the performance statistic. - The second argument
`Ref`

must have a default value. This is used as reference for the performance statistic, and will be demonstrated shortly. - The third argument
`Yrs`

can have a default value of`NULL`

or specify a numeric vector of length 2 with the first and last years to calculate the performance statistic, or a numeric vector of length 1 in which case if it is positive it is the first`Yrs`

and if negative the last`Yrs`

of the projection period.

The first line of a `PM`

function must be `Yrs <- ChkYrs(Yrs, MSEobj)`

. This line updates the `Yrs`

variable and makes sure that the specified year indices are valid. For example:

```
MSEobj <- runMSE(silent=TRUE) # example MSE object
ChkYrs(NULL, MSEobj) # returns all projection years
```

`## [1] 1 50`

`ChkYrs(c(1,10), MSEobj) # returns first 10 years `

`## [1] 1 10`

`ChkYrs(c(60,80), MSEobj) # returns message and last 20 years`

```
## Yrs[2] is greater than MSEobj@proyears. Setting Yrs[2] =
## MSEobj@proyears
```

```
## Yrs[1] is greater than MSEobj@proyears. Setting Yrs[1] = Yrs[2] -
## Yrs[1]
```

`## [1] 30 50`

`ChkYrs(5, MSEobj) # first 5 years`

`## [1] 1 5`

`ChkYrs(-5, MSEobj) # last 5 years`

`## [1] 46 50`

`ChkYrs(c(50,10), MSEobj) # returns an error`

`## Error: Yrs[1] is > Yrs[2]`

When the default value for `Yrs`

is `NULL`

, the `Yrs`

variable is updated to include all projection years:

```
Yrs <- ChkYrs(NULL, MSEobj)
Yrs
```

`## [1] 1 50`

#### Creating a PM function

We’ll go through the steps of creating the `P50`

function.

First, we create a new object of class `PMobj`

, and populate the `Name`

slot with a short but descriptive name:

```
PMobj <- new("PMobj")
PMobj@Name <- "Spawning Biomass relative to SBMSY"
```

The next line populates the `Caption`

slot with a brief caption including the years over which the performance statistic is calculated. The *if* statement is not crucial, but avoids the redundant `SB > 1 SBMSY`

in cases where `Ref=1`

.

Next we store the value of the `Ref`

argument in the `PMobj@Ref`

slot so that information is contained in the function output.

`PMobj@Ref <- Ref`

The `Stat`

slot is an array that stores the variable which we wish to calculate the performance statistic; an output from the `runMSE`

function with dimensions `MSE@nsim`

, `MSE@nMPs`

, and `MSE@proyears`

(or fewer if the argument `Yrs != NULL`

).

In this case we want to calculate a performance statistic related to the biomass relative to BMSY, and so we assign the `Stat`

slot as follows:

`PMobj@Stat <- MSEobj@SB_SBMSY[, , Yrs[1]:Yrs[2]]`

Note that we are including all simulations and MPs and indexing the years specified in `Yrs`

.

Next we use the `calcProb`

function to calculate the mean of `PMobj@Stat > PMobj@Ref`

over the years dimension. This results in a matrix with dimensions `MSE@nsim, MSE@nMPs`

:

`PMobj@Prob <- calcProb(PMobj@Stat > PMobj@Ref, MSEobj)`

Note that in order to calculate a probability the argument to the `calcProb`

function must be a logical array, which is achieved using the `Ref`

slot.

Also note that in this case `PMobj@Stat > PMobj@Ref`

is equivalent to `MSEobj@SB_SBMSY[, , Yrs[1]:Yrs[2]] > 0.5`

.

The `PM`

functions have been designed this way so that in most cases the `PMobj@Prob <- calcProb(PMobj@Stat > PMobj@Ref)`

line is identical in all `PM`

functions and does not need to be modified. The exception to this is if we don’t want to calculate a probability but want the actual mean values of `PMobj@Stat`

, demonstrated in the example below.

In the next line we calculate the mean of `PMobj@Prob`

over simulations using the `calcMean`

function:

`PMobj@Mean <- calcMean(PMobj@Prob)`

Similar to the previous line, this line is identical in all `PM`

functions and can be simply copy/pasted from other `PM`

functions without being modified. The `Mean`

slot is a numeric vector of length `MSEobj@nMPs`

with the overall performance statistic, in this case the probability of \(\text{SB} > 0.5\text{SB}_\text{MSY}\) across all simulations and years.

Finally, we store the names of the MPs and return the `PMobj`

.

#### Creating Example `PMs`

and Plot

As an example, we will create another version of `DFO_plot`

using some custom `PM`

functions and a customized version of `TradePlot`

.

First we create the plot using `DFO_plot`

:

`DFO_plot(MSEobj)`

From the help documentation (`?DFO_plot`

) we can see that this function plots mean spawning biomass relative to \(\text{SB}_\text{MSY}\) and fishing mortality rate relative to \(F_\text{MSY}\)
over the final 5 years of the projection.

First we’ll develop a `PM`

function to calculate the mean \(\text{SB}_\text{MSY}\) for the last 5 years of the projection period. Notice that this is very similar to `P50`

described above, with the modification of the `Caption`

and the `Prob`

slots, and the `Yrs`

argument. We are calculating a mean here instead of a probability and are not using the `Ref`

argument:

```
MeanB <- function(MSEobj = NULL, Ref = 1, Yrs = -5) {
Yrs <- ChkYrs(Yrs, MSEobj)
PMobj <- new("PMobj")
PMobj@Name <- "Spawning Biomass relative to SBMSY"
PMobj@Caption <- paste0("Mean SB/SBMSY (Years ", Yrs[1], " - ", Yrs[2], ")")
PMobj@Ref <- Ref
PMobj@Stat <- MSEobj@SB_SBMSY[, , Yrs[1]:Yrs[2]]
PMobj@Prob <- calcProb(PMobj@Stat, MSEobj)
PMobj@Mean <- calcMean(PMobj@Prob)
PMobj@MPs <- MSEobj@MPs
PMobj
}
```

We develop a `PM`

function to calculate average F/FMSY in a similar way:

```
MeanF <- function(MSEobj = NULL, Ref = 1, Yrs = -5) {
Yrs <- ChkYrs(Yrs, MSEobj)
PMobj <- new("PMobj")
PMobj@Name <- "Fishing Mortality relative to FMSY"
PMobj@Caption <- paste0("Mean F/FMSY (Years ", Yrs[1], " - ", Yrs[2], ")")
PMobj@Ref <- Ref
PMobj@Stat <- MSEobj@F_FMSY[, , Yrs[1]:Yrs[2]]
PMobj@Prob <- calcProb(PMobj@Stat, MSEobj)
PMobj@Mean <- calcMean(PMobj@Prob)
PMobj@MPs <- MSEobj@MPs
PMobj
}
```

Similar to developing custom MPs we need to tell R that these new functions are `PM`

methods:

```
class(MeanB) <- "PM"
class(MeanF) <- "PM"
```

Now we can test our performance metric functions:

```
data.frame(MP=MeanB(MSEobj)@MPs,
SB_SBMSY=MeanB(MSEobj)@Mean,
F_FMSY=MeanF(MSEobj)@Mean)
```

```
## MP SB_SBMSY F_FMSY
## 1 curEref 0.3258801 2.186392e+00
## 2 FMSYref 0.8950520 1.008874e+00
## 3 FMSYref50 1.9335515 5.044369e-01
## 4 FMSYref75 1.3170624 7.566554e-01
## 5 NFref 4.2925589 4.017549e-15
```

How do these results compare to what is shown in `DFO_plot`

?

We could also use the `summary`

function with our new `PM`

functions, but note that these results are not probabilities:

`summary(MSEobj, 'MeanB', 'MeanF')`

`## Calculating Performance Metrics`

```
## Performance.Metrics
## 1 Spawning Biomass relative to SBMSY Mean SB/SBMSY (Years 46 - 50)
## 2 Fishing Mortality relative to FMSY Mean F/FMSY (Years 46 - 50)
##
##
## Performance Statistics:
## MP MeanB MeanF
## 1 curEref 0.33 2.2e+00
## 2 FMSYref 0.90 1.0e+00
## 3 FMSYref50 1.90 5.0e-01
## 4 FMSYref75 1.30 7.6e-01
## 5 NFref 4.30 4.0e-15
```

Finally, we will develop a customized plotting function to reproduce the image produced by `DFO_plot`

.

We can produced something fairly similar quite quickly using the `TradePlot`

function:

`TradePlot(MSEobj, 'MeanB', 'MeanF', Lims=c(0,0))`

```
## MP MeanB MeanF Satisificed
## 1 curEref 0.33 2.2e+00 TRUE
## 2 FMSYref 0.90 1.0e+00 TRUE
## 3 FMSYref50 1.90 5.0e-01 TRUE
## 4 FMSYref75 1.30 7.6e-01 TRUE
## 5 NFref 4.30 4.0e-15 TRUE
```

Adding the shaded polygons and text requires a little more tweaking and some knowledge of the `ggplot2`

package. We will wrap up our code in a function:

```
NewPlot <- function(MSEobj) {
# create but don't show the plot
P <- TradePlot(MSEobj, 'MeanB', 'MeanF', Lims=c(0,0), Show='none')
P1 <- P$Plots[[1]] # the ggplot objects are returned as a list
# add the shaded regions and the text
P1 <- P1 + ggplot2::geom_rect(ggplot2::aes(xmin=c(0,0,0), xmax=c(0.4, 0.8, Inf),
ymin=c(0,0,1), ymax=rep(Inf,3)),
alpha=c(0.8, 0.6, 0.4), fill="grey86") +
ggplot2::annotate(geom = "text", x = c(0.25, 0.6, 1.25), y = Inf,
label = c("Critical", "Cautious", "Healthy") ,
color = c('white', 'darkgray', 'darkgray'),
size=5, angle = 270, hjust=-0.25)
# Re-order layers so text and points are not covered
P1$layers <- c(P1$layers, P1$layers[[3]], P1$layers[[4]])
P1$layers[[3]] <- P1$layers[[4]] <- NULL
P1
}
NewPlot(MSEobj)
```