Approximate maximum likelihood estimation for population genetic inference

2.1 Maximization by stochastic approximation

Here, we summarize stochastic approximation in its classic form. In Subsection 2.2., we discuss how stochastic approximation can be adopted to obtain maximum likelihood estimates.

Let Y ∈ ℝ be a random variable depending on Θ∈ℝp. The function L(Θ)=E(Y∣Θ) is to be maximized in Θ, but L(Θ) as well as the gradient ∇L(Θ) are unknown. If realizations y(Θ) of Y∣Θ can be obtained for any value of Θ, argmaxΘ∈ℝpL(Θ) can be approximated by stochastic approximation methods (see Spall, 2003, sections 6 and 7, for an overview).

Similar to gradient algorithms in a deterministic setting, the stochastic approximation algorithms are based on the recursion

starting with some Θ0∈ℝp.

The unknown gradient can be substituted by an approximation based on observed finite differences. For Θ∈ℝ this algorithm was first described in Kiefer and Wolfowitz (1952), followed by a multivariate version in Blum (1954), where the gradient is approximated by finite differences in each dimension. In iteration n, each element l = 1, …, p of the gradient is approximated by

where e_l is the lth unit vector of length p, cn∈ℝ+ is a decreasing sequence and y(Θn+cnel) and y(Θn−cnel) are independent realizations of Y∣Θn+cnel and Y∣Θn−cnel, respectively.

Thus, for each iteration of this algorithm, 2p observations of Y are needed. A computationally more efficient method was introduced by Spall (1992): The finite differences approximation of the slope is only obtained along one direction that is randomly chosen among the vertices of the unit hypercube, so only two observations per iteration are necessary. More specifically, in iteration n, a random vector with elements

for l = 1, …, p is generated. Then, the gradient is approximated by

Spall showed that for a given number of simulations this algorithm reaches a smaller asymptotic mean squared error in a large class of problems than the original version. Phrased differently, fewer realizations of Y are usually necessary to reach the same level of accuracy. See Spall, 2003, Section 7, for a more detailed overview.

2.2 Approximate maximum likelihood (AML) algorithms

Suppose, data D_obs are observed under model ℳ with unknown parameter vector Θ∈ℝp. Let L(Θ;Dobs)=p(Dobs∣Θ) denote the likelihood of Θ. For complex models often there is no closed form expression for the likelihood, and for high dimensional data sets, the likelihood can be difficult to estimate. As in ABC, we therefore consider L(Θ;sobs)=p(sobs∣Θ), an approximation to the original likelihood that uses a d-dimensional vector of summary statistics s_obs instead of the original data D_obs. The more informative s_obs is about D_obs, the more accurate is this approximation. In the rare cases where s_obs can be chosen as a sufficient statistic for Θ, L(Θ;Dobs)∝L(Θ;sobs).

An estimate L^(Θ;sobs) of L(Θ;sobs) can be obtained from simulations under the model ℳ⇐Θ⇒ using kernel density estimation. With L^(Θ;sobs) as our objective function, we approximate the maximum likelihood estimate Θ^ML of Θ using the following algorithm:

The approximate maximum likelihood estimate is obtained at the final step of this algorithm, i.e. Θ^AML-FD,N:=ΘN. Notice however, that more sophisticated versions of the algorithm involve averaging over the last s steps.

Alternatively, an algorithm based on Spall’s simultaneous perturbations method can be defined as follows:

Here, the approximate maximum likelihood estimate is denoted Θ^AML-SP,N:=ΘN. In more general statements, we denote both approximate maximum likelihood estimates by Θ^AML.

Instead of the likelihood, one might want to equivalently maximize the log-likelihood, which can be preferable from a numerical point of view. This can be done both with the AML-FD and the AML-SP algorithm by replacing L^(Θ+;sobs) and L^(Θ−;sobs) by their logarithms in step 1d.

In a Bayesian setting with a prior distribution π(Θ), the algorithms can be modified to approximate the maximum of the posterior distribution, Θ^MAP, by multiplying L^(Θ−,sobs) and L^(Θ+,sobs) with π(Θ−) and π(Θ+), respectively.

