2.1. Empirical linguistic data

The number of tones for each language was obtained from the PHOIBLE database (Moran & McCloy, Reference Moran and McCloy2019). Missing data were filled in from the World Phonotactics Database (Donohue et al., Reference Donohue, Hetherington, McElvenny and Dawson.2013). A random sample of 20 other languages was obtained manually from language descriptions (see the Supporting Information). The number of vowels and consonants in the phoneme inventory of each language was also obtained from Moran and McCloy (Reference Moran and McCloy2019). The vowel ratio in basic vocabulary was estimated from data in the Automated Similarity Judgement Program (ASJP version 20; CLLD, 2022) using a method similar to Everett (Reference Everett2017), except that we use the proportion of vowels used in the basic vocabulary compared with the total number of segments. Where there were multiple sources available for a language, the mean value was used.

2.2. Estimated ancestral states

The geographic–phylogenetic tree of Bantu was taken from Grollemund et al. (Reference Grollemund, Branford, Bostoen, Meade, Venditti and Pagel2015). This is a Bayesian phylogenetic reconstruction of the Bantu subfamily of languages which provides estimated ages and geographic locations for each ancestral node. These estimates were only available for the maximum clade credibility tree. Following the documentation, a reference year of 2000 CE was used.

Ancestral states for language data were estimated using Random-Walk Markov Chain Monte-Carlo Bayesian ancestral state reconstruction in BayesTraits (Meade & Pagel, Reference Meade and Pagel2022). A two-step procedure was used, first estimating a distribution of models from existing data, then estimating unknown values. Although the number of tones is a discrete measure, we used continuous ancestral state reconstruction since there were too many possible unique states for reliable estimate of a discrete state model. Proto-Bantu was fossilised as having two tones, based on traditional reconstructions from historical linguists (Nurse & Philippson, Reference Nurse and Philippson2006; Stevick, Reference Stevick1969). A uniform prior from zero to 12 tones was used. The Bayesian process ran for a million iterations of burn-in and a million more iterations of sampling.

The same process was used to estimate the number of vowels and consonant segments in the phoneme inventory, and the vowel ratio in basic vocabulary. The prior biases were set according to the empirical distribution of all sources in PHOIBLE: a gamma distribution for the number of vowels (k = 3.3, θ = 0.3) and number of consonants; and a normal distribution for the vowel ratio (mean = 0.4659, SD = 0.0556). While the vowel ratio is technically bounded, the empirical distribution is far from either 0 or 1, making the bounds unreachable in practice. No nodes were fossilised for these alternative linguistic measures.

These processes resulted in a distribution of estimations for the linguistic measures at each node of the geo-phylo tree.

2.3. Humidity data

Humidity estimates for Sub-Saharan Africa for the last 10,000 years were taken from geophysical simulations (Davies-Barnard et al., Reference Davies-Barnard, Ridgwell, Singarayer and Valdes2017; Valdes et al., Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix and Williams2017). These estimate the Stephenson screen measurement of specific humidity (1.5 m above ground) for regularly spaced points in time and space across the world. A reference year of 1850 CE was used. To obtain estimates for specific points in time and space, the data were interpolated using a general additive model that predicted humidity given the year, latitude and longitude. A knot was specified for each unique year, latitude and longitude, and latitude and longitude had interactions by year. The model also included a tensor product smooth for all three variables. This effectively allows smooth interpolation between empirical observations in a three-dimensional space. In informal testing, the fit to the data was better for the general additive model (R ² = 0.81) than bilinear interpolation (R ² = 0.74).

2.4. Summary of data

The steps above resulted in a geo-phylo tree, where each leaf and ancestral node had estimates for their age, geographic location, humidity and the number of tones (and the alternative linguistic measures of vowels). This comprised data for 210 leaf languages for number of tones, 161 leaf languages for number of vowels and consonants, and 383 leaf languages for the ratio of vowels in basic vocabulary.

