The Evolutionary Dynamics of Negative Existentials in Indo-European

Where in earlier work diachronic change is used to explain away exceptions to typologies, linguistic typologists have started to make use of explicit diachronic models as explanations for typological distributions. A topic that lends itself for this approach especially well is that of negation. In this article, we assess the explanatory value of a specific hypothesis, the Negative Existential Cycle (NEC), on the distribution of negative existential strategies (“types”) in 106 Indo-European languages. We use Bayesian phylogenetic comparative methods to infer posterior distributions of transition rates and parameters, thus applying rational methods to construct and evaluate a set of different models under which the attested typological distribution could have evolved. We find that the frequency of diachronic processes that affect negative existentials outside of the NEC cannot be ignored—the unidirectional NEC alone cannot explain the evolution of negative existential strategies in our sample. We show that non-unidirectional evolutionary models, especially those that allow for different and multiple transitions between strategies, provide better fit. In addition, the phylogenetic modeling is impacted by the expected skewed distribution of negative existential strategies in our sample, pointing out the need for densely sampled and family-based typological research.


INTRODUCTION
The negative existential domain, the expression of negated existential statements, may appear to be a simple, unremarkable, area of grammar. In many languages, it simply involves the deployment of the usual means used to express the affirmative existential with the standard verbal negation marker. This leads to clauses such as (1), from Swedish (Indo-European, Germanic), where the structural coding means involved in expressing existence are deployed alongside the Swedish standard verbal negation marker inte.
A thorough look at the negative existential domain in the languages of the world, however, suggests that it is expressed by a variety of construction types. In many languages, one finds a dedicated negative existential marker used as a negative existential copula. This marker is not used to negate verbal predicates, but may be used to negate other domains of nominal predication such as predicate location. This is illustrated in (2) from Turkish, where the negative existential marker yok is deployed, but not the Turkish standard verbal negation marker, a suffix. In another construction type, the standard negation marker is used as the only marker of negative existence, without the existential marker used in affirmative existential clauses. This is illustrated by the Tongan examples in (3a-b). The negation marker 'ikai is used in (3a) to negate the main verbal predicate, and in (3b) as the sole marker of negative existence, without another existential marker.
Turkish (Turkic; own knowledge) Tongan (Austronesian, Polynesian, Veselinova 2014: 1342 This cross-linguistic variation led Croft (1991) to propose a typology for the negative existential domain with synchronic and diachronic components (fully explained in the following section). In the synchronic portion, Croft identifies the construction types illustrated by (1-3) above as well as three intermediate construction types. The six types are defined based on a comparison of the negation marker(s) used to negate verbal predicates and the expression of negation in negative existential clauses: are they identical? Distinct? Is the verbal negation marker used as the negative existential copula, or do we find some intermediate situation? Croft's synchronic typology has been successful as a cross-linguistic taxonomy of negativeexistential constructions and fits well with other variables in the typology of negation (e.g., articles in Veselinova and Hamari, Forthcoming).
The dynamic component of Croft's typology, the Negative Existential Cycle (NEC), connects these six construction types in a cycle where each construction is the source for another. Elaborating on the NEC, Veselinova (2013Veselinova ( , 2014Veselinova ( , 2016; see also Verkerk and Shirtz, forthcoming) showed that while the diachronic transitions proposed by Croft are indeed attested, other transitions are also involved in the rise of innovative negative existential constructions or in changes to the typological classification of old constructions. It is unclear, however, how widespread these non-NEC transitions are, and how much of the attested cross-linguistic variation in the negative existential domain arises out of, or can be explained by, the transitions in Croft's NEC. This article directly tackles this by asking: how likely is it that the attested variation in the negative existential domain resulted from the transitions that compose the NEC? To do so, we use Bayesian phylogenetic comparative methods to analyze data from 106 Indo-European languages, collected by consulting grammars and published texts, as well as questionnaires filled by language experts. This article tests the likelihood that the NEC is the main set of transitions behind the cross-linguistic variation in the Indo-European negative existential domain, and compares it to the likelihood of other potential sets of transitions. We use the results of this modeling to illustrate our answer to another question: what is the relationship between rational quantitative and statistical modeling approaches and "traditional," analytic approaches in studies of morphosyntactic change and diachronic typology? In a way, we additionally explore the feasibility of phylogenetic diachronic typology; how many languages, or how big of a language family does one need to investigate a complex typological hypothesis? Of course, our answer to these questions is limited to the Indo-European family and to our current sample, as Section 2 elaborates.
The rest of this section further defines the negative existential domain, describes Croft's typology and NEC in more detail, and sketches some of the major issues with the NEC. In Section 2 we turn to describe the data and the methods used in this study, and turn to present some of our findings in Section 3. There, we give an overview of the negative existential domain in the different subfamilies of Indo-European, and illustrate two types of transitions that are not included in the NEC. Our phylogenetic modeling of the Indo-European negative existential domain is presented in Section 4, and our findings are summarized and discussed in Section 5.

Negative Existentials: Definitions and Synchronic Typology
A cross-linguistic study of the negative existential domain requires that it be defined without reference to any language-specific property, i.e., as a comparative concept (Haspelmath 2010(Haspelmath , 2016Croft 2016). Furthermore, the different types of negative existential constructions need to be defined as what Croft (2016) calls "hybrid comparative concept": a combination of functional and formal properties defined without reference to any language-specific properties so they are identifiable in different languages. Following Croft (1991) and Veselinova (2013Veselinova ( , 2014, we view existential constructions as expressing the existence or presence of a particular figure constituent relative to a specific location (the ground) or generally "in the world." Its negated counterpart expresses the fact that a particular figure constituent does not exist or is not present generally "in the world" or relative to some ground location. As this definition does not refer to any language-specific grammatical devices, it can be and has been successfully deployed cross-linguistically, and qualifies as a functional comparative concept.
The definition adopted here is largely compatible with the approach of Creissels (2013Creissels ( , 2019; see also Clark 1978) who views what is usually referred to by the term existential predication as an inverse-locative predication 1 . In this type of predication, the important information is the presence of the figure against some locative ground, rather than the location of the figure, and thus it is the inverse of clauses expressing a predicate location. Creissels' (2019) definition of inverse-locative predication, however, focuses on constructions that are not formally related to constructions expressing locative predication. Instances where the difference between clauses expressing locative and inverselocative semantics has to do with information packaging devices such as word order or topic/focus markers (as in some Indo-Iranian languages; Shirtz 2019), are not included in Creissels' inverse-locative.
This article's focus is on clauses expressing the negative existential domain regardless of their structural similarity to clauses expressing locative predication. We do include here instances where the main distinction between predicate location and existential constructions has to do with the relative order of the figure and the ground or with other information-packaging devices such as articles or so-called topic/focus markers. In this sense, our approach to the negative existential domain, as well as the approach of Croft (1991) and Veselinova (2014), is compatible with Creissels' criteria for inverselocative predication, except for his requirement that it be structurally unrelated to locative predication. Croft's (1991) typology of the negative existential domain is composed of six language types. Their identification rests on comparing the negation marker in clauses expressing the negative existential domain to the negation markers used in standard verbal negation, and marginally on the existence of other negative existential constructions in the same language. The constructions on which the typology rests, then, are bundles of functional and abstract formal properties and as instances of Croft's (2016) "hybrid comparative concepts" are cross-linguistically identifiable. We classify our languages in terms of Croft's (1991) typology, but instead of classifying in terms of language types, we use construction types (as in Veselinova's approach and in Verkerk and Shirtz, forthcoming). This implies that languages may have more than one type of negative existential construction, and that attested constructions may undergo change that is in part independent of other constructions that may be attested in that language.
The six construction types in Croft's typology are divided into three major types, called Type-A, Type-B, and Type-C, and three transitional types, called Type-A∼B, Type-B∼C, and Type-C∼A 2 . In Type-A, the same marker used in standard verbal negation accompanies the affirmative existential marker in clauses expressing negative existence. This was illustrated by the Swedish clause in (1) above, where the standard verbal negation marker inte is used to negate the verbal locative copula finns "be.at." In Type-B, a special marker distinct from the standard verbal negation marker is used in negative existential clauses. This was illustrated by the Turkish clause in (2) above, where the negative existential marker yok is used.
In the intermediate Type-AB one finds instances of Type-A and instances of Type-B that may be diachronically related, each in its own functional niche. This situation is very common in Iranian languages, where an innovative negative copula often emerges from a reduction of the verbal negation marker and the present-tense copula or some other copular element, thus leading to the innovation of a Type-B construction. But a similar reduction does not occur with the past-tense copula, thus a conservative Type-A construction is retained. This is illustrated by (4a-b), from Sivandi. In (4a) the past tense copula, also used in affirmative existential clauses, is negated by the Sivandi Standard Verbal Negation marker na . The negative existential in (4b) is expressed by nund, a negative copula that resulted from the reduction of the standard verbal negation marker na with some other element.
Sivandi (Iranian, Lecoq, 1979: 89, 150) In Type-C, the standard verbal negation marker is used as a negative existential marker, without an affirmative existential marker. This was illustrated above by the Tongan examples in (3a-b), where the standard verbal negation marker 'ikai is also used as the negative existential marker in (3b). In the intermediate Type-BC a special negative existential marker is also used as a verbal negation marker under some circumstances, but other verbal negation markers also exist. That is, the domain of verbal negation includes several markers, one of which also functions as the negative existential marker. This is illustrated by (5a-c), from Darai (Indo-Aryan). The clauses in (5a-b) illustrate two Darai verbal negation constructions: the nai-prefix in (5a) and in (5b) the particle nidzə. This particle is also used in (5c) as the negative existential marker, without the Darai affirmative existential copulas. Thus, in Darai a special negative existential marker is also used as one of the verbal negation markers.
Darai (Indo-Aryan;Dhakal 2012: 134, 134, 137) The sixth type in Croft's typology, Type-CA, includes a negative existential marker which also functions as the verbal negation marker, but is optionally used alongside an affirmative existential marker. This is illustrated here by (6a-b), from Marathi, where the verbal negation marker nahi in (6a) is also used as the negative existential marker in (6b), optionally cooccurring with the existential marker ahe.
Marathi (Indo-Aryan, Croft 1991: 12; his glosses and translation) While using Croft's (1991) typology to classify negative existential constructions is not always straightforward, it often is so, and it is an illuminating taxonomy for the IE data analyzed here. What this paper sets out to test, then, is how well does the diachronic component of Croft's typology, the negative existential cycle (NEC), explain the attested variation in negative existential construction types across the family.

Negative Existentials: Dynamic Typology
The diachronic component of Croft's typology arranges the six construction types in a unidirectional cycle such that each type is the source of one other type of negative existential. The cycle is presented in (7), arbitrarily beginning with Type-A, and cycling through the different types until we return to Type-A.
Croft's dynamization of his typology is appealing. It is simple and unidirectional, and each transition on the cycle is described and illustrated by Croft as an instance of an internal mechanism of morphosyntactic change: reanalysis + actualization 3 or extension. The emergence of Type-B negative existential markers, for example, often involves a reanalysis of the relationship between a negation marker and a copula as a single unit, actualized by a phonological reduction of the two or changes in the distribution of the negated copula. Whether the reanalysis occurs with all copular forms or only with some, it is clear that Type-AB will likely occur at some point in such transitions from Type-A to Type-B, leading to a Type-A > Type-AB > Type-B pathway.
Croft illustrates the transition from Type-B to Type-C with two processes. First, an extension of the negative existential marker into the domain of verbal negation results in competition between the negative existential marker and the standard verbal negation marker. This can occur, for example, when a new main clause construction emerges from the reanalysis of a nonfinite form of the verb and an existential marker. When the two forms are in competition, or when each form is used in its own functional niche, the system is an instance of Croft's Type-BC and if the negative existential marker overtakes the entire domain of verbal negation, we arrive at Type-C. The second type of process involves an extension of the negative existential marker to reinforce the standard verbal negation marker under some conditions. As predicted in Jespersen's cycle (Jespersen 1917;van der Auwera 2009, see also van Gelderen forthcoming), where novel negation markers arise out of older negation-reinforcing elements which end up replacing the older markers, the negative existential marker may be reanalyzed as the main verbal negation marker. At first, the complete replacement will occur only in certain situations, and the system will be an instance of Type-BC, but after the complete loss of the old verbal negation marker, the system will be best classified as Type-C. The negative existential marker of Type-C, then, may optionally combine with the affirmative existential marker, often for information management purposes, innovating a Type-CA construction. When the combination of the old negative existential marker and the affirmative existential marker becomes obligatory, we arrive back at Type-A.

Issues With Croft's Cycle
Croft's proposal, then, includes both a synchronic typology of six types and a set of diachronic transitions between them that results in a cycle. Veselinova (2013Veselinova ( , 2014, 2016 see also Verkerk and Shirtz, forthcoming) further explored the different negative existential types in Slavic and Oceanic languages, and the different transitions between these types as proposed by Croft. Doing so, she identified several issues with the dynamic portion of Croft's proposal.
First, Veselinova notes that languages may have two (or more) negative existential constructions of different types. She illustrates this with the co-existence of Type-B and Type-C in Tahitian and the coexistence of Type-B and Type-BC in Kapingamarangi. Verkerk and Shirtz (forthcoming) further illustrate these patterns with the data from the Eastern Indo-Aryan languages Kupia and Standard Oriya. They also note that pairs of negative existential types may differ in whether their coexistence is at all possible given the way Croft's typology is set up. There is nothing prohibiting the coexistence of multiple Type-A or multiple Type-B negative existential constructions, but multiple Type-C constructions cannot coexist. In such a situation, by definition, the deployment of each particular negation marker will be limited in some way (as there are two verbal negation markers), and hence the two constructions are an instance of multiple Type-BC constructions. This entails that in situations where a Type-B and a Type-C negative existential constructions are found in the same language, predictable changes to one construction, Type-B > Type-BC, entail changes to the classification of the other against the NEC direction (Type-C > Type-BC). Veselinova (2014Veselinova ( , 2016 also identifies transitions that are not represented on Croft's cycle. These include a transition from Type-AB directly to Type-BC (without an intermediate Type-B), potentially documented in Russian and Hawai'ian, and a transition from Type-B directly to Type-C (without an intermediate Type-BC). The directionality of these changes is the same overall direction of the NEC, but a stage is "skipped." These transitions, together with the Type-C > Type-BC proposed by Verkerk and Shirtz (forthcoming), form a set of transitions that are outside of the set of transitions proposed by Croft. This suggests that there is more to the diachronic changes that negative existential constructions undergo than the NEC. The diachronic processes involved are summarized by Veselinova (2016: 155) as follows: "They include (i) subordination processes; (ii) the reanalysis of an external negator into a negator external to the proposition; (iii) a direct inheritance of a construction; (iv) the use of negative existentials with nominalized verb forms." van Gelderen (forthcoming) discusses the NEC in relation to two other negative cycles, the Jespersen Cycle (see the previous section) and the Givón Cycle 4 (Jespersen 1917, Givón 1978; see also van Gelderen forthcoming, van der Auwera et al., forthcoming) as well as the Copula Cycle 5 (Katz 1996), demonstrating how other diachronic processes impact the NEC.
The attested transitions that are not a part of the NEC, by virtue of being in the opposite direction to the NEC or by virtue of "skipping" an NEC stage or two, seem to depend on very specific configurations in the grammar of individual languages, such as the coexistence of Type-B and Type-C in a single language, and these configurations may be cross-linguistically rare. As a result, it seems that transitions that are not a part of the NEC are infrequent. But without a wider cross-linguistic survey, we simply do not know how rare these situations are, and it could very well be that their rarity is a result of diachronic instability.

MATERIALS AND METHODS
The previous section suggests that while the transitions in the NEC seem to be common, other transitions are attested as well, including transitions in the same general direction of the NEC and transitions in the opposite direction. The question, then, is how much of the attested variation in the negative existential domain does the NEC explain? This article answers this question for the Indo-European language family. This section describes the data, the theoretical assumptions, and the methods used in this article.

Data
The data used in this article is the classification of negative existential constructions in 106 Indo-European languages. The Indo-European language family is well suited for a study such as the one done here. First, it is a large family that includes many subfamilies with a wide geographic dispersal and deep historical records. While not a requirement, this dispersal raises the likelihood for variation in the typological classification of negative existential constructions. This variation is required for a meaningful quantitative testing of the NEC. In families with little to no variation to explain, it will be difficult to reject transition sets based on their low explanatory power. Further, the documentation of Indo-European is quite extensive, and while there are still several lacunae in the documentation of the family (e.g., the Indo-European languages of Pakistan and the Pamir region), many branches and sub-branches are well documented.
The data collection for this article relied on three types of data sources. First, as with many typological surveys, we made extensive use of published grammatical descriptions. The domain of existential clauses, affirmative or negative, however, is often not directly mentioned in such sources. This may be because of their marginal nature or low frequency, or because they sometimes do not have any unique or unpredictable grammatical properties (e.g., in compositional Type-A constructions). To have as wide a coverage as possible, then, we also used translation questionnaires and analyzed published textual data. The questionnaire is based on the one used by Veselinova (2014; see Supplementary Information S1) and was filled by language experts. It includes questions about affirmative and negative verbal clauses, affirmative and negative existential clauses, and other types of nonverbal predication. We also made use of published textual data (not necessarily computer-readable corpora) which accompanied documentation projects. This was required when the reference or sketch grammars did not include an explicit discussion or illustration of existential and negative existential predication, but were clear and enabled us to go through textual data published as a part of a documentation project.

Typology-Historical Morphosyntax-Phylogenetic Modeling?
This article is situated at the intersection of linguistic typology, historical morphosyntax, and phylogenetic modeling. This requires the article to adhere to the main assumptions of each of the three fields. As this is a typological study, we approach the negative existential domain here as a functional comparative concept, and define it without reference to any language-specific construction. Further, Croft's definition of the six negative existential types are instances of hybrid comparative concepts, based on both form and function (Croft 2016, see Section 1.1).
More controversial is the relationship between the more traditional, analytic approaches to historical morphosyntax and newly adopted statistical, phylogenetic approaches. These two approaches are often viewed as competing, or even contradictory in their assumptions (see, for example, replies to Dunn et al., 2011 in Linguistic Typology 15.2;especially Dryer 2011 andPlank 2011). We argue, however, that when it comes to the study of morphosyntactic change, the two do not contradict each other, but rather highlight different aspects of language change. As such, they are best viewed as complementing each other by answering slightly different questions so that each of them may inform the hypotheses and the work done in the other (see again Linguistic Typology 15.2 for Levinson et al., 2011 andCroft et al., 2011; as well as Gray 2012 andDunn 2014: Sect. 5.2). We use the domain of the negative existential in Indo-European to illustrate this point.
The field of historical morphosyntax in general, and morphosyntactic reconstruction in particular, has been rife with controversy about the very plausibility of its goals. The identification of the mechanisms of morphosyntactic change (Harris and Campbell 1995) and the introduction of a constructional interpretation for morphosyntactic change (Barðdal and Eythórsson 2012;Barðdal 2013;Barðdal and Gildea 2015), however, enable an explicit statement of the methods used and assumptions required in the analysis of morphosyntactic change. These assumptions include the identification of cognate constructions using a set of principles that are parallel to those used in the identification of lexical cognates, and the identification of the plausible mechanism of change involved (see Gildea et al., 2020 for a detailed survey and discussion). These principles are also applied in diachronic typology (Bybee 1988;Bybee et al., 1994;Hendery 2012;Sansò 2017).
The goal of phylogenetic comparative modeling of the type pursued here, focusing on change in a typological variable, involves the estimation of the likelihood of a set of transitions, or historical changes between construction types, in a set of observed data (Pagel 1999). This estimation depends on the topology of the family tree, which should be arrived at independently (e.g., using the Comparative Method), and the length of its branches, estimating time elapsed since the diversification of two languages (Pagel 1999: 878, Dunn 2014.2). The cognate status of the attested constructions, as well as the specific mechanisms involved in the rise of each of these constructions, matters less for such phylogenetic modeling. That is, the modeling pursued here treats transitions between construction types with no regard to the actual process of change "on the ground." Several different processes can often lead from one construction type to another, but which of these actually occurred is not a part of the model. The fact that analytic, "traditional" methods in historical morphosyntax and phylogenetic comparative modeling of morphosyntactic change highlight different aspects of the data may be taken to suggest that these are competing methods. We however believe the exact opposite: the fact that different aspects of the historical record are highlighted by these methods allows them to complete and inform each other. Croft's NEC and Veselinova's critique of the NEC were arrived at using analytic morpho-syntactic methods. Both were arrived at without taking into account a specific family tree topology (although family relations must have been taken into account implicitly), and without testing the NEC against other plausible pathways using quantitative tools. Both of these obviously motivate and inform the current study. Testing how much of the cross-linguistic variation the NEC can explain will either fortify it as the main set of diachronic transitions in the negative existential domain, or propose other (sets of) transitions active in this domain that can then be further explored by a more direct analysis of language data.

Phylogenetic Comparative Methods
We model the type of negative existential strategy that each language in our sample has in terms of an explicit phylogenetic process, i.e., as the outcome of evolutionary processes that take place on the branches of a phylogenetic tree. The phylogenetic tree set is given, we are not inferring phylogenies, but rather using them to do quantitative diachronic typology: testing an influential hypothesis using phylogenetic methods. Phylogenetic comparative models have been used to estimate what typological strategy the ancestors of sampled languages must have had (Maurits and Griffiths 2014), and the rates of evolutionary change (Cathcart et al., 2020). We will focus on which transition parameters are most relevant for explaining the distribution of strategies attested in our sample (see also Dunn et al., 2017). Here, we test whether the transition parameters associated with the NEC are essential for explaining the diachrony of the distribution of negative existential strategies in the current sample of 106 Indo-European languages.
Doing this requires three components: data on negative existential strategies, a tree sample of phylogenies of the languages under investigation, and a set of models, grounded in a particular way of thinking about evolutionary change. This section describes the sample of phylogenetic trees we use, and sketches the relevant model of evolutionary change. Our dataset is covered in Section 3. More details regarding specific phylogenetic comparative testing are given in Section 2.3.
Since none of the currently available Bayesian tree sets (Bouckaert et al., 2012;Chang et al., 2015;Heggarty et al. in review) sample all of the languages in our dataset, we use trees from Glottolog (Hammerström et al., 2014) which have been given necessary branch lengths by Dediu (2018). Dediu (2018) takes cladogram-like trees from four different sources, and adds branch length, vital for phylogenetic comparative analysis, using nine different methods. We describe in Supplementary Information S2 how we opted for two of these trees. We used the function multi2di() in the R package ape (Paradis et al., 2004; R Core Team 2020) to create 250 trees in which the polytomies (nodes in the tree which lead to more than two clades or languages) present in Glottolog were resolved in a random fashion. Subsequently, branch lengths were added to these newly created branches in a random fashion corresponding to the distribution that the branch lengths for each of these trees have. This resulted in a sample of 500 phylogenetic trees.
There are different models for the evolution of different types of characters (Pagel 1999;Meade and Pagel 2019): binary characters (a language either has a characteristic, like having one or more click consonants, or it does not), continuous characters (a real number, such as the entropy of object-verb word order in a parallel corpus; Levshina et al. in review;or Greenberg's 1960 morpheme-word ratio), or multistate characters in (Comrie's 2013 WALS chapter on the alignment of case marking of full noun phrases, a language may have one of six possible alignment patterns). In this study, we are concerned with a multistate character with exactly six states that detail the interaction between existential and standard negation (see Section 1).
The standard model to account for multistate characters is the continuous-time Markov process of character evolution (Pagel et al., 2004). This model describes the probability of change between states of a character (here, negative existential strategies) in terms of a set of transition rate parameters. Here, "continuous-time" implies the character can change its state at any instant of time rather than at fixed intervals; "Markov," from "Markov chain," indicates that the probability of changing from one state to another depends only on the current state, and not on any earlier states. The changes that take place in the continuoustime Markov model are summarized using a transition rate matrix or a Q matrix, where the individual transition rate parameters are designated by q, followed by codes for two states. Figure 1 illustrates this matrix for the negative existential domain. It is the set and values of transition rate parameters captured by the Q matrix, as well as the likelihoods that are associated with different states at the internal nodes of the tree, that are tracked during analysis.
Because there are six types in Croft's (1991) typology, there are 6×6-6 30 transition rate parameters. States cannot change into themselves (these dependencies are marked in Figure 1 by "-"), hence these represent the probabilities that the state stays the same. In the Q matrix, the diagonal "no change" probabilities and the off-diagonal transition rate parameters (q's) sum to zero. Croft's (1991) NEC proposes a diachronic typology using only six of these changes between types, as indicated by Figure 1. Croft's (1991) NEC is very ambitious given the possible transition-rate parameter space: modeling change between six types using only six diachronic pathways between types.
The large number of types and correspondingly large number of transition rate parameters, together with the rarity of some types (see Section 3), pose a practical problem. The "one in ten" rule in statistics also applies to phylogenetic comparative analysis, i.e., having ten data points (species or languages) per free parameter is an aim during data collection, despite the fact that actual sample size is reduced through phylogenetic dependencies (Mundry 2014). This implies needing a sample size of 300+ languages to run the model in which all 30 transition rate parameters are included. Ideally, construction types would be distributed evenly across that ideal sample, but this is not realistic (see Section 3.1, Croft 1991; Veselinova 2016). To have a reliable estimate, for instance, whether change to Type-CA is more likely to come from Type-C (as predicted by the NEC) or any other type, we would probably need even more data. In other words, our Indo-European dataset is still too small, and the distribution of types is too skewed, to comprehensively test Croft's (1991) NEC. However, we will try regardless of these issues and report on the results in Section 4. 6 We aim at model optimization that will 1) test which set of transition-rate parameters explains best the distribution of negative existential types in our Indo-European dataset and 2) compare the best fit models to the NEC. The Bayesian model used here allows several options for model optimization. The first option we have is simply excluding certain transition-rate parameters manually. This is also how we test the NEC model, by excluding the 24 transition-rate parameters in black and red typeface in Figure 1 by setting them to zero. The second option, Reverse Jump MCMC (RJ MCMC, Green 1995;Pagel and Meade 2006), automatically turns on and off transition-rate parameters while at the same time estimating and reducing the number of different rates. In the posterior models, transitionrate parameters that do not contribute to the model are excluded and the number of individual transition-rate parameters is typically reduced such that a small number of rates is shared across parameters, optimizing the model and making it "more elegant." Excluding transition rate parameters manually and doing RJ MCMC can also be done at the same time.
Again, the large number of types and correspondingly large number of transition rate parameters poses a practical problem. If, for example, we want to compare a model with ten transition rate parameters (perhaps the six FIGURE 1 | Q matrix of the six states negative existential constructions may be in. Type-A has been coded as 1, AB as 2, B as 3, BC as 4, C as 5, and CA as 6. The first number refers to the state that is left (rows), and the second number refers to the state that is entered (columns). Thus qij is the transition parameter from Type-i to Type-j (example: q12 is Type-A > AB; q21 is Type AB > A, etc). The changes between the states described by Croft's (1991) Negative Existential Cycle have been marked in green color. The set of opposite transitions is marked in red color. Note that, with the exception of q61 and q16, changes from smaller to bigger numbers (for example, q25; Type-AB > C) designate changes in line with the direction of the NEC, while changes from bigger to smaller numbers (for example, q52; Type C > AB) go against the direction of the NEC.
6 One solution to the large transition rate parameter space problem that we have tried for an earlier version of this paper is to exclude Croft's (1991) transitional types AB, BC, and CA from the dataset and model. The three possible types are then A, B, and C, AB would be re-coded as A&B, BC as B&C, and CA as A&C. This reduces the transition rate parameter space to 3 × 3-3 6; and three of these transition rate parameters, i.e., change from A > B, B > C, and C > A are implied by Croft's (1991) NEC; the other three transition rate parameters, A < B, B < C, and C < A, are the reverse of Croft's (1991) NEC. However, while this effectively eliminates the large transition rate parameter space problem, it is not viable because 1) it does not do justice to the data, for example, split languages have to be coded as A&B&C; 2) because of the resulting prevalence of having two states (we have especially many AB languages which would be coded as A&B), including/excluding specific transition rate parameters is not informative, as any combination of transition rate parameters becomes equally likely. NEC parameters plus four more) with ten random transitionrate parameters picked out of the 30 in our Q matrix, we have to face the fact that there are 30,045,015 possible sets of 10 parameters out of 30 transition-rate parameters. It is very likely that some combinations of ten parameters fit our data better than others; but identifying these combinations is exactly the problem. Testing so many models is not feasible: a single model takes four to seven hours to run on a normal desktop computer, and the number of models grows exponentially when we add models with 11, 12, 13 transition-rate parameters out of our set of 30. To improve our chances of finding the model with the best fit, then, we have to rely on RJ MCMC. In Section 4, we introduce the models we tested and discuss their fit. We used BayesTraits V3.0.2 (Meade and Pagel 2019) to conduct phylogenetic Bayesian MCMC analysis, specifically its component MultiState (Pagel, Meade, and Barker 2004). We construct various models we want to test, focusing on the transition-rate parameters; these are covered in Section 4. We conducted a single MCMC analysis for each model, which was run for 2×10 7 iterations, with a burn in of 1×10 7 , sampling every 10 5 iteration, resulting in a sample of 1,000 posterior estimates. Convergence was assessed by checking the absence of a correlation between the posterior likelihood and the iteration number. Lack of autocorrelation between samples was assessed visually. When used, Reverse Jump was used on all transition-rate parameters, with a default exponential prior with mean 50. When Reverse Jump was not used, the default uniform prior with distribution 0-100 was used for the transition rate parameters. BayesTraits V3.0.2's built-in stepping stone sampler was used to estimate log marginal likelihoods after the MCMC analysis was concluded. The log marginal likelihoods were used to assess the fit of the various models in Section 4.

NEGATIVE EXISTENTIALS IN INDO-EUROPEAN: A TYPOLOGICAL SURVEY
This study surveys the expression of the negative existential domain in 106 Indo-European languages. Data on 42 languages come from Verkerk and Shirtz (forthcoming), data on 13 Slavic languages come from Veselinova (2014), and data on 51 additional languages are added to this article (see Supplementary Information S3). We aimed to sample as extensively as possible, but were constrained by both available resources and time during a global pandemic. This section first briefly summarizes the results of our survey, and highlights some noteworthy areal tendencies (see Verkerk and Shirtz, forthcoming for a more detailed discussion): the relative typological stability of the negative existential in some branches or areas and its relative instability in other branches or areas. Then, we provide a brief historical and phylogenetic overview of the data, and following these two sections we discuss morphosyntactic innovations in the negative existential domain that do not involve a change in the typological classification of a construction and innovations that lead to transitions that are not included in the NEC.

General and Areal Overview of the Typology
The grammatical expression of the negative existential domain in Indo-European is varied, with each of Croft's six types attested somewhere in the family. This variation is not homogenous and some types are quite frequent while others are rare. Furthermore, the variation is not equally distributed across Indo-European, and some families exhibit a rather uniform typology of the negative existential domain while other families are diverse. The diagram in Figure 2, constructed based on the NEC itself, indicates the raw counts of each attested construction type. 7 In our Indo-European sample, Type-A (37.5%) and Type-AB (31%) are the most common types. The biggest difference between Figure 2 and the world-wide surveys in Veselinova (2014Veselinova ( , 2016 is that in our Indo-European sample, Type-AB (current paper: 31%; Veselinova: 8.9%) is far more common than type-B (current paper: 12.5%; Veselinova: 29.7%), but our Indo-European sample resembles Veselinova's findings for Berber and Uralic. Aside from these differences, the Indo-European data, just like Veselinova (2014Veselinova ( , 2016 worldwide sample, mostly confirms Croft's (1991) remark that types A and B are more common than Type-C, and the transitional types AB, BC, and CA are uncommon, with the caveat that Type-AB is quite common across Indo-European. However, we found six CA languages, much more than the two instances in Veselinova (2016: 150) world sample and her family-based studies (1 CA language out of 109 languages).
The areal distribution of the different construction types is illustrated by the map in Figure 3. It illustrates how the unequal distribution of construction types in general is magnified when focusing on certain areas and certain families. For a first impression, one can simply contrast the relative color FIGURE 2 | Frequency diagram of attested types in our 106 language sample. Note that six languages have more than one type; Oriya is Type-A and Type-BC, Kumzari is Type-A and Type-C, Kupia, Chitpavani Goan Konkani, Goan Konkani, and Varhadi-Nagpuri are Type-B and C. 7 It is not unheard of that languages have more than one type, Veselinova (2016: 154) discusses 9 other cases. Interestingly, the types do not have to be (sometimes cannot be) consecutive in the NEC, showing that different parts of the negative existential domain in an individual language can undergo different transitions.
Frontiers in Communication | www.frontiersin.org September 2021 | Volume 6 | Article 661862 8 uniformity across Europe, especially in the Romance and Germanic-speaking areas, to the diversity found in the Indo-European languages of Western, Central, and South Asia.
The relative uniformity across the European part of the map in Figure 3 is the result of two larger Indo-European families, Romance and Germanic, exhibiting little to no typological variation in the expression of the domain. The Romance languages in our sample uniformly negate existential statements using the same negation marker used in standard verbal negation. Thus, they consistently exhibit a Type-A negative existential construction. This is illustrated by (8a-b) below from Piedmontese (Turinese): Piedmontese (Turinese: Romance; Emanuele Miola p.c.) Typological uniformity is also attested across Germanic, where alongside Type-A constructions (illustrated in (1) above, from Swedish), existential statements may also be negated using a negative indefinite pronoun or determiner. In (9), from Swedish, the existential statement is negated by the Swedish Negative indefinite pronoun inget, and in the English translation of the example, the statement is negated by the English marker no.
Swedish (Bordal 2017: 6; their glosses and translation) The uniformity of the Romance and Germanic families stands in contrast to the diversity attested in Iranian and Indo-Aryan, and to a lesser degree in Slavic. The factors involved in this difference may include borrowing, diachronic replication, substrate factors, or universal tendencies (e.g., Nichols 1992). It is beyond the scope of the current paper to argue which of these factors led to the typological uniformity of negative existential construction in Romance and Germanic on the one hand, and to the typological variation in Slavic, Iranian, and Indo-Aryan on the other hand. For now, suffice it to mention that the languages of Western Europe form a Sprachbund (e.g., Haspelmath 2001; van der Auwera 2011), and propose that the uniformity across Germanic involves some sort of diachronic replication. Finally, note that the pattern whereby the Iranian and Indo-Aryan families exhibit much more typological variation than the Germanic and Romance families is not limited to the negative existential domain. A September 2021 | Volume 6 | Article 661862 9 similar pattern is attested, for example, with the alignment of core arguments, where Iranian and Indo-Aryan are very diverse (e.g., Haig, 2008;Verbeke, 2013), while the Indo-European languages of Western Europe are rather uniform.

Historical and Phylogenetic Overview
Not surprisingly, the uneven areal distribution of variation in construction types goes hand-in-hand with uneven distribution across different Indo-European subfamilies. We illustrate this on a randomly chosen phylogenetic tree from our tree sample in Figures 4 and 5. As described in the methodology, the trees we used for phylogenetic comparative analysis were built by making the polytomies binary in a random way, and assigning branch lengths to these newly created branches on the basis of the distribution of existing branch lengths from Dediu (2018). This leads to sometimes unrealistic higher order groupings, such as the one we find here relating Hittite and Celtic. We do not argue that this is how the Indo-European languages actually evolved; this is simply one of many possibilities that was selected for display purposes only. Given the size of the sample and tree, we split it such that the non-Indo-Iranian part of the family is displayed in Figure 4, and Indo-Iranian is displayed in Figure 5. The pie plots on the internal nodes of the tree represent marginal ancestral state reconstructions conducted in the R package corHHM (Beaulieu et al., 2013;R Core Team 2020). These are again illustrations on a single tree; the analyses were conducted on the full tree sample and reconstructions will differ across trees (see Section 4). A simple parsimony reading can be misleading. For example, in Bihari (Chitwania Tharu, Ranna Tharu, Darai, Sadri, and Bhojpuri), Darai and Sadri have not changed Type-A > Type-BC or Type-B, but rather, the Tharu languages and Bhojpuri are likely to have finished the NEC cycle and reached Type-A again. Hence, the reconstructions are partly realistic, partly a consequence of the tree structure coupled with gaps in the diachronic record as intermediate stages are often not present in the dataset (or not recorded at all).
While much of this paper focuses on transitions that are not a part of the NEC, it should be mentioned that many instances of FIGURE 4 | An overview of negative existential construction types in non-Indo-Iranian languages, overlaid on a modified Indo-European Glottolog tree (Hammarström et al., 2014;Dediu 2018 Figures 4 and 5 can be deceptive due to missing intermediary stages and languages we did not sample. The transition from Type-A to Type-AB can be clearly seen in Palula (Indo-Aryan). The verbal negation marker in Palula is na (Liljegren 2016). Existential predicates may be negated by na, or by the special negative existential copula náinu, a reduction of na NEG + hínu COP (Liljegren 2016: 413). The Palula data, then, illustrates the NEC's A > AB transition. Two Kurdish languages in our sample illustrate change from Type-AB to Type-B. In Mukri (Central) Kurdish, standard negators are nā-for present tense, ne-for past tense (Opengin 2016: 74). Negative existential strategies show another tense-based split, with standard negation used in the past tense, and in non-past tense, the negative copula negation nī is used. Bahdini Kurdish (Kurmanji) is Type-B, with standard negation using the clitic or prefix na-, ne-, and a special negative existential copula tun- (Thackston 2006).
While there are several BC languages in our sample, it is not easy to find a clear example of the type B > BC transition. One tentative example is Darai (Type-BC), where the special negative existential marker nidze is used as one of two nonexistential negation markers (Dhakal 2012). Closely related Sadri has a Type-B construction with the potential cognate special negative existential verb nʌkh (nʌkhe, Jordan-Horstmann 1969). Perhaps an earlier stage of these languages was Type-B, with Sadri being conservative and all other languages in this subgroup being innovative (see below on the CA > A transition). A tentative example of the BC > C transition can be found between Dhivehi and Sinhalese. Dhivehi has a special negative existential copula net (<OIA nāsti, Fritz 2002) and is Type-B, Sinhalese is Type-C, with the free standing negative existential naeae that also functions as postverbal predicate negator (Chandralal 2010).
The transition between type C > CA is attested twice in the Indo-Aryan Southern zone. In Standard Goan Konkani the verb-like FIGURE 5 | An overview of negative existential construction types in Indo-Iranian languages, overlaid on a modified Indo-European Glottolog tree (Hammarström et al., 2014;Dediu 2018 negator nay is combined with existential predicates (Ghatage 1966). In closely related Chitpavani Goan Konkani the likely cognate nãy ∼ naĩ is used as standard negator and as special negative existential, without an existential predicate (Bhide 1982). Similarly, Marathi (Type-CA, the special negative existential is nāhi, Croft 1991) can be contrasted with Katkari (Type-C), where the cognate negator nahĩ ∼ nay does not yet combine with an existential verb (Kulkarni 1969). Evidence for a CA > A transition is found in Chitwaniya Tharu (Indo-Aryan), for example, where the verbal negation marker hoyne, a combination of the old h-copula and a negation marker, is used to negate existential statements alongside an innovative, and obligatory, existential verb. This negation marker, however, is still deployed as a nonverbal negative copula, without a synchronic verbal copula, in some conservative nonverbal predication constructions, where it is often followed by an emphatic clitic marker (Paudyal, 2014). Varli (Abraham and Abraham 2012) has likewise undergone the CA > A transition, as the negator nahĩː (likely cognate with Marathi (CA) nāhi and Katkari (C) nahĩ ∼ nay) can no longer be used without an existential predicate (as is still optional in Marathi).

Illustrations of Innovations That do Not Involve a Change of Construction Type
The typological stability in Romance and Germanic, as well as in some sub-branches of other families, may lead one to believe in some extreme conservatism in the verbal and existential negation in these families. Reality, however, is more complex and across these families there are several instances of innovation in these domains, as well as innovations in the domain of existential constructions in general. These innovations do not lead to a change in the typological classification of the expression of negative existence in these languages.
Across Romance, innovations are attested in the expression of existential predication and in the expression of negation. Existential predicates in Catalan, illustrated in (10), are expressed by a combination of the locative adverbial clitic hi "there" and a third person form of haver "have." In Romanian, on the other hand, existential predicates can be expressed by a se găsi, the middle form of the verb "to find," as illustrated by (11), or by the verbs a exista "to exist" or a fi "to be." Catalan (Romance, Wheeler et al., 1999: 460) Romanian (p.c. Andreea Calude) As the innovative expression of existential predication involves verbs in both languages, it is only natural for them to be negated by the standard verbal negation marker (at least initially). Thus, an innovative existential verb in a language which already had a Type-A negative existential construction, the common situation in Romance, leads to a novel negative existential construction without a change in the typological classification of the domain.
The use of a negative indefinite pronoun or determiner to negate existential predication in Germanic also hides instances of innovation. This involves the rise of innovative negative indefinite forms, such as German kein, Dutch geen, and Swedish ingen when compared to English no. While these forms are related, they involve different types of syntactic and lexical innovations (German kein from *nih "neither" and *aina-"one"; Dutch geen from neh "and not" and ein "one" (Philipa et al., 2003(Philipa et al., -2009), Swedish ingen from einn "one" +-gi, privative suffix; the use of negative indefinite pronouns across different nonverbal predicates differs quite radically, see Verkerk and Shirtz, forthcoming). Thus, once a negative existential construction with a negative indefinite pronoun as the negation marker exists, an innovation in the domain of this marker would not alter the typological classification of the construction itself.
Similar innovations can be found in Greek, where innovations in the expression of negation occurred from time to time (e.g., Kiparsky and Condoravdi 2006). Across Indo-Iranian, locative verbs have often been co-opted into existential predication, and by extension also negative existential predication. This often results in an innovative Type-A construction that may or may not lead to a change in typological classification.

Illustrations of Innovations That Are
"Outside" the NEC Figures 4 and 5 use a phylogenetic tree to illustrate the synchronic variation in negative existential constructions in Indo-European. More than that, these figures also illustrate one proposal for the reconstruction of the type of negative existential in ancestral nodes on the tree. A closer look at the tree would suggest that there are several transitions that cannot be explained in terms of the NEC. We describe such transitions in this section on the basis of our own analysis, as the tree can mislead through data gaps. These transitions are of two types. First, we find transitions where a development within the domain of negative existence itself leads to a change in the classification of a negative existential construction. Second, we find transitions that involve innovations that occur outside of the negative existential domain but affect it. This includes innovative negative existential constructions entering the domain and innovations in the realm of verbal negation.
The first type of transition was illustrated by Macedonian and Bulgarian, where Veselinova (2014) shows a transition from Type-A (as illustrated by Old Church Slavonic data) directly to Type-BC in Bulgarian and Macedonian, without moving through the intermediary types. Another example can be found in Kumzari where we find a split between Type-A and Type-C, the latter evolving directly from an older Type-A construction. Verbal negation in Kumzari is expressed by a post-verbal na (van der Wal-Anonby 2015: 211-213; see also the main clause in (12a)). Affirmative existentials are expressed by clauses containing the figure NP, an optional locative ground, and a copula. Now, the source of the Kumzari Frontiers in Communication | www.frontiersin.org September 2021 | Volume 6 | Article 661862 enclitic copula is the Old Iranian *h-copula, and this copula underwent a great deal of phonological reduction that resulted in its complete deletion in the Kumzari 3SG form (van der Wal-Anonby 2015). This resulted in clauses such as the subordinate clause in (12a) or the unipartite clause in (12b) where only the figure NP is expressed. The negative existential is simply expressed by clauses composed of the figure NP followed by the verbal negation marker na.
Kumzari (Iranian, van der Wal-Anonby, 2015: 184;164;140) The Kumzari construction illustrated in (12c), then, is an instance of Type-C, as the standard verbal negation marker is used as the sole marker of existential negation. This Type-C construction did not arise out of a previous Type-BC, but given the source of the enclitic copula in the Old Iranian *hverbal copula with its subsequent extreme phonological lenition, arose out of an older Type-A construction (other instances of a Type-A > Type-C transition are illustrated by Croft, 1991).
The second type of change involves innovations outside the negative existence domain that affect it. One type of such an innovation has been proposed by Verkerk and Shirtz (forthcoming), where the rise of an innovative verbal negation marker may change the classification of Type-C negative existential construction in the opposite direction to the NEC into Type-BC, as now only one of several verbal negation markers is used in clauses expressing negative existence. An opposite scenario also occurs, where a loss of an older verbal negation marker may affect the classification of an existing construction. In Eastern Indo-Aryan, for example, there were (at least) a preverbal negation marker and a postverbal negation marker. In Kupia, the older preverbal negation marker was lost, but was fossilized in some lexical verbs including "not know," "be unable," and in the negative copula nenj-which is still used in negative existential clauses (Christmas and Christmas 1973: 310). Thus, the Kupia negative existential construction with nenj-went from being a Type-A construction to being a Type-B construction, without an intermediate Type-AB. Similarly, changes in the domain of verbal negation in Assamese seem to have led to a change from a previous Type-BC to a Type-B construction. As these constructions coexisted with a conservative Type-A construction, Assamese is now classified as Type-AB.
Another type of innovation involves a novel verbal existential marker that is negated by the standard verbal negation marker. Thus the novel negative existential construction in this case is a Type-A construction. This was illustrated above by Catalan and Romanian existential constructions. Such innovations are common in Iranian, where verbs translatable as "be.at" evolved to express predicate location and existence. This is illustrated by the Sivandi example below, where dar "be.at" is used as the existential copula and is negated by the Sivandi verbal negation marker, a preverbal na. When such a construction evolves in a language that already has a Type-B construction, and each construction settles in its own functional niche, the result would be Type-B > Type-AB transition, in the opposite direction of the cycle.
Sivandi (Iranian, Lecoq 1979: 15) We conclude this section with an overview of changes discussed here and elsewhere that do not fit the NEC: Kumzari, see also Croft (1991)  The transitions that form the NEC, then, may account for many attested changes in the Indo-European negative existential domain, but there are other attested or plausible changes that are not a part of the NEC set of transitions. This was already mentioned in Croft's original description of the NEC (1991), developed by Veselinova (2013Veselinova ( , 2014, and further systematized here.

PHYLOGENETIC COMPARATIVE MODELING
This section deals with diachronic change in strategies of negative existentials as modeled on the branches of a phylogenetic tree set. There are two main options for designing the set of transition rate parameters: leaving both the selection of parameters and the estimation of their (communal) rate to the RJ MCMC analysis, and selecting transition rate parameters to be included or excluded manually. We use both approaches, and additionally try to combine them. The RJ models and results are presented in Section 4.1, the manual models are presented in Section 4.2. All code and results can be found in Supplementary Information S5.

Reverse Jump MCMC Models
In this section, we report on a set (or rather, a chain) of Reverse Jump MCMC models where we exclude, step by step, transition rate parameters that were infrequently turned on in the set of posterior models. Our starting model is a full RJ model (exponential prior, mean 50), which generates posterior models by allowing transition rate parameters to be excluded and/or their rate to be set equal across parameters. The posterior distribution of models gives a sense of which transition rate parameters are of vital importance (i.e., often turned on) and whether rates differ across rate parameters. Note that BayesTraits (Meade and Pagel 2019) does not allow for a large number of free transition rate parameters (30) to be estimated without RJ. Similarly, too low an included number of free transition rate parameters, especially where one of the states ends up unreachable, are also disallowed by BayesTraits. The first model, RJ.FULL, is an RJ model with no prior restrictions on the transition rate parameters, so all transition rate parameters (also called qs, see Figure 1 for the Q matrix) can be turned on. The number of times in which each transition rate parameter is indeed turned on in RJ.FULL's set of posterior models is given in Figure 6.
In 989/1,000 posterior models, there was a single rate estimated for all transition rate parameters that were turned on. This implies that there is little evidence for different types of changes occurring at different rates in models; however, this may also be due to the lack of constraints on excluded parameters. We could hypothesize, for example, that change from transitional types AB, BC, CA to nontransitional types A, B, and C would occur at a faster rate than vice versa, but there is no evidence for that in this model. In the RJ.FULL model, the mean transition rate is 0.39 (median 0.32). Figure 6 shows that very few transition rate parameters are consistently turned off in the posterior models (only q15 (A>C) and q25 (AB>C)). Conversely, very few transition rate parameters are consistently turned on in posterior models (q31 (B>A) has the largest frequency, featured in over 80% of posterior models). Most transition rate parameters are turned on about 40-60% of the time. Hence, we do not observe a clear pattern of which transition rate parameters are relevant and which are not. There are several explanations for this pattern: 1) there is a multimodal distribution of well-fitting transition rate parameters that is dependent on the characteristics of the phylogenetic trees; 2) dependencies between transition rate parameters, such that they replace each other across models; we could, for instance, imagine that in some model, q12 (A > AB) and q23 (AB > B) are turned on, while in another model, these two are not needed, but only q13 (A > B) is turned on; 3) the sheer amount of transition rate parameters allows for a multitude of likely models, all of about an equal good fit. Unfortunately, it is difficult to tease apart the cause for this mixed pattern of turning on and off transition rate parameters. The lower and upper bounds for the number of transition rate parameters that were turned on for the RJ.FULL model were 9 and 26, compared to the six transitions of the NEC. However, there are millions and millions of options to create models with 9-26 parameters, so this information is not useful. It would further be an immense task to find out if there are indeed correlations between characteristics of the trees and the transition rate parameters that are turned on or off in subsets of the posterior models. Therefore, we have to conclude that RJ.FULL does not immediately point us toward an elegant, clear model of diachronic change. Nevertheless, we find the model informative because 1) it demonstrates how Reverse Jump MCMC models work in the context of character evolution; FIGURE 6 | Frequency of transition rate parameters being turned on in the RJ.FULL model. The y axis marks the number of posterior models in which a transition parameter is turned on, with the horizontal line marking being turned on in 50% of the posterior. Type-A has been coded as 1, AB as 2, B as 3, BC as 4, C as 5, and CA as 6. Hence, q12 refers to change A > AB; q13 refers to change A > B, etc. See Figure 1 for the Q matrix.
2) it demonstrates the intrinsic difficulty in modeling a feature with six states; and 3) RJ.FULL gives us at least some sense of which transition rate parameters are relevant and which are not, albeit limited. In the remainder of this section, we build on RJ.FULL.
The way forward is to manipulate the set of transition rate parameters, such that we exclude from the prior those parameters that were not turned on in RJ.FULL often, in the hope of making the model more decisive and obtain a higher log marginal likelihood (log mLh). The mLh of a model "is the integral of the model likelihoods over all values of the models parameters and over possible trees, weighted by their priors" (Meade and Pagel 2019: 14). It is the main mechanism used by BayesTraits (and other software) to assess model fit. The mLh is computationally expensive, and is therefore estimated using stepping stone sampling (Xie et al., 2011) in BayesTraits, which provides an estimated log mLh. We can compare the log mLh of the two models, and see if there is evidence for a significantly better fit of the better fitting model by calculating log Bayes Factors (BF). The better fitting model is the one with a higher log mLh (because log likelihoods are negative, it makes sense to think about the better-fitting model being the one that is closer to zero).
(15) Log Bayes Factor 2(log marginal likelihood better fitting model − log marginal likelihood worse fitting model) The log Bayes Factor can be interpreted such that a BF > 2 constitutes positive evidence against the null hypothesis, the bigger the BF the more convincing the evidence (Kass and Raftery 1995: 777, their two log e (B 10 )). 8 We build smaller RJ models by excluding parameters that were not often turned on in RJ.FULL. In the first model, RJ_17PAR, all parameters that were absent in 50% or more of the posterior models of RJ.FULL (q12,q14,q15,q16,q23,q24,q25,q26,q51,q52,q53,q54,and q62) were excluded in the prior. What follows are models in which further infrequent transition parameters are excluded, each constructed on the basis of the preceding one. Model names such as RJ_17PAR, RJ_16PAR, etc. are constructed for these models such that they indicate the number of transition rates PARameters that is included. The models are displayed in Figure 7, and the results of all of them are reported in Table 1.
As indicated by the ordering in Figure 7, the RJ models were calculated consecutively, i.e., model RJ_16PAR was constructed by excluding an infrequent transition rate parameter from RJ_17PAR, etc. The figures that detail how often the transition rate parameters are turned on in each model, which served to frame this successive exclusion of parameters, are included in Supplementary Information S5. Note that many different choices could have been made in this successive exclusion of transition rate parameters and that we have not exhaustively sampled the set of possible models in any sense. 9 Doing this, we ultimately arrive at RJ_9PAR_CA and can no longer exclude any parameters from the RJ model (models with eight parameters out of the nine in RJ_9PAR_BC, _C, or _CA do not run). Hence, any RJ model has a minimum of nine transition rate parameters. RJ_9PAR_CA has two epicenters of change: cyclical change between Types-A, AB, and B, and then chance centering around Type-CA, with movement between the two epicenters through B > CA and CA > A. To investigate this specific RJ model with a "hub" FIGURE 7 | Overview of RJ models. Shaded transition rate parameters are excluded; nonshaded transition rate parameters are included (but not necessarily turned on in the posterior models). Green indicates the six transition rate parameters that model the NEC; red indicates the transition rate parameters that revert the NEC. 8 log Bayes Factors (Kass and Raftery 1995: 777). 0 to 2 weak evidence against null hypothesis. 2 to 6 positive evidence against null hypothesis. 6 to 10 strong evidence against null hypothesis. > 10 very strong evidence against null hypothesis. 9 As an ad-hoc test, we constructed a RJ parallel to the RJ_16PAR model, excluding the transition rate parameters that were well attested in the RJ.FULL, such that only the following were left in the model: q12 q14 q15 q16 q23 q24 q25 q26 q35 q45 q51 q52 q53 q54 q62 q64. This model performed much worse than the RJ_16PAR, its log mLh was −159.45; log BF 2(−143.99-−159.45) 30.92, providing decisive evidence for RJ_16PAR over this alternative model with 16 transition rate parameters. In addition, RJ is of critical importance for the good fit of the models reported in Table 1. Without RJ, the alternative RJ_16PAR model has a log mLh was −171.99, so again much worse than the alternative RJ_16PAR model (log BF 25.08).
Frontiers in Communication | www.frontiersin.org September 2021 | Volume 6 | Article 661862 more fully, we constructed RJ_9PAR_BC and RJ_9PAR_C, with the same amount of parameters, but change out of Type-B toward Type-BC (RJ_9PAR_BC) and toward Type-C (RJ_9PAR_C). These models are depicted in Figure 8.
Model fit assessment using log Bayes Factors is given in Table 1, where each model is compared to RJ.FULL. We can identify several "zones" of model fit. The best fitting RJ model is RJ_15PAR, with a log mLh of -142.95. However, the fit of RJ_14PAR is not significantly different from RJ_15PAR (log BF < 2), and the fit of RJ_13PAR, RJ_16PAR, and RJ_17PAR is only marginally worse than RJ_15PAR's. These models score much better than RJ.FULL, suggesting some manual restrictions on transition rate parameters help model fit. The next "zone" of model fit is that of RJ_12PAR through all RJ_9PARs, with log mLh between −147.02 and −149.23 (log BF 4.42). These scores are significantly worse than those of the first "zone," suggesting that RJ prefers a larger number of free parameters to choose from. Last, we can compare the three models with nine parameters, which differ in the "hub" through which change in types BC, C, and CA is directed. Here, RJ_9PAR_CA and RJ_9PAR_BC have similar fit (BF < 2), whereas RJ_9PAR_C performs significantly worse, log BF 4.42 and 2.92, respectively.
As described above, RJ MCMC estimation can set the transition rate(s) to be equal or shared across q-parameters. This happened in RJ.FULL in 989/100 posterior models, showing there is no evidence for multiple rates even in the models including only 9 or 10 transition rate parameters. This is true for all other RJ models in Table 1, with RJ_12PAR being possibly the only exception. One might have expected some evidence for two or more rates as more transition rate parameters were excluded, and the space for model optimization shrank; this is not borne out by the results reported in Table 1.
The inference of the ancestral state for Proto-Indo-European is detailed in Table 2. 10 From RJ.FULL at the top to RJ_9PAR_C at the bottom, we observe a distinct tendency for Type-A to be reconstructed for Proto-Indo-European with increasing certainty. Type-AB, despite its frequency in the data set, is not estimated to be ancestral in the best-fitting models. Throughout the consecutive exclusion of the transition rate parameters, parameters leading away from Type-B are excluded, decreasing the probability that Proto-Indo-European was Type-B. Note that this shows that the differences in ancestral state estimation across models depend directly on the transition rate parameters that are included. The ancestral state estimations of the three models with nine parameters match the "hub" out of which change between Type-BC, C, and CA is directed.

Manual Models
Alongside using RJ MCMC to establish which transitions are relevant, we also tested models that we constructed manually, inspired by Croft's (1991) NEC and by the changes observed in our dataset (see Section 3). The minimum number of parameters that has to be included is six; otherwise, there are states/types that cannot be reached and BayesTraits will not run. Therefore, two minimal models are Croft's (1991) Figure 9, illustrating the variable rates in the two models, shows that neither model makes a lot of sense given what we know about diachronic change in negative existentials. The NEC model suggests more or less comparable rates toward A, AB, B, BC, and CA, hardly any change toward Type-C, and a lot of change toward Type-A. Croft (1991) and Veselinova (2016) have pointed out that type CA is rare in the languages of the world, which is also true of our sample (6 instances of CA out of 106 languages). In the NEC model, the root of the tree, Proto-Indo-European, is estimated to be Type-CA with 0.99 probability, which explains the massive change away from CA. Given the results in Section 4.1 (Proto-Indo-European was likely Type-A), the cross-linguistic rarity of CA, and its unstable nature, this result is probably false. REV.NEC shows even larger rate disparity across parameters, with change toward Type-C, BC, and B being far more common than toward the other parameters. Proto-Indo-European is estimated to be Type-A with 0.96 probability in REV.NEC, which does not explain this disparity in rates.
We further tested a range of models informed by the NEC, the set of changes outside NEC, presented in Section 3, and the prevalence of attested change toward Types A and AB: 1. NEC + extra: parameters included in the NEC plus other attested changes q13 (A > B), q15 (A > C), q16 (A > CA), and q24 (AB > BC) (see (14)). 2. REV.NEC + extra: exact reverse of NEC plus extra parameters from NEC (not reversed). 3. NEC + ALL_X: parameters included in the NEC and four parameters that lead to Type-X, with separate analyses for each type. For instance, model NEC + ALL_A includes the NEC + q21 (AB > A), q31 (B > A), q41 (BC > A), q51 (C > A). 4. REV.NEC + ALL_X: parameters included in REV.NEC and four parameters that lead to Type-X, with separate analyses for each type. 5. PARSIMONY: only parameters that can be observed on the tree when a strict parsimony analysis is conducted. This implies looking at the tree presented in Figure 4 and Figure 5, and observing changes leading to languages we have data on, ignoring uncertainty in the ancestral state estimation.  CA)) and five parameters back into type A (q21 (AB > A), q31 (B > A), q41 (BC > A), q51 (C > A), q61 (CA > A)).
We include in Figure 10 and Table 3 only the best fitting models out of the sets above (see for a full description Supplementary Information S4). For convenience, models NEC and REV.NEC are also included in Figure 10, and Table 3 lists RJ.FULL again. We ran the manual models twice: once using a uniform prior and once using RJ. Using uniform priors, each transition rate parameter is included and takes its own, individual rate (initially sampled from the uniform prior, 0-100). Using RJ, we again allow transition rate parameters to be turned off, and allow for a unified rate of change across included transition rate parameters. We use both uniform priors and RJ because we want to test 1) the manual models informed by our typological and historical analysis directly using the uniform prior, i.e., without further optimization by RJ; and 2) if the fit of the manual models improves by using RJ, most importantly considering whether a unified rate of change is supported. The RJ manual models also provide us with a better comparison to the RJ models discussed in Section 4.1. Table 3 shows that RJ manual models outperformed models with a uniform prior by a positive to a large margin (except for NEC and REV.NEC, as discussed above).
The RJ models presented in Table 1 perform better than the manually constructed models reported on in Table 3, regardless of the latter's prior. Three RJ manual models, REV.NEC + ALL_X, REV.NEC + extra, and PARSIMONY, did not converge or did not run. Most manual models outperform NEC and REV.NEC. The NEC + extra and the PARSIMONY models did not fit better than the NEC model (log BF > 2). This probably has to do with how both emphasize change away from Type-A. Type-A is the most common type attested, and including transition rate parameters toward it, especially q21 (AB < A, such as included in RJ models in Section 4.1, REV.NEC, REV.NEC + extra), improved model performance.
Out of all the models where we add parameters to those involved in the NEC (NEC + extra, NEC + ALL_X, NEC + PARSIMONY), the NEC + ALL_B model performs best (uniform: log mLh −158.49, RJ: log mLh -155.53). However, NEC + ALL_A and NEC + ALL_AB perform equally well ( Supplementary Information S4). Figure 11 illustrates the transition rate parameter settings for the most common posterior model of NEC + ALL_B (attested in 594/1,000 posterior models, with two rates). It shows that the transition rate parameters of the NEC are not all turned on as q23 (AB > B) is turned off. Type-AB becomes an endpoint type, where languages get stuck. Nevertheless, cyclicity still moves from A > B > BC > C > CA > A in this model, with additional parameters leading to type B, out of which only one (BC > B) takes the fast rate.
Out of all models compared in Table 3, the best performing one is ALL_THROUGH_B (log BF > 2 with all other models, an illustration is given in Figure 12). Again, however, Supplementary Information S4 states that ALL_THROUGH_B, ALL_THROUGH_C, and ALL_THROUGH_CA perform equally well. We believe that this can be at least partially explained by the distribution of negative existential types in our sample. There is a very skewed distribution toward Type-A and AB and fewer instances of Type-B, BC, C, and CA (see Figure 2). Hence, it makes (mathematical) sense to have change leading out of an infrequent type to more frequent types, especially Type-A and AB. REV.NEC + ALL_A performs equally well as ALL_THROUGH_B, showing again that models which allow for transitions toward the commonly attested types are preferred. This result may be distinctive for Indo-FIGURE 10 | Overview of the manual models. Shaded transition rate parameters are excluded; non-shaded transition rate parameters are included (but not necessarily turned on). Green indicates the six transition rate parameters that model the NEC; red indicates the transition rate parameters that reverse the NEC. European, with its marked prevalence of Type-A, AB, and B; or it may tentatively suggest that there is something different regarding Type-A, AB, and B, through which a special negative existential marker arises, and Type-BC, C, and CA, through which a negative existential marker replaces in part and takes over standard negation. This hypothesis is further fleshed out in Section 5.
None of the manual models we thought to be most relevant outperform the best-performing RJ models presented in Section 4.1. In the beginning of this section we already had to conclude that the NEC as formulated by Croft (1991), with six changes, does not suffice in a the modeling context. Regardless of the approach, RJ or manual plus RJ, the best fitting models are those that allow for a "cycle within a cycle," i.e., to have several ways to move between parameters, and not the unidirectional way implied by the NEC.

