Monte Carlo Evaluation of the Likelihood for N_e From Temporally Spaced Samples

Corresponding author: Eric C. Anderson, Department of Statistics, University of Washington, Box 354322, Seattle, WA 98195., eriq{at}cqs.washington.edu (E-mail)

Communicating editor: G. A. CHURCHILL

	ABSTRACT

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

A population's effective size is an important quantity for conservation and management. The effective size may be estimated from the change of allele frequencies observed in temporally spaced genetic samples taken from the population. Though moment-based estimators exist, recently Williamson and Slatkin demonstrated the advantages of a maximum-likelihood approach that they applied to data on diallelic genetic markers. Their computational methods, however, do not extend to data on multiallelic markers, because in such cases exact evaluation of the likelihood is impossible, requiring an intractable sum over latent variables. We present a Monte Carlo approach to compute the likelihood with data on multiallelic markers. So as to be computationally efficient, our approach relies on an importance-sampling distribution constructed by a forward-backward method. We describe the Monte Carlo formulation and the importance-sampling function and then demonstrate their use on both simulated and real datasets.

REDUCTIONS in population size can lead to inbreeding, which increases the probability of population extinction in typically outbreeding species (FRANKHAM 1995 ). Reductions in population size also lead to a loss of genetic diversity, which may restrict a population's ability to adapt to changing conditions (SOULE 1986 ). To predict the risk to a population from these types of genetic factors, biologists are often interested in knowing the effective population size, N_e. An effective size is defined by comparison to an ideal population model, the Wright-Fisher model. The Wright-Fisher model assumes discrete, nonoverlapping generations of constant size, and it assumes that the gametes that unite to form adults in one generation are randomly sampled with replacement from the previous generation. The variance effective size of a natural population is the size of a Wright-Fisher population that would experience a comparable increase in variance of gene frequency over time. The inbreeding effective size is defined similarly, but is based on the increase in gene identity by descent over time.

It is possible to estimate the variance effective size from observed changes in allele frequencies in a population over time. Moment-based estimators using F-statistics have been developed for this purpose (KRIMBAS and TSAKAS 1971 ; NEI and TAJIMA 1981 ; POLLAK 1983 ; WAPLES 1989 ; JORDE and RYMAN 1995 ). Recently, WILLIAMSON and SLATKIN 1999 described a method to estimate N_e by the method of maximum likelihood. To find the maximum-likelihood estimate _e of N_e, given allele frequencies observed in samples taken from a population at different times, one models the population underlying the samples as a Wright-Fisher population. _e is then the size of that underlying, ideal population for which the observed data are most probable. In simulation studies WILLIAMSON and SLATKIN 1999 showed that the maximum-likelihood estimator outperformed the moment-based estimators, and they also demonstrated how a likelihood approach may be extended to estimate parameters in more complex population models.

This likelihood method has been restricted to data on diallelic loci, because, with data on multiallelic loci, evaluating the likelihood for N_e exactly is computationally intractable. Here we describe the problem as one of inference from a hidden Markov chain (BAUM et al. 1970 ). We develop an algorithm for importance sampling, which makes it possible to compute the likelihood by Monte Carlo.

	FORMULATION OF THE MODEL AND MONTE CARLO

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

The model:
The data are genetic samples collected at different generations. The first sample is collected at generation 0 and the last sample at generation T. Any samples drawn at intervening generations may be evenly or irregularly spaced in time. For notational simplicity, we assume for now that individuals are genotyped at a single locus, though we describe later the extension to multiple, independently segregating loci. The data include K different allelic types, indexed by k = 1, ... , K. The allele frequencies observed in samples taken from different generations will differ due to genetic drift and sampling variation.

Let Y_t = (Y_t,1, ... , Y_t,K) be the counts of the K different allelic types in the sample at generation t, and let S_t denote the number of diploid individuals in the sample. We assume that the samples were taken from a Wright-Fisher population of size N_e and denote the unobserved population allele counts at generation t by X_t = (X_t,1, ... , X_t,K), with ^K_k=1 X_t,k = 2N_e. By the formulation of the Wright-Fisher model, the X_t form a Markov chain in time, with transitions defined by multinomial probabilities depending on N_e,

(1)

The genetic sample at a time t is assumed to be drawn with replacement from the copies of alleles present in the population at time t. This is equivalent to drawing the sample Y_t from a very large gamete pool produced by the population at time t: sampling plan II of WAPLES 1989 . This type of sampling applies to many organisms, especially those species with high fecundity that may be sampled as juveniles, or those that may be sampled (preferably noninvasively) as adults in populations having census sizes considerably larger than their effective sizes (WAPLES 1989 ). The sample allele counts Y_t, given the latent variable X_t, are conditionally independent of all the other variables and follow the multinomial distribution depending on the parameter N_e, the sample size S_t, and X_t,

(2)

when S_t > 0. If there is no sample taken from the population at generation t, then S_t 0, and we define P_Ne(Y_t|X_t) 1.

