On stochastic logistic population growth models with immigration and multiple births

ARTICLE IN PRESS

Theoretical Population Biology 65 (2004) 89–104

http://www.elsevier.com/locate/ytpbi

On stochastic logistic population growth models with immigration and multiple births James H. Matisa,
and Thomas R. Kiffeb a

Department of Statistics, Texas A&M University, 3143 TAMU, College Station, TX 77843-3143, USA b Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA Received 20 May 2003

Abstract This paper develops a stochastic logistic population growth model with immigration and multiple births. The differential equations for the low-order cumulant functions (i.e., mean, variance, and skewness) of the single birth model are reviewed, and the corresponding equations for the multiple birth model are derived. Accurate approximate solutions for the cumulant functions are obtained using moment closure methods for two families of model parameterizations, one for badger and the other for fox population growth. For both model families, the equilibrium size distribution may be approximated well using the Normal approximation, and even more accurately using the saddlepoint approximation. It is shown that in comparison with the corresponding single birth model, the multiple birth mechanism increases the skewness and the variance of the equilibrium distribution, but slightly reduces its mean. Moreover, the type of density-dependent population control is shown to influence the sign of the skewness and the size of the variance. r 2003 Elsevier Inc. All rights reserved. Keywords: Birth–death process; Population growth model; Differential equations; Moment closure methods; Saddlepoint approximations

1. Introduction Population growth models are central in modern ecological theory. The classic Verhulst–Pearl logistic equation incorporates density-dependent growth and has a rich history dating back to 1838 (see, e.g., Renshaw, 1991; Nasell, 2001). Its full solution, and especially its equilibrium value called the carrying capacity, are simple mathematical expressions. This simplicity, no doubt, has contributed to its widespread use in practice, particularly in models for managing ecological populations (see, e.g., Pielou, 1977). Two specific examples used subsequently in this paper concern managing fox populations (Anderson et al., 1981) and badger populations (Anderson and Trewhella, 1985) in Europe; other examples include Pech and Hone (1988), Coyne et al. (1989), and Roberts (1996). The development of a stochastic logistic model is reviewed in Nasell (2001). This well-known model


Corresponding author. Fax: 979-845-3144. E-mail address: [email protected] (J.H. Matis).

0040-5809/$ - see front matter r 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.tpb.2003.08.003

(see, e.g., Pielou, 1977; Renshaw, 1991) is also considered in Matis and Kiffe (1996), which gives differential equations for approximate cumulant (e.g., mean, variance, skewness) functions of the population size variable. Recently, a simple analytical saddlepoint approximation for the equilibrium distribution, based on these cumulant approximations, was proposed (Matis et al., 2003), and it is expected that the simplicity and accuracy of this approximation will lead to a much wider usage of equilibrium size distributions in ecological modeling. Arguably, the logistic growth model would be more realistic biologically with the inclusion of two additional features. The first of these is immigration. The solution for the deterministic logistic growth model with immigration is much more involved mathematically (see, e.g., Banks, 1994), though a relatively simple solution still exists for its carrying capacity. The stochastic model with immigration is investigated by Matis and Kiffe (1999), which gives differential equations for its cumulant functions. Though the case for immigration is often compelling ecologically, its more involved mathematical solutions have restricted its use

ARTICLE IN PRESS 90

J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

in practice. For example, in the previously cited applications immigration is logically present, however the simpler, tractable ordinary logistic growth model without immigration was sufficient for the broader objectives of the studies. The second feature that would add realism, particularly for small populations, is the inclusion of a multiple birth mechanism. The aforementioned deterministic models are based on the notion of continuous fractional increases, and the present stochastic models assume unit births. Multiple births have previously been included in some stochastic models with linear dynamics (see, e.g., Matis and Kiffe, 2000), but we are not aware of the inclusion of this feature in logistic growth models. The principal motivation of this paper is to develop and investigate the stochastic logistic model with immigration and with multiple births. Section 2 reviews the deterministic ordinary logistic model with immigration, and illustrates its use with the classic badger and fox population models. Section 3 outlines the corresponding stochastic model with single births, and illustrates the use of relatively simple and accurate approximations for its population size distributions. Section 4 develops the new model with multiple births, and also illustrates the user-friendly distribution approximations. Section 5 then investigates the effect of the multiple birth generalization on the equilibrium size distribution. Concluding discussion is given in Section 6.

independent of population size, i.e., IðNÞ ¼ I;

ð3Þ

also called constant immigration. Defining a ¼ a1  a2

and

b ¼ b1 þ b2 ;

ð4Þ

the differential equation for N which combines (1)–(3) is N’ ¼ I þ aN  bN 2 : ð5Þ The closed form solution of (5) is NðtÞ ¼ fa þ bð1  debt Þ=ð1 þ debt Þg=2b;

ð6Þ

where b ¼ ða2 þ 4bIÞ1=2 ; g ¼ ð2bN0  aÞ=b; d ¼ ð1  gÞ=ð1 þ gÞ with N0 denoting the initial population size. This solution is given in an equivalent form in Banks (1994, p. 117). The equilibrium value, or carrying capacity, for this special case is K ¼ ða þ bÞ=2b:

ð7Þ

In the special case of no immigration, with I ¼ 0; solution (6) simplifies considerably, and the carrying capacity reduces easily from (7) to K ¼ a=b

ð8Þ

which is a well-known result for the ordinary logistic model (Pielou, 1977). 2.2. Examples

2. Deterministic logistic growth models with immigration 2.1. Assumptions and model solution Consider a population of size N: For subsequent convenience, we describe changes due to reproduction and mortality separately, and assume population birth and death rates lðNÞ ¼ a1 N  b1 N 2 ;

ð1Þ

mðNÞ ¼ a2 N þ b2 N 2

ð2Þ

respectively, where a1 ; a2 ; b1 ; b2 40: The ai are interpreted as the ‘‘intrinsic rates’’ and the bi as the ‘‘crowding coefficients.’’ For simplicity, the assumed birth and death rate functions are quadratic, but the theory would be easily extended to any polynomial rate function, as outlined in Matis and Kiffe (1999). Nasell (2001, 2003) uses an alternative parameterization, however the above model formulation is useful for subsequent purposes in which b1 and b2 are involved directly. The population is also assumed to have an immigration rate (number of individuals/time) which is

We previously applied these models to describe the spread of muskrats in the Netherlands (Matis and Kiffe, 1999), with parameter values estimated from data. For present purposes, however, we consider applications of the theory to fox and badger populations, for which simple deterministic logistic models without immigration were sufficient previously for the intended purpose. We will subsequently extend each model to include immigration and multiple births. 2.2.1. Badgers The population biology for badgers in Europe is described in an extensive monograph by Anderson and Trewhella (1985). Their intrinsic rate estimates for southwest England are a1 ¼ 0:61 and a2 ¼ 0:37: A density of five to eight adults=km2 is regarded as typical in England, which for convenience we scale and set to K ¼ 60=ð10 km2 Þ: The authors set b2 ¼ 0; noting that ‘‘the principal density-dependent check on growth appears to act on fecundity as opposed to y mortality’’ (Anderson and Trewhella, 1985, p. 338). Their model has no immigration, hence it follows from (4) and (8) that b1 ¼ 0:004:

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

