Stock-production functions resembling Beverton-Holt

Researchers dissatisfied with the performance of the Beverton-Holt model, in contexts where “Beverton-Holt-like” behavior is expected, have introduced a plethora of alternative model forms. This paper presents a formalization of what constitutes “Beverton-Holt-like behavior” which includes many of these forms, and shows that the class of functions so defined has a coherent and non-trivial mathematical theory. Data from the stock production database assembled by Ransom Myers is used to illustrate why such generalizations have been sought in the first place, and to highlight the difficulties in choosing between model forms on purely empirical grounds. Special attention is given to a parametric family of functions within this class, here called “θ-BH” functions. These functions cover a broad range of shapes, including both the Beverton-Holt and hockey stick functions, and share useful properties with these two widely-used models.

from, a deep analysis of the ecology of an organism. In practice, it is more common to 11 use a stock production function as a substitute for such an analysis. That is, one argues 12 on qualitative grounds that the "true" relationship should be at least approximately of 13 a certain algebraic form, and selects parameters for this form by some process that 14 side-steps the need for detailed theories of individual behavior, bioenergetics, and so on. 15 One commonly used form is the "Beverton-Holt" function [1]: This form is appealing for a number of reasons, not least of which is that it has only two 17 parameters, both of which can be given plausible interpretations as "real" quantities 18 (this will be considered more carefully in Sections 4 and 5). It has become one of the 19 most commonly used of all stock-production models. 20 The Beverton-Holt function has also been frequently criticized, however. Aside from 21 reservations about the usefulness of stock-production theory in general, a recurrent 22 complaint is that theory or empirical data suggest that Beverton-Holt has the wrong 23 qualitative behavior for the situation in hand. 24 For example, while it is trivially true that there must be some theoretical upper 25 limit on the population, this might not have any practical meaning for the population 26 under consideration, which might be better described by the Shepherd function 27 F (X) = rX/(1 + bX γ ) with γ < 1 [2]. Or one might expect production to attain a 28 maximum at some finite stock and then decrease, like the Ricker function 29 F (X) = rX exp(− rX eK ) [3]. If Allee effects are a concern, the Thompson function 30 F (X) = rKX δ /(rX δ + K) might be appropriate [4]. Some of the possibilities are  In many situations, biological considerations lead to the expectation that production 44 should be nearly proportional to stock when habitat is freely available, and nearly 45 constant when some dimension of habitat is fully utilized. In particular, this is the case 46 when population dynamics is driven principally by contest competition [5,6]. These 47 qualitative expectations can be formalized as constraints on the mathematical form of a 48 stock-production function Y = F (X):

49
(I) For some 0 < r < ∞, the graph of F approaches the line Y = rX as X goes to 0. 50 (II) For some 0 < K < ∞, the graph of F approaches the line Y = K as X goes to ∞. 51 Conventionally, r is interpreted as an "intrinsic rate of increase" or a 52 "density-independent survival," and K as the ultimate "carrying capacity" of the habitat 53 (this will be considered more carefully in Section 4). 54 These are typically are the only properties that are justifiable a priori, so F should 55 be otherwise unremarkable. A simple formulation of this is: 56 (III) F (X) is non-decreasing and concave.

57
A more stringent notion of "unremarkability" is explored in S1 Appendix. 58 For the rest of this paper, continuous functions satisfying the conditions (I), (II), and 59 (III) will be called Beverton-Holt-like. It is mathematically convenient to exclude the

Examples of Beverton-Holt-like functions 88
Some examples of Beverton-Holt-like functions are collected in Table 1 in standard form, 89 together with their defects. Fig 3 shows a sampling of these. Table 1. Examples of Beverton-Holt-like functions in standard form. To express these in the original units for stock and production, apply the substitutions x = rX/K, y = Y /K, and multiply the defects by K 2 /r.
Since all Beverton-Holt-like functions look more or less alike to the human eye, it is 91 reasonable to wonder if there is any point to this proliferation of algebraic forms. This 92 will be addressed in Sections 4 and 6; for now, it will just be noted that many of the 93 functions in Table 1 were introduced specifically to address practical problems 94 encountered with Beverton-Holt.

