Modeling linguistic evolution: a look under the hood

3.1 Computing the likelihood of a tree under an evolutionary model: the pruning algorithm

At the outset, we do not know the latent parameters α and β (even if we can assume a fixed tree topology), and must find values which maximize the likelihood of the tree (under a Maximum Likelihood approach), or values in regions of high posterior probability (under a Bayesian approach). This process involves computing the probability of the tree, which requires summing over all possible state combinations of the interior nodes of the tree; this number increases exponentially with the number of interior nodes in the tree, making it computationally intractable for large trees. However, the independence structure of the tree — the fact that the state of a node depends only on the state of its direct ancestor and the length of the branch along which it evolves — can be exploited to efficiently compute this likelihood according to the pruning algorithm (Felsenstein 2004:253–5), which starts at the tips of the tree (i.e., the terminal nodes, which represent languages with observed data) and recursively computes the likelihood that each interior node of the tree is in state 0 or 1 (in the case of a binary character), given a tree topology and some values for α and β (which I write ℓα,β(n=i); this notation varies widely across the literature on biology and graphical statistics), given the likelihood that each of its children are in either state 0 or 1 (branchn,child denotes the length of the branch connecting n to one of its children).

At the tips of the tree, ℓ(state(n)=i) is equal to either 0 or 1 (provided that data are not missing). The overall likelihood of the tree is computed via the sum of the probabilities that the root is in each state (given its direct descendants) weighted by the stationary probability of each state^[2]:

The pruning algorithm is also a crucial ingredient in searching the space of possible tree topologies for regions of higher likelihood.

Yang (2014:133–143) gives an extensive survey of maximum likelihood methods for finding rates of evolution. In the case of a single, symmetric rate of evolution, the maximum likelihood estimate can be found via a simple line search, since only one variable is involved. For more than one rate, the math is slightly more complicated. MCMC approaches do not require much in the way of complicated math; different values for the rates in question are iteratively proposed, and accepted stochastically if they improve the likelihood of the tree (this procedure is described in greater detail below). MCMC approaches are generally Bayesian, and return posterior distributions over rates.

Via Bayes’ Rule, the relationship between the posterior probability p(θ|D), the likelihood p(D|θ), the prior p(θ), and the marginal likelihood is the following (θ denotes the unknown parameter[s], while D denotes the observed data):^[3]

The marginal likelihood can be difficult to compute, as it requires summing over or integrating out the entire hypothesis space (p(D)=∑θp(D,θ)=∑θp(D|θ)p(θ), if θ is discrete); since it does not vary from hypothesis to hypothesis, it can be ignored, as the relative probability of a hypothesis is as important (for most practical purposes) as its absolute probability.

While inferring distributions for α and β via MCMC, we care about the relative posterior probabilities (which, in this case, are equivalent to the relative likelihoods) of different values of α and β. It is standard to randomly initialize α and β; the likelihood of the tree under the current values of α and β (ℓα,β(tree)) is computed via the pruning algorithm, and new values α′,β′ are proposed, drawn from Gaussian distributions centered around α,β:

The likelihood under α′,β′, ℓα′,β′(tree), is computed as well. The proposed values α′,β′ are accepted as new current values if a randomly generated number a (drawn from the Uniform distribution between 0 and 1) is less than ℓα′,β′(tree)ℓα,β(tree).^[4] Above, σα,σβ denote the step sizes for the two parameters. These values are nontrivial, as they have an effect on the extent to which the space of possible states is explored; in general the state space is adequately traversed by a Markov chain when roughly 23.4% of proposed states are accepted (Rosenthal 2011). After a period of time, the algorithm can be expected to draw samples from the posterior distributions of α and β.

Inference is generally run on three or four chains, all of which are initialized separately at different values, but under optimal circumstances are ultimately expected to “mix,” i.e., to sample from the same posterior distribution. An initial quantity of samples (in conservative practice, the first half) is generally discarded (the so-called “burn-in”). The chains’ convergence can be monitored using the Potential Scale Reduction Factor R^ for each posterior sample, which is the square root of the ratio of the between-chain variance and the average within-chain variance, under the assumption that the posterior distribution is normal (Gelman and Rubin 1992). Values below 1.1 indicate good convergence.

