Ukrainization and the Effect of Russian Language on the Web: The Google Trends Case Study

ABSTRACT

We consider the question of language diffusion and examine the latest attempt to popularize the Ukrainian language. This effort has been undertaken by the Ukrainian government since 1991. While sociological studies report positive dynamics in Ukrainian language dissemination, in practice, the language adoption progress might be slower than reported. We test this hypothesis using google trends web search data. We apply a Bayesian beta regression model and show that the proposed model is appropriate. The results suggest that the Ukrainian language popularization policy is successful and that the 2022 Russian invasion has considerably intensified the usage of the Ukrainian language.

Introduction

The Russian Federation has enormous military power. According to the global fire power ranking¹, its army is placed second on our planet. Moscow also possesses a powerful cultural weapon that is frequently overlooked – the Russian language (Decker 2021). The Russian language has an enormous unifying power for 193 ethnic groups that live in Russia, with more than 95% speaking Russian (Gladkova 2015). Despite the dissolution of the Soviet Union, which led to the creation of many new independent states, Moscow is more than willing to use both cultural and linguistic ties with a view to promoting its geopolitical interests (Decker 2021; Feklyunina 2016; Forsberg and Smith 2016). The usage of this (soft) power is especially applicable (but is not limited) to a range of former Soviet republics including Ukraine, Belarus, and Kazakhstan. In this study, we examine the case of Ukraine, and in particular the most recent Ukrainization attempt that followed the 2014 Euromaidan revolution. Since 1991, the Ukrainian government has opposed the promotion of Russian cultural and linguistic establishments on Ukrainian territory. These attempts have become more comprehensive since 2014. The extensive Ukrainization policy (see Kiss 2022) was perceived by Donetsk and Luhansk regions as an effort to prohibit the Russian language. This could be an important element that contributed to the development of separatist movements in eastern Ukraine. In addition, this situation was a perfect opportunity for Russian president Vladimir Putin to pursue his geopolitical aims. The Ukrainian case gives us a chance to rigorously examine the effect of de-Russification policies in the post Soviet space; the latter being the main objective of this study.

Several studies suggest that language is an important aspect of society’s self-identification (Arel 2002; Kulyk 2011; Poses and Revilla 2021; Reznik 2018; Ushchyna 2020). An implementation of Ukrainization policy implies increasing the usage of the Ukrainian language, promoting the local culture in various domains, and discouraging both Russian cultural influence and the Russian language. For a comprehensive historical overview of Ukrainian language development, we refer to Flier and Graziosi (n.d.). While such policy can be beneficial for Ukrainian identity and selfhood, it is important to remember that for some regions, especially for those that share a border with the Russian Federation, the corresponding ties to Russia are extremely strong. In addition, many people who live in these regions consider Russian as their native language and have (historically) held pro-Russian views (Bureiko and Moga 2018). Such a regional diversity effect has often been exploited by politicians. As Reetta Toivanen (2007) observes, “Language boundaries, real or imagined, can easily become exploited politically.” Not surprisingly, Ukrainian politicians have also capitalized on linguistic and cultural differences to accumulate electoral support (Kulyk 2011).

Language planning development as part of Ukrainization policy has always been an important concern in Ukraine (Kiss 2022). In 1989, the Ukrainian Soviet Socialist Republic adopted the Law of Languages; this law declared the Ukrainian language as the only official language in the state. It is convenient to divide the 1991–2021 time frame into three periods: the early independence years (1991–2004); the Orange revolution, which took place in 2004; and the post-Euromaidan era, which started in 2014. Under President Leonid Kravchuk (1991–1994), a de-Russification of schools was begun. President Viktor Yushchenko (2005–2010) started to de-Russify the media landscape. However, only mild Ukrainization attempts were made by president Yushchenko’s government. One example of such actions is the requirement that television and radio should have a quota of 75% minimum Ukrainian-language programs, and that there should be an audio-dubbing in Ukrainian for programs that were originally broadcast using other languages. In 2010, when the election was won by Viktor Yanukovych, Ukrainization policy implementation experienced a significant slowdown. A law sponsored by Kivalov and Kolesnichenko (both members of the Verkhovna Rada), called Principles of State Language Policy from 2012, granted the Russian language the status of a regional language. Effectively, the law allowed the usage of minority languages (and Russian in particular) in government institutions such as schools and courts in regions where national minorities exceeded 10% of the population (Elder 2012). The most important event that happened in 2014 on the linguistic front was when parliament repealed (February 23, 2014) the Principles of State Language Policy law. Although President Oleksandr Turchynov (February 23, 2014–June 7, 2014) and President Petro Poroshenko (2014–2019) refused to sign the removal of the law (so the law remained in force until February 2018), the parliament’s decision provided a pretext for Moscow to militarily annex Crimea and promote separatist movements in the east (Reznik 2018).

After the Euromaidan events, Ukrainian society experienced major changes. The loss of Crimea and the war in the eastern Ukraine led to a dramatic increase in national identity and influenced the growing usage of the Ukrainian language. For example, in 2018, the Kivalov–Kolesnichenko language law was declared unconstitutional by the Constitutional Court of Ukraine. According to the Ukrainian Center for Economic and Political Studies (Centr Razumkova 2016), 69%, 27%, and 2% of Ukrainians consider the Ukrainian language, the Russian language, or another language, respectively, as their native language. Nevertheless, certain care should be taken when discussing the “native language” term (Hentschel and Palinska 2022; Zeller 2021). Specifically, it does not necessarily mean that this language is used in practice but instead, it might correspond to nationality, heritage, or the country of residence.

The de-Russification efforts of the Ukrainian government included Ukrainian language quotas for television and radio broadcasting (Ogarkova 2018) and a de-communization law under which the majority of geographical names with reference to Soviet era were changed. To further support de-Russification laws, in February 2017 the Ukrainian government banned the commercial importation of books from Russia. According to Ukrainian sources, the Ukrainization policy is very fruitful. The 2021 sociological service poll of the Ilko Kucheriv Democratic Initiatives Foundation and Razumkov Center (2021), which considered 2,019 respondents aged 18 and older, showed that 78% of Ukrainians identify the Ukrainian language as their native language as opposed to 18 who stated that their native language is Russian. As expected, Russia-oriented media sources like Russia Today (RT) report that the majority of the population in Ukraine is actively using the Russian language and in fact prefers to use Russian. For example, based on polls from the Social Monitoring Center, RT reports that more than 50% of Ukrainian citizens are willing to consume books and media that are delivered in Russian. Moreover, they also report that less than a third of the population supports the usage of the Ukrainian language only. RT thus arrives at the conclusion that the “forceful” Ukrainization of the population, which started in 2014, is not effective overall (for details, please see Latyshev, Onischuk, and Medvedeva 2021). We would like to note that we could not independently verify the Social Monitoring Center poll, since this resource is no longer available.