100

60

91

I=10

50 40

80 Population Size

Population Size

I=6 I=2

30 20

I=3

60 40 20

10 10

20

30 Time (yrs)

40

50

10

20 30 Time (yrs)

40

50

Fig. 1. Solutions for deterministic population size function, NðtÞ; for family of badger models with two levels of immigration, I:

Fig. 2. Solutions for deterministic population size function, NðtÞ; for family of fox models with two levels of immigration, I:

We investigate two cases which add immigration to the model. Consider first the case of I ¼ 6; which, at 10% of K; we call ‘‘high immigration’’. We propose a model where the intrinsic rates and the carrying capacity are fixed, and the density-dependent rate b is adjusted to compensate for immigration. Solving for b from (5), the equilibrium size, K; satisfies

3. Stochastic logistic growth model with immigration

b ¼ ðI þ aKÞ=K 2 :

ð9Þ

3.1. Assumptions and general approach Let the population size N now be a random variable. We assume that the probabilities of a single birth, single death, and single immigrant, respectively, for a population of size N ¼ n in an infinitesimal time interval, Dt; are ln Dt; mn Dt and IDt where
ln ¼

a1 n  b 1 n2

for noa1 =b1 ;

0

otherwise;

ð10Þ

For I ¼ 6 with the assumed intrinsic rates and K ¼ 60; Eq. (9) yields b ¼ 0:005667; which under our assumption (of b2 ¼ 0) is also b1 : For a second case, consider I ¼ 2 which, at 3.3% of K; we classify as ‘‘low immigration’’. Substitution into (9) yields b1 ¼ 0:004556: The solutions, from (6), for the transient population sizes of these two cases, assuming initial value N0 ¼ 6; are given in Fig. 1.

mn ¼ a2 n þ b2 n2

2.2.2. Foxes The parameter estimates for foxes in Europe come from Anderson et al. (1981). The intrinsic rates are a1 ¼ 1:00 and a2 ¼ 0:50: Typical carrying capacities are noted to range from 0.1 to 4 foxes=km2 ; which we scale and set to K ¼ 100 foxes=ð50 km2 ). The authors set b1 ¼ 0; assuming that ‘‘density dependence acts primarily on mortality’’ (Anderson et al., 1981, p. 766). This illustration with densitydependent deaths provides a nice contrast to the badger illustration with density-dependent births. It follows from (8) that b2 ¼ 0:5=100 ¼ 0:005 for the model without immigration. Consider again two immigration cases. For high (10%) immigration, with I ¼ 10; and a and K fixed at 0.5 and 100, one has from (9), b2 ¼ 0:006: Low (3%) immigration yields I ¼ 3 and b2 ¼ 0:0053: The transient solutions using (6) for the population sizes of these two cases, with N0 ¼ 10; are given in Fig. 2.

Clearly, assumptions (10)–(12) are analogs in the stochastic model of assumptions (1)–(3) in the deterministic model. We review briefly some general approaches for investigating this model. First, methods for obtaining the exact population size distribution are reviewed. The Kolmogorov equations are used to find the transient distributions, and the special case of obtaining the equilibrium distributions directly is outlined. Then, the ‘‘moment-closure’’ procedure, which yields approximate solutions for low-order cumulant functions, is presented. These approximations are shown to be very accurate for the muskrat population illustration with single births in Matis and Kiffe (1999). This paper extends the methodology by using the cumulant approximations to form Normal and saddlepoint approximations of the population size distributions. These methods are illustrated for the assumed badger and fox populations.

ð11Þ

and In ¼ I:

ð12Þ

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

92

3.2. Direct solutions for population size distributions

respectively, are

Consider first the equilibrium size distributions. In the model without immigration (i.e., with I ¼ 0), ultimate extinction is certain, as the population size is bounded above and n ¼ 0 is an absorbing barrier. However, Nasell (2001) outlines the calculation of a ‘‘quasi-stationary’’ distribution of population size conditional upon N40: As a contrast, the mode with immigration has a true equilibrium distribution. It may be obtained directly, numerically, from the recurrence relationships

k1 ðtÞ ¼ m1 ðtÞ;

mn Pn ¼ ðln1 þ IÞPn1

ð13Þ

for nX1: These equations, also called the ‘‘balance equations,’’ imply that the probability of population size n with a decrease to n  1 is equal to the probability of size n  1 with an increase to n; thereby ruling out any ‘‘drift’’ in mean size. Nasell (2001) shows that this approach is valid only for the case I40; for which P0 is positive. One can find the transient distributions at any time t; and of course the equilibrium distributions as well, of the population size NðtÞ from the Kolmogorov equations (see, e.g., Bailey, 1964). This set of basic differential-difference equations is given for the present model in Matis and Kiffe (1999). The population size NðtÞ with immigration is unbounded and hence the system of Kolmogorov equations is infinite. Approximate solutions for the transient distributions may be obtained by ‘‘truncating’’ the system at some large upper value, say U: In practice, U is chosen to be so large that the approximation error is negligible, and the resulting probability functions are deemed to be exact for all practical purposes. 3.3. Solutions using cumulant functions Consider now an alternative approach. The (uncorrected) moments, denoted for simplicity as mi ðtÞ for i ¼ 1; 2; y; of the distribution at time t are defined as the expectations of the powers of NðtÞ; i.e., N X mi ðtÞ ¼ E½N i ðtÞ ¼ ni Pn ðtÞ; ð14Þ n¼0

where Pn ðtÞ denotes Prob½NðtÞ ¼ n : It is convention to use m to denote both a moment, as in (14), and a death rate, as in (11). However, in this paper only the moments, mi ðtÞ; are time-dependent functions, and hence it should always be clear in the context of the problem which interpretation is intended. The moment functions may be transformed to cumulant functions, denoted ki ðtÞ; (Smith, 1995; see also Johnson et al., 1992). In particular, the first three cumulant functions, which are the mean, variance and skewness functions

k2 ðtÞ ¼ m2 ðtÞ  m21 ðtÞ; k3 ðtÞ ¼ m3 ðtÞ  3m1 ðtÞm2 ðtÞ þ 2m31 ðtÞ:

ð15Þ

In order to derive direct differential equations for the cumulant functions, consider first the ratio r ¼ a1 =b1 : We define u as the smallest integer greater than or equal to r: Hence it follows from (10) that lu1 40 and lu ¼ 0: We define now a slightly modified birth–death process with the quadratic function ln ¼ a1 n  b1 n2

ð16Þ

for all n: Under this modification, the birth rate ln is negative for n4u: However, the modification would have no measurable effect on the distribution of NðtÞ if the probability of NðtÞ4u is negligible. We suggest in Matis and Kiffe (1999) using this modification only for models which satisfy the condition Prob½NðtÞ4u

o0:01 for all t40: This may be checked in practice by calculating the standard Normal variate z ¼ ðu  1=2 KÞ=k2 ; and using the modification only for models with z42:33: It is easy to verify subsequently that all of our proposed badger and fox models easily satisfy this criterion and hence are well below the 1% threshold. Hence the modification in (16) has negligible effect. Consider now the moment generating function (mgf), defined as Mðy; tÞ ¼

