An explicit solution for calculating optimum spawning stock size from Ricker’s stock recruitment model

Introduction

One of the most difficult problems in the assessment of fish stocks is establishing the relationship between the spawning stock and subsequent recruitment (Hilborn & Walters, 1992). Stock-recruitment models have been used for decades in fisheries management as a means of formalizing this relationship (Beverton & Holt, 1957; Ricker, 1954). Over time, a variety of functional forms have emerged to capture varying assumptions about depensatory and compensatory mortality (Hilborn & Walters, 1992). In a classroom setting, deterministic versions of the models provide useful constructs for teaching about management reference points such as maximum sustained yield (MSY).

In particular, Ricker’s stock recruitment model (Ricker, 1954; Ricker, 1975) is one of the most widely used models to describe the population dynamics of fishes, such that (1)${\displaystyle R=\alpha S{e}^{-bS},}$ R is the number of recruits produced, S is the number of spawners, α is the dimensionless number of recruits per spawner produced at very low spawner density, and b is the strength of density dependence (units: spawner⁻¹). It is common to substitute α = e^a into Eq. (1) and rewrite it as (2)${\displaystyle R=S{e}^{a-bS}.}$ To make the model reflect a stochastic process, Eq. (2) is typically multiplied by a log-normal error term, so that (3)${\displaystyle R=S{e}^{a-bS}{e}^{\epsilon },}$ and ε is a normally distributed error term with a mean of −1∕2 σ and variance σ. This non-zero mean ensures that a is interpreted as the mean recruits per spawner rather than the median (Hilborn, 1985). Part of the model’s popularity is due to the relative ease with which its parameters are estimated. After log transformation, Eq. (3) is typically rewritten as (4)${\displaystyle ln\left(R∕S\right)=a-bS+\epsilon ,}$ and the parameters are estimated via simple linear regression. I note here that estimation of the parameters via a simple observation-error model like (4) can lead to substantial biases in a and b if the sample size is low (n ≤ 10) due to autocorrelation in the residuals ε (Walters, 1985).

Once the model has been fit to data and any necessary bias corrections made, the parameters can be used to derive various biological reference points of interest to fisheries managers. Some of these metrics are rather trivial to compute. For example, the spawning stock size leading to maximum recruit production (S_MSR) is simply 1/b. However, other reference points are much less straightforward to calculate. In particular, the spawning stock expected to produce the maximum sustainable yield (S_MSY) under deterministic dynamics is of common interest.

To find S_MSY, I express the yield (Y) as (5)${\displaystyle Y=R-S=S{e}^{a-bS}-S,}$ and then take the derivative of Y with respect to S: (6)${\displaystyle \frac{dY}{dS}=\left(1-bS\right)e^{a-bS}-1.}$ S_MSY is then determined by setting Eq. (6) to zero and solving for S. Upon initial inspection, however, there does not appear to be an explicit solution to this equation in terms of S, and therefore S_MSY is typically solved “by trial” (Ricker, 1975) with some form of gradient method (e.g., Newton’s as in Hilborn, 1985).

To simplify this issue for common applications, Hilborn (1985) developed a simple model whereby the ratio of spawning stock size at MSY to that at the unfished equilibrium (S_MSY∕S_r) is a linear function of the parameter a. Specifically, for 0 < a ≤ 3 he estimated that (7a)${\displaystyle \frac{S_{\mathrm{MSY}}}{S_r}=\frac{S_{\mathrm{MSY}}}{\left(a∕b\right)}=0.5-0.07a,\text{and}}$ (7b)${\displaystyle S_{\mathrm{MSY}}=\frac{a\left(0.5-0.07a\right)}{b}.}$ Although this approximation is very useful due to its simplicity, there is no underlying fundamental support for the statistical form of the relationship.

Methods

Here I make use of the Lambert W function, W(z), to demonstrate an explicit solution to Eq. (4) that precludes the need to estimate S_MSY via numerical methods or Hilborn’s (1985) linear approximation. This function has been used for explicit solutions to Roger’s random predator equation in ecology (McCoy & Bolker, 2008) and susceptible-infected-removed (SIR) models in epidemiology (Reluga, 2004; Wang, 2010). Specifically, W(z) is defined as the function that satisfies (8)${\displaystyle W\left(z\right)e^{W\left(z\right)}=z}$ for any complex number z (Lambert 1758 and Euler 1783 as cited in Corless et al., 1996). Here we are interested only in real values, however, so I replace z with x and note that W(x) is only defined for x ≥ − 1∕e (Corless et al., 1996). Furthermore, this function is not injective and has two values for −1∕e ≤ x ≤ 0, but as I show below, we are concerned only with the region where x > 0 and W(x) is a singular, non-negative value.

