Change in stock-recruit relationship

While the stock-recruit relationship is constant in openMSE, it is possible to change the effective stock-recruit relationship by changing the mean \(\mu\) of the recruitment deviations through OM@cpars$Perr_y. This article describes how reference points implicitly change with \(\mu\). The internal code used to calculate annual reference points in MSEtool::runMSE is not altered when the mean of OM@cpars$Perr_y is altered, so it is up to the user to decide whether reference points (see code below) need to be re-calculated.

Stock-recruit parameters

With the Beverton-Holt stock recruit relationship, the recruitment in year \(y\) is:

\[ R_y = \dfrac{\alpha SB_y}{1 + \beta SB_y} \delta_y\]

where \(\delta_y\) are recruitment deviates. For stationary productivity around the recruitment predicted by the S-R curve, we sample \(\delta_y\) from a lognormal distribution with \(E[\delta_y] = 1\) and \(V[\delta_y] = \exp(\sigma^2-1)\). These values are passed to the operating model in OM@cpars$Perr_y. Parameters \(\alpha\) and \(\beta\), along with \(\phi_0\), can be used to calculate the corresponding unfished recruitment \(R_0\) and steepness \(h\):

\[ h = \dfrac{\alpha \phi_0}{4 + \alpha \phi_0}\]

\[ R_0 = \dfrac{1}{\beta}\left(\alpha - \dfrac{1}{\phi_0}\right)\]

We define the unfished biomass as \(SB_0 = R_0 \phi_0\). For easier reading, we have dropped the superscript \(\textrm{SR}\) associated with \(R_0\), \(h\), and \(\phi_0\).

If we want to change the productivity via recruitment deviations, we can re-scale the recruitment deviations. Set \(\delta'_y = \mu \delta_y\), with \(\mu < 1\) for a less productive system. In effect, we will have created a different stock-recruit relationship:

\[ R_y = \dfrac{\alpha SB_y}{1 + \beta SB_y} \delta'_y = \dfrac{\alpha' SB_y}{1 + \beta SB_y} \delta_y\]