2.3 Parametric bootstrap

Confidence intervals and estimates of the bias and standard error of Θ^AML can be obtained by parametric bootstrap: B bootstrap datasets are simulated from the model ℳ⇐Θ^AML⇒ and the AML algorithm is run on each dataset to obtain the bootstrap estimates Θ^AML,1∗,…,Θ^AML,B∗. These estimates reflect both the error of the maximum likelihood estimator as well as the approximation error of the algorithm.

The bias can be estimated by

where Θ¯∗=(∑i=1BΘ^i∗)/B. The corrected estimator is

The standard error of Θ^AML can be estimated by the standard deviation of Θ^1∗,…,Θ^B∗,

We compute separate bootstrap confidence intervals for each parameter that are based on the assumption that the distribution of Θ^AML−Θ can be approximated sufficiently well by the distribution of Θ^AML∗−Θ^AML, where Θ^AML∗ is the bootstrap estimator. Then, a two-sided (1 − α)-confidence interval is defined as

where q(β)(Θ^AML∗) denotes the β quantile of Θ^AML,1∗,…,Θ^AML,B∗ (Davison and Hinkley, 1997, sections 2.2.1, 2.4).

3.1 Theoretical guidelines for choosing tuning parameters

The stochastic approximation part of our approximate maximum likelihood algorithm requires two tuning parameter sequences, the stepsize a_n and the parameter c_n responsible for the bias-variance trade-off of the finite difference approximations to the derivative. According to Spall (2003, chap. 7), tuning constants satisfying an,cn→0 for n → ∞, as well as ∑n=1∞an<∞ and ∑n=1∞an2cn2<∞ ensure convergence of the algorithm. (Notice that there are further technical conditions for convergence, such as a unique optimum and a uniformly bounded Hessian matrix.) It can be shown that for

optimum rates of convergence can be obtained when α = 1 and γ = 1/6. As discussed in Spall (2003, p. 187), however, smaller values for α and γ often lead to a better practical performance. The constant A is usually taken as a small fraction of the maximum number of iterations N and guards against too large jumps in the initial steps.

Another important factor affecting the algorithm is the precision of the kernel density estimate used. Under smoothness conditions on the underlying likelihood, the kernel density estimate used for estimating the likelihood is consistent, if kn→0, Hn→0 and nHnd→∞ in the case of d summary statistics. Notice however that for stochastic approximation to work, the variance of the estimates does not need to tend to zero. In particular large sample sizes k_n are not needed in the individual steps. Nevertheless good estimates of the gradient are obviously helpful. If we measure the quality of the pointwise density estimates by the the mean squared error (MSE), optimizing the MSE [see Rosenblatt (1991)] suggests to choose Hn=cn−(4+d), (for each diagonal enrtry of H_n) with a constant c^∗ depending on the unknown likelihood function.

As c^∗ is unknown in practice, we use a data driven bandwidth selection rule with our algorithm. Many different methods to estimate bandwidth matrices for multivariate kernel density estimation have been proposed (see Scott, 2015, for an overview), and in principle any method that fulfills the convergence criteria is valid to estimate H_n. For some types of densities, non-sparse bandwidth matrices can give a substantial gain in efficiency compared to e.g. diagonal matrices (Wand and Jones, 1995, section 4.6). However, as the bandwidth matrix is computed for each likelihood estimate, computational efficiency is important. In our examples, we estimate a diagonal bandwidth matrix H_n using a multivariate extension of Silverman’s rule (Härdle et al., 2004, equation 3.69). Using new bandwidth matrices in each iteration introduces an additional level of noise that can be reduced by using a moving average of bandwidth estimates.

In density estimation, the choice of a specific kernel function is usually considered less important than the bandwidth choice. The kernel choice plays a more important role when used with our stochastic gradient algorithm however. Indeed, to enable the estimation of the gradient even far away from the maximum, a kernel function with infinite support is helpful. The Gaussian kernel is an obvious option, but when the log-likelihood is used, the high rate of decay can by chance cause very large steps leading away from the maximum. Therefore, we use the following modification of the Gaussian kernel:

In degenerate cases where the likelihood evaluates to zero numerically, we replace the classical kernel density estimate by a nearest neighbor estimate in step 1(c):

If L^⁢(Θ−;sobs)≈0 or/and L^⁢(Θ+;sobs)≈0 (with “≈” denoting “numerically equivalent to”), find

and recompute the kernel density estimate L^⁢(Θ−;s) and L^⁢(Θ+;s) in step 1c using Smin− and Smin+, respectively, as the only observation.

3.2 Heuristic tuning guidelines

Here, we summarize some practical recommendations that turned out to be useful in our simulations.

As the parameters we want to estimate may often lie in regions of different length, and to reflect the varying slope of L in different directions of the parameter space, a good choice of a and c can speed up convergence. We therefore replace a∈ℝ+ in eq. (3) by a p-dimensional diagonal matrix a=diag(a(1),…,a(p)) with a(i)∈ℝ+ for i = 1, …, p (this is equivalent to scaling the space accordingly). The optimal choice of a and c depends on the unknown shape of L. Therefore, we propose the following heuristic based on suggestions in Spall (2003, sections 6.6, 7.5.1, 7.5.2) and our own experience to determine these values as well as A∈ℕ, a shift parameter to avoid too fast decay of the step size in the first iterations.

Let N be the number of planned iterations and b∈Rp a vector that gives the desired stepsize in early iterations in each dimension. Choose a starting value Θ0.

1. Set c to a small percentage of the total parameter range (or the search space, see Section3.4), to obtainc₁.

2. SetA=⌊0.1*N⌋.

3. Choosea:

   1. Estimate∇L(Θ0;sobs)by the median ofK₁finite differences approximations (step 1 of the AML algorithm),∇^¯c1L^(Θ0;sobs).

   2. Set

      a(i)=|b(i)(A+1)α(∇^¯c1L^(Θ0;sobs))(i)| for i=1,…,p.

As a is determined using information about the likelihood in mn_mn only, it might not be adequate in other regions of the parameter space. To be able to distinguish convergence from a too small step size, we simultaneously monitor the growth of the likelihood function and the trend in the parameter estimates to adjust a if necessary. Every N₀ iterations the following three tests are conducted on the preceding N₀ iterates:

• Trend test (too small a⁽ⁱ⁾):For each dimension i = 1, …, p a trend inΘk(i)is tested using the standard random walk model

  Θk(i)=Θk−1(i)+β+ϵk,

  whereβdenotes the trend andϵk∼N(0,σ2). The null hypothesisβ = 0 can be tested by a t-test on the differencesΔk=Θk(i)−Θk−1(i). If a trend is detected, a⁽ⁱ⁾is increased by a fixed factorf∈ℝ+.

• Range test (too large a⁽ⁱ⁾):For each dimension i = 1, …, p, a⁽ⁱ⁾is set toa(i)/fif the trajectory ofΘk(i)spans more thanπ_a% of the parameter range.

• Convergence test:Simulate K₂likelihood estimates atΘn−N0and at Θ_n. Growth of the likelihood is then tested by a one-sided Welch’s t-test. (Testing for equivalence could be used instead, seeWellek, 2010.)

We conclude that the algorithm has converged to a maximum only if the convergence test did not reject the null hypothesis three times in a row and at the same time no adjustments of a were necessary.

The values of the tuning parameters we used in our examples presented in Section 4 are summarized in Table 1.

3.3 Starting points

To avoid starting in regions of very low likelihood, a set of random points should be drawn from the parameter space and their simulated likelihood values compared to find good starting points. This strategy also helps to avoid that the algorithm reaches only a local maximum. More details can be found in Section 4.

3.4 Constraints on the parameters