I begin my explicit solution of S_MSY by setting Eq. (6) to zero, such that (9)${\displaystyle \left(1-b{S}_{\mathrm{MSY}}\right)e^{a-b{S}_{\mathrm{MSY}}}=1.}$ After rearranging terms and multiplying both sides by e, we arrive at (10)${\displaystyle \left(1-b{S}_{\mathrm{MSY}}\right)e^{1-b{S}_{\mathrm{MSY}}}=e^{1-a}.}$ At this point I note the relationship between Eqs. (10) and (8), with 1 − bS_MSY = W(z) and e^1−a = z. Therefore, we can write (11)${\displaystyle 1-b{S}_{\mathrm{MSY}}=W\left(e^{1-a}\right),\text{and hence}}$ (12)${\displaystyle S_{\mathrm{MSY}}=\frac{1-W\left(e^{1-a}\right)}{b}.}$ We now have an explicit solution for S_MSY that depends only on the parameters a and b from Eq. (2). As mentioned above, W(x) is only defined for x ≥ − 1∕e, which does not pose any problems here because x = e^1−a > 0∀a ∈ ℝ. For visualization purposes, I show a plot of W(e^1−a) versus a in Fig. 1.

We can also derive an explicit formula for calculating the fraction of the return harvested at S_MSY, which I call U_MSY. As Ricker (1975) shows, (13)${\displaystyle U_{\mathrm{MSY}}=b{S}_{\mathrm{MSY}},}$ and therefore substituting (12) into (13) gives (14)${\displaystyle U_{\mathrm{MSY}}=1-W\left(e^{1-a}\right).}$

In practice W(x) may be approximated numerically using some form of gradient method. Corless et al. (1996) recommend Halley’s method, with the update equation given by (15)${\displaystyle w_{j+1}=w_j-\frac{w_j{e}^{w_j}-x}{e^{w_j}\left(w_j+1\right)-\frac{\left(w_j+2\right)\left(w_j{e}^{w_j}-x\right)}{2w_j+2}}.}$ I used an initial guess of w₀ = 3∕4ln(x + 1) based on the shape of W(x) over the range of a typically considered in fisheries research (i.e., 0 < a < 3 as in Hilborn, 1985). If, however, one must estimate W(x) numerically, then one should ask whether doing so is, in fact, computationally faster. Therefore, as a test I randomly selected 1,000 values each for 0 < a ≤ 3 and 10⁻⁵ ≤ b ≤ 10⁻³, and then solved for S_MSY using both Newton’s method as suggested by Ricker (1975), and Halley’s method as in Eq. (15).

Results and Discussion

Recent analyses have relied on estimating S_MSY via Hilborn’s (1985) linear approximation when calculating optimal yield profiles (Fleischman et al., 2013) or the effects of observation error on biases in parameter estimates (Su & Peterman, 2012). On the other hand, solving for S_MSY using W(x) and Halley’s method is not only convenient; it also offers an appreciable computational advantage over the standard Newton method. Although both methods converged in less than 10 iterations during my test, Halley’s method was always faster and less variable overall (Fig. 2). Therefore, estimating S_MSY via Halley’s method might save significant time in applications such as management strategy evaluations that are much more computationally intensive than a simple one-case solution.

Although implementing Eq. (15) may seem a bit daunting to individuals less familiar with numerical methods, a variety of contemporary software packages (e.g., MATLAB, R) include built-in functions to calculate W(x) directly. This means that anyone using a personal computer to estimate the parameters in a Ricker model can easily estimate S_MSY from Eq. (12) as I demonstrate in Table 1; I show the results from my R implementation for a range of a and b in Fig. 3. For those preferring to use Microsoft Excel, there is no built-in function to calculate W(x), but I have implemented Eq. (15) as the VBA function ‘LAMBERTW’ in the Microsoft Excel add-in file ‘LambertWfunc.xlam’ (see Fig. S1 for download and install instructions).

> library("gsl")
> a = 1
> b = 5e-4
> Smsy = (1 − lambert_W0(exp(1 - a))) / b

>> a = 1
>> b = 5e-4
>> Smsy = (1 - lambertw(exp(1 - a))) / b

Here I have outlined a new method to easily calculate S_MSY from the productivity (a) and density-dependent (b) parameters in a Ricker model using readily available functions in several software packages. This method is much more straightforward than trying to solve for S_MSY using numerical methods and should be useful in many classroom settings. Although there could be some utility in actually going through the exercise of numerically deriving the answer, it is rare nowadays, for example, for anyone to code a random number generator because of their ubiquitous implementation in standard software. In addition, the explicit analytical solution is closed-form with respect to the special functions, and therefore precludes the need to estimate S_MSY via Hilborn’s (1985) approximation. Thus, due to the speed and ease with which these new equations are calculated, I recommend that practitioners use them for the estimation of S_MSY and U_MSY in lieu of those listed in Appendix III of Ricker (1975) and Table 7.2 of Hilborn & Walters (1992).

Supplemental Information