where \(\alpha' = \mu \alpha\), \(E[\delta'_y] = \mu\) and \(V[\delta'_y] = \exp(\sigma^2-1)/\mu^2\), and \(\beta\) is unchanged. In other words, by setting the mean of OM@cpars$Perr_y to \(\mu\), we create a stock-recruit relationship defined by \(\alpha'\) and recruitment variance of \(V[\delta'_y]\).

Implied reference points

Annual unfished and MSY reference points can be re-calculated accordingly with \(\alpha'\), \(\beta\), and \(\phi_0\). The implications on the reference points can be seen in “new” parameters \(R_0'\) and \(h'\) corresponding to “new” \(\alpha'\):

\[ h' = \dfrac{\alpha' \phi_0}{4 + \alpha' \phi_0} =\dfrac{\mu h}{1 + h (\mu - 1)}\]

\[ R_0' = \dfrac{1}{\beta}\left(\alpha' - \dfrac{1}{\phi_0}\right) = R_0\dfrac{h(4\mu + 1) - 1}{5h-1}\] Then, \(SB_0' = R_0' \phi_0\).

In the simplest case with constant \(\phi_0\), we expect \(h' < h\) and \(R_0' < R_0\) when \(\mu < 1\), and also expect \(SB_0\), \(F_{\textrm{MSY}}\), \(\textrm{MSY}\), and \(SB_\textrm{MSY}\) to decrease.

The new \(SPR_{\textrm{crash}} = (\alpha'\phi_0)^{-1} = (\mu \alpha \phi_0)^{-1}\) increases (\(F_{\textrm{crash}}\) decreases) when \(\mu < 1\).

Examples

Contour plots

Let’s explore the ratio of the new and old values of \(R_0\) and \(h\) as a function of \(\mu\) and the original steepness, with \(\beta = 1\) and \(\phi_0 = 1\):

Interestingly, there can exist a combination of low steepness and low \(\mu\) such that the population crashes, as implied by \(R_0' \le 0\) in red.

Forward projections

Let’s explore the dynamics in a forward projecting model in MSEtool. We’ll create an operating model with 4 simulations with varying \(\mu\) but constant stock-recruit parameters among simulations. Then we’ll project the model using the NFref (no fishing) MP.

library(MSEtool)
## Loading required package: snowfall
## Loading required package: snow
OM <- MSEtool::testOM

# Generate four simulations with same alpha, beta, and phi_0
OM@nsim <- 4
OM@Linfsd <- OM@Msd <- OM@Ksd <- c(0, 0)
OM@Linf <- mean(OM@Linf) |> rep(2)
OM@K <- mean(OM@K) |> rep(2)
OM@t0 <- mean(OM@t0) |> rep(2)
OM@M <- mean(OM@M) |> rep(2)
OM@L50 <- mean(OM@L50) |> rep(2)
OM@L50_95 <- mean(OM@L50_95) |> rep(2)
OM@h <- mean(OM@h) |> rep(2)

OM@Prob_staying <- OM@Frac_area_1 <- OM@Size_area_1 <- c(0.5, 0.5)

# With four simulations, let's set the projection Perr_y to these values:
mu_vec <- c(0.01, 0.75, 1, 1.5)
OM@cpars$Perr_y <- matrix(1, OM@nsim, OM@maxage + OM@nyears + OM@proyears)
OM@cpars$Perr_y[, OM@maxage + OM@nyears + 1:OM@proyears] <- mu_vec

Hist <- Simulate(OM, silent = TRUE)
MSE <- Project(Hist, MPs = "NFref", silent = TRUE)

Let’s see what happens with the NFref MP:

Despite the same stock-recruit parameters, the spawning biomass all diverge with \(\mu\):

This terminal SB should be the new \(SB_0\), so let’s create a function to calculate these values analytically instead of projecting:

recalc_ref_pt <- function(Hist, mu = rep(1, Hist@OM@nsim)) {
  StockPars <- Hist@SampPars$Stock
  FleetPars <- Hist@SampPars$Fleet
  R0 <- StockPars$R0
  hs <- StockPars$hs
  
  if(StockPars$SRrel[1] == 1) {
    R0_new <- R0 * (hs * (4 * mu + 1) - 1)/(5 * hs - 1)
    hs_new <- mu * hs /(1 + hs * (mu - 1))
  } else {
    R0_new <- 1 + 0.8 * R0 * log(mu)/log(5 * hs)
    hs_new <- hs * mu^0.8
  }
  
  MSY_ref_pt <- lapply(1:Hist@OM@nsim, function(x) { # Internal function for calculating unfished and MSY reference points
    sapply(Hist@OM@nyears + 1:Hist@OM@proyears, function(y) {
      MSEtool:::optMSY_eq(x, 
                          M_ageArray=StockPars$M_ageArray, 
                          Wt_age=StockPars$Wt_age, 
                          Mat_age=StockPars$Mat_age,
                          Fec_age=StockPars$Fec_Age, 
                          V=FleetPars$V_real, 
                          maxage=StockPars$maxage, 
                          R0=R0_new,
                          SRrel=StockPars$SRrel, 
                          hs=hs_new, 
                          SSBpR=StockPars$SSBpR,
                          yr.ind=y, 
                          plusgroup=StockPars$plusgroup)
    })
  })
  
  crash_ref_pt <- local({ # Internal function for calculating SPR crash
    boundsF <- c(1E-3, 3)
    F_search <- exp(seq(log(min(boundsF)), log(max(boundsF)), length.out = 50))
    
    lapply(1:Hist@OM@nsim, function(x) {
      sapply(Hist@OM@nyears + 1:Hist@OM@proyears, function(y) {
        Ref_search <- MSEtool:::Ref_int_cpp(F_search, 
                                            M_at_Age = StockPars$M_ageArray[x, , y],
                                            Wt_at_Age = StockPars$Wt_age[x, , y], 
                                            Mat_at_Age = StockPars$Mat_age[x, , y],
                                            Fec_at_Age = StockPars$Fec_Age[x, , y], 
                                            V_at_Age = FleetPars$V_real[x, , y],
                                            maxage = StockPars$maxage,
                                            plusgroup = StockPars$plusgroup)
        SPR_search <- Ref_search[2, ]
        RPS <- Ref_search[3, ]
        
        if(StockPars$SRrel[x] == 1) {
          CR <- 4 * StockPars$hs[x]/(1 - StockPars$hs[x])
        } else if(StockPars$SRrel[x] == 2) {
          CR <- (5 * StockPars$hs[x])^1.25
        }
        alpha <- mu[x] * CR/StockPars$SSBpR[x, 1] # New alpha
        
        if(min(RPS) >= alpha) { # Unfished RPS is steeper than alpha
          SPRcrash <- min(1, RPS[1]/alpha) # Should be 1
          Fcrash <- 0
        } else {
          SPRcrash <- MSEtool:::LinInterp_cpp(RPS, SPR_search, xlev = alpha)
          Fcrash <- MSEtool:::LinInterp_cpp(RPS, F_search, xlev = alpha)
        }
        c(SPRcrash = SPRcrash, Fcrash = Fcrash)
      })
    })
  })
  ref_pt <- abind::abind(simplify2array(MSY_ref_pt), simplify2array(crash_ref_pt), along = 1)
  dimnames(ref_pt)[[1]] <- c("MSY", "FMSY", "SBMSY", "SBMSY_SB0", "BMSY_B0", "BMSY", "VBMSY", "VBMSY_VB0", "RMSY", "SB0", "B0", "R0", "h",
                             "N0", "SN0", "SPRcrash", "Fcrash")
  return(ref_pt)
}

Let’s compare the ratio of the calculated \(SB_0\) and the projected terminal SSB:

new_ref_pt <- recalc_ref_pt(Hist, mu = mu_vec)
SSB0_prime <- new_ref_pt["SB0", OM@proyears, ]
plot(mu_vec, SSB0_prime/terminal_SB, xlab = expression(mu), ylab = "SB0/Terminal SB")

Looks like our new reference points match the values from the projections. Note with \(\mu = 0.01\), the population crashes so \(SB_0 = 0\).

Now let’s map out \(h'/h\) and \(R_0'/R_0\) as a function of \(\mu\):

If \(\mu < 0.08\), our population crashes because \(R_0' < 0\), and \(SPR_{\textrm{crash}} = 1\):

SPRcrash <- new_ref_pt["SPRcrash", OM@proyears, ]
plot(mu_vec, SPRcrash, xlab = expression(mu), ylim = c(0, 1))

Ricker function

With a Ricker function where \(R_y = \alpha SB_y \exp(-\beta SB_y) \delta_y\), steepness is \(h = 0.2(\alpha\phi_0)^{0.8}\) and \(R_0 = \log(\alpha\phi_0)/(\beta \phi_0)\).

The new SR relationship is \(R_y = \alpha' SB_y \exp(-\beta SB_y) \delta_y\), with

\[h' = 0.2 (\alpha' \phi_0)^{0.8} = \mu^{0.8} h\]

\[R_0' = 1 + \dfrac{0.8 R_0\log(\mu)}{\log(5h)}\]