The first group of factors for the use of ART are factors that influence predisposition to infertility. Based survey data analysis some studies reveal that obesity, low education, higher age of the woman and higher age at marriage are associated with higher risks of infertility (Righarts et al. 2015; Sarac, Koc 2017). Environmental factors of infertility have been identified in medical studies: for example, air pollution affects female infertility (Carre et al. 2017; Conforti et al. 2018), and water and soil pollution affect male reproductive health (Di Nisio, Foresta 2019).

The second group of factors for the use of ART is technology awareness factors: the impact of information on the effectiveness of the procedure is discussed in the reports of ART clinics in the United States (Wu 2019), and in Russia the impact of youth awareness of the existence and principles of the technology is estimated within a pilot study in Saransk (Dadaeva, Baranova 2019).

The third group of factors for the use of ART are accessibility factors, or socio-economic factors: the level of development of the health care system, the work of public funds (Adamson 2009).

The fourth group of factors in the use of ART is cultural factors: religious composition of population, traditions in the division of gender roles, level of gender discrimination in the society (Bennett 2016).

Several publications highlight the impact of reproductive technology policy, its effectiveness and benefits to the State. Many existing studies are based on the US data, as measures to regulate the application of the technology differ across the states and were introduced at different times. In the US, health insurance mandates – legislatively enshrined requirements for the employer or individual to acquire health insurance – are arranged differently. This creates an opportunity to use the difference-in-difference method and compare the effectiveness of various regulations. In the work of Marianne P. Bitler and Lucie Schmidt (Bitler, Schmidt 2012), the comparative effectiveness of different formats of health insurance for IVF use is evaluated on the basis of survey data. Basing on the logit model, the authors find that the introduction of insurance mandates did not affect inequalities in access by race and socioeconomic status.

An important question in research on the US is whether the state policy affects treatment programme choices. The choice is compared between a more “aggressive” but more effective treatment option, where more than two embryos are transplanted during an IVF procedure, and a safer one when only one embryo is transplanted. The impact of the health insurance options is analyzed using the least squares method (OLS) with two-way fixed effects at clinic and market levels (Hamilton, McManus 2011). Insurance mandates increase the availability of the technology and can reduce the number of multiple births, but their effectiveness varies widely depending on the type of insurance. In the study by Kasey S. Buckles (Buckles 2013), a similar result is obtained using the difference-in-differences method. Analysis of the effectiveness of insurance programmes for the use of ART was also carried out on the basis of individual patient data with the help of pooled regression (Baker et al. 2007). A cohort analysis of the data of the ART register in the United States has also been applied to the subject of the study (Crawford et al. 2016). This study identified no change in the effectiveness of the use of ART, but revealed the impact of insurance mandates on the general prevalence of ART.

A number of studies consider ART as a factor influencing fertility rates: number of births, eventual cohort fertility, total fertility rate.

Based on the simulation model of population reproduction, built using the Monte Carlo method on historical data for France, Henri Leridon (Leridon 2004) shows that the development of the technology existing at the time of the study made it possible to make up for only 50% of the loss associated with postponement of births in the age group 30 to 35 years old. In his works Lone Schmidt (Schmidt 2005, 2007) on the basis of micro and macro data using difference-in-differences and difference-in-difference-in-differences methods shows that insurance mandates have increased the number of births of first children born by women over 35 years of age. Local studies have been conducted in Russia, for example, a study at the municipal level (Vasilyeva, Peregontseva 2016), in which the authors consider the impact of the development of ART on fertility in the Bryansk oblast based on the analysis of patient cards. The authors concluded that the opening of the ART office had a positive impact on fertility rates in the region.

Studies show the presence of a significant positive relationship between the number of ART procedures and the total fertility rate (TFR) in Europe (Burcin et al. 2014). Another result was achieved when analyzing the impact of the growth of government support for infertility treatment on TFR and the age of the mother at birth of the first child in the United States (Machado et al. 2015). The probit model shows that insurance mandates for infertility treatment in the 1980-1990s did not affect the total birth rate of cohorts, rather, they led to a shift in the births calendar and postponement of motherhood.

Basing on the literature review results, the author of this study systemized factors affecting the use of ART and composed a list of relevant control variables. Figure 1 shows a scheme of influence of these factors. The factors are divided into three levels according to the decision steps on the use of ART: 1) biological need; 2) social need; 3) technology availability.