In this work, we do not intend to make any political claims in favor of Ukrainization or Russification of Ukraine. We do want to note that due to different polls, there is a considerable lack of clarity and uncertainty about the actual usage of the Ukrainian and Russian languages in Ukraine. In addition, some responders might be reluctant about stating their true preferences in both the controlled and the uncontrolled territories. With this in mind, we aim to rigorously investigate the hypothesis that the usage of the Ukrainian language is actually growing by examining an independent and self-sufficient data source. In order to accomplish this, we utilize the google trend data from 2011 to 2021 in order to study the dynamic of change in percentage of Ukrainian language usage. In order to understand the effect of Ukrainization on different regions – namely, regions that are geographically (historically and culturally) closer to Russia or those having such proximity to the “West”² – we apply a Bayesian beta regression model which can take into account the effect of regions (noting that significant parts of two regions were annexed by Russian Federation in 2014). We verify that the model fits the data well and that the Ukrainization policy is successful in all regions (except for Crimea and Sevastopol). However, we also show that in practice, Ukrainian language adoption by the Ukrainian population might be slower than reported.

In addition, we provide several conceptual and methodological contributions. First, we show that the annexed regions of Crimea and Sevastopol should have experienced a severe deterioration in Ukrainian language usage. A possible reason for this decline is the fact that usage of the Ukrainian language is discouraged and that many Ukrainians, including Crimean Tatars, have fled from Crimea since 2014. (For additional details about linguistic conflicts, refer to Müller and Wingender 2021; in particular, see “Characterisation of the Language Situation in the Republic of Crimea from the Perspective of Geolinguistics” by Yuri Dorofeev, in Müller and Wingender 2021, Part III, where the situation in Crimea is discussed.) Regions that are under the control of the Ukrainian government show a statistically significant increase of Ukrainian language usage. Our model also suggests that some regions (Donetsk and Luhansk), in Ukraine will show very slow adoption of the Ukrainian language, which is not very surprising since Donetsk and Luhansk regions have a strong pro-Russian agenda and are partly controlled by separatists. Moreover, with respect to these regions, the question of multiethnicity of Ukraine is raised. The question of Ukraine being a multiethnic country was recently examined by Volodymyr Kulyk (2022b), where the author focuses on the disappearing differentiation between the two largest population groups in Ukraine, Ukrainians and Russians, in Soviet times.

The rest of the paper is organized as follows. In the next section, we formally define the methods used in this paper. We show that one can use the Bayesian beta regression models to fit a language usage data extracted from google trends, and also validate the model and perform efficient prior selection. The results are discussed in the following section. The next section is dedicated to the full-scale Russian invasion in 2022, where we show that the proposed model can be easily adjusted to account for large-scale changes. In the final section we summarize our findings and discuss both the limitations of the proposed method and the possible directions for future research.

Methods

Data Collection

Google and other media companies collect useful data about their customers’ preferences. Here, we take advantage of publicly available Google trends data and, in particular, of the interest by sub-region option which is attainable for specific search terms. Under our setting, a search term is translated into two languages – Ukrainian and Russian – while making sure that these terms are written differently, in order to distinguish between them. Effectively, we now have two search terms that have the same meaning but, are treated as different search labels by Google. Then, the sub-region data provide the percentage of queries recorded for these two terms for a specific region and for a given time frame. Using the Leipzig Corpora Collection (Quasthoff, Goldhahn, and Eckart 2014)³, a list of 50 popular terms was created (for the full list, see Appendix A). The recorded terms are associated with the News, Web, and Wikipedia domains. To ensure a fair comparison, frequent Russian and Ukrainian words were recorded. Finally, for each term, region, and time frame (year), the data from the google trends website was extracted. Eventually, for each region and for years from 2011 to 2021, we calculated the average proportion of searches in the Ukrainian language. Formally, we are working with a quantity:

(1) $proportion\stackrel{def}{\hspace{0.17em}=}\frac{percentageofsearchesinUkrainian}{\begin{array}{c}percentageofsearchesinUkrainian\\ +percentageofsearchesinRussian\end{array}}.$(1)

For the rest of the paper we refer to (1) as the proportion. In order to ensure that the sample of size 50 is indeed representative, we took several random samples (without replacement), of sizes 30 and 40 out of 50. Using these reduced data-sets, we performed a statistical analysis. The obtained results were similar to the ones presented in this manuscript.

Exploratory Analysis

Figure 1 depicts the map of Ukraine divided into administrative districts (regions); the detailed region mapping is given in Table 1.

Table 1. 27 Ukrainian Administrative Divisions; The $Crime{a}^{\ast }$ Peninsula and the City of $Sevastopol{ }^{\ast }$ Were Annexed in 2014 by Russian Federation; The Regions of $Donets{k}^{\mathrm{\Delta }}$ and $Luhans{k}^{\mathrm{\Delta }}$ are Partially Controlled by Separatists

Region id	Region name	Admin. division	Region id	Region name	Admin. division
1	Chernivtsi	West	15	Sumy	Center
2	Ivano-Frankivsk	West	16	Vinnytsia	Center
3	Khmelnytskyi	West	17	Zhytomyr	Center
4	Lviv	West	18	$\mathrm{C}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}{\mathrm{a}}^{\ast }$	South
5	Rivne	West	19	Dnipropetrovsk	South
6	Ternopil	West	20	Kherson	South
7	Volyn	West	21	Mykolaiv	South
8	Zakarpattia	West	22	Odessa	South
9	Cherkasy	Center	23	$\mathrm{S}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}{\mathrm{l}}^{\ast }$	South
10	Chernihiv	Center	24	Zaporizhzhia	South
11	Kyiv (city)	Center	25	$\mathrm{D}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{s}{\mathrm{k}}^{\mathrm{\Delta }}$	East
12	Kyiv (region)	Center	26	Kharkiv	East
13	Kirovohrad	Center	27	$\mathrm{L}\mathrm{u}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{s}{\mathrm{k}}^{\mathrm{\Delta }}$	East
14	Poltava	Center

Figure 1. Ukraine map with regions divided into western, central, southern, and eastern regions. The western part of Ukraine consists of regions that are geographically close to “Western” counties, namely, to Poland, Slovakia, Hungary, and Roumania. The eastern part of Ukraine contains regions that are geographically close to the Russian Federation.

The average proportion of Ukrainian language usage for each region between 2011 and 2021 is depicted in Figure 2. While Figure 2 shows positive dynamics for the Ukrainian language, it is important to rigorously investigate the phenomena using an appropriate regression analysis. With a view to modeling proportions and in order to allow a natural interpretation of the obtained results, we propose to utilize the beta regression model (Ferrari and Cribari-Neto 2004; Figueroa-Zúñiga, Arellano-Valle, and Ferrari 2013), which is discussed below.

Figure 2. The dynamic of the proportion of the Ukrainian language used from 2011 to 2021.

The Proposed Beta Regression Model

A continuous random variable $Y$ is said to have a $Beta$ distribution if its probability density function is given by

