Reliable detection and quantification of selective forces in language change

2.1 Maximum likelihood estimation methods

Maximum likelihood estimation is a conceptually simple yet powerful technique for estimating parameter values in a model and selecting between multiple candidate models. The basic setup is that we have both some empirical data, denoted X, and a probabilistic model that tells us how likely the observation X is given some choice of parameters Θ. We then estimate the values of the parameters by determining the combination that maximises the likelihood of the data. This procedure lies at the heart of many statistical methods, including linear regression. In such a model, parameters are chosen to maximise the likelihood of the data given a statistical model of the residuals (Severini 2000; Silvey 1970), for example, that the residuals are drawn from a normal distribution. It can also be viewed as a special case of Bayesian inference with a uniform prior.

In this work we are concerned with frequency time series, that is, a sequence of measurements X = {(x _t , t)} = {(x₁, t₁), (x₂, t₂), …, (x _m , t _m )} where x _i is the fraction of instances of use of a linguistic variable (e.g., all past-tense forms of a specific verb) in which a particular variant (e.g., the regular form) was used during a short time window centred on time t _i . Thus, the dataset X = {(0.2, 1), (0.5, 2), (0.75, 5)} would imply a proportion of usage of the regular form of 20 % at time 1, 50 % at time 2 and 75 % at time 5 (and no frequency data at any other time).

The underlying evolutionary model of language change determines a set of transition probabilities, Prob(x_i+1, t_i+1|x _i , t _i , Θ), that tell us how likely it is that, given a proportion x _i at time t _i and parameters Θ, the proportion will be x_i+1 at time t_i+1. In the previous example, the dataset X would determine the transition probabilities Prob(0.5, 2|0.2, 1, Θ) and Prob(0.75, 5|0.5, 2, Θ), whose exact numerical values would depend on the choice of model parameters Θ. We assume that contributions to changes in variant frequencies at different points in time are uncorrelated, which means that we can write the likelihood of the entire frequency time series as the product of the transition probabilities for each interval:

It is this likelihood function that we will maximise to determine the set of parameters Θ that best describes the cultural evolutionary dynamics, and that we will use to compare different models.

There are two main ways to choose the form of the transition probabilities Prob(⋅|⋅, Θ), a choice that is crucial to parameter estimation and subsequent interpretation. One could simply assume that the frequencies follow a prescribed trajectory between the two points, subject to some fluctuations around them. For example, in linear regression, frequencies would be assumed to vary linearly with time with normally-distributed residuals. Logistic regression is similar, but instead assumes that frequencies vary following a nonlinear logistic function that is commonly used as a model for S-shaped language change (Kauhanen and Walkden 2018; Kroch 1989; Reali and Griffiths 2010). A weakness of this approach is that without an underlying model of language production and transmission that may be operationalised as frequency time series, it is difficult to relate the parameters obtained to the behaviour of individuals or speech communities.

The alternative is to derive the transition probabilities from an explicit agent-based model of language change, many of which can be understood as a variant of the Wright–Fisher model of evolution (Blythe and Croft 2021). As noted in the introduction, the transfer of a model from genetics is justified on theoretical grounds (Croft 2000; Hull 2010) and one can interpret the parameters by appealing to models of language use (Baxter et al. 2006) or iterated Bayesian learning (Reali and Griffiths 2010). A drawback of the Wright–Fisher model is that exact expressions for the transition probabilities (Crow and Kimura 1970) are complex and difficult to work with computationally, as their associated transition matrices may become numerically intractable (Paris et al. 2019). This has motivated many different approximation schemes (Tataru et al. 2016). In this work, we apply a selfcontained Beta-with-Spikes approximation scheme that was developed and tested by Montero and Blythe (2023) and found to provide reliable estimates for parameter values without incurring undue computational cost. In the following we overview the conceptual components of this approach that are most relevant to linguistic applications, directing the reader to Montero and Blythe (2023) for technical details.

2.2 The Wright–Fisher model

The Wright–Fisher model describes a population of N replicating individuals of different types, each of which directly corresponds to a different variant of a linguistic variable. While any number of mutually exclusive types can be included in the model, we will focus on the case with two distinct types. Its extension to three or more distinct types is straightforward (Tataru et al. 2016). The quantity x _t is the proportion of individuals of a specific variant in generation t, as described above. The process of replication has the effect that, in generation t + 1, each individual is of the variant of interest with probability g(x _t , s), which depends both on the composition of the population and a measure of selection strength that we denote s. It is assumed that each of the N individuals in the new generation has its type assigned independently. This replication process is illustrated in Figure 1.

Figure 1:

Schematic representation of the transition from generation t to generation t + 1 in a Wright–Fisher process with N = 10.