At the first level, the patient’s medical diagnosis is determined: under the influence of environmental factors and previous reproductive behaviour, the reproductive health status is defined, and it determines the occurrence of the infertility diagnosis.

The awareness of the need for using ART occurs at the second level – the level of social need formation. Here, the decision is influenced by the woman’s age (with age the need to use ART increases if the woman does not have children), the marital status of the woman (the decision of the woman can be influenced by the partner’s opinion and the need to realize oneself as a mother). Through attitude to the technology the formation of social need is influenced by education (how much a woman trusts the technology) and religion (the acceptability of medical intervention in the process of conception).

The third level is the accessibility level. Factors of the third level have an impact on the ability to take advantage of the technology if a woman has already decided that she would like to use ART. These factors can be influenced by state policy measures, which include: level of development of the mother and child health care system which show how safe and effective the procedure can be for a woman; economic factors including level of income and income inequalities showing the financial availability of the technology; social factors including information and territorial accessibility of ART services.

We are particularly interested in the impact of economic and social factors on the use of ART in the regions of Russia and the combination of the influence of these factors with the influence of fertility policies increasing the availability of ART.

State support for IVF in Russia began with the inclusion of this procedure into the programme of state guarantees with funding from the federal budget and regional budgets in 2013. In fact, this change took its effect in 2014, which is reflected by the frequent mention of this year in news reports. The inclusion of IVF in the programme of state guarantees is due to the fact that, firstly, since 2014 the MHI programme has begun to include the most commonly used methods of high-tech assistance to enhance their accessibility, and secondly, the national project Demography identified the number of IVF procedures for the future. Analysis of the information published at the websites of specialized clinics and departments showed that the IVF programme funded by MHI has not been launched at the same time in all Russian regions: in some places it started in 2013, while in others – in 2014. This might be an additional aspect worth of studying while estimating the impact of IVF inclusion in MHI, but it requires a detailed analysis of territorial programs of state guarantees in all regions.

Since 2016 IVF is included in the basic MHI programme with funding from the FMHI, i.e. the procedure can be arranged by prescription of the gynecologist and covered by the MHI policy. There are no restrictions on the number of IVF attempts via MHI, the only exception is due to the frequency of attempts within one year: no more than two attempts per year are allowed, which is due to medical reasons.

This study examined the change in IVF availability after incorporating it into the state guarantee Programme in 2014. It is not currently possible to analyze the change in the availability of IVF after its inclusion in the basic MHI programme in 2016, since as of carrying out this study, data from the ART register are only available for 2017, i.e. one year after the inclusion of IVF into MHI. The author of the study decided to consider 2014 the year of incorporation of IVF into the MHI due to the fact that in 2013-2014, along with other adjustments in MHI due to the expansion of the high-tech medical care, the programme itself underwent changes, so in 2013 it could not yet be seen as fully operational.

The IVF procedure is arranged so that more than one embryo is often transplanted to the patient to increase the likelihood of pregnancy. Therefore, the author supposes that the distribution of the use of the technology in regions strongly correlates with the distribution of multiple births, which we know from the Federal State Statistics Service data.

Before verifying this assumption, let’s look at the dynamics of multiple births in Russia (Figure 4). Since 2004, there has been an increase in the proportion of multiple births, and for the urban population the increase is faster than for the rural population. 2004 is the year since which the number of children born with the use of ART exceeds 100 per 1 million population. In 2012, the growth of the proportion of multiple births slowed down, which may be due to the recommendation of the Ministry of Health to transplant no more than two embryos.

The relevance of this assumption might also be verified based on the US data, which are more complete. The aggregated data on the number of live births using ART and the number of multiple births across the US is plotted in Figure 5. The data is aggregated form for the period from 2015 to 2017, firstly, to level random effects in the data, and secondly, to i account for the specifics of the available data. Pearson’s correlation coefficient is 0.93 and is significant at a 1% level, which supports the proposed hypothesis.

Having confirmed our assumption, we can construct a dependent variable – the proportion of children born using ART in the total number of live births in the region as follows:

total_IVF_children_i = pct_mult_birth_i . total_IVF_children, (1)