Such a system forms a hidden Markov chain with the dependence structure shown in the directed graph of Fig 1. The allele counts in the population when the first sample is drawn, X₀, are nuisance parameters. To avoid estimating X₀ and the associated consistency problems with the maximum-likelihood estimator (NEYMAN and SCOTT 1948 ), we consider the integrated likelihood, assuming a uniform prior distribution, (X₀), on the population allele counts at time 0. The likelihood for N_e is the probability of the data Y = (Y₀, ... , Y_T) given the parameter N_e. The probability of Y is the sum of the joint probability of Y and the latent variables X = (X₀, ... , X_T) over the space of all X:

(3)

View larger version (22K): In this window In a new window Download PPT slide   Figure A1. Figures representing and for 2N* = 20. The normal curve is the density for . (a) Reflections and translations as described in the Appendix Long-dashed lines represent the curve after reflection through L or H, while the short-dashed lines represent the reflected curve after one or more successive translations. (b) If = 1 and = 1, then X(i) is constrained to be in {1 ... , 2N* - 1}. The shaded regions correspond to those values of for which X(i) = 13 by . The total shaded area is equal to µ,2(X = 13; 10, 1, , 1). (c) If = 0 then X(i) may take the value 0. The shaded area shows µ,2(X = 0; 10, 0, , ). (d) If = 0 then X may take the value 2N*. The shaded area shows µ,2(X = 2N*; 10, , , 0).

For the case of K = 2 and N_e small, the likelihood in (3) may be computed exactly. WILLIAMSON and SLATKIN 1999 effected the summation in (3) in terms of multiplication of transition probability matrices. The dimension of the square matrices is (N_e - 1)!/[(N_e - K)!(K - 1)!], which increases rapidly with N_e and K. We note that the hidden Markov form of the system allows a more efficient computation of the likelihood using the algorithm of BAUM 1972 . Nonetheless, exact evaluation for multiple alleles would still require prohibitively large amounts of computation and storage. An alternative is to estimate P_N^e(Y) by Monte Carlo.

Monte Carlo evaluation:
For likelihood inference, we must evaluate P_N^e(Y) for a number of different values of N_e. Expressing this probability as an expectation with respect to the distribution of X gives

(4)

In this form the expectation would be taken over the marginal probabilities of X, and it could be estimated by Monte Carlo as

(5)

for large m, with X⁽ⁱ⁾ being the ith realization from the marginal distribution of X. Such a naive scheme fails, however, because P_Ne(Y|X⁽ⁱ⁾) varies greatly over the values of X realized from their marginal distribution, resulting in enormous Monte Carlo variance.

Instead, we pursue a more efficient Monte Carlo approximation by using importance sampling (HAMMERSLEY and HANDSCOMB 1964 ). We express P_N^e(Y) as an expectation with respect to a different distribution P^*_Ne(X) having the property that P^*_Ne(X) > 0 for all X such that P_Ne(Y, X) > 0. Thus, we have

(6)

where ^*_Ne indicates that the expectation is over the space of X weighted by the distribution P^*_Ne(X). The expectation (6) may be estimated by Monte Carlo, giving

(7)

for large m, where X⁽ⁱ⁾ is the ith realization of X drawn from P^*_Ne(X). The Monte Carlo variance of _Ne(Y) is made small when varies little across the possible values of X and would be minimized if P^*_Ne(X) were exactly proportional to P_Ne(Y, X). Such a distribution of X would, by definition, be the conditional distribution P_Ne(X|Y). Unfortunately, for the same reasons that P_Ne(Y) cannot be computed exactly, it is infeasible to compute P_Ne(X|Y). Nonetheless, the Monte Carlo variance of _Ne(Y) will be reduced to the extent that P^*_Ne(X) resembles P_Ne(X|Y). We now describe a method for rapid simulation of X⁽ⁱ⁾'s from a distribution P^*_Ne(X) that is close to P_Ne(X|Y). As is required for the importance sampling, it is also possible to compute P^*_Ne(X⁽ⁱ⁾) quickly for each X⁽ⁱ⁾ generated.

Sampling from P^*_Ne(X) by a forward-backward method:
BAUM et al. 1970 describe a method applicable to general, hidden Markov chains for realizing latent variables, such as X = (X₀, ... , X_T), from their exact conditional distribution given the observed variables, such as Y = (Y₀, ... , Y_T). Their algorithm first employs a "forward step" in which the conditional probability distributions of each X_t, given the observed variables up to and including Y_t, are recursively computed and stored using the relation

(8)

which is normalized by the sum of that quantity over all the values of X_t. The last such conditional distribution computed is P(X_T|Y₀, ... , Y_T). The "backward step" begins with simulating a value X⁽ⁱ⁾_T from this distribution (where, as before, the superscript (i) indicates a realized value of a random variable). One then proceeds backward, realizing X⁽ⁱ⁾_T-1 from its conditional distribution given all of the observed variables, Y, and X⁽ⁱ⁾_T. In similar fashion, one realizes X⁽ⁱ⁾_T-2 and so forth back to X⁽ⁱ⁾₀. In this backward phase, each X⁽ⁱ⁾_t is simulated from its conditional distribution given all the data Y and all of the components of X that have been realized so far. That is, X⁽ⁱ⁾_t is drawn from

(9)