The parameter space will usually be subject to constraints (e.g. rates are positive quantities). They can be incorporated by projecting the iterate to the closest point such that both Θ− as well as Θ+ are in the feasible set (Sadegh, 1997). Even if there are no imperative constraints, it is advisable to restrict the search space to a range of plausible values to prevent the algorithm from trailing off at random in regions of very low likelihood.

To reduce the effect of large steps within the boundaries, we clamp the step size at π_max% of the range of the search space in each dimension.

4.1 Multivariate normal distribution

One dataset, consisting of i.i.d. draws from a 10-dimensional normal distribution, is simulated such that the maximum likelihood estimator for Θ is Θ^ML=X¯~N(5·110,I10) where I₁₀ denotes the 10-dimensional identity matrix and 1₁₀ the 10-dimensional vector with 1 in each component. To estimate the distribution of Θ^AML-SP, the AML-SP algorithm is run 1000 times on this dataset with summary statistics S=X¯.

At start, 1000 points are drawn randomly on (−100, 100) × … ×(−100, 100). For each of them, the likelihood is simulated and the 5 points with the highest likelihood estimate are used as starting points. On each of them, the AML-SP algorithm is run with k_n = 100 for at least 10,000 iterations and stopped as soon as convergence is reached (for ≈90% of the sequences within 11,000 iterations). Again, the likelihood is simulated on each of the five results and the one with the highest likelihood is considered as a realization of Θ^AML-SP. Based on these 1000 realizations of Θ^AML-SP, the density, bias and standard error of Θ^AML-SP are estimated for each dimension (Table 2, Figure 1).

The results are extremely accurate: The densities of the 10 components of Θ^AML-SP are symmetric around the maximum likelihood estimate with a negligible bias. Compared to the standard error of Θ^ML, which by construction equals 1, the standard error of Θ^AML-SP is nearly 20 times smaller. Bootstrap confidence intervals that were obtained from B = 100 bootstrap samples for 100 datasets simulated under the same model as above show a close approximation to the intended coverage probability of 95% (Table 2).

To investigate the impact of the choice of summary statistics on the results, we repeat the experiment with the following set of summary statistics:

Figure 1:

Density of the components of Θ^AML-SP obtained with S=X¯ in one dataset estimated from 1000 converged sequences with a miminum length of 10,000 iterations by kernel density estimation. Vertical dashed line: Θ^ML.

For ≈90% of the sequences, convergence was detected within 14,000 iterations. Bootstrap confidence intervals are obtained for 100 simulated example datasets. Their coverage probability matches the nominal 95% confidence level closely (Table 2). The components behave very similar to the previous simulations, except for the estimates of the density for components 9 and 10 (Figure 2). Compared to the simulations with S=X¯, the bias of components 9 and 10 is considerably increased, but it is still much smaller than the standard error of Θ^ML. To investigate how fast the bias decreases with the number of iterations, we re-run the above described algorithm for 100,000 iterations without earlier stopping (Figure 3). Both bias and standard error decrease with the number of iterations.

Figure 2:

Density of the components of Θ^AML-SP obtained with S^∗ (eq. 4) in one dataset estimated from 1000 converged sequences with a miminum length of 10,000 iterations by kernel density estimation. Vertical dashed line: Θ^ML.

Figure 3:

Absolute bias and standard error of the AML estimator of μ₉ and μ₁₀ using the summary statistics in eq. (4) estimated from 1000 runs of the algorithm on the same dataset.

For the comparison of the AML-FD and AML-SP algorithm, 1000 datasets are simulated under the same normal distribution model as above. S^∗ (eq. 4) is used as our vector of summary statistics.

For each dataset, one random starting value is drawn on the search space (0,10)×…×(0,10) and the two approximate maximum likelihood algorithms AML-FD and SP are run on each of them with k_n = 50. Different values of c are tested in short preliminary runs and a small set of good values is used. Here, the convergence diagnostics and the methods to make the algorithms more robust are not used (except from shifting the iterates back into the search space if necessary).

