Spatio-temporal impact of self-financed rotavirus vaccination on rotavirus and acute gastroenteritis hospitalisations in the Valencia region, Spain

Background Several studies have shown a substantial impact of Rotavirus (RV) vaccination on the burden of RV and all-cause acute gastroenteritis (AGE). However, the results of most impact studies could be confused by a dynamic and complex space-time process. Therefore, there is a need to analyse the impact of RV vaccination on RV and AGE hospitalisations in a space-time framework to detect geographical-time patterns while avoiding the potential confusion caused by population inequalities in the impact estimations. Methods A retrospective population-based study using real-world data from the Valencia Region was performed among children aged less than 3 years old in the period 2005–2016. A Bayesian spatio-temporal model was constructed to analyse RV and AGE hospitalisations and to estimate the vaccination impact measured in averted hospitalisations. Results We found important spatio-temporal patterns in RV and AGE hospitalisations, RV vaccination coverage and in their associated adverted hospitalisations. Overall, ~ 1866 hospital admissions for RV were averted by RV vaccination during 2007–2016. Despite the low-medium vaccine coverage (~ 50%) in 2015–2016, relevant 36 and 20% reductions were estimated in RV and AGE hospitalisations respectively. Conclusions The introduction of the RV vaccines has substantially reduced the number of RV hospitalisations, averting ~ 1866 admissions during 2007–2016 which were space and time dependent. This study improves the methodologies commonly used to estimate the RV vaccine impact and their interpretation.

The World Health Organization (WHO) recommended including RV vaccination worldwide. The recommended schedule is two (RV1) or three (RV5) oral doses and should be completed between 6 and 32 weeks of age. Currently, 98 countries have introduced RV vaccines into their national immunisation programs [3]. This measure has had a major impact on the burden of AGE, decreasing RV outpatient visits and hospitalisations by 60-90% in Europe [4][5][6][7], Although in Spain RV vaccines are recommended by the Spanish Paediatric Association but not funded by the National Health System (NHS), several post-authorization studies have also shown their effectiveness and impact on AGE and RV-AGE hospitalisations [8][9][10][11][12]. The Valencia Region of Spain could show a specific coverage-related impact of RV vaccines on AGE and RV-AGE hospitalisations and costs, despite the low-medium vaccine coverage (40-50%) [8].
Following WHO recommendations, most postauthorization studies usually estimate impact of the RV vaccine by comparing trends of RV or AGE hospitalisations in pre-and post-vaccination periods [7,13,14]. However, this ecological design is highly prone to bias and confounding [15][16][17].
In fact, a number of key studies have shown that the spread of infectious diseases are heterogeneously distributed in space because places differ in their environmental and population characteristics [18,19]. Consequently, epidemiological studies are often confounded by complex and dynamic spatio-temporal processes [18,20]. RV vaccine uptake and hospitalisations could, therefore, vary from time to time and between places for different reasons, including complex interaction of population demographics, socioeconomic inequalities, environmental factors, circulation of RV strains and their interactions across space and time [21]. Spatial variation in RV vaccination coverage [22] and in RV hospitalisations has been previously shown in the USA, Germany, Brazil, New Zealand [23][24][25].
A previous study in Spain showed strong variability in both vaccination coverage and RV/AGE hospitalisation rates over time and between health departments [8]. Thus, it would be important to evaluate variations in the RV/AGE hospitalisation risk and the impact of RV vaccination in a space-time framework to detect geographicaltime patterns while avoiding the potential confusion caused by population inequalities in the impact estimates [7,8,12,18,24,26].
Our aim is to assess the spatio-temporal impact of RV vaccines on RV and AGE-associated hospitalisations in children under 3 years of age in the Valencia Region using real-world data. In this study, real space-time rotavirus vaccination impact is predicted in terms of number of averted hospitalisations.