The two evolutionary parameters, N and s, can be afforded a linguistic interpretation as follows. The effective population size, N, quantifies the scale of fluctuations in the variant frequencies around a smooth trajectory of change. The smaller the effective population size, N, the larger the fluctuations. When selection is weak (s is close to zero), the time taken for a variant to go extinct is proportional to N (Crow and Kimura 1970). Its interpretation as a population size, effortless in population genetics, does not work as well when studying language change. Through agent-based models of language learning and use (Baxter et al. 2006; Blythe and Croft 2021), we understand that N generically correlates with the size of the speech community. However, heterogeneous social network structures can result in N correlating only weakly with the size of the human population (Blythe and Croft 2021; Bromham et al. 2015; Wichmann et al. 2008). The population size N is also related to the total usage of all variants of a linguistic variable in a given generation, and how long it is retained in memory (Baxter et al. 2006; Reali and Griffiths 2010). Intuitively, speakers will be more consistent in their usage of specific variants the more they encounter them, and the longer they recall these encounters. This increased individual consistency will be reflected as smaller fluctuations in time series data. Although our analysis will not allow these different contributions to N to be distinguished, we will be able to determine which linguemes are subject to greater or lesser uncertainty in transmission between speakers.

The selection strength parameter, s, represents a tendency for the variant of interest to increase in relative usage (s > 0) or decrease (s < 0). Here, s subsumes all factors that could lead to a variant systematically increasing or decreasing in frequency over time, whether they originate in cultural, cognitive or language internal factors (Croft 2000; Labov 1994, 2001, 2010). Similarly to the various factors that may influence the effective population size, we will not be able to distinguish them from the value of s alone. However, as we show below, we can gain useful information by looking for common features of variants which are found to have similar selection strengths.

The parameter s specifies the probability g(x, s) that an individual in generation t + 1 is an offspring of an individual with frequency x and selection strength s in the previous generation. Here, we take

which has been commonly used in the theoretical characterisation of language change (Kauhanen and Walkden 2018; Yang 2000). In Figure 2 we plot the transition probability Prob(x|x₀, s) that results from this definition for the case where x0=12 and for different values of s. We see that larger values of s shift the peak of this distribution towards higher values of this frequency x, consistent with the notion of a bias towards the corresponding variant.

Figure 2:

Probability distribution of a variant frequency x after one generation of evolution in the Wright–Fisher model, starting from a frequency x0=12. As the selection strength s increases, the distribution becomes more sharply peaked on larger values of x.

In the literature, one can find relationships between the selection strength s and the probability g(x, s) different to that specified above (Crow and Kimura 1970; Tataru et al. 2016). Our chosen formula has the useful property that g(x, s) + g(1 − x, −s) = 1, which means that if one of two variants in a population has a selection strength s, the other one implicitly has a selection strength −s. This choice thus lends a symmetry between positive and negative selection strengths of the same magnitude, which aids the interpretation of the results. The strength s = 0 represents pure drift, where any changes in usage over time are due to the stochasticity of replication alone, and not the presence of selective forces. Under pure drift, g(x, s) reduces to g(x, 0) = x.

2.3 Beta-with-Spikes approximation

For a single generation of evolution in the Wright–Fisher model, the transition probability is the binomial distribution

because there are N individuals and a success probability of g(x _t , s). If the effective population size N is known, and the time between two frequency measurements corresponds to a single generation, one can use this expression for the transition probability in Equation (1) to construct the overall likelihood of a series of measurements. In the present application to linguistic corpus data, neither of these requirements hold. N is a parameter that we need to estimate, and measurement times are not in general separated by a fixed interval that constitutes a single Wright–Fisher generation. The Beta-with-Spikes approximation, introduced by Tataru et al. (2015) and extended by Montero and Blythe (2023), is designed to deal with these complexities.

For two observations made at times x _t and x_t+k (i.e., separated by k generations) the BwS approximation is

Here, P_0,k, P_1,k, α _k and β _k are parameters that determine the shape of the distribution. These parameters have the following interpretation. P_0,k is the probability that the variant has gone extinct by the kth generation, and P_1,k is the probability that it has driven the other variant to extinction in that time. α _k and β _k control the shape of the distribution of variant frequencies, conditioned on neither of them having gone extinct. These parameters can be determined from the mean and variance of this conditional distribution (see Montero and Blythe 2023). Note that all four parameters depend on N and s, as well as the sequence of observed frequencies xti. Therefore, they need to be recalculated for each time series and combination of model parameters.

A crucial advantage of the BwS approximation is that it accounts for the fact that changes in variant frequencies cannot be arbitrarily large. If a variant has a low frequency (x close to zero), then a downward fluctuation should cause it to become extinct, rather than attain a negative frequency. It is the spikes (represented by the delta functions) in the Beta-with-Spikes expression (4) that incorporate this constraint. By contrast, a normal approximation to the same transition rates (as used by Feder et al. 2014; Newberry et al. 2017) allows, in principle, arbitrarily large or negative x, instead of being constrained to the range 0 ≤ x ≤ 1. This difference is illustrated in Figure 3 which shows the statistical distances between the BwS and normal approximation and the exact WF transition probability, for different values of the initial frequency x₀ and two values of the selection strength s. We see that for both pure drift (s = 0) and strong selection (s = 0.5), the BwS approximation stays consistently closer to the exact distribution for all values of x₀, which is reflected in lower values of the statistical distance. In particular, the BwS approximation is significantly better than the normal approximation for initial frequencies x₀ close to the edges of the interval, and for strong selection.

Figure 3:

Comparison of the statistical distances of the BwS and normal approximations to the exact WF distribution as a function of the initial frequency x₀. Left: statistical distance for pure drift (s = 0). Right: statistical distance for strong selection (s = 0.5). The Beta-with-Spikes approximation has lower statistical distance to the exact distribution (meaning it approximates it more accurately) for every value of s and x₀, but especially for extreme values of x₀ close to 0.0 or 1.0 and for strong selection.