For each algorithm, 5 million datasets are simulated. This results in different numbers of iterations (N): For the FD algorithm, N = 5000, for the SP algorithm, N = 50,000. The results of the different algorithms are compared after the same number of simulations. The corresponding runtimes are shown in Table 3. As expected, one iteration of the FD algorithm takes approximately 10 times longer than one iteration of the SP algorithm: the likelihood is estimated at 20 points in each iteration in the FD algorithm compared to 2 points in the SP algorithm.

In the worst case (c = 1, SP), more than 50% of the runs enter regions of approximately zero likelihood.

For c = 2, both algorithms behave very similar (Figure 4). For c = 1, the AML-SP algorithm converges slightly faster and for c = 0.5 it converges considerably faster and the AML-FD algorithm needs many more iterations to reach the same accuracy. This reflects its theoretical properties derived in Spall (1992).

However, the AML-SP algorithm is less stable than AML-FD and can more easily trail off randomly to regions of very low likelihood. The SP algorithm even enters regions of zero likelihood (at standard double precision) considerably more frequently (Table 3), causing problems with the computation of the log-likelihood. In our simulations we could easily overcome these problems by the adjustments proposed in Sections 3.2–3.4 that make the algorithm more robust. But in regions of very low likelihood, small differences in the likelihood values will still lead to large gradients of the log likelihood, and thus to large steps. Then, it can help to reduce π_max.

In the following examples, we use the AML-SP algorithm and to simplify the notation, define Θ^AML:=Θ^AML-SP.

Figure 4:

Euclidean distance of Θ^n to Θ^ML. The distance is plotted only every percent of the iterations. Right panel: same as left panel, but with both axes on a log-scale.

4.2 The evolutionary history of Bornean and Sumatran orang-utans

Pongo pygmaeus and Pongo abelii, Bornean and Sumatran orang-utans, respectively, are Asian great apes whose distributions are exclusive to the islands of Borneo and Sumatra. Recurring glacial periods led to a cooler, drier, and more seasonal climate. Consequently the rain forest might have contracted and led to isolated populations of orang-utans. At the same time, the sea level dropped and land bridges among islands created opportunities for migration among previously isolated populations. However, whether glacial periods have been an isolating or a connecting factor remains poorly understood. Therefore, there has been a considerable interest in using genetic data to understand the demographic history despite the computational difficulties involved in such a population genetic analysis. We will compare our results to the analysis of the orang-utan genome paper (Locke et al., 2011) and a more comprehensive study by Ma et al. (Ma et al., 2013); these are analyses that have been performed genome-wide. Both studies use ∂a∂i (Gutenkunst et al., 2009), a popular software that has been widely used for demographic inference. ∂a∂i is based on a numerical solution of the Wright-Fisher diffusion approximation. We use the orang-utan data set of Locke et al. (2011), consisting of SNP data of four-fold degenerate (synonymous) sites taken from 10 sequenced individuals (five individuals each per orang-utan population, haploid sample size 10). See the Appendix for more details on the dataset.

As in Locke et al. (2011) and Ma et al. (2013), we consider an Isolation-Migration (IM) model where a panmictic ancestral population of effective size N_A splits τ years ago into two distinct populations of constant effective size N_B (the Bornean population) and N_S (the Sumatran population) with backward migration rates μ_BS (fraction of the Bornean population that is replaced by Sumatran migrants per generation) and μ_SB (vice versa; Figure 5A).

Figure 5:

(A) Isolation-migration model for the ancestral history of orang-utans. Notation: N_A, effective size of the ancestral population; μ_BS (μ_SB), fraction of the Bornean (Sumatran) population that is replaced by Sumatran (Bornean) migrants per generation (backwards migration rate); τ, split time in years; N_B (N_S), effective population size in Borneo (Sumatra). (B) Binned joint site frequency spectrum (adapted from Naduvilezhath et al., 2011). (C) Folded joint site frequency spectrum of biallelic SNPs at four-fold degenerate sites in the Bornean and Sumatran orang-utan samples.

N_A is set to the present effective population size that we obtain using the number of SNPs in our data set and assuming an average per generation mutation rate per nucleotide of 2⋅10−8 and a generation time of 20 years (Locke et al., 2011), so N_A = 17,400.