where total_IVF_children_i is the number of children born with the use of ART in the ith region; pct_mult_birth_i is the proportion of multiple births in the ith region in total multiple births; total_IVF_children is the total number of children born using ART in Russia.

The following regression equation is proposed to evaluate models:

Pct_IVF_children_it = Ins_{t –1} + Ins_t–1 . E_it–1 + E_it–1 + S_it–1 + K_it–1 + ε_it , (2)

where t is the period of time; i is the region; Pct_IVF_children_it is the proportion of children born using ART; Ins_t–1 is a dummy variable of inclusion of IVF in the MHI programme; E_it–1 is the matrix of economic factors; S_it–1 is the matrix of social factors; K_it–1 is the matrix of control variables; ε_it is random errors¹.

Three methods for evaluating models based on panel data were used for the analysis: a model with fixed effects of time and the region, a model with random effects, and pooled regression. We estimate four following model specifications: model 1 – pooled regression; model 2 – fixed effects of time; model 3 – fixed effects of the region; model 4 – random effects model.

According to the criterion of adjusted R-square(R²_adj), the best model among all fixed effects models is the model with the inclusion of fixed effects of the region: 0.8 versus 0.34. A model with simultaneous incorporation of effects for time and region was not constructed because of insufficient data to account for such a number of dummy variables. The results of the test for linear constraints, Breusсh–Pagan and Hausman tests, indicate the choice of the fixed effects model. This model is used for further analysis.

Table 3 provides estimates of coefficients for variables of interest. Significant factors among the variables of interest are the inclusion of IVF in MHI, the inclusion of IVF in MHI in combination with income levels, the level of income in the region, the level of education and the availability of information. In two out of four specifications, the level of income inequality is also significant.

The directions of influence for most variables correspond to the expectations of the author. The level of education has a significant positive impact on the proportion of births using ART: more educated citizens are more positive about the use of high-tech medical care. Also, the positive impact of education may be attributed to the tendency to postpone births to a later age due to the increase in the number of years of study, and as age increases, so does the need for the application of ART. The use of the Internet also has a significant positive impact on the dependent variable: greater awareness of both the service itself and the possibility of obtaining government support increases the use of ART. The level of income has a significant positive effect on the proportion of births: the appearance of a child in the family is in any case a cost, and birth using IVF requires additional costs, even if the procedure is funded within the MHI programme. The resulting signs of explanatory variables are consistent with significant valuations at the market/region level in the papers devoted to the US situation.

The key result of the analysis is an assessment of the impact of IVF inclusion in the HMI programme. Two of the three coefficients illustrating this impact were proved to be significant. Based on the model with fixed effects of the region, which was chosen when comparing the models, we see that when IVF is included in MHI (∆Ins = 1) other things being equal:

∆pct_IVF_children_i = –2,05 + 0,22 . ln(Aver_real_income_i), (3)

where ∆pct_IVF_children_i is the change in the proportion of children born using ART in the region; ln(Aver_real_income_i) is the logarithm of real average per capita income in the region.

Then ∆pct_IVF_children_i > 1 will be true for regions where ln(Aver_real_income_i) >9,32, or in translation into rubles – more than 11,159 rubles.

The findings show that inclusion in the programme had a significant positive effect on the proportion of births using ART only in wealthier regions. In poor regions, there is a negative effect from the programme (in the sample for 2017 there were 14 such regions). With the growing wealth of the region, the positive effect increases. Consequently, the inclusion of IVF in the MHI programme helped increase the use of the technology, but exacerbated inequalities in access to it.

Here we will consider how the use of alternative variables characterizing the same phenomenon will affect the results of model construction. To verify this relation we estimate six additional models using: 1) data on RAHR clinics instead of data on clinics collected by the author; 2) total incidence rate of infertility instead of the indicator of new diagnoses per year; 3) inclusion of maternal and child health rating as characteristics of the quality of system functioning; 4) and 5) environmental factors instead of infertility indicators, which is due to the influence of types of environmental pollution on male and female factors of infertility identified in medical studies; 6) level of urbanization as a proxy variable for territorial accessibility of services instead of the proportion of the population of cities with ART clinics in the total population of the region. The model with fixed effects of regions has been used as a specification to include alternative factors.