The main task in applying the BwS approximation is to estimate the parameters P_0,k, P_1,k, α _k and β _k for successive generations k = 1, 2, 3, …. The strategy of Tataru et al. (2015) is to match up the moments of the BwS distribution to those of the Wright–Fisher model after k generations have elapsed. This method works well when the selection strength s small, but less so when it is large. Montero and Blythe (2023) have improved on the method, particularly in the large s regime, by iterating (3) one generation at a time, and reading off the extinction probabilities, mean and variance at each. The code that implements this procedure, and generates parameter estimates is available here.

In the context of cultural evolution, it is not obvious what period of time counts as a generation in the Wright–Fisher model. In principle, this is a free parameter which would also need to be optimised by maximum likelihood estimation (and furthermore demand an interpretation). Fortunately, this is unnecessary. Montero and Blythe (2023) further show that the optimum values of 1/N and s are both proportional to the chosen generation time. In other words, the generation time serves only to set the units in which the parameters N and s are measured. It is however important to use the same generation time across multiple time series when one wishes to compare the values of N and s that are obtained: otherwise, they would be in different units and not comparable. In this work we generally take the shortest time between successive observation points as the generation time. If one makes it shorter than this, the computational effort increases without any improvement in the quality of the estimates obtained. If one makes it longer, one must then aggregate multiple data points which then entails a loss of temporal resolution. However, as we discuss below, it is sometimes beneficial to combine data points to reduce sampling error that is not accounted for in the present maximum likelihood analysis.

2.4 Distinguishing selection from drift

As established in the introduction, the social, linguistic and cognitive forces driving language change are very diverse. Still, their measurable effects can be broadly characterised as belonging to one of two types. Systematic biases drive the evolutionary process in a specific direction, and can be modelled as selective forces. Frequency effects and stochasticity in transmission produce random, unbiased drift whose effects are always present, albeit not always sufficient to explain the behaviour of the data. Quantitative, empirical analyses benefit from the simple yet powerful and flexible characterisation of language change afforded by this binary description.

By using the transition probabilities (4) in the likelihood function (1), we can find the maximum likelihood values of the effective population size N* and selection strength s* via

In practice, we find that the likelihood function L(X|N, s) has a single maximum, which can be located by successively optimising on N at fixed s and vice versa.

It is important to establish whether the selection strength s* is significantly different to zero: otherwise, the null model of stochastic drift (s = 0) would be sufficient to explain the behaviour of the data without the need for selection (Blythe 2012; Newberry et al. 2017). In order to do this, we compare the maximal likelihood under selection, L(X|N*, s*), with the maximal likelihood under pure drift. That is, we first restrict to s = 0 and determine the optimal effective population size N0*:

Then we compare the models with and without selection by computing the likelihood-ratio

This quantity can be compared to a reference distribution to find a p-value, an estimation of the probability that the observed time series could have arisen from drift alone^[1] (Severini 2000; Silvey 1970). To achieve this, we follow the procedure outlined by Feder et al. (2014) and generate 1,000 artificial time series spanning the same time period as the empirical data X with parameter values s = 0 and N=N0*. For each of these we compute the maximum likelihood values N*, s* and N0*, using the same sequence of steps as for the original empirical time series. We then compute the likelihood ratio λ and determine what fraction of the artificial time series has a larger λ than the one that was observed. This provides an empirical p-value for the null hypothesis of drift.

3.1 Drift versus selection in past-tense English verbs

A simple example of competition between two variants is provided by English verbs with an irregular past tense form which in many cases coexists with a regular form. This competition has been studied from a variety of quantitative perspectives (Cuskley et al. 2014; Karjus et al. 2020; Karsdorp et al. 2020; Lieberman et al. 2007; Newberry et al. 2017; Ringe and Yang 2022). Of greatest relevance to the present work are those studies that aimed to distinguish drift from selection as the mechanism behind changes in the relative frequencies of the regular and irregular forms over time.

Newberry et al. (2017) applied the Frequency Increment Test (FIT, Feder et al. 2014) to a set of verbs from the Corpus of Historical American English (COHA, Davies 2010). This is a maximum-likelihood method that rests on a normal approximation to the Wright–Fisher transition probabilities. Like the Beta-with-Spikes maximum-likelihood method in Section 2.4, this method yields estimates of the effective population size and selection strength, along with a p-value for the null hypothesis of pure drift. However, there are situations where results are flagged as unreliable due to the frequency increments failing a normality test (Newberry et al. 2017). Karjus et al. (2020) further noted that the results can also be sensitive to the size of the window over which frequencies are estimated.