$f(y|\alpha ,\beta )=\left\{\begin{array}{cc}\frac{y^{\alpha -1}{(1-y)}^{\beta -1}}{\mathcal{B}\left(\alpha ,\beta \right)}& 0\le y\le 1,\\ 0& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\end{array}\right.$

where α > 0, β > 0, and $\mathcal{B}\left(\alpha ,\beta \right)$ is the beta function (Grimmett and Stirzaker 2001). For $Y\sim Beta\left(\alpha ,\beta \right)$, it holds that $\mathbb{E}Y=\frac{\alpha }{\alpha +\beta }$, and therefore, for convenience, we consider the reparametrisation: $\alpha =\mu \varphi ,\mathrm{a}\mathrm{n}\mathrm{d}\beta =\left(1-\mu \right)\varphi$.

Under this reparametrisation, we write $Y\sim Beta\left(\mu ,\varphi \right)$, and arrive at $\mathbb{E}Y=\mu$. We further assume that the Ukrainian language proportion in region $i\in \left\{1,\dots ,27\right\}$ and year $j\in \left\{2011,\dots ,2021\right\}$ is $y_{ij}$, and that $y_{ij}\sim Beta({\mu }_{ij},\varphi )$. In order to ensure that ${\mu }_{ij}\in \left[0,1\right]$ and with a view to providing a natural interpretation via odds ratios, we utilize the logit link function and define:

$\mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{{\mu }_{ij}}{1-{\mu }_{ij}}\right)=\left({\beta }_0+{\beta }_0^{\left(i\right)}\right)+\left({\beta }_1+{\beta }_1^{\left(i\right)}\right)x_j\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{\eta }_{ij},$

where the effect ${\beta }_0$ is the intercept that characterizes baseline state of proportion, and ${\beta }_1$ is the baseline rate of proportion growth. The regional effects ${\beta }_0^{\left(i\right)}$ and ${\beta }_1^{\left(i\right)}$ are associated with region $i$, and $x_j$ is the covariate, specifically, $x_j$ is a function of a year $j\in \left\{2011,\dots ,2021\right\}$. Under the proposed model, the slopes $\left({\beta }_1+{\beta }_1^{\left(i\right)}\right)$, have an appealing interpretation as the change of log-odds that corresponds to a one unit increase in $x_j$, namely

$\begin{array}{rl}& \mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{{\mu }_{ij+1}}{1-{\mu }_{ij+1}}\right)-\mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{{\mu }_{ij}}{1-{\mu }_{ij}}\right)={\eta }_{ij+1}-{\eta }_{ij}=\\ & =\left({\beta }_0+{\beta }_0^{\left(i\right)}\right)+\left({\beta }_1+{\beta }_1^{\left(i\right)}\right)\left(x_j+1\right)-\left({\beta }_0+{\beta }_0^{\left(i\right)}\right)\\ & +\left({\beta }_1+{\beta }_1^{\left(i\right)}\right)x_j=\left({\beta }_1+{\beta }_1^{\left(i\right)}\right).\end{array}$

For computational efficiency, we center the year covariate and define $x_j=j-2016$, where 2016 is the mean of the {2011,…,2021} set. The available data size is not large, and thus we believe that the latter justifies the usage of the Bayesian approach. In addition, we aim to explore a general machinery for future research, which might include a good prior knowledge about the model parameters. Specifically, we propose to use a Bayesian model which is defined via:

(2) $\begin{array}{rl}& y_{ij}|{\mu }_{ij},\varphi \sim Beta\left({\mu }_{ij},\varphi \right),i\in \left\{1,\dots ,27\right\},j\in \left\{2011,\dots ,2021\right\},\\ & {\mu }_{ij}=\frac{{\mathrm{e}}^{{\eta }_{ij}}}{1+{\mathrm{e}}^{{\eta }_{ij}}},i\in \left\{1,\dots ,27\right\},j\in \left\{2011,\dots ,2021\right\},\\ & \varphi \sim U\left(0,{10}^4\right),{\beta }_0,{\beta }_1,{\beta }_0^{\left(i\right)},{\beta }_1^{\left(i\right)}\sim N\left(0,{\sigma }^2\right),i\in \left\{1,\dots ,27\right\}.\end{array}$(2)

The proposed model can be potentially extended to include change-points (regime switching), with a view to specifying different behaviors of the proportion time series, and by specifying distinct parameters ${\varphi }_{ij}$, instead of the single parameter $\varphi$. Our experiments imply that the model in (2) (with σ = 1 parameter for the ${\beta }_0,{\beta }_1,{\beta }_0^{\left(i\right)},{\beta }_1^{\left(i\right)}$ for $i\in \left\{1,\dots ,27\right\}$ coefficients prior), fits the data well. A more detailed discussion regarding the choice of the σ parameter (prior sensitivity) and the goodness of fit is provided below. We proceed with the computational aspects.

Computation and Validation

Computation

Using No U-Turn Sampling scheme (Carpenter, Andrew Gelman, and Hoffman et al. 2017), we generated three chains with 10,000 iterations per chain, where the first 5,000 iterations were used as a warm-up. The 5,000 remaining samples for each chain were thinned by a factor of 10. Therefore, we had the total of 1,500 samples to perform inference. The Gelman-Rubin diagnostics (Brooks and Gelman 1998; Gelman and Rubin 1992) shows good convergence characteristics with the corresponding statistic values around 1.0 for all model parameters.

Convergence Diagnostics

The MCMC sampler shows good convergence results. Figure 3 shows the convergence of the GelmanRubin statistic (Brooks and Gelman 1998; Gelman and Rubin 1992) of three independent MCMC runs of the No U-Turn sampler and Figure 4 depicts a graphical summary. The first, the second, and the third column of Figure 4 correspond to trace, sample auto correlation function, and density plots, respectively. The first, the second, and the third row of Figure 4, correspond to the first, the second and the third independent MCMC run, respectively. For additional typical convergence results, please see Online Appendix B.

Figure 3. Gelman-Rubin diagnostic for parameter ϕ.

Figure 4. Summary of three Markov Chain Monte Carlo runs for the parameter ϕ.

Prior Sensitivity Analysis for Models with Same Parameter Vector

Note that since our attention is restricted to models with priors $\varphi \sim U\left(0,{10}^4\right)$, and ${\beta }_0,{\beta }_1,{\beta }_0^{\left(i\right)},{\beta }_1^{\left(i\right)}\sim N\left(0,{\sigma }^2\right)$ for $i\in \left\{1,\dots ,27\right\}$, for different values of σ, all competitive models have the same parameter vector $\theta =\left(\varphi ,{\beta }_0,{\beta }_1,{\beta }_0^{\left(1\right)},\dots ,{\beta }_0^{\left(27\right)},{\beta }_1^{\left(1\right)},\dots ,{\beta }_1^{\left(27\right)}\right)$. Let $M_i$ and $M_j$ be two models that have the same parameter vector. Then, the Bayesian Factors (BF) for models $M_i$ and $M_j$, is given by:

$\mathrm{B}{\mathrm{F}}_{ij}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\frac{p(\mathbf{y}|M_i)}{p(\mathbf{y}|M_j)}.$

Under the same parameter vector assumption, it holds that

$\begin{array}{rl}& p\left(\mathbf{y}|M_i\right)=\underset{\mathrm{\Theta }}{\int }p\left(\mathbf{y}|\theta ,M_i\right)p\left(\theta |M_i\right)d\theta \\ & =\underset{\mathrm{\Theta }}{\int }p\left(\mathbf{y}|\theta ,M_i\right)p\left(\theta |M_i\right)\frac{p\left(\theta |\mathbf{y},M_j\right)}{p\left(\theta |\mathbf{y},M_j\right)}d\theta \\ & =\underset{\mathrm{\Theta }}{\int }\frac{p\left(\mathbf{y}|\theta ,M_i\right)p\left(\theta |M_i\right)}{p\left(\theta |\mathbf{y},M_j\right)}p\left(\theta |\mathbf{y},M_j\right)d\theta \\ & =\underset{\mathrm{\Theta }}{\int }\frac{p\left(\mathbf{y}|\theta ,M_i\right)p\left(\theta |M_i\right)}{p\left(\mathbf{y}|\theta ,M_j\right)p\left(\theta |M_j\right)p{\left(\mathbf{y}|M_j\right)}^{-1}}p\left(\theta |\mathbf{y},M_j\right)d\theta \\ & =p\left(\mathbf{y}|M_j\right)\underset{\mathrm{\Theta }}{\int }\frac{p\left(\mathbf{y}|\theta ,M_i\right)p\left(\theta |M_i\right)}{p\left(\mathbf{y}|\theta ,M_j\right)p\left(\theta |M_j\right)}p\left(\theta |\mathbf{y},M_j\right)d\theta .\end{array}$

