A parameter estimation routine for the vitality-based survival model
Introduction
Anderson, 1992, Anderson, 2000 presents the vitality-based survival model; a parametric model for relating stressors and environmental properties to organism survivorship in terms of three parameters: the mean rate of vitality loss, the variability in that rate and a rate of accidental mortality. In that work, parameter estimates were obtained via a maximum likelihood estimator (MLE) using the ms( )—minimum sum—routine in S-Plus® (see Chambers and Hastie, 1992). However, numerical overflow errors and over-sensitivity to initial values caused the routine to be unstable. In this paper, we present a robust MLE routine, still based on the ms( ) routine in S-Plus®, for this three-parameter decaying inverse Gaussian model. The routine contains an initial parameter estimation scheme and numerical safeguards as well as a treatment of right censored data. In addition, the issue of when to define the start time is discussed. For access to the parameter estimation code, see the Appendix A.
In an editorial, Jorgensen (1997) states that “parameter estimation is still one of the weak points in ecological modelling ….” And it seems that each modelling situation brings its own set of parameter estimation difficulties. The estimation routine can be highly sensitive to initial or intermediate parameter values. Givens and Poole (2002) show how the population dynamics models used by the International whaling commission for assessments of bowhead whale stocks leads to troublesome likelihood surfaces which, with the data available, leads to a nearly chaotic dynamical system. Multiple routines may be necessary to first estimate and then refine parameters. Stollenwerk et al. (2001) simplify the nonlinear parameter estimation for a stochastic phytoplankton biomass time series model by performing the estimation in several steps; first linearizing around a deterministic fixed point to get an initial MLE estimate of some parameters, then by simulating the state-dependent noise to estimate the remaining autocorrelation parameter, and finally by adjusting the autocorrelation parameter by numerically minimizing the deviation of the cumulative spectrum from the empirical one. Estimation routines must be be as flexible as the models they are built for. Sielken (1985) chronicles the development of the MLE routine for the Hartley–Sielken family of dose–response models; one that can handle a variety of data situations including right censored data. The Hartley–Sielken model, like the vitality model, is derived for either constant continuous dosing or a single exposure at one point in time. Our parameter estimation routine was constructed with difficulties like those mentioned above, in mind. There are, of course, many other numerical issues which can be encountered.
The vitality model (Anderson, 2000) is based on the premise that mortality occurs through two processes; one dependent on the organism’s history and the other independent of that history. The history-dependent component is defined by a directed random walk (Weiner process) of the abstract variable, vitality, to its absorbing boundary (death). This age-dependent process requires two parameters—r, the a rate of vitality loss and s, a variability in that rate—and can be formulated as one minus an inverse Gaussian survival distribution. The history-independent, or accidental, component of the mortality acts independently of age or condition of the organism and is described by an exponential model (Poisson process) characterized by the average time between mortality events 1/k.
Combining these mortality components gives the cumulative survivalwhere Φ is the cumulative normal distribution.
The model characterizes the effect of natural and xenobiotic stressors on organism survival. For many stressors the rate of vitality loss, r, can be linearly related to stressor level. By fitting the model to survival curves produced under different levels of a stressor, empirical relationships between the model parameters and stressor level can be established. The parameters not only characterize the survival curve, but quantify the temporally changing distribution of vitality within a population and provide information on the selective mortality of weaker population members over the span of the survival observations. (See Anderson, 2000 for further explanation.)
Section snippets
Maximum likelihood estimator
Typically, survival study data is interval-censored; that is the mortalities are counted at the end of each time period rather than continuously. To define a maximum likelihood estimator (MLE) of the model parameters using interval censored survival data, the incremental mortality probability (density) is writtenfor the T time intervals defined by [ti−1,ti] for i∈[1,T]. The likelihood function is thenwhere ni is the observed mortality in time period i. The maximum
Numerical safeguarding
As with any numerical routine, care must be taken to safeguard the routine from numerical problems. The three most important cases in the parameter search algorithm, with their remedies, are discussed below.
In the original routine, the main numerical difficulty arose from the exponential e2r/s2 in the survival function Eq. (1), which in the parameter search could result in an overflow error if 2r/s2 became too large. To remedy this as well as to reduce the numerical error when the extremely
Goodness-of-fit test
A goodness-of-fit test is included as an indicator of whether the data could reasonably have originated from the probability distribution Eq. (2) with parameters as estimated by the method of maximum likelihood. For this case (where the distribution parameters are estimated from the data) we use an analog of the χ2 (or Pearson’s C) test (cf. Larsen and Marx, 1986, Section 9.4)where the s are computed by replacing the unknown parameters in Eq. (2) with their MLEs. Here, n
Examples
To examine the ability of the routine to fit both poorly and well-conditioned data sets, we present two examples. One from a study of the mortality effects of various levels of gas supersaturation on rainbow and brown trout (White et al., 1991). In the example used here, 75 rainbow trout of mean weight 3.4 g were subjected to gas supersaturation at a level of 125%. This was chosen as an example of right censored data and to show the ability of the routine to deal with this particularly scattered
Acknowledgements
This work was supported by the Bonneville Power Administration under contract number 00000154-00001, project number 89-108.
References (12)
- et al.
Testing nonlinear stochastic models on phytoplankton biomass time series
Ecol. Model.
(2001) - Anderson, J.J., 1992. A vitality-based stochastic model for organism survival. In: DeAngelis, D.L., Gross, L.J. (Eds.),...
A vitality-based model relating stressors and environmental properties to organism survival
Ecol. Monographs
(2000)- Chambers, J.M., Hastie, T.J. (Eds.), 1992. Statistical Models in S. Wadsworth & Brooks/Cole, Pacific Grove, California,...
Editorial: ecological modelling in 100 volumes
Ecol. Model.
(1997)- Kendall, M., Stuart, A., 1997. The Advanced Theory of Statistics, 4th ed., vol. 2. MacMillan, New York,...
Cited by (16)
The vitality model: A way to understand population survival and demographic heterogeneity
2009, Theoretical Population BiologyJoint modelling of breeding and survival in the kittiwake using frailty models
2005, Ecological ModellingCitation Excerpt :Such models can prove difficult to fit using classical approaches (Link et al., 2002), but very few attempts have been made. Parameter estimation has long been recognized as central to ecological modelling [(Jorgensen, 1997; Williams et al., 2002; Salinger et al., 2003)]. Our main objective is to develop a new parametrization for a model to estimate the survival and breeding probability jointly using data from a long-lived species (the kittiwake).
Contaminants as viral cofactors: Assessing indirect population effects
2005, Aquatic ToxicologyIntrinsic units in growth modeling
2004, Ecological Modelling