Karsdorp et al. (2020) avoid these issues by taking the rather different approach of training a neural network on simulated time series generated by the Wright–Fisher transition probabilities (3) for different values of s (but fixed N = 1, 000). Each time series in the training set is labelled according to whether it was generated purely by drift (s = 0) or if selection was operating (s ≠ 0). Once trained, the network yields a binary classification of empirical time series, according to whether they are more similar to the examples of drift or selection in the training data. We refer to this as the Time Series Classification (TSC) approach. The advantage of TSC is that no approximation to the Wright–Fisher transition probabilities is made. Moreover, one can manipulate the training data so that it displays artifacts of binning or finite sample sizes that are features of real time series, which in turn should improve the reliability of the classification. This approach does however come with some drawbacks. Whilst the output from the classification algorithm is a value between 0 and 1, it does not have an obvious interpretation as a probability. Karsdorp et al. (2020) used a threshold of 0.5 to label timeseries as arising from drift or selection. The method further does not provide an estimate of the strength of s, and since N was fixed in the training set, this amounts to an assumption that this single value of N was appropriate for all empirical time series. This could be an issue since Newberry et al. (2017) report a wide range of values of N for this data set (from around 80 to around 22,500).

In Section A of Appendix we report the maximum likelihood estimates of N and s, along with the p-value for the drift hypothesis, obtained using the BwS method for the same set of verbs that were considered by Newberry et al. (2017) and Karjus et al. (2020) using FIT and by Karsdorp et al. (2020) using TSC. We perform the analysis by extracting annual frequency data of the variants of interest from COHA and aggregating it into 10-, 20- and 40-year bins. The reason for this is a trade-off between the more precise frequency estimates that derive from larger bins and the greater temporal resolution obtained from a larger number of bins over the relevant historical period. By employing different binning strategies, we can gain insights into the consequences of this trade-off. Variable-width binning strategies have also been successfully applied in previous studies (Newberry et al. 2017). In these, the number of tokens per bin is kept roughly constant at an arbitrarily chosen value, at the expense of varying their temporal width. For the purpose of comparing the different methods, we have chosen to look only at fixed-binning strategies, although the BwS method could be combined with variable-width binning.

We focus first on the role played by selective forces, which we quantify by appealing to the p values associated with the null hypothesis of pure drift as described in Section 2.4. In Figure 4 we compare the results obtained from the three different methods by ordering the verbs from left to right by decreasing BwS p-value, averaged over the three temporal binnings. Each panel corresponds to a different analysis method, and indicates the p-value for the hypothesis of pure drift for each verb and binning protocol. We recall that higher p-values are more suggestive of the historical changes being due to drift: these are represented with colours ranging from light to dark blue, with darker colours representing higher p-values. Meanwhile, low p-values point towards other forces (such as selection) being present and are represented with different shades of red. While we use the standard p-value threshold of 0.05 in the transition between blue (drift) and red (selection) in this representation, we acknowledge that these mechanisms lie in a continuum by making the transition between these extremes smooth.

Figure 4:

Results for the detection of selective forces in 36 COHA verbs, with three different methods and for three different temporal binnings of 10, 20 and 40 years. Results for both the FIT and BwS likelihood-ratio algorithms produce a p-value for the pure drift hypothesis. Blue shades represent higher p-values (i.e., higher likelihood of the data under drift), while red shades represent p-values under the traditional 0.05 threshold of significance for selection. Time series where the normal approximation that FIT relies on is inaccurate are crossed out. Results for the TSC method from Karsdorp et al. (2020) are classified in a binary way as either drift or selection. The average p-value across the three bins widths obtained through the BwS algorithm is shown along the horizontal axis. We note that the BwS method gives results consistent with TSC when FIT is unreliable.

We see from Figure 4 that the three distinct methods give broadly consistent results, with those verbs towards the left being more compatible with change through pure drift, and those to the right with change from selection. More precisely, the correlation coefficients between the p-values obtained with different methods are 0.63 (Pearson) between FIT and BwS, 0.68 (biserial) between TSC and FIT, and 0.62 (biserial) between BwS and TSC. Analyses producing high p-values for selection (i.e. implying that drift alone can explain the behaviour of the data) are indicated with blue colours, whereas those where selection is more significant are red. Results obtained through the FIT method are generally consistent with those obtained with the BwS method. However, 30 of the FIT results (27.8 % of the total) are flagged as ‘unreliable’ due to a failed normality test. These reliability issues are designed out of the BwS method, as it does not require normally-distributed increments (Montero and Blythe 2023; Tataru et al. 2016). Confidence in the method’s reliability is also gained by benchmarking with synthetic and genetic data (Montero and Blythe 2023) and through the consistency with the independent TSC results. The higher precision of the BwS at high selection strengths leads to higher significance (lower p-values) in its detection of selective forces when compared to the normal approximation, leading to redder colours in Figure 3.

The TSC appears to give a cleaner classification of verbs according to drift and selection, and greater consistency with different choices of bin size. This is likely due in part to the training data being subjected to the same binning protocol as the empirical time-series, but also because a strict threshold was applied to the neural network’s output value to partition into the two classes. While the TSC neural networks produce a value between 0 and 1 as their output, making it more nuanced than this binary classification would suggest, this number is not a probability or a p-value like those produced by BwS or FIT. Thus, an arbitrary threshold is necessary in order to classify time series as driven by drift or selection. A higher or lower threshold would put the boundary between the two classes in a different place. This hinders the interpretability of the result and the estimation of significance levels.

Our results further demonstrate that variation in p-values under different binning strategies, previously observed within the FIT analysis (Karjus et al. 2020), remains evident under the less restrictive BwS analysis presented here. We consequently regard this variability as an inherent feature of the time series data: that is, some changes are harder to classify than others. That is, this uncertainty need not be a failure of the method, but a reflection of linguistic reality. For example, it could reflect different variants being used less predictably by speakers, or by the constraints on variation changing over short timescales (Tagliamonte 2011).