N X

mi ðtÞyi =i!:

ð17Þ

i¼0

We show in the Appendix that a partial differential equation for the mgf of the present birth–immigration– death process with the rates in (11), (12) and the modification in (16) is @M @M ¼ ½ðey  1Þa1 þ ðey  1Þa2

@t @y þ ½ðey  1Þðb1 Þ þ ðey  1Þb2

@2M @y2

þ ðey  1ÞIM:

ð18Þ

The cumulant generating function (cgf) is defined as Kðy; tÞ ¼ log Mðy; tÞ;

ð19Þ

which upon substitution into (18) yields @K @K ¼ ½ðey  1Þa1 þ ðey  1Þa2

@t @y

 2 # 2 @ K @K þ ½ðey  1Þðb1 Þ þ ðey  1Þb2

þ 2 @y @y þ Iðey  1Þ:

"

ð20Þ

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

The cumulants are defined as the coefficients in the power series expansion X Kðy; tÞ ¼ ki ðtÞyi =i!: ð21Þ iX1

Substituting (21) into (20) and equating powers of y yields the following set of differential equations for loworder cumulants: k’ 1 ðtÞ ¼ I þ ½a  bk1 ðtÞ k1 ðtÞ  bk2 ðtÞ; k’ 2 ðtÞ ¼ I þ ½c  dk1 ðtÞ k1 ðtÞ þ ½2a  d  4bk1 ðtÞ k2 ðtÞ  2bk3 ðtÞ; k’ 3 ðtÞ ¼ I þ ½a  bk1 ðtÞ k1 ðtÞ þ ½3c  b  6dk1 ðtÞ  6bk2 ðtÞ k2 ðtÞ þ ½3a  3d  6bk1 ðtÞ k3 ðtÞ  3bk4 ðtÞ

mean is close to the unique deterministic solution. The simple equilibrium solutions derived in this manner are denoted k
i ; i ¼ 1; y; 3: The accuracy of distribution approximations based on these cumulant approximations will be explored. If the skewness is small, the distribution will be approximated at any time of interest using the Normal distribution with the approximate value of the first two cumulants. If, on the other hand, the skewness indicates a noticeable lack of symmetry, the approximations for the first three cumulants may be substituted into the saddlepoint point approximation based on third-order truncation, expressed in analytical form (see Renshaw, 1998) as f ðxÞ ¼ ð4p2 cÞ1=4 expfðk32  3k2 c þ 2c3=2 Þ=ð6k23 Þg;

ð22Þ

with c ¼ a1 þ a2

93

ð24Þ where

and

d ¼ b1  b2 :

ð23Þ

Special cases of (22) have been given previously. In particular, the equations for the mean and variance, k’ 1 ðtÞ and k’ 2 ðtÞ; were derived in a different way in Matis and Kiffe (1999), and all three equations are given for the case I ¼ 0 in Matis et al. (2003). Exact solutions of these equations are not available, as the equation for any specific jth order cumulant function involves terms up to the ð j þ 1Þth order. Approximate solutions are available using ‘‘moment-closure’’ methods (see, e.g., Whittle, 1957) in which the differential equations are solved with higher-order cumulants set equal to 0. Our experience indicates that accurate approximations of cumulant functions of any specific order, say i; may be obtained by solving the set of ði þ 1Þ equations with all cumulants of order above i þ 1 set to 0. Therefore, in order to solve for the first three cumulant function approximations, we obtain an equation for k’ 4 ðtÞ; using Mathematica (Wolfram, 1991), by further expansion of (20). The set of four equations are then solved setting k5 ðtÞ ¼ 0: These approximations are denoted k# i ðtÞ; for i ¼ 1; y; 3: Often the solution for only the equilibrium distribution is sufficient for the objectives at hand. Equilibrium cumulant approximations are easily available by setting the derivatives equal to 0, and solving the system of equations numerically. The accuracy of this procedure for the case I ¼ 0 has been previously investigated by Nasell (2001) and Matis and Kiffe (2000). The application of the approach for I40; based on only the three equations in (22) with k4 ðtÞ ¼ 0; will follow. This approach of finding the roots for (22) yields a number of vector solutions, however most are inadmissible due to a negative mean or variance. In the event of multiple admissible solutions, our approach is to choose the vector solution (in our experience the only one) whose

c ¼ k22 þ 2k3 ðx  k1 Þ:

ð25Þ

One limitation of this simple form of the approximation is that f ðxÞ is defined only for x values for which cX0: However, this saddlepoint approximation yields a completely general (i.e., distribution-free) approximation, at any time of interest, for at least the tail of the population size distribution in the direction of the skewness. This tail is usually the one of primary interest in ecology. 3.4. Examples In the examples that follow, approximate cumulant functions obtained from (22) using moment-closure methods will be compared with the corresponding exact functions obtained from the Kolmogorov equations. Two measures will be used. Max is defined for the mean and variance as the maximum percent approximation error over time, i.e., Max ¼ maxf100jk# i ðtÞ  ki ðtÞj=ki ðtÞg over the range 0oto100: The skewness functions have a zero, hence for this cumulant Max is defined as the larger of the error rates at the two extremes of the skewness function. Eq is defined as the percent approximation error of an equilibrium cumulant, i.e., Eq ¼ 100jk
i  ki ð100Þj=ki ð100Þ: The accuracy of the Normal and saddlepoint approximations of the population size distributions will also be investigated. 3.4.1. Badgers The comparative solutions for the first three exact cumulant functions to time t ¼ 50 years of the badger population sizes, with starting value N0 ¼ 6 are illustrated in Fig. 3 for the two cases. Table 1A lists the numerical values for these cumulants at t ¼ 100; by which time both cases are in effective equilibrium, as apparent in Fig. 3. The index of skewness, a relative

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

94

I=6

50

300

120

200

I=2

Variance

I=2

30 20

I=2

80 60

I=6

40

10

Skewness

100

40 Mean

400

140

60

0

10

20 30 Time (yrs)

40

I=6

-300 0

50

0 -100 -200

20

(A)

100

10

20 30 Time (yrs)

(B)

40

50

0

(C)

10

20 30 Time (yrs)

40

50

Fig. 3. Exact cumulant functions for family of single birth badger models with two immigration levels: (A) mean value, (B) variance, and (C) skewness.

Table 1 Equilibrium sizes of cumulants and of indices of skewness for single birth stochastic models Immigration Level Cumulant

10%

3.3%

ðAÞ Badgers k1 k2 k3 g

59.33 51.17 31.32 0.09

58.87 74.36 57.26 0.09

Cumulant

10%

3%

ðBÞ Foxes k1 k2 k3 g

98.64 157.3 183.6 0.09

98.22 185.11 198.52 0.08

Table 2 Percent errors of cumulant approximations at equilibrium (Eq) and maximum error over time (Max) of cumulant function approximations for single birth stochastic models Immigration Level Cumulant

10%

3.3%