Therefore,

(3) $\begin{array}{c}\mathrm{B}{\mathrm{F}}_{ij}=\frac{p(\mathbf{y}|M_i)}{p(\mathbf{y}|M_j)}={\int }_{\mathrm{\Theta }}\frac{p(\mathbf{y}|\theta ,M_i)p(\theta |M_i)}{p(\mathbf{y}|\theta ,M_j)p(\theta |M_j)}p(\theta |\mathbf{y},M_j)d\theta .\end{array}$(3)

Suppose further that for any two models, $M_i$ and $M_j$ and for any $\theta$, it holds that

$p(\mathbf{y}|\theta ,M_i)=p(\mathbf{y}|\theta ,M_j).$

Note that this condition corresponds to the beta regression model in the manuscript, since the only difference between two models is the prior parameter $\sigma$. In this case (3), simplifies to

$\mathrm{B}{\mathrm{F}}_{ij}={\int }_{\mathrm{\Theta }}\frac{p(\theta |M_i)}{p(\theta |M_j)}p(\theta |\mathbf{y},M_j)d\theta ={\mathbb{E}}_{\theta |\mathbf{y},M_j}\frac{p(\theta |M_i)}{p(\theta |M_j)}.$

If one has access to samples from the posterior distribution that corresponds to the $M_j$ model, it is possible to compare $M_j$ to any model $M_i$ without even fitting the $M_i$ model. Specifically, it holds that

${\widehat{\mathrm{B}\mathrm{F}}}_{ij}=\frac{1}{N}\sum _{i=1}^N\frac{p({\theta }_i|M_i)}{p({\theta }_i|M_j)},$

where ${\theta }_1,\dots ,{\theta }_N$ are samples from the posterior distribution that corresponds to $M_j$.

For the proposed model in the manuscript, let $M_{{\sigma }_1}$ and $M_{{\sigma }_2}$ be two competitive models. Then, an estimator ${\widehat{\mathrm{B}\mathrm{F}}}_{M_{{\sigma }_1},M_{{\sigma }_2}}$ for $\mathrm{B}{\mathrm{F}}_{M_{{\sigma }_1},M_{{\sigma }_2}}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}p(\mathbf{y}|M_{{\sigma }_1})p(\mathbf{y}|M_{{\sigma }_2}{)}^{-1}$, can be obtained via

${\widehat{\mathrm{B}\mathrm{F}}}_{M_{{\sigma }_1},M_{{\sigma }_2}}=\frac{1}{N}\sum _{i=1}^N\frac{p({\theta }_i|M_{{\sigma }_1})}{p({\theta }_i|M_{{\sigma }_2})},$

where ${\theta }_1,\dots ,{\theta }_N$ are samples from the posterior distribution that corresponds to $M_{{\sigma }_2}$, and $p({\theta }_i|M)$ is the joint prior distribution that corresponds to model $M$ (Chan and Eisenstat 2015). We used posterior samples that correspond to two models, $M_1$ and $M_5$ to produce Figure 5. The left panel of Figure 5, shows the logarithm of Bayes factor $\left({\widehat{\mathrm{B}\mathrm{F}}}_{M_{\sigma },M_5}\right)$ as a function of σ and one can observe that the largest ${\widehat{\mathrm{B}\mathrm{F}}}_{M_{\sigma },M_5}$ is located around σ = 1. The right panel of Figure 5 depicts the logarithm of Bayes factor $\left({\widehat{\mathrm{B}\mathrm{F}}}_{M_{\sigma },M_1}\right)$ as a function of σ. Combining the observations from the left and the right panel with BF interpretation (see Kass and Raftery 1995 for details), we conclude that $M_1$, namely σ = 1, constitutes an appropriate prior.

Figure 5. Logarithms of Bayes factors as a function of σ.

Remark 1 (A hierarchical model) As an alternative, it is possible to consider a hierarchical model:

(4) $\begin{array}{rl}& y_{ij}|{\mu }_{ij},\varphi \sim Beta\left({\mu }_{ij},\varphi \right),i\in \left\{1,\dots ,27\right\},j\in \left\{2011,\dots ,2021\right\},\\ & {\mu }_{ij}=\frac{e^{{\eta }_{ij}}}{1+e^{{\eta }_{ij}}},i\in \left\{1,\dots ,27\right\},j\in \left\{2011,\dots ,2021\right\},\\ & {\beta }_0,{\beta }_1,{\beta }_0^{\left(i\right)},{\beta }_1^{\left(i\right)}\sim N\left(0,{\sigma }^2\right),i\in \left\{1,\dots ,27\right\}.\\ & \varphi \sim U\left(0,{10}^4\right),\\ & \sigma \sim U\left(0,100\right),\end{array}$(4)

The obtained results are similar to the ones reported above; please see Online Appendix C for additional details. For example, the posterior summary for ${\varphi }_{ij}$, ${\beta }_0$, ${\beta }_1$, and σ parameters that are associated with the hierarchical model (7) provided in Table C1. Indeed, the σ parameter is around 1.0 as expected. Figure C1 shows the slopes and the intercept of the hierarchical model.

Goodness of Fit

In order to validate that the proposed model fits the data, we use the extension of the classical ${\chi }^2$ test proposed by Johnson (Johnson 2004). Johnson showed that if one draws parameters samples from the posterior distribution and evaluate the Parson’s goodness-of-fit statistic at these values, then, regardless of the dimension of the parameter vector, the Parson’s goodness-of-fit statistic is asymptotically distributed as a ${\chi }^2$ random variable with K – 1 degrees of freedom.

Let $y_1,\dots ,y_n$ be a scalar-valued, continuous, identically distributed, conditionally independent observations drawn from probability density function $f(y|\theta )$ ($\theta$ is a multidimensional parameter vector). Let $p(\theta |\mathbf{y})$ be the posterior density of $\theta$ based on the data $\mathbf{y}$ and let $\tilde{\theta }$ be a sample from the posterior distribution, namely, from $p(\theta |\mathbf{y})$. The procedure of constructing the Bayesian ${\chi }^2$ test for goodness of fit is as follows.

1. Let $0=a_0<a_1<\cdots <a_{K-1}<a_K=1$ be (user-predefined), quantiles from a uniform distribution. In addition, let $p_j\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{a}_j-a_{j-1}$.

2. Define the vector ${\mathbf{z}}_j(\tilde{\theta })$ for $1\le j\le K$ to be a K-length vector such that its jth element is 1 and all other elements are zero if

(5) $\begin{array}{c}F(y_i|\tilde{\theta })\in \left[a_{j-1},a_j\right],\end{array}$(5)