Such observations motivate a more detailed investigation of the classifiability of individual time series. A time series that shows limited variation in parameter values under different temporal binnings is more classifiable than one that shows more variation. With our interest in selection, the two most relevant parameters are s, the selection strength, and the p-value associated with the drift hypothesis. We can visualise the variation in these parameters by performing a Principal Components Analysis (Jolliffe 2002) on combinations obtained through different binning strategies (in this case, bins of 10, 20 and 40 years). The interior of the resulting ellipses indicates the range of variation of the two parameters over different binning strategies. This way, they provide a visualization of not just the average, but the uncertainty and covariance of s and the p-value under different binning strategies. We show these ellipses for the COHA verbs in Figure 5. The upper panel contains the full range of p and s values obtained through the analysis, while the lower panel zooms in on the region where the drift p-value is smaller than 0.05 (i.e., the conventional threshold for rejecting the null hypothesis). We see a correlation between the maximum likelihood value of s and the p-value (both through the positions and rotation angles of each ellipse).

Figure 5:

Variability in the selection strengths s and p-values for the null hypothesis pure drift for the COHA verbs. Each cross shows the mean value of the two parameters for each verb obtained when aggregating frequencies into temporal bins of different lengths. Each ellipse indicates the variability in the parameters at the level of one standard deviation. The vertical axis is an indicator of selection, defined as one minus the p-value associated with the drift hypothesis. The lower panel shows those verbs that fall within the range of p-values that is conventionally used to reject the null hypothesis for a single observation. In this panel we see a clear split into those that are regularising (negative s) and are irregularising (positive s).

The ellipses that lie entirely within the lower panel correspond to the verbs that are most likely to be driven by selection. We see a clear split between four verbs with positive selection (catch, light, wake and quit), which corresponds to them becoming more irregular over time, and six verbs (learn, lean, burn, smell, dwell and spill) with negative selection, and thus regularising over time. In this analysis, the frequency x is the fraction of irregular forms used in the relevant context. Across the entire plane, there is evidence of both regularisation and irregularisation, although in most cases it is difficult to rule out drift as an explanation for the changes, as was observed by Newberry et al. (2017).

In interpreting these results, it is important to recognise that the presence of fluctuations around a smooth change trajectory will tend towards a higher drift p-value, since in the analysis drift is the sole source of fluctuations. It is possible that fluctuations in the corpus derive from other sources, such as sampling effects associated with a finite corpus. Some methods for estimating parameters in the Wright–Fisher model attempt to account for such fluctuations separately to drift (Tataru et al. 2016). These are, however, typically difficult to implement, and instead we sidestep the issue by ensuring sufficiently many tokens in each temporal bin that the frequency is well estimated. As such, we might expect to see stronger evidence for selection as the bin width is increased, which appears to be true for some (but not all) of the verbs with intermediate p-values. This suggests that some language changes may be dominated by the random effects of drift and therefore exhibit strong fluctuations even in very large corpora.

To summarise, we have shown in this section that the BwS method can be readily applied to historical corpus data for changes in the frequencies of linguistic variants. It provides estimates of parameters in the Wright–Fisher model that do not rest on an assumption that frequency increments are drawn from a normal distribution, and we find broad consistency in the strength of support for a drift hypothesis with complementary methods.

3.2 Competing linguistic motivations in English verbs

In the previous section we observed a split between some verbs that were regularising and some that were irregularising. While the extension of regular inflection at the expense of irregulars seems to be the norm (e.g. Bybee 1995; Sims-Williams 2016), irregularisation is however an attested phenomenon. Cuskley et al. (2014) found that the processes of regularisation and irregularisation tend to take place with similar frequency, something that is also perhaps suggested by Figure 5, which shows a similar density of verbs along the branch with positive s (towards irregularity) and negative s (towards regularity). Ringe and Yang (2022) suggest that irregularisation may occur if the number of verbs within an irregularity class is high enough to surpass a productivity threshold. Following Bybee (2001) and Prasada and Pinker (1993), both Cuskley et al. (2014) and Newberry et al. (2017) propose phonological analogy as a potential mechanism for irregularisation. Couched in the terms of the present work, this would correspond to the general rule (adding -ed) contributing a negative value to s whilst rules that apply only to a specific subset of verbs contribute a positive value to s. Note that we do not necessarily imply that these contributions are additive: for example, in optimality theory (Hayes 1999; Prince and Smolensky 1997), higher-ranked rules take precedence over lower-ranked rules. In general, we may regard opposing forces on linguistic variation as arising from competing motivations which have been discussed in a variety of language change contexts (e.g., Bates and MacWhinney 1987, 1989; DuBois 1985; Haiman 1983; Hawkins 2004; Kirby 1997). By whatever mechanism this opposition is resolved, an overall positive s value here indicates that the irregularising rule is dominant.

