Challenges of sampling and how phylogenetic comparative methods help: with a case study of the Pama-Nyungan laminal contrast

2.1 Phylogenetic autocorrelation across the sciences

Different fields have their own lines of literature grappling with phylogenetic autocorrelation extending back many decades. In comparative anthropology, this issue was noted as early as 1889 by Sir Francis Galton in the context of cross-cultural datasets, which lack independence due to shared histories of cultural innovation and exchange between societies (Naroll 1961: 15). This phenomenon, known as Galton’s Problem, is now more precisely understood as a form of statistical autocorrelation, i.e., similarity between observations that correlates with their proximity, in this case, their proximity in evolutionary time. The same phenomenon has been recognised in comparative biology too. A seminal study concerning comparative studies of phenotypes, Felsenstein (1985) demonstrates that data from species cannot be assumed to be independently drawn from the same distribution, because species are related to one another via a branching, hierarchical phylogeny, thus, statistical methods that assume independent, identically-distributed observations will inflate the significance of the test (discussed further in Section 2.3 below). Linguists, it was argued, had been somewhat slower than those in other fields to acknowledge exposure to Galton’s problem, or phylogenetic autocorrelation (Perkins 1989: 293). However, this is a central concern of Dryer (1989: 259) and has been addressed in a considerable body of linguistic typological literature since then.

Statistical non-independence due to shared history is thus no new revelation, not in comparative anthropology, nor in comparative biology, nor in linguistic typology. However, there are many possible approaches to dealing with its challenges and a sizeable body of literature on the topic. As we will see, although precise strategies are varied, a notable commonality to all fields is a history of first attempting to address phylogenetic autocorrelation through the development of sampling methods for the creation of phylogenetically independent—or phylogenetically balanced—samples. The most striking differences between disciplines emerges only later, following the uptake in comparative biology of phylogenetic comparative methods.

2.2 Phylogenetic autocorrelation in linguistics