There are no sufficient summary statistics at hand, but for the IM model the joint site frequency spectrum (JSFS) between the two populations was reported to be a particularly informative summary statistic (Tellier et al., 2011). However, for N samples in each of the two demes, the JSFS has (N+1)2−2 entries, so even for small datasets it is very high-dimensional. To reduce this to more manageable levels we follow Naduvilezhath et al. (2011) and bin categories of entries (Figure 5B). As the ancestral state is unknown, we use the folded binned JSFS that has 28 entries (see Figure 5C for the folded JSFS of the two population samples). To incorporate observations of multiple unlinked loci, mean and standard deviation across loci are computed for each bin, so the final summary statistics vector is of length 56.

The AML-SP algorithm with the JSFS as summary statistics was used together with the coalescent simulation program msms (Ewing and Hermisson, 2010). This allows for complex demographic models and permits fast summary statistic evaluation without using external programs and the associated performance penalties.

4.2.1 Simulations

Before applying the AML-SP algorithm to the actual orang-utan DNA sequences, we tested it on simulated data. Under the described IM model with parameters NB=NS=17,400, μBS=μSB=1.44⋅10−5 and τ=695,000, we simulated 25 datasets with 25 haploid sequences per deme, each of them consisting of 75 loci with 130 SNPs each. We define loci as unlinked stretches of DNA sequence, which are so short that recombination within a locus can be disregarded.

For each dataset, 25 AML estimates were obtained with the same scheme: 1000 random starting points were drawn from the parameter space; the likelihood was estimated with n = 40 simulations. Then, the five values with the highest likelihood estimates were used as starting points for the AML algorithm. The algorithm converged after 3000–25,000 iterations (average: ≈8000 iterations; average runtime of our algorithm: 11.7 h on a single core).

For the migration rates μ_SB and μ_BS, the per dataset variation of the estimates is considerably smaller than the total variation (Table 4). This suggests that the approximation error of the AML algorithm is small in comparison to the error of Θ^ML. For the split time τ and the population sizes N_S and N_B, the difference is less pronounced, but still apparent. For all parameters, the average bias of the estimates is either smaller or of approximately the same size as the standard error. As the maximum likelihood estimate itself cannot be computed, it is impossible to disentangle the bias of Θ^ML and an additional bias introduced by the AML algorithm.

Table 4:

Properties of Θ^AML in the IM model with short divergence time (τ=695,000 years).

μBS		True	1.44e−05
		Space	(1.44e−06, 0.000144)
		Mean	1.42e−05
		Median	1.39e−05
		Bias	−1.74e−07
		Mean se	2.1e−06
		Total se	4.09e−06
μSB		True	1.44e−05
		Space	(1.44e−06, 0.000144)
		Mean	1.26e−05
		Median	1.22e−05
		Bias	−1.78e−06
		Mean se	2.02e−06
		Total se	3.81e−06
τ		True	695,000
		Space	(139,000, 6,950,000)
		Mean	1,030,000
		Median	920,000
		Bias	334,000
		Mean se	310,000
		Total se	437,000
NS		True	17,400
		Space	(1740, 174,000)
		Mean	22,000
		Median	21,200
		Bias	4610
		Mean se	3060
		Total se	4130
NB		True	17,400
		Space	(1740, 174,000)
		Mean	21,600
		Median	20,600
		Bias	4230
		Mean se	3090
		Total se	4790

1. Figures: marginal densities of components of Θ^AML, estimated by kernel density estimation. Black solid line: density of Θ^AML. Coloured dotted and dashed lines: density of Θ^AML in three example datasets. Vertical black line: true parameter value. Summaries: true, true parameter value; space, search space; mean, mean of Θ^AML; median, median of Θ^AML; bias, bias of Θ^AML; mean se, mean standard error of Θ^AML per dataset; total se, standard error of Θ^AML.