In this section, we investigate a distinct motivation that may favour irregularisation, namely the phonological simplicity that is afforded by omitting a sound repetition that would occur under application of the regular rule. Specifically, we consider verbs whose infinitives end in alveolar stops (/d/ or /t/) and have an irregular past form where the regular -ed termination is omitted. Examples include I bled instead of I bleeded or she bet instead of she betted. Verbs where devoicing of final /d/ or changes in the root vowel take place on top of the omission of the termination are also considered. Thus, we hypothesise that the regular form is preferred from the point of view of inflectional simplicity (i.e. using the regular everywhere leads to a simpler inflectional system), while the irregular form is favoured by phonological simplicity. By applying the BwS algorithm to estimate the s parameter (and in particular, its sign), we can assess how these competing motivations play out.

For this investigation we switch to the 2019 English Google Books corpus (Michel et al. 2011), as the number of verbs falling into this category and whose past tense forms are both sufficiently frequent and can be reliably identified is relatively small. The larger size of Google Books relative to COHA allows more examples to be included. We identified 19 English verbs whose irregular and regular forms both show usage above 1 % at least in one 5-year bin in the Google Books corpus in the considered time frame (1809–2009). These verbs are: bend, bet, bite, blend, build, fit, glide, knit, light, pat, plead, quit, slide, speed, spit, thrust, tread, wed, and wet. A difficulty in the analysis is that the irregular past-tense form can coincide with certain present-tense forms. A major exception is when the verb is preceded by a third-person singular pronoun (e.g., the present he bets vs. the irregular past she bet), which can easily be distinguished in the bigram dataset. We recognise that this separation is not perfect: for example, certain English varieties do not use the third person marker -s, but we consider the effect of these contributions to be negligible in the corpus. We also kept only those cases where the pronoun was judged to appear at the start of a sentence (by virtue of capitalisation), so as to exclude contexts where the pronoun is followed by the infinitive in a question or an inversion. Again there are situations where capitalised pronouns can appear mid-sentence, but these are also rare. With this, total counts of usage for verbs with non distinct irregular past tense forms range roughly between 2,600 (knit) and 120,000 (pat), while counts for verbs whose irregular past tense is distinct from their base form range between 600,000 (glide) and 40,000,000 (build).

In order to formally test whether a potential bias towards irregularisation is significant, a similar analysis was carried out on a baseline set of 34 English verbs whose base form does not end in /d/ or /t/. Data was extracted from Google Books and all verbs satisfy the same conditions on minimal usage in the time frame of interest (1809–2009). The chosen verbs are: awake, blow, burn, catch, cleave, creep, dive, dream, dwell, freeze, grow, hang, heave, hew, kneel, lean, leap, learn, shake, shear, shine, slay, slink, smell, sneak, spell, spill, spoil, strew, string, strive, swell, wake and weave. Total usage for these verbs in the Google Books corpus for the specified period ranges between 211,000 tokens (slink) and 31,900,000 (learn), in the same orders of magnitude as the /d/,/t/ set.

The maximum likelihood parameters for these 53 verbs are given in Section B of Appendix. Here, we visualise our findings by plotting ellipses in the plane spanned by the selection strength and the indicator of selection, following the same procedure as previously described for the COHA verbs, albeit with the addition of a 5-year temporal binning strategy. With this, each ellipse in the s-p plane for each verb is produced by averaging the results of the analyses of at most four temporal binnings. We recall that these ellipses characterise the variability in these parameters as the temporal binning is varied. The upper panel in Figure 6 shows the results for all 53 verbs.

Figure 6:

Parameter estimates for verbs ending in alveolar stops (red) and verbs in the baseline set (blue) in the Google Books data set. The top panel shows the entire range of drift p-values and includes all 53 verbs. The bottom panel is restricted to p < 0.05, thus focusing on verbs that are likely to be undergoing directed selection. The distribution of verbs in the alveolar stop set seems to be skewed to the region where s > 0 and p < 0.05, suggesting they are more likely to be irregularising than the other verbs.

For the purpose of comparing the two sets of verbs, we partition the s-p plane into four regions: those with positive or negative selection strengths; and those where the p-value falls above or below 0.05. The lower panel of Figure 6 zooms in on this latter region, which we may regard as showing evidence of selection. In both panels, red crosses and ellipses correspond to verbs ending in alveolar stops, while blue crosses and ellipses correspond to verbs in the baseline set. Given our interest in irregularisation, three groups of verbs can be identified. Sixteen verbs (awake, bend, bet, bite, catch, fit, hang, light, quit, shake, slide, sneak, spit, strew, wake, wed) have their confidence regions (ellipses) completely contained in the region of likely selection of the irregular form (p < 0.05 and s > 0, lower-right panel). Of those, 9 are in the alveolar stop set and 7 are in the baseline set. Six verbs (freeze, kneel, leap, plead, swell, thrust) have confidence regions only partially contained in this region of the s-p plane, indicating that, while selective forces towards the irregular form are a plausible explanation to their dynamics, the pure drift hypothesis cannot be confidently ruled out. The remaining 31 verbs (8 in the alveolar stop set, 23 in the baseline set) have confidence regions contained entirely outside this region of likely irregularisation.

These results suggest that verbs in the alveolar stop set are more likely to be selected towards their phonologically simpler irregular form than their counterparts in the baseline set. To test the significance of these findings, we construct the 2 × 3 contingency table shown in Table 1, where one dimension expresses belonging to the alveolar stop or the baseline sets, while the other dimension expresses whether the verbs’ ellipse falls in the irregularisation region in the bottom panels of Figure 6. The p-value for the null hypothesis that the baseline and alveolar stop verbs are drawn from the same distribution is 0.031, as obtained by applying the G-test of goodness-of-fit to the contingency table (McDonald 2014). This indicates that the specific rule favouring phonological simplicity likely outcompetes a general tendency towards regularity.

