Computational historical linguistics

3.1 Demonstration of genetic relationship

In (Jäger 2013), a dissimilarity measure between ASJP word lists is developed. Space does not permit to explain it in any detail here. Suffice it to say that this measure is based on the average string similarity between the corresponding elements of two word lists while controlling for the possibility of chance similarities. Let us call this dissimilarity measure between two word lists the PMI distance, since it makes crucial use of the pointwise mutual information (PMI) between phonetic strings.

To demonstrate that all Romance languages and dialects used in this study are mutually related, I will use the ASJP word lists from Papunesia, i.e. “all islands between Sumatra and the Americas, excluding islands off Australia and excluding Japan and islands to the North of it” (Hammarström et al. 2016) as training data and the ASJP word lists from Africa as test data.^[6] Input for inference are PMI distances between pairs of languages/dialect, and the output is the classification of this pair as related or unrelated, where two doculects count as related if they belong to the same language family according to the Glottolog classification. The distribution of PMI distances are shown in Figure 2. The graphics illustrates that all doculect pairs with a PMI distance ≤ 0.75 are, with a very high probability, related. The largest PMI distance among Romance dialects (between Aromanian and Nones) is 0.65.

Figure 2:

PMI distances between related and unrelated doculects from Papunesia, and between the Romance doculects.

A statistical test confirms this impression. I fitted a cumulative density estimation for the PMI distances of the unrelated doculect pairs from the training data, using the R-package logspline (Kooperberg 2016). If a pair of doculects has a PMI distance d, the value of the cumulative density function for d can then be interpreted as the (one-sided) p-values for the null hypothesis that the doculects are unrelated.

Using a threshold of α = 0.0001, I say that a doculect pair is predicted to be related if the model predicts it to be unrelated with a probability ≤ α. In Table 2, the predictions are tabulated against the Glottolog gold standard.

These results amount to ca. 0.3% of false positives and ca. 84% of false negatives. This test ensures that the chosen model and threshold is sufficiently conservative to keep the risk of wrongly assessing doculects to be related small. Since the method is so conservative, it produces a large amount of false negatives.

In the next step, I compute the probability of all pairs of Romance doculects to be unrelated, using the model obtained from the training data. Using the Holm–Bonferroni method to control for multiple tests, the highest p-value for the null hypothesis that the doculect pair in question is unrelated is 1.8 × 10^- 5, i.e. all adjusted p-values are < α. We can therefore reject the null hypothesis for all Romance doculect pairs.

3.2 Pairwise string comparison

All subsequent steps rely on a systematic comparison of words for the same concept from different doculects. Let us consider as an example the words for water from Catalan and Italian from ASJP, “aigua” and “acqua”. Both are descendants of the Latin “aqua” (Meyer-Lübke 1935: 46). In ASJP transcription, these are aixw~3 and akwa. The sequence w~ in the Catalan word encodes a diacritic (indicating labialization of the preceding segment) and is removed in the subsequent processing steps.

A pairwise sequence alignment of two strings arranges them in such a way that corresponding segments are aligned, possibly inserting gap symbols for segments in one string that have no correspondent in the other string. For the example, the historically correct alignment would arguably be as follows:

In this study, the quality of a pairwise alignment is quantified as its aggregate pointwise mutual information (PMI). (See List 2014 for a different approach.) The PMI between two sound classes a,b is defined as

Here s(a,b) is the probability that a is aligned to b in a correct alignment, and q(a),q(b) are the probabilities of occurrence of a and b in a string. If one of the two symbols is a gap, the PMI score is a gap penalty. I use affine gap penalties, i.e. the gap penalty is reduced if the gap is preceded by another gap. The captures the heuristics that both addition and deletion of phonetic material often targets several consecutive segments.

If the PMI scores for each pair of sound classes and the gap penalties are known, the best alignment between two strings (i.e. the alignment maximizing the aggregate PMI score) can efficiently be computed using the Needleman–Wunsch algorithm (Needleman and Wunsch 1970).