3.2 Ancestral state reconstruction and stochastic character mapping

Once maximum likelihood estimates or posterior distributions for α and β are obtained, it becomes possible to reconstruct the probability that a feature is present at internal nodes of the tree, and also to simulate mutations along branches, a process known as stochastic character mapping (Nielsen 2002; Bollback 2006). First, the probability that the feature is present at the root is calculated:

A Bernoulli variate parameterized by this probability is drawn, yielding a value i = {1, 0}. Subsequently, for each child n of the root, the probability that the feature is present at n given the sample value i for the root can be computed:

Again, a Bernoulli variate is drawn; this procedure carries on down the tree until the nodes directly above the tips are reached. After state values have been simulated at the interior nodes of the tree, the birth and death locations can be simulated along branches by drawing values from exponential distributions parameterized by α and β, respectively.

3.3 Extensions and linguistic applications

The CTM model of character evolution and the pruning algorithm serve as the backbone for a large number of models of evolution which have enjoyed use in linguistics. The model described above is among the most simple; it assumes a “strict clock,” i.e., that across all branches in the tree, α and β are drawn from the same global distributions. This simplicity and restrictiveness make for tractable and fast inference and prevent against overfitting, but the assumptions involved may not be entirely realistic. A “relaxed clock” assumes that each branch has its own individual rate parameter, which can be multiplied by the global change rates.^[5] The covarion family of models assumes that there are regions of slow and fast change over the tree; in some formulations, the slow rates are usually equivalent to the non-slow gain and loss rate multiplied by some constant between 0 and 1 (Wang et al. 2006). Hidden rates approaches extend this model to more than two rate classes (Beaulieu and O’Meara 2014).

Methodologies which infer evolutionary rates have enjoyed some use to date in diachronic linguistics. Dediu (2010) uses (among other methods) a two-rate model to estimate the stability of typological features in a number of families, quantifying stability according to the gain rate relative to the loss rate. Chang (2014) introduces a latent-state character model suited to the evolution of root-meaning characters. Under this model, a character is born at a certain rate, after which point it remains “latent” for a period of time; it becomes “active” at a certain rate, or dies; if active, it dies at a certain rate.

Pagel (1994) and Pagel and Meade (2006) present methods designed to measure co-evolutionary dependencies between a pair of traits. Two inference procedures are carried out: first, an independent model is run, where each trait is gained or lost independently according to a CTM process with four rates (two per trait). This is followed by a dependent model where the joint evolution of the two traits is modeled via a CTM process transitioning between four states representing all 2 × 2 combinations of presence and absence; the dependent model has eight rates, which treat all possible transitions between combinations differently (with the exception of the changes 0,0→1,1 and 1,1→0,0, the rates of which are set to zero). In Bayesian approaches to this procedure, evidence in favor of the dependent hypothesis versus the independent hypothesis is assessed according to the Bayes Factor of the two models (the ratio of their marginal likelihoods). This method provides an explicit approach to testing hypotheses regarding evolutionary associations between two features, provided that all possible combinations are observed at the tips of the tree. At the same time, Maddison and FitzJohn (2015) caution against making inferences regarding the causal relationship between the features under study; it could be the case that they are associated with each other because they are caused by a third factor (cf. Pearl 2009); furthermore, the authors note that this methodology will infer a correlation between features even if one of the features evolves only once. This technique is applied by Dunn et al. (2011) to the question of whether pairs of features involved in Greenbergian word order harmonies exhibit dependent evolution, with a negative result.

We may additionally desire to uncover correlated evolution between more than two traits. However, it is impossible to run the procedure described above for all possible combinations of a large number of traits. Reversible-jump MCMC provides a solution to this problem, allowing parameters pertaining to pairwise correlation between traits to be added or deleted over the course of the inference procedure. Haynie and Bowern (2016) utilize this method to analyze trajectories in the gain and loss of color terms over a posterior tree sample taken from Bowern and Atkinson 2012, finding broad support for Berlin and Kay (1969) theory that color terms are gained in a partially fixed order, but find also evidence for alternative pathways (e.g., green gained in the absence of red) and loss of color terms.