95
The hockey stick, Beverton-Holt, and Skellam functions have all been used by 96 multiple researchers. Moreover, less familiar functions seem invariably to include one or 97 more of these as a special or limiting cases. It is obvious that any Beverton-Holt-like 98 function is bounded above by the hockey stick of the same r and K, and it is shown in 99 S1 Appendix that any Beverton-Holt-like function which is "sufficiently boring" (in a 100 technical sense explained in the Appendix) is bounded above by the Skellam of the same 101 r and K. differential equation with which it is closely related (Section 3.3). [9] calls it the 104 "δ-power B&H", and derives it from a size-based recruitment model. It makes a fleeting 105 appearance in [12] without acquiring a name there. Despite its simplicity, it has not 106 seen much use as a stock-production model. The θ-BH functions can be obtained from 107 Beverton-Holt by a simple change of variables: Y is a θ-BH function of X if and only if 108 Y θ is a Beverton-Holt function of X θ .

109
The negative power distribution is derived in [6] as a mixture of Skellam functions. 110 Mixtures of Beverton-Holt-like functions will be discussed more generally in Section 2.5. 111 The positive power distribution is a trivial variant of this, but I haven't been able to 112 find a case in which it has been used as a stock-production relation (perhaps because it 113 is not easily fitted by the usual least-squares method).

114
The logistic hockey stick is presented in [7] as a function in three parameters α, µ, 115 and θ, most naturally expressed as a cumulative integral: The standard form of this turns out to depend only on θ, so this is a natural shape 117 parameter; when working with the standard form, it is convenient to use λ = e 1/θ + 1, 118 resulting in the (admittedly bizarre) expression in Table 1. The natural domain of λ 119 from this derivation is 1 < λ < ∞, with the limiting cases λ → 1 and λ → ∞ yielding 120 the Skellam and hockey stick, respectively. The algebraic expression in the table   121 continues to describe a Beverton-Holt-like function for the larger interval 0 < λ < ∞, 122 with a removeable singularity at λ = 1.

123
The bent hyperbola is the specialization of the more general bent hyperbola of [13] 124 to the context of Beverton-Holt-like functions. It is presented in [10] as a function in 125 three parameters β, S * , and γ: The standard forms of these turn out to depend only on the quantity 127 λ = S * / S * 2 + γ 2 /4, which is taken as the shape parameter here. The natural domain 128 the larger interval −∞ < λ ≤ 1.

132
The bent hyperbola reveals a limitation of the defect as a measure of distance from 133 the hockey stick: although λ = 1 is the hockey stick, and the convergence as λ → 1 is 134 uniform, the defect is infinite for all λ < 1. and differentiable at all but countably many points [14].

145
If F is in standard form, ϕ is a probability density on (0, ∞). This gives an 146 interesting interpretation of the defect: if ϕ is the density function associated with an f 147 in standard form, then with the special cases p = θ, q = 1 and p = 1, q = λ, respectively, of the generalized 157 t-distribution of [17].

169
Considerations of habitat heterogeneity lead naturally to such mixtures. This will be 170 pursued a bit further in S1 Appendix.

171
The "negative power" function, which includes both Skellam and Beverton-Holt as 172 limiting cases, can be exhibited as a continuous mixture of Skellam functions: where ν is a gamma distribution with expected value 1: If the term 1 + x/λ is modified to (1 + x/λ) + , the negative power function continues 175 to be an r-K function for negative values of the parameter λ, and is Beverton-Holt-like 176 when λ ≤ −1. This "positive power" function has some charm: it has the interesting (and biologically natural?) property of actually attaining the production capacity at a 178 If Y is a Beverton-Holt-like function of X, and Z is a 182 Beverton-Holt-like function of Y , then Z is also a Beverton-Holt-like function of X.

183
That is, compositions of Beverton-Holt-like functions are again Beverton-Holt-like. In 184 particular, models obtained by iterating Beverton-Holt-like functions cannot produce 185 "interesting" population dynamics, as iteration of Ricker functions famously can [19].
This property is sometimes convenient for simulation modeling [20]. It is exploited 191 systematically, for example, in the EDT framework of [21].