The increase in the indicator of the explanatory strength of the R²_adj model is observed only in models with the inclusion of environmental factors: an increase of 0.01. Thus, it can be concluded that environmental variables are a little better at describing the real rate of infertility incidence in regions than official statistics. In other models, alternative factors are insignificant. The following variables of interest remained significant in all specifications: the logarithm of real average per capita income (+), inequality (-), percentage of employment with higher education (+), share of households with access to the Internet (+), in 4 specifications out of 6 the inclusion of IVF in the MHI and its impact associated with the income level is significant. The results of the evaluation of the basic specification are generally stable regarding the use of alternative factors.

An analysis of information provided on the websites of regional ministries and health care departments has shown that there was mobility between nearby regions to obtain ART services: some patients are sent to other regions for treatment. Also, on the sites of some clinics there is an IVF programme according to the transport scheme. When using this programme, part of the treatment takes place in the city where the patient lives, and the IVF procedure is performed directly at a clinic in another city. In addition, there are regions where there are no IVF centers in the sample, but we still calculate the proportion of children born using ART based on the method of forming a dependent variable. It follows that the rate of births using ART also depends on where the region is located geographically and with which regions it borders. Misappropriation of spatial effects can lead to a shift in coefficient estimates.

To take into account these factors, we will use methods of modelling spatial relationships. To do this, we shall carry out the following steps:

1. We shall choose a method to form a matrix of spatial weights.
2. We shall carry out statistical tests for the presence of spatial dependence in the data.
3. If there is a significant relationship, we shall choose a method for estimating spatial panel regression.

A key element in the analysis of spatial relationships is the matrix of spatial weights. Several classes of methods for creating similar matrices are given in the literature.

The basis of the first group of methods is geographic distances between objects. Distances between the centers of spatial objects, calculated as the distance between points on the plane with latitude and longitude coordinates (for objects in one region/district), as the distance between points on a circle (for countries/regions), as the distance by railway (Vakulenko et al. 2012) or by roads (Demidova 2013) or the length of the common borders can be used as distances.

For this group of methods, several variations in calculating diagonal elements of the spatial weights matrix (W) are used after obtaining the distance matrix between objects. The first option is calculation of the values inverse to distance (d_ij) : or inverse to the square of the distance: (Ivanova 2019; Getis, Aldstadt 2004). The second option is the use of the following metric: where γ > 0 is the degree of difference between the influence of near and remote objects (Harris et al. 2011). The closer the value of parameter γ is to zero, the stronger the influence of close objects. In any case, each element of the matrix of spatial weights i ≠ j reflects the proximity of objects: the larger it is, the greater the similarity of the location of the objects. We assume diagonal elements w_ij = 0, i = j as we are interested in the relationships with other objects.

The second group of methods is the creation of an adjacency matrix, for example, a weighted matrix by the number of neighbours. The adjacency matrix is created in different ways based on the use of graph analysis methods. A variant of the matrix, in which the presence of a boundary is denoted as 1 (Kolomac 2010) is common, after which the weighted matrix is calculated.

The third group of methods is to create a proximity matrix in the characteristic space. Distance matrices in physical space and space of characteristics are combined for analysis.

Two sets of distance data have been used to analyze spatial influence. The first array is data on distances between the centers of the regions, the second is data on distances by railway between the centers of the regions (if there are ART clinics in the region, they are certainly available in the regional center) collected by the author. If there is no railway station in the city, we used the distance to the nearest station, to which we added the distance of travel by car. For pairs of regions between which there is no railway connection, we used the distance of flight by plane. Such distance would be less than by railway or by road, and it is assumed that a rational person would rather take a flight for these pairs of regions, that is, if deciding to travel so far at all. Finally, for those regions, between which there is no air communication, a transfer flight is used: a person travels to the nearest city, which has air communication with the region that is difficult to access, then travels there by plane.

To construct a matrix of spatial weights, five options for calculating elements are used: and with values γ = {1; 0,5; 0,2; 0,1}. To improve the accuracy of the matrix, the analysis includes an adjustment matrix based on the location of clinics, where w_ij = 0 if there is no IVF clinic in either of the two regions. To obtain the final matrix of spatial weights, the balance matrix and the adjustment matrix are multiplied elementally. After multiplying, the values along the lines are weighted: the sum of the values by line is equal to 1.