where F is the cumulative distribution function which corresponds to $f(y|\theta )$.

• (3) Using the definition in (5), let: $m\left(\tilde{\theta }\right)={\sum }_{i=1}^n{z}_i\left(\tilde{\theta }\right)$. Essentially, the jth component of $m(\tilde{\theta })$ (let us call it $m_j(\tilde{\theta })$, is the number of observations that fell into the jth bin. Note that the bins are determined by the quantiles of the inverse distribution function evaluated at $\tilde{\theta }$.

• (4) Finally, we define

(6) $\begin{array}{c}{R}^B(\tilde{\theta })={\sum }_{j=1}^K\frac{(m_j(\tilde{\theta })-n{p}_j){}^2}{n{p}_j},\end{array}$(6)

where $n{p}_j$ is the expected number of points that should lend in bin j.

Theorem 2.1. Under some regularity conditions (Johnson 2004), $R^B$ converges to the ${\chi }^2$ distribution with K – 1 degrees of freedom as $n\to \mathrm{\infty }$.

Under this setting, the null hypothesis is that there is no significant difference between the observed and the expected values.

Practical considerations: It was shown that the number of bins $K=⌈n^{0.4}⌉$ works well in practice. In principle, it is preferred to base the goodness-of-fit statistic on more than a single sampled value from the posterior distribution. This means that in practice we should aim to calculate an average with respect to samples from the posterior distribution, namely, we use

${\hat{R}}^B=\frac{1}{N}\sum _{i=1}^N{R}^B({\tilde{\theta }}_i),$

where ${\tilde{\theta }}_i$ for $i\in \left\{1,\dots ,N\right\}$ are N samples from the posterior distribution.

In order to test the adequacy of the proposed model, we perform the Bayesian ${\chi }^2$ goodness-of-fit test (Johnson 2004). Following Johnson’s recommendation (Johnson 2004), we define the number of bins $K=n^{0.4}$, where n is the sample size. In our case, there are 27 regions with 11 observations for each district and thus K = 10. The corresponding ${\chi }^2$ test statistic estimator ${\hat{R}}^B$ for goodness of fit, was calculated based on 1,500 posterior samples. For the σ = 1 model, the point estimator ${\hat{R}}^B$ is 12.649 and the corresponding 95% confidence interval is (12.734,13305). Since it holds that ${\chi }_{K-1,0.95}^2\approx 16.919$, we conclude that ${\hat{R}}^B<{\chi }_{K-1,0.95}^2$, so this suggests that the proposed model indeed provides an adequate fit to the data.

Results

As mentioned above, the inference is based on 1,500 posterior samples. In this section, all the results are reported with respect to σ = 1 prior parameter. For the posterior distribution summary tables with respect to σ = 5 prior parameters, please refer to Online Appendix D. Tables 2–4, show the full posterior distribution summaries associated with model (2) and with respect to σ = 1 prior parameter.

Table 2. Posterior Distribution Summaries for $\varphi$ and Baseline Effects ${\beta }_0$ and ${\beta }_1$; the Summary is with Respect to $\sigma =1$ Prior Parameter

			Quantiles
Parameter	Mean	Std. Dev.	0.025	0.25	0.5	0.75	0.975
$\varphi$	369.0	34.6	305.7	344.5	368.1	391.9	437.9
${\beta }_0$	−1.218	0.182	−1.569	−1.342	−1.224	−1.097	−0.854
${\beta }_1$	0.066	0.181	−0.289	−0.063	0.069	0.190	0.414

Table 3. Posterior Distribution Summaries for Intercepts; the Summary is with Respect to $\sigma =1$ Prior Parameter

				Quantiles
Region	Parameter	Mean	Std. Dev.	0.025	0.25	0.5	0.75	0.975
Chernivtsi	${\beta }_0+{\beta }_0^{\left(1\right)}$	−0.265	0.033	−0.328	−0.287	−0.265	−0.243	−0.2
Ivano-Frankivsk	${\beta }_0+{\beta }_0^{\left(2\right)}$	0.592	0.033	0.528	0.571	0.591	0.614	0.66
Khmelnytskyi	${\beta }_0+{\beta }_0^{\left(3\right)}$	−0.246	0.033	−0.31	−0.268	−0.245	−0.225	−0.183
Lviv	${\beta }_0+{\beta }_0^{\left(4\right)}$	0.324	0.032	0.263	0.302	0.323	0.345	0.387
Rivne	${\beta }_0+{\beta }_0^{\left(5\right)}$	0.166	0.032	0.103	0.146	0.167	0.188	0.23
Ternopil	${\beta }_0+{\beta }_0^{\left(6\right)}$	0.608	0.034	0.543	0.585	0.609	0.631	0.674
Volyn	${\beta }_0+{\beta }_0^{\left(7\right)}$	0.248	0.031	0.187	0.228	0.248	0.269	0.309
Zakarpattia	${\beta }_0+{\beta }_0^{\left(8\right)}$	−0.259	0.033	−0.322	−0.28	−0.258	−0.236	−0.196
Cherkasy	${\beta }_0+{\beta }_0^{\left(9\right)}$	−0.783	0.034	−0.85	−0.806	−0.783	−0.76	−0.716
Chernihiv	${\beta }_0+{\beta }_0^{\left(10\right)}$	−1.239	0.037	−1.31	−1.265	−1.239	−1.214	−1.168
Kyiv (city)	${\beta }_0+{\beta }_0^{\left(11\right)}$	−1.039	0.036	−1.107	−1.064	−1.039	−1.016	−0.972
Kyiv (region)	${\beta }_0+{\beta }_0^{\left(12\right)}$	−0.987	0.036	−1.056	−1.01	−0.987	−0.963	−0.917
Kirovohrad	${\beta }_0+{\beta }_0^{\left(13\right)}$	−1.246	0.039	−1.324	−1.271	−1.245	−1.221	−1.168
Poltava	${\beta }_0+{\beta }_0^{\left(14\right)}$	−1.15	0.037	−1.225	−1.174	−1.15	−1.126	−1.08
Sumy	${\beta }_0+{\beta }_0^{\left(15\right)}$	−1.336	0.039	−1.413	−1.361	−1.336	−1.309	−1.259
Vinnytsia	${\beta }_0+{\beta }_0^{\left(16\right)}$	−0.556	0.033	−0.618	−0.579	−0.556	−0.532	−0.492
Zhytomyr	${\beta }_0+{\beta }_0^{\left(17\right)}$	−0.7	0.034	−0.766	−0.722	−0.701	−0.679	−0.629
Crimea	${\beta }_0+{\beta }_0^{\left(18\right)}$	−4.986	0.18	−5.363	−5.102	−4.976	−4.865	−4.653
Dnipropetrovsk	${\beta }_0+{\beta }_0^{\left(19\right)}$	−1.751	0.044	−1.837	−1.78	−1.75	−1.72	−1.668
Kherson	${\beta }_0+{\beta }_0^{\left(20\right)}$	−1.627	0.043	−1.709	−1.655	−1.627	−1.599	−1.546
Mykolaiv	${\beta }_0+{\beta }_0^{\left(21\right)}$	−1.595	0.043	−1.68	−1.624	−1.596	−1.567	−1.511
Odessa	${\beta }_0+{\beta }_0^{\left(22\right)}$	−2.007	0.049	−2.104	−2.039	−2.005	−1.973	−1.912
Sevastopol	${\beta }_0+{\beta }_0^{\left(23\right)}$	−5.162	0.179	−5.519	−5.275	−5.16	−5.036	−4.825
Zaporizhzhia	${\beta }_0+{\beta }_0^{\left(24\right)}$	−1.861	0.045	−1.951	−1.89	−1.86	−1.831	−1.774
Donetsk	${\beta }_0+{\beta }_0^{\left(25\right)}$	−2.599	0.06	−2.719	−2.639	−2.596	−2.559	−2.489
Kharkiv	${\beta }_0+{\beta }_0^{\left(26\right)}$	−1.976	0.048	−2.069	−2.007	−1.974	−1.943	−1.887
Luhansk	${\beta }_0+{\beta }_0^{\left(27\right)}$	−2.654	0.065	−2.78	−2.695	−2.654	−2.61	−2.53

Table 4. Posterior Distribution Summaries for Slopes; the Summary is with Respect to $\sigma =1$ Prior Parameter

				Quantiles
Region	Parameter	Mean	Std. Dev.	0.025	0.25	0.5	0.75	0.975
Chernivtsi	${\beta }_1+{\beta }_1^{\left(1\right)}$	0.108	0.011	0.087	0.101	0.108	0.115	0.128
Ivano-Frankivsk	${\beta }_1+{\beta }_1^{\left(2\right)}$	0.11	0.011	0.089	0.103	0.11	0.117	0.132
Khmelnytskyi	${\beta }_1+{\beta }_1^{\left(3\right)}$	0.138	0.011	0.117	0.13	0.137	0.145	0.159
Lviv	${\beta }_1+{\beta }_1^{\left(4\right)}$	0.124	0.011	0.103	0.118	0.124	0.131	0.145
Rivne	${\beta }_1+{\beta }_1^{\left(5\right)}$	0.124	0.01	0.104	0.117	0.124	0.131	0.144
Ternopil	${\beta }_1+{\beta }_1^{\left(6\right)}$	0.131	0.011	0.11	0.124	0.131	0.139	0.154
Volyn	${\beta }_1+{\beta }_1^{\left(7\right)}$	0.137	0.011	0.117	0.13	0.137	0.145	0.159
Zakarpattia	${\beta }_1+{\beta }_1^{\left(8\right)}$	0.146	0.01	0.125	0.139	0.146	0.153	0.165
Cherkasy	${\beta }_1+{\beta }_1^{\left(9\right)}$	0.115	0.011	0.094	0.107	0.114	0.122	0.136
Chernihiv	${\beta }_1+{\beta }_1^{\left(10\right)}$	0.085	0.013	0.06	0.076	0.085	0.094	0.11
Kyiv (city)	${\beta }_1+{\beta }_1^{\left(11\right)}$	0.033	0.011	0.012	0.026	0.033	0.041	0.055
Kyiv (region)	${\beta }_1+{\beta }_1^{\left(12\right)}$	0.07	0.011	0.048	0.063	0.07	0.077	0.094
Kirovohrad	${\beta }_1+{\beta }_1^{\left(13\right)}$	0.105	0.013	0.08	0.097	0.105	0.114	0.131
Poltava	${\beta }_1+{\beta }_1^{\left(14\right)}$	0.095	0.012	0.072	0.087	0.095	0.103	0.119
Sumy	${\beta }_1+{\beta }_1^{\left(15\right)}$	0.089	0.013	0.065	0.08	0.088	0.097	0.114
Vinnytsia	${\beta }_1+{\beta }_1^{\left(16\right)}$	0.139	0.011	0.117	0.131	0.138	0.147	0.162
Zhytomyr	${\beta }_1+{\beta }_1^{\left(17\right)}$	0.117	0.011	0.096	0.11	0.117	0.124	0.138
Crimea	${\beta }_1+{\beta }_1^{\left(18\right)}$	−0.287	0.054	−0.396	−0.321	−0.286	−0.249	−0.186
Dnipropetrovsk	${\beta }_1+{\beta }_1^{\left(19\right)}$	0.042	0.015	0.011	0.032	0.041	0.052	0.073
Kherson	${\beta }_1+{\beta }_1^{\left(20\right)}$	0.082	0.014	0.057	0.072	0.082	0.091	0.11
Mykolaiv	${\beta }_1+{\beta }_1^{\left(21\right)}$	0.077	0.013	0.051	0.068	0.078	0.086	0.104
Odessa	${\beta }_1+{\beta }_1^{\left(22\right)}$	0.07	0.016	0.039	0.059	0.069	0.081	0.1
Sevastopol	${\beta }_1+{\beta }_1^{\left(23\right)}$	−0.254	0.054	−0.361	−0.291	−0.254	−0.218	−0.146
Zaporizhzhia	${\beta }_1+{\beta }_1^{\left(24\right)}$	0.052	0.015	0.022	0.042	0.051	0.062	0.081
Donetsk	${\beta }_1+{\beta }_1^{\left(25\right)}$	0.013	0.022	−0.028	−0.002	0.013	0.027	0.056
Kharkiv	${\beta }_1+{\beta }_1^{\left(26\right)}$	0.045	0.016	0.015	0.034	0.045	0.056	0.076
Luhansk	${\beta }_1+{\beta }_1^{\left(27\right)}$	0.046	0.022	0.003	0.031	0.045	0.06	0.088

Figure 6 depicts interval estimates from posterior draws associated with model (2).

Figure 6. Left panel shows the regional intercepts and right panel shows the regional slopes.

The left plot shows the intercepts. The intercepts of the western part of Ukraine are generally larger as compared to the eastern part, thus indicating that the usage of the Ukrainian language was always stronger in the western part as compared to other regions. On the other hand, the intercepts of Crimea, Sevastopol, and eastern regions such as Luhansk and Donetsk that have geographic (historical and cultural) proximity to the Russian Federation indicate a lower usage of the Ukrainian language as compared to the Russian language.

While the left panel of Figure 6 aligns well with our intuition, it is even more instructive to consider the right panel of Figure 6. The latter corresponds to the regression slopes and thus shows the dynamics of the proportion, namely, the growth of the proportion of Ukrainian language usage. By observing the slopes, we can see that all the government-controlled territories have positive slopes (the results are statistically significant; please see Table 4 for details). This indicates that the Ukrainization policy is indeed quite successful. Of course, we can also see that the slopes in Crimea and Sevastopol (two regions that were annexed by Russia in 2014) have negative slopes (this result is also statistically significant).

The mean and the corresponding 95% credible intervals (CI) for slopes in Crimea are –0.287 and (–0.396,–0.186), respectively. The situation in Sevastopol is similar: the mean and the corresponding 95% credible intervals (CI) for slopes are –0.254 and (–0.361,–0.146), respectively. Namely, we see a strong Russification of the annexed regions. On the other hand, there are regions in Ukraine that show a rapid growth of Ukrainian language usage. Two such regions are Zakarpattia and Khmelnytskyi, that introduce the slopes of 0.146, (0.125,0.165) and 0.138, (0.117,0.159), respectively. The results suggest that the effect of the Russification of Crimea is stronger than the corresponding Ukrainization effects in other Ukrainian regions.

The regional effect slopes in the Ukrainian territory indicate that the western regions in Ukraine show the biggest increase in the proportion, which is not surprising. We acknowledge that Donetsk and Luhansk are partially controlled by separatists, so negative slopes are expected there. Surprisingly enough, the Luhansk region shows slow growth, with mean slope of 0.046 and 95% CI of (0.003,0.088). However, for the Donetsk region, we observe the mean slope of 0.013 and the corresponding 95% CI of (–0.028,0.056). The latter might indicate that the situation in the Donetsk region is more radical as compared to the Luhansk region.

Overall, the obtained regression slopes indicate that the Ukrainization policy is working, although the progress might be slower than reported by Ukrainian sources. Using the posterior samples and the 2021 population estimates⁴, we were able to approximate the overall population percentage that perform Google searches in Ukrainian language. Our data show that about 35% (the 95% credible interval is (32.56%, 37.48%), see Table 5), perform their search in Ukrainian.

Table 5. Posterior Summary of the Number of People Who Perform Google Search Using Ukrainian Language. The Mean UA Search Column Corresponds to the Estimator of the Number of People that Search Using the Ukrainian Language

Region	Population estimate	Mean UA search	95 Cred Interval
Chernivtsi	8.92 × 10^5	5.07 × 10^5	(4.81 × 10^5, 5.33 × 10^5)
Ivano-Frankivsk	1.35 × 10^6	1.03 × 10^6	(9.95 × 10^5, 1.06 × 10^6)
Khmelnytskyi	1.23 × 10^6	7.50 × 10^5	(7.14 × 10^5, 7.84 × 10^5)
Lviv	2.48 × 10^6	1.79 × 10^6	(1.72 × 10^6, 1.85 × 10^5)
Rivne	1.14 × 10^6	7.85 × 10^5	(7.57 × 10^5, 8.12 × 10^5)
Ternopil	1.02 × 10^6	7.98 × 10^5	(7.74 × 10^5, 8.23 × 10^5)
Volyn	1.02 × 10^6	7.35 × 10^5	(7.09 × 10^5, 7.60 × 10^5)
Zakarpattia	1.25 × 10^6	7.67 × 10^5	(7.32 × 10^5, 8.02 × 10^5)
Cherkasy	1.16 × 10^6	5.22 × 10^5	(4.88 × 10^5, 5.57 × 10^5)
Chernihiv	9.63 × 10^5	2.96 × 10^5	(2.67 × 10^5, 3.25 × 10^5)
Kyiv (city)	2.95 × 10^6	8.72 × 10^5	(7.94 × 10^5, 9.51 × 10^5)
Kyiv (region)	1.80 × 10^6	6.22 × 10^5	(5.73 × 10^5, 6.75 × 10^5)
Kirovohrad	9.08 × 10^5	2.97 × 10^5	(2.70 × 10^5, 3.23 × 10^5)
Poltava	1.36 × 10^6	4.58 × 10^5	(4021 × 10^5, 4.99 × 10^5)
Sumy	1.36 × 10^6	4.58 × 10^5	(2.75 × 10^5, 3.33 × 10^5)
Vinnytsia	1.51 × 10^6	810 × 10^5	(7.65 × 10^5, 8.54 × 10^5)
Zhytomyr	1.18 × 10^6	5.57 × 10^5	(5.24 × 10^5, 5.93 × 10^5)
Crimea
Dnipropetrovsk	3.11 × 10^6	5.49 × 10^5	(4.76 × 10^5, 6.30 × 10^5)
Kherson	1.01 × 10^6	2.30 × 10^5	(2.04 × 10^5, 2.57 × 10^5)
Mykolaiv	1.10 × 10^6	2.52 × 10^5	(2.24 × 10^5, 2.80 × 10^5)
Odessa	2.36 × 10^6	3.77 × 10^5	(3.27 × 10^5, 4.31 × 10^5)
Sevastopol
Zaporizhzhia	1.65 × 10^6	2.76 × 10^5	(2.39 × 10^5, 3.14 × 10^5)
Donetsk	4.07 × 10^6	3.01 × 10^5	(2.36 × 10^5, 3.74 × 10^5)
Kharkiv	2.61 × 10^6	3.86 × 10^5	(3.31 × 10^5, 4.43 × 10^5)
Luhansk	2.11 × 10^6	1.72 × 10^5	(1.13 × 10^5, 2.11 × 10^5)
Summary	4.13 × 10^7	1.44 × 10^7	(1.34 × 10^5, 1.55 × 10^5)
Percentage	100.00%	34.98%	(32.56%, 37.48%)

The posterior summary of the number of people who perform search in Ukrainian language for each region is given in Table 5. From Table 5, we arrive at the conclusion that the percentage of the Ukrainian population that performs Ukrainian language searches is about 35%; this might not align well with the reported 78% who reported that they consider Ukrainian as their native language (Tables D1–D3).

Nevertheless, our model predicts that the number of regions that will exceed the 50% usage threshold of the Ukrainian language grows. Specifically, in 2021, only 9 regions (out of 27), exceeded the 50% threshold. However, in 2026 and 2031, we predict that 11 and 15 regions will exceed this threshold. The latter is an indication of the overall success of the Ukrainization policy, at least with respect to Ukrainian language usage.

The full-scale Russian invasion in 2022 requires an additional validation of the proposed model. In the next section we examine the model and the consequences of recent events on Ukrainian language popularization. We show that the proposed model is still valid subject to the introduction of a change point that occurs in 2022.

The 2022 Russian Invasion

First, we examine how well the current 2011–2021 model fits the data when including the results from 2022. In this case, there are again 27 regions with 12 observations for each district and thus K = 11. The corresponding ${\chi }^2$ test statistic estimator ${\hat{R}}^B$ for goodness of fit, was calculated based on 1,500 posterior samples. For the σ = 1 model, the point estimator ${\hat{R}}^B$ is $19.921$ and the corresponding 95% confidence interval is $\left(19.570,20.271\right)$. Since it holds that ${\chi }_{K-1,0.95}^2\approx 18.307$, we conclude that ${\hat{R}}^B>{\chi }_{K-1,0.95}^2$, so this suggests that the proposed model is not adequate for the 2022 data. While the ${\chi }^2$ test statistic is important from the mathematical point of view, Figure 7, which depicts the 2022 data and the 2011–2021 model prediction intervals for every region is quite instructive. It is interesting to note that the 2011–2021 model is adequate for the western regions. Nevertheless, the majority of data points for the central, southern, and eastern regions are above the prediction intervals. This might be due to the fact that the western regions were less exposed to military actions.

Figure 7. The 2022 data and the 2011–2021 model 95% prediction intervals for all regions.

The above findings are not very surprising and are basically supported by recent studies of Rating 2022 and Kulyk 2022a. As noted by Kulyk 2022a, many Ukrainians tend to blame the Russian population for the war and the associated crimes of the occupying forces. Specifically, according to Kulyk’s study, Russian is now considered to be the language of the enemy and thus, many Ukrainians refuse to utilize Russian and instead try to use Ukrainian, as it is now considered to be the language of the resistance. Our findings support this claim.

The obtained results and in particular, a careful observation of Figure 7, indicates that the 2011–2021 model should be adjusted. We propose to extend the original model by introducing a change point (Rizzo 2019, Chapter 11). The extension of the (2) is as follows.

(7) $\begin{array}{rl}& y_{ij}\sim \left\{\begin{array}{c}Beta({\mu }_{ij},{\varphi }_1)i\in \left\{1,\dots ,27\right\},j\in \left\{2011,\dots ,2021\right\}\\ Beta({\mu }_{ij},{\varphi }_2)i\in \left\{1,\dots ,27\right\},j=2022\end{array}\right.\\ & {\mu }_{ij}=\frac{{\mathrm{e}}^{{\eta }_{ij}}}{1+{\mathrm{e}}^{{\eta }_{ij}}},i\in \left\{1,\dots ,27\right\},j\in \left\{2011,\dots ,2021,2022\right\},\\ & {{\mu }^{\prime }}_{ij}=\frac{{\mathrm{e}}^{{\eta }_{ij}}}{1+{\mathrm{e}}^{{\eta }_{ij}}},i\in \left\{1,\dots ,27\right\},j=2022,\\ & {\varphi }_1,{\varphi }_2\sim U\left(0,{10}^4\right),\\ & {\beta }_{0,}{\beta }_{1,}{\beta }_0^{\left(i\right)},{\beta }_1^{\left(i\right)},\beta 0\left(i\right),\beta 1\left(i\right)\sim N\left(0,{\sigma }^2\right),i\in \left\{1,\dots ,27\right\}.\end{array}$(7)

Here, similarly to the original model, we define:

$\begin{array}{rl}& \mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{{\mu }_{ij}}{1-{\mu }_{ij}}\right)={\eta }_{ij}={\beta }_0+{\beta }_0^{\left(i\right)}+\left({\beta }_1+{\beta }_1^{\left(i\right)}\right)x_j\mathrm{a}\mathrm{n}\mathrm{d}\\ & \mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{{{\mu }^{\prime }}_{ij}}{1-{\mu }_{ij}}\right)={{\eta }^{\prime }}_{ij}=={\beta }_0+{{\beta }_0}^{\prime \left(i\right)}+{\beta }_1+{{\beta }_1}^{\prime \left(i\right)}{x}_j.\end{array}$

