## Abstract

A maximum-likelihood method for demographic inference is applied to data sets consisting of the frequency spectrum of unlinked single-nucleotide polymorphisms (SNPs). We use simulation analyses to explore the effect of sample size and number of polymorphic sites on both the power to reject the null hypothesis of constant population size and the properties of two- and three-dimensional maximum-likelihood estimators (MLEs). Large amounts of data are required to produce accurate demographic inferences, particularly for scenarios of recent growth. Properties of the MLEs are highly dependent upon the demographic scenario, as estimates improve with a more ancient time of growth onset and smaller degree of growth. Severe episodes of growth lead to an upward bias in the estimates of the current population size, and that bias increases with the magnitude of growth. One data set of African origin supports a model of mild, ancient growth, and another is compatible with both constant population size and a variety of growth scenarios, rejecting greater than fivefold growth beginning >36,000 years ago. Analysis of a data set of European origin indicates a bottlenecked population history, with an 85% population reduction occurring ∼30,000 years ago.

PATTERNS of genetic variation in contemporary populations can be used to make inferences about past population size changes. Ideally, likelihood methods using the full data would be applied to make such inferences. For the case of DNA sequence polymorphism and where no recombination occurs between the variable sites, methods are available for carrying out such inferences (Beerli and Felsenstein 2001; Kuhner *et al.* 1998; Nielsen 1999). With incomplete linkage between sites, such approaches are frequently computationally infeasible. An exception is the case in which only two chromosomes are sampled at each locus, where Marth* et al.* (2003) have shown that maximum-likelihood methods are feasible. These computational difficulties have led to the use of summary statistics such as Tajima's *D* (Tajima 1989) for making inferences about past demography. For example, Wall and Przeworski (2000) and Pluzhnikov* et al.* (2002) tested compatibility between observed values of Tajima's *D* and values observed in simulations under constant-size and alternative demographic scenarios. Weiss and Von Haeseler (1998) also focused on summaries of the data by implementing a likelihood approach based on mean pairwise differences and segregating sites for a model of complete linkage.

With free recombination between sites, the problem is greatly simplified. In this case, sites can be considered to be statistically independent of each other and the data are completely characterized by the number of polymorphic sites and the frequency spectrum. That is, we can represent the full data by **m** = (*m*_{0}, *m*_{1}, *m*_{2}, *m _{n}*

_{−1}), where

*m*

_{0}is the number of sites monomorphic in the sample, and, for

*i*> 0,

*m*is the number of polymorphic sites in which the derived allele is present

_{i}*i*times in the sample of

*n*chromosomes. We assume all polymorphic sites are biallelic. Then, is

*L*, the number of sites surveyed, and is the total number of segregating sites in the sample,

*S*. Also note that

*m*

_{0}=

*L*−

*S*. In this case of free recombination between sites, full-likelihood approaches are computationally undemanding. This case has been examined by Wooding and Rogers (2002), Polanski and Kimmel (2003), and Marth

*et al.*(2004) and is also the focus of our study. We examine the statistical properties of demographic inferences based on

**m**, using maximum likelihood and assuming sites are independent. By utilizing the entire frequency spectrum,

**m**, rather than a summary statistic such as Tajima's

*D*, this approach captures all available information in data sets consisting of unlinked polymorphic sites.

With linkage between sites, there is a statistical nonindependence between polymorphic sites, and thus **m** for a set of linked sites would contain less information than that for a set of unlinked sites. It follows that our results for unlinked sites give an idea of the best one can do using **m** or summary statistics, such as Tajima's *D*, which can be calculated from it. It is important to note that for linked sites **m** does not completely characterize the data, and that full likelihood, which incorporates information about linkage disequilibrium, for instance, might result in better inferences.

The models examined here consist of either exponential growth or an instantaneous decrease followed by exponential growth, which requires simultaneous estimation of either two or three parameters, respectively. We illustrate that, particularly for recent growth scenarios, data sets consisting of large numbers of segregating sites are required to produce good estimates based solely on frequency spectrum data. Our results provide a theoretical perspective on the feasibility of frequency spectrum-based parameter estimation with a modest amount of data, and we present methods to determine the approximate variance and covariance associated with such estimators under any demographic scenario of interest. The maximum-likelihood method is also applied to three human data sets. The first is an African data set consisting of the original data set of Frisse *et al.* (2001) as well as 40 additional locus pairs (A. Di Rienzo, unpublished data). The Seattle single-nucleotide polymorphisms (SNPs) data, consisting of both African-American and European data sets (http://pga.gs.washington.edu), is also examined. Each of these data sets consists of linked segregating sites within effectively unlinked loci, and a procedure is outlined by which one can extend this method to such data and construct the appropriate confidence regions associated with the estimators.

## MODEL AND METHODS

### Model:

The demographic model considered is that of a population of constant effective size *N*_{0} until time *T* when there was an instantaneous decrease to an intermediate population size (*N*_{int}) followed by exponential growth to the current population size (*N*_{rec}). As illustrated in Figure 1, this model involves four demographic parameters: *N*_{0}, *N*_{int}, *N*_{rec}, and *T*, where *T* is the time at which the instantaneous size change occurred. *T* is measured in units of 4*N*_{0} generations before the present. We assume the mutation rate per site, *u*, is small, so that the occurrence of more than one mutation occurring in the history of the sample at a single site can be ignored. We find it convenient to introduce the parameters *f*_{int} = *N*_{int}/*N*_{0} and *f*_{rec} = *N*_{rec}/*N*_{0} and specify the model by 4*N*_{0}*u*, *f*_{int}, *f*_{rec}, and *T*. This demographic model is flexible and can be generalized to the case of exponential growth with no bottleneck by setting *f*_{int} equal to one or to the case of a population reduction with no recovery by setting *f*_{rec} equal to *f*_{int}. We also assume that the population is unstructured (panmictic) and that the polymorphic sites are unlinked.

### Maximum-likelihood method:

The maximum-likelihood approach followed here is that of Wooding and Rogers (2002) and Polanski and Kimmel (2003). Our analyses require a population survey of variation at a set of *L* unlinked sites. For *L* unlinked sites, **m** is multinomially distributed, 1where *P*_{0} is the probability that a site is monomorphic in the sample, and, for *i* > 0, *P _{i}* is the probability that a site is polymorphic with

*i*copies of the derived allele. The

*P*'s are functions of the four parameters of the demographic model (θ

_{i}_{0},

*f*

_{int},

*f*

_{rec}, and

*T*) and the sample size

*n*.

To obtain the maximum-likelihood estimates of the parameters one maximizes the right-hand side of (1). We note, however, that we can write the probability of the data as 23where *p _{i} *(=

*P*/(1 −

_{i}*P*

_{0})) is the probability that a site is polymorphic with