It is possible that other effects may be responsible for this subset of verbs tending to irregularise. For example, it is well understood that higher frequency items tend to tolerate greater irregularity (Bybee 2007). Given the selection criteria imposed to arrive at the set of 12 verbs in this analysis, it is possible that the sample is skewed towards higher frequency and more irregular forms. However, as noted, the total token counts for both the baseline set and the alveolar stop set span similar ranges, and also have similar averages (of around 5 million for both sets). Therefore we consider this alternative explanation unlikely.

This is not the only phonological conditioning on irregularisation that can be inferred from Figure 6. The subset of verbs ending in a short vowel plus a lateral (dwell, smell, spell, spill, swell) seem to be a lot more likely to regularise under selection than other verbs in the study. A similar G-test to the one performed on the alveolar stop test on Table 1 reveals that this tendency is significant with p < 0.003. The origin of this tendency is, however, unclear.

To summarise, in this section we have shown that by focussing on a subset of verbs that are subject to specific combination of competing motivations, the Wright–Fisher model combined with the BwS approximation can be used to determine the net effect of this competition. Specifically, we have acquired evidence that phonological simplicity dominates inflectional simplicity in this competition, suggesting perhaps that this is an instance of an OCP constraint (Obligatory Contour Principle, Leben 1973; Stemberger 1981). OCP constraints disfavour pairs of identical or near-identical consonants from being in close proximity to each other. In particular, the constraint here appears to be an OCP-place constraint (Frisch et al. 2004; McCarthy 1986; Pozdniakov and Segerer 2007), meaning that it does not just affect identical consonants, but all alveolar stops independently of voicing.

3.3 Spanish spelling reforms

So far we have assumed that the evolutionary parameters (the effective population size N and the selection strength s) have been constant over time. In the case of competition between regular and irregular verbs this is a reasonable assumption, due to the factors favouring one over the other likely being cognitive or linguistic in origin. By contrast, social pressures like prestige, taboo, or language contact (Hernández-Campoy and Conde-Silvestre 2012; Labov 2001; McMahon 1994) are inherently time-dependent, and we may expect the selection strength in particular to change over time. Here, we investigate this possibility in the context of a purposeful change made by a regulating institution through prescriptive grammar and spelling rules (Anderwald 2012; Rubin et al. 2013), the acceptance or rejection of which we expect to be reflected by a change in the value of s. While well established algorithms like change-point analysis (Taylor 2000) exist for the detection of change in time series, these suffer from shortcomings that make them inadequate for a more nuanced analysis of change in language and culture. First, change-point analysis is based on the assumption that the data is distributed around a constant average before and after a change, which changes the value of said average instantaneously. This makes this methodology only fit for the detection of rapidly occurring S-shaped curves of language change, where the usage frequency of a variant quickly changes and stabilises. Secondly, change-point analysis provides no extra linguistic information, as it does not assume a model of the underlying evolutionary dynamics. Montero and Blythe (2023) solve this issue by setting out a procedure for estimating times at which the parameters N and s change, thus measuring changes in the evolutionary dynamics of the data rather than its average. We briefly recapitulate and then apply this method below.

The specific changes of interest are spelling reforms in Spanish that were introduced by the Real Academia Española (RAE), the central regulatory institution of the standard Spanish language. Since its creation in 1713, the RAE has regulated Spanish orthography following the phonemic principle over etymological or conservative approaches (Baddeley and Voeste 2012). We study words affected by one of the following reforms: (A) The simplification of the <ss>  digraph to a single <s>  in 1763, due to the different sounds that both spellings represented having merged in the 16th century (Real Academia Española 1763); (B) The replacement in 1815 of etymological <x>  with <j>  in all non word-final contexts where it represented the phoneme /x/ (Real Academia Española 1815); (C) The replacement, also in 1815, of <y>  with <i>  in all non word-final closing diphthongs; (D) The reversal of accentuation rules for words ending in <n>, introduced in 1881. This reform stipulated that words ending in <n>  with a tonic last syllable had to be accentuated, while words ending in <n>  with a tonic penultimate syllable lost their previously prescribed accent (Real Academia Española 1881). We treat words that gain an accent and words that lose an accent as independent sets (D.1 and D.2, respectively).

We now seek to estimate the time at which each reform occurred by appealing only to the time series data and no external information. The basic idea (see also Montero and Blythe 2023) is to allow different parameter combinations (N, s) to apply before and after a time T. That is, for t < T the Wright–Fisher model with parameters (N₁, s₁) applies, and for t > T the parameters (N₂, s₂) apply. The data likelihood, obtained by combining Eqs. (1) and (4), is then maximised with respect to all five parameters (i.e., N and s each before and after the change, and the time T of the change itself).