At the time of the publication of Dunn 2015, relatively little linguistic research had made use of ancestral state reconstruction and stochastic character mapping, a fact noted in the author, with some exceptions such as Maurits and Griffiths 2014, a cross-family study of the evolution of word order, in which the authors discuss how choice of prior affects ancestral state reconstruction at the root, among other issues. The state of affairs has changed since then.

A handful of papers have investigated character histories over tree topologies using parsimony. Jäger and List (2016) use parsimony to reconstruct root-meaning traits to Proto-Indo-European and Proto-Austronesian, and evaluate the performance of their model explicitly using metrics that are commonplace in other branches of computational linguistics. The authors adopt a strategy of comparing computationally reconstructed root-meaning traits against the “list of cognate classes present in the proto-language according to expert assessments” (p. 3). From this wording, it is not clear to me whether an automatically reconstructed root-meaning trait is evaluated against whether the expert considers the root to exist in the proto-language with a particular meaning, or whether the expert simply considers the root to exist in the proto-language. If the former, then all caveats regarding the difficulty of semantic reconstruction apply. If the latter, a root-meaning reconstruction could potentially be judged correct simply because the expert believes the root to be present in the proto-language, but with some other meaning. This could potentially lead to inflated recall values. Given the fact that the recall values reported by the authors never exceed 0.734 (still a relatively high number), it seems unlikely that the latter potentially anti-conservative metric was employed.

It is useful to have a baseline established via expert assessments against which to measure reconstruction methodologies. However, experts often disagree in their views; for this reason, it is as crucial to know how a qualitative conclusion was reached as it is how a quantitative result was achieved. It may be the case that a philologist carries out reconstruction effectively implements the parsimony algorithm by hand, i.e., by applying Occam’s Razor and attempting to avoid recurrent changes (though the scores that the authors report show that these processes are not 100% identical), such that this accuracy might not be particularly surprising. It is worth noting that the authors seem to treat this exercise as more of a sanity check — a step along the way to future work detecting large-scale recurrent change between branches.

Additionally, some research has carried out ancestral state reconstruction using likelihood-based methodologies. Recent work by Widmer et al. (2017) uses stochastic character mapping over a tree sample of Indo-European languages to demonstrate that at least one strategy of NP recursion was available over the course of Indo-European prehistory, according to their evolutionary model; ancestral state reconstruction is a crucial ingredient in this analysis.

Dunn et al. (2017) investigate the evolution of morphological case-marking on the subjects of a large number of verbs in Germanic languages; verbs in many Germanic languages co-occur with morphologically dative and accusative as well as nominative subjects. The study, based on a large amount of carefully curated data, is the first in linguistics to make use of (Pagel 1999) MULTISTATE methodology. Each verb is treated as a non-binary character with values denoting how the subject is marked, i.e., with N(ominative), A(ccusative), or D(ative) case. The authors use a handcrafted tree topology which they claim represents the consensus of historical linguists; however, a subgroup like Germanic, where internal language contact is well-documented, is an area where a tree sample, rather than a fixed tree, would be highly welcome, since there is arguably a degree of phylogenetic uncertainty. The authors compare an unconstrained model, where transitions between all three strategies occur at non-zero rates, against a set of models which set certain transition rates to zero in order to constrain the evolutionary pathway between states. They find that while a model representing “dative sickness” (where transitions A → D, A → N, and D → N are permitted, but not transitions D → A, N → D, or N → A) is not favored in terms of likelihood, it reconstructs ancestral state values which agree completely with qualitative judgments made by the authors and other specialists; they interpret this apparent success as evidence in favor of the dative sickness model.

