An adaptive feedback methodology for determining information content in stable population studies.

We develop statistical and mathematical based methodologies for determining (as the experiment progresses) the amount of information required to complete the estimation of stable population parameters with pre-specified levels of confidence. We do this in the context of life table models and data for growth/death for three species of Daphniids as investigated by J. Stark and J. Banks [17]. The ideas developed here also have wide application in the health and social sciences where experimental data are often expensive as well as difficult to obtain.


1.
Introduction. Pesticides protect crops and improve their productivity by killing or repelling pest species. However, they also present a risk to nontarget species and ecosystems. For example, processes such as spray-drift, leaching, and runoff can contaminate aquatic ecosystems [14]. Thus an ecological risk assessment (ERA) is often used to evaluate the effects of a chemical on the environment. These ERAs are mainly based on measurements of individuals, such as acute mortality, even though the goals for protecting the environment are often in terms of a population [13,15]. Thus some ecotoxicologists are proposing and presenting supporting research for the use of demographic data with the population growth rate as the endpoint of interest [10,11,12,13]. However, the collection of demographic data can be both time-consuming and costly. Therefore an important question is whether partial demographic data can be used instead of full demographic data while still providing an accurate picture of the impact of a toxicant on a population. This question was recently investigated by Stark and Banks [17] using life table data for three species of Daphniids, Ceriodaphnia dubia, Daphnia magna, and D. pulex. Life tables were developed weekly, and several population endpoints, including the stable population growth rate [9] (λ), were compared. The authors statistically compared the population growth rate for each week to the population growth rate calculated from the complete life table to determine the earliest time that gave statistically insignificant additional results.
We approach this problem from a modeling perspective in order to develop an adaptive feedback mathematical methodology to determine the minimum amount of data needed to estimate stable model parameters with pre-specified confidence levels. We first present the parameter estimation or inverse problem methodology that we use in our analysis in Section 2. In Section 3, we investigate the problem and explain our ideas using the classic and well-known Verhulst-Pearl logistic growth model. Using our approach developed and illustrated with this simple example, we then in Section 4 consider the data used by Stark and Banks [17]. Finally, in Section 5 we present our conclusions and directions for further research.
2. Inverse problem methodology. We first describe the mathematical and statistical model used in the inverse problem (See [8,7] for further details). Let x represent the N -dimensional vector of state variables and let θ θ θ = [q, x 0 ] (assuming initial conditions x 0 are unknown) represent the parameters to be estimated in the mathematical model dx dt = g(t, x(t), q) (1) x(t 0 ) = x 0 .
We now consider n longitudinal data {y j } n j=1 , collected at discrete time points {t j } n j=1 . Since it is possible for a state variable to not be observed in the collected data, we define the m-dimensional observation process f (t; θ θ θ) = Cx(t; θ θ θ), with m ≤ N . Then f (t j ; θ θ θ) = Cx(t j ; θ θ θ), for j = 1, . . . , n. It should be noted that the data {y j } will not be exactly f (t j ; θ θ θ) since measurement procedures contain error or uncertainty. Thus it is necessary to discuss the statistical model as well as the mathematical model.
The vector θ θ θ 0 represents the p×1 "true" or "nominal" parameter vector that generates the observations Y j and we assume this "true" or "nominal" parameter vector exists. The random error (or measurement/observational error) is represented by the m × 1 random vector E j . We assume E j , for j = 1, . . . , n, are independent and identically distributed with zero mean and covariance matrix given by V 0 = Var(E j ) = diag(σ 2 1 , . . . , σ 2 m ), for j = 1, . . . , n. In order to estimate the "true" parameter vector from a set Ω θ of possible parameters, since we are assuming absolute error, we use an ordinary least squares (OLS) approach and minimize the OLS cost functional Standard errors (SE) are then calculated as SE k ≈ Σn kk . We test use of this mathematical/statistical framework for problems with constant parameters.
3. The Verhulst-Pearl logistic growth population model. We first investigate the problem of determining the minimum amount of data that is needed to collect using the simple and well-known Verhulst-Pearl logistic growth population model [16] dx The state variable x(t) represents the population size at time t, parameter K represents the carrying capacity, and parameter r represents the intrinsic growth rate. The exact solution of this model is given by where x 0 represents the scaled (relative to the carrying capacity) initial population size, i.e., x(0) = x 0 . We first created simulated data {y j } n j=1 from equation (4) using the "true" parameter set θ θ θ 0 = [K, r] = [17.5, .7] and with x 0 = 0.1. We then added Gaussian noise to the simulated data with mean 0 and standard deviation 0.01. Using the inverse problem methodology described above, we carried out estimation procedures for the parameters θ θ θ = [K, r] using the ordinary least squares method. Note that we did not try to estimate the initial population size x 0 . This is because when we consider the daphnia population data from [17], the initial daphnia population is known. Since the model is very simple and only two parameters are being estimated, a Nelder-Mead optimization algorithm is sufficient. Thus we used the MATLAB function fminsearch to minimize the cost functional with respect to θ θ θ, where f (t j , θ θ θ) = x(t j ). We provided the MATLAB function with an initial guess θ θ θ 0 = [8,1] and obtained the minimized cost J n (θ θ θ), whereθ θ θ =θ θ θ n represents the optimized value of θ θ θ using n number of points. We then calculated the standard errors for each estimated parameter. Since the model is simple, we calculated the sensitivity matrix by computing the partial derivatives The model sensitivities are plotted in Figure 1. 3.1. Length of time for data collection. We are interested in determining the smallest amount of data that is needed to still provide an adequate parameter estimationθ θ θ. Thus we performed multiple inverse problems with an increasing number of data points. We first use only three data points {y 1 , y 2 , y 3 } to estimatê The simulated data began on time t = 0 and was "collected" each time unit ("day") until time t = 24. Figure 2a contains the graphs of the estimated curves and the data and Figure 2b contains the same curves and data over a smaller time period.
The estimated parameter valuesθ θ θ i and corresponding standard errors obtained for each increasing dataset are plotted in Figure 2c and Figure 2d, respectively.
We can see that parameter r can be estimated with very few data points and the standard errors remain small for parameter r with increasing n. However, K is under estimated while using data up to time t = 3 and then is over estimated. Around day 6 or 7, K becomes close to its true value and the model estimation is adequate. Similarly, the standard error for parameter K is large using data up to time t = 3 and then decreases as more data are used in the estimation procedure. Around day 6 or 7, the standard errors for parameter K begin to level off. It should be noted that this was to be expected due to the model sensitivities with the "true" parameters. We can see that the model is sensitive to r early on, while the model does not become sensitive to K until around day 6 or 7. Thus we can expect to be able to estimate r with very few data points but the data will not contain information about parameter K until at least day 6 or 7.
In this simple example, it is easy to look at all of the data and see when the parameter values and standard errors level off. However, a more important problem is developing an adaptive feedback methodology that can be used as the data are being collected (i.e., establishing "real-time" stopping rules). To do this, we first calculate the percent change in parameter values for an increasing number of data points n, using the following P n − P n−1 P n−1 100, where P n represents the estimated parameter value using n number of data points and P n−1 represents the estimated parameter values using n − 1 number of data points. The percent change in parameter values for increasing n is plotted in Figure  3a. For each parameter (K, r), we then determine the time when the percent change in parameter value is less than a threshold of 5%, i.e., P n − P n−1 P n−1 100 < 5.
We note that 5% was chosen somewhat arbitrarily, but represents a sufficiently small relative error to suggest that the parameter values have leveled off to reflect a stable population. We then choose the maximum of the times associated with each of the two parameters. We plotted a vertical dashed line at this time in Figure  3a. We can see that the percent change in parameter r is above 5% when using up to 4 data points but then below 5% when using up to 5 data points. However the percent change in parameter K is less than 5% once 8 data points are used. Thus a vertical dashed line is plotted at n = 8 (max{5, 8}) in Figure 3a.
We then calculate the percent change in the standard error using the following where s n and s n−1 represent the standard error using n and n − 1 number of data points, respectively. The percent change in standard error for increasing n is plotted in Figure 3b. Again, for each parameter, we determine the time when the percent change in standard error is less than a threshold of 5%, i.e., We then choose the maximum of the times associated with each of the two parameters. We plotted a vertical dashed line at this time in Figure 3b. We can see that the percent change in standard error for parameter r is less than 5% using 7 data points, but the percent change in standard error for parameter K requires 8 data points to reach below the threshold. Thus, a vertical dashed line is plotted at n = 8 (max{7, 8}) in Figure 3b. Note that we want the parameter values to converge to the "true" or nominal parameter value as more and more data are collected and thus we use the percent change in parameter values equation (equation (5)). When using the standard error, we want to find the time when the standard errors are decreasing and begin to level off, even if they are still changing relative to the previous standard error. Thus we calculate the percent change in standard error using equation (6). Using these two methods, we have two different proposed last data points (proposed last data P and proposed last data SE, Figure 3). Since we want both the parameter values to converge and the standard errors to level off, we take the maximum of these times as our final proposed last data point. Thus, in Figure 2, a vertical dashed line is plotted at time t = 7, or equivalently n = 8. Note that in this case, both proposed last data points were equivalent (n = 8). As the data are being collected, the percent change in parameter values and standard errors can be calculated and once both values are below their thresholds, we can be confident that sufficient data has been collected to obtain reliable parameter values with acceptable confidence levels. Our adaptive feedback methodology is summarized in the following steps. Each time data are collected 1. Perform an inverse problem using all the available data to obtain the estimated parameter(s). 2. Calculate the associated standard error(s). 3. Calculate the percent change in parameter values. 4. Calculate the percent change in standard errors. 5. Stop collecting data once the percent change in parameter values decreases below a threshold and the percent change in standard errors decreases below a threshold.