Now, for the new 2011–2022 (σ = 1) model, the point estimator ${\hat{R}}^B$ is $15.348$ and the corresponding 95% confidence interval is $\left(14.180,16.254\right)$. It holds that ${\chi }_{K-1,0.95}^2\approx 18.307$, and we conclude that ${\hat{R}}^B<{\chi }_{K-1,0.95}^2$, so this suggests that the proposed model fits the data well. Figure 8 shows a comparison of slopes for the original 2011–2021 model and the new 2011–2022 model.

Figure 8. The 2011–2021 vs 2011–2022 models slope comparison.

Figure 8 is very instructive. The majority of slopes in the 2011–2022 model are much higher when the corresponding slopes in the 2011–2021 model. This further indicates that the Russian invasion contributes to the development and the acceptance of the Ukrainian language and supports the findings of Rating 2022 and Kulyk 2022a. In Crimea and Sevastopol, the change in the slope is actually negative, namely, from –0.2871 to –0.7409 in Crimea and from –0.2537 to –0.9045 in Sevastopol.

Discussion and Conclusion

In this paper, we have demonstrated that google trend data combined with spatial information can provide important insights into language dissemination trends. The proposed beta regression model fits the data well and is able to explain spatial variations. It is important to note that the data are open and that no costly experiments are required. The proposed model combined with google trends data can potentially serve as a verification mechanism to language poll experiments.