Setting and study population
This is a retrospective, population-based study using real-world data from the Valencia Region, including all children less than 3 years old living in the Region between 2005 and 2016.
The Valencia Region of Spain has approximately 4, 900,000 inhabitants. Of them, around 3% (~150,000 children) are younger than 3 years old. The regional health system is divided into 34 public hospitals (24 of them with paediatric emergency rooms) and 241 health care districts structured into 24 health departments. As RV vaccines are administered to infants from six weeks of age, children with the first dose of RV vaccine recorded before six weeks of age were excluded from the study.

Data sources
The Valencia Region has a set of multiple electronic databases collecting health and sociodemographic data from 98% of the population [27]. The population information system (SIP) was used to determine the population and their socio-demographics characteristics (sex, date of birth, health department, and health care district). Health care district and department are assigned by place of residence. Hospitalisations were collected from the minimum basic data set (MBDS). The vaccine information system (SIV) was used to obtain the vaccinated population; this source captures the immunisation history of each individual. Population, hospitalisation, and vaccination data were linked at individual level through a unique personal identification number [28].
Vaccination status was assessed as a time-varying variable. Children were considered vaccinated from the date of the first dose of RV5 or RV1 and unvaccinated before that date. Children with no recorded rotavirus vaccination in SIV were considered as unvaccinated.
Vaccination coverage was calculated as the proportion of the children < 3 years old vaccinated with at least one dose of RV1 or RV5.

Spatio-temporal analyses
The database for the analysis gathered population and hospitalisations aggregated by vaccination status, sex, age, health department, biennial periods (two-years period), and health care district.
A Bayesian spatio-temporal ecological model was constructed to analyse RV and AGE hospitalisation rates and to estimate the impact of vaccination on hospitalisations.
The model assumed that the number of hospitalisations (for RV or AGE) in the different observation units, Y = {y 1 , …, y vsadtm , …, y n }, followed a binomial distribution, where "v" indexes the two vaccination status, "s" the two sexes, "a" the 3 age groups (0, 1 and 2 years old), "d" the 24 health departments, "t" the 6 (biennial) periods, and "m" the 241 health districts. From now on, we will index y by y i instead of y vsadtm where i spans all the values of the sub-indexes v, s, a, d, t and m to make the notation shorter. Thus, the model assumed proceeds as follows: Where θ i is the hospitalisation rate and N i the population for each observation unit. θ i was modelled considering the logit link as follows: where logð δ m 1 − δ m Þ acts as an offset term to control for the hospital attraction (people who live near the hospital are more frequently admitted to it than those who live far from hospital, (see Additional file 1)), where δ m is the estimated hospitalisation rate for all causes measured in each health care district. This rate was estimated using the spatial Besag-York-Mollié model [29] on hospital admissions for any cause. This offset makes that if no other term in the linear predictor had an effect, the corresponding risk, θ i , would be that corresponding to general hospital admissions for that health care district. β 0 is the intercept term and β j are the parameters associated with the categories of the covariates, X j : vaccination status, sex and age. The health department random effect, α d , was considered to fit the differences in admission policies between hospitals. α d was considered to have the following distribution where σ is also estimated within the model. No spatial dependence was considered for this term because it is expected to fit the admission policies of each hospital, which should not follow any spatial pattern. The biennial period effect, u t , was introduced to control the expected temporal variability in RV and AGE incidence. It was modelled as a random effect considering correlation between adjacent periods by a first order random walk modelled as an intrinsic conditional autoregressive (ICAR) prior distribution. Besides the temporal and spatial (health department) terms already mentioned, it was considered appropriate to include a spatio-temporal term that could jointly vary in time and space. The random effect v tm reproduces this effect. This term is assumed to follow a spatio-temporal autoregressive model [30]. Thus, the spatio-temporal effect for the first period was formulated as and for the following periods where W tm follows a spatial Besag, York and Molliè model [29] for each time period t inducing spatial dependence on v t m . On the other hand, ρ controls the temporal dependence in v t m . This parameter is assumed to follow a uniform prior distribution between − 1 and 1. Non-informative flat prior distributions were considered for β j ( j = 0, . . , 3) parameters. Uniform prior distributions between 0 and 5 were considered for the standard deviations of all the random effects in the model. Predictive distributions were used to estimate the number of rotavirus hospitalisations averted in order to assess the impact of rotavirus vaccination by health care district and time period. The number of cases averted by vaccination was calculated as the difference between the hospitalisations predicted by the adjusted model without the vaccine effect and the hospitalisations predicted by the model explained above. R (Foundation for Statistical Computing, Vienna, Austria) and WinBUGS (Cambridge Biostatistics Unit and the Imperial College School of Medicine, London) software were used to perform the analysis using MCMC methods. A total of 2000 initial iterations were used as burn-in period of the MCMC. Subsequently, 10, 000 iterations were run and only 1 in every 10 of them was saved. Three chains were simulated in total. MCMC convergence was assessed by visual inspection of history plots of posterior samples, the Brooks-Gelman-Rubin scale reduction factor, and the effective sample size implemented in the R2WinBUGS package of R. All statistical analyses conducted for this study are completely reproducible, and the data and the R code used for statistical analysis can be found as supplemental digital content to the paper.