192
The Beverton-Holt family is not the only family with a "composition law", however. 193 If F and G are θ-BH for the same value of θ, G • F is also θ-BH for this θ, with Moreover, if F and G are hockey sticks, G • F is also a hockey stick, with Since the limit as θ → ∞ of the θ-BH functions having given r and K parameters is 196 the hockey stick with these parameters, hockey sticks can be thought of as ∞-BH 197 functions; with this convention, Eqs (9) and (11) are both special cases of Eq (10). suppose that F is closed under composition. Then there is some 0 < θ ≤ ∞ such that F 202 is the family of θ-BH functions.

203
Since this is a mathematical fact, rather than a biological one, the demonstration is 204 relegated to S2 Appendix. This section will describe this connection, and show that the θ-BH functions are 213 associated in the same way with θ-logistic differential equations.
for the abundance P , where X is the initial population and ϕ is a continuous function 220 It is shown in standard textbooks on differential equations that, under mild 222 assumptions on ϕ, Eq (12) has a unique solution for any X > 0, at least on some 223 non-trivial interval containing t = 0 [22] . If there is some interval [0, T max ) which works 224 for all X > 0, then any T strictly between 0 and T max gives rise to a stock-production 225 relation F between X = P (0) and Y = P (T ).  Perhaps the oldest population model of all is the differential equation with a constant. This appears (implicitly) already in the pioneering work of Graunt on 232 human demographics [24]. The solution, P (t) = P (0)e at , gives rise to the 233 stock-production function Y = rX, where r = e aT .

234
Malthus observed that exponential growth cannot persist indefinitely, and must 235 therefore be modified by some kind of density dependence [25]. An early mathematical 236 formulation of this is the logistic differential equation of Verhulst [26],  Typical derivations of Eq (14) have an ad hoc flavor, starting from the desired 242 qualitative behavior of ϕ and simply taking the "simplest" form that works [26,27]. alternatives have been considered from a very early date. In particular, the models for θ > 0 already appear in Verhulst's 1838 paper [26, page 116].

246
Eq (17) is sometimes called the "θ-logistic" equation. In population-modeling 247 contexts, it has been called the "Richards equation" after its appearance in [28].