Under our model, the prediction of the proportion for the forthcoming years is straightforward. Using the posterior samples associated with the model, the prediction and confidence intervals of the proportion for region i and year j ≥ 2023 can be derived from (7). From Figure 9, which shows the prediction of the Ukrainian language usage proportion in the annexed Crimea and Sevastopol, we conclude that the situation is quite distressing for the Ukrainian language.

Figure 9. Crimea and Sevastopol; proportion of Ukrainian language usage prediction until 2040. The graph shows the 0.025, 0.5, and the 0.975 quantiles of the proportion.

It is also very important to consider the situation in partly controlled territories of Donetsk and Luhansk. The prediction of the proportion for Donetsk and Luhansk is depicted in Figure 10. Our analysis indicates that the situation in both regions looks dreadful for the Ukrainian language. The prediction of the Ukrainian language proportion usage for all regions is depicted in Figures 11 and 12.

Figure 10. Donetsk and Luhansk; proportion of Ukrainian language usage prediction until 2040. The graph shows the 0.025, 0.5, and the 0.975 quantiles of the proportion.

Figure 11. Prediction of Ukrainian language usage proportion until year 2040 (part 1).

Figure 12. Prediction of Ukrainian language usage proportion until year 2040 (part 2).

There are several limitations of this study that should be discussed. First, due to the geopolitical situation in Donetsk and Luhansk regions, we cannot distinguish between the territories that are under control of the Ukrainian government or under control of separatists. That is, the data are combined for these two regions and therefore we cannot observe the corresponding effect directly. It is also important to consider the effect of the Russian invasion of Ukraine that started on February 24, 2022. Specifically, one should keep in mind that many people died or fled these regions and a considerable number of cities and villages were devastated. Moreover, if we consider the annexed regions of Crimea and the city of Sevastopol, there might be an additional effect that is related to service availability. Namely, an individual located in Crimea and seeking, say, government assistance, will need to use the Russian language. However, this logic also applies to the territories that are under the control of the Ukrainian government. An additional limitation of this study is web content availability. Specifically, there exists more content in the Russian language. Furthermore, by using the open google trend data, we cannot distinguish between individuals with respect to say age, education, and so forth. However, this work can be potentially extended by designing appropriate experiments with human subjects; the statistical machinery would remain almost identical. Finally, it is important to note that many Ukrainian citizens use Surzhyk (Hentschel and Palinska 2022; Hentschel and Taranenko 2021), a mixed language that contains both Russian and Ukrainian words. The proposed method can not distinguish between pure Ukrainian, pure Russian, and Surzhyk speakers. However, the method can provide evidence regarding the proportion of Ukrainian and Russian words used.

Despite the above limitations, this work demonstrates the value of the spatial google trend language data availability. Moreover, it lays a foundation for various extensions and future work. For example, starting from February 2022, additional queries and keywords might become popular. It will be of interest to develop a model that both considers the available 2011–2022 data and, takes into account the new set of war-related queries. One possible direction is to consider the Beta rectangular distribution link function as suggested by Bayes et al. (Bayes, Bazán, and García 2012), since it can provide a more robust modeling of proportions with respect to outliers. Applying such a model will be increasingly important as time passes, and additional data from google trends becomes available. To conclude, we conjecture that the application of the proposed model is of great value since it can help to examine the effectiveness of government policies with respect to language dissemination.

Acknowledgments

We are thoroughly grateful to the editor and to anonymous reviewers for their valuable and constructive remarks and suggestions.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Supplemental data

Supplemental data for this article can be accessed on the publisher’s website at https://doi.org/10.1080/10758216.2023.2224568

Notes

1. https://www.globalfirepower.com

2. Poland, Slovakia, Hungary, and Roumania

3. https://wortschatz.uni-leipzig.de/en

4. http://database.ukrcensus.gov.ua/PXWEB2007/eng/news/op_popul_e.asp

Appendix Appendix A.

Google trends search terms