Because of the conditional independence structure in a hidden Markov chain, (9) reduces to P(X_t|Y₀, ... , Y_t, X⁽ⁱ⁾_t+1), which may be computed using the distributions stored during the forward step by the relation

(10)

At the end of the backward step, it is thus clear that the resulting realization (X⁽ⁱ⁾₀, ... , X⁽ⁱ⁾_T) is from the conditional distribution of X given Y.

An approximation for multiple alleles: In our application, with multiple alleles at a locus, since there are so many possible states that each X_t may take, the above procedure is computationally infeasible. However, we make use of the BAUM et al. 1970 algorithm in spirit, employing two alterations to make it feasible to simulate from P^*_Ne(X) and to compute its value. We emphasize that although the method described below involves a series of "approximations" by which P^*_Ne(X), differs from P_Ne(X|Y), the final sampling and computation of P^*_Ne(X) is exactly from the P^*_Ne(X) as constructed, so its use in (7) gives a true Monte Carlo estimate.

The first approximation is to perform the forward-backward cycle separately for each allele. To describe this, we introduce some more notation. Denote by X_(k) the vector (X_0,k, ... , X_T,k) of latent counts of the kth allele from time t = 0 to t = T. Similarly we define Y_(k) = (Y_0,k, ... , Y_T,k). To do the forward-backward cycle separately over alleles we first focus on allele 1, simulating X⁽ⁱ⁾₍₁₎ by the forward-backward mechanism as if the data were on a diallelic locus with observed counts Y₍₁₎ from samples of size S₀, ... , S_T through time. Once we have realized X⁽ⁱ⁾₍₁₎ we update the sizes of the population and the sample. Thus we define the updated population size vector 2N^*₍₂₎ = (2N_e - X⁽ⁱ⁾_0,1, ... , 2N_e - X⁽ⁱ⁾_T,1) and an updated sample size vector 2S^*₍₂₎ = (2S^*_0,2, ... , 2S^*_T,2) = (2S₀ - Y_0,1, ... , 2S_T - Y_T,1), in effect removing the first allelic type from the remainder of the data and the population. We then use the forward-backward mechanism again to simulate X⁽ⁱ⁾₍₂₎ as though the data were counts Y₍₂₎ from a diallelic locus drawn from a population with sizes that change over time N^*₍₂₎ and sample sizes S^*₍₂₎. This continues sequentially over alleles, updating population sizes and sample sizes as above: 2N^*_(k) (2N^*_(k-1) - X⁽ⁱ⁾_(k-1)) and 2S^*_(k) (2S^*_(k-1) - Y⁽ⁱ⁾_(k-1)), until X_(K-1) has been realized, which also determines that X_(K) (2N^*_(K-1) - X⁽ⁱ⁾_(K-1)). (Here and later we use the notation A B to mean "the value B is assigned to the variable A.") At the end one has obtained a realized value X⁽ⁱ⁾, which may be used in (7).

P^*_Ne(X) using a continuous approximation: Although realizing alleles sequentially, as above, greatly reduces the number of terms required to use (8) and (10), the method would still involve a prohibitive amount of summation over binomial probabilities. Thus we construct P^*_Ne(X), employing a normal approximation to binomial probabilities, which replaces all such sums by analytically tractable integrals. Recall that if W Binomial(n, p), then the transformed variable sin^-1(W/n)^1/2 is approximately normally distributed with variance 1/(4n). Note that this quantity does not depend on p. Hence we use this transformation to define the quantities _t,k = sin^-1[] when S^*_t,k > 0, and _t,k = sin^-1[]. By realizing the continuous values ⁽ⁱ⁾_t,k in a forward-backward framework within a continuous setting, the computational demands are greatly reduced. And then, by transforming each ⁽ⁱ⁾_t,k back into the appropriate discrete X⁽ⁱ⁾_t,k we have a way to realize X⁽ⁱ⁾ from P^*_Ne(X) and to compute the probability P^*_Ne(X⁽ⁱ⁾). The details of this procedure are given in the Appendix. We use it to compute the Monte Carlo estimate _Ne(Y) using (7).

Monte Carlo variance and multiple loci:
The quantity _Ne(Y) is only an estimate of the true value P_Ne(Y). By the central limit theorem, for large m, P_Ne(Y) will be approximately normally distributed (HAMMERSLEY and HANDSCOMB 1964 ) with mean _Ne(Y) and a variance that may be approximated without bias by the quantity

(11)

These facts may be used to obtain a confidence interval estimate around each _Ne(Y).

The ability to estimate N_e typically requires data from many loci. The extension to data on J independently segregating loci, indexed by j = 1, ... , J, is straightforward—each locus is treated separately, and the estimated likelihoods from each locus are multiplied together. Thus, let _Ne,j(Y) be the Monte Carlo likelihood estimate from the data on the jth locus. The Monte Carlo likelihood estimate using all the loci is then

(12)

This requires that the initial allele counts have independent prior distributions, (X₀). An implicit assumption of (12) is consequently that the loci used are in linkage equilibrium at t = 0. ^J_Ne(Y) will also have an approximately normal distribution. An unbiased estimator for its Monte Carlo variance (derived in Equation A9 in the Appendix) is

