The study uses data of the Federal State Statistics Service of the Russian Federation (Rosstat) (https://rosstat.gov.ru/). The sample included data on 82 regions of Russia, including the Tyumen and Arkhangelsk regions, without distinguishing districts. The observation period is limited to 10 months from March to December 2020, that is, from the beginning of the epidemic up to the initiation of mass vaccination in the country. In addition to Rosstat data, the study used information from Yandex’s public data visualization and analysis service, Yandex DataLens², and Weather Archive website³.

The indicator of excess mortality per 100 000 population in the region – excess mortality – was used as a dependent variable in the construction of models. Following the majority of researchers, we have deliberately refrained from using the indicators of morbidity and/or mortality from COVID-19 provided by official statistics.

Already in early stages of the pandemic, experts pointed out obvious advantages of using excess mortality rate as the most objective indicator for comparisons (Leon et al. 2020). Indeed, morbidity statistics have an obvious drawback: the registered number of cases significantly depends upon intensity of testing, which may vary by country/region of one country / time periods. They more often test vulnerable groups, for example, the elderly, which can distort the detected morbidity structure (Danilova 2020). Many people get sick asymptomatically, the share of undetected cases is not known exactly and depends, in turn, upon availability of medical services. In addition, different criteria for determining COVID-19 as the official cause of death can be used (Ivanov 2020). It is also important that the coronavirus infection increases the likelihood of death from other diseases, acting as an indirect cause of death. Excess mortality is an excess of the actual mortality rate over the expected one, therefore it includes both official mortality from COVID-19 and unregistered cases as well as cases that are indirectly accounted for by the pandemic. Estimates of excess mortality over the two years of the pandemic, from January 1, 2020 to December 31, 2021, conducted by specialists in 191 countries, showed that excess mortality exceeds the official COVID-19 mortality rates more than three times (COVID-19 Excess Mortality Collaborators 2022).

There are different ways of evaluating excess mortality. This article used the following approach: first, the indicator of relative (per 100 000 population) mortality in each region was calculated for each month of each year – from 2017 to 2019. Then, for each region, a three-year average relative mortality for each month was calculated. Further, for each region, the average relative mortality in the same month in 2017-2019 was subtracted from the relative mortality for each month (from March to December) in 2020. Thus, we have obtained excess mortality for each region of Russia for every month in 2020:

(1)

M – actual deaths;

EM – excess mortality;

i – region index;

j – month index;

t – year index, t = 2017-2019.

In some papers, only the previous year 2019 is used as a basis for comparison (Druzhinin and Molchanova 2021), however, this method seems inaccurate, because random fluctuations in mortality can be observed in one individual year. Using the average for several years as a basis makes it possible to even out such fluctuations. In their study, Goroshko and Patsala use the average indicator for 5 years as a basis (Goroshko and Patsala 2021), however, in recent years mortality in Russia has been significantly reducing, in particular, due to cardiovascular component, so the average for 5 years can turn out to be higher than the 2020 actual indicators. Even using a three-year average as a basis, we have identified some negative values in monthly indicators of excess mortality by region. Since the linear-logarithmic form of the mortality function will be used later in the study, the variable of excess mortality per 100 000 population was transformed by adding to all its values a constant equal to the sum of the modulus of the highest negative value and figure of one.

Independent variables for modeling were selected in accordance with the groups of morbidity and mortality factors identified above (Table 1). To take into account demographic factors, the following three variables were used: the number of pensioners per 100 000 population, the share of internal migrants calculated as net in-migration per 100 000 population, and the share of population living in the capital of the subject of the Russian Federation and other cities with population over 300 000 people. A pensioner is a citizen of the Russian Federation who has realized the right to receive a pension.

The share of pensioners reflects, first, the age structure of the region, and second, what is important for this study, roughly shows the share of the non-employed elderly. Indicator of the number of pensioners per 100 000 population is available on a quarterly basis.

The number of internal migrants was calculated as the sum of the number of interregional and intraregional migrants, this indicator is available by month. Indicator of the share of population living in the capital of the subject of the Russian Federation and other cities with population over 300 000 people was manually calculated for each region as the sum of the number of all residents living in large cities of the region, attributed to the population of the region. This is an annual indicator.

To assess environmental factors, the following variables were used: average temperature and relative air humidity, and the share of green spaces within the city limits per 1 hectare. Data on average temperatures and humidity were taken for each region by month from WeatherArchive website. In winter months, the average temperature in most regions is below zero, therefore, to meet purposes of further logarithmation, the indicator had to be converted. A constant was added to each observation equal to the sum of the modulus of the highest negative value (-37.3 in December 2020 in the Republic of Sakha (Yakutia)) and figure of one. In order to obtain the share of green spaces within the city limits per 1 hectare, for each region, indicator of the total area of green spaces within the city limits was divided by the total area of urban land within the city limits.

Political factors are presented in the analysis as Yandex self-isolation index variable. Self-isolation of population is one of the major non-medical measures to control spread of the epidemic. In order to measure severity of the applied policy, researchers use various indicators. Some of them collect official data on all measures introduced or canceled and compile integral indexes for individual regions of the country or entire countries. This principle is used, in particular, to construct a well-known OxCGRT Severity Index (Oxford Coronavirus Government Response Tracker). Others use population geolocation indicators, which are usually collected by cameramen and other specialized companies, for example, Google Mobility in the U.S. or Baidu Maps in China (Brodeur et al. 2021). In Russia, to assess mobility of population, Yandex self-isolation index has been developed – an integral indicator that is calculated daily based on data on the use of various Yandex applications and services from the very first days of the epidemic. It compares the level of urban activity on a particular day and a regular day before the epidemic. If the activity level is the same as during rush hours of a regular weekday, it means that the self—isolation index is low, 0 point. If the city is quiet as at night, the index equals to 5 points.

Naturally, the Yandex index shows how the anti-infection measures are actually implemented rather than a formal set of anti-infection measures in effect in the region. In this sense, it can be considered as a proxy variable of the measures applied, adjusted for “compliance” of the population. At the same time, it is precisely this approach to taking into account policy measures that seems to be correct – after all, it is actual mobility that affects morbidity/mortality rather than decrees of the governor or resolutions of Rospotrebnadzor (Federal Service for the Oversight of Consumer Protection and Welfare). And vice versa, integral indices compiled on the basis of a set of measures taken require adjustments to the level of their implementation. The Yandex index has already been used in studies on prevalence of COVID-19 on the basis of the Russian data (Dokhov and Topnikov 2021; Egorov et al. 2021).

Yandex website only presents graphical distribution of the index by day and city, its quantitative values by region are not publicly available. To receive quantitative data, the authors requested the Research Department of Yandex LLC on an individual basis. Since this is a daily index, its monthly average for each month was calculated for each region.

Population per capita monetary income, total average area of residential premises per inhabitant, and unemployment among population aged 15 years and older were used as economic indicators. Monthly per capita monetary incomes were adjusted to constant prices in December 2020 using the basic consumer price index (CPI) for goods and services of the Federal State Statistics Service. Rosstat data on the total average area of residential premises per inhabitant of the region are available by month. Data on the unemployment rate among population aged 15 and older are available on a monthly basis.

Of all factors characterizing the health system, the study uses the number of doctors and nurses for one very important reason. The number of hospital beds, as well as expenditures of the healthcare system during the pandemic significantly increased, precisely as a response to the increased morbidity and mortality. Therefore, inclusion of these indicators in the model as independent variables does not seem well-grounded due to obvious endogeneity. Unlike beds and monetary expenses, it is impossible to increase the number of medical personnel in a short period of time, since medical professions require special and long-term education.

As outliers, we have removed data on Moscow, St. Petersburg and Sevastopol from the sample. The outliers are values of the indicator “the share of population living in the capital of the subject of the Russian Federation and other cities with population over 300 000 people” for Moscow, St. Petersburg and Sevastopol, because the indicator in the subject cities equals to 100%. In addition, the Chukotka District was excluded from the sample due to lack of data on the self-isolation index. In order to avoid duplication of information, data on autonomous districts within regions were not used. After removing the outliers, a balanced panel sample was obtained with a total of 780 observations in 78 regions over 10 months.

Distribution of excess mortality by region is shown on the map (Appendix). For visualization, we have selected October with one of the highest excess mortality rates. The value volatility of the dependent variable is high, excess mortality is unevenly distributed across the country. The leaders include the Krasnodar Krai, Volgograd, Saratov and Belgorod regions and some other regions. If we consider dynamics in the country’s average excess mortality by month (Figure 1), we see that the growth is observed throughout 2020 with August as the only exception.

Dynamics in Yandex self-isolation index shows peak values across the country in April, which is quite expected, since it was April when the lockdown was in effect in most regions with the strictest restrictions on population movement (Figure 2).

Table 2 shows descriptive analysis of data used in the study. There is a significant value variation in variables by month and region: the minimum gap in indicators such as the number of medical specialists, share of pensioners or air humidity is about two times, while the spread of independent variable values – excess mortality – ranges from minus 30 to plus 126. Thus, all selected variables demonstrate a significant volatility and can be included in the regression analysis. The correlation matrix constructed for the selected variables did not show any high pair correlation (above 0.6), which indicates absence of multicollinearity and possibility to use all variables in the model simultaneously.

Based on theoretical models of health economics and empirical studies conducted in other countries or earlier in Russia and with due regard to the pandemic complex nature and its both medical and social specific features, we expect that excess mortality will be associated with five groups of factors, namely: demographic characteristics of the region; peculiar features of the natural environment; restrictive policy in force in the region; and economic factors and state of the health system.

The study used regression models of panel data as the main method. For the ease of result interpretation, a linear-logarithmic form of the model has been selected. Therefore, at the stage of data preparation, we have transformed variables of excess mortality and average temperature so that values of all observations became positive. The mortality function looks as follows (2):

Three regression models have been evaluated: ordinarily least squares method (OLS or pooled regression), and fixed and random effect models (RE and FE). The models were evaluated using the R statistical package, version 4.0.2.

A consistent comparison of the estimated models using the Hausman specification test, Breusch–Pagan test and F-test suggests that a fixed effect model is more preferable. However, in this model we lose a number of important variables for which monthly values are not available. Due to peculiar features of the panel model with fixed effects, impact of all variables, which are annual, “goes” into effects of the regions. Therefore, we further interpret the coefficients of both the FE model and the OLS model, which turned out to be statistically significant.

As Table 3 shows, estimates of both models indicate a statistically significant positive relationship between excess mortality in the region and relative number of migrants. Excess mortality is positively associated with air humidity and negatively with its average temperature, with humidity factor being more important. The estimates of both models also confirm a positive correlation of the dependent variable with indicators of per capita income and unemployment. In addition, the FE model reveals a negative relationship between excess mortality and relative number of pensioners in the region and a positive relationship with the self-isolation index. The OLS model re-confirms a significant and negative dependence of the mortality variable upon relative number of medical specialists.

To test results’ robustness of the linear econometric models, we have used the Random Forest model, a powerful algorithm based on the method of decision trees. It allows to establish nonlinear relationships between the selected variables. The models included the same factors that we used for regression models. At the first step, the sample was divided into a training and a test one in the proportion of 75:25. The number of trees was set to equal to 500. After testing the model on both samples, at the next stage, individual factors were checked for importance. The estimates obtained by the Random Forest method do not have an exact quantitative interpretation, like coefficients in linear regression models, but they make it possible to understand which factors are more important in the model. Figures 3 (a) and 3 (b) present results of the modeling carried out in different ways. In both cases, the horizontal axis shows significance of each factor in the model. The modeling results showed that the most important variables associated with excess mortality are as follows: relative air humidity (humidity); average temperature (temp); self-isolation index (self.isolation.aver); number of internal migrants per 100 000 population (migrants); per capita monetary income (income). Just as in the regression models, the share of urban population, housing provision and share of green spaces in cities turned out to be least significant.

Using the “Random Forest” method, we proved that most of the factors that turned out to be significant in the regression panel model with fixed effects were also significant in the nonlinear model. Thus, the modeling results turned out to be quite sustainable.

In general, the study hypotheses have been confirmed: in each of the five groups of the factors identified during literature review and descriptive data analysis, the modeling has identified variables that were statistically significantly associated with excess mortality rates.

Among demographic factors, a positive relationship between excess mortality and number of migrants has been confirmed, which is consistent with studies conducted both in Russia and other countries (Nakada and Urban 2021; Mikhailova and Valsecchi 2020). As expected, the coronavirus infection, like other similar diseases, spreads faster in those regions where population mobility is higher.

An interesting result is a relatively low excess mortality in regions with higher relative number of pensioners. It would seem that mortality among the elderly should be higher (Singh et al., 2021). However, people of this age group, first, do not work and objectively have fewer social contacts, and second, they are more careful and voluntarily comply with self-isolation requirements, even in the established restrictions are invalid (Kalinin et al. 2020). Therefore, disease incidence in this group is lower compared to younger population reducing average mortality in the region. This result is consistent with earlier studies (Goroshko and Patsala 2021; Zemtsov and Baburin 2020; Pilyasov et al. 2021).

The paper confirms significance of climatic factors: other things being equal, high temperatures reduce excess mortality, while high humidity, on the contrary, increases excess mortality. This conclusion confirms results of the study (Pramanik et al. 2022) regarding COVID-19 prevalence in Russia, although it disagrees with findings of another study (Sabgayda and Zubko 2021) based on data on two cold months of 2020.

The self-isolation index turned out to be significantly and positively associated with excess mortality. This means that in those regions and in those months when the index is higher and mobility is lower, mortality rates are higher as well. Apparently, there is an inverse relationship out here: poor epidemiological situation makes people stay home regardless of formal restrictions, either because of illness and quarantine, or simply because of the fear of getting infected. This fact is consistent with conclusion of the study (Maloney and Taskin 2020). However, such conclusion obviously contradicts the claims about effectiveness of the compulsory large-scale isolation of population, which was used by authorities in many countries, including Russia, in the first months of the pandemic.

Among significant economic factors, per capita income especially noteworthy. Its positive association with excess mortality was confirmed by all models used in the work. It is natural that in richer regions with well-developed industry, transport and wholesale trade, the level of business activity is higher and more people continued to work even during the lockdown, therefore both morbidity and mortality were relatively higher. The obtained conclusion contradicts results of many foreign and one Russian work (Khalatbari-Soltani et al. 2020; Caul 2020; Ettensperger 2021; Dokhov and Topnikov 2021) but does comply with results of a cross-country study (Chaudhry et al. 2020) and studies based on data of the Russian regions (Zemtsov and Baburin 2020; Pilyasov et al. 2021).

Another economic factor – unemployment –turned out to be positively associated with excess mortality as well. One can assume that people who lost their job due to closure of enterprises were made search for vacancies or part-time jobs and could not work on a remote basis. For the first time ever such conclusion was made based on the Russian data and does correspond to findings of studies conducted in other countries (Sun et al. 2021).

Finally, a negative relation between excess mortality in the region and number of medical specialists is another important result of the study. This conclusion is consistent with findings of many foreign studies (Buja et al. 2022; Cifuentes-Faura 2021) and one Russian study (Stepanov 2020). Unlike number of hospital beds, which follows development of the epidemic and therefore positively correlates with both morbidity and mortality, the number of doctors and nurses is exogenous, at least in the short term. Thus, one can argue that regions with high number of medical specialists better coped with morbidity and, subsequently, had relatively low excess mortality, other things being equal.

Findings of this study may be used for developing public health policy in case of repeated waves of the epidemic or emergence of new viruses. Vast extent of the Russian territory and significant variation in climatic and socio-economic characteristics of its regions require differentiated policy measures, rather than its standardization. In particular, priority should be given to epidemiological situation in regions with humid climates and low temperatures, high incomes, intensive migration, and high unemployment – for example, in the framework of the next vaccination program. The COVID-19 pandemic has clearly demonstrated the need for significant investments in medical education to increase the number of medical specialists, which is extremely unevenly distributed across regions nowadays. This approach turns out to be more effective in terms of reducing mortality rather than restrictions on mobility.

Limitations of the study are mainly related to shortcomings of information. Thus, individual statistical indicators are presented in Rosstat databases by region only on an annual or quarterly basis, significantly limiting capabilities of the analysis. Information on anti-infection policy measures implemented by regions is either unified nor systematized. There are no reliable data on the scale of population testing by region and in dynamics, making it impossible to include this important factor into modeling.

Furthermore, in this study, we have intentionally used information for 2020 only, when mass vaccination was yet to be initiated in the country. However, a further analysis, starting from 2021, does require inclusion of this factor. Meanwhile, there are no official data on the number of vaccinated by region and in dynamics yet. It is possible to use findings of sociological surveys of population and employers along with such data accumulation.