ðAÞ Badgers k1 k2 k3

Max

Eq

Max

Eq



















1.6

0.2 1.4




0.2

Cumulant

10%

1.0 3%

ðBÞ Foxes k1 k2 k3 n

Max

Eq

Max

Eq

























0.3

0.7

0.8

0.8

Denotes o0:1%:

(dimensionless) measure defined as g ¼ k3 =ðk2 Þ3=2 ;

ð26Þ

is also listed for each case. Table 2A gives the error rates comparing the exact solutions to the corresponding

approximate solutions based on the Max and Eq measures. Notable among the findings from these tables are the following: 1. The mean value functions increase monotonically, with equilibrium values of k1 ¼ 59:3 and 58.9 for high and low immigration, respectively. The two variance functions rise rapidly initially to respective peaks, and then taper off to equilibrium values of k2 ¼ 51:2 and 74.4. The skewness functions asymptote to k3 ¼ 31:3 and 57:3; both of which are relatively small with indices jgjo0:1: 2. It is apparent in Table 2A that the moment-closure approximations from (22) are very accurate, with Maxo0:2 for the mean and variance functions and Maxo1:5 for the skewness functions. Hence the approximate solutions are visually indistinguishable from the exact solutions in Fig. 3, and the equilibrium approximations are, for all practical purposes, those given in Table 1. 3. Eq is larger in magnitude than Max for the skewness in the high immigration case. This anomaly is possible because Max for skewness functions is based on the error at only two extreme skewness values rather than the whole range. Fig. 4 illustrates the exact equilibrium size distributions, and also the corresponding distribution approximations based on the Normal distribution with the cumulant approximations. The accuracy is remarkable for both cases, due to their low indices of skewness. The small negative skewness is just barely visible in each graph. The saddlepoint approximations were also calculated, and in each case they fit the exact distribution perfectly from a visual standpoint. For the sake of simplicity, however, the Normal approximation seems adequate for all practical purposes. These approximations may also be used for the transient size distributions of the immigration models. The maximum magnitudes for the indices of skewness

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

0.1

0.1 I=2

Probability

95

I=6

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02 20

40

60

80

100

20

40

60

80

100

Size Fig. 4. Exact and Normal approximation of equilibrium size distributions for single birth badger models with two immigration levels.

3.4.2. Foxes Fig. 6 illustrates the comparative exact cumulant functions of the fox population sizes. Table 1B gives the solutions for the cumulants at t ¼ 100; by which time the distributions are in effective equilibrium. Table 2B lists the Max and Eq error rates comparing the exact solutions in Fig. 6 to the corresponding approximate solutions. The results are similar qualitatively to those for the badger models. The following are of interest.

0.05 0.04 0.03 0.02

Probability

0.01

20

40

60

80

(A)

Fig. 5. Exact and approximate size distribution at time of maximum relative skewness for single birth badger model with low immigration, I ¼ 2: (A) Normal approximation and (B) saddlepoint approximation.

1. Both cases have increasing mean value functions, with the equilibrium values of 98.64 for high immigration and 98.22 for low immigration. The variance functions have the same shape as before, rising rapidly to some peak, and then declining gradually to equilibrium values, which are 157.3 for high immigration and 185.11 for low immigration. The skewness values are relatively small ðgo0:1Þ for both cases. 2. The approximations for all three cumulant functions are very accurate ðMaxo1:0Þ for both cases and hence for all practical purposes their exact curves in Fig. 5 also represent the approximate solutions. The approximations for the equilibrium cumulant values are also very accurate ðEqo1:0Þ for both cases.

over time are 0.25 for low and 0.18 for high immigration. This ‘‘worst-case’’ scenario for low immigration occurs at t ¼ 4:6: Fig. 5 illustrates the exact size distribution at this time, together with the Normal and saddlepoint approximations based again on cumulant approximations. The graphs for this extreme case suggest that though the Normal approximation may be adequate for most practical purposes at any given time, the saddlepoint approximation represents a definite improvement and is quite accurate over its entire admissible range.

Fig. 7 illustrates the exact equilibrium size distributions together with their Normal approximations. These Normal approximations are very accurate even in the tails due to the small indices of skewness. The small positive skewness is visible, but tiny. The maximum magnitudes for the indices of skewness over time are 0.28 for high immigration and 0.29 for low immigration. This value for low immigration occurs at t ¼ 3:1; and Fig. 8 illustrates the Normal and saddlepoint approximations of the population size distribution at this time. As before with the badgers, the Normal approximation may be adequate at

0.05 0.04 0.03 0.02 0.01

20 (B)

40 Size

60

80

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

96 100 I=10 80

350

1000

300

750

I=3

500

40

Skewness

60

Variance

Mean

250 I=3

200 150

I=10

0

10

20 30 Time (yrs)

40

-250 -500

50

-750

50

0

10

20 30 Time (yrs)

(B)

40

I=10

0

100 20

(A)

250

50

I=3

0

10

(C)

20 30 Time (yrs)

40

50

Fig. 6. Exact cumulant functions for family of single birth fox models with two immigration levels: (A) mean value, (B) variance, and (C) skewness.

0.05

0.05

Probability

I=3

I=10

0.04

0.04

0.03

0.03

0.02

0.02

0.01

0.01 80

100

120

80

140

100

120

140

Size Fig. 7. Exact and Normal approximation of equilibrium size distributions for single birth fox models with two immigration levels.

any given time, yet the saddlepoint approximation is a clear improvement when the index of skewness is relatively large.

0.05

3.5. Summary for single birth models

0.03

Previous research in Matis and Kiffe (1999) has shown that the cumulant function approximations based on moment-closure procedures may be very accurate for the single birth model with immigration. The present additional real-world examples support this key finding. The results also show that the simple saddlepoint approximation of the population size distribution, at any given time, is very accurate over its entire range. When the skewness is relatively small, e.g., for the equilibrium distributions, the even simpler Normal approximation seems sufficiently accurate for all practical purposes.

0.02

0.04

Probability

0.01

(A)

20

40

60

80

100

20

40

60 Size

80

100

0.05 0.04 0.03 0.02

4. Stochastic logistic growth model with immigration and multiple births 4.1. Assumptions In mimicking the classic deterministic growth models, virtually all stochastic analogs have assumed unit births.

0.01

(B)

Fig. 8. Exact and approximate size distribution at time of maximum relative skewness for single birth fox models with low immigration, I ¼ 3: (A) Normal approximation and (B) saddlepoint approximation.

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

A principal objective of this paper is to enhance the biological realism of the stochastic model by incorporating a multiple birth mechanism, and then to investigate the properties of this model. Multiple births may be incorporated into continuous time stochastic models as follows. Let the probability of a litter of size i at the time of birth be denoted pðiÞ: We assume that immigration still occurs as single, independent units. It would be easy to generalize that, but unit immigration seems more natural for the present and most other biological examples. The birth assumptions in time interval Dt are hence assumed to change from (10) and (12) to: Probfincrease of single individual given N ¼ ng
½I þ pð1Þða1 n  b1 n2 Þ Dt for noa1 =b1 ; ¼ ð27Þ IDt otherwise: Probfincrease of i41 individuals given N ¼ ng
pðiÞða1 n  b1 n2 ÞDt for noa1 =b1 ; ¼ 0 otherwise:

ð28Þ

The death assumption restated from (11) is Probfdecrease of single individual given N ¼ ng ¼ ða2 n þ b2 n2 ÞDt:

ð29Þ

4.2. Solutions using cumulant functions A partial differential equation for the mgf Mðy; tÞ follows using the rules in the Appendix. Modifying the birth assumptions to be quadratic equations as in (16), the equation for Mðy; tÞ is
 @M X iy @M @2M ¼  b1 ðe  1ÞpðiÞ a1 þ ðey  1Þ 2 @t @y @y i
 @M @2M þ b2 a2 ð30Þ þ ðey  1ÞIM @y @y2 which generalizes (18). Transforming to the cgf, Kðy; tÞ; using (19) one has ( "  2 #) @K X iy @K @2K @K ¼  b1 ðe  1ÞpðiÞ a1 þ 2 @t @y @y @y i ( " #)   @K @2K @K 2 y þ b2 þ ðe  1Þ a2 þ @y @y @y2 þ ðey  1ÞI:

ð31Þ

Substituting the power series in (21) yields the following equations for low-order cumulants: k’ 1 ðtÞ ¼ I þ ½a1 m1  a2  b1 m1 k1  b2 k1 k1  ½b1 m1 þ b2 k2 ;

97

k’ 2 ðtÞ ¼ I þ ½a1 m2 þ a2  b1 m2 k1 þ b2 k1 k1 þ ½2a1 m1  2a2  b1 m2 þ b2  4b1 m1 k1  4b2 k1 k2  2½b1 m1 þ b2 k3 ; k’ 3 ðtÞ ¼ I þ ½a1 m3  a2  b1 m3 k1  b2 k1 k1 þ ½3a1 m2 þ 3a2  b1 m3  b2  6b1 m2 k1 þ 6b2 k1  6b1 m1 k2  6b2 k2 k2 þ 3½a1 m1  a2  b1 m2 þ b2  2b1 m1 k1  2b2 k1 k3  3½b1 m1 þ b2 k4 ;

ð32Þ

where mj denotes the jth moment of the litter size distribution, i.e., X mj ¼ i j pðiÞ ð33Þ i

for integer j: Note that in the special case of unit births, for which mj ¼ 1 for all j; (32) reduces to (22). Approximations to these cumulant functions will be obtained using moment closure procedures, as before. The accuracy of these cumulant function approximations and of the Normal and saddlepoint approximations for population size distributions will be investigated with the two examples. 4.3. Examples 4.3.1. Badgers (a) Litter size distribution and model: Neal (1977) reviews data on the litter size distribution of badgers in Europe. He notes that the number of cubs in a litter varies from one to five, and gives the following litter size distribution, which is based on a relatively large sample ðn ¼ 97Þ of litters ‘‘seen above the ground some 8–10 weeks after birth’’ (Neal, 1977, p. 222): i pðiÞ

1 0.11

2 0.51

3 0.28

4 0.08

5 0.02

The moments for these data are m1 ¼ 2:39; m2 ¼ 6:45 and m3 ¼ 19:37; which using relationships equivalent to (15) give a variance of s2 ¼ 0:74 and a skewness of k3 ¼ 0:43: For subsequent comparative purposes, the coefficient of variation, cv ¼ s=m; is 0.36 and the index of skewness using (25) is 0.68. The mean of 2.39 takes into account initial mortality during development below ground, and hence is consistent with the higher estimated mean of 2.7 cubs/female at birth given by Anderson and Trewhella (1985). In order to yield comparable mean size functions, it is apparent by comparing k’ 1 ðtÞ in (22) and (32) that the a1 and b1 parameters in the multiple birth models should be reduced by a factor of m1 from those in the single birth model. This reparameterization slows down the

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

98 60 I=6

50

300

1500

250

1000

I=2 200

30 20

100

10

50 0

(A)

10

20 30 Time (yrs)

40

I=2

150

50

0

10

20 30 Time (yrs)

0 I=6 -500

I=6

(B)

I=2

500

Skewness

Variance

Mean

40

40

50

0

10

(C)

20

30 Time (yrs)

40

50

Fig. 9. Exact cumulant functions for family of multiple birth badger models with two immigration levels: (A) mean value, (B) variance, and (C) skewness.

birth rate to compensate for the mean litter size of m1 41 when birth events occur. Accordingly, the single birth model parameters of a1 ¼ 0:61 and b1 ¼ 0:005667 and 0:009; for the two levels of immigration, were reduced by a factor of m1 ¼ 2:39 to a01 ¼ 0:255 and b01 ¼ 2:37 103 and 3:77 103 ; respectively. (b) Model solutions: The exact solutions for the first three cumulant functions of these two models are illustrated in Fig. 9, the solutions at t ¼ 100 are listed in Table 3A, and the error rates are tabulated in Table 4A. The following results are notable: 1. The exact cumulant functions have the same qualitative shape characteristics as those of the single birth model. Both cases reach an effective equilibrium, with means of k1 ¼ 58:9 and 57.9 for high and low immigration respectively. The corresponding variances are k2 ¼ 82:9 and 133.6, and the skewnesses are k3 ¼ 46:6 and 114:7; which give small indices of skewnesses ðjgjo0:1Þ: 2. The accuracy of the cumulant approximations is exceptional for the means ðMaxo0:1%Þ and very good for the variances and skewnesses ðMaxo2:5%Þ: Hence for all practical purposes, these approximate solutions are given in Fig. 9. The cumulant approximations for the equilibrium distributions also have small errors, except for the skewness approximations which have a roughly 10% error in both cases. This larger error is consistent with our experience, that often the approximation for the highest order cumulant (in this case the third order one) has large errors. Fig. 10 gives the exact equilibrium distributions and their approximations based on the Normal distribution. The Normal approximations are excellent, with barely visible negative skewness. The maximum magnitudes for the indices of skewness are 0.38 for high immigration and 0.49 for low immigration. The latter case occurs at t ¼ 5:0; at which time the population size distribution is markedly skewed to the right as illustrated in Fig. 11. The Normal approximation does poorly. The saddlepoint approximation with the cumulant approximations does

Table 3 Equilibrium sizes of cumulants and of indices of skewness for multiple birth stochastic models Immigration Level Cumulant

10%

3.3%

ðAÞ Badgers k1 k2 k3 g

58.86 82.94 46.61 0.06

57.85 133.6 114.7 0.07

Cumulant

10%

3%

ðBÞ Foxes k1 k2 k3 g

96.27 441.8 2247 0.24

94.72 552.4 2453 0.19

Table 4 Percent errors of cumulant approximations at equilibrium (Eq) and maximum error over time (Max) of cumulant function approximations for multiple birth stochastic models Immigration Level Cumulant

10%

3.3%

ðAÞ Badgers k1 k2 k3

Max

Eq

Max

Eq



















0.8

10.6

2.1 2.4

0.1 9.8

Cumulant

10%

3%

ðBÞ Foxes k1 k2 k3 n

Max

Eq







1.4 11.7

0.2 4.3

Max 2.7 21.3 58.3

Eq


0.2 4.6

Denotes o0:1%:

considerably better, though it is not admissible for most of the lower tail. 4.3.2. Foxes (a) Litter size distribution and model: Lloyd (1980) gives a litter size distribution for foxes at a different

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

0.05

0.05

Probability

0.04

I=2

0.03

0.02

0.02

0.01

0.01 40

I=6

0.04

0.03

20

99

60

80

20

100

40

60

80

100

Size Fig. 10. Exact and Normal approximation of equilibrium size distributions for multiple birth badger models with two immigration levels.

The moments for these data are m1 ¼ 4:45; m2 ¼ 22:55 and m3 ¼ 126:982; for which the variance is s2 ¼ 2:7475 and skewness is k3 ¼ 2:18: The coefficient of variation is 0.37, which is very close to that of badger litters, and the index of skewness is 0.49, which is somewhat smaller (28%) than for the badgers.

0.05 0.04 0.03 0.02

Probability

0.01

20

40

60

80

(A)

0.05 0.04

The previous population size models with unit births had only one rate parameter, a1 ¼ 1:00; which was reduced by m1 ¼ 4:45 to a01 ¼ 0:2247: This was substituted into (31), so that the models with multiple births would have mean value functions comparable to the single birth models. (b) Model solutions: The exact solutions for the three low-order cumulant functions are given in Fig. 12, the cumulants at t ¼ 100 in Table 3B, and the error rates in Table 4B. The following results stand out:

0.03 0.02 0.01

20 (B)

40 Size

60

80

Fig. 11. Exact and approximate size distribution at time of maximum relative skewness for multiple birth badger model with low immigration, I ¼ 2: (A) Normal approximation and (B) saddlepoint approximation.

qualitative stage in the reproductive cycle than the badgers, namely during the last 10 days of gestation. This distribution represents our best available data and is used for illustrative purposes. The distribution is (Lloyd, 1980, p. 280): i 1 2 3 4 5 6 7 8 9 10 pðiÞ 0.035 0.070 0.176 0.211 0.317 0.087 0.0070 0.017 0 0.017

1. Both cases reach an effective equilibrium by t ¼ 100; with equilibrium means of k1 ¼ 96:3 and 94.7 for high and low immigration, respectively. Both of these are substantially different from the deterministic carrying capacity of K ¼ 100: The corresponding variances are k2 ¼ 441:8 and 552.4 with skewnesses of k3 ¼ 2247 and 2453. The indices of skewness, g ¼ 0:24 and 0.19, are substantial. 2. The cumulant function approximations for the two cases have substantial errors for a short period of time in their variance and skewness functions around their peak values. The variance for the low immigration case has an error rate of 21% at t ¼ 5:0; the time of peak variance, but it reduces rapidly and is less than 1% by t ¼ 17: Similarly, the maximum error in the skewness approximation is very large, roughly 60%, in this case and is sizeable also for the high immigration case. However the errors diminish rapidly after the times of peak values, which in both cases is t ¼ 3:4: The error rates for the equilibrium values, which are small, are representative of most of the time range. Therefore, the exact solutions in

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104 100

I=10 I=3 Variance

Mean

80 60 40 20

1000

10000

800

8000 Skewness

100

600

I=3

400

I=10

4000 I=10 2000

200

0

10

20 30 Time (yrs)

(A)

40

0 0

50

I=3

6000

10

(B)

20

30 Time (yrs)

40

50

0

10

20

(C)

30 Time (yrs)

40

50

Fig. 12. Exact cumulant functions for family of multiple birth fox models with two immigration levels: (A) mean value, (B) variance, and (C) skewness.

0.025

0.025

0.025

Probability

I=3

I=10 (Saddle)

I=10

0.02

0.02

0.02

0.015

0.015

0.015

0.01

0.01

0.01

0.005

0.005

0.005

40 60 80 100 120 140 160 Size

40 60 80 100 120 140 160 Size

40 60 80 100 120 140 160 Size

Fig. 13. Exact and Normal approximation of equilibrium size distributions for multiple birth fox models with two immigration levels, I ¼ 3 and 10, and also saddlepoint approximation for I ¼ 10:

Fig. 12 for the models with immigration are indistinguishable from the approximate solutions except for short periods of time around the peak values. The errors for the equilibrium mean and variances approximations are all small, however the skewness approximations have a roughly 5% error for both cases.

0.025 0.02 0.015 0.01 0.005 Probability

Fig. 13 illustrates the exact equilibrium distributions with their Normal approximations. In both cases, there is small but visibly noticeable lack of symmetry. As expected, the saddlepoint approximations are more accurate in each case. As an illustration, Fig. 13 also contains the very accurate saddlepoint approximation for the high immigration case, which at g ¼ 0:24 has the largest index of skewness. The maximum magnitudes for the indices of skewness over time are g ¼ 0:43 for high immigration and g ¼ 0:72 for low immigration. The Normal and saddlepoint approximations for the low immigration distribution at t ¼ 3:4; the time of maximum skewness, are illustrated in Fig. 14. The saddlepoint approximation is a clear improvement over its admissible range. These approximation results seem remarkable given that this is a worst-case scenario.

0.03

The findings for the multiple birth models are similar qualitatively to those for the single birth models. All

40

60

80

100

120

140

20

40

60

80 Size

100

120

140

0.03 0.025 0.02 0.015 0.01 0.005

(B)

4.4. Summary for multiple birth models

20 (A)

Fig. 14. Exact and approximate size distribution at time of maximum relative skewness for multiple birth fox models with low immigration, I ¼ 3: (A) Normal approximation and (B) saddlepoint approximation.

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

three cumulant functions are approximated well, except sometimes for a short period of time in the neighborhood of their peak values. The mean and variance approximations for the equilibrium distributions are very accurate. The approximation for the equilibrium skewness has substantial proportional errors in some cases, however the skewness in equilibrium is relatively small. The equilibrium distributions all have small though visible lack of symmetry, hence the Normal approximations with the cumulant approximations seem adequate for most practical purposes. Nevertheless, the saddlepoint approximation provides a visible improvement in accuracy in every case, and it provides a good general approximation even at the transient time points of large relative skewness.

5. Effect of multiple births on equilibrium distributions Two types of comparisons on the cumulants of the equilibrium distributions are of interest. The first is a straightforward comparison of cumulants for the simple- vs. multiple birth mechanisms. The other is a comparison of the effects of the two types of nonlinear mechanisms in the paper, one with density-dependent births and the other with densitydependent deaths. The two ‘‘real-world’’ examples used in this paper are based on two litter size distributions which fortuitously are quite similar. The two coefficients of variation are very close, with cv ¼ 0:36 and 0.37, and the two indices of skewness are roughly comparable, at g ¼ 0:68 and 0.49. Though not exactly the same, the similarity in litter size distributions facilitates the qualitative comparisons of the two density-dependent mechanisms. The subsequent comparisons are based for simplicity on the exact cumulants of the equilibrium distributions given in Tables 1 and 3. However the cumulant approximations are very close in each case and therefore have the same qualitative comparisons.

101

equilibrium is (see, e.g., Matis and Kiffe, 1996): k
3 ¼ ðb1  b2 Þk
2 =ðb1 þ b2 Þ:

ð34Þ

The population control for the foxes comes through density-dependent mortality and the assumed model has a linearly increasing birth rate, with b1 ¼ 0 in (10). Hence, without immigration, k
3 ¼ k
2 with resulting positive skewness. The positive skewness carries over in the present fox models even with immigration. Conceptually, the larger the population size, the more frequent are the birth events which would tend to produce an upper tail in the population size distribution, yielding positive skewness. Multiple births would accentuate outlying large population sizes, and increase the skewness as indicated in the observed results. On the other hand, population control for the badgers comes through density-dependent births, with b2 ¼ 0 in (10). Hence, without immigration, k
3 ¼ k
2 in (34), with resulting negative skewness. The negative skewness is also manifested in our present badger models with immigration. The assumed quadratic birthrate generates a population size, u; above which the birth rate is zero. Conceptually, though the population size, N; is unbounded due to immigration, the quadratic birthrate function slows down the birth rate drastically as N approaches u; which in turn dampens the upper tail. This constraining force on the upper population sizes tends to yield negative skewness. For sure, the skewness effects are very tiny with the present examples. However the observation that a quadratic birthrate tends to yield negative skewness on the population size distribution is supported by all our previous experience. In particular, the indices of skewness for our models with immigration and unit density-dependent births in Matis and Kiffe (2000) are g ¼ 0:62 and 1:06 for our two African bee models and range from 0:154 to 0.034 for our six muskrat models, the single positive skewness occurring for the case of massive immigration in which the immigration rate exceeded the birth rate for any population size.

5.1. Comparisons of skewnesses

5.2. Comparisons of variances

Two results are apparent. First, the magnitude of the skewness for the multiple birth case always exceeds that of the corresponding single birth case. For the badgers, the increases associated with multiple births are 49% and 100% for the two immigration cases, whereas for the foxes both cases have a substantial, roughly 12-fold increase. Second, it is noteworthy that the equilibrium skewness values for all badger models are negative, whereas for all fox models, they are positive. The result is expected in the absence of immigration, as the Bartlett et al. (1960) approximation for the skewness in

It seems self-evident that as one introduces an additional source of variation to the model through the random litter size, the variability of the population size in the resulting process would increase, as with the linear processes in Matis et al. (1995). This property is apparent in the results in Tables 1 and 3. The incorporation of multiple births increases the variances of the equilibrium distributions by 62% and 80%, respectively, for the badger models and by 181% and 198%, respectively, for the fox models. These differences are illustrated in Fig. 15, which compares the Normal approximations of the equilibrium distributions.

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

102

0.1

0.1

Probability

Badgers, I=2

Foxes, I=3

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02 20

40

60

80

100

0.1

20 40 60 80 100 120 140 0.1

Badgers, I=6

Foxes, I=10

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02 20

40

60

80

100

20 40 60 80 100 120 140 Size

Fig. 15. Comparative Normal approximations of equilibrium size distributions for single (—) vs. multiple (- -) birth models for two immigration cases of badger and fox population growth models.

The relationship of the increase in variance to the type of density-dependent mechanism may be investigated by obtaining the difference of k’ 2 ðtÞ in (22) and (32). Let m1 ¼ 1 in (32), which adjusts only for scale differences in the a1 and b1 coefficients. The difference in these rates is then Dk’ 2 ðtÞ ¼ ðm2  1Þ½a1 k1  b1 k21  b1 k2 :

ð35Þ

This simple difference ignores the fact that k1 and k3 are not equal in the two models, due to the effect of the higher-order litter size moments, i.e., for mi with i41: Therefore, Dk’ 2 ðtÞ is only an approximation of the true difference. Note that the coefficients in the second factor on the right are the birth rate parameters, a1 and b1 ; and in fact this factor is equal to k’ 1 ðtÞ in (22) with I ¼ a2 ¼ b2 ¼ 0: Hence (35) simplifies to Dk’ 2 ðtÞ ¼ ðm2  1Þk’ 1 ðtÞ;

ð36Þ

in which k’ 1 ðtÞ represents the mean for a model without immigration nor death. This population would be monotonically increasing, i.e., k’ 1 ðtÞ40: Moreover m2 41 for multiple births (except in the trivial case j ¼ 0; 2 with pð jÞ ¼ 0:5). Hence Dk’ 2 ðtÞ40 for all t in (36) implying that the variance of population size is larger for the multiple birth model. In our two examples, with m2 ¼ 6:45 for badgers and m2 ¼ 22:55 for foxes, the increases average 70% for the former and 190% for the latter. It is also clear from (35) that the size of the difference is related to the type of density-dependent mechanism. Because k2 is positive, population control through density-dependent births, with b1 40 as with badgers,

would tend to have a smaller increase in equilibrium variance than the control through density-dependent deaths, with b1 ¼ 0 as with foxes. This is consistent with the observed results, e.g., in Fig. 15, in which the fox models have large proportional increases in variance due to multiple births than the badger models. 5.3. Comparisons of means Setting m1 ¼ 1 in (32) to adjust for scale differences, k’ 1 ðtÞ for both the single and multiple birth models, from (22) and (32), is k’ 1 ðtÞ ¼ I þ ½a  bk1 ðtÞ k1 ðtÞ  bk2 ðtÞ; where b40: Because as noted previously, the variance k2 ðtÞ for the multiple birth model exceeds that of the corresponding single birth model, it follows that the equilibrium mean of the former will be less than that of the latter, as is the case in all present models. It also follows that, for comparable immigration, the model with density-dependent births will have a smaller decrease in the mean, due to its smaller variance. Consistent with this observation, the reductions in equilibrium sizes are 0.8% and 1.7% for the badger models, compared to 2.4–3.6% for the comparable fox models.

6. Conclusions Multiple births are common for numerous animal species. This paper presents the differential equations

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

for the mean, variance, and skewness functions of the transient population size distributions of stochastic logistic growth models with multiple births. These equations generalize equations previously given assuming single births during birth episodes. Previous research has shown that one can use moment closure methods to find accurate approximate solutions for the single birth stochastic logistic model, whether without immigration (Nasell 2001, 2003) or with immigration (Matis and Kiffe, 1999). Outstanding among the findings of this research is that one may also use such methods to obtain very accurate approximations of the equilibrium mean and variance of the stochastic logistic models with immigration and realistic multiple birth distributions. The skewness approximation has larger proportional error, but the magnitude of the error is small. Consequently, the equilibrium size distribution is quite accurately approximated by the Normal, and even more accurately by the saddlepoint approximation, for the cases investigated. This research was motivated by the results in Matis et al. (2003) in which gives accurate analytical approximations of the equilibrium size distributions for the simple logistic model without immigration. The current methods for finding equilibrium distributions for the model with immigration are based directly on the balance equations in (13) or on the stationary solution of the Kolmogorov equations. Both of these methods are formidable numerically for ecological researchers, and especially so for the generalization to the multiple birth case. The method in this paper for finding the equilibrium distribution approximation is far easier, consisting of finding the stationary solutions for the three equations in (22) or (32), and then using the Normal or saddlepoint, in (24), approximation. The transient distribution approximations are found by solving the set of three differential equations in (22) or (32), together with an equation for the fourth-order cumulant. We are hopeful that the simplicity of the new methodology will lead to a much wider utilization in practice of the biologically more realistic multiple birth model. An important reason for using the multiple birth model, when appropriate, instead of the simpler single birth model is that the inclusion of this additional source of stochasticity to the model increases the variance of the resulting population size variable at any given time, including of course the equilibrium distribution, as in Fig. 15. The underestimation of the natural variability of population size resulting from improper use of the single birth model may have serious consequences when the resulting distributions are used in animal management models. We recognize that more research is required to investigate the characteristics of the models and their properties, such as accuracy, of the approximation

103

methods for litter size distributions that are very different, for example those encountered in many nonmammalian species. Research is in progress to develop user-friendly procedures for generalizing the saddlepoint approximation to include the full range of possible population sizes. The methods to do so are outlined in Wang (1992), and we seek simple ways to implement them in practice. Research is also required to incorporate the multiple birth mechanism into other models. One such area is into stochastic models without immigration, as in Nasell (2001, 2003), which is outside the scope of our present application areas. Research is in progress to extend the methods to include such multiple birth mechanisms in multi-species models, such as in the host–parasite model for bees and mites in Matis and Kiffe (2002). We contend that stochasticity from numerous factors is inherent in all natural systems, and that practitioners will utilize such enriched stochastic models when userfriendly methodology becomes available.

Acknowledgments The suggestions of the associate editor (K.D.) and two anonymous referees have improved the paper substantially, and we are grateful to them for their input.

Appendix Eq. (18) could be derived by direct but tedious manipulation of the Kolmogorov equations. However it is easily obtained using the ‘‘random-variable’’ (operator) technique. Let the probabilities of possible changes in the NðtÞ count from t to t þ Dt be denoted as ProbfNðtÞ changes by j unitsg ¼ fj ðNÞDt þ oðDtÞ: The possible changes assumed in (11), (12) and (16) may then be written in terms of the fj functions as f1 ¼ a1 N  b1 N 2 þ I;

ð37Þ

f1 ¼ a2 N þ b2 N 2 :

ð38Þ

Bailey (1964, p. 73) gives the following operator equation for the moment generating function, Mðy; tÞ:   @M X jy @ ¼ ðe  1Þ fj Mðy; tÞ; ð39Þ @t @y ja0 where for X f ðNÞ ¼ ak N k one has   X @k M @ f : M¼ ak @y @yk

ARTICLE IN PRESS J.H. Matis, T.R. Kiffe / Theoretical Population Biology 65 (2004) 89–104

104

Eq. (18) follows by substituting (37) and (38) into (39). The technique is also illustrated in Matis and Kiffe (2000). Similarly, the changes assumed for the multiple birth models are (38) and f1 ¼ ða1 N  b1 N 2 Þpð1Þ þ I; 2

fi ¼ ða1 N  b1 N ÞpðiÞ

and

for i41:

Substitution of these into (39) yields (30). References Anderson, R.M., Trewhella, W., 1985. Population dynamics of the badger (Meles meles) and the epidemiology of bovine tuberculosis (Mycobacterium bovis). Philos. Trans. R. Soc. London B 310, 327–381. Anderson, R.M., Jackson, H.C., May, R.M., Smith, A.M., 1981. Population dynamics of fox rabies in Europe. Nature 289, 765–771. Bailey, N.T.J., 1964. The Elements of Stochastic Processes. Wiley, New York. Banks, R.B., 1994. Growth and Diffusion Phenomena. Springer, Berlin. Bartlett, M.S., Gower, J.C., Leslie, P.H., 1960. A comparison of theoretical and empirical results for some stochastic population models. Biometrika 47, 1–11. Coyne, M.J., Smith, G., McAllister, F.E., 1989. Mathematic model for the population biology of rabies in raccoons in the mid-Atlantic states. Am. J. Vet. Res. 50, 2148–2153. Johnson, N.L., Kotz, S., Kemp, A.W., 1992. Univariate Discrete Distributions, 2nd Edition. Wiley, New York. Lloyd, H.G., 1980. The Red Fox. B. T. Batsford, London. Matis, J.H., Kiffe, T.R., 1996. On approximations the moments of the equilibrium distribution of a stochastic logistic model. Biometrics 52, 980–991. Matis, J.H., Kiffe, T.R., 1999. Effects of immigration on some stochastic logistic model: a cumulant truncation analysis. Theoret. Popul. Biol. 56, 139–161.

Matis, J.H., Kiffe, T.R., 2000. Stochastic Population Models—A Compartmental Perspective. Lecture Notes in Statistics, Vol. 145. Springer, New York. Matis, J.H., Kiffe, T.R., 2002. On interacting bee mite populations: a stochastic model with analysis using cumulant truncation. Environ. Ecol. Statist. 9, 237–258. Matis, J.H., Zheng, Q., Kiffe, T.R., 1995. Describing the spread of biological populations using stochastic compartmental models with births. Math. Biosci. 126, 215–247. Matis, J.H., Kiffe, T.R., Renshaw, E., Hassan, J., 2003. A simple saddlepoint approximation for the equilibrium distribution of the stochastic logistic model of population growth. Ecol. Modelling 161, 239–248. Nasell, I., 2001. Extinction and quasi-stationarity in the Verhulst logistic model. J. Theoret. Biol. 211, 11–27. Nasell, I., 2003. Moment closure and the stochastic logistic model. Theoret. Popul. Biol. 63, 159–168. Neal, E.G., 1977. Badgers. Blandford, Poole, Dorset, UK. Pech, R.P., Hone, J., 1988. A model of the dynamics and control of an outbreak of foot and mouth disease in feral pigs in Australia. J. Appl. Ecol. 25, 63–77. Pielou, E.C., 1977. Mathematical Ecology. Wiley, New York. Renshaw, E., 1991. Modelling Biological Populations in Space and Time. Cambridge University Press, Cambridge, UK. Renshaw, E., 1998. Saddlepoint approximations for stochastic processes with truncated cumulant generating functions. J. Math. Appl. Med. Biol. 15, 1–12. Roberts, M.G., 1996. The dynamics of bovine tuberculosis in possum populations, and its eradication or control by culling or vaccination. J. Anim. Ecol. 65, 451–464. Smith, P.J., 1995. A recursive formulation of the old problem of obtaining moments from cumulants and vice versa. Am. Statist. 49, 217–218. Wang, S., 1992. General saddlepoint approximations in the bootstrap. Statist. Probab. Lett. 13, 61–66. Whittle, P., 1957. On the use of the Normal approximation in the treatment of stochastic processes. J. Roy. Statist. Soc. B 19, 268–281. Wolfram, S., 1991. Mathematica. A System for Doing Mathematics by Computer, 2nd Edition. Addison-Wesley, Redwood City, CA.