3.2.
Frequency of data collection. When designing an experiment, it might be difficult to collect data daily. We now consider the amount of data needed to collect in order to provide reasonable parameter estimations when the data are collected, for example, every other day. To investigate this question, we repeated the process described above using data from time t = 0 to t = 24 with data being "collected" every other day. The results are plotted in Figures 4 and 5. Again, we calculated a proposed last data point using the percent change in parameter values by determining the time when the percent change for each parameter reached below the threshold and then by choosing the maximum of the times associated with each of the two parameters (using 5 data points, Figure 5a). We then calculated a proposed last data point using the percent change in standard error by determining the time when the percent change in standard error for each parameter reached below the threshold and then by choosing the maximum of the times associated with each of the two parameters (using 6 data points, Figure 5b). Finally, we chose the maximum of these two threshold times as our final proposed last data point (using 6 data points). Thus, we plotted a vertical dashed line in Figure 4 at time t = 10, or equivalently n = 6, to represent our proposed time when enough data has been collected. A benefit here is that only 6 data points are needed whereas about 8 data points were needed when collecting the data daily. However, for this example, by collecting the data daily, we would have sufficient data a few days earlier. 4. Daphnia population growth. Using our knowledge from this simple example, we now consider the question of whether partial demographic data can replace full demographic data while still providing a reasonable estimation of the stable population growth rate λ [9, p.81-89] using the daphnia population data from Stark and Banks [17].  represents the stable population growth rate [9], or multiplication rate of the stable population, i.e., where x(t) represents a daphnia population at time t and x 0 represents the initial daphnia population (x(0) = x 0 ). Since we know that each replicate of each species began with 10 daphnia, then the model becomes The life tables that Stark and Banks developed are based on the stable population model [9], which assumes the following: 1. The only changes in the population are from births and deaths (i.e., there is a closed population) 2. The population is growing at a constant rate over time (i.e., the birth and death rates are constant) If we tacitly assume 1. and 2., as done in [17], then we can let b represent the constant per capita birth rate, d represent the constant per capita death rate, and r represent the constant intrinsic rate of increase, where r = b − d. Then and so  Note that equation (9) is equivalent to equation (7) with λ = e r . For further information on the stable population model and the role of this model in developing the life tables, see Appendix A in [4].