CONCLUSION
The results presented above further underline some of the claims put forth by Veselinova (2013Veselinova ( , 2014Veselinova ( , 2016; see also Verkerk and Shirtz, forthcoming): the NEC does not represent the entire set of historical changes in the domain of negative existence, and other transitions do occur. It is obvious, however, from the results presented here as well as from the work cited above that the six transitions of the NEC do occur both in Indo-European and in other language families (there is indication for a complete, or a nearly complete, cycle occurring across several subfamilies of Indo-Aryan). This means, we believe, that the NEC is neither "false" nor "unhelpful" for understanding historical changes in the negative existential domain. The NEC is simply not the complete story, and this explains, at least in part, the results of phylogenetic modeling presented in Section 4. Unlike what a simple, unidirectional, cyclic model would imply, the domain of negative existence in Indo-European is not easily modeled as a closed subsystem of grammar. Some transitions involve innovations in pre-existing negative existential constructions, but many transitions that were identified here and elsewhere (see also Veselinova 2016: 151ff) involve innovations in constructions outside the domain of negative existence that either lead to innovative (negative) existential constructions or influence the classification of already present negative existential constructions.
The changes from outside the domain of negative existence affect the different negative existential construction types in different ways. The nine parameter RJ models from Section 4.1 suggest cyclical change AB > A > B > AB, with changes involving BC, C, and CA modeled differently. This may be because renewal of standard negation strategies outside of the negative existential domain and the emergence of new existential verbs impacts negative existential types A, AB, and B more directly than types BC, C, or CA. The six possible changes between types A, AB, and B are all attested, suggesting that change between A, AB, and B may in fact be bidirectional. The tipping point seems to be B > BC, marking whether the special negative existential makes its move into standard negation or not, but transitions to BC, C, and CA not predicted by the NEC are attested as well, caused by different types of diachronic changes (A > BC in Macedonian/Bulgarian, AB > BC in Russian/Hawaiian, Veselinova 2014; A > C in Kumzari as described above). This suggestion hence remains tentative at this point, because we have limited information on transitions outside of the NEC, and because types BC and CA (C less so) are rare in Indo-European. The skewed distribution of construction types poses a problem both for analytical work, as we do not yet have enough data to count transitions outside the NEC and categorize them in a sensible way, as well as for the phylogenetic models we constructed in Section 4. This is most clear from the ALL_THROUGH_X models (Supplementary Information S4), but also from the rest of the results: the rarer constructions (BC, C, and CA) are modeled as ancestral, with change toward common types A and AB. Extending the data set is an obvious solution here, both in terms of Indo-European languages and including other large language families (given Veselinova 2014 study of Polynesian, Oceanic/Austronesian seems to be an obvious candidate).
A separate issue for studying the NEC and negative existentials at large is the occurrence of multiple strategies in the same language (see Section 1.3). This is not a very common issue, but frequent enough across the languages of the world (Veselinova 2014(Veselinova , 2016) that we cannot ignore it, as is usually done in typology (see Dryer, 2013 andComrie, 2013 for two different strategies to "do away with" this issue). Further analytical work should be devoted to finding common diachronic pathways in how multiple strategies arise, coexist, and resolve in the negative existential domain.
There is no simple historical scenario that explains the synchronic variation in the Indo-European domain of negative existence. The reason we propose here, following Veselinova (2016), is that functional domains such as the negative existential domain are not always closed ecosystems of constructions, and innovative constructions of different types may enter these domains. Constructions resulting from these different processes may coexist, each deployed in its own functional niche, or replace each other after some period of time. When there are many pathways leading into the domain from "outside," the source of many constructions will not be a construction "inside" the domain. The more pathways leading into the domain, the messier the historical process may seem and the more difficult it is to model.
Our proposal entails, then, that processes whose origin is "outside" the domain of negative existence result in transitions of a different nature from processes whose origin is "within" the domain, and this leads to difficulties in modeling changes with a unidirectional model. Innovations whose origin is "within" the negative existential domain involve a reanalysis (+actualization) or an extension of an older negative existential construction, which may lead to a change in the typological classification of a construction. This innovation type is the one assumed in the NEC and other unidirectional models. The (re-)classification of older constructions is less central for innovative constructions involving reanalysis (+actualization) or extension of some material which is "outside" the negative existential domain. The result, if these novel constructions end up replacing the older constructions, is a set of transitions with the same endpoint but with different starting points. Such transitions, when frequent enough, mean that unidirectional models are unlikely to have adequate explanatory power.
Despite the fact that some processes with an "outside" source have been illustrated above and elsewhere, confirming or rejecting our proposal requires a more direct analysis and systematic collection of such instances. If this interpretation is on the right track, we should be able to identify negative existential constructions whose source is clearly outside the negative existential domain. These constructions may be innovative on leaf level (e.g., Romanian Type-A construction with "to be found") or in some ancestral stage (e.g., innovative locative copulas in Iranian). Testing this proposal, then, would require further analytic work, highlighting the complementary relationship between phylogenetic comparative and analytic methods in historical morphosyntax. Subsequent phylogenetic modeling could be used to test the hypothesis that change "outside" and "inside" of the domain is dependent on the construction type.
As we have already mentioned, while our sample is sizable and non-sparse in a diachronic typological context, it is still not comprehensive enough to test the hypotheses here to the fullest extent. This will have to wait for a more comprehensive sample of Indo-European languages or another big language family. Further, our results should not be imposed onto other families: different language families may involve different tendencies and differ in the common transitions between construction types (e.g., Dunn et al., 2011). Our results suggest that the NEC is not an accurate general typological hypothesis, as it does not fully explain the distribution of negative existential construction types in Indo-European (see again Veselinova (2014) for problems raised for the NEC from a Polynesian perspective). We do not believe Indo-European is in any way special, and suspect the patterns we found here are attested throughout the languages of the world. Further investigation of different families will bear out this hypothesis.
Finally, we proposed that one reason for this may be the frequency of state transitions arising from "outside" the negative existential domain. Testing this hypothesis would involve "traditional," analytic, studies of language change in specific subgroups of Indo-European. This illustrates the relationship we envision between phylogenetic comparative studies in historical morphosyntax and more "traditional," analytic studies: they inform and complement one another. More generally, the results here suggest that the more the expression of a functional domain interacts with other domains, the more likely are changes that depend less on the current typological classification of the domain, and the more difficult it will be to model changes in it by a unidirectional model.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.