If one wishes to see whether a particular model makes reconstructions which agree with a human expert (in order to confirm that such a model is appropriate to assume for orthogonal work, as do Jäger and List), that is indeed reasonable. However, the degree to which ancestral state reconstruction agrees with a philological reconstruction should have no bearing on model selection. This accuracy no way serves as independent evidence in favor of the dative sickness model, as Dunn et al. (2017) seem to imply; the authors have built their a priori assumptions regarding directionality of change in case assignment into the CTM model they use. When we look at reconstructions made by the unconstrained model (p. e15), we see that N subjects are erroneously reconstructed with high probability for a number of verbs, but for the dative sickness model, N subjects are reconstructed with zero probability; this appears to account for a large part of the model’s accuracy, as assessed by the authors. These facts should be unsurprising: in the dative sickness model, the rates for N → A and N → D change are set to zero. This means that N can be reconstructed with nonzero likelihood at an interior node of the tree only if N is present at all of the tips that descend from the node in question. In the dative sickness model, the D → A transition rate is also set to zero; hence, if D, A, and N are attested at the tips of the tree, A will be reconstructed at the root with 100% probability. If only D and A are present, the same is true. If only D and N are present, D is most likely to be reconstructed, but if the rates A → D and A → N are positive, there is still a small chance of A being reconstructed. And indeed, D is reconstructed only for verbs where there are no unambiguous instances of A at the tips of the tree; otherwise, A is reconstructed with 100% probability. It is pleasing when the result of a quantitative methodology dovetails with a philological analysis — evaluating the results of ancestral state reconstruction according to expert assessments provides a necessary sanity check for linguists — but in many situations, it behooves us to investigate exactly why this is the case. A simpler example is as follows: if, for instance, we wished to use a similar model to reconstruct word order for a family of languages exhibiting SOV and SVO word order, but chose to set the rate of SVO → SOV change to zero, we could expect SOV word order to be reconstructed with 100% probability. It may be that we believe SOV to be the ancestral word order, but this is not an observable fact against which to evaluate whether the model we have chosen is “correct.”

This discussion raises, it is hoped, some important questions regarding what it is exactly that we wish to do with CTM methods, ancestral state reconstruction, and related methodologies. It is important to recall that there is more than one way to allocate credibility to a particular hypothesis. The authors’ dative sickness model, which set certain transition rates to zero, is not favored in terms of likelihood. But it may be the case that transitions of the type D → A, N → A, and N → D occur with a probability that is greater than zero, but considerably lower than the transitions A → D, A → N, D → N. After all, the authors do provide some examples of dative to accusative change (p. e6), and it is therefore puzzling why they favor a model which fixes the probability of this type of change at zero. It could be the case that evidence can be adduced in favor of a weak formulation of the authors’ dative sickness model, one which predicts that D → A, N → A, and N → D changes happen at a less frequent rate than other changes. In theory, this could be done by running the unconstrained model for all of the verbs in the dataset, and then inspecting the rate estimates to see whether this pattern holds (the authors do not make the resulting MLEs available).

Figure 2:

Median values for log rates, aggregated across verbs.

I attempted to investigate this question using the data and tree from the paper, which the authors have made publicly available. I used MCMC to sample from the log posterior distributions of all six rates, constraining parameter space such that evolutionary rates did not exceed ≈ 2.7 mutations per 1000 years. The code used for inference can be found at https://github.com/chundrac/vanguard. In general, I found that this procedure did not always lead to convergence, and that chain acceptance rates often exceeded their recommended values (see above); this was undoubtedly an artifact of carrying out MCMC for a large number of rates over a relatively small tree, with a somewhat constrained state space. For this reason, we cannot necessarily expect to obtain the same rate values across different random number seeds. Median values for each rate, aggregated across verbs, are shown in the boxplot in Figure 2. The mean of the median log values of the “favored” rates (qAD, qDN, qAN) is slightly higher than that of the “disfavored” rates (−23.8 vs. −28.1), but this difference is not significant (Mann-Whitney U=698,p=.58). This is undoubtedly due to the fact that values for the rate qDA, included among the disfavored rates, are relatively high; while lower than that of qAD (−30.5 vs. −24.8), it is not significantly so. The rate qAN is also higher than qNA (−21.3 vs. −25.5), but this difference is insignificant as well. The same is true for qDN as opposed to qND (−25.2 vs. −28.2); like the others, this difference is insignificant. All in all, the rates for trajectories favored by the sickness model are higher than their counterparts in the opposite direction, but there is no significant evidence in favor of these asymmetries. Given the fact that other quantitative studies have found relatively strong support for the dative sickness hypothesis (del Prado Martín and Brendel 2016), evolutionary methods may not be the approach best suited to this question.