Tones had a weak phylogenetic signal (Pagel's λ = 0.35), while humidity had a very strong phylogenetic signal (Pagel's λ = 1.00, although we are not suggesting that humidity is inherited culturally). The number of vowels had a strong phylogenetic signal (Pagel's λ = 0.98), while the number of consonants had a slightly weaker phylogentic signal (Pagel's λ = 0.69). The ASJP vowel ratio had a very strong phylogenetic signal (Pagel's λ = 1.00)

All raw data, processing scripts and results are available in an online repository (https://github.com/seannyD/TonesClimateGeoPhylo_Public).

3.1. The null model of environmental effects

As a starting point of this investigation, we simulate the evolution in Bantu tones according to a null model with random changes. To do this, we iterate over the tree nodes and apply the following function:

$$T_{desc} = T_{anc} + \epsilon B$$

Here, the tones at the descendant node T _desc are equal to the ancestor node T _anc plus the product of a normally distributed error $\epsilon$ with standard deviation of 1 and the branch length B. The humidity information for every node is retained from the original dataset. According to this process, we simulate 100 random datasets and for each of these, we applied a simple Frequentist univariate regression to the leaf nodes (i.e. the extant languages) with tones as the outcome and humidity as the predictor variable. Thus, for each of these datasets, we obtain a slope of the predictor humidity along with a significance value at the level α = 0.05. In a dataset without significant bias, we would expect the type-I error rate for the predictor variable slope to be ~0.05, i.e. the slope in around five out of 100 datasets would be incorrectly significant.

The result of this simulation shows that the calculated slopes range widely from strongly negative to strongly positive slopes. Figure 2 shows a histogram of those predictor slopes.

Figure 2. Histogram of the distribution of humidity predictor slopes for the simulated datasets under the null model.

Recall that the basis of those datasets is a null model in which humidity does not have any effect. Nevertheless we find that in 74 out of 100 datasets the univariate linear model detects a significant slope (either positive or negative) which is an extraordinarily high type-I error rate. This is indicative of a severe bias in the data generation process associated with the genealogical tree. However, what is it about this process that makes it exceptionally susceptible to these errors?

To explore this issue further, we need to investigate the original dataset specifically regarding the distribution of humidity and tones along the structure of the phylogenetic model: Figure 3 is a visualisation of the diachronic development of humidity and tones along the nodes of the Bantu family tree. Here, each line segment represents a link between an ancestor node and a descendant node plotted against their respective humidity or tone values and the respective node age. The thickness of the lines is proportional to the number of terminal nodes that descend from this link.

Figure 3. Evolution of humidity (left) and tones (right) along the Bantu language family tree.

Especially salient is that in the plot of the humidity development, we find that over time, two clusters seem to form, one with high humidity values (> 0.016) and one with lower humidity values (< 0.016), the latter of which exhibits more variance. This clustering arises if there is one branch early in the tree that moves into more or less humid regions. We can demonstrate this phenomenon when we select two clusters and plot them separately. In Figure 4, we selected two representative clusters for each of the variables by selecting a node (represented by black dots) and plotting all descendant nodes from this point. Note that for this visualisation, different nodes were selected for tones and humidity.

Figure 4. Two examples of clustering of decedents of two nodes in the phylomorpho space.

In Figure 4, we can see that the two clusters differ considerably in the distribution of their terminal nodes. The line segments on the left of the graphs indicate the 80% highest density credible intervals for the distribution of the extant languages regarding humidity and tones. This means that, in the humidity graph, 80% of the languages descended from the blue lineage fall into a narrow range between approximately 0.0118 and 0.0142 humidity. A similar distribution can be observed for tones. This, however, is not necessarily problematic; to reveal the biasing effect of this distribution, we need to consider two major points.

Firstly, these clusters can arise owing to small changes early in the tree. When we consider the humidity graph in Figure 4, we see that the difference between the blue and red cluster is due to the fact that the ancestral community strongly and rapidly changed their environment between the most recent common ancestor of the red and blue lineage at about 3,500 BP and the ancestor node at 2,500 BP. After that, there are further no strong shifts. In effect, most of the contrast between the groups is due to a sudden change in environment by the ancestor node of the blue cluster. The environment in the history of this cluster is approximately homeostatic before 3,500 and after 2,500 BP. A similar, albeit weaker effect can be found for tones as well (see Figure 4).

The second driver of bias which builds upon this comes into play when the genealogically highly clustered tones and humidity interact. Specifically when clusters align or partially align (i.e. many tones – high humidity, or many tones – low humidity and vice versa), strong positive and negative correlations are easily created. Further, owing to the cluster contrasts resulting from changes early in the tree, smaller variations early on can result in strong founder effects for subsequent nodes. This is the reason why we find alternatingly strong positive or negative effects, despite the simulation datasets being generated by a null model.

3.2. Node importance bias

To define even more precisely the driving factor of this process, we can investigate node importances, specifically, which nodes have the strongest impact on the outcome. Put differently: in which nodes are small variations likely to cause a bias in the synchronic data?

To calculate which nodes are most influential in shaping the outcome, we create a new dataset for each node in the tree where a fixed value is added to the tone values of all descendants of this node (in this case, we add 1). This means that we create a large number of datasets that each differ in the tone values of one lineage. In a second step, we take the predictor slope of humidity on tones in the new dataset and calculate the absolute distance to the predictor slope in the original, unaltered, dataset. Subsequently, we obtain for every dataset a value of how strongly changing the tone value of a certain node changes the calculated predictor slope. This results in a per-node value of importance since nodes that are more influential cause a greater shift in the regression slope. Figure 5 shows the Bantu family tree with every node coloured according to its importance.

Figure 5. Bantu family tree with node importance ranging from blue (less important) to red (more important).

We can observe that the most important nodes are those that (1) are located somewhat early in the tree and that (2) have many, closely clustered descendants. For this tree, eight out of 499 nodes can be viewed as highly important which means that around 1.6% of the observations in the dataset determine whether we see a positive/negative/no correlation between humidity and tones synchronically. we expect that such biases exist in other language families and across other phonological features. This is a major problem for correlational studies on synchronic data, especially for studies which only control at the level of language family rather than a full phylogeny, since this example demonstrates that intra-family biases can lead to distorted findings.

4.1. Accounting for genealogical bias

Given the aforementioned issues of using synchronic data that are the result of a genealogical process, we aim for a diachronic model to remove these biases. Rather than modelling the synchronic relationship between humidity and tones, we model the difference in tones between an ancestor node and a descendant node.

If humidity indeed had an influence on the diachronic development of tones, we would expect to see that languages in more humid environments tend to develop more tones. In other words, if this conjecture is true, then we would see that languages located in more humid climates would tend to acquire more tones over time than languages in more arid environments would. This focus on changes between nodes removes the genealogical bias since the clusters that form in synchronic data are not present.

4.2. Causal relationships

There are multiple causal relationships associated with humidity and tones. Here, we use tools from causal inference to map out these causal relationships in order to arrive at the form of our statistical model. Figure 6a shows some complex relationships between hypothetical cultures and climates as a causal network (or ‘Bayesian network’; see Pearl, Reference Pearl1985). Imagine a group of cultures that disperses over time (T0, T1, T2) and moves across a geographic space (locations L0, L1, L2). We are interested in estimating the effect of humidity X on the number of tones in the language spoken by culture Y. However, the number of tones in Y is also affected by cultural inheritance from an ancestor culture Z and borrowing from geographically proximal cultures N ₁ and N ₂. These associations need to be partialled out of the effect of X on Y. Furthermore, the local humidity is part of a complex network of climatic dependencies in time and space. This creates formal confounds for estimating the effect of X on Y, for example non-causal pathways via climatic dependencies and inheritance or borrowing. Therefore, the geospatial dependencies also need to be controlled for in order to accurately estimate the effect of humidity on the number of tones.

Figure 6. (a) Causal relationships between cultural variables (white circles) at different points in time (T0, T1, T2) and the environment (blue circles) at different geographic locations (L0, L1, L2) via inheritance (black and grey lines, e.g. from Z to Y), borrowing (red lines, e.g. from N ₁ to Y), and a target effect of the environment on the cultural variable (blue line, e.g. X to Y). (b) A simplification of the graph focused on node Y, with the environmental effects collapsed into one node.

In order to take account of these multiple confounds in a statistical model, we collapse the causal graph to focus on local effects (Figure 6b). In this causal model, T is the change in number of tones between an ancestor node and a descendant node (see Section 4.3.1). We model the change in tones from the ancestor node to the descendant node to be able to model the relative change with a Gaussian outcome rather than the absolute number of tones. This also removes the ancestor node in the predictor (mathematically, the two formulations Tones _descendant − Tones _ancestor ≈ predictors and Tones _descendant ≈ Tones _ancestor + predictors are equivalent).

The effect of the number of tones in phylogenetically close languages on tone change (P → T) represents confounds associated with genealogical relationships beyond the effects that are already controlled for (see Section 4.1). These residual effects include but are not limited to:

• • Rate of change – some lineages might exhibit a higher propensity to lose/acquire tones, i.e. a higher change rate.

• • Drift – some closely related languages might develop more tones independently (after the split from a common ancestor), but owing to a common cause inherited from their ancestor. In effect, if tones are not inherited from their ancestor, but a certain phonological situation is, that leads to the development of tones in the future.

It is calculated in the model as the average tones in a language's closest relatives or, in the second model, as a Gaussian Process (see discussion in Section 4.3.4).

Further, the effect of geographically close languages on the number of tone change (G → T ) represents geographical confounds associated with the (non-)acquisition of tones from geographically proximate groups. These influences include but are not limited to:

• • Intra-family contact – borrowing of tones from related languages that are geograpically close.

• • Contact – borrowing of tones from unrelated languages that are geographically close.

• • Isolation – effects from a lack of contact with other groups.

These effects can cause a change in tones since, for example, languages may increase their number of tones through borrowing when surrounded by languages with many tones. If, conversely, a language is surrounded by languages without tones, no tones can subsequently be borrowed. G is calculated in the model as the average tones in a language's vicinity or, in the second model, as a Gaussian Process (see discussion in Section 4.3.4). In the current work, we calculate this based only on Bantu languages in our sample. While there are neighbouring languages from other families that may have tone, we leave this for future work, since we do not have reliable data on the geo-phylo histories of those languages.

The Local Humidity node represents the climate network of humidity across time and space. If humidity influences tones, then we would assume it to do so on every level. Higher humidity would increase the change in tones between ancestor and descendant nodes but it would also increase the number of tones in the vicinity of a language and the tones in its related languages. In the main model only the effect of H → T is modelled. In order to model the assumption that the effect of humidity on language takes place at some point between the ancestor node and the child node, the model measures H as the average humidity of the location at the ancestor node and the descendant node. That is, the time points of changes from one humidity level to another between the ancestor and the descendant are unknown, but we assume that they average out across the sample. We assume consistency (Rehkopf, Glymour, & Osypuk, Reference Rehkopf, Glymour and Osypuk2016) and temporal stability: that the effect of humidity is not a composite measure and has the same effect in any situation. This might be violated by, for example, modification of ambient humidity by modern air conditioning, though this would only have very recent effects on the data.

The statistical model predicts T. We make a causal Markov assumption, meaning that three effects are not represented in the statistical model: H → P, H → G and P → G. This is because these effects are not directly causally relevant for predicting T.

We also assume that the causal network above is complete and that no major unmeasured confounds exist that causally influence both change in tones and humidity. This may not be realistic, but is a simplifying assumption for the sake of this study. Potential candidates for additional confounding effects include, for example, the effect of population size on rates of cultural change (Greenhill, Hua, Welsh, Schneemann, & Bromham, Reference Greenhill, Hua, Welsh, Schneemann and Bromham2018), or the relationship between climate and population dispersal (Honkola et al., Reference Honkola, Vesakoski, Korhonen, Lehtinen, Syrjänen and Wahlberg2013).

Finally, we assume that specific humidity is the only climatic variable that affects tone. While some studies of linguistic adaptation use temperature as their environmental variable, Everett (Reference Everett2017, p.4) explains that ‘temperature is an (imperfect) proxy for specific humidity because air at colder temperatures can “hold” less water’.

4.3. The inference models

For the computational modelling of the causal model outlined in Section 4.2, we employ the probabilistic programming language STAN (Stan Development Team, 2022) to design two different Bayesian models (named LIN and GP) which differ in their handling of the geography and phylogeny variables (see Section 4.3.4). The following sections go into detail about the individual components of the models and what parameters are used, while at the end, the full model equations are given in Section 4.3.5.

4.3.1. The model core

The cores of both models are the outcome and the main linear model:

$$\eqalign{T_{Desc, i}-T_{Anc, i}\;&\sim Normal( {\mu_i, \;\sigma_i} ) \cr {\mu _i} &= ( {{\alpha_\mu + \eta H_i + \gamma Gt_i + \rho Pt_i} ) {\rm d} A_i}}$$

Here, the change from the number of tones at the ancestor node T _Anc to that of the descendant node T _Desc is modelled as a normal outcome with mean μ and variance σ. The mean for each observation i is calculated as the sum of the population intercept $\alpha _\mu$, and the products of the linear coefficients η, γ, ρ, the respective humidity H and tones in geographically (Gt) and phylogenetically (Pt) close languages at observation i. This sum is multiplied by the temporal distance dA between the ancestor node and the descendant node, measured in units of 1000 years. This last term serves as a normalising factor which normalises the effects to effects per 1000 years. Doing this is important as we expect that the effects scale with time, i.e. languages that only recently split from their ancestor are expected to have less time to undergo any changes. This prevents the model from finding a null effect for languages with a small time distance from their respective ancestor.

4.3.2. Modelling time distance variance increases

The mean μ, however, is not the only part of the model affected by time distance between ancestor and descendant node. We also need to account for the possibility that with increasing time difference between the node pair, variance in the observation also increases.

On theoretical grounds, we can make such an assumption since we expect the unmeasured effects to compile over time such that with increasing temporal distance, the outcome is affected by increasing variance. In practical terms, the more time difference there is between the observation at the descendant node and the ancestor node, the longer other effects outside of this model can act on the process and increase the variance. The equation for this part of the model is as follows:

$$\sigma _i = \alpha _\sigma + \beta {\rm d}A_i$$

Here, for every datapoint i, the model calculates the individual variance to be equal to a global baseline variance $\alpha _\sigma$ plus the product of the coefficient β with the temporal distance dA between the node pair. Note that β is inferred as an unconstrained variable with a standard normal prior which results in an unconstrained slope. This relaxes the assumption about whether a variance increase exists in the data; the specific shape of a potential difference in variance dependant on dA is therefore inferred from the data.

Doing this prompts the model to assume the possibility that more distant descendant–ancestor node pairs fit the data less well. This again is intended to prevent the model from finding a null effect solely because of the unequal distribution of variance along the temporal axis. Since β is unconstrained, however, if there is not enough evidence for this in the data, the model will not assume such a variance increase.

4.3.3. Modelling uncertainty in reconstruction

As discussed in Section 2, the information about the tones at each node is drawn from a Bayesian phylogenetic model and the interior nodes of the tree therefore are reconstructions. However, since the tone values for interior nodes are posterior draws, we need to account for reconstruction uncertainty for these nodes. Therefore, we use this mechanism in the model to incorporate this reconstruction uncertainty:

$$T_i\;\sim Normal( {T_{REC, i}, \;T_{SE, i}} ) $$

Each interior node's tone value T is drawn from a normal distribution with a mean equal to the posterior mean of that node in the phylogenetic inference model T _REC and a variance equal to the standard error of this posterior distribution T _SE. Doing this ensures that the model recognises different degrees of uncertainty in the tone values of interior nodes.

4.3.4. Modelling contact and phylogenetic effects

The treatment of the potentially confounding variables of phylogeny and geography is where the two models (model LIN and model GP) differ. There are two ways of incorporating both variables in such models, each with their own drawbacks.

In model LIN, the geographical control is, simplified, the average tone count in the geographically/phylogenetically close languages. Model LIN is therefore a model with fewer variables that is additionally sensitive to which other languages are contemporary to the node in question. The drawback here is that it contains stronger assumptions about what constitutes geographic and phylogenetic closeness. Model GP foregoes these assumptions and models geography and phylogeny as Gaussian processes based on pairwise distance matrices. This makes the model better at modelling gradation in these variables, but at the same time it strongly increases the model complexity and, owing to the constraints of the kernel functions, cannot account for which languages are contemporaneous.

The specific implementation of these two different architectures is as follows. For Model LIN, the geography variable represents the per-node mean and variance of the number of tones in the languages that are contemporaneous to the node and less than 500 km away. Since the tone data are drawn from a phylogenetic tree where the splits occur at different times, the tones in the contemporaneous lineages are linearly interpolated according to the following function:

$$T_{Interpolated} = T_{Ancestor} + p( {T_{Descendant}-T_{Ancestor}} ) $$

where p is the proportional time point where the node in question is temporally located. For example, let N be a node for which the contemporaneous tone value of a neighbouring lineage is calculated. Assume that N occurs temporally exactly half way between the ancestor (five tones) and descendant (two tones) nodes of this neighbouring lineage. Figure 7 visualises this example.

Figure 7. Visualization of the interpolation of tones in neighbouring lineages at a given split time.

Here, the interpolated tone value of the neighbouring lineage at the time of node N is 3.5 since 5 + 0.5(2 − 5) = 3.5.

The phylogeny variable is devised in much the same way where for every node, the mean and standard deviation were calculated of all phylogenetically close nodes with a cophenetic distance of less than 0.02, which corresponds to a small phylum (see Figure A1 in the Appendix for an illustrative example). Note that, here, the nodes were not filtered for contemporaneity in accordance with the intended purpose of this variable as representing residual phylogenetic influences of a lineage (see Section 4.2).

Owing to how the two variables are designed, there is variance in the tone means, which is accounted for in the model:

$$\matrix{ {Gt_i\;\sim Normal( {Gt_{MEAN, i}, \;Gt_{SD, i}} ) } \cr {Pt_i\;\sim Normal( {Pt_{MEAN, i}, \;Pt_{SD, i}} ) } \cr } $$

Here, we estimate the average tones in the geographically/phylogenetically close languages (Gt, Pt) of node i as a variable drawn from a normal distribution with mean equal to the tone mean of the observed data (Gt _MEAN, Pt _MEAN) and variance equal to the standard deviation of the observed data (Gt _SD, Pt _SD). Doing this accounts for the inherent variation in the data aiming to incorporate the variance in the distribution of tones in close languages; taking only the mean could obscure important patterns.

In Model GP, the geography and phylogeny variables represent residual expected geographical/phylogenetic similarity between individual languages. The model architecture is built to infer these similarities from Gaussian processes on the basis of geographical and phylogenetic distance (as detailed above). The model structure is as follows:

$$\matrix{ {\left[{\matrix{ {\gamma_1} \cr {\gamma_2} \cr \ldots \cr {\gamma_n} \cr } } \right]} & {\sim MVNormal\left({\left[{\matrix{ 0 \cr 0 \cr \ldots \cr 0 \cr } } \right], \;G} \right)} \cr {\left[{\matrix{ {\rho_1} \cr {\rho_2} \cr \ldots \cr {\rho_n} \cr } } \right]} & {\sim MVNormal\left({\left[{\matrix{ 0 \cr 0 \cr \ldots \cr 0 \cr } } \right], \;P} \right)} \cr }$$

Here, the expected confounds for every node ($\gamma _1 \ldots \gamma _n,\ \rho _1 \ldots \rho _n$) are drawn from multinormal distributions with means of 0 and covariance matrices G and P as the kernel function.

The covariance matrix for the geographic confound (G) is a variation of a squared quadratic kernel whereas P represents a variation of the Ornstein–Uhlenbeck kernel:

$$\matrix{ {G_{tj} = \nu _G{\rm e}^{-\lambda _G\delta _{G, tj}^2 }} \cr {P_{rk} = \nu _P{\rm e}^{-\lambda _P\delta _{P, rk}}} \cr } $$

Here, δ = |x − x′| where x − x′ is the phylogenetic/geographical distance between two languages in a distance matrix. This modelling approach is similar to what is described in McElreath (Reference McElreath2018, chapter 14.5). As mentioned before, the geographical variable in this implementation cannot account for which languages are contemporaneous. This is due to the constraints of the Gaussian process kernels which can only process positive definite matrices, which requires all languages to be simultaneously present.

4.3.5. Full model summaries

After the detailed discussion of the model components, this section shows the models in full with all components put together. The first model (Model LIN) is a model where the effects of phylogeny and geography are modelled linearly. Recall that the model structure indicated here refers to the model using analysing the variable tones. The models for all four variables are equivalent with the exception that the priors used for the modelling of the vowel ratio variables are more informative: for those parameters that in the following equations show a prior of Normal(0, 1), we used Normal(0, 0.5) instead in these cases and instead of Exponential(1), we used Exponential(5).

(Equation of model LIN)$$\eqalign{T_{Desc, i}-T_{Anc, i}&\sim Normal( \mu _i, \;\sigma _i) \cr \mu _i &= ( \alpha _\mu + \eta H_i + \gamma Gt_i + \rho Pt_i) dA_i \cr \sigma _i &= \alpha _\sigma + \beta dA_i \cr T_i&\sim Normal( T_{REC, i}, \;T_{SE, i}) \cr Gt_i&\sim Normal( Gt_{MEAN, i}, \;Gt_{SD, i}) \cr Pt_i&\sim Normal( Pt_{MEAN, i}, \;Pt_{SD, i}) \cr \alpha _\sigma &\sim Exponential( 1) \cr \beta , \;\alpha _\mu , \;\eta , \;\gamma , \;\rho &\sim Normal( 0, \;1) } $$

The second model (Model GP) is a variation of the first in which the effects of phylogeny and geography are modelled using Gaussian processes:

(Equation of model GP)$$\eqalign{ T_{Desc, i}-T_{Anc, i}&\sim Normal( \mu _i, \;\sigma _i) \cr \mu _i &= ( \alpha _\mu + \eta H_i + \gamma _i + \rho _i) dA_i \cr \sigma _i &= \alpha _\sigma + \beta dA_i \cr T_i &\sim Normal( T_{REC, i}, \;T_{SE, i}) \cr \alpha _\sigma & \sim Exponential( 1) \cr \beta , \;\alpha _\mu , \;\eta &\sim Normal( 0, \;1) \cr \left[{\matrix{ {\gamma_1} \cr {\gamma_2} \cr \ldots \cr {\gamma_n} \cr } } \right]&\sim MVNormal\left({\left[{\matrix{ 0 \cr 0 \cr \ldots \cr 0 \cr } } \right], \;G} \right)}$$

$$\eqalign{\left[{\matrix{ {\rho_1} \cr {\rho_2} \cr \ldots \cr {\rho_n} \cr } } \right]&\sim MVNormal\left({\left[{\matrix{ 0 \cr 0 \cr \cdots \cr 0 \cr } } \right], \;P} \right)\cr G_{tj} &= \nu _G{e}^{-\lambda _G\delta _{G, tj}^2 } \cr P_{rk} &= \nu _P{e}^{-\lambda _P\delta _{P, rk}} \cr \lambda _G, \;\lambda _P&\sim Exponential( 0.3) \cr \nu _G, \;\nu _P&\sim Exponential( 0.5) } $$

In all models, the input variables were z-scaled (not centred) such that one unit in the model corresponds to one standard deviation in the data.

4.4. Model comparison

Since both ways of including geographical/phylogenetic confounds have advantages and disadvantages, we compared the two model approaches along with the individual confounds. For this, we used the R-package loo (Vehtari et al., Reference Vehtari, Gabry, Magnusson, Yao, Bürkner, Paananen and Gelman2022) to calculate the pairwise expected log predictive density (ELPD) for each model to obtain the most predictive model given the individual parameters. Table 1 shows the variables and model type that achieved the best scores.

Table 1. Variables in best models

Here, the checks indicate a control variable being included and GP indicates whether the model uses the linear or the Gaussian Process approach. We can see that in the models for the two vowel ratio models, phylogeny is more predictive than in the raw tones or raw vowels models. Further, the only Gaussian Process model is that of vowel ratio. The detailed coefficients of the model comparison analysis can be found in the appendix. Going forward, we only report the results from these models in the following sections, meaning that whenever we refer to the tones model, we refer to the model with the specifications indicated in Table 1.