Time-varying parameters 252
The models (13), (14), and (17) where p(t) is the particular solution with p(0) = 1 and r(t) = exp( The stock-production functions associated with a given ϕ as in Section 3.1 for different 264 choices of T , form a one-parameter family.

265
More generally, let ϕ be a continuous function from (0, ∞) × [0, ∞) into (−∞, ∞), 266 and suppose that there is some 0 < T max ≤ ∞ such that, for each 0 ≤ s < T max and 267 each 0 < X < ∞, the initial value problem 268 dP dt = ϕ(P, t), P (s) = X has the unique solution F s,t (X), s ≤ t < T max . Then F = {F s,t | 0 ≤ s < t < T max } is 269 a two-parameter family of stock-production functions, degenerating to a one-parameter 270 family in the autonomous case (since then F s,t = F 0,t−s ).

271
In the cases considered in Sections 3.2-3.4, this family was homogeneous. This turns 272 out to be a very special property:

273
Theorem 3. Let F be as above. If the F s,t are all r-K functions having the same 274 standard form, this form is the standard θ-BH for some 0 < θ < ∞.

275
This result is closely related to Theorem 2. A proof is given in S2 Appendix.

Beverton-Holt-like functions as population models 277
It is hard to avoid the suspicion that all this mathematics is out of proportion to the 278 original problem-that the proliferation of algebraic forms for the stock-production 279 function is more a matter of scholastic hair-splitting than practical biology.

280
Such doubts are only exacerbated by the fact that such functions are typically used 281 in a crudely empirical way. The usual goal is simply to describe a data set, with no 282 pretense that the fitted form is anything other than a conveniently simple 283 approximation to an inconveniently complex reality.

284
This section and the next will try to explain why so many "alternatives to" or

Physical interpretation of model parameters 288
In the conceptual model of Section 2, the r and K parameters correspond to "real" 289 quantities: 290 r is the "intrinsic rate of increase" one would expect to see in the absence of 291 crowding-a measure of habitat quality.

292
K is the "carrying capacity" or maximum production potential-a measure of 293 habitat quantity.

294
There are thus two sets of "r" and "K" parameters present when a Beverton-Holt-like 295 model is considered: values describing the habitat, and values describing the population 296 dynamics.

297
Some applications of stock-production modeling rely on identifying the two. This 298 can be done in either direction:

299
One can attempt to obtain information about physical habitat from population 300 data. For example, parameters obtained by fitting a stock-production model to 301 population data may be used as estimates of physical values, in the course of 302 setting harvest levels [29] or estimating extinction risks [30].

303
One can attempt to obtain information about population dynamics from habitat 304 data. For example, stock-production models parameterized with physical values 305 (from survival experiments, habitat mapping, etc.) may be used to game 306 management or restoration alternatives [20].

307
What makes this identification dangerous is that the r and K parameters by 308 definition concern properties of the fitted curve at the fringes of the data (strictly  Of course, it is to be expected that using the wrong model will produce incorrect 320 results-the Beverton-Holt function is not special in this regard. However, if it is typical 321 for "true" stock-production curves to sit above the Beverton-Holt curve having the same 322 asymptotes, then fitting Beverton-Holt functions will typically over-estimate r and K. 323 And this is precisely what is asserted in in the cited papers.

324
At the root of these charges is the very slow rate at which the Beverton-Holt 325 function rises toward its asymptote, which corresponds to a seemingly inefficient use of 326 resources. For example, Eq (1) predicts that even when the habitat is 100% over-seeded 327 (in the sense that X = 2K/r, the stock which would yield a production of 2K in the 328 absence of density-dependence), fully one-third of the productive potential of the 329 habitat will remain unexploited (in the sense that Y = 2K/3). Whether this is 330 reasonable or not of course depends on how organisms actually interact with the 331 habitat, and with one another, in the situation at hand. But as a generic assumption, to 332 be used in the absence of population-specific detail, it is at least open to challenge.

333
One measure of just how slowly the Beverton-Holt function approaches the 334 horizontal asymptote is that the "wedge" between the curve and the asymptote has 335 infinite area. In the notation of Table 1 This section will make the discussion of Section 4 more concrete by considering some 345 actual data. Because the papers of Myers et al. are so cogent, and because these "meta-analytic methods to combine results across many populations" [33], they seem 348 particularly well suited to the present purpose. The focus will be on populations of 349 Coho salmon Oncorhynchus kisutch discussed in [7]. salmon, Oncorhynchus kisutch, in which "stock" is female spawners and "production" is 353 outmigrant smolts, arranged very subjectively by shape.

354
The fitted curves will be discussed below. Disregarding these for the moment, 355 several features of these data are worth remarking.

356
First, all the individual datasets are quite small. The largest consists of 26 points.

357
This is typical of stock-production data: monitoring programs that follow a consistent 358 methodology for a quarter of a century are unfortunately (albeit understandably) rare. 359 Next, the data are quite noisy. Again, this is typical of stock-production data: 360 environmental drivers (flow, water temperature, food availability) can vary dramatically 361 from year to year, and the estimates of both stock and production usually come from 362 sampling, with large uncertainties.

363
Because the datasets are small and noisy (or as Ransom Myers put it, "nasty, 364 brutish, and short" [33]), it is unrealistic to expect the data to point unambiguously to 365 a particular algebraic form. One might decide that a dataset seems "Beverton-Holt-ish", 366 or "Ricker-ish", but it seems hopeless to choose between, say, the logistic hockey stick 367 and bent hyperbola, on empirical grounds alone. It will be seen that there are good 368 reasons for considering models with a degree of freedom beyond r and K, but one more 369 "shape" parameter is probably all that can be accommodated. The observed stock-production pairs (X i , Y i ), 1 ≤ i ≤ n, are assumed to be related via a 377 stock-production function F from some parametric family F (·, p). Since the purpose 378 here is simply illustrative, a simple multiplicative error structure is assumed: where the ε i are i.i.d. normal deviates with mean 0 and variance σ 2 .

380
(Stock-production datasets are typically time-series, with X i in turn a function of 381 earlier Y j , j < i, and autocorrelation should be taken into account. Furthermore, both 382 X i and Y i are usually themselves estimated from sampling data, and hence have error 383 structures of their own. A more general framework here is state-space modeling [35].)

384
Eq 22 gives rise to the log-likelihood The maximum-likelihood estimate for the parameters is (p,σ 2 ), wherep minimizes the 386 residual sum of squares andσ 2 = 1 n RSS(p).

388
The Bayesian Information Criterion is then 389 BIC(F ) = n logσ 2 + n log(2π) + (k + 1) log(n) (25) where k is the length of p. Given an alternate model form The quantity exp(− 1 2 ∆ BIC(F, G)) is interpretable as an evidence ratio.

392
The models considered will all be θ-BH: The present concern is how the performance of this model depends on θ. That is, rather 394 than treating θ as a fitting parameter, a range of values for θ are considered, and the  Dependence of fitted models on the shape parameter θ. By at least one widely-used criterion, a broad range of θ values are more-or-less equally compatible with the data. The fitted r and K parameters however, vary considerably between models (note the logarithmic scale).

404
In most cases, the fitted curves look very similar to one another. This subjective  To put it more dramatically, it is not possible to estimate r and K very well without 411 committing to a particular model form, and the data themselves are of little help in Essentially the same phenomenon is analyzed in [36], where a discretized version of 416 the θ-logistic differential equation discussed in Section 3 is fitted to real and simulated 417 population time-series. These authors show that treating θ as a fitting parameter on par 418 with the (appropriate analogs of) r and K is very unstable, in the sense that small 419 perturbations of the data can result in large changes to the parameter estimates, and in 420 the case of simulated series, that model fitting cannot reliably recover the "true" values 421 used to generate the series in the first place.  Evidence-ratio curves for θ-BH fits to stock-recruitment time-series for all stocks from selected orders from the Myers database.

429
All but a handful of these curves are monotone on the entire interval 0.5 ≤ θ ≤ 4.

430
That is, if θ is treated as a parameter on a par with r and K, to be estimated by 431 maximum likelihood, the fitting routine wants to drive the model to one of the two 432 degenerate forms θ → ∞ or θ → 0.

433
The curves which do have a local maximum are very flat: in only eleven cases is the 434 evidence ratio at some intermediate θ more than 20% higher than at both endpoints.

435
For 170 (31%) of the stocks, the curves are "extremely flat," in the sense that they vary 436 no more than 1% across the full range of θ. distributions are used to explore other kinds of data.

460
Almost any data analysis these days is likely to include some linear regressions, 461 complete with R 2 statistics, even when there is no particular reason to expect data to 462 come from a truly linear relationship, or for the errors to be independent fits and stock-production fits are typically put. The mean of a normal fit, or the slope 469 of a linear fit, are attempts to capture some property of the "center" or "main body" of 470 the data. The r and K parameters of a Beverton-Holt-like function, however, concern 471 properties of the fitted curve at the fringes of the data, or in many cases well beyond 472 them.

473
It is well-known that interpolation is a much safer process than extrapolation, and 474 when a fitted model is to be used in an "interpolatory" way, for example, to make 475 short-term forecasts, or to explore the implications of modest changes to habitat under 476 hypothetical conditions broadly similar to historical conditions, the precise form of the 477 fitted model is unlikely to be of great importance.

478
However, to estimate a true intrinsic productivity or ultimate carrying capacity from 479 passively-observed stock-prodution data, by fitting any kind of stock-production model, 480 is always an extrapolation. Statistics is not magic, and information which is not present 481 in the data to begin with cannot be extracted from it by mathematical manipulations. 482 In the case of stock-production:

483
If the data do not include cases of low stock densities, the data contain no 484 intrinsic information about productivity at low densities. Any estimate of intrinsic 485 productivity will be driven primarily by a priori assumptions about the model 486 form.

487
If the data do not include cases of production near the carrying capacity, the data 488 contain no intrinsic information about carrying capacity. Any estimate of carrying 489 capacity will be driven primarily by a priori assumptions about the model form. 490 Even if the data do include low or high stock densities, the "best" fit of a model 491 form to the overall data need not reflect the actual behavior at these extremes