## QLy nyf11 n

Solving the above equation for p, and checking that the second derivative is negative, yields the maximum likelihood estimator p = y/ n (the specific value of p-hat for a given value of y and n is the maximum likelihood estimate). In general, however, it is rare to find the maximum likelihood estimators analytically. For these situations, the maximum n Survival probability ( p)

iv at

Survival probability ( p)

Figure 2 Likelihood functions for a binomial distribution under three sets of observed data: (a) shows the likelihood profiles, and (b) shows the negative log-likelihood profiles. Observed data are (X, N) - (2, 3) (solid line), (X, N) - (10, 15) (dashed line), and (X, N) = (20, 30) (dotted line). For ease of comparison, the likelihood and negative log-likelihood profiles are scaled by the maximum or minimum value, respectively. The maximum likelihood estimate of p, or equivalently the minimum negative log-likelihood, occurs at p — X/N (bold vertical line).

likelihood estimates are found by numerically maximizing the likelihood function (Figure 2). Equivalently, it is common practice to find the maximum likelihood estimates by minimizing the negative log-likelihood function:

The above example of cohort survival is relatively straightforward; it illustrates the general approach of using maximum likelihood to fit biologically motivated statistical models to population dynamics. Two of the biggest challenges encountered when fitting more complicated statistical models are constructing the likelihood function, and numerically finding the most likely parameter estimates. To discuss these details for the different types of statistical models in eqns -, it is useful to return to the algal dynamics example.

Models with observation error only

Equation  shows a statistical model with only observation error. To start, assume that the observation error is normally distributed:

where Nt is the true population dynamics at time t and —2 is the variance. A normally distributed observation error means that observations can take on negative values, which may be appropriate for our example if algal density was estimated indirectly through an index, such as relative flourometric units. The corresponding likelihood and negative log-likelihood function for the entire data set are

L(r, K, u21 y) - f (yi | r, K, u2) ■■■ f {y„\r, K,u2)

where Nt (r, K) is the deterministic prediction given by the solution to eqn  that depends on parameters r, K, and the initial algal abundance N0.

The likelihood function in eqn  depends on three variables, which means that the function requires minimization over three dimensions. To help visualize the likelihood function, Figure 3 shows the negative log-likelihood contours for different values of r and K, but for a fixed value of — . The minimum value of the negative log-likelihood function, and thus the most likely set of parameter estimates, are obtained by numerically searching over parameter space.

A number of search algorithms can be used to minimize the negative log-likelihood function (Table 2). Each takes a different approach, but all search algorithms attempt to iteratively move through parameter space toward lower values of the function, while striking a balance between robustness (in the sense of moving in the correct direction) and avoiding excessive function evaluations. The problem is that most search algorithms are not guaranteed to find the global minimum. This can be seen in the search paths of the algal dynamics example (Figure 3). Depending on where the search algorithm starts in parameter space, the search paths either arrive at the global minimum or end prematurely at the incorrect parameter estimates. While o o o o o

CS O

 1 \ — • Maximum per capita growth rate (r) Maximum per capita growth rate (r) Figure 3 Negative log-likelihood surface from eqn  for the algal population dynamics in Figure 1. The surface is shown for different parameters values of r and K, and a fixed value of u2. Three example search paths are shown (dashed lines). Each shade of gray represents a different starting point in parameter space. The start point is shown with an open circle. One of the search paths finds the minimal value, but two terminate prematurely at incorrect parameter estimates. search algorithms are programmed as automated routines, the process of finding the most likely parameter estimates requires user interaction to ensure that the global minimum has been found. Some good practices to keep in mind are to start searches from multiple locations in parameter space, experiment with multiple search routines, and plot the likelihood profile in the area surrounding the solution. Evolutionary algorithms (Table 2), and a new methodology called data cloning (see below for more details) offer the possibility of global minimization. The most likely parameter estimates for the statistical model given by eqn  are r = 1 .70, K = 9.90 x 1 05, and u2 = 1 .82 x 1 09. So far, we have assumed a normal distribution for observation error. To develop statistical models for other distributions, the parameters of the distribution describing observation error need to be linked to the deterministic model. To illustrate the process, Table 3 summarizes the steps for fitting statistical models with lognormal distributed observation error, gamma distributed observation error, and Poisson distributed observation error. The most likely parameter estimates for the different distributions are shown in Table 4, and the predicted dynamics for the most likely parameter estimates are superimposed over the raw data in Figure 1. 