The quantities s(a,b) and q(a),q(b) must be estimated from the data. Here I follow a simplified version of the parameter estimation technique from (Jäger 2013). In a first step, I set

Also, I set the initial gap penalties to - 1. (This amounts to Levenshtein alignment.) For instance, the alignment for the example given above would come out as

This alignment has one identity (the initial a) and three non-identities, so the PMI₀ score is - 3.

Using these parameters, all pairs of words for the same concept from different doculects are aligned.

From those alignments, s(a,b) is estimated as the relative frequency of a and b being aligned among all non-gap alignment pairs, while q(a) is estimated as the relative frequency of sound class a in the data. The PMI scores are then estimated using equation (1). For the gap penalties I used the values from (Jäger 2013), i.e. - 2.49 for opening gaps and - 1.70 for extending gaps (i.e. for gaps preceded by another gap). Using those parameters, all synonymous word pairs are realigned.

In the next step, only word pairs with an aggregate PMI score ≥ 4.45 are used. (This threshold is taken from Jäger 2013 as well.) Those word pairs are realigned and the PMI scores are reestimated. This step is repeated ten times.

The threshold of 4.45 is rather strict; almost all word pairs above this threshold are either cognates or loans. For instance, for the language pair Italian/Albanian, the only translation pair with a higher PMI score is Italian peSe/Albanian peSk (“fish”), where the former is a descendant and the latter a loan from Latin piscis (cf. http://ielex.mpi.nl). For Spanish/Romanian, two rather divergent Romance languages, we find eight such word pairs. They are shown alongside with the inferred alignments in Table 3.

The aggregate PMI score for the best alignment between two strings is a measure for the degree of similarity between the strings. I will call it the PMI similarity henceforth.

3.3 Cognate clustering

Automatic cognate detection is an area of active investigation in CHL (Dolgopolsky 1986; Bergsma and Kondrak 2007; Hall and Klein 2010; Turchin et al. 2010; Hauer and Kondrak 2011; List 2012; List 2014; Rama 2015; Jäger and Sofroniev 2016; Jäger et al. 2017, inter alia). For the present study, I chose a rather simple approach based on unsupervised learning.

Figure 3:

PMI similarities for synonymous and non-synonymous word pairs.

Figure 3 shows the PMI similarities for words from different doculects having different or identical meanings. Within our data, synonymous word pairs are, on average, more similar to each other than non-synonymous ones. The most plausible explanation for this effect is that the synonymous word pairs contain a large proportion of cognate pairs. Therefore “identity of meaning” will be used as a proxy for “being cognate”. The implicit assumption underlying this procedure is that cognate words always have the same meaning. This is evidently false when considering the entire lexicon. There is a plethora of examples, such as English deer vs. German Tier “animal”, which are cognate cf. (Kroonen 2013: 94) without being synonyms. However, within the 40-concept core vocabulary space covered by ASJP, such cross-concept cognate pairs are arguably very rare.

For each concept, a weighted graph is constructed, with the words denoting this concept as vertices. Two vertices are connected if the predicted probability of these words to be synonymous (based on their PMI similarity and the logistic regression model) is ≥ 0.5. The weight of each edge equals the predicted probabilities. The nodes of the graph are clustered using the weighted version of the Label Propagation algorithm (Raghavan et al. 2007) as implemented in the igraph software (Csardi and Nepusz 2006). As a result, a class label is assigned to each word. Non-synonymous words never carry the same class label. Table 4 illustrates the resulting clustering for the concept “person” and a subset of the doculects. A manual inspection reveals that the automatic classification does not completely coincide with the cognate classification a human expert would assume. For instance, the descendants of Latin homo are split into classes 1, 2, 5, and 7. Also, Gheg Albanian 5eri and Sardinian omini have the same label but are not cognate.

Based on evaluations against manually assembled cognacy judgments for different but similar data (Jäger and Sofroniev 2016; Jäger et al. 2017), we can expect an average F-score of 60%–80% for automatic cognate detection. This means that on average, for each word, 60%–80% of its true cognates are assigned the same label, and 60%–80% of the words carrying the same label are genuine cognates.

3.4 Phylogenetic inference

3.4.2 Application to the case study

For the case study, doculects were represented by two types of binary characters:

1. Inferred class label characters (cf. Subsection 3.3). Each inferred class label is a character. A doculect has value 1 for such a character if and only if its word list contains a word carrying this label.

2. Soundclass-concept characters. There is one character for each pair (s,c) of a sound class s and a concept c. A doculect has value 1 for that character if and only if its word list contains a word w for c that contains s in its transcription.

^[11]

Both types of characters carry a diachronic signal. For instance, the mutation 0→1 for class label 6/concept person (cf. Table 4) represents a lexical replacement of Latin “homo” or “persona” by descendants of Latin “christianus” in some Romance dialects (Meyer-Lübke 1935: 179). The mutation 0→1 for the soundclass-concept character k/person represents the same historical process. Soundclass-concept characters, however, also capture sound shifts. For instance, the mutation 0→1 for b/person reflects the epenthetic insertion of b in descendants of Latin “homo” in some Iberian dialects.^[12]

I performed Bayesian phylogenetic inference on those characters. The inference was carried out using the Software MrBayes (Ronquist and Huelsenbeck 2003). Separate rate models were inferred for the two character types. Rate variation across characters was modeled by a discretized Gamma distribution using four rate categories. I assumed no rate variation across edges. Root probabilities were identified with equilibrium probabilities. An ascertainment correction for missing all-0 characters was used.

I assumed rates to be constant across branches. This entails that the fitted branch lengths reflect the expected amount of change (i.e. the expected number of mutations) along that branch.

In such a model, the likelihood of a phylogeny does not depend on the location of the root (the assumed Markov process is time reversible.) Therefore, phylogenetic inference provides no information about the location of the root. This motivates the inclusion of the Albanian doculects. Those doculects were used as outgroup, i.e. the root was placed on the branch separating the Albanian and the Romance doculects.

I obtained a sample of the posterior distribution containing 2,000 phylogenies. Figure 4 displays a representative member of this sample (the maximum clade credibility tree). The labels at the nodes indicate posterior probabilities of that node, i.e. the proportion of the phylogenies in the posterior sample having the same sub-group.

Figure 4:

Representative phylogeny from the posterior distribution. Labels at the internal nodes indicate posterior probabilities.

These posterior probabilities are mostly rather low, indicating a large degree of topological variation in the posterior sample. Some subgroups, such as Balkan Romance or the Piemontese dialects, achieve high posterior probabilities though.

Notably, branch lengths carry information about the amount of change. According to the phylogeny in Figure 4, for instance, the Tuscan dialect of Italian (ITALIAN_GROSSETO_TUSCAN) is predicted to be the most conservative Romance dialect (since its distance to the latest common ancestors of all Romance dialects is shortest) and French the most innovative one.

These results indicate that the data only contain a weak tree-like signal. This is unsurprising since the Romance languages and dialects form a dialect continuum where horizontal transfer of innovations is an important factor.

Phylogenetic trees, like traditional family trees, only model vertical descent, not horizontal diffusion. They are therefore only an approximation of the historical truth. But nevertheless, they are useful as statistical models for further inference steps.

3.5 Ancestral state reconstruction

If a fully specified model is given, it is possible to estimate the probability distributions over character states for each internal node.

Let M=〈T,θ⃗〉 be a model, i.e. a phylogeny T plus further parameters θ⃗ (rates and root probabilities, possibly specifications of rate variation). Let i be a character and n a node within T.

The parameters θ⃗ specify a Markov process, including rates, for the branch leading to n. Let 〈π0,π1〉 be the equilibrium probabilities of that process. (If n is the root, 〈π0,π1〉 are directly given by θ⃗.)

Let M(n_i = x) be the same model as M, except that the value of character i at node n is fixed to the value x. L(M) is the likelihood of model M given the observed character vectors for the leaves.

The probability distribution over values of character i at node n, given M, is determined by Bayes’ rule:

Figure 5:

Ancestral state reconstruction for character person:3.

Figure 5 illustrates this principle with the Romance part of the tree from Figure 4 and the character person:3 (cf. Table 4). The pie charts at the nodes display the probability distribution for that node, where white represents 0 and red 1.

This kind of computation was carried out for each class label character and each tree in the posterior sample for the latest common ancestor of the Romance doculects. For each concept, the class label for that concept with the highest average probability for value 1 at the root of the Romance subtree was inferred to represent the cognate class of the Proto-Romance word for that concept.^[13] For the concept person, e.g. character person:3 (representing the descendants of Latin “persona”), comes out as the best reconstruction.

3.6 Multiple sequence alignment

In the previous step, for the concept eye, the class label 6 was reconstructed for Proto-Romance. Its reflexes are given in Table 5.

A multiple sequence alignment (MSA) is a generalization of pairwise alignment to more than two sequences. Ideally, all segments within a column are descendants of the same sound in some common ancestor.

MSA, as applied to DNA or protein sequences, is a major research topic in bioinformatics. The techniques developed in this field are mutatis mutandis also applicable to MSA of phonetic strings. In this subsection, one approach will briefly be sketched. For a wider discussion and proposals for related but different approaches, see (List 2014).

Here I will follow the overall approach from (Notredame et al. 2000) and combine it with the PMI-based method for pairwise alignment described in Subsection 3.2. (Notredame et al. 2000) dub their approach T-Coffee (“Tree-based Consistency Objective Function For alignment Evaluation”), and I will use this name for the method sketched here as well.

In a first step, all pairwise alignments between words from the list to be multiply aligned are collected. For this purpose, I use PMI pairwise alignment. Some examples would be

The last row shows the score of the alignment, i.e. the proportion of identical matches (disregarding gaps).

In a second step, all indirect alignments between a given word pair are collected, which are obtained via relation composition with a third word. Some examples for indirect alignments between okiu and oky would be:

The direct pairwise alignment matches the i in okiu with the y in oky. Most indirect alignments pair these two positions as well, but not all of them. In the last columns, the i from okiu is related to the k of oky, and the y from oky with a gap. For each pair of positions in two strings, the relative frequency of them being indirectly aligned, weighted by the score of the two pairwise alignments relating them, is summed. They form the extended score between those positions.

The optimal MSA for the entire group of words is the one where the sum of the pairwise extended scores per column is maximized. Finding this global optimum is computationally not feasible though, since the complexity of this task grows exponentially with the number of sequences. Progressive alignment (Hogeweg and Hesper 1984) is a method to obtain possibly sub-optimal but good MSAs in polynomial time. Using a guide tree with sequences at the leaves, MSAs are obtained recursively leaves-to-root. For each internal node, the MSAs at the daughter nodes are combined via the Needleman–Wunsch algorithm while respecting all partial alignments from the daughter nodes.

For the words from Table 5, this method produces the MSA in Table 6. The tree in Figure 4, pruned to the doculects represented in the word lists, was used as guide tree.

Using this method, MSAs were computed for each inferred class label that was inferred to be present in Proto-Romance via Ancestral State Reconstruction.

3.7 Proto-form reconstruction

A final step toward the reconstruction of Proto-Romance forms, Ancestral State Reconstruction is performed for the sound classes in each column, for each MSA obtained in the previous step.

Consider the first column of the MSA in Table 6. It contains three possible states, v, w, and the gap symbols -. For each of these states, a binary presence–absence character is constructed. For doculects which do not occur in the MSA in question, this character is undefined.

The method for ancestral state reconstruction described in Subsection 3.5 was applied to these characters. For phylogeny in the posterior sample, the probabilities for state 1 at the Proto-Romance node was computed for each character. For each column of an MSA, the state with the highest average probability was considered as reconstructed.

The reconstructed proto-form for a given concept is then obtained by concatenating the reconstructed states for the corresponding MSA and deleting all gap symbols. The results are given in Table 7.

3.8 Evaluation

To evaluate the quality of the automatic reconstructions, they were compared to the corresponding elements of the Latin word list. For each reconstructed word, the normalized Levenshtein distance (i.e. the Levenshtein distance divided by the length of the longer string) to each Latin word (without diacritics) for that concept was computed. The smallest such value counts as the score for that concept. The average score was 0.484. The extant Romance doculects have an average score of 0.627. The most conservative doculect, Sardinian, has a score of 0.502, and the least conservative, Arpitan, 0.742. The evaluation results are depicted in Figure 6.

These findings indicate that the automatic reconstruction does in fact capture a historical signal. Manual inspection of the reconstructed word list reveals that, to a large degree, the discrepancies to Latin actually reflect language change between Classical Latin and the latest common ancestor of the modern Romance doculects, namely Vulgar Latin.^[14] To mention just a few points: (1) Modern Romance nouns are mostly derived from the (Vulgar-)Latin accusative form (Herman 2000: 3), while the word lists contains the nominative form. For instance, the common ancestor forms for “tooth” and “night” are dentem and noctem. The reconstructed t in the corresponding reconstructed forms are therefore historically correct. (2) Some Vulgar Latin words are morphologically derived from their Classical Latin counterparts, such as mons→montanea “mountain” (Meyer-Lübke 1935: 464) or genus→genukulum “knee” (Meyer-Lübke 1935: 319). This is likewise partially reflected in the reconstructions. (3) For some concepts, lexical replacement by non-cognate words took place between Classical and Vulgar Latin, such as via→strata “path”,^[15]ignis→focus “fire” (Meyer-Lübke 1935 293), or iecur→ficatum “liver” (Herman 2000: 106). Again, this is reflected in the reconstruction.

On the negative side, the reconstructions occasionally reflect sound changes that only took place in the Western Romania, such as the voicing of plosives between vowels (Herman 2000: 46), as in figat “liver” (←ficatum).

Let me conclude this section with some reflections on how the reconstructions were obtained and how this relates to the comparative method.

A major difference to the traditional approach is the stochastic nature of the workflow sketched here. Both phylogenetic inference and ancestral state reconstruction are based on probabilities rather than categorical decisions. The results shown in Table 7 propose a unique reconstruction for each concept, but it would be a minor modification of the workflow only to derive a probability distribution over reconstructions instead. This probabilistic approach is arguably an advantage since it allows to utilize uncertain and inconclusive information while taking this uncertainty properly into account.

Figure 6:

Average normalized Levenshtein distance to Latin words: reconstruction (dashed line) and extant Romance doculects (white bars).

Another major difference concerns the multiple independence assumptions implicit in the probabilistic model sketched in Subsection 3.4. The likelihood of a phylogeny is the product of its likelihoods for the individual characters. This amounts to the assumptions that the characters are mutually stochastically independent.

The characters used here (and generally in computational phylogenetics as applied to historical linguistics) are mutually dependent in manifold ways though. For instance, the loss of a cognate class makes it more likely that the affected lineage will acquire another cognate class for the same semantic slot and vice versa.

This problem is even more severe for phonetic change. Since the work of the Neogrammarians in the nineteenth century, it is recognized that many sound changes are regular, i.e. they apply to all instances of a certain sound (given contextual conditions) throughout the lexicon. Furthermore, both regular and irregular sound changes are usually dependent on their syntagmatic phonetic context, and sometimes on the paradigmatic context within inflectional paradigms as well. Bouchard-Côté et al. (2013) and Hruschka et al. (2015) propose more sophisticated probabilistic models of language change than the one used here to take these dependencies into account.^[16]

Last but not least, the treatment of borrowing (and language contact in general) is an unsolved problem for computational historical linguistics. Automatic cognate clustering does not distinguish between genuine cognates (related via unbroken chains of vertical descent) and (descendants of) loanwords. This introduces a potential bias for phylogenetic inference and ancestral state reconstruction, since borrowed items might be misconstrued as shared retentions.