(13)

This can be used to compute a confidence interval estimate around ^J_Ne(Y).

When displaying the Monte Carlo likelihood curve it is preferable to plot the log-likelihood values, log ^J_Ne(Y), for different values of N_e. In this case, the endpoints of the confidence intervals may be similarly log-transformed.

	SIMULATED AND REAL DATASETS

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

We demonstrate the method by computing log-likelihood curves for N_e from three different datasets. First, to verify that the method gives correct results we apply it to a simple simulated dataset (dataset 1) for which it is possible to compute the likelihood exactly. We simulated samples of 20 diallelic loci from 100 diploid individuals at generations 0, 6, and 12 from a Wright-Fisher population of 25 diploid individuals. For each locus, the initial allele frequency in the population at time zero was an independently drawn uniform real number between 0 and 1. The log-likelihood for N_e given these data was estimated for values between 10 and 52, in steps of 2, using m = 20,000 realizations of X from P^*_Ne(X) for each locus and each N_e.

We simulated a second dataset (dataset 2) to see how the method performed with multiallelic markers taken from a Wright-Fisher population. The dataset included three samples of 100 diploids for 12 five-allele loci at generations 0, 4, and 8 from a population of 50 diploids. The allele frequencies at each locus in generation 0 for these simulations were independently drawn from a uniform Dirichlet density with five components. For this dataset, the log-likelihood was computed for values of N_e between 20 and 100 in increments of 4 using m = 50,000 realizations of X for each locus and each value of N_e.

Finally, we computed a log-likelihood curve for N_e given data on a population of Drosophila in BEGON et al. 1980 . These data were analyzed using F-statistics by BEGON et al. 1980 as well as by POLLAK 1983 . They observed allele frequencies in three samples at each of nine enzyme loci. The first two samples were taken a little more than 1 yr apart, and the third sample was taken some 8 mo later. Though the natural populations do not have discrete generations, they have been modeled previously by Begon et al. and Pollak as populations with discrete generations. Because of the different growth rates of flies during different seasons, seven generations separate the first two samples, while only two generations separate the second two samples (BEGON et al. 1980 ). The sample sizes for all loci were the same, with larger sample sizes taken in the latter sampling periods. The sample sizes were S₀ = 190, S₇ = 250, and S₉ = 335 flies. POLLAK 1983 notes that since BEGON et al. 1980 sampled adult flies, their sampling scheme is closer to what is known in the literature as sampling scheme I than it is to sampling scheme II. However, as discussed by WAPLES 1989 , the probability models underlying the two different sampling schemes are very similar when the actual size of the population is much larger than the effective size of the population. This is the case with these Drosophila. BEGON et al. 1980 report census sizes in the tens of thousands of flies, while the estimated N_e is orders of magnitude smaller. Because of this, it is still reasonable to analyze the data using the likelihood method we have developed here.

The data appear as allele frequencies in Table 1 of BEGON et al. 1980 . Unfortunately the allele frequencies at the Pgm locus are misreported there and fail to sum to one. We thus used only the remaining eight loci. Of these eight, three had three alleles, two had four alleles, two had five alleles, and one had six alleles. We evaluated ^J_Ne(Y) at values of N_e between 200 and 1200 in increments of 50, with two more points (N_e = 425 and N_e = 475) included near the peak of the likelihood curve. For each point we used m = 500,000 realizations of X.

	RESULTS

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

For each of the datasets, we were able to use our importance-sampling method to compute a log-likelihood curve. Using a program we wrote in C, the runs for datasets 1 and 2 each took 10 hr on a 266-Mhz laptop computer. The log-likelihood curves from datasets 1 and 2 appear as solid lines in Fig 2. The estimated 90% confidence intervals around each value of log ^J_Ne(Y) appear as two dashed lines bordering the log-likelihood curve. Despite the fact that few Monte Carlo replicates (m = 20,000 and 50,000) were used, the Monte Carlo variance is minimal, as indicated by the fact that the dotted lines practically lie on top of the estimated log-likelihood curve. In both cases, the true values of N_e (25 and 50, respectively) are well within 2 units of log-likelihood from the maximum-likelihood estimates, which may be read from the graph as 24 and 56. Since dataset 1 consists only of diallelic loci, it is possible to compute the exact log-likelihood curve. This exact curve has been plotted as a dotted line in Fig 2A. It is impossible to distinguish the exact curve because the Monte Carlo estimate is very accurate in this case.

View larger version (16K): In this window In a new window Download PPT slide   Figure 2. Log-likelihood curves estimated by Monte Carlo from datasets 1 and 2. The values of Ne at which the likelihood was computed are indicated by vertical lines above the horizontal axis in each figure. The log-likelihood values are connected by a solid line. Vertical bars intersecting the solid line indicate 90% confidence intervals around log JNe(Y) computed using the Monte Carlo variance estimate (13). The endpoints of the confidence intervals are connected by dashed lines. (a) Dataset 1 is simulated data from 20 diallelic loci. (b) Dataset 2 is simulated data from 12 loci with five alleles each.

The log-likelihood curve computed for the data of BEGON et al. 1980 is shown in Fig 3. It took 54 hr on a 450-Mhz desktop computer to produce the results. As before, the 90% confidence intervals around the Monte Carlo estimates appear as dotted lines. With this dataset, even with m = 500,000 realizations of X, the Monte Carlo variance is not negligible. It is, however, small enough that reliable inferences may be made from the log-likelihood curve. The maximum-likelihood estimate of N_e is 500. Using the values of N_e at which the log-likelihood has decreased 2 units from its maximum gives an estimate of a 95% confidence interval for the true N_e. These points are 250 and 975. By contrast, POLLAK 1983 , using an F-statistic method, estimated N_e to be 251 with a standard error of 115. We discuss the discrepancy between the two estimates in the next section. Our results are not comparable to the N_e estimated by BEGON et al. 1980 because, at the time, those authors were unable to make a single estimate of N_e using the samples at all three time points.

View larger version (24K): In this window In a new window Download PPT slide   Figure 3. Log-likelihood curves from the data of BEGON et al. 1980 estimated by Monte Carlo. The format of the plot is as for Fig 2.

	DISCUSSION

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

As discussed in WILLIAMSON and SLATKIN 1999 , estimating N_e by maximum likelihood has advantages over estimating N_e using F-statistics. Until now, it was impractical to compute the likelihood for N_e using all the data when loci with more than two alleles were available. While it has been suggested that one may bin low-frequency alleles together to turn multiallelic loci into apparently diallelic loci and then apply exact likelihood calculation methods to such reduced data, this invariably throws away some information. Furthermore, different binning strategies lead to different results. Allowing full use of the data, the Monte Carlo likelihood procedure described here is a preferable way to analyze temporal data on multiallelic loci. The method is suitable for multiallelic loci such as the microsatellite markers becoming available in a wide variety of species.

Monte Carlo methods use realizations of random variables to estimate an expectation by a sample average. There are a number of ways one can express the likelihood of N_e as an expectation, and then estimate it by Monte Carlo, but few of those schemes will be successful, because most will have high Monte Carlo variance. We attempted several different schemes before settling on the importance-sampling method presented here. Although these less sophisticated Monte Carlo estimators produced reasonable estimates in very small problems, when applied to data involving loci with many alleles these methods failed to converge reliably, even after many days of computation (E. C. ANDERSON and E. G. WILLIAMSON, unpublished data).

The importance-sampling method we use is successful because our importance sampling function, P^*_Ne(X), closely resembles P_Ne(X|Y), the conditional probability of X given Y. This is achieved by recognizing the hidden Markov chain structure of the problem and using the forward-backward algorithm of BAUM et al. 1970 . In doing so, we have developed a Monte Carlo estimator with demonstrably small Monte Carlo variance. Though the computational demands of this procedure are substantial, the reduction in Monte Carlo variance obtained makes it worthwhile. Nonetheless, it may be possible to improve the estimates by making additional changes to P^*_Ne(X) so that it more closely resembles P_Ne(X, Y), especially in the tails of the distribution. This would further reduce the Monte Carlo variance.

It should be pointed out that while many Monte Carlo problems involving high-dimensional random variables like X make use of Markov chain Monte Carlo (MCMC) methods, our method is not an MCMC method. In MCMC, successive realizations are correlated. In our method each X⁽ⁱ⁾ realized from the distribution P^*_Ne(X) is independent of the other realized values. As a result, our method does not have the same problems of convergence assessment as does Markov chain simulation (GELMAN 1996 ).

It is interesting that our maximum-likelihood estimate differs so much from the estimate given by POLLAK 1983 for the same data. Though perhaps some of this is attributable to the fact that we chose not to use the incorrectly reported data at the Pgm locus, there are differences between the two estimation methods that could also account for some of the discrepancy. The most notable differences occur when combining information from multiple samples in time. Consider the fact that a better estimate of N_e may be made from two samples taken many generations apart than from two samples separated by fewer generations. Likewise two large samples will yield a better estimate than two small samples. When there are many samples, the relative information content in different intersample intervals will depend on the relative sample sizes and the number of generations between the samples. By its nature, the maximum-likelihood approach will appropriately weight information from different intervals. In contrast, Pollak's F-statistic, F_Kr, neither includes terms for sample size nor interval length between samples, and, Ñ_Kr, his estimate of N_e based on F_Kr, includes a term for only the number of generations between the first and the last sample and is invariant to permutations of the sample sizes at different times. Since the data from BEGON et al. 1980 span sampling intervals of different lengths and include different sample sizes at different times, differences between our results and those in POLLAK 1983 should be expected.

The Monte Carlo variance of our estimate of the likelihood given the data from a natural population of Drosophila was higher than the variance associated with our estimates from simulated data. Although a good estimate was achieved after sufficient computation, it may still be that data generated under a model that differs from the Wright-Fisher model present difficulty for the Monte Carlo likelihood method. For example, it may be that the effective size of the natural Drosophila population was different during the two different sampling intervals. We note that one could extend the likelihood framework to allow for N_e changing over time. For example, if the estimated census size of the population were available and was known to change over time it would be more sensible to estimate directly the ratio, , of the effective size of the population to the census size of the population. This ratio would be more useful for the purposes of modeling genetic change in the population than a single estimate of N_e over the entire time period between the first and last samples. The forward-backward approach implemented here could easily be modified to estimate this parameter, .

It would also be possible to extend the present approach to compute likelihoods from different stochastic models of allele frequency change. We suspect there would seldom be enough information in the data to estimate accurately models in which allele frequency change is due jointly to genetic drift and some other genetic mechanism such as mutation or selection. However, our methods could be modified to handle different demographic models. For example, consider an alteration of the Wright-Fisher model such that the current generation is formed by sampling genes without replacement from a gamete pool in which each parental allele is represented by M gametes. In such a case, allele frequency changes occur due to multivariate hypergeometric sampling determined by the parameters N and M. For values of M that are not very small, the normal distribution is still a reasonable approximation to the hypergeometric distribution, and importance sampling could proceed as before. The importance-sampling distribution P^*_N,M(X) should then account for the decrease in variance (by a factor of 2N[M - 1]/[2N M - 1]) of the hypergeometric distribution relative to the binomial distribution. The same sort of sampling model and adjustment in the importance-sampling distribution could be used to accommodate scenarios in which one's genetic samples (the data Y) were assumed drawn without replacement from the population.

Another possible extension would be to stochastic models involving populations of organisms with more complex life histories, for example, overlapping generations or age-structured populations. As it becomes easier to gather extensive genetic data on populations, and as understanding of the structure within those populations improves, it will be possible to specify much richer probability models for temporal population genetic data. The forward-backward method here, and variations of it, should be useful in estimating population parameters from such models using Monte Carlo likelihood.

A software package, MCLEEPS, implementing the algorithm described in this article is available for free download from http://www.stat.washington.edu/thompson/Genepi/Mcleeps.shtml.

	ACKNOWLEDGMENTS

We thank an anonymous referee for helpful comments on the manuscript and for suggesting the extension to the population model with hypergeometric sampling. This study was supported by National Science Foundation grant BIR-9807747 to E.T. Additionally, E.A. was supported by National Science Foundation grant BIR-9256537, the University of Washington QERM program, and a Burroughs Wellcome Fund PMMB training fellowship. E.W. was supported by National Institutes of Health grant GM40282 to M. Slatkin.

Manuscript received December 27, 1999; Accepted for publication August 7, 2000.

	APPENDIX

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Using _t,k and _t,k in a continuous setting:
We define the random variables _t,k = sin^-1[] when S^*_t,k > 0, and _t,k = sin^-1[ when N^*_t,k > 0. These quantities have an approximate normal distribution, which is independent of their means. We use them in our construction of the importance-sampling function P_Ne(X). Below, we concentrate on their use for realizing X⁽ⁱ⁾_(k), keeping in mind that if k > 1 then we will have already realized X⁽ⁱ⁾_(k-1), and we use the updated population and sample sizes N^*_(k) and S^*_(k). If k = 1 then N^*₍₁₎ = (N_e, ... , N_e) and S^*₍₁₎ = (S₀, ... , S_T), respectively.

The forward step:
Following CAVALLI-SFORZA and EDWARDS 1967 , if _t-1,k is normally () distributed with mean µ_t-1 and variance ²_t-1, then, after a generation of genetic drift in a population of N^*_t,k diploids, _t,k has an approximate normal distributions with mean µ_t-1 and variance ²_t = ²_t-1 + . If there are data Y_t,k from a sample of size S_t,k at time t, then _t,k has an approximate normal distribution with mean _t,k and variance , so, given that _t,k (µ_t, ²_t), the conditional distribution of _t,k given _t,k is also normal. These relations form the basis of a continuous approximation for doing the forward step. For the purpose of realizing X we assume that the uniform prior on X₀ is equivalent to a diffuse prior on _0,k. Therefore _0,k |_0,k (µ₀, ²₀) with µ₀ = _0,k and ²₀ = . With that as a starting point, we work iteratively forward in time, assigning values

(A1)

(A2)

if S^*_t,k = 0. If S^*_t,k > 0, however, then one first computes µ_t and ² as in (A1) and (A2), but then further updates the values to reflect the information in the sample at time t:

(A3)

(A4)

This is analogous to computing a posterior distribution from a normal prior and normal data (see, for example, GELMAN et al. 1996 , p. 43).

Carrying this out until t = T gives values for the mean and variance of _T,k given _0,k, ... , _T,k, assuming they follow a normal distribution. In fact, for each t, it gives us the parameters for the normal distribution of _t,k conditional on _r,k for all r t. We are thus in a position to realize ⁽ⁱ⁾_t,k's in the backward step and transform those ⁽ⁱ⁾_t,k's back into the X⁽ⁱ⁾_t,k's that we need.

The backward step:
The backward step is more complicated than the forward step, because after realizing each value of ⁽ⁱ⁾_t,k we must transform it into the discrete value X⁽ⁱ⁾_t,k that we require. This transformation process requires some extra bookkeeping to ensure that we do not waste time realizing X⁽ⁱ⁾'s that are incompatible with the data. This is described in the next section of the APPENDIX. We first realize the value ⁽ⁱ⁾_T,k from a (µ_T, ²_T) distribution. Then we transform that to the realization X⁽ⁱ⁾_T,k by a many-to-one map , which has two effects: the first is that of folding and translating the distribution of _T,k so that it is bounded between 0 and /2, mapping ⁽ⁱ⁾_T,k (-, ) to a value ^*_T,k [0, ]. The second is transforming that ^*_T,k into the appropriate value X⁽ⁱ⁾_T,k (see the next section).

Working backward, each ⁽ⁱ⁾_T,k for t = T - 1 down to t = 0, is realized from a (µ_t, ²_t) distribution and then transformed into the corresponding ^*_T,k and X⁽ⁱ⁾_t,k by . In keeping with (10), before ⁽ⁱ⁾_t,k is realized, µ_t and ²_t must be appropriately updated, on the basis of the values of µ_t and ²_t stored during the forward step and the realized value ⁽ⁱ⁾_t+1,k. This involves making the assignments

(A5)

(A6)

in the order as written.

Computing the probability P ^*_Ne(X⁽ⁱ⁾):
By carrying out the forward-backward steps above on the first allele, the realization X⁽ⁱ⁾₍₁₎ is obtained. Then, N^*₍₂₎ and S^*₍₂₎ are computed and used in the forward-backward steps to obtain X⁽ⁱ⁾₍₂₎. Executing these steps for all the alleles yields the realization X⁽ⁱ⁾, which is used in (7). P_Ne(Y, X⁽ⁱ⁾) in (7) is easily computed using the expansion shown between the large parentheses in (3).

It remains only to compute P^*_Ne(X⁽ⁱ⁾), which can be done by recording the probability of realizing each component X⁽ⁱ⁾_t,k. Although this probability depends on the values of µ_t, ²_t, N^*_(k), and several bookkeeping variables, we denote it here simply by (X⁽ⁱ⁾_t,k). (The actual function is described later in this APPENDIX.) As long as the realization of X⁽ⁱ⁾_(k) over alleles occurs in the same order over k (k = 1, 2, ... , K) for each i, then

(A7)

Details of :
The fact that we are realizing X⁽ⁱ⁾_(k)'s one allele at a time requires that we do some extra bookkeeping to keep our importance-sampling scheme efficient. Primarily, we must avoid realizing X⁽ⁱ⁾'s for which P_Ne(Y, X⁽ⁱ⁾) = 0. Potential problems arise because by the method we use to realize values from P^*_Ne(X), X_t,k may only take values between 0 and 2N^*_t,k, inclusive. If 2N ^*_t,k = 0 at any value of t, then for any s > t, X⁽ⁱ⁾_s,k must also be 0. To avoid situations in which this leads to P_Ne(Y, X⁽ⁱ⁾) being 0 (as when X⁽ⁱ⁾_t,k = 0 and Y_t,k > 0) we introduce the following scheme and additional notation:

(A8)

Knowing the above quantities, we can define the function . In the remainder of this section and in the following one we drop the t and k subscripts for clarity.

With N* and positive integers, {0, 1}, and {0, 1, ... , min(2N* - , - )}, let (; N*, , , ): ¹ {, ... , 2N* - } x [0, /2] be the many-to-one map that takes a realization of (-, ) to the ordered pair (X, *), where X is an integer such that X 2N* - , and * is a real number between 0 and /2, inclusive. may be described by the following pseudocode. We first define the quantities L = sin^-1() and

Then,

if ( = 2N* - or = - = 0) then * 0 else if (L < < H) then * else if ( < L)

and if ( = 0) then *

else if ( = 1) then _[L] 2L - (this is reflection around = L), and then

if (L _[L] < H) then * _[L]

else we know _[L] H, and we consider the sequence _[i] = i(L - H) + _[L], i = 1, 2, ... , and we assign * _[i*], where i* is the least i such that L < _[i] < H. (The sequence _[i] represents successive translation leftward.)

else if ( > H)

and if ( = 0) then * /2

else if ( > 1) then _[H] 2H - (this is reflection around = H), and then

if (L < _[H] < H) then * _[H]

else we know _[H] < L and we consider the sequence _[j] = j(H - L) + _[H], j = 1, 2, ... , and we assign * _[j*], where j* is the least j such that L _[j] < H. (The sequence _[j] represents successive translation rightward.)

finally we use *, making the assignment X 2N*sin²* + 0.5,

where x denotes the largest integer x. The reflections and translations are depicted graphically in Fig 1A.

The probability _µ,²(X = X⁽ⁱ⁾; N*, , , ) of realizing X = X⁽ⁱ⁾:
If (µ, ²), and (X, *) = (; N*, , , ), then we denote by _µ,²(X = X⁽ⁱ⁾; N*, , , ) the marginal probability that X = X⁽ⁱ⁾. The value of _µ,²(X = X⁽ⁱ⁾; N*, , , ) can be expressed using the notation from the above section. First, _µ,²(X = X⁽ⁱ⁾; N*, , , ) = 0 if X⁽ⁱ⁾ < or X⁽ⁱ⁾ > 2N* - , though such values of X⁽ⁱ⁾ should never occur from anyway. Second, there are cases when constrains X⁽ⁱ⁾ to be either 0 or 1 with probability 1. Hence if 2N* - = or - = = 0 then _µ,²(X = ; N*, , , ) = 1.

If, on the other hand, 2N* - > and - > 0, then for X⁽ⁱ⁾ = 0 and X⁽ⁱ⁾ = 2N* we have

while for 0 < X⁽ⁱ⁾ < 2N* - we define a = sin^-1[] and b = sin^-1[] and have

where I{·} is the indicator function and P(a < b) is the probability that a (µ, ²) random variable is between a and b, namely ^b_a(2²) exp{}d. We compute this probability by numerical integration in our programs. In practice, the infinite sums are approximated by summing the first several terms of the series, until the contribution of the next term is very small (e.g., <10^-7). Values of for different values of and appear as shaded regions in Fig 1B&NDASH;D.

This folding and translating might seem to be a very involved process, but it is computationally much faster than realizing from a truncated normal distribution and computing the probability of X⁽ⁱ⁾ when is from such a distribution.

Multilocus Monte Carlo variance calculation:
We derive an expression for the variance of a product of J independent random variables, W_j, j = 1, ... J:

Denoting _Ne,j(Y) in (13) by W_j and taking the expectation gives the same result, verifying that the expression in (13) is unbiased for Var(^J_Ne(Y)).

	LITERATURE CITED

TOP ABSTRACT FORMULATION OF THE MODEL... SIMULATED AND REAL DATASETS RESULTS DISCUSSION APPENDIX LITERATURE CITED

BAUM, L. E., 1972 An inequality and associated maximization technique in statistical estimation for probabilistic functions of Markov processes, pp. 1–8 in Inequalities—III: Proceedings of the Third Symposium on Inequalities Held at the University of California, Los Angeles, September 1–9, 1969, edited by O. SHISHA. Academic Press, New York.

BAUM, L. E., T. PETRIE, G. SOULES, and N. WEISS, 1970  A maximization technique occurring in the statistical analysis of probabilistic functions on Markov chains. Ann. Math. Stat. 41:164-171.

BEGON, M., C. B. KRIMBAS, and M. LOUKAS, 1980  The genetics of Drosophila subobscura populations. XV. Effective size of a natural population estimated by three independent methods. Heredity 45:335-350.

CAVALLI-SFORZA, L. L. and A. W. F. EDWARDS, 1967  Phylogenetic analysis: models and estimation procedures. Evolution 21:550-570.

FRANKHAM, R., 1995  Inbreeding and extinction: a threshold effect. Conserv. Biol. 9:792-799.

GELMAN, A., 1996 Inference and monitoring convergence. pp. 131–143 in Markov Chain Monte Carlo in Practice, edited by W. R. GILKS, S. RICHARDSON and D. J. SPIEGELHALTER. Chapman & Hall, New York.

GELMAN, A., J. B. CARLIN, H. S. STERN and D. B. RUBIN, 1996 Bayesian Data Analysis. Chapman and Hall, New York.

HAMMERSLEY, J. M., and D. C. HANDSCOMB, 1964 Monte Carlo Methods. Methuen & Co. Ltd., London.

JORDE, P. E. and N. RYMAN, 1995  Temporal allele frequency change and estimation of effective size in populations with overlapping generations. Genetics 139:1077-1090[Abstract].

KRIMBAS, C. B. and S. TSAKAS, 1971  The genetics of Dacus oleae. V. Changes of esterase polymorphism in a natural population following insecticide control—selection or drift? Evolution 25:454-460.

NEI, M. and F. TAJIMA, 1981  Genetic drift and estimation of effective population size. Genetics 98:625-640[Abstract/Free Full Text].

NEYMAN, J. and E. L. SCOTT, 1948  Consistent estimates based on partially consistent observations. Econometrica 16:1-32.

POLLAK, E., 1983  A new method for estimating the effective population size from allele frequency changes. Genetics 104:531-548[Abstract/Free Full Text].

SOULÉ, M. E. (Editor), 1986 Conservation Biology: The Science of Scarcity and Diversity. Sinauer and Associates, Sunderland, MA.

WAPLES, R. S., 1989  A generalized approach for estimating effective population size from temporal changes in allele frequency. Genetics 121:379-391[Abstract/Free Full Text].

WILLIAMSON, E. G. and M. SLATKIN, 1999  Using maximum likelihood to estimate population size from temporal changes in allele frequencies. Genetics 152:755-761[Abstract/Free Full Text].

This article has been cited by other articles:

				J. P. Bollback, T. L. York, and R. Nielsen Estimation of 2Nes From Temporal Allele Frequency Data Genetics, May 1, 2008; 179(1): 497 - 502. [Abstract] [Full Text] [PDF]

				P. E. Jorde and N. Ryman Unbiased Estimator for Genetic Drift and Effective Population Size Genetics, October 1, 2007; 177(2): 927 - 935. [Abstract] [Full Text] [PDF]

				R. S. Waples and M. Yokota Temporal Estimates of Effective Population Size in Species With Overlapping Generations Genetics, January 1, 2007; 175(1): 219 - 233. [Abstract] [Full Text] [PDF]