4.2.
Adaptive feedback methodology on the parameter λ. Based on the results from our simple example using the logistic population growth model, we investigate the minimum amount of daphnia population data needed by estimating λ, the parameter of interest, and by calculating the standard errors for an increasing number of data; once the estimated parameter values and standard errors level off, we will deem that sufficient data will have been collected. We can determine this by using the methodology described above. We first calculate the percent change in parameter values for increasing n using equation (5) and determine when this reaches below a threshold. We then repeat this process for percent change  in standard error using equation (6). Our proposed last data point collection is then the maximum of these two times. However, there is a difference between our simulated data in the logistic model example and the data from Stark and Banks. In the logistic model example, the state variable represented total population size and the data was total population size at given times. For the daphnia population model described above, the state variable also represents total population size, however the data are not total population size. Since the offspring are removed daily, the data consists of only two generations, not the total population. Since the life tables were built on the data (see Appendix A in [4] for how to develop the life table), we used the calculated λ from each week of the life table to determine the total population for each week using equation (8). We repeated this process for each replicate of each species. Thus, the datasets {y j } n j=1 now consist of the total population for each week.
We use the inverse problem methodology to carry out estimation procedures for θ = λ using the ordinary least squares method. Again, we did not try to estimate x 0 since we know x 0 = 10. Using the MATALB function fminsearch, we minimized the cost functional with respect to θ, where f (t j , θ) = x(t j ). We provided the MATLAB function with an initial guess which was equivalent to the λ calculated from the life table for that week, and obtained the minimized cost J n (θ), whereθ =θ n represents the optimized value of θ using n number of points. To compute the standard errors for λ, we used the partial derivative Following the procedure in the logistic model example, we performed multiple inverse problems to estimate θ = λ for an increasing number of data points. We first  Figure 6. C. dubia number of offspring and survival data were used to build a life table [17]. The calculated λ from the life table for each week was then used to calculate the total population each week. Using this total population data, parameter valuesθ n are estimated each week and the resulting curves are potted along with the total population data (1st row). The estimated parameter valuesθ n are plotted as well as the standard errors forθ n for each replicate (2nd row). The dashed vertical line represents the proposed last data by Stark and Banks. The dash-dot vertical line represents our proposed last data.
used only two data points {y 1 , y 2 } to estimateθ 2 and then repeated this process until all n data points had been used, resulting in the set {θ 2 , . . . ,θ n }. We then plotted the total population data (calculated from the λ for each week from the life tables) along with the estimated curves over the entire time period using eacĥ  errors, and percent change in standard errors were also plotted. This process was then repeated for each replicate and for each species. Appendix B in [4] contains all of the graphs for each replicate of each species. Figure 6a contains the graphs of the estimated curves and the total population data for the first replicate for C. dubia (the other replicate curves and data are similar, see Appendix B in [4]). The estimatedθ n values for each replicate of C. dubia are plotted in Figure 6b, and the associated standard errors are then plotted in Figure 6c. Using the Student-Newman-Keuls test, the results from Stark and Banks indicated that the λ calculated using up to 5 weeks of C. dubia data was not statistically different from the λ calculated using all 9 weeks of data and thus they suggested that 5 weeks of data was sufficient. We plotted a vertical dashed line in Figure 6 at week 5 to represent their proposed last data point.
To determine our proposed last data point, we first calculated when the percent change in parameter values (equation (5)) was less than 5% and when the percent change in standard errors (equation (6)) reached below 25% (Figure 7). Note that we used a larger threshold for the percent change in standard error since we are interested in when the standard errors are decreasing and begin to level off, even if they are still changing relative to the previous standard error. Our proposed last data point for each replicate was calculated as the maximum of these two times. This process was repeated for each replicate of each species (see Appendix B in [4]). To determine the proposed last data for each species, we chose the proposed time for each replicate that occurred the most number of times. For example, for C. dubia, the proposed last data for replicates 1, 2 and 3, occurred at 8 weeks, 8 weeks, and 7 weeks, respectively (Appendix B in [4]). Thus, we chose 8 weeks as the proposed last data point for species C. dubia using our proposed methodology.
When considering the D. pulex data, Stark and Banks again concluded that using 5 weeks of data provided a sufficient estimate of the population growth rate.  [17]. The calculated λ from the life table for each week was then used to calculate the total population each week. Using this total population data, parameter valuesθ n are estimated each week and the resulting curves are potted along with the total population data (1st row). The estimated parameter valuesθ n are plotted as well as the standard errors forθ n for each replicate (2nd row). The dashed vertical line represents the proposed last data by Stark and Banks. The dash-dot vertical line represents our proposed last data.
Our results indicate that 8 weeks of data are needed ( Figure 10). Thus, we again propose about 3 more weeks compared to Stark and Banks' proposed time. When considering the D. magna data, Stark and Banks concluded that using 7 weeks of data provided a sufficient estimate of the population growth rate. Our results again indicate that 8 weeks of measurements are needed (Figure 8). For each species, our results consistently indicated that 8 weeks of data were needed in order to provide a reasonable estimate for the population growth rate λ.  We should note that the percent change in parameter values and standard errors do not always decrease below the threshold and remain there, but sometimes increase back above the threshold. To handle this situation, one must use all of the information available to decide when enough information has been collected. For C. dubia, although the percent change in parameter values begins below the threshold, 2 data points is clearly not enough data (Figure 7). For the percent change in standard error using C. dubia data, replicate 3 drops below the threshold at week 4, however replicates 1 and 2 are still above the threshold at this time. Thus, we would not propose to stop collecting data at week 4. Although all of the replicates are below the threshold at week 5, replicate 3 increased from about 1.5 to 23 percent, which is quite a big change. Thus, one more week of collecting data might be recommended. At week 6, again replicate 3 increases to about 42 percent, another large change, and replicates 1 and 2 also increase. Thus, our proposed last week of collecting data occurred on week 8, rather than an earlier week. If data had been collected for only replicates 1 and 2, then we would have proposed to stop collecting data on week 5. Similarly, when considering D. magna percent change in parameter values, the levels reach below the threshold on week 3 ( Figure 9). However, there is no indication that the estimated λ values have converged (Figure 8), thus another week would be recommended. When considering species D. pulex, the percent change in standard error decreased below the threshold on week 4 for replicates 2 and 3 ( Figure 11). However, the standard errors at this time are still relatively high. Thus, another week of collecting data would be recommended. Thus, in summary, we believe our proposed methodology can be extremely useful if employed with proper caution and a measure of conservatism when in doubt.
Stark and Banks developed the life tables weekly, even though the data was collected daily. We decided to investigate the results using daily life tables. Once  Figure 10. D. pulex number of offspring and survival data were used to build a life table [17]. The calculated λ from the life table for each week was then used to calculate the total population each week. Using this total population data, parameter valuesθ n are estimated each week and the resulting curves are potted along with the total population data (1st row). The estimated parameter valuesθ n are plotted as well as the standard errors forθ n for each replicate (2nd row). The dashed vertical line represents the proposed last data by Stark and Banks. The dash-dot vertical line represents our proposed last data.
we created the daily life tables, we repeated the same process using the daily λ values from the life tables ( Figure 12). Again we plotted a dashed vertical line to represent the last data point that was suggested by Stark and Banks. Although we did not calculate the percent change in parameter values and percent change in standard errors using the daily life tables, we plotted a dash-dot vertical line to represent our proposed last data point calculated from the weekly life tables. Using the daily data, we can see that the changes in parameter values and standard errors are much