As an alternative measure of performance, we compare the likelihood estimated at the true parameter values and the AML estimate with the highest estimated likelihood on each dataset (Table 5). In none of the 25 simulated datasets the likelihood at the true value is higher, whereas it is significantly lower in 11 of them (significance was tested with a two-sided Welch test using a Bonferroni-correction to account for multiple testing, i.e. for the 100 tests in Tables 5 and 7, a test with p-value <0.00025 is considered significant at an overall significance level of 5%). This suggests that the AML algorithm usually produces estimates that are closer to the maximum likelihood estimate than the true parameter value is.

Similarly, we compare our results to estimates obtained with δaδi, Θ^δaδi (Gutenkunst et al., 2009). As δaδi is based on the diffusion approximation, whereas we are using the coalescent, we do not expect the maximum likelihood estimators to be equal, and the coalescent likelihood should be higher at Θ^AML, if the algorithm performs well. In all the 25 simulated datasets the coalescent likelihood at Θ^δaδi is lower than at Θ^AML with the highest likelihood, and the difference is significant in 13 datasets (Table 5). To make sure that the results are not biased due to possible convergence problems of δaδi, we use the true Θ as starting point. In none of the runs, δaδi reported convergence problems.

Table 5:

Comparison of likelihood estimates at Θ, Θ^AML and Θ^δaδi in the IM-model with short divergence time (τ=695,000 years).

1. Rank: rank of L^(Θ) or L^(Θ^δaδi), respectively, among L^(Θ^AML,1),…,L^(Θ^AML,25), where Θ^AML,j denotes the j’th realization of Θ^AML in the dataset at hand; p-value: p-value of Welch’s t-test for H0:L(Θ)=L(Θ^AML,max) vs. H1:L(Θ)≠L(Θ^AML,max) with Θ^AML,max=argmax{L^(Θ^AML,j):j=1,…,25}; boxplot: boxplot of logL^(Θ^AML,1),…,logL^(Θ^AML,25); grey ×, L^(Θ); grey +, L^(Θ^δaδi).

To investigate the impact of the underlying parameter value on the quality of the estimates, simulation results were obtained also for 25 datasets simulated with the divergence time τ twice as large. Here, all parameter estimates, especially τ, N_B and N_S had larger standard errors and biases (Table 6). Apparently, the estimation problem is more difficult for more distant split times. This may be caused by a flatter likelihood surface and by stronger random noise in the simulations. Only the migration rates are hardly affected by the large τ: a longer divergence time allows for more migration events that might facilitate their analysis. The higher level of difficulty shows up also when comparing the likelihood at the true value and at Θ^AML: in 13 out of 25 datasets, the likelihood at Θ is significantly lower than at Θ^AML, whereas it is significantly higher in three cases (Table 7). The likelihood at Θ^δaδi is significantly lower in 15 datasets and never significantly higher.

Table 6:

Properties of Θ^AML in the IM model with long divergence time (τ=1,390,000 years).

μBS		True	1.44e−05
		Space	(1.44e−06, 0.000144)
		Mean	1.24e−05
		Median	1.12e−05
		Bias	−1.97e−06
		Mean se	3.38e−06
		Total se	4.99e−06
μSB		True	1.44e−05
		Space	(1.44e−06, 0.000144)
		Mean	1.25e−05
		Median	1.17e−05
		Bias	−1.88e−06
		Mean se	3.28e−06
		Total se	4.38e−06
τ		True	1,390,000
		Space	(139,000, 6,950,000)
		Mean	2,360,000
		Median	2,280,000
		Bias	967,000
		Mean se	876,000
		Total se	1e+06
NS		True	17,400
		Space	(1740, 174,000)
		Mean	25,700
		Median	24,600
		Bias	8350
		Mean se	6620
		Total se	7560
NB		True	17,400
		Space	(1740, 174,000)
		Mean	26,000
		Median	25,200
		Bias	8670
		Mean se	6680
		Total se	7870

1. Figures and summaries as in Table 4.

Table 7:

Comparison of likelihood estimates at Θ, Θ^AML and Θ^δaδi in the IM-model with τ=1,390,000 years.

1. Figures and summaries as in Table 5.