The Parameters of the Barry and Hartigan General Markov Model Are Statistically NonIdentifiable

Phylogenetic estimation using model-based distance, maximum likelihood (ML), and Bayesian approaches have exploded in popularity in the last decade. Markov models employed by these approaches (e.g., the GTR model) typically assume that nucleotide or amino acid sequences evolve through a homogeneous, stationary, or reversible process over edges of the Tree of Life (Felsenstein 2004). However, in some cases, molecular sequence evolution will fail to satisfy these assumptions; nucleotide or amino acid frequencies can sometimes vary greatly between sequences under examination, and the use of Markov models that fail to account for this property can positively mislead phylogenetic estimation (Galtier and Gouy 1995, Foster and Hickey 1999). Jayaswal et al. (2005) have showed that, for such cases, a more general Markov model proposed first by Barry and Hartigan (1987) (the “BH” model) is a useful alternative to the GTR family of models because it makes many fewer assumptions, namely that 1) evolution follows a Markov process across each edge and, 2) the data patterns at each site are independently and identically distributed. As a result, the BH model allows the character frequencies to differ across nodes in the tree and the evolutionary process to differ not only across edges but also along the two directions on an edge. Here, we show that despite the fact that the BH model has been shown to improve accuracy of phylogenetic estimation under some conditions, there are potentially serious problems with the identifiability of its parameters. An identifiability problem arises under the BH model, because, as shown in detail below, different sets of parameters can lead to the same probability distribution of data patterns under the model and hence the same likelihood of the data. Specifically, we will show that there always exist multiple different sets of transition probabilities along edges and different sets of character base frequency vectors at internal nodes that will lead yield the same likelihood. If the BH model is being used to infer the tree topology alone, this will not create a problem. This follows from the main result of Steel (1994) which shows that LogDet distances are tree additive for the BH model, implying that sets of pairwise distributions for different topologies are always different. However, if researchers are interested in the nucleotide transition probabilities along edges of the phylogeny, the nucleotide frequencies at internal nodes, or approximations to edge lengths as, for example, proposed by Jayaswal et al. (2005), nonidentifiability becomes a major problem. Below, we explain the nature of this nonidentifiability and the effect it has on estimation of parameters. The three taxa tree used in simulations. The Q( e ) are oriented so that rows refer to the state at the internal node. The simulating model has the same Q( e ) for each edge determined from a GTR model with edge length 0.5, stationary frequencies 0.1,0.35,0.4,0.15, and exchangeabilities 0.1,1,0.5,0.5,1 and 0.1 for AC,AG,AT,CG,CT. and GT. Equation 2 implies that the probabilities of data patterns with and without switching the rows of the conditional probabilities along three edges and the base character frequencies at internal node i remain the same. Therefore, this model is not identifiable. Note that the same permutation was applied to each of the transition matrices. More generally, for sequence data with m different states and three taxa, there are m! permutations that have the same site pattern distribution. Furthermore, for any topology that has more than one internal node, the problem is compounded because multiple permutations of transition matrices can be performed at each internal node without a change in the likelihood of the data. To demonstrate how and why the nonidentifiability of the BH model may create problems, we simulated large sequence length data sets under the widely used GTR model, which is a special case of the BH model. For these simulated data sets, we know the true generating parameters, and thus we can make a comparison between the true values and the estimates of parameters. With long sequence lengths, uncertainty due to estimation is minimal, and therefore large departures of estimates from the true values can be attributed to nonidentifiability effects. For the simulation, we selected a very simple three-branch tree (Fig. 1) and simulated data with the Seq-Gen program (Rambaut and Grassly 1997) under a GTR model with arbitrarily chosen parameters: as a stationary nucleotide frequency vector of π = {0.10,0.35,0.4,0.15}, state exchangeabilities of r = {0.1(A↔C),1(A↔G),0.5(A↔T),0.5(C↔G),1(C↔T),0.1(G↔T)}, and a sequence length of 100,000 sites. For this data set, we estimated parameters of the BH model using the software described in Jayaswal et al. (2005). Here, S is the set of 24 permutations of rows, E contains the three branches, and Q^i,j(e)s is the sth permutation of Q^(e). SSQs gives the minimum departure of estimates from the true matrices over all permutations of the rows of the estimated matrices; the same permutation is applied for each edge. If minimum sum of squares is attained without requiring permutation, the estimated permutation is the “correct” permutation, otherwise it is a “wrong” permutation. For the large sequence length data sets we considered, only one permutation of the matrices gave a small sum of squares. Here, we show the results for a simulated data set with 100,000 sites. To avoid difficulties with local maxima and to investigate whether the correct permutation tends to get estimated in practice, we used 100 sets of randomly generated matrices along the three branches in Figure 1 as initial values for ML estimation and obtained 100 sets of estimates for the single data set. Of these 100 sets of estimates, 97 converged on the global optimum with log-likelihoods within 0.5 of − 349,111.0, whereas three of them were local maxima with log-likelihoods less than -349,366.7. Below, we restrict attention to the 97 cases where the global optimum was estimated. For these 97 estimates, SSQss ranged from 0.000067 to 0.00018 indicating that the actual values of the entries of estimated matrices were very close to the entries in the true generating matrices once they were permuted appropriately. However, for 93 of the cases, the estimated matrices did not correspond to the original permutation of the matrix; the estimated permutations covered 23 of the 24 possible permutations (Table 1). Thus, the nonidentifiability problem with the BH model arises because the correct permutation of the estimated matrices cannot be recovered in practice. The numbers of times that each of the 24 permutations were the estimated permutation among the 97 runs The numbers of times that each of the 24 permutations were the estimated permutation among the 97 runs One of the potential benefits of the BH model is that it allows the researcher to both accommodate changing state (nucleotide or amino acid) frequencies during phylogenetic estimation as well as estimate the ancestral character state frequency vector at internal nodes of the tree. The latter parameters are of general interest to molecular evolutionists as they may give insights into the physical properties of ancestral molecules and environments (see, e.g., Boussau et al. 2008). Unfortunately, one result of the nonidentifiability problem we describe is that the ancestral base character frequency vector of the internal node will also only be accurately estimated up to permutation. To illustrate this, we estimated the base character frequency vectors of internal node i for the five most frequent permutations recovered above (Table 2). Because incorrect permutations were estimated, the estimated base frequencies are drastically different from the true values with an averaged sum of square error of approximately 0.13 as compared with 0.002 for the correct permutations. The mean vectors of estimates of internal node frequency vectors for the five most frequent permutations The mean vectors of estimates of internal node frequency vectors for the five most frequent permutations Another serious problem can occur if researchers wish to estimate edge lengths for their phylogenetic tree based on the estimated BH model parameters. For the BH model, the entries in the joint probability matrices along edges are optimized directly. No edge length estimation is involved as the evolutionary processes along edges need not correspond to a continuous-time Markov process. However, Jayaswal et al. (2005) described a way to extract an approximate branch length estimate from the BH model parameters by assuming a continuous-time Markov process. A difficulty arises because the method requires logarithms of all the eigenvalues of the conditional probability matrices of the two evolutionary directions along a branch to exist. Unfortunately, we find that for most estimated permutations this is not the case. For example, for the eigen decompositions of the transition matrices along branches of 24 permutations of the example above we found that typically only one permutation has all positive eigenvalues, whereas the others all have complex and/or negative-valued eigenvalues. For idealized simulation cases where the generating model is GTR and there is sufficient data, the correct permutation of the BH estimate will correspond to a case where there the conditional probability matrices have positive eigenvalues, as in our example. However, for real data that have evolved under heterogeneous historical conditions there is no guarantee that any of the permutations of these matrices will meet this condition or, if they do, that they correspond to the correct permutations. For example, we fit the BH model to the Plasmodium species phylogenomic data set and tree described by Davalos and Perkins (2008) that had eight taxa and 77,313 sites after gaps were removed. We estimated the parameters using the BH model with random initial starting values. For each of the 13 edges of the tree in Figure 1(A) of Davalos and Perkins (2008), we found that the transition matrices always had negative eigenvalues even after ignoring the complex-valued portions of the eigenvalues. Thus, the approximate method of Jayaswal et al. (2005) could not be used to obtain edge lengths for this case. From the foregoing discussion, it should be clear that for BH model parameter estimates to be useful to researchers, a method for estimating the correct permutation is necessary. One possibility is given in the implementation of Jayaswal et al. (2005), where they recommended initial parameter values for BH optimization be 18 on the diagonal and 124 on the off-diagonal elements. This satisfies the condition of “diagonal largest in column” (DLC) for transition matrices in both evolutionary directions. The DLC condition is discussed in Chang (1996) as a potential condition for identifiability. If the true evolutionary process yields transition matrices with this property, then starting the optimization with initial parameters that satisfy DLC could help with the estimation of the correct permutation. We did an experiment randomly generating 100 data sets each with 1000 sites under our generating GTR model which does satisfy the DLC condition. For 100 sets of estimates of transition probability matrices obtained using Jayaswal et al. recommended initial values, we always obtained the correct permutation. However, in practice, part of the reason for considering a BH model is to fit nonstationary processes; in the latter, there is no guarantee that the DLC condition will hold for the correct permutation. To demonstrate this, we simulated under a nonstationary model using the tree in Figure 1. For this generating model, the frequency vector was {0.45,0.05,0.05,0.45} for a GTR model along edge (a,i); the frequency vector was {0.05,0.45,0.45,0.05} for edges (i,b) and (i,c). All three edges shared the same exchangeability vector {1.0, 5.6, 1.0, 1.0, 5.6, 1.0} for A↔C,A↔G,A↔T,C↔G,C↔T,G↔T and the root frequency vector is {0.1,0.4,0.4,0.1}. The model was nonstationary because different GTR models were used on the other edges. Here, although the transition matrices do satisfy the DLC condition, the joint probability matrices do not. In this experiment, the estimates obtained based on optimization from Jayaswal et al. recommended initial parameter settings for the joint probability matrices did not correspond to the correct permutations and the eigenvalue vectors of the estimated transition matrices along the three edges had negative values. Testing our idea, we used the estimates obtained from our previous simulation experiments and estimated the permutation using Equation 4. We found that in all cases examined the sums of squares of the estimated permutations were smaller than others. Secondly, we checked whether the s obtained by this strategy was the true correct permutation and verified that for the cases where the estimation converged on global optimal likelihood, the correct permutation was always selected by our method. In general, we expect that this method will perform well for real data sets where the nucleotide composition of sequences is changing gradually (rather than abruptly) over the tree topology. To extend our proposal to more than three taxa, one could minimize the sum of Equation (4) over all internal nodes. One approach to doing so is as follows. Given an ordering of the internal nodes, successively to determine the permutations minimizing Equation (4) for each of the internal nodes. For each internal node, when minimizing, hold the frequency vectors at all other internal nodes fixed. Iterate the process until no further improvement in the sum is possible. Iteration is required here because the best permutation for node i depends on the frequencies at a. Consequently, if the permutation at node a changes when it is considered, the permutation for node i minimizing Equation (4) might change as well. Although the phylogenetic tree estimated by the BH model is an identifiable parameter, we have shown that there is a problem with identifiability of the joint probability matrix parameters. If ignored, this nonidentifiability problem can mislead researchers interested in the ancestral character state compositions estimated by the BH method or approximate edge lengths generated by existing implementations. We have proposed a solution that will select the correct permutation when the data evolves in a “close to stationary” manner. However, if this assumption is not correct, then researchers should be aware of the potential pitfalls stemming from the identifiability problem we have discussed. This work was supported by Discovery grants from the Natural Sciences and Engineering Research Council of Canada awarded to C.F., E.S., and A.J.R. We would like to thank Dr V. Jayaswal for helpful advice and discussions concerning the BH model and software implementation. We would like to thank Drs J. Robinson, P. Foster, and an anonymous reviewer for many useful suggestions.

Fetched live from OpenAlex and de-inverted. Abstracts are not stored in this database: the inverted indexes are 8.6 GB of the frame’s 9.3 GB of text, and the host has 13 GB free.