In linguistic typology, the use of phylogenetically balanced language samples remains the predominant way of accounting for phylogenetic autocorrelation and literature on this topic extends back several decades. Bell (1978: 145–149) argues that common strategies which simply ensure equally-weighted representation of “all major families” or all continents is inadequate due to differing rates of divergence among families. He estimates the number of language groups separated by more than 3,500 years of divergence and uses it as a heuristic for estimating genealogical biases in a selection of proposed language samples. He concludes that European languages tended to be overrepresented and Indo-Pacific languages underrepresented in typological language samples at his time of writing. He attributes this to a corresponding over/under-representation among quality language resources, which is a persistent problem for comparative linguistics. Perkins (1980, 1988 created a sample of 50 languages, later adapted by Bybee (1985), which attempts to account for both genealogical and areal biases by selecting no more than one language from each language phylum following Voegelin and Voegelin (1966) and no more than one language from each cultural and geographic area following work in comparative anthropology (Kenny 1975; Murdock 1967). This method attempts to account for non-independence due to areal spread, unlike Bell’s heuristic measure which accounts only for genealogical bias, however it does not account for differing ages of divergence and size of language phyla in the way Bell does.

Balanced sampling methods seek to produce linguistic samples that are independent, by selectively excluding the vast majority of attested languages, as necessitated by their extensive, inherent non-independence. As typologists have developed these methods, they have confronted two main complications.

The first complication is that it may be difficult to find criteria for the inclusion/exclusion of languages which truly remove all dependencies, or which are uncontroversial. Dryer (1989: 261) refers to the example of the inclusion of three languages in Perkins’ sample (Ingassana (also known as Gaam, gaam1241), Maasai (masa1300) and Songhai (koyr1242)) which potentially are related as part of the Nilo-Saharan family, and thus non-independent, although these relationships are remote and subject to debate. One aspect of this problem is that the maximal extent of presently established language families is partially a product of the extent of adequate documentation and scholarly attention, rather than a reflection of the fullest extent to which the family may be reconstructed (Levinson et al. 2011). Two languages which are presently understood to be unrelated, and therefore statistically independent, may in fact belong to a shared larger grouping, which has not yet been identified due to poor documentation or lack of historical-comparative study. A second aspect is that language families undoubtedly share deep-time relationships that are currently beyond the reach of the comparative method, even if all extant languages were documented and compared completely. Both challenges can lead to languages being deemed as independent when in reality they are not. Dryer (1989: 263) raises a related concern, which is that languages selected on the basis of genealogical independence may nonetheless share characteristics due to non-genealogical processes—language contact and borrowing. This motivates the use of areal criteria in addition to genealogical ones when constructing an independent sample. As Dryer (1989: 284) acknowledges however, linguistic areas may be also subject to the same concerns about undetected historical non-independence and it is possible that the whole world may, in effect, function as a single linguistic area, such that the distribution of certain linguistic features may reflect extremely remote areal or genealogical patterns rather than some true tendency of human language.

The second complication is that once all genealogical and areal criteria are adhered to, the resulting sample may be too small for use in statistical analysis (Cysouw 2005; Jaeger et al. 2011; Piantadosi and Gibson 2014). In response, linguists have proposed various procedures for constructing samples which, if not fully independent, at least have a high degree of independence. Dryer’s proposed solution is to build a sample of languages of approximately equal relative independence (at the level of major subfamilies within Indo-European, such as Romance, Germanic, and so on) for each of five large linguistic areas which are assumed to be independent, or at least sufficiently independent for statistical purposes. Any statistical test can then be applied to each of the five areas and only if the same result is replicated in all five areas is it considered statistically significant. If the same result is replicated in four of five areas, this falls short of statistical significance, although Dryer (1989: 272–273) considers such cases to be evidence of a “trend”. Nichols (1992: 41) uses Dryer’s area-by-area testing method as part of a three-pronged approach. For any given question, Nichols first conducts a chi-square test of the world sample and then re-tests the significance of the finding using either Dryer’s method or by running the same test on only the sample of “New World” languages (comprising North, Central and South America). Rijkhoff et al. (1993) and Rijkhoff and Bakker (1998) develop another approach to account for the possibility of non-independence across large linguistic areas and large, as-yet-undetected families. They permit multiple languages within a family to be included but develop a measure, based on the density of nodes in a known language phylogeny, to determine how many languages should be included. In this way, they also aim to account for the fact that some language families will have greater internal diversity than others (see also Bakker 2011; Miestamo et al. 2016).

Another approach is to include/exclude languages based on their typological profile. Following the logic that historical relatedness and interactions tend to result in elevated similarity, these methods bias their sample in favour of typological diversity, as a proxy for independence. Dryer and Haspelmath (2013) propose setting a minimum threshold of typological distance between languages, calculated from the World Atlas of Language Structures (WALS), such that languages must be sufficiently typologically distinct from others in the sample to warrant inclusion. Bickel (2009) develops an alternative algorithm based on Dryer (1989), which allows all uniquely-valued data points within a family to be included in the sample, but then reduces the weighting of data points in the final analysis where a particular value is over-represented within a family. In other words, if all the languages in a particular family share the same value for a typological variable of interest, those observations may be reduced to a single data point.

In these ways, developments in typological methodology have treated historical non-independence between languages as a challenge to be addressed through sampling. Earlier researchers sought to maintain the independence of their sample by maximising the genealogical distance between the languages in their sample, such that no two languages were known to belong to the same family. Later, with subsequent acknowledgement of the possibility of non-independence from very large language families, as well as large-scale areal diffusion and effects from as-yet undetected or unconfirmed historical relations, it became apparent that it may be impossible to create a sample which is simultaneously independent and sufficiently large to generate statistical significance. As discussed above, typologists have primarily responded to this dilemma by developing a variety of robustness checks, even bootstrapping-like processes, whereby languages are sampled at an approximately equal relative level of independence and the sample is then subdivided in some way and a statistical test replicated over each subdivision. More recent years have seen the continued evolution of statistics and robustness checking methods (for an overview, see Roberts 2018), although balanced sampling remains a common element of modern, large-scale comparative linguistic studies (for example, Blasi et al. 2017; Everett 2017; Everett et al. 2015).

Before turning to biology, it is worth underscoring how linguistic typology has arrived at its current mode of response to phylogenetic autocorrelation. The starting point is that many conventional statistical methods require observations that are independent, yet languages are non-independent. For four decades, the response has been to change the dataset, by means of balanced sampling, so that it better corresponds to the requirements of the statistics. Doing so requires excluding the vast majority of documented languages from the dataset and hence from the analysis, and even then, the result is still not truly independent. In the next section, we will see that biology initially followed the same path. The key breakthrough, though, was to invert the response to the original problem that phylogenetic autocorrelation posed: to change not the dataset to suit the statistics, but the statistics to suit the dataset. Those changed statistics are phylogenetic comparative methods.

2.3 Phylogenetic autocorrelation in comparative biology

Comparative biology faces the same issue of phylogenetic autocorrelation as comparative linguistics. Many conventional statistical methods assume that observations are independent, which is problematic since observations come from species, which are related to one another through shared evolutionary histories.

Earlier approaches to phylogenetic autocorrelation in biology are in a similar vein to the sampling methods in linguistic typology discussed in the previous section. Harvey and Mace (1982: 346–347) seek to find a taxonomic level to sample from, which strikes the right balance in terms of being sufficiently statistically independent without being so conservative that sample sizes become prohibitively small, an aim similar to Dryer (1989). Their proposed solution is to identify and sample from the lowest taxonomic level which can be “justified on statistical grounds”. One method of doing this is suggested by Clutton-Brock and Harvey (1977: 6–8), who conduct a nested analysis of variance and then select the taxonomic level containing the greatest level of variation. Similar to the methods of Dryer and Haspelmath (2013) and Bickel (2009), this approach makes reference to diversity in the traits of the species (cf. diversity in typological traits) to guide the sampling procedure.^[1]

As in linguistics, areality is also an issue in biology. Geographical and ecological proximity can lead to similarities in taxa (i.e., species or languages) which is causally separate from the effects of genealogy. Two distinct, causal scenarios can be distinguished. In the first scenario, material is passed directly between taxa, such as lateral transfer of genetic material between species, especially but not exclusively in prokaryotic life forms such as bacteria (Keeling and Palmer 2008), or borrowing between languages. In the second scenario there is no direct transfer of material, rather a shared environment leads to similar developments in taxa, such as parallel dwarfism on islands or, in some cases more contentiously, parallel conditioning of language by its environment (Blasi et al. 2017; Everett 2017, 2021; Everett et al. 2015). In both kinds of scenario, there is a causal, areally-correlated contribution to similarity which is separate from the contribution due to shared genealogy. While it is true that modern, genomic studies can circumvent some of the difficulties due to the second scenario in biology, it should be noted that phylogenetic comparative methods in biology predated the emergence of widespread genomic sequencing, and for many species including those attested only as fossils, genetic data is still not available. Consequently, the problem of convergent evolution due to areality was and still is a genuine, hard problem that comparative biology has faced, and should not be misunderstood as a problem specific to linguistics. In an approach with strong conceptual similarities to the area-by-area robustness checking of Dryer (1989) and Nichols (1992), Baker and Parker (1979: 85–86) discussed how the causal effects of ecological areas might be addressed while constructing a sample which is genealogically balanced. To do so, Baker and Parker (1979) replicate their analysis within individual families as well as within different ecological areas, with the assumption that if the same associations are observed within different areas as across the dataset as a whole, then one can discount the possibility that the full analysis is simply picking up differences between different families or different ecological areas.

In essence, both linguistics and biology face the same phenomenon of phylogenetic autocorrelation including the complication of areality, and for several decades explored strikingly similar methodological responses based on sampling (for further discussion, see Bromham 2017). However, in recent decades the primary methods in linguistic typology and biology have diverged as biology has undergone a fundamental shift. While typologists continue to focus on sampling procedures as the response to phylogenetic autocorrelation, comparative biologists have moved to a more direct, statistical solution. Since the solution addresses phylogenetic autocorrelation, not areality, our focus will narrow now to the genealogical aspects of taxon relatedness. We return to the separate and additional problem of areality in Section 7.2.

2.4 Phylogenetically independent contrasts

Felsenstein (1985) demonstrates that it is possible to account for phylogenetic non-independence in a statistical model without the need to remove data or compromise the unit of analysis (for example, by collapsing or averaging observations within a subgroup). Felsenstein’ breakthrough insight is that this can be achieved not by directly comparing non-independent observations but by comparing phylogenetically independent contrasts (PICs) between observations. His method has become, by one estimate, the most widespread in comparative biology (Nunn 2011: 162). The essential insight is relatively straightforward. Consider the tree in Figure 1. Any traits of A and B will be non-independent observations, since much of their evolutionary history is shared: all of the evolutionary change between points I and H, and between H and G, has contributed equally to both A and B. However, any differences (or in biological parlance, contrasts) between A and B have the particular status that they must have arisen after the split at point G. That period of development, after split G until the modern species (or languages) A and B is not shared with any other part of the tree. It is independent. Felsenstein’s insight is that by examining phylogenetic contrasts such as this, one can obtain observations that truly are independent. It is then possible to apply standard statistical tests to the phylogenetically independent contrasts (rather than directly to observed values) without phylogenetic autocorrelation introducing bias into the results.

Figure 1:

A phylogeny of six species or languages.

In the remainder of this subsection we discuss some finer technical points of Felsenstein’s notion for readers who are interested. Others may wish to skip ahead directly to the next subsection.

In order to calculate PICs not only between sister tips of a tree such as A and B, but also between sister interior nodes such as H and K, or node-tip sisters such as G and C, one requires in addition to a phylogeny, a model according to which the variable evolves. As a starting point, Felsenstein assumes a Brownian motion model of evolution, since Brownian motion is one of the simplest and most fundamental of all stochastic processes. In a Brownian motion model, an evolving quantitative trait can wander positively or negatively with equal probability, and each new time step is independent from the last, with the resulting effect that displacement of the variable over time will be drawn from a normal distribution with a mean of zero and variance proportional to the amount of elapsed time (Felsenstein 1985: 8). An observed contrast can be scaled by dividing it by the standard deviation of its expected variance. This gives a statistically independent contrast of expectation zero and unit variance (i.e. variance equal to 1). This process can be repeated for all adjacent tips in the tree. Contrasts can then be extracted from adjacent nodes in the tree, where the value of the node is an average of the observed values of the tips below it. In the end, there will be a collection of phylogenetically independent contrasts, all of expectation zero and unit variance, to which statistical analysis can be applied.

One drawback of Felsenstein’s initial method is the reliance on the assumption of Brownian motion as a model of variable evolution. Grafen (1989) subsequently devises a similar method, the phylogenetic regression, which has the flexibility to incorporate models of evolution other than Brownian motion. Further, Grafen’s method is able to be applied in situations where phylogenetic information is incomplete (for example, where the phylogeny is an incomplete work-in-progress rather than an accepted gold-standard). This method is a phylogenetic adaptation of generalised least squares (GLS). In this model, the value of a dependent variable, y _i , is predicted by the equation y _i  = α + βx _i  + ϵ, where α is the intercept, β is the regression slope, x is the independent variable and ϵ is an error term (Nunn 2011: 164). Phylogenetic information can be incorporated into the error term, in the form of a variance-covariance matrix of phylogenetic distances between tips in a tree. PICs and GLS are mathematically equivalent when a Brownian motion evolutionary model is assumed and the reference tree is fully bifurcated, so PICs are essentially a special case of GLS where these assumptions are met (Nunn 2011).

2.5 Phylogenetic comparative methods beyond biology

Linguistic typology and comparative anthropology have long faced the same essential problem of phylogenetic autocorrelation that comparative biology contends with. Initially, all three disciplines followed similar trajectories, responding to phylogenetic autocorrelation through the development of increasingly elaborate methods of balanced sampling. By historical accident it was in biology that the breakthrough of examining PICs occurred, but the breakthrough is a solution to an inherent problem that transcends disciplinary boundaries. Anthropologists, recognising the same problem in kind, followed this breakthrough in biology with their own uptake of phylogenetic comparative methods around 10–20 years later (e.g. Holden and Mace 2003, 2009; Jordan et al. 2009; Mace et al. 1994; Nunn 2011), and recently there has been growing interest in the application of phylogenetic comparative methods in linguistics (e.g. Bentz et al. 2018; Birchall 2015; Calude and Verkerk 2016; Cathcart et al. 2020; Dunn et al. 2011; Dunn et al. 2017; Jäger and Wahle 2021; Macklin-Cordes et al. 2021; Maurits and Griffiths 2014; Maslova 2000a, 2000b; Verkerk 2014, 2017; Zhou and Bowern 2015).

One of our motivations for this paper, however, is that despite the increasing uptake of phylogenetic comparative methods in linguistics, there has been little attempt until now to explain why phylogenetic comparative methods can best be understood as a continuation of a tradition of inquiry that typology is greatly invested in. Previously, that tradition of inquiry, whether in comparative biology, comparative anthropology or linguistics, had led to methods of balanced sampling. Like balanced sampling methods, phylogenetic comparative methods are a response to phylogenetic autocorrelation. Methodologists working on balanced sampling have striven to generate samples that come as close as possible to phylogenetic independence, but the goal cannot be fully attained even with the most elaborate sampling procedures, and in the meantime procedures of balanced sampling require the exclusion of the vast majority of documented languages from the dataset and hence from the analysis. As it turns out, the solution is to be found not in phylogenetically independent samples, but in phylogenetically independent contrasts (PICs). By focussing on PICs, Felsenstein unlocked a method for obtaining truly independent observations, without excluding data. This is why typologists have every reason to be keenly interested in phylogenetic comparative methods: they solve a problem which has stood at the centre of our discipline for decades.

In the sections that remain, we shift our focus away from theory and onto practicality: how can typologists begin making use of phylogenetic comparative methods? In Sections 3–5, we introduce key phylogenetic concepts and techniques that typologists can employ, followed by a phylogenetic typological case study in Section 6. In Section S1 of the Supplementary Materials, we provide an extended practical introduction to a suite of computational tools that have been designed with the typologist in mind (Round 2021a, 2021b), enabling phylogenetic comparative methods to be used in everyday typological research. In Section 7, we return to the topic of areality.

3.1 Phylogenetic signal in continuous variables

Blomberg et al. (2003) provide a suite of tools for quantifying phylogenetic signal, which have become somewhat of a standard in the field (cited 3,780 times as of September 2021, according to Google Scholar).^[3] Recent comparative studies using these tools include Balisi et al. (2018), Hutchinson et al. (2018) and Macklin-Cordes et al. (2021). Blomberg et al. (2003) present a descriptive statistic, K, which is generalisable across phylogenies of different sizes and shapes. In addition, they provide a randomisation test for checking whether the degree of phylogenetic signal for a given dataset is statistically significant. K can be calculated using either phylogenetically independent contrasts (PICs) (Felsenstein 1985) or generalised least squares (GLS) (Grafen 1989) (see Section 2.4). In a Brownian motion model, where variable values can wander up and down with equal probability through time, PIC variances are expected to be proportional to elapsed time. Among more closely related languages, where there has been less divergence time for variable values to wander, the variance of PICs is expected to be low. The randomisation test works by comparing whether observed PICs are lower than the PIC values obtained by randomly permuting the data across the tips of the tree. The process of permuting data across tree tips at random is repeated many times over. If the real variances, with data in their correct positions on the tree, are lower than 95% of the randomly permuted datasets, then the null hypothesis of no phylogenetic signal can be rejected at the conventional 95% confidence level. In other words, closely related languages resemble one another to a statistically significantly greater degree than would be expected by chance.

The descriptive statistic, K, quantifies the strength of phylogenetic signal. As with the randomisation procedure above, the input is a set of observed values, where each observation is associated with a tip of the reference tree. Blomberg et al. (2003: 722) give an explanation of the calculation of the K statistic. To summarise briefly, K is calculated by, firstly, taking the mean squared error (MSE₀), as measured from a phylogenetic mean,^[4] and dividing it by the mean squared error (MSE) calculated using a variance-covariance matrix of phylogenetic distances between tips in the reference tree (the same variance-covariance matrix of phylogenetic distances incorporated into the error term in GLS-based phylogenetic regression, as discussed in the previous section). This latter value, MSE, will be small when the pattern of covariance in the data matches what would be expected given the phylogenetic distances in the reference tree, leading to a high MSE₀/MSE ratio and vice versa. Thus, a high MSE₀/MSE ratio indicates higher phylogenetic signal. Finally, the observed ratio can be scaled according to its expectation under the assumption of Brownian motion evolution along the tree. This gives a K score which can be compared directly between analyses using different tree sizes and shapes. Where K = 1, this suggests a perfect match between the covariance observed in the data and what would be expected given the reference tree and the assumption of Brownian motion evolution. Where K < 1, close relatives in the tree bear less resemblance in the data than would be expected under the Brownian motion assumption. K > 1 is also possible—this occurs where there is less variance in the data than expected, given the Brownian motion assumption and divergence times suggested by the reference tree. In other words, close relatives bear greater resemblance than would be expected, given the overall phylogenetic diversity.

As discussed, the assumption of a Brownian motion model of evolution, where a variable is free to wander up or down in value, with equal probability, as time passes, is central to quantification of phylogenetic signal with the K statistic. Blomberg et al. (2003: 726–727) extend their approach to cover two different modes of evolution as well. This is achieved by incorporating extra parameters into the variance-covariance matrix to reflect different evolutionary processes. The first evolutionary model alternative is the Ornstein–Uhlenbeck (OU) model (Felsenstein 1988; Garland et al. 1993; Hansen and Martins 1996; Lavin et al. 2008) whereby variables are still free to wander up or down at random, but there is a central pulling force towards some optimum value. The second alternative is an acceleration–deceleration (ACDC) model, developed by Blomberg et al. (2003) where a variable value moves up or down with equal probability (like Brownian motion) but the rate of evolution will either accelerate or decelerate over time.

Other statistics for quantifying phylogenetic signal have been proposed and warrant mention. Freckleton et al. (2002) propose using the λ (lambda) statistic, based on earlier work by Pagel (1999). As for Blomberg et al. (2003) this approach works with a variance-covariance matrix showing the amount of shared evolutionary history between any two tips in the tree (the diagonal of the matrix, the variances, will indicate the total height of the tree; the off-diagonals, the covariances, will indicate the amount of shared evolutionary history between two given entities, before they diverge in the tree). The statistic, λ is a scaling parameter which can be applied to this variance-covariance matrix. Scaling the values in the matrix by λ transforms the branch lengths of the tree, from λ = 1, where branch lengths are left unscaled, to λ = 0, where all covariances in the matrix will be zero, in other words, no covariance through shared evolutionary history is indicated between any tips, thus all tips will be joined at the root by branches of equal length (a star phylogeny). Freckleton et al. (2002) present a method for finding the λ parameter that maximises the likelihood of a set of observations arising, given a Brownian motion model of evolution. If λ is close to 1, this indicates high phylogenetic signal, where the data closely fit expectation given the shared evolutionary histories in the tree and a Brownian motion model of evolution. Further measures which have been proposed are I (Moran 1950), a spatial autocorrelation measure which was adapted for phylogenetic analyses by Gittleman and Kot (1990), and C _mean (Abouheif 1999), which is a test for serial independence (for an overview, see Münkemüller et al. 2012). In an evaluation of different methods Münkemüller et al. (2012) find that, assuming a Brownian motion model of evolution, C _mean and λ generally outperform K and I. However, C _mean considers only the topology of the reference tree (i.e., the order of the branches from top to bottom), but not branch length information, and the value of the C _mean statistic is partially dependent on tree size and shape, so it lacks comparability between different studies. In addition, λ shows some unreliability with small sample sizes (trees with <20 tips).

3.2 Phylogenetic signal in binary variables

The methods so far described concern continuously-valued data. Other methods have been proposed for quantifying phylogenetic signal in binary and categorical variables too. Abouheif (1999) presents a simulation-based approach for testing whether discrete values along the tips of a phylogeny are distributed in a phylogenetically non-random way. Although this method is useful for testing whether the phylogenetic signal in a set of discretely-valued data is statistically significant, it does not provide a quantification of the level of phylogenetic signal which is comparable between different datasets. Although specific to binary data only, Fritz and Purvis (2010) present a statistic, D, which quantifies the strength of phylogenetic signal for binary variables.

The D statistic is based on the sum of differences between sister tips and sister clades, Σd. To summarise, following Fritz and Purvis (2010), differences between values at the tips of the tree are summed first (all tips will either share the same value, 0 or 1, with 0 difference; or one will be 0 and the other will be 1, for a difference of 0.5). Nodes immediately above the tips are valued as an average of the two tips below (either 0, 0.5 or 1) and the differences between sister nodes is summed. This process is repeated for all nodes in the tree, until a total sum of differences, Σd, is reached. At two extremes, data may be maximally clumped, such that all 1s are grouped together in the same clade in the tree and likewise for all 0s, or data may be maximally dispersed, such that no two sister tips share the same value (every pair of sisters contains a 1 and a 0, leading to a maximal sum of differences). Lying somewhere in between will be both a phylogenetically random distribution and a distribution that is clumped to a degree expected under a Brownian motion model of evolution. A distribution of sums of differences following a phylogenetically random pattern, Σd _r , is obtained by shuffling variable values among tree tips many times over. A distribution of sums of differences following a Brownian motion pattern, Σd _b is obtained by simulating the evolution of a continuous trait along the tree, following a Brownian motion process, many times over. Resulting values at the tips above a threshold are converted to 1, values below the threshold are converted to 0. The threshold is set to whatever level is required to obtain the same proportion of 1s and 0s as observed in the real data. Finally, D is determined by scaling the observed sum of differences to the means of the two reference distributions (the expected sums of differences under a phylogenetically random pattern and under a Brownian motion pattern).

Scaling D in this way provides a standardised statistic which can be compared between different sets of data, with trees of different sizes and shapes, as with K for continuous variables. One disadvantage of D, however, is that it requires quite large sample sizes (>50), below which it loses statistical power, increasing the chance of a false positive result (type I error).

Although we have restricted our focus to continuous and binary data here, some recent developments in testing for phylogenetic signal in other kinds of data warrant brief mention also. For example, Borges et al. (2019) have developed a statistic, δ, for quantifying phylogenetic signal in multivalued categorical variables. Other developments concern multivariate and multidimensional data. Zheng et al. (2009) present a multivariate version of the K statistic discussed in Section 3.1, for measuring phylogenetic signal in groups of related variables. Their statistic also incorporates measurement error. Finally, Adams (2014) presents K _mult , a statistic for detecting phylogenetic signal in multivariate traits, i.e. conceptually unitary evolutionary traits that are defined by multiple values (e.g. in biology, a set of measurements that together define skull shape).

In this section we have introduced the fundamental notion of phylogenetic signal—the degree to which the distribution of synchronic diversity reflects the shape of a phylogeny—and some key methods for estimating it. Of course, doing this requires a phylogeny to begin with, and typologists may have questions about the suitability of current linguistic trees for such purposes. It is to this important topic that we turn next.

6.1 Phylogenetic signal in the binary laminal contrast

In Figure 5b, we saw that the distribution of the laminal contrast in Pama-Nyungan hews closely to major subgroup boundaries, so we will not be surprised if a D test returns a strong confirmation of phylogenetic signal. Figure 6, which plots the presence and absence of a laminal contrast against the Pama-Nyungan tree, reinforces this expectation. We tested a set of 216 Pama-Nyungan languages (Round 2019), each coded for the binary presence/absence of the phonemic laminal contrast. To account for phylogenetic uncertainty, the statistic is calculated using 100 individual reference phylogenies.

Figure 6:

The distribution of the presence (light) and absence (dark) of a laminal place of articulation contrast across Pama-Nyungan, displayed on a maximum clade credibility (MCC) tree. An MCC tree is a single tree within a tree sample which most adequately represents the highest-probability subgroups in the trees of the sample. This MCC tree is taken from the sample of 100 highly-likely Paman-Nyungan phylogenies used in the current study.

The 100 results are summarised in Table 1. The mean D statistic obtained is low, at −0.439, indicating that the data is phylogenetically clumped to an even greater degree than expected under a Brownian motion model of evolution. Results like this can emerge when the variable under study has changed only rarely, and the changes have mostly been deep within the tree. This is the case in Pama-Nyungan, where variation in the presence/absence of the laminal contrast is mainly between major subgroups rather than within them. Returning to the statistical results, the hypothesis of randomness is rejected (p < 0.001) and the hypothesis of phylogenetic signal is not rejected (p = 0.987 ± 0.005). The values of the D statistic have a small standard deviation (0.019), indicating that a similar result is obtained for all 100 reference trees. In sum, the D test results confirm, in a quantitative manner and taking into account our uncertainty in the Pama-Nyungan phylogenetic tree, what our inspection of the map in Figure 5b could only suggest: that the binary presence/absence of the laminal contrast in Pama-Nyungan has strong phylogenetic signal.

6.2 Phylogenetic signal in continuously-valued phonotactic variables

Languages vary not only in what contrastive segments they have but also in how frequently they use them (Frisch et al. 2004; Hall 2009; Macklin-Cordes and Round 2020; Wedel et al. 2013). For example, Pitta Pitta (pitt1247, Blake 1990) and Burduna (burd1238, Burgman 2007) are similar in that they both contrast laminal stops, nasals and laterals in word-initial position. However, a closer examination reveals notable differences. In word-initial position before /u/, 29% of the consonantal laminals in Pitta Pitta are pre-palatal while 71% are dental, whereas in Burduna the frequencies are reversed, with 68% pre-palatal and just 32% dental. Frequency measures such as these can be viewed as continuous variables that can be investigated for phylogenetic signal (Macklin-Cordes et al. 2021). In this section we examine continuous variables of this kind, which describe the relative predominance of pre-palatals versus dentals in nine phonotactic positions, across 76 languages that possess the contrast. Data is from a phonemicised lexical database of Australian languages, which is under development (Round 2017), and which extends and enhances the Chirila database (Bowern 2016). Raw data tables and details of the primary language documentation sources are provided in Section S2 of the Supplementary Materials.

Our choice of nine variables is informed by the typological literature on Australian phonology. One long-established characteristic of Australian laminals is that their relative frequencies are sensitive to the quality of neighbouring vowels (Dixon 1970, 1980).^[7] Most Australian languages have three contrastive vowel qualities (Round 2022), with /i/ contexts favouring the laminal pre-palatal, /u/ contexts favouring the dental, and /a/ contexts somewhere in between. Here we examine the relative predominance of pre-palatals in word-initial position before /i, a, u/ and in intervocalic position before /i,a,u/ and after /i,a,u/.^[8] We apply the randomisation test described in Section 3.1 and then calculate a K statistic. As in our D test, we address phylogenetic uncertainty by applying the statistical tests using a sample of 100 reference trees.

Results are summarised in Table 2. The randomisation test finds phylogenetic signal to be statistically significant (p < 0.05) in all 9 variables and 100 reference trees except in two cases: these were the a_V and V_a contexts, for the same, one tree. Given that both contexts are judged to have significant phylogenetic signal in all other 99 trees in the 100-tree sample, we conclude that phylogenetic signal is present at a statistically significant level in all nine phonotactic variables.

The findings for the K statistic differ among the variables. For the word-initial variables, K is high, ranging from 0.783 to 1.322, whereas for the intervocalic variables it is uniformly lower, ranging from 0.337 to 0.696. In all cases, the standard deviation is low, indicating that similar results are obtained for all 100 reference trees. To put these K values in perspective, Blomberg et al. (2003) examined 121 biological traits of a wide variety of plant and animal organisms, finding mean K of 0.35 for behavioural traits, 0.54 for physiology and 0.83 for traits related to body size. Macklin-Cordes et al. (2021) estimated K for biphones (sequences of two adjacent phonemes) in Pama-Nyungan and found mean K of 0.52 for biphones of individual segments, and K of 0.63 when segments are binned into groups by place or manner of articulation. This suggests that our laminal phonotactic variables exhibit a level of phylogenetic signal at least as high as many evolved, biological traits, as well as the Pama-Nyungan biphone variables investigated in Macklin-Cordes et al. (2021).

The highest K value, at 1.322, is for laminals in word-initial position before /i/. A K value well above 1 is consistent with a scenario in which a linguistic property varies between deep branches of the tree, but much less so within the subgroups below those branches. This is true of Pama-Nyungan laminals word-initially before /i/. In the western half of the family, this position favours pre-palatals, reflecting a typical effect of the neighbouring vowel, whereas in the eastern half, the initial position in a word is one which favours dentals, irrespective of the following vowel.

A novel and consistent finding is that laminals exhibit stronger phylogenetic signal in word-initial position than intervocalically. There are many reasons why this might be so, and here we consider just one. Pertinacity (Dresher and Lahiri 2005) refers the perpetuation of linguistic patterns even as the items that instantiate them change. For instance, though a borrowed word may be new, its phonology is often reshaped to conform to the existing patterns in the recipient language (Hyman 1970), which then perpetuates the phonological patterns even as the set of items instantiating them changes. Similarly, if neologisms conform to existing statistical patterns in the lexicon, they too will contribute to pertinacity. Because our phonotactic variables are based on whole lexicons, and not merely a basic vocabulary list, lexical turnover will have been an important contributor to their historical dynamics. If it is the case that word-initial laminals have been subject to more-pertinacious changes than intervocalic laminals, such as more reshaping of borrowed words, or neologism which more closely replicates existing statistical patterns in the lexicon, then this could potentially lead to the difference in phylogenetic signal that we find. Whether there is additional evidence to support this hypothesis remains a question for future research, however the fact that such a hypothesis is able to emerge, illustrates how phylogenetic analysis can supplement the typologist’s existing toolkit for generating theoretically interesting hypotheses from the analysis of cross-linguistic data.

6.3 Genealogically-sensitive proportions of languages with a laminal contrast

We turn now to examine the phylogenetically-weighted proportion of Pama-Nyungan languages that have a laminal contrast. We know already, just by counting, that the simple proportion of Pama-Nyungan languages with a laminal contrast is 157/216 = 0.727. Our question here is, what is the proportion when genealogy is taken into account? As discussed in Section 5, there are different methods available for calculating this phylogenetic quantity, just as there are different kinds of non-phylogenetic averages. Here we compare the ACL and BM methods introduced earlier. We account for phylogenetic uncertainty by calculating them with respect to a sample of 100 reference trees. Table 3 reports the results. In this case the answer is broadly similar according to all three methods: the simple proportion is 0.727, the ACL-weighted proportion is somewhat higher, at 0.761 (SD 0.009) and the BM-weighted proportion marginally lower, at 0.705 (SD 0.003). The standard deviations of the phylogenetically weighted proportions are low, indicating that a similar result is obtained for all 100 reference trees. As mentioned in Section 5, an ACL proportion is more sensitive to genealogical structure deep within the tree than the BM method is, thus if we wish to remain conservative about our confidence in deep tree structure, we could conclude that a figure of around 71% (but perhaps as high as 76%) provides a good representation of the proportion of Pama-Nyungan languages that possess a laminal contrast. Note that unlike for balanced sampling, we did not need to discard any data, meaning that our results provide a faithful reflection of the evidence provided by all 216 languages and they do so while taking phylogenetic autocorrelation, including our uncertainty about Pama-Nyungan genealogy, into account.

Our case study has illustrated the application of methods and principles introduced in earlier sections. We have confirmed that the presence/absence of a laminal contrast in Pama-Nyungan has significant phylogenetic signal, notwithstanding a long history in the literature of emphasising its apparent areality. An examination of phylogenetic signal in continuously-valued phonotactic variables prompted us to notice a major east-west split in the treatment of word-initial laminals before /i/ and suggested a potential difference in the pertinacity of laminals and their statistical frequencies in word-initial versus intervocalic positions. Finally, having first confirmed the presence of phylogenetic signal, we then calculated genealogically-weighted proportions of the Pama-Nyungan languages which have the laminal contrast. This was done taking into account phylogenetic uncertainty in the Pama-Nyungan tree, and using two weighing methods which allow us to compare the consequences of investing a more conservative or less conservative degree of confidence in the deep-time branching structure of the trees.

7.1 Comparison across families and deep-time genealogy

Throughout this paper we have discussed phylogenetic comparative methods primarily within the scope of a single family. In this single-family, single-tree context we have examined phylogenetic uncertainty, testing for phylogenetic signal and the estimation of genealogically-sensitive averages and proportions. However, at the end of our discussion of uncertainty in phylogenetic trees in Section 4, we mentioned the problem of comparing across language families. We noted that logically, if it is believed that multiple families ultimately are related genealogically, then it is not possible to compare them without implicating a grand phylogeny that links them all. Methods which place all families on an equal footing merely do this by positing a rake tree. Thus, as radical as it may sound to say that we must hypothesise a deep-time tree which links currently-distinct families together, this is in fact something linguists have been doing for decades, covertly. Consequently, the question is not whether to use a grand, supra-familial tree but instead, which grand tree to use. Until now, linguists have generally declined to engage in positing grand trees that span beyond the reach of the comparative method, for the eminently good reason, that such trees cannot be demonstrated to be correct. However, as we have emphasised, trees do not need to be verifiably correct to be gainfully used with phylogenetic comparative methods. Instead, trees are hypotheses. Even if we do not, or cannot, know what the correct tree is, we surely can distinguish between more or less plausible hypotheses. Once we view the creation of grand trees as a matter of hypothesis generation, then there is every reason to begin working with them earnestly. For readers who find themselves still sceptical, consider the issue presented in the form of this question: Is a rake tree truly the best hypothesis that linguists could come up with about deep-time relatedness, entailing that every language family everywhere in the world is exactly equally related to every other? If our answer is anything other than an unequivocal yes, then we are effectively, tacitly entertaining the existence of other, more plausible grand trees.

To summarise so far, in order to apply phylogenetic comparative methods not only within but also across known families, we join the families in a grand tree. If the grand tree is a rake, then we are effectively continuing current practice in supra-familial language sampling. If the grand tree is otherwise, then we are beginning to explore alternative hypotheses for deep-time relatedness. As with the examples discussed earlier in the paper, phylogenetic comparative methods can be applied to multiple, alternative grand trees in order to reflect phylogenetic uncertainty and to investigate its implications.

Given this state of affairs, it strikes us that an important task for linguistic typology in coming years will be to establish an inventory of deep-time genealogical hypotheses, represented as phylogenies, as key ingredients for phylogenetic typological research, much in the way that the field in previous decades developed a variety of sampling techniques. Hypotheses within this inventory might come from many sources, whether from detailed interdisciplinary studies such as Matsumae et al. (2021) or novel linguistic attempts such as Jäger (2018), or more prosaically in the form of random samples of plausible hypotheses that meet certain constraining assumptions. There is ample scope for innovation. In Section S1 of the Supplementary Materials, we provide an extended description of a set of tools (Round 2021a) designed specifically with linguists in mind, for generating hypotheses about linguistic genealogy either within or across families, by creating and adjusting explicit linguistic phylogenies (see also Dediu 2018 for constructing within-family trees).

7.2 Areality

In scientific discussions with colleagues, we have encountered the concern that phylogenetic comparative methods cannot work, because they do not take into account the effects of areality (similarly, in published work see e.g. Blench 2015; François 2014). We believe that this concern may follow from a partial misapprehension about what phylogenetic comparative methods ought to be able to achieve. By way of comparison, it would be amiss to argue that a good model of gender should not be incorporated into a sociolinguistic analysis, merely because it does not account for geography. One could argue with good justification that we also desire an account of geography, but that is not the same thing as rejecting the successful model of gender. Similarly, we should not dismiss the breakthrough that Felsenstein achieved, dealing with genealogy far more effectively than in previous methods, merely because areality remains as difficult a problem as it always was. Here we briefly discuss why areality remains a hard problem and what can be done about it.

Viewed in mathematical and statistical terms, phylogenies are rather simple geometric objects. One consequence of their simplicity is that PICs can also be defined in a simple and effective manner. In contrast, the relationships implied by thousands of years of areality, including interactions with languages that have left no direct descendants, are significantly more complex. As mentioned in Section 4, comparative biology is also confronted with similarities shaped by areality, including in high-stakes fields such as bacteriology. Thus it is not for lack of motivation or interest that mathematical biologists are yet to produce methodological solutions to areality that match the solutions for phylogeny. The work is well underway, but the mathematics of historical networks, which such phenomena imply, is truly challenging (Elworth et al. 2019).

In this context, it is imperative for typologists to continue grappling with the problem of areality, though not by rejecting phylogenetic comparative methods, but instead by supplementing them. Recent methodological work that addresses areality in concert with phylogenetic comparative methods includes Cathcart et al. (2018) on areality in grammatical change, and Verkerk (2019) on estimating areality effects in relation to phylogenetic uncertainty. Similarly, it will be important to continue to learn more about the empirical facts of areality and its typological implications, to better understand its expected quantitative impact on the performance of phylogenetic comparative methods. For example, in the domain of lexical phylogenetic inference, Bowern et al. (2011) clarified empirical levels of lexical borrowing among hunter-gatherer and small-scale agriculturalist societies, providing crucial empirical knowledge about areality which could then be compared with the results of robustness studies (Greenhill et al. 2009), to suggest that at known empirical rates of borrowing, quantitative inference of phylogenies from lexical data should not suffer from significant impairment.

In all likelihood, areality will remain a tough challenge for linguistic typology, as it is for comparative biology, for some time to come. The problems that areality presents are different to and more complex than those of phylogeny. However, the mere fact that areality is hard is no sound reason to reject the advances offered by phylogenetic comparative methods. Instead, as always, the best available methods for handling genealogy must be supplemented with the current best attempts at handling areality.