Inverse matrices for distances on railways and between centers in a straight line are very similar. Matrices with large values of the gamma parameter (1 and 0.5) are very sparse, since they do not give more value to closer regions than farther ones. Therefore, in further analysis, these matrices are not used.

Several tests have been developed to verify spatial dependence (Millo 2017). We use the Moran test, or Moran index, as it is the most commonly applicable to spatial relationship analysis tasks. For the Moran test, the null hypothesis is the absence of a spatial relationship between the values of the indicator. The formula for calculating the global spatial autocorrelation index is as follows:

(4)

where w_ij is the element of the spatial weights matrix located at the intersection of the ith row and the jth column; x_i,j is the value of the variable for which the presence of spatial relationships is verified, in our case – the proportion of those born with the use of ART in the region.

Table 4 presents the results of the assessment of Moran indices by year using the weight matrices proposed in the previous section.

In 2013, there is no significant spatial relationship in the values of the proportion of births with the use of ART in the regions of Russia. In other years, a significant positive relationship is observed for all years when using inverse distances and for almost all years when using another metric for evaluating elements of the weight matrix. Moran indexes values for railway-based matrices are greater than the values for the matrix for normal distances between centers. This suggests that matrices based on railway data show a closer spatial relationship in the data. Consequently, they can be considered more accurate for capturing spatial effects. For further analysis, we shall use Matrix 2 with inverse distances across railways, since Moran index estimates for this matrix are positive for all the years under consideration which is consistent with intuition about the nature of spatial relationships.

Models on panel data with spatial effects are evaluated either using the maximum likelihood method or generalized moment method (can only be used to estimate models with spatial interaction in errors). We will apply the maximum likelihood method: it is more versatile with respect to spatial evaluations of models because it can be used to evaluate all types of models.

Spatial effects can be taken into account in three ways: including explanatory variables in the spatial lag model, spatial lag of the dependent variable, and calculation of spatial errors. According to these options, there are four main types of spatial models with fixed effects (Kosfeld 2006; Fischer, Getis 2009):

1. Spatial cross model with fixed effects (Spatial cross- regressive (SLX) model with fixed effects):

(5)

where α_i is here and further the constant for region i; β_j is here and further the vector-coefficient with explanatory variable j; x_jit is here and further an element of the matrix of explanatory variables; γ_j is the coefficient at the spatial lag of the explanatory variable; W (x_it) is spatial lag of explanatory variable j in region i at time t ; ε_it is here and further random errors.

1. Spatial lag model (SAM) with fixed effects (Spatial lag model (SAM) with fixed effects):

(5)

where is the coefficient at spatial lag; W (x_it) is the spatial lag of the dependent variable for region i at time t.

1. Spatial error model (SEM) with fixed effects:

(6)

where ρ is the spatial factor; W (ε_it) is spatial lag of the standard error in region i at time t; v_it is a random error.

1. Spatial Durbin model with fixed effects:

(7)

For analysis, we decided to build three specifications of spatial models. The first specification is a model with a spatial lag of the dependent variable: the high rate of births using IVF in the region may be, firstly, a signal of the use of the technology, resulting in a network effect (the more people use the technology, the more new people learn about it) and increased demand from those who had previously been in doubt both within the region and in neighbouring regions; it may be a signal of increased use of the technology in the region and attract “reproductive migration”, thirdly, it may reflect hidden characteristics of quality of service delivery that were not considered by us in the basic model specification. The second specification is a model with a spatial lag of the number of clinics: the more clinics in nearby regions, the higher the proportion of births using technology in the region. The third specification is the Darbin model with the inclusion of both spatial lag of the dependent and explanatory variable.

Note that the evaluation of spatial models in RStudio requires a balanced panel, so we had to exclude data for five regions from the analysis: for the Ivanovo oblast there are no data on female infertility for 2016, for the Republic of Crimea and Sevastopol the data are available until 2015, for the Republic of Ingushetia the data on female infertility for 2014 and for the Chechen Republic – Gini Index data for 2010-2011 and Internet access data for 2010-2013 are unavailable. As a result, a balanced panel for 2011-2017 in 80 regions of the Russian Federation was obtained for analysis.