Results
The study included 721,471 children < 3 years old. Of these, 189,247 were vaccinated against RV. There were a total of 17,482 AGE hospitalisations, of which 28% (4871) were codified as RV. AGE and RV hospitalisations accounted for 8.4 and 2.4% respectively of all hospitalisations (207,014 hospitalisations for any cause). Vaccinated children accounted for 2248 AGE and 200 RV admissions.

Spatio-temporal hospitalisation rate and relative risk
Risk of RV and AGE hospitalisations decreased with RV vaccination (Table 1). RV and AGE hospitalisation rates were 86% (95% CI: 84-88) and 47% (95% CI: 45-50) lower in vaccinees, respectively. Risk of RV and AGE hospitalisation also decreased with increasing age, by 72% (95% CI: 70-74) and 58% (95% CI: 56-60) respectively in two-year-old children as compared to those aged less than one year old. Risk of RV and AGE-hospitalisation was respectively 19% (95% CI: 15-23) and 15% (95% CI: [12][13][14][15][16][17][18] lower in girls as compared to boys. A strong variability in both RV and AGE hospitalisation rates was found between health departments (Additional file 2). Risk of AGE hospitalisation showed a downward trend during the study (Additional file 2), while the RV rate only declined between 2005 and 2010. Once controlled the vaccine effect, RV peaked in 2013-2014, with an 8% (95% CI: 6-14) higher rate than the average risk for the whole study period (Additional file 2). Additional structured spatio-temporal interaction was found for both outcomes. The spatio-temporal effect maps (Additional file 2) showed spatial clusters after adjusting for confounders.

Spatio-temporal RV vaccination coverage
Rotavirus vaccination coverage varied considerably across the Valencia Region during the study period, with pockets of undervaccination (lower coverages) in many health care districts. Vaccination rates increased over the years in the districts. In 2016, 50% of the health care districts had a coverage higher than 53% (IQR: 35-64%) (Fig. 1). The overall RV vaccination coverage increased from 0 to 49% during the study period.

Spatio-temporal RV vaccination impact
The number of hospitalisations averted by vaccination was coverage-dependent (   These results are in accordance with the vaccine effectiveness estimated in the Valencia Region previously [9]. Regarding the spatio-temporal results, substantial variability was seen in RV vaccine coverage and hospitalisation risk for RV and AGE among health departments and health care districts. Spatio-temporal clusters were clearly distinguished. These patterns could be explained by climatic, environmental, sociodemographic, or economic differences, or by the different admission policies of health departments. Although other impact studies reported relevant reductions in both RV and AGE hospitalisations in children < 5 years following RV vaccination [4,6,7,13,14,[31][32][33], only two of them showed a coverage-dependent response [8,34]. Moreover, many of them were timetrend ecological studies comparing hospitalisation data in pre and post-vaccine populations and a historical prevaccine group [7,13,14,33]. Even though this is the most commonly used method, it has been associated with potential confusion bias [15,16]. The reported impact of the vaccination could be due to other secular trends caused by, changes in reporting, in medical practices, in health seeking behaviour, etc. [35]. Besides, vaccine impact based on hospitalisation data is prone to confounding, because hospitalisations rates are closely related to changes in the quality, access and use of the health care system which often occur simultaneously with introduction of new vaccines [17]. On the other hand, few spatial and spatio-temporal models have studied RV and AGE dynamics and none of them included the vaccination status of the population.  Spatial variation in RV hospitalisations explained by sociodemographic characteristics of the population has previously been shown in studies conducted in Germany and New Zealand [23,24]. Other studies in the USA and Brazil found that spatio-temporal variation in birth rate can lead to secular changes in the RV pattern [21,25]. Finally, a study conducted in Bhutan showed that rainfall and temperature explain much of the spatio-temporal dynamics of diarrhoea (possibly due to RV infection in approximately 23% of cases) [31]. The studies developed in Germany and New Zealand were based in aggregated data over time, however, caution should be taken when interpreting this analysis because the area-specific risk may be overestimated or underestimated. Furthermore, none of these standard models considered spatiotemporal dependence; however, what occurs in a health care district is intimately related to what occurs in the adjacent one and is also related to what happened previously [36]. The present study analysed the impact of RV vaccination on RV and AGE hospitalisations from a different point of view. We developed a sophisticated spatiotemporal model that allowed us to estimate the RV vaccination impact in terms of adverted hospitalisations according to the number of children vaccinated. The spatio-temporal approach improves the commonly used methodologies to estimate the RV vaccine impact and its interpretation as follows. First of all, this analysis showed the evolution of the impact of RV vaccination and the risk of hospitalisation for RV/AGE in the Valencia Region at the health care district level over time. Second, adjusting by spatial variables such us health care district and health department in the analysis, several potentially attributable biases can be controlled. Those biases could have been caused by economic inequalities, environmental factors, socio-demographic differences or even possible changes in hospitalisations-admission policies [21,[37][38][39]. Moreover, the hospitalisation rate for any cause of each health care district was included to adjust the confusion caused by hospital attraction or other secular trends [17]. Finally, the Bayesian approach used allowed us to adequately capture dependencies among health areas and the potential relationship of data over time that cannot be easily modelled in classical statistics [40,41].
Nevertheless, some limitations of our study should be highlighted. First of all, RV vaccines are not included in the official immunisation schedule, which may suggest differences between rotavirus vaccinees and nonvaccinees in terms of socioeconomic conditions and health-seeking behaviour. Therefore, socioeconomic factors might be an important confounder of our results and admissions at private hospitals should also be considered in future studies.
Secondly, although the positive predictive value of the rotavirus ICD-9-CM code identifying acute gastroenteritis attributable to rotavirus using MBDS resulted in 90% [9], different immunochromatographic methods with different sensitivities and specificities could have been used in the different hospitals during the study period [42]. In fact, based on the difference found in the number of hospitalisations prevented for AGE and RV (1866 vs. 1142),~40% of underdiagnosis in RV hospitalisations was detected in the present study. Thirdly, health care district and health department could have varied over time; but only the last updated information was available. Fourthly, children who were unable to receive RV vaccines according to manufacturer recommendations (i.e. immunocompromised children) were not excluded from the analysis due to the lack of information.
Finally, it should be noted that both vaccines (RV1 and RV5) were used concurrently until 2010. But, RV5 was the only rotavirus vaccine available in Spain between 2010 and 2016. Therefore, results will have a limited value for estimating the impact of RV1.

Conclusions
In summary, the introduction of the RV vaccines has substantially reduced the number of RV hospitalisations. The sophisticated spatio-temporal analysis allows us to show the impact of different vaccine coverage rates in terms of avoided hospitalisations in a geographical-time framework. Interestingly, our study predicted that8 606 RV hospitalisations could have been adverted with all children vaccinated. This study improves the methodologies commonly used to estimate the RV vaccine impact and its interpretation. The spatio-temporal model avoided the potential confusion caused by population inequalities in the impact estimations. It also detects spatial clusters of the RV and AGE-hospitalisation risk attributable to common environmental, demographical, or cultural effects shared by neighboring regions.