After identifying the time T that maximises the data likelihood, one needs to determine if the additional complexity of the five-parameter model is compensated by a sufficiently improved description of the data. To achieve this, we obtain an empirical p-value for the null hypothesis that the selection strength s was constant over the entire time period by following a procedure similar to that described in Section 2.4. Specifically, we determine the maximum likelihood values of N and s without a change point, and generate 500 synthetic time series that match the length of the observed series with these parameter values. For each of these time series, we then optimise the five-parameter likelihood that applies when the selection strength changes at a single point in time. An empirical p-value is then given by the fraction of such time series whose five-parameter likelihood exceeds that of the real trajectory. Although computational constraints limit the number of synthetic time series that can be analysed this way, we find that situations where the five-parameter fit has a high likelihood are extremely rare, and there is little to be gained by estimating their rarity to greater precision. One can then apply a threshold, e.g., p < 0.05, to decide whether to accept the more complex model. Having split the time series once, one can apply the method again to each sub-series, thereby identifying secondary change points. This procedure terminates when none of the sub-series admits a subdivision that yields a sufficiently improved description of the data according to the threshold that has been imposed.

To apply this method to the Spanish Spelling reforms, we identify a set of commonly used words that are affected by each one, and average the relative frequencies of usage of their old spellings over all members of each set. The number of words in each set ranges from 16 to 27. The exact sets are specified in Section C of Appendix. This procedure generates a single effective time series for each of the reforms, and has been found effective in related corpus analyses (Amato et al. 2018).

While this averaging over sets of words decreases the sampling noise in the data and increases the inferential power of the analysis, cultural data still suffers from issues that may affect the applicability of the method. Particularly, corpora tend to contain lower token counts in earlier time periods. When translated to frequency time series, this leads to greater sampling noise fluctuations that may be misidentified as changes in the effective population size parameter N. This issue can be remedied by applying a sampling error equalisation algorithm, as laid out by Montero and Blythe (2023). This method creates subsamples of the larger token counts in the data set, in such a way that sampling effects are of equal magnitude throughout the data. In this way, any significant changes in N detected by the method must be due to changes in effective population size parameter, and not a consequence of unequal sampling noise.

Our results are shown in Figure 7. Despite the aggregation of words within each category (to improve the inferential power) and 5-year bins (to reduce computational effort), we find that the resulting trajectories are still subject to considerable fluctuations. The frequency plotted is that of the deprecated variant, which we find is eliminated in all five cases—this highlights the acceptance and influence of the Real Academia Española amongst the literate population. We show with a red line and dot the time at which the reform was introduced, and with a black dot and solid line the first time T at which subdividing the time-series improves the fit to the data, with a p-value threshold of p = 0.05 applied.

Figure 7:

Application of the BwS algorithm for the detection of changing forces to the reference data set of Spanish spelling reforms in the 2019 Spanish 1-g Google Books corpus, with temporal binning of the frequency data of 5 years. For each set of words that undergo a rule change, the ratio of usage of the old form is plotted over time. The ratio of usage of all old forms converges to zero after each reform. Red dots with solid vertical lines represent the year of publication of the RAE spelling reforms (Real Academia Española 1763, 1815, 1881). Dark blue dots with solid vertical lines represent the year at which selection strengths changed as detected by the maximum likelihood method with a p-value below 0.05. These fall within a period ΔT of 12 years or less relative to the date of the reform. Note that the temporal resolution of the time series is of 5 years, so an error of 10 years is equivalent to just two data points. Dashed vertical lines represent secondary points of change in evolutionary parameters, also detected with a p-value below 0.05. The number of such secondary points depends on the time series.

In all five cases, we find evidence that the selection s changed significantly over time. In each case, the first detected change point falls within twelve years of the reform being introduced, even when the trajectory is strongly fluctuating. We note that the algorithm does have a tendency to detect the reform after it has occurred, rather than at its inception. This is due to the algorithm not distinguishing past from present, making both the beginning and the end of the sharp decline following a reform are considered equivalent.

Further iterating the algorithm, we can further subdivide the time series, as described above. In doing so, we detect secondary time divisions (dashed lines in Figure 7), whose p-values are below 0.05. In time series (B), the earlier secondary point detects the beginning of the rapid decline in usage that was deemed less significant than the end by the first application of the algorithm. The later secondary point in (B) and the secondary point in (A) are not associated to documented reforms, and may reflect slight changes in social attitudes or simply be quirks of the data.

Table 2 further records the s-values before and after the main change point. All s, N and p-values for every main and secondary point detected by the algorithm can be found in Section D of Appendix. For categories (A), (D.1) and (D.2), we find that the s value decreases after the detected year of the reform, corresponding to an acceptance of the reform by the speech community. The other two categories however show the opposite trend, with the s value becoming less negative across the reform. We note from Figure 7 that both categories (B) and (C) feature a rapid elimination of the deprecated form, and that this change was in progress before the reform was introduced. It has been suggested that in many cases, language reforms tend to reflect pre-existing trends, as opposed to actuating the change (Rutten and Vosters 2021). Our analysis provides further evidence of this, and further suggests that the impact of the reform on the speech community may be limited in such cases.

In summary, this analysis indicates that the BwS method can be used successfully to characterise evolutionary forces that change over time from time series data alone. As an unsupervised method, it does not rely on any prior knowledge as to when the change may have occurred, although it does benefit from a large sample size being available, obtained here by aggregating multiple instances of a change together. We have found that the estimated time at which the selection strength changed corresponds well with the time at which the corresponding reform was introduced, and comparing these strengths before and after the reform allows us to assess its impact on the speech community.