Statistical Tests of Gamma-Distributed Rate Heterogeneity in Models of Sequence Evolution in Phylogenetics

Likelihood ratio tests (LRTs) for comparing models of sequence evolution have become popular over the last few years (Goldman 1993 ; Yang, Goldman and Friday 1994, 1995 ; Huelsenbeck and Crandall 1997 ; Huelsenbeck and Rannala 1997 ). In their simplest form, such tests compare a simpler null hypothesis (H0) with a more complex alternative hypotheses (H1) which is a generalization of H0. H0 can be derived from H1 by fixing one or more of its free parameters at particular values, and the hypotheses are described as nested. Although it is also possible to test non-nested models (Goldman 1993 ), nested models are often preferred, as statistical tests are simpler to perform and their results can be easier to interpret. The test statistic for an LRT can be written as 2δ = 2 ln(L̂H1/L̂H0) = 2(ln(L̂H1) − ln(L̂H0)), where L̂H0 and L̂H1 are the maximum-likelihood (ML) scores under hypotheses H0 and H1, respectively. This statistic measures how much improvement H1 gives over H0, and when the hypotheses are nested, 2δ will always be nonnegative. For these nested hypotheses, and under certain regularity conditions, the asymptotic distribution of 2δ (i.e., for large amounts of data) will be χ2k. Here, k is the number of degrees of freedom by which H0 and H1 differ, that is, the number of free parameters of H1 whose values must be fixed to derive H0 (Wald 1949 ; Silvey 1975 ; Felsenstein 1981;Goldman 1993 ; Yang, Goldman, and Friday 1994, 1995 ). (In effect, each free parameter contributes a χ21 variate to the distribution of 2δ, with the sum of k independent χ21 variates being distributed as χ2k.) Statistical tests assessed using such χ2 distributions have now become a widespread and useful tool in phylogenetics (Huelsenbeck and Crandall 1997 ; Huelsenbeck and Rannala 1997 ). Recently, there has been renewed interest in testing whether the predicted χ2 distribution gives a reliable estimate of the true distribution of 2δ under realistic conditions (e.g., with finite sequence lengths). Whelan and Goldman (1999) investigated cases in which the competing hypotheses were different models of nucleotide substitution. Under three specimen experimental designs (representing realistic phylogenies and nucleotide substitution processes), we found that the χ2 distribution was acceptable for performing tests of the significance of parameters describing the relative rate of transition and transversion substitutions (typically denoted κ; one free parameter) and describing equilibrium base frequencies (typically denoted π = (πA, πC, πG, πT); three free parameters, because of the constraint that ΣX πX = 1) when these parameters are estimated by ML (Whelan and Goldman 1999 ). Combinations of these parameters give most of the commonly used models of nucleotide substitution. The parameter κ is present in the models K2P of Kimura (1980) and HKY of Hasegawa, Kishino, and Yano (1985) , which reduce to the models JC of Jukes and Cantor (1969) and FEL of Felsenstein (1981) , respectively, when κ is fixed at 1. The parameters π are present in FEL and HKY, which reduce to JC and K2P, respectively, when π is fixed at (¼, ¼, ¼, ¼). In contrast, for tests of the significance of the incorporation of the shape parameter of a gamma distribution used to describe among-sites variation in the rate of nucleotide substitution (Yang 1993, 1994, 1996; typically denoted α; one free parameter), Whelan and Goldman (1999) consistently found an unacceptable fit between the predicted χ2 distributions and estimates of the true distributions of 2δ. We provisionally attributed this effect to the fact that reducing a hypothesis H1 incorporating a gamma distribution (denoted, e.g., JC+Γ) to a hypothesis H0 without a gamma distribution (in this example, JC) requires α to be fixed at ∞, which is on the boundary of the set [0, ∞) of values permitted under H1. This contravenes the conditions required for the asymptotic distribution of 2δ to be χ2 (see, e.g., Self and Liang 1987 ). Ota and colleagues (Ota, Waddell, and Kishino 1999 ; Ota et al. 2000 ) have also investigated effects caused by parameters lying on their permitted boundaries. They illustrated that ML estimates of these parameters cannot be asymptotically normally distributed and how, as a consequence, the ML estimates must be biased and the distribution of 2δ (2lnLR in their notation) cannot be χ2. Andrews (2000) showed that parameters which lie on their permitted boundaries can also adversely affect results of nonparametric bootstraps. We refer to a parameter which under H0 is fixed at a value on the boundary of the set of values permitted under H1 as a boundary parameter. Self and Liang (1987 , p. 608, case 5) showed that (in the absence of any other boundary parameters not being tested) the contribution to 2δ of a single boundary parameter is (asymptotically) a distribution which is a 50:50 mixture of χ20 and χ21 distributions. The χ20 distribution takes the value 0 with probability 1; consequently, this mixture distribution takes the value 0 with probability ½, and takes a value drawn from a χ21 distribution with probability ½. If the statistical test being performed also tests for the significance of k − 1 other parameters which are not on their boundaries, these contribute χ2k−1 to the distribution of 2δ, and the complete distribution of 2δ is a 50:50 mixture of χ2k−1 and χ2k (Self and Liang 1987 , p. 609, case 9). We denote this distribution by χ̄2k. These are the situations of relevance to the phylogenetic tests described here. The χ̄2k distribution has a mean of k − ½; significance levels of observed statistics X may be calculated according to P(χ̄2k > X) = (P(χ2k−1 > X) + P(χ2k > X))/2, where the terms to the right of the equality can be found using standard tables or computer programs (e.g., the chi2 program of the PAML package; Yang 1997 ). Other features of χ̄2 distributions are more difficult to derive and were calculated using locally written programs. If two or more parameters lie on their boundaries, the situation is more complicated (Self and Liang 1987 ; Ota et al. 2000 ), and if these parameters are not independent, the distribution of 2δ may not be expressible as any sum of χ2 distributions (e.g., Self and Liang 1987 , pp. 608–609, case 8). Models incorporating more than one boundary parameter are beginning to appear in phylogenetic analyses (e.g., Huelsenbeck and Nielsen 1999 ). Using simulated data, Whelan and Goldman (1999) performed a number of model comparisons to investigate the true distribution of 2δ when the hypotheses being compared differed by various combinations of the parameters κ, π, and α. We now report the results of further simulations designed to consider specifically whether the χ̄2 distributions described by Self and Liang (1987) and Ota et al. (2000) are acceptable under realistic conditions for statistical comparisons of nested models where H1 incorporates a gamma distribution for rate heterogeneity and H0 does not. We again use the models JC, K2P, FEL, and HKY (using ML estimators of π in FEL and HKY, denoted FELMLE and HKYMLE in Whelan and Goldman 1999 ) and their counterparts, which, in addition, include the gamma distribution (JC+Γ, K2P+Γ, FEL+Γ, and HKY+Γ). Nine model comparisons (labeled A–J in table 1 ) represent all possible combinations of parameters which may then be present in H0 and H1, subject to the condition that α always be included in H1 but not H0. Thus, each comparison is always testing for the significance of including α, either singly (comparisons A–D) or in combination with other parameters (E–J), and either with (B–D, H–J) or without (A, E–G) other (untested) parameters present. Model trees used for data simulation are those described in Whelan and Goldman (1999) , labeled according to the original data set whose analysis gave the phylogeny and parameters used for data simulation: cytochrome b, ψη-globin, and mtDNA (erroneously labeled “D-loop” by Whelan and Goldman 1999 , although, in fact, it comprises two protein-coding regions and three tRNAs). Table 1 summarizes our results on the fit of the values of 2δ obtained from the analysis of simulated data sets to the predicted χ̄2 distributions. For all but one of the 27 comparisons in table 1 (test of HKY vs. HKY+Γ for the cytochrome b model phylogeny), we observe a good fit between the distributions. Note that this is true irrespective of the combination of parameters being tested and of the combination of parameters present in both models and thus untested. Ota et al. (2000) comment on practical difficulties that can arise in iterative ML optimizations. These are due to likelihood surfaces which are almost flat, and may be particularly pronounced in cases where optimal parameter values are on or near the boundaries permitted. Possible effects are premature termination of ML iterations either before a parameter estimate reaches a boundary at which it should lie or before it escapes from a boundary at which it should not lie. We encountered some such problems in a few of the analyses reported in table 1 . Typically, difficulty in achieving an accurate ML score is greater for more complex models (e.g., H1 rather than H0) and for analyses with parameters estimated near permitted boundaries. These effects lead to slightly depressed estimates of 2δ when it is near 0. We find that these effects account for our example in which the fit between our estimate of the distribution of 2δ and the predicted χ̄2 distribution is rejected as inadequate (results not shown), although even in this case the estimated mean and 95% point of the distribution of 2δ are close to those of the predicted χ̄2 distribution (table 1 , comparison D for cytochrome b). We conclude that the theoretical (asymptotic) prediction that the parameter α contributes a χ̄21 distribution (50:50 mixture of χ20 and χ21) to the distribution of 2δ is accurate in these examples, even for the finite-sized data sets considered, and that χ̄2k distributions (50:50 mixtures of χ2k−1 and χ2k) are appropriate for performing tests involving the boundary parameter α and k − 1 additional parameters whose fixed values under H0 do not lie on the boundary permitted under H1. Combining the results described by Ota and colleagues (Ota, Waddell, and Kishino 1999 ; Ota et al. 2000 ), those described by Whelan and Goldman (1999) , and those described above, it seems certain that the use of a χ21 distribution for testing the shape parameter α of the gamma distribution or other boundary parameters is inappropriate and potentially misleading. As explained by Ota et al. (2000) , the past use of a χ21 distribution instead of the correct χ̄21 distribution (or the use of a χ2k distribution instead of a χ̄2k distribution in tests involving k − 1 other parameters) may have led to overly conservative conclusions (i.e., too great a tendency to reject H1). Although confidence in past results regarding rate heterogeneity will often be maintained (typically rejecting H0 when only H1 uses a gamma distribution to allow for such rate heterogeneity; Yang 1996 ), this need not always be the case. As more complex and realistic models are compared, including some which use gamma distributions to model factors other than rate heterogeneity, observed values of 2δ may be a lot smaller. For example, the failure to use appropriate χ̄2 distributions leads to two false negative results in Huelsenbeck and Nielsen (1999 , table 2 , data set COI, second codon position, test of HKY85 vs. HKY85+Γκ, and NADH6, third codon position, HKY85+Γκ vs. HKY85+Γr+Γκ—although we note that these comparisons are in no way crucial to those authors' findings). From the theory presented by Self and Liang (1987) and the results shown by Ota and colleagues (Ota, Waddell, and Kishino 1999 ; Ota et al. 2000 ) and those shown in this paper, we conclude that the use of χ̄2 distributions is appropriate for LRTs of boundary parameters in phylogenetics. These modified distributions result in lower thresholds for rejection of H0 in favor of H1, and appropriate values for significance testing are presented in table 2 . in Whelan and Goldman (1999) , we suggested that simulation always be used for tests incorporating the parameter α modeling rate heterogeneity. We are now able to recommend instead the (time-saving) use of table 2 for all future LRTs concerning the presence of one boundary parameter (and in the absence of other boundary parameters not being tested). There have been few cases to date of false negatives due to inappropriate distributions being used to test significance, but as models of sequence evolution become more complex and the relevant values of 2δ decrease, the use of the correct distribution for significance testing becomes increasingly important. Mike Hendy, Reviewing Editor Keywords: likelihood ratio tests gamma distribution maximum likelihood model comparison molecular evolution phylogenetics Address for correspondence and reprints: Nick Goldman, University Museum of Zoology, Department of Zoology, University of Cambridge, Downing Street, Cambridge CB2 3EJ, United Kingdom. E-mail: n.goldman@zoo.cam.ac.uk. Table 1 Probabilities and Other Information Relating to the Fit of Predicted χ̄2 Distributions to Estimates of the True Distribution of 2δ Table 1 Probabilities and Other Information Relating to the Fit of Predicted χ̄2 Distributions to Estimates of the True Distribution of 2δ Table 2 Percentage Points of χ̄2 Distributions Table 2 Percentage Points of χ̄2 Distributions N.G. is supported by a Wellcome Trust Fellowship in Biodiversity Research. S.W. is supported by a BBSRC Research Studentship. We are grateful to R. Ota and P. J. Waddell for allowing us to see prepublication versions of Ota et al. (2000), and to Z. Yang for alerting us to the method of calculating P(χ̄2k > X).