5.
Discussion. An important question in the field of ecotoxicology is determining whether or not partial demographic data can replace full demographic data while still obtaining an accurate estimate of the endpoint of interest. To investigate this problem, we used a simple logistic model in order to develop an adaptive feedback methodology that can be used as the data are being collected for determining when enough data has been collected to obtain reasonable parameter estimates. By carrying out an inverse problem each time data are collected, one can plot the estimated parameters over time, standard errors over time, and models using each estimated parameter. Once the parameter estimates converge (measured by percent change in parameter estimates decreasing below a threshold) and the standard errors level off (measured by percent change in standard error decreasing below a threshold), one can have confidence in concluding that a sufficient amount of data has been collected to still provide reasonable parameter estimates. We note that the methodologies developed here are meant to provide "real-time" stopping rules for data collection for stable population models. In addition, they suggest a means for optimal design of experiments when investigating known stable assemblages of species. Although our results supported the results by Stark and Banks, the conclusions reached should be taken with reservations since a simple model using constant birth and death rates was assumed in both studies. Previous works that support using time-varying coefficients when modeling daphnia populations [2,1] and other populations [6,5,3] suggest future efforts with time-varying coefficients should be carried out to further develop our ideas. We note that our methodological   developments are founded on having both a correct mathematical as well as a correct statistical model for the data collection. Current efforts include testing models with better fits to data to verify that the methodology proposed here is also useful in models that describe the data sets more closely.