*i*copies of the derived allele, conditional on the site being polymorphic in the sample.

*P*

_{0}and the

*p*'s can be written in terms of θ

_{i}_{0}and mean properties of sample gene trees. For example, for small θ

_{0},

*P*

_{0}can be expressed in terms of the mutation parameter and the mean total length of the gene genealogy, 4where τ(

*n*) is the mean total length of the gene tree of a sample of

*n*chromosomes measured in units of 4

*N*

_{0}generations (Hudson 1990). We define an

*i*-branch to be a branch of the gene tree such that a mutation that occurs on the branch results in

*i*copies of the mutation in the sample. The mean total length of

*i*-branches in units of 4

*N*

_{0}generations we denote by τ

*(*

_{i}*n*). Then,

*P*is ∼θ

_{i}_{0}τ

*(*

_{i}*n*), and 5

When time is measured in units of 4*N*_{0} generations, τ* _{i}*(

*n*) and τ(

*n*) are functions of

*f*

_{int},

*f*

_{rec}, and

*T*, but do not depend on

*N*

_{0}or θ

_{0}. Thus, to find the maximum-likelihood estimates of the four parameters, we can first find the maximum-likelihood estimates,

*f̂*

_{int},

*f̂*

_{rec}, and

*T̂*, by maximizing . The maximum-likelihood estimate of θ

_{0}, if desired, can then be obtained as θ̂

_{0}= (

*S/L*)/τ̂(

*n*), where τ̂(

*n*) is the mean gene tree length with

*f*

_{int},

*f*

_{rec}, and

*T*set equal to the maximum-likelihood estimates. In this article, we focus on estimation of the parameters

*f*

_{int},

*f*

_{rec}, and

*T*, which requires maximizing and does not require specifying

*L*or

*m*

_{0}. Equivalently, we can consider estimation of the parameters based on the probability of

**m**conditional on

*S*, which is 6which does not depend on

*L*or

*m*

_{0}.

To find the maximum-likelihood estimates of *f*_{int}, *f*_{rec}, and *T*, we estimate Prob(**m**|*S*) at a set of points on a rectangular grid of values in the three-dimensional space of *f*_{int}, *f*_{rec}, and *T* values. For each point in the grid, we estimate the τ* _{i}*(

*n*)'s by generating 100,000 replicate gene trees with simple one-site coalescent simulations. The τ

*(*

_{i}*n*)'s can also be obtained as described elsewhere (Griffiths and Tavaré 1998; Wooding and Rogers 2002; Polanski and Kimmel 2003), and this method can be generalized to any demographic model for which the relevant τ

*(*

_{i}*n*)'s can be calculated or estimated. From the estimated τ

*'s, the*

_{i}*p*'s are calculated, which in turn are used to calculate the product, . Since we have ignored the problem of estimating θ

_{i}_{0}, we do not require

*L*, and the results are all given conditional on specified numbers of polymorphic sites.

### Required data:

Our analyses require data in the form of unlinked polymorphic sites. Ascertainment bias is not considered in this article, so we assume that sites are randomly chosen with no prior knowledge of polymorphism and sequenced in each sampled chromosome. Frisse* et al.* (2001) sequenced ∼25 kb and found 120 segregating sites in an African Hausa sample of 30 chromosomes. If we consider genome-wide polymorphism levels to be similar to that data, then ∼104,000 sites would have to be sequenced in 30 chromosomes to assemble a data set consisting of 500 segregating sites. These 104,000 sites could be sequenced in small unlinked segments throughout the genome to obtain frequency spectrum data from unlinked polymorphic sites. We consider only biallelic polymorphic sites and assume that the ancestral/derived status of each allele is known. However, not having knowledge of the ancestral state makes only a minimal difference to our results (data not shown). We also assume that all sites are surveyed in a sample of *n* chromosomes (or *n*/2 diploid individuals), but more general sampling is easily accommodated. For example, one can separate a data set into a series of frequency spectra, each with a different *n*. A global likelihood may then be obtained for the entire data set by multiplying the result of Equation 6 for sets of sites with different sample sizes.

### Sample size comparison:

We compare the effect of sample size and number of unlinked polymorphic sites on both the power to reject the null hypothesis of constant population size and the quality of estimates of specific demographic parameters. The number of segregating sites, when indicated, is scaled on the basis of the average total branch length of a random gene genealogy (τ), which will vary according to sample size and demographic scenario. For example, suppose we wish to compare sample sizes of 50 and 100 chromosomes under a growth scenario where 40-fold expansion occurred beginning 10,000 years ago. In this case τ(50) (in units of 4*N*_{0} generations) for a sample size of 50 chromosomes is 4.93, τ(100) for a sample size of 100 chromosomes is 6.06, and τ(100)/τ(50) is 1.23. In words, the average total branch length of a random gene genealogy is 1.23 times greater for a sample size of 100 than for a sample size of 50. This indicates that, for every 500 segregating sites discovered in a sample size of 50, ∼615 segregating sites would be found in a sample size of 100 if the same number of sites were sequenced. Thus, when we compare sample size 50 to sample size 100, we compare *n* = 50, *S* = 500 to *n* = 100, *S* = 615. Normalizing the number of sites in this way serves to facilitate comparisons between analyses of different sample size because it takes into account the expected number of polymorphic sites in the samples of each size.

### Asymptotic properties:

To investigate whether asymptotic approximations of confidence regions and variances of estimators are applicable to data sets of modest size, we first determined whether 95% confidence regions obtained from the log-likelihood ratio have the expected coverage properties. Confidence regions include those points on the grid with a log-likelihood ratio <3 or 3.9 for simultaneous estimation of two or three parameters, respectively. We also compared the observed variances and covariances of the estimators with the approximate variances and covariances calculated by estimating the inverse of the information matrix, 7

We estimate the expected log-likelihood with 8

with the *P*_{i} estimated by coalescent simulation for a set of points on a narrow grid of *f*_{rec}, *f*_{int}, and *T* values around the true parameter values of *f*_{rec0}, *f*_{int0}, and *T*_{0}. The second partial derivatives relevant to the information matrix are then approximated from best-fit second-degree polynomial curves. The variance and covariance observed from simulation can then be compared to the appropriate terms of the inverse of *I _{ij}* to determine whether the asymptotic approximations apply to data sets of modest size.

## RESULTS

### Accuracy and precision of estimated τ_{i}'s:

_{i}

As described in model and methods, we estimate the relevant τ* _{i}*(

*n*)'s from 100,000 replicate gene trees generated by one-site coalescent simulations. We find that our τ

*(*

_{i}*n*)'s, estimated from simulation, are in very close agreement to τ

*(*

_{i}*n*)'s calculated numerically by the method of Polanski and Kimmel (2003). For the case of

*f*

_{rec}= 2.0,

*f*

_{int}= 0.15, and

*T*= 0.0375 for a sample size of 46, the maximum-likelihood parameters for the Seattle SNPs European data set, we find that our simulated τ

*(46)'s differ, at most, by 0.05% from the calculated τ*

_{i}*(46)'s. Additionally, we calculate the log-likelihood of the Seattle SNPs European data set using 10 independent τ*

_{i}*(46) estimates, each resulting from 100,000 replicate gene trees, and find that the likelihoods calculated from our simulated τ*

_{i}*(46)'s differ little between trials, ranging from −11987.278 to −11987.302. Since the log-likelihood ratio critical values relevant for our construction of confidence regions range from 3.86 to 9.1, such a negligible fluctuation would not affect our inferred acceptance regions.*

_{i}### Power curves:

Power analyses were conducted using a chi-square test with *n −* 2 d.f. for a sample size of *n* chromosomes to determine the degree of growth that would be required to reject the null hypothesis of constant population size using only frequency spectrum information. For smaller numbers of segregating sites, the degrees of freedom may vary slightly, as frequency categories with expected site counts of less than five are collapsed. The expected frequency spectrum under the null hypothesis is calculated from 9

(Ewens 1979) by multiplying each *p _{i}* by the number of segregating sites. The observed frequency spectrum is obtained by estimating the

*p*'s from 100,000 replicates for each combination of

_{i}*f*

_{int},

*f*

_{rec}, and

*T*values and then multinomially sampling from these simulated

*p*'s. For each sample size, the number of segregating sites at each

_{i}*f*

_{rec}value is scaled on the basis of 500 polymorphic sites in a sample size of 20, as described in model and methods.

#### Recent growth beginning 10,000 years ago:

Figure 2a shows power curves for the scenario of recent growth beginning at *T* = 0.0125 (which, for humans, would correspond to 10,000 years ago on the basis of a generation time of 20 years and *N*_{0} of 10,000, roughly corresponding to the advent of agricultural society). We consider sample sizes of 10, 20, 50, 100, and 250 chromosomes with 500 polymorphic sizes (scaled as described above), which, for a sample size range of 10–250, corresponds to a range of 398–859 sites at constant population size (*f*_{rec} = 1) to 390–1137 sites at the most extreme growth scenario considered (*f*_{rec} = 250). With sample sizes of ≤20, the power to reject the null hypothesis of constant population size never exceeds 0.15, even with 500-fold growth. With a sample size of 50, power reaches ∼0.5 with 50-fold growth, but barely rises above 0.6 at the largest magnitude of growth considered. As sample size reaches 100, one can reliably detect 20-fold growth, and a sample size of 250 allows for a power near 1 to reject the null hypothesis with only 5-fold growth (Figure 2a). It is clear that recent rapid growth can be reliably detected with frequency spectrum data only with fairly large samples (>100 chromosomes), and the most modest growth scenarios may be detected only with samples consisting of at least 250 chromosomes when data sets consist of only 500 (scaled) polymorphic sites.

#### More ancient growth onset:

Power to reject the constant size hypothesis is also dependent upon the time that exponential growth begins, as illustrated in Figure 2b. While power is minimal for small sample sizes with growth beginning 10,000 years ago, power increases dramatically with more ancient growth. For example, while a sample size of 20 with 500 polymorphic sites yields virtually no power to detect any magnitude of growth beginning 10,000 years ago, if growth instead began 50,000 years ago, a sample size of 20 with the same number of sites would be sufficient to reliably detect 10-fold growth.

### Asymptotic properties:

We evaluate the distribution of our maximum-likelihood estimates to determine whether asymptotic theory provides an adequate approximation of the 95% confidence regions and variance associated with our parameter estimates. Table 1 illustrates the proportion of maximum-likelihood estimates for which the true value of the parameters lies outside the asymptotic 95% confidence region. Our simulations indicate that for large amounts of data, asymptotic theory does provide a good approximation of the 95% confidence region for the demographic scenario examined. For smaller amounts of data, the asymptotic approximation appears to be conservative, with the true parameter values lying within the 95% confidence region in ∼97–98% of the runs. We also examine a specific case corresponding to the Hausa data set, which consists of 597 sites in a sample size of 30. For a data set of this size [simulated from the Hausa maximum-likelihood estimate (MLE)], we find that asymptotic approximation is especially conservative, rejecting the true parameter values in only 2.04% of the 5000 simulated data sets.

We also determine the variance and covariance of our maximum-likelihood estimators in two dimensions by both asymptotic theory and simulation (Table 2). In this analysis, we assume that *N*_{rec} is known and is fivefold >*N*_{0} (*f*_{rec} = 5), while *f*_{int} and *T* are jointly estimated. For this demographic scenario, asymptotic theory provides a good approximation for the simulated variance and covariance only when the data set consists of a large number of segregating sites.

### Quality of estimators:

We evaluate the quality of our maximum-likelihood estimators by examining the distribution of the estimates under both two-dimensional and three-dimensional models.

#### Two-dimensional estimators:

A recent growth scenario was examined in which the population size was constant until exponential growth occurred beginning 10,000 years ago. Data sets were simulated with *f*_{rec} ranging from 10- to 320-fold growth, *f*_{int} fixed at 1, and *T* = 0.0125 (10,000 years ago based on a generation time of 20 years and *N*_{0} of 10,000). For this scenario of recent growth, a sample size of at least 250 chromosomes with ∼16,000 segregating sites is required for 90% of the distribution of *f̂*_{rec} to lie within a factor of four of the true *f*_{rec} value for all magnitudes of growth examined. This is illustrated in Figure 3, which compares the *f̂*_{rec} distribution under this recent growth scenario for sample sizes of 50 and 250 for 20-fold growth. As the magnitude of growth increases, *f̂*_{rec} becomes biased more severely upward. This result is similar to that obtained for growth beginning at *T* = 0.0625 (Table 4).

The estimates of *T*, however, are not subject to the upward bias seen in *f̂*_{rec}. Instead, estimates of *T* are improved as the magnitude of growth increases (Table 3). Sample size also has a dramatic effect on the distribution of *T̂*, as illustrated in Figure 4, which compares the *T̂* distribution for sample sizes of 50 and 250. With 5000 sites in a sample size of 50 chromosomes, 95% of *T̂* estimates lie within a factor of 3 of the true *T* value for 10-fold growth, with 95% of the distribution lying within a factor of 1.5 for 320-fold growth, the most severe growth scenario examined.

We explored another growth scenario in which the onset of growth was more ancient, beginning 50,000 years ago. In this case, the quality of the estimators improved, and 90% of the *f̂*_{rec} distribution was within a factor of four of the true *f*_{rec} value for a sample size of 50 with 5000 sites (as opposed to a sample size of 250 and ∼16,000 sites under the more recent growth scenario). Figure 3 reveals the improvement in the *f̂*_{rec} estimator with the more ancient time of growth onset. As in the recent growth scenario, increasing the degree of growth both increased the bias and widened the quantiles of the *f̂*_{rec} distribution (Table 4). The *T* estimates under the more ancient growth scenario were also improved over the analogous recent growth estimates, with 95% of the *T̂* distribution within a factor of two of the true *T* value for all magnitudes of growth examined with data sets as small as a sample size of 50 with 1000 sites.

#### Three-dimensional estimators:

We consider a three-dimensional model of a constant-sized population that experienced an instantaneous decrease to 0.05 times its initial size 100,000 years in the past, followed by exponential growth until the present to a final size of 5 times the initial population size (*f*_{rec} = 5; *f*_{int} = 0.05; *T* = 0.125). All three parameters were estimated for 5000 simulated data sets. Under this model, 90% of the distribution of each of the three estimators falls within a factor of 4 of the respective true values with data sets as small as 500 sites in a sample size of 30 (Table 5). If a large data set consisting of 10,000 polymorphic sites in 50 chromosomes were available, 95% of the estimates of all three parameters would lie within a factor of 1.5 of the true parameter values. As would be expected, estimates of any of the three parameters are improved by fixing one of the parameters at its true value (data not shown), indicating that incorporation of prior knowledge of one of the parameters would be beneficial.

### Applications:

We apply the maximum-likelihood method to data obtained from an African Hausa population (A. Di Rienzo, unpublished data) as well as to the African-American and European (CEPH) samples of the Seattle SNPs data set (http://pga.gs.washington.edu).

#### Hausa data:

The Hausa data set consisted of the data of Frisse* et al*. (2001) in conjunction with additional unlinked locus pairs (A. Di Rienzo, unpublished data), which resulted in a data set consisting of 30 chromosomes and 597 polymorphic sites in an African sample, the Hausa of Cameroon. The sites in this data set include linked polymorphic sites within 50 effectively unlinked loci, but in the maximum-likelihood analysis we treat each site as though it provides independent information. As seen in Table 6, *f̂*_{rec} = 3.1, *f̂*_{int} = 1, and *T̂* = 6.1 for this data set.

We perform a χ^{2} goodness-of-fit test on the Hausa data set to determine whether the maximum-likelihood parameters can be accepted as an explanation of the Hausa data. However, this test assumes that each site is independent, which is not the case for this data set. Because the linkage between sites will affect the 95% critical value of the χ^{2} test statistic, we determine the critical value of the test statistic distribution for this data set by coalescent simulation with recombination (Hudson 1983, 2002). We simulate 5000 data sets, each consisting of 30 chromosomes and 50 unlinked loci. The input parameters for the simulation included Watterson's estimate of θ (estimated to be 0.0012/bp), the recombination rate (estimated to be 5.99 × 10^{−4}/bp), the average locus length (10,286 bp), and a gene-conversion to crossing-over ratio of 2. Polymorphic sites within the middle 8000 bp were ignored to mimic the locus pair data (Frisse *et al.* 2001). Because the ancestral/derived status of each allele was not considered, each simulated frequency spectrum was folded at frequency 0.5 prior to performing the χ^{2} goodness-of-fit test. On the basis of these simulations, we found the 95% critical value of the χ^{2} test statistic to be 39.39, as opposed to a critical value of 23.68 (14 d.f.) if all sites were independent. The χ^{2} goodness-of-fit test statistic for the Hausa data set under its maximum-likelihood estimate of *f̂*_{rec} = 3.1, *f̂*_{int} = 1, and *T̂* = 6.1 is 26.30 (*P* = 0.304), indicating that this scenario cannot be rejected at the 0.05 significance level. Note, however, that this demographic scenario would have been rejected without properly accounting for the linkage within the data set. We also considered the equilibrium model of constant population size for this data set and obtained a χ^{2} goodness-of-fit test statistic of 29.24 (*P* = 0.204). On the basis of an analysis of 10 locus pairs (a subset of the 50 locus pairs examined here), Frisse* et al*. (2001) also concluded that the Hausa data are consistent with the equilibrium model.

#### Seattle SNPs:

We also examined both the African-American and the European samples of the Seattle SNPs data, which are located at the University of Washington-Fred Hutchinson Cancer Research Center (UW-FHCRC) Variation Discovery Resource (http://pga.gs.washington.edu). These data sets consist of 12,587 and 7712 total SNPs across 138 loci in the African-American and European samples, respectively. We considered only those SNPs that were sequenced in the entire panel of 48 (African-American) or 46 (European) chromosomes to facilitate evaluation of confidence regions and goodness-of-fit by simulation. Additional analyses incorporating more of the SNPs are described in the discussion. The frequency spectrum of nonsynonymous SNPs has been shown to differ from that of synonymous SNPs (Cargill *et al.* 1999; Fay *et al.* 2001; Wooding and Rogers 2002), so we also removed all SNPs that result in an amino-acid coding change to minimize the inclusion of those SNPs subject to nonneutral evolutionary processes. This resulted in a final data set of 5892 SNPs for the African-American data set and 4211 SNPs for the European data set. We applied our maximum-likelihood method to these data sets, treating all SNPs as unlinked, and found that the three-dimensional maximum-likelihood estimates of *f*_{rec}, *f*_{int}, and *T* are *f̂*_{rec} = 1.9, *f̂*_{int} = 1, and *T̂* = 0.27 for the African-American data set and *f̂*_{rec} = 2, *f̂*_{int}= 0.15, and *T̂* = 0.0375 for the European data set (Table 6). These estimates suggest a scenario of very slow growth over a long period of time with no bottleneck for the African-Americans and a fairly recent population bottleneck with ∼13-fold recovery for the Europeans.

To determine whether the demographic model we consider is compatible with the Seattle SNPs data, we simulate the distribution of the goodness-of-fit test statistic for this data set as described for the Hausa data set. For these simulations, each data set consisted of 48 chromosomes and 138 loci. The input parameters were that of the Hausa data set, except substituting the average length of a Seattle SNPs locus. In this case, each locus was simulated, fixing the number of segregating sites to be the average number of segregating sites per locus in the African-American or European Seattle SNPs data set, so each simulated data set contained the same total number of segregating sites as our observed Seattle SNPs African-American or European data set. Because the Seattle SNPs data sets do not specify the ancestral/derived status of each allele, each simulated frequency spectrum is again folded. The 95% critical values of the distribution were found to be 48.68 (African) and 137.36 (European) as opposed to 35.17 (23 d.f.) and 33.92 (22 d.f.) if all sites were unlinked.

Using these simulated critical values, a χ^{2} goodness-of-fit test indicates that the maximum-likelihood parameters produce an expected frequency spectrum that is not significantly different from the observed Seattle SNPs European data (χ^{2} = 122.98; *P* = 0.2015). Therefore, we can accept our simple bottleneck model as a reasonable explanation for this data set. The same test indicates that a constant population size model is not compatible with the European data (χ^{2} = 207.286; *P* < 1 × 10^{−4}).

However, the χ^{2} goodness-of-fit test on the African-American data set reveals that the frequency spectrum predicted by the maximum-likelihood estimates of *f*_{rec}, *f*_{int}, and *T* is significantly different from the empirical Seattle SNPs African-American frequency spectrum (χ^{2} = 86.64; *P* = 2 × 10^{−4}); therefore, our simple demographic model cannot be accepted as a complete explanation of the African-American data set, although the fit is better than that predicted by the constant population size model (χ^{2} = 268.66; *P* ≪ 1 × 10^{−4}). Figure 5 provides a visual comparison of the observed Seattle SNPs frequency spectrum to the frequency spectra predicted by both the maximum-likelihood parameters and constant population size parameters, indicating that the lack of fit of the maximum-likelihood parameters does not seem to be confined to any particular nonsingleton frequency class. However, our demographic model with the maximum-likelihood parameters appears to provide a better fit to the data than the equilibrium model, particularly in the singleton class. In addition, we note that a χ^{2} goodness-of-fit test shows that the Hausa data are compatible with *f*_{rec} = 1.9, *f*_{int} = 1, and *T* = 0.27, the estimates obtained from the Seattle SNPs African-American data set (χ^{2} = 33.46; *P* = 0.113).

Because the sites in the Hausa and Seattle SNPs data sets are not entirely unlinked, asymptotic approximation of confidence intervals is not appropriate. We simulated 10,000 data sets as described above for both the Hausa and Seattle SNPs data sets, using their respective maximum-likelihood estimates for input parameters, and applied the maximum-likelihood method to the folded frequency spectrum of each simulated data set. For the Hausa and Seattle SNPs African-American data sets, we estimated both *f*_{rec} and *T*, fixing *f*_{int} at 1, which was the maximum-likelihood estimate for both data sets. All three parameters were estimated for the data sets simulated from the Seattle SNPs European parameters. The ratio of the log-likelihood at the maximum-likelihood parameters to the log-likelihood at the parameters from which the data set was simulated could then be calculated. From the log-likelihood ratio distribution, we determined the 95% critical value to be 3.86 for the Hausa data set and 4.85 for the Seattle SNPs African-American data set as compared to the asymptotic critical value of 3.0 for two-dimensional maximum-likelihood estimates. The 95% critical value of 9.1 was found for the European data set as compared to the asymptotic critical value of 3.9 for three-dimension estimation. Using the critical values from simulation, we can easily reject the constant-size population model for the Seattle SNPs African-American and European data sets since the log likelihood ratios are 96 and 35, respectively (Table 6).

Figure 6a provides a visual representation of the 95, 99, and 99.9% confidence regions of the Hausa data set obtained by including all parameter values for which the log-likelihood ratio is ≤3.86, 6.38, and 10.08, respectively. Likewise, Figure 7 illustrates the analogous confidence regions for the Seattle SNPs African-American (Figure 7a) and European (Figure 7b) data sets.

## DISCUSSION

Our power analyses on models with exponential growth beginning 10,000 years ago illustrate that the frequency spectrum does not provide sufficient information to reject the null hypothesis of constant population size when either small sample sizes (<50 chromosomes) or small numbers of unlinked sites (<1000) are available. This result should serve as a cautionary note to researchers interested in demographic models involving expansion as recent as 10,000 years ago. Prior knowledge of the model of interest should also be considered when determining whether the frequency spectrum retains the requisite information for demographic inference, as the power to detect departures from constant size increases with both the extent of growth (Figure 2, a and b) and the time since the onset of growth (Figure 2b).

Application of the maximum-likelihood method on recent growth scenarios reveals that data sets consisting of at least 250 chromosomes with at least 10,000 scaled segregating sites (15,838–17,140 segregating sites depending upon the true *f*_{rec} value) are required for the *f̂*_{rec} distribution to have 95% critical values that fall within a factor of four of the true *f*_{rec} value if growth began as recently as 10,000 years ago (*T* = 0.0125). Unless large sample sizes and many unlinked sites are surveyed, the frequency spectrum alone provides little information about the magnitude of growth that has occurred relatively recently. As *f*_{rec} increases, the frequency spectrum becomes more distinct from what would be expected under a constant size scenario. However, with increasingly extreme recent growth, the frequency spectrum becomes less distinguishable from that of other severe growth scenarios, and it becomes more difficult to estimate the *f*_{rec} parameter with frequency spectrum information alone.

While it is difficult to accurately estimate *f*_{rec} for scenarios of recent growth, *T* can be estimated with more modest amounts of data. The distribution of *T̂* has 95% critical values that fall within a factor of four of the true *T* value for sample sizes as small as 50 chromosomes and 5000 sites, as compared to a sample size of 250 and 15,838–17,140 sites required to estimate *f*_{rec} to the same accuracy. Additionally, *T̂* is not subject to the upward bias seen in *f̂*_{rec} and estimates of *T* actually improve with increasing *f*_{rec}. Estimates of both *f*_{rec} and *T* improve as the onset of growth becomes more ancient. This observation is consistent with our observation that power to reject the null hypothesis of constant population size with frequency spectrum data increases with scenarios of more ancient growth (Figure 2b).

Simultaneous estimation of all three parameters results in estimator distributions where 90% of the estimates lie within a factor of four of the true parameter values with data sets as small as 500 segregating sites in a sample size of 30 for a model where the population decrease and subsequent expansion began 100,000 years ago and the present population is only five times the initial population size. These estimates benefit from both a more ancient time of growth onset and a modest magnitude of growth that is not subject to the upward bias seen in more severe growth scenarios. The ability of the frequency spectrum alone to elucidate the time and magnitude of population size change events is, therefore, greatly dependent upon the underlying demographic model. While ancient demographic events may be inferred relatively accurately from contemporary frequency spectrum patterns, more recent and severe episodes of growth are problematic for this method and require exceedingly large amounts of unlinked data. For these recent growth scenarios, it is possible that more informative estimates could be obtained by using a method that uses linked polymorphic sites and considers additional aspects of the data such as levels of linkage disequilibrium.

Evaluation of the asymptotic properties of our maximum-likelihood estimators indicates that asymptotic theory provides a reasonable approximation of the confidence intervals associated with the estimators. As we illustrate with the Hausa and Seattle SNPs data sets, it is also possible to construct these confidence intervals around a maximum-likelihood estimate through simulation. By simulating data sets that closely match the properties of the observed data set, one can estimate the critical value of this log-likelihood ratio distribution and construct corresponding confidence regions. This procedure is particularly relevant when asymptotic approximation is not appropriate, such as when the segregating sites in a data set are not unlinked.

We apply the maximum-likelihood method to both the African Hausa data set and the African-American and European samples of the Seattle SNPs data set. In both the Hausa and the Seattle SNPs data sets, the segregating sites are not entirely unlinked, but the maximum-likelihood analysis treats them as though each site provides independent information. However, we illustrate how one may use simulation to construct confidence regions and use goodness-of-fit tests that take into account the linkage between sites.

In the Seattle SNPs African-American data set, the simulated 95% confidence interval clearly allows for rejection of the constant population size model, since the log-likelihood of observing the data is almost 100 units less with the constant size parameters than with the estimated parameters. The maximum-likelihood estimates of *f̂*_{rec} = 1.9, *f̂*_{int} = 1, and *T̂* = 0.27 based on the Seattle SNPs African-American data correspond to a slow, ancient growth scenario where growth began >200,000 years ago to a present size of approximately two times the initial population size.

The simulated 95% confidence region around the Seattle SNPs African-American maximum-likelihood estimates, shown in Figure 7a, includes only a very narrow range of *f*_{rec} values within 1.6–2.5. However, the confidence region includes a wide range of *T* values ranging from as recent as 80,000 years ago (*T* = 0.1) to the most ancient time examined, 800,000 years ago (*T* = 1), assuming a generation time of 20 years and an *N*_{0} of 10,000. Even with the most recent compatible *T* value, it is not surprising that this data set allows for rejection of the constant size hypothesis with an estimate of only twofold growth. Our power analyses show that a data set consisting of 50 chromosomes has a power of 0.9 to reject the constant size hypothesis with only 1000 unlinked sites for twofold growth beginning 100,000 years ago. While the Seattle SNPs data set does not consist of entirely unlinked sites, our analysis included >5000 polymorphic sites across 138 loci, which should allow for comparable power.

Despite the compact confidence region (Figure 7a), visually reasonable fit of the frequency spectrum under the maximum-likelihood parameters to the observed data (Figure 5), and compatibility with the Hausa data set, a χ^{2} goodness-of-fit test indicates that our simple three-dimensional demographic model with the maximum-likelihood estimates obtained from the Seattle SNPs African-American data set is incompatible with the Seattle SNPs data (*P* = 2 × 10^{−4}), even when linkage is taken into account. The visual comparison between the observed and maximum-likelihood frequency spectra in Figure 5 seems to indicate that the number of singletons expected under the maximum-likelihood parameters is very close to the observed value, and therefore the incompatibility must be due to some combination of the other frequency categories. The loci in the Seattle SNPs data set were chosen because of their role in inflammatory pathways and may reflect the action of evolutionary forces other than population size changes. Although we removed those SNPs that result in amino-acid coding changes, which are more apt to be subject to natural selection, it is still probable that this data set includes SNPs that are mildly deleterious and may influence the frequency spectrum toward greater numbers of low-frequency variants and mimic evidence of growth.

The African-American population sampled for the Seattle SNPs data set may also be subject to population structure and admixture, which could affect the frequency spectrum and confound our inference about demographic history (Ptak and Przeworski 2002). To determine the effect of European admixture on maximum-likelihood estimates obtained from an African data set, we randomly combined 6 Italian chromosomes (A. Di Rienzo, unpublished data) with the 30 Hausa chromosomes at each of the 50 locus pairs of the Hausa data set, which resulted in a total data set of 657 polymorphic sites in 36 chromosomes. This represents ∼17% European admixture, which is consistent with admixture estimates obtained from African-American populations (Parra *et al.* 1998). Admixture of this proportion had virtually no effect on the maximum-likelihood estimates or confidence intervals based on the original Hausa data set (data not shown). Regardless, that does not eliminate the possibility that either the true population structure could involve admixture in different proportions or admixture in a larger data set such as the Seattle SNPs would produce a more prominent effect. The frequency spectrum of the Seattle SNPs data set is certainly not consistent with an equilibrium model of constant population size, although the degree of growth predicted is less than that of some previous reports based on African populations (Aris-Brosou and Excoffier 1996; Pritchard *et al.* 1999). However, our estimate of twofold growth beginning as recently as 80,000 years ago is consistent with a recent study based on the frequency spectrum in an African-American population (Marth *et al.* 2004).

The maximum-likelihood parameters estimated from this data set are consistent with the Hausa data set, which contains noncoding loci that are less likely to be subject to confounding factors such as selection. However, this analysis does not preclude population structure within Africa as a potential influence on the maximum-likelihood estimates of the Hausa data set. The maximum-likelihood estimates from the Hausa data set (*f̂*_{rec} = 3.1, *f̂*_{int} = 1, and *T̂* = 6.1) correspond to a scenario of slow, ancient threefold growth, beginning several million years ago. However, the confidence region associated with this data set (Figure 6a) is consistent with a wide range of growth scenarios, including both the demographic history estimated from the Seattle SNPs data set (*f̂*_{rec} = 1.9, *f̂*_{int} = 1, and *T̂* = 0.27) and *f*_{rec} = 1, which corresponds to constant population size. Additionally, Figure 6b provides a close-up view of the acceptance region of Figure 6a, considering only more recent values of *T* where the onset of growth occurs no more than 80,000 years ago. If we focus on these *T* values, it is clear that our confidence regions on this data set do not exclude scenarios of 20-fold or more growth, provided that the time of onset is correspondingly more recent. For example, if we believe that the Hausa population has undergone growth >5-fold, then our analysis indicates that the growth must have begun no earlier than *T* = 0.045 (36,000 years ago if *N*_{0} is 10,000 assuming a generation time of 20 years). Growth of that magnitude or larger is rejected at the 1% level (Figure 6b) for values of *T* > 0.045 and <3 (∼36,000–2.4 million years ago).

Our analysis of the Seattle SNPs European data set reveals an estimated demographic history of *f̂*_{rec} = 2.0, *f̂*_{int} = 0.15, and *T̂* = 0.0375, which corresponds to an 85% reduction in population size at *T* = 0.0375 (30,000 years ago assuming *N*_{0} = 10,000 and a 20-year generation time) and then ∼13-fold exponential growth to a current population size of twice the ancestral size. In constructing our data set, we exclude all SNPs that are not successfully typed in every chromosome to facilitate construction of appropriate confidence regions and estimation of χ^{2} critical values through simulation. However, we note that if all SNPs that were typed in at least half of the sampled chromosomes (7410 SNPs) were included in this analysis, we get only a slightly different estimate (*f̂*_{rec} = 1.25, *f̂*_{int} = 0.2, and *T̂* = 0.05) that differs by <2 log-likelihood units from our maximum-likelihood estimate based on the filtered data (−20,668.96 *vs.* −20,670.93 when all SNPs are included). Since the 95% confidence region includes all parameter values within 9.1 log-likelihood units from the maximum, it is not likely that this filtering of the data would result in a significant shift in our acceptance region.

A χ^{2} goodness-of-fit test indicates that frequency spectrum produced by the estimated parameters (*f̂*_{rec} = 2.0, *f̂*_{int} = 0.15, and *T̂* = 0.0375) is a reasonable match to the observed Seattle SNPs European frequency spectrum with a *P*-value of 0.2015. The constant size model is both rejected by the goodness-of-fit test and excluded by the simulated likelihood-ratio confidence region for this data set (Table 6, Figure 7b). These results implicate a bottlenecked history for this European data set, which is consistent with previous studies (Marth *et al.* 2003, 2004) and the “Out of Africa” model for human population history (Harpending *et al.* 1998). Since the Seattle SNPs European data set is composed of the same coding loci as the Seattle SNPs African-American data set, it seems reasonable that the lack of agreement between the frequency spectrum predicted by the maximum-likelihood parameters and the observed African-American frequency spectrum is more likely to be due to population structure than to the presence of slightly deleterious variants in the data set. Of course, the good fit of the European maximum-likelihood parameters does not preclude the possibility of population structure or selection within the European data set as well.

As we have indicated earlier, there are a number of confounding factors to consider when attempting to infer demographic history based on frequency spectrum information, including population structure (past or present) and selection. An additional complication that is not considered by these analyses is ascertainment bias and genotyping error. It has been shown that ascertainment bias can lead to large errors in maximum-likelihood-based demographic inference (Kuhner *et al.* 2000; Wakeley *et al.* 2001). Polanski and Kimmel (2003) have also shown that exclusion or misclassification of low-frequency SNPs can result in estimated growth rates that are significantly lower than the true value. Note, however, that the sites represented in the Di Rienzo and Seattle SNPs data sets were chosen without prior indication of polymorphism status. Therefore, the analyses on these data sets would not be influenced by ascertainment bias due to using a discovery sample for SNP identification. However, we cannot exclude the possibility that genotyping errors have biased our inferences.

### Conclusions:

Analysis of this maximum-likelihood method indicates that demographic inferences can be drawn from frequency spectrum data when sufficient amounts of data are available. Asymptotic theory or simulation can be used to determine the variance and covariance associated with these estimators to determine whether the maximum-likelihood estimates would be meaningful for a particular demographic model and amount of data that may be available. However, our results show that very large amounts of data may be required to obtain practical confidence regions, particularly in models involving recent growth. For growth beginning as recently as 10,000 years ago, the power to reject the hypothesis of constant population size is very low with sample sizes of <20 chromosomes. To make accurate inferences under this type of recent-growth model using the frequency spectrum alone, both large samples (>100 chromosomes) and a large number of unlinked sites (>5000 sites) are required, although estimators improve as the time of onset of growth becomes more ancient. In scenarios of extreme growth, there is also a severe bias in *f̂*_{rec}, even with large amounts of data. However, *T* can be estimated with more modest amounts of data, and *T̂* is not subject to the bias seen in *f̂*_{rec}, indicating that one may obtain reasonable estimates of the time of population size-change events, even if the magnitude is biased. This maximum-likelihood method incorporates all available information contained in unlinked polymorphic sites, and parameter estimation methods based on summaries of the frequency spectrum require even larger amounts of data to be equally as informative. Therefore, for scenarios where the entire frequency spectrum of modest data sets does not provide an adequate amount of information, it may be necessary to incorporate additional aspects of linked data to improve estimates of demographic parameters.

Application of the maximum-likelihood method to three human data sets implicates differing demographic histories for African *vs.* European data sets. The African Hausa data set is compatible with a wide range of growth scenarios, ranging from slow, ancient growth to some scenarios of very recent, rapid growth. However, we can reject episodes of greater than 5-fold growth beginning >36,000 and <2.4 million years ago on the basis of this data set. The Seattle SNPs African-American data set also supports a model of growth, although a goodness-of-fit test indicates that the best-fit model of ancient, slow growth is not sufficient to explain the observed frequency spectrum. Maximum-likelihood analysis of the Seattle SNPs European data set reveals that the best-fit model is one of a population bottleneck occurring ∼30,000 years ago, reducing the population to 15% of the ancestral size, followed by 13-fold growth to a current population size that is twice the ancestral size.

## Acknowledgments

We thank A. Di Rienzo for access to the Hausa data set prior to publication. Construction of the Hausa data set was supported by National Institutes of Health (NIH) grant HG-02098. A. Adams was supported by NIH/National Institute of General Medical Sciences grant 5 T32 GM07197.

## Footnotes

Communicating editor: M. Nordborg

- Received April 17, 2004.
- Accepted July 13, 2004.

- Genetics Society of America