Casting the Population Model as a Statistical Model

Statistical models comprise two interwoven components: a deterministic model describing the deterministic biological process and a stochastic component describing how random variation influences the process. The deterministic component typically, but not always, forms the biological process of interest. For example, the deterministic model may be relatively straightforward, such as exponential population growth; could include a description of multiple populations, such as coupled predator-prey models; or may explicitly describe within-population structure, such as size-structured dynamics. The possibilities are endless - and really only limited by the numerical methods available to simulate population dynamics. The stochastic component(s) of a statistical model are described by probability distributions that can take on wide range of functions. For example, probability distributions can be used to represent observation error, to describe demographic stochasticity, to account for environmental stochasticity in model parameters, or to account for variation among individuals in the population from structure not accounted for in the deterministic model (e.g., size structure). The stochastic component of a statistical model is often thought of as an 'error' term that allows for variation around a deterministic model. However, while true for pure observation error, stochasticity can introduce randomness into the parameters and variables that generates population dynamics not present in deterministic models. As a result, the statistical model that describes the biological process necessarily includes a careful statement of both the stochastic and deterministic components.

As a practical example, consider a series of statistical models to describe the algal dynamics illustrated in Figure 1. Based on our understanding of algal growth in a limited medium, we may feel that the population dynamics can be described by continuous logistic growth:

model depends on how the process error is introduced, an example is

where Nt is the predicted algal density at time t, r is the maximum per capita growth rate, and K is the algal carrying capacity. While the deterministic model for logistic growth has an analytical solution, many of the statistical models based on logistic growth must be solved numerically. To reflect this, and the general fact that biologically motivated statistical models rarely have analytical solutions, the discussion below focuses on numerical methods.

A number of statistical models can be derived from eqn [ 1 ] depending on how stochasticity is introduced. If stochasticity influences how the states are observed, but not their inherent dynamics, it is referred to as 'observation error'. For our example of algal growth, the statistical model with only observation error has the following structure:

where Yt is the observed algal density, and f (N, 8) is a probability distribution with parameter(s) 8 that describe the observation error surrounding the true density N? Ifstochasticity directly influences the population dynamics, then it is referred to as 'process error'. While the final dNt f Nt —-— = rNA 1--

where Wt is the random variation through time, and g(Nfi 0) is a probability distribution with parameter(s) 0 (for short, written as Wt ~ g(Nt, 0)). The stochastic component Wt is typically expanded as Wt = u(Nt)£f, where u(Nt) is the intensity of the stochasticity and is a probability distribution with mean zero and unit variance describing the process noise (abbreviated as Wt in the equations below). The statistical models shown in eqns [2] and [3] are just two examples; many others are possible depending on how stochasticity is thought to influence the observed dynamics. For example, a model with both observation and process error would be dNt ( NÄ - = rN^1-T] +

where Wt is the process error and Yt is the observed algal density.

The final step to creating a statistical model is to define the probability distribution(s) that describe the stochastic process. As with the deterministic model, there are many choices available (e.g., Table 1). If the response variable (i.e., the event that the probability distribution describes) is continuous and takes on values between negative infinity and positive infinity (referred to as the domain), then a normal distribution may be an appropriate choice. If the continuous response variable can take on only positive values, then a lognormal or gamma distribution may be a good choice. When the response variable is discrete, such as count data, then the probability distribution could be described with a Poisson or negative binomial distribution. The choice of probability distribution should reflect whether the response variable is continuous or discrete, the domain of the stochastic process, and how well the shape of the distribution reflects the stochastic process (e.g., relationship between the mean and the variance). In some situations, probability distributions emerge directly from the assumptions of stochastic process. For example, if we assume that the change in population numbers results from a sufficiently large number of independent births and deaths, then the stochastic component of eqn [3] is a Wiener increment written as Wt = , where is a normal distribution with mean zero and unit variance (£f ~ N(0, 1 )) and u is the process error variance.

Table 1 Common probability distributions

Probability distributiona

J2ira2

E{Y} = ex^ ß + — Var{Y} = exp(2ß + 2a2) - exp(2ß + a2) Gamma f (Y = y|a,ß) =

Y is a continuous random variable that can take on values between negative and positive infinity. The parameter p is the mean, and a2 is the variance.

Y is a strictly positive continuous random variable.

Y is a strictly positive continuous random variable. The parameter a is the shape parameter, ¡3 the scale parameter, and r(a) is the gamma function.

Y is continuous random variable strictly between 0 and 1. The function r() is the gamma function.

Y is a discrete random variable, where y is the number of successes in n trials, and p is the probability of success for each trial.

Y is the discrete number of occurrences over a fixed period of time, and A is the intensity parameter.

Y is the discrete number of Bernoulli trials required to get the rth success. p is the probability of success for each trial.

Y is a discrete random variable, and r() is the gamma function. The beta-binomial is referred to as a hierarchical model because it is a binomial distribution where the rate parameterp is itself a beta distribution. It is often used when the variance in the observed data is greater than binomial.

(Continued)

Table 1 (Continued)

Probability distributiona

Multivariate distributions Multinomial f (Y = yi,..., Yk = yk\n, pi,..., Pk -1) n! fk -1 \f k -1

Y is the discrete number of successes of the /'th outcome (of k) in n trials, and Pi is the probability of success for each trial. Similar to the binomial distribution, overdispersion in the multinomial can be included by assuming that the parameters p have a Dirichlet distribution.

aE{Y} is the expectation (average value) of the distribution, and Var{Y} is the variance.

The equations are shown using statistical notation as follows: uppercase letters denote a random variable (Y), and lowercase letters denote a realized value of the random variable (y); the function f() denotes either the probability mass function (for discrete random variables) or the probability density function (for continuous random variables); the symbol j' represents the phrase 'given', and the lowercase letters denote the parameters of the distribution. Thus, the notation f(Y=y I ^,a2) is read as 'the probability mass (or density) that the random variable Ywill take on a particular realization y, for the parameters ^ and a2