Elizabeth E. Holmes and William F. Fagan. 2002. Validating population viability analysis for corrupted data sets. Ecology 83:2379-2386.

A pdf version of Appendix B is also available for viewing.

Appendix B. Derivations of the distributions of the estimated parameters.

Distribution of the estimate, Eq. 2

Consider for the moment, an idealized estimate using subsampled data to eliminate overlap in the N_t+/N_t ratios and L = 1. Let's call it . We can derive the distribution of by observing that the slope of var(ln(N_{t +} ) - ln(N_t)) vs. ( = 1,2,…,') is basically

since the var(ln(N_t+)-ln(N_t)) vs. line is generally straight. Using

where (,) is a gamma distribution with shape and scale , the distribution of is straight-forward to derive.

(B.1)

Note that the sequential N_t+'/N_t ratios are chosen so that there is no overlap thus each ratio is independent. Assume for the moment, that the two gamma distributions are independent -- which they are not. In this case, we can show as follows that the limiting distribution of Eq. B.1 as df_' and df₁ become large is ² with df_' degrees of freedom.

The moment generating function of is . Thus the moment generating function for the distribution in Eq. B.1 is

.

Take the natural log of this to get,

.

Using the Taylor expansion for ln(1+x) and multiplying the second element by ,

.

Ignoring higher order terms, the ln(mgf) has the form:

which is the ln(mgf) of the following ² distribution:

As noted, the gamma distributions for the variances of ln(N_t+'/N_t) and var(ln(N_t+1/N_t)) are actually correlated. The effect of the correlation, as seen from numerical experiments, is to cause the distribution in Eq. B.1 to approach the limiting distribution faster (i.e., when the dfs in the gamma distributions are smaller).

The used in the Dennis-Holmes method is somewhat different than the idealized used in this derivation. First, the
N_t+'/N_t ratios cannot generally be subsampled due to short time series. This means the ratios are correlated and df_' is substantially less than the number of ratios minus one; additionally the lack of subsampling makes biased. The data are running sum transformed (L > 1); this leads to further bias. These are trade-offs that improve estimation for short corrupted time series by reducing the number of negative variance estimates (percent errors column in Table B1). Despite the differences, understanding the limiting distribution for the idealized helps us understand why when we estimated a non-idealized (L > 1 and data not subsampled) from simulated data, we observed that showed a distribution of the form for a wide range of time series lengths, non-process to process error ratios, and filter lengths (Table B1).

Monte Carlo estimation was used to numerically estimate the ² distributions for the estimates used in the Dennis-Holmes method (= the slope of ln(R_t+ /R_t) vs. for = 1,2,3,4). Monte Carlo estimation uses parameter estimates from samples of data generated with simulations to calculate the distribution of the parameter estimate (this is akin to parametric bootstrapping). We generated 5000 time series of length n using the model, N_t+1 = N_t exp(+ _p), O_t = N_t exp(_np) where the process error, _p ~ Normal(0, _p), and the non-process error, _np ~ Normal(0, _np). Let mean( ) denote the mean of all 5000 estimates. For each simulation, we calculated the statistic = . We then found the best fitting df_slp parameter such that

.

This was done by finding the df_slp that maximized the P value from a Kolmogorov-Smirnov goodness of fit test. The fitting process was repeated for different time series lengths (n), filter lengths (L), ratios of process to non-process error (_p/_np) and and _p. The best fitting df_slp values for different n, L and (_p/_np) are given in Table B1 with the P values for the fitted distribution. The observed bias and parameters from the simulations are given in Table B2. The degrees of freedom depended mainly on the length of the time series, n, and the length of the filter, L. There was an approximately linear relationship between n, L and the df_slp values in Table B1. The following formula gives a close approximation of the numerically calculated df_slp:

df_slp = 0.333 + 0.212 n - 0.387 L   for   n > 15.

Variance of , Eq. 3

Given n observations, O₁, O₂, O₃ … O_n, of the true population size, N₁, N₂, N₃ … N_n, the O_t series is transformed into a running sum, R₁, R₂, R₃ … R_r where r = n-L+1 and .

Denote as the mean of the N's that comprise the running sum, R_t: , and recall O_t = _np,t N_t.

Note that is the mean of the ln(R_t+1/R_t) ratios from the time series; however, for corrupted time series, the variance of the is not 1/(n-L) times the variance of the ln(R_t+1/R_t) ratios, as it would be the case for uncorrupted time series:

Using the variance of the ln(R_t+1/R_t) ratios would lead to high overestimation of the variance of . This overestimation is greater for smaller L.

Estimate of the distribution of from data, Eq. 4

If and were known, it would be straight-forward to specify the distribution of , (i.e., Normal() however, instead we have to use estimates of and which themselves have some distribution. Below is outlined an estimate of the distribution of which uses only . Deriving a distribution based on both and appears problematic given the nature of the distribution of (see below) and given that the estimate is not independent of .

By simple algebra, we can rewrite as

Point estimate of

A point estimate of can be calculated by noting that , thus

where from the data.

Tables

Table B1. Numerically calculated degrees of freedom for the ² distribution describing the ratio and percent errors (negative variance estimates). P values give the fit of the empirical distribution from a Kolmogorov-Smirnov goodness of fit test. For the degrees of freedom calculations, negative were removed from the sample. For the simulations, = -0.01 and = 0.1. The results were not sensitive to alternate parameter values in the ranges: –0.2<<0.2 or 0.01<<0.3.

Table B2. Numerically calculated mean bias between and , expressed as a percentage of , and parameter (Eq. 11) describing the relationship between and . Negative estimates were removed from the sample before calculations.  For the simulations, = -0.01 and = 0.1. The results were not sensitive to alternate parameter values in the ranges: –0.2<<0.2 or 0.01<<0.3.