Spatial prediction of air pollution levels using a hierarchical Bayesian spatiotemporal model in Catalonia, Spain ================================================================================================================= * Marc Saez * Maria A. Barceló ## Abstract Our objective in this work was to present a hierarchical Bayesian spatiotemporal model that allowed us to make spatial predictions of air pollution levels in an effective way and with very few computational costs. We specified a hierarchical spatiotemporal model, using the Stochastic Partial Differential Equations of the integrated nested Laplace approximations approximation. This approach allowed us to spatially predict, in the territory of Catalonia (Spain), the levels of the four pollutants for which there is the most evidence of an adverse health effect. Our model allowed us to make fairly accurate spatial predictions of both long- term and short-term exposure to air pollutants, with a low computational cost. The only requirements of the method we propose are the minimum number of stations distributed throughout the territory where the predictions are to be made, and that the spatial and temporal dimensions are either independent or separable. **Highlights** We show a hierarchical Bayesian spatiotemporal model. Our model provides predictions of both long-term and short-term exposure. The computational cost is low. The model only needs a minimum number of stations being distributed throughout the territory. The other requirement of our model is that the spatial and temporal dimensions are either independent or separable. ![Figure1](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/06/07/2021.06.06.21258419/F1.medium.gif) [Figure1](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/F1) **Graphical abstract** **Software and Data availability** We used open data with free access using these sources. Air pollutants Departament de Territori i Sostenibilitat, Generalitat de Catalunya [Available at: [https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Qualitat-de-l-aire-als-punts-de-mesurament-autom-t/tasf-thgu](https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Qualitat-de-l-aire-als-punts-de-mesurament-autom-t/tasf-thgu), last accessed on March 14, 2021]. Meteorological variables METEOCAT, Generalitat de Catalunya. Meteorological data from XEMA [Available at: [https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Dades-meteorol-giques-de-la-XEMA/nzvn-apee](https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Dades-meteorol-giques-de-la-XEMA/nzvn-apee), last accessed on March 14, 2021]. AEMET. AEMET Open Data [in Spanish] [Available at: [http://www.aemet.es/es/datos\_abiertos/AEMET\_OpenData](http://www.aemet.es/es/datos_abiertos/AEMET_OpenData), last accessed on March 14, 2021]. Digitized cartography of the ABS Departament de Salut. Cartography [Available at: [https://salutweb.gencat.cat/ca/el\_departament/estadistiques\_sanitaries/cartografia/](https://salutweb.gencat.cat/ca/el_departament/estadistiques_sanitaries/cartografia/), accessed on March 14, 2021]. Code will be available at [www.researchprojects.es](http://www.researchprojects.es) Key words * Spatial predictions * Hierarchical Bayesian spatiotemporal model * Stochastic Partial Differential Equations (SPDE) * Integrated Nested Laplace Approximations (INLA). ## 1. Introduction In studies assessing the health effects of exposure to air pollution, there is the problem of how to estimate that exposure. Air pollution monitoring station locations do not usually coincide with where the majority of the subjects exposed to such pollution are found. In fact, the air pollution monitoring stations are not often distributed homogeneously in the territory under study, and it is quite usual that large areas, even some densely populated one, do not have any stations at all. Many studies use the measurements observed in the geographical region of the study to estimate, by means of point estimators, exposure levels for that entire region. The estimators most widely used are the inverse-distance weighted average and the arithmetic mean of the values of the air pollutant observed in several monitor stations, although sometimes the values of the pollutants observed in the nearest monitoring station are also used as estimators. The problem, as Wannemuehler et al. (2009) pointed out, is that when air pollution levels exhibit spatial variation across the study region, using these point estimators leads to a bias, as a consequence of ignoring the spatial structure (i.e., spatial dependence) of the data. Furthermore, when that biased estimated level is related to a health variable, this leads to an underestimation of the health effect of interest (Wannemuehler et al., 2009). There are numerous studies that propose models to estimate the levels of air pollutants, explicitly incorporating both spatial and temporal dependence (Cameletti et al., 2011, 2013; Pirani et al., 2013; Shaddick et al., 2013; Liang *et al.,* 2015, 2016; Calculli et al., 2015; Cheam et al., 2017; Mukhopadhyay and Sahu, 2018; Chen et al., 2018; Clifford *et al.,* 2019; Nicolis *et al.,* 2019; Wan *et al.,* 2021) (to refer to only some of those that have appeared in the last ten years). However, we must point out that very few studies attempt to predict air pollution levels in locations where there is no monitoring station (ie, spatial prediction) (Cameletti et al., 2011, 2013; Pirani et al., 2013; Shaddick et al., 2013; Mukhopadhyay and Sahu, 2018; Nicolis et al., 2019), or, having them, to perform out-of-sample temporal predictions (Wan et al., 2021). The spatial domain of these studies ranges from cities (Santiago de Chile - Nicolis *et al.,* 2019-; Beijing - Wan *et al.,* 2021-) to countries (EU-15 countries - Shaddick et al., 2013-), passing through metropolitan areas (Greater London - Pirani et al., 2013-) and regions (Po valley, northern Italy - Cameletti et al., 2011, 2013-; England and Wales - Mukhopadhyay and Sahu, 2018-). The pollutants that are predicted in these studies are coarse particles, PM10, those with a diameter of 10 micrometres (μm) or less (Cameletti et al., 2011, 2013; Pirani et al., 2013; Mukhopadhyay and Sahu, 2018), fine particles, PM2.5, those with a diameter of 2.5 μm or less (Mukhopadhyay and Sahu, 2018; Nicolis *et al.,* 2019; Wan *et al.,* 2021), nitrogen dioxide, NO2 (Shaddick et al., 2013; Mukhopadhyay and Sahu, 2018) and ozone, O3 (Mukhopadhyay and Sahu, 2018). Regarding the frequency at which pollutants are observed, daily data (Cameletti et al., 2011, 2013; Pirani et al., 2013; Mukhopadhyay and Sahu, 2018) dominate, although hourly data (Nicolis *et al.,* 2019; Wan *et al.,* 2021) and annual data (Shaddick et al., 2013) are also used. The models used in most of these articles, in addition to incorporating spatial and temporal dependencies, include explanatory variables among which appear, in decreasing order of the number of studies, meteorological variables (Cameletti et al., 2011, 2013; Pirani et al., 2013; Shaddick et al., 2013; Nicolis *et al.,* 2019; Wan *et al.,* 2021), other pollutants different from the one predicted (Cameletti et al., 2011, 2013), topographical variables (altitude – Cameletti et al., 2013; Wan *et al.,* 2021- and distances to sea and roads - Shaddick et al., 2013 -, and to mountains - Wan *et al.,* 2021-), site types (Pirani et al., 2013; Mukhopadhyay and Sahu, 2018), and land use variables (Shaddick et al., 2013). With one exception (Wan *et al.,* 2021), the studies use a Bayesian approach since it is the one that best allows the uncertainty of complex space-time data to be incorporated. Most of the studies that use the Bayesian approach perform the inference using the Monte Carlo Markov Chain (MCMC) (Cameletti et al., 2011; Pirani et al., 2013; Shaddick et al., 2013; Mukhopadhyay and Sahu, 2018; Nicolis *et al.,* 2019). Only one uses the Stochastic Partial Differential Equations (SPDE) representation of the INLA approximation (Cameletti et al., 2013). Using MCMC implies a high computational model complexity that, in some cases, prevents the practical application of the methods proposed by these studies. As an exception, it is worth mentioning Nicolis et al. (2019), who use the *spTimer* package (Bakar and Sahu, 2015). This package, which uses MCMC, allows large space-time data sets to be handled with fast computation and very good data processing capacity. The INLA approach is much more computationally effective than MCMC, producing accurate approximations to posterior distributions, even for very complex models (Lindgren and Rue, 2015). These few studies that provide methods for spatial prediction use a relatively large number of monitoring stations. In this study we intend to present an equally effective model that allows the use of information from a small number of monitoring stations. Furthermore, we intend to make spatial predictions with a computational cost much lower than existing methods. Specifically, our objective in this work was to present a hierarchical Bayesian spatiotemporal model that allowed us to make spatial predictions of air pollution levels in an effective way and with very few computational costs. In this work, we used the SPDE representation of the INLA approximation to spatially predict, in the territory of Catalonia (Spain), the levels of the four pollutants for which there is the most evidence of an adverse health effect: PM10, NO2, O3 and PM2.5. We performed the spatial predictions at a point level (defined by its UTM coordinates), allowing them to be aggregated later in any spatial unit required. We were especially interested in the long-term exposure to air pollutants. That is, by living in a certain area an individual is exposed to a mix of pollutants that have lasting effects on their health. We also considered the performance of our method to spatially predict short-term exposure to air pollutants, which has more temporary effects on health. ## 2. Methods ### 2.1. Data We obtained information on the hourly levels of air pollution for 2011-2020 from the 143 monitoring stations from the Catalan Network for Pollution Control and Prevention (XVPCA) (open data) (*Departament de Territori i Sostenibilitat, Generalitat de Catalunya*, 2021), located throughout Catalonia (Figure 1), and that were active during that period. The pollutants we were interested in for making spatial predictions were PM10, NO2, O3 and PM2.5 (all of them expressed as μm/m3) (air pollutants of interest, hereinafter). However, the monitoring stations also measured other pollutants: nitrogen monoxide (NO), sulphur dioxide (SO2), carbon monoxide (CO), benzene (C6H6), hydrogen sulphide (H2S), dichloride (Cl2), and heavy metals (mercury, arsenic, nickel, cadmium and lead). We have used these other pollutants as covariates. ![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/06/07/2021.06.06.21258419/F2.medium.gif) [Figure 1.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/F2) Figure 1. **Location of Catalonia (Spain).** Not all pollutants of interest were measured at all the monitoring stations. Thus, during the entire period 2011-2020, PM10 was measured at 122 stations, NO2 at 77 stations, O3 at 62 stations and PM2.5 at 42 stations. As can be seen in Figure 2, most of the monitoring stations were located in the city of Barcelona and in its metropolitan area. In the rest of the territory, the stations were located in cities (especially those that measure NO2 and PM2.5) and, in the case of O3, also in rural areas. On the other hand, in 2020 (which we used to spatially predict short- term exposure), the number of air pollution monitoring stations dropped considerably, from 143 to 78. In particular, those stations that measured particles dropped dramatically (PM2.5 from 42 to 3, 92.88% less; PM10 from 122 to 36, 70.49% less). The number of stations that measured O3 went from 62 to 50 (19.35% less stations) and NO2 from 77 to 67 stations (12.99% less) (Table 1). ![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/06/07/2021.06.06.21258419/F3.medium.gif) [Figure 2.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/F3) Figure 2. **Distribution in the territory of Catalonia (Spain), of the air pollution monitoring stations, according to where air pollutants (PM10, NO2, O3 and PM2.5) are measured.** View this table: [Table 1.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/T1) Table 1. Description of the air pollutants. As we said, our primary interest was in spatially predicting long-term exposure to air pollutants. In this case, we used the monthly averages, after obtaining the daily averages from the hourly data, from January 2011 to December 2019. To make the spatial predictions of the short-term exposure, we used the daily averages from January 1, 2020 to November 29, 2020. We carried out the spatial predictions at a point level, with the centroids being Basic Health Areas (ABS, for its acronym in Catalan from here on). Catalan health planning defines an ABS as the elementary territorial unit through which primary health care services are organized (*Atenció Primària Girona. Institut Català de la Salut*, 2021). The ABSs are either made up of neighbourhoods or districts in urban areas or by one or more municipalities in rural areas. Their delimitation is determined by geographical, demographic, social and epidemiological factors and, in particular, based on the accessibility the population has to services and the efficiency in the organization of health resources [18]. Catalonia is divided into 376 ABSs, with a population between 371 and 72,321 inhabitants (mean 20,266 inhabitants, standard deviation 13,391, median 18,457 inhabitants, first quartile -Q1- 10,554, third quartile -Q3- 27,529). The population density was in the range of 0.31-34,590.58 inhabitants/km2 (mean 3,486.36, standard deviation 6,719.23, median 309.18, Q1 44.83, Q3 3,752.54). In Catalonia, 769 of the 947 of the municipalities belong to a single ABS. Of the 178 remaining, 46 were divided into more than one ABS, 37 of them into a maximum of five ABSs, eight between six and 14 ABSs and one, (the city of Barcelona) into 67 ABSs (Idescat, 2021). Less than a third of ABSs have at least one air pollution monitoring station (105 from a total of 376). An ABS has five monitoring stations, six ABSs have three stations, 22 ABSs have two stations and the remaining 76 have only one station. As covariates, we included the altitude of the air pollution monitoring station (in m) and the area of the ABS (in km2). The altitude (as well as other information related to the monitoring station, such as its latitude and longitude) were obtained from the *Departament de Territori i Sostenibilitat* (2021). We transformed the geographic coordinates (latitude and longitude) to UTM coordinates (in km) using the R package *rgdal* (Bivand *et al*., 2021). The areas of the ABS, as well as the UTM coordinates of their centroids, were calculated using QGIS (version 2.18) from the digitized cartography of the ABS (information of 2018) (open data) (*Departament de Salut*, 2021). It is known that, at least in the short term, exposure to air pollution is correlated with various meteorological variables. For this reason, in the case of spatial prediction of short-term exposure, we also included several meteorological variables as covariates. Most of them, such as temperature (in °C), relative humidity (in %), wind speed at 10m (in m/s) and atmospheric pressure (hPA) influence the dispersion of the pollutant; although some also influencing its formation, for instance, global solar radiation (W/m2) (O3 is a secondary pollutant, formed when the two atoms that make up oxygen gas dissociate under the action of light solar). The sources of the data were the stations of the Network of Automatic Meteorological Stations (XEMA) of the Meteorological Service of Catalonia (METEOCAT) (open data). We also used the daily data from the State Meteorological Agency’s (AEMET) automatic stations. Albeit not as much as the air pollutants monitoring stations, the meteorological stations are also dispersed throughout the territory. Catalonia has 217 meteorological stations, 188 belonging to METEOCAT and 29 to AEMET. All of them measured all the meteorological variables every day. We used the same model (explained in this work) for short-term exposure to carry out the spatial prediction of the daily values of the meteorological variables at the ABS level, for the year 2020. Further details concerning this, can be found in Ribas et al. (2021). ### 2.2. Model specification We specified a hierarchical spatiotemporal model as follows: At the top of the hierarchy: ![Formula][1] where i denoted the air pollution monitoring station where the pollutant was observed; t was the time unit; *si* was the location of the station; *Y*(. , . ) the spatiotemporal process, the realization of which corresponded to the pollutant measurements (at station i and time unit t); and *ɛ*(. , . ) was the measurement error defined by a Gaussian white-noise process (i.e., spatially and temporally uncorrelated) (![Graphic][2] was the nugget effect). At the next level, we specified the following measurement equation: ![Formula][3] where *y*(. , . ), is the realization of the spatiotemporal process; *μ*(. , . ) denoted the large-scale component, depending on the covariates; and *η*(. , . ) was a spatiotemporal process. The spatiotemporal process was an independent in time Gaussian field (GF) with zero mean and a Matérn covariance function: ![Formula][4] where Κν is the modified Bessel function of the second type and order *ν* > 0. *ν* is a parameter controlling the smoothness of the GF, *σ*2 is the variance and *κ* > 0, is a scaling parameter related to the range, *ρ*, the distance to which the spatial correlation becomes small. We used ![Graphic][5], where *ρ* corresponded to the distance where the spatial correlation is close to 0.1 for each *ν* (Lindgren et al., 2011). ![Graphic][6], where *ϕ* is a parameter controlling the rate of decay of the spatial correlation as the distance ![Graphic][7] increases. Due to its computational problems, we chose to represent the GF as a Gaussian Markov Random Field (GMRF) (Rue *et al.,* 2009). GMRFs are defined by a precision matrix with a sparse structure allowing inference to be performed in a computationally effective way. We linked the GF and GMRF through the Stochastic Partial Differential Equations (SPDE) approach (Lindgren et al., 2011). The SPDE allowed us to find a GMRF, with local neighbourhood and sparse precision matrix (instead of spatiotemporal covariance function and the dense covariance matrix of a GF, respectively), that best represented the Matérn field. Further details can be found in Lindgren et al. (2011) and in Cameletti et al. (2013). We specified the large-scale component, *μ*(. , . ), as a generalized linear mixed model (GLMM) with response from the Gaussian family. Specifically, for each of the pollutants of interest (PM10, NO2, O3 and PM2.5) we specified two GLMMs: one for long-term exposure and the other for short-term exposure. Long-term exposure: ![Formula][8] Short-term exposure: ![Formula][9] where i denoted the air pollution monitoring station where the pollutant was observed (i=1,2,…143); t was the time unit (month in the case of long-term exposure, day in the case of short-term exposure); *μi,t* = *E*(*yit*), *yit* denoted the air pollutants of interest, PM10, NO2, O3 and PM2.5; *pollutantsj,it* corresponded to the pollutant j measurements at station i and time unit t. Pollutants considered were, first, the pollutants of interest other than the pollutant for which the spatial prediction was made and, second, the rest of the pollutants (i.e., NO, SO2, CO, C6H6, H2S, Cl2, mercury, arsenic, nickel, cadmium and lead); *areai* was the area of the ABS i; *sd_yi,.*, *ηi* and *τ* denoted random effects. In the models, we included *sd_yi, year*, *sd_yi,week* structured random effects, indexed on a standard deviation of the air pollutant that was being predicted, in the ABS i, during a particular year (2011 to 2018) and a particular week of 2020 (weeks 1 to 37), respectively. We chose, a random walk of order one (rw1) as the structure of the random effect. In the integrated nested Laplace approximations (INLA) approach (Rue et al., 2009, 2017), the random walk of order 1 for the Gaussian vector x is constructed assuming independent increments (R INLA project, 2021a): ![Formula][10] Following the INLA approach, when, as in our case, the random effects are indexed on a continuous variable, they can be used as smoothers to model non-linear dependency on covariates in the linear predictor. *ηi* denoted a random effect indexed on the air pollution monitoring station. This random effect was unstructured (independent and identically distributed random effects) and captured individual heterogeneity, that is to say, unobserved confounders specific to the station and invariant in time. We also included *τmonth* and *τday*, structured random effects indexed on time, in order to control the temporal dependency associated to possible seasonal effects throughout the year (long-term exposure) and throughout the week (short-term exposure). In this case, a model for seasonal variation with periodicity *m* (12 for long-term exposure, seven for short-term exposure), for the random vector (x1, x2,…,xn) (n>m) was obtained assuming that the sums were independent Gaussian with a precision Y. The density for x is derived from the n-m+1 increments (R INLA project, 2021b): ![Formula][11] ### 2.3. Inference Inferences for GMRFs were made following a Bayesian perspective, using the INLA approach (Rue et al., 2009, 2017). We started from the SPDE representation, which uses a finite element representation to define the Matérn field as a linear combination of basis functions defined on a triangulation of the domain (mesh, hereinafter). This consists of subdividing the domain into a set of non-intersecting triangles meeting in, at most, a common edge or corner (Lindgren *et al.,* 2011; Cameletti *et al.,* 2013). Then, instead of projecting the subsequent mean of the random field onto mesh nodes to target locations where we do not have observed data, we performed the spatial prediction of the random field jointly with the parameter estimation process. For this, we projected the mesh into those locations with no air pollutants observed and then we jointly computed the posterior means at all the locations (with observed and unobserved air pollutants measurements) (Krainski *et al.,* 2020). We separately estimated each year (long-term exposure) and each week (short- term exposure) and then merged every year and every week. We used priors that penalize complexity (called PC priors). These priors are robust in the sense that they do not have an impact on the results, and furthermore, they have an epidemiological interpretation (Simpson *et al.,* 2017). All analyses were carried out using the free software R (version 4.0.3), through the INLA package (Rue et al., 2009, 2017; R INLA project, 2021c). The maps were represented using the *leaflet* package (Cheng *et al.,* 2019). ### 2.4. Measures of predictive performance The predictive performance of each model was assessed by cross-validation, considering a training set (2011 to 2018 for long-term exposure, weeks 1 to 36 - January 1 to September 8, 2020 -, for short-term exposure) and a test set (2019 for long-term exposure and weeks 37 - September 9 - to 48 - November 29, 2020- for short-term exposure). The prediction accuracy was assessed by: - Mean absolute percentage error (MAPE) ![Formula][12] where N was the total number of available observations in the test set; *y*(*si,t*) were the pollutant measurements (at station i and time unit t) at the test set; and ![Graphic][13] were the posterior means. - Root mean square error (RMSE) ![Formula][14] - Correlation coefficient ![Formula][15] - Actual coverage of the 95% prediction intervals ### 2.5. Sensitivity analysis We conducted two sensitivity analyses. In the first place, we carried out a new cross-validation and, in the second place, we changed the spatiotemporal model. In both cases, we consider the spatial prediction of long-term exposure to NO2. As regards cross-validation, we considered, as training sets, five random samples from the monitoring stations in which NO2 was measured during the entire period 2011-2019. Specifically, we considered random samples of, approximately, 75% of the stations (58 out of a total of 77 stations), of 70% (55 stations), of 50% (41 stations), of 45% (35 stations), and of 20% (18 stations). As a test set, we considered the rest of the stations (19, 22, 36, 42 and 59 remaining stations, respectively). Next, we calculated the measures’ prediction accuracy (explained previously). With respect to the spatiotemporal model above, we considered an independent in time Gaussian field (GF), following Camelleti *et al*. (2013) we assumed a spatiotemporal Gaussian field that changes in time according to an autoregressive of order one (AR(1)). Returning the measurement equation {2}: ![Formula][16] the realization of the spatiotemporal process, *η*(. , . ), was specified as, ![Formula][17] where |*ϕ*| < 1. Here, it was *ω*(*si,t*) what was assumed to be zero mean Gaussian and a Matérn covariance function: ![Formula][18] In the addition, in the GLMM specification of the large-scale component, *μ*(. , . ), in the linear predictor we included structured random effects indexed on year, *τyear*, in order to capture the long-term trend. ![Formula][19] With this analysis, our objective was to compare not only the predictive performance of the model {1-2}, {4-5} with the one specified above {1-3}, but, above all, to compare the computation time in the inference of both specifications. ## 3. Results Descriptive results are shown in Table 1. Regarding long-term exposure, we observed that, with the exception of O3, the daily averages of pollutants decreased in 2019 (PM2.5 22.39% less, NO2 12.88% less and PM10 8.48% less). In contrast, the daily average of O3 increased by 3.97% in 2019 compared to 2011-2018. With regard to short-term exposure, the levels of NO2 and of the PM10 were higher from September 9 (week 37) (23.11% and 4.87%, respectively). Conversely, the levels of O3 and of PM2.5 (although in this case only measured in three stations) were lower than the levels before September 9 (19.34% and 9.81%). When we excluded the lockdown (which took place in Spain from March 14 - week 11- to June 21 - week 25 -, both 2020), the variation from September 9 changed sign for PM10, it was 7.67% lower, they were moderated for NO2 (which was 3.64% higher) and O3 (12.48% lower), while they were increased in the case of PM2.5 (13.74% lower). The measures of predictive performance are shown in Table 2. With the exception of PM2.5, the results for long-term exposure were quite good. Achieved coverages of the 95% credibility intervals for predictions were greater than 90%, correlation coefficients were greater than 0.80, and MAPEs less than 10%. Furthermore, if Table 2 is compared with Table 1, it is observed that the reduction in the variability of the spatial prediction, measured between the ratio of the RMSE and the standard deviations of the pollutants observed, was, at most, one third of the standard deviations of the pollutants during the period 2011-2018 (19.40% for NO2, 25.71% for O3 and 33.89% for PM10), again with the exception of PM2.5 (the RMSE in this case was 59.92% of the standard deviation in the period 2011-2018). Therefore, except for PM2.5, our method managed to significantly reduce the variability of the spatial prediction around fairly accurate predictions. Although quite good, note that, in relative terms, the results for PM10 were somewhat worse than for gaseous pollutants (NO2 and O3). View this table: [Table 2.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/T2) Table 2. Measures of predictive performance. Spatiotemporal process independent in time Gaussian field The poor results obtained for PM2.5 are because of its smaller sample size. Although it is true that in the period as a whole up to 42 stations measured PM2.5, the year with the lowest number of active stations was 2018 (31 stations), with the rest of the years ranging between 33 and 35 active stations. No year fell below 40 active stations for the rest of pollutants (the year with the lowest number of stations measuring PM10 was 2018 with 94 stations, while the other years ranged between 100 and 107 stations; in the cases of NO2 and O3 it was 2015 with 59 and 44 stations, respectively, with the other years oscillating between 62 and 66, and 45 and 57, respectively). Regarding the short-term exposure, first, predictive performance was worse when we did not exclude the lockdown period (which took place in Spain from March 14 to June 21, 2020) than when we did. In fact, note that predictive performance measures were much better for gaseous pollutants (NO2 and O3). The results for the coarse particles, PM10, were quite poor (we did not interpret the results for PM2.5 as it was measured in only three stations). This was likely due to the lower number of stations where PM10 was measured (36 stations, versus 67 for NO2 and 50 for O3, see Table 1). The variability of the spatial prediction was reduced much less than in long-term exposure, especially for NO2. The RMSEs were between 33.83% for O3 and 43.49% for NO2, of the standard deviations of the pollutants (excluding lockdown). The results of the sensitivity analyses, when the number of stations in the training set was greater than 40 and when the spatiotemporal Gaussian field changed in time according to an AR(1) (model {1-2}, {4-5}), were quite similar to the results for the spatiotemporal process independent in time Gaussian field (model {1-3}) and all the stations were included in the training set (Table 3). View this table: [Table 3.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/T3) Table 3. Sensitivity analyses. Measures of predictive performance When we varied the number of stations in the training set, but used every year (2011 to 2019), the predictive performance seemed to depend on the size of the sample. The more stations the training set had, the better the results, while dramatically deteriorating with a small sample size. In particular, the cut-off appears to be 40 stations. Below this, the predictive performance measures were poor. Although the predictive performance of the spatiotemporal Gaussian field model changed in time according to an AR (1) (model {1-2}, {4-5}) was very similar to that of the spatiotemporal process independent in time Gaussian field (model {1- 3}) (perhaps somewhat worse, in relative terms), the computation time was much longer. Using a 6-core Intel Core i9 (2.9 GHz 32 GB RAM), while the model inference {1-3} required on average 0.05 seconds per observation (a total of 569 seconds on average), the model {1-2}, {4-5} required 0.354 seconds (a total of 3,947 seconds), that is, seven times more computing time. The maps of the posterior means and the posterior standard deviations for 2019 (in quintiles) of the spatiotemporal process independent in time Gaussian field (model {1-3}) for the long-term exposure of PM10, NO2 and O3 are shown in Figures 3. We decided not to represent the posterior means for PM2.5 because of its poor predictive performance. The spatial distributions of the subsequent means of PM10 and NO2 were quite similar, although in the high levels of NO2, (fourth and fifth quintiles) there was somewhat more spatial variation. Note that, unlike PM10 and NO2, the lowest levels of O3 (first and second quintiles) occurred in the urban areas. As expected, the uncertainty, as measured by the posterior standard deviations, was, in general, higher in those areas with few (or no) monitoring stations. Note, however, that higher levels of air pollutants do not always coincide with higher standard deviations. ![Figure 3a.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/06/07/2021.06.06.21258419/F4.medium.gif) [Figure 3a.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/F4) Figure 3a. **Posterior mean and posterior standard deviation of PM10 for 2019. Spatiotemporal process independent in time Gaussian field** ![Figure 3b.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/06/07/2021.06.06.21258419/F5.medium.gif) [Figure 3b.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/F5) Figure 3b. **Posterior mean and posterior standard deviation of NO2 for 2019. Spatiotemporal process independent in time Gaussian field** ![Figure 3c.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/06/07/2021.06.06.21258419/F6.medium.gif) [Figure 3c.](http://medrxiv.org/content/early/2021/06/07/2021.06.06.21258419/F6) Figure 3c. **Posterior mean and posterior standard deviation of O3 for 2019. Spatiotemporal process independent in time Gaussian field** ## 4. Discussion Our results were quite good in terms of predictive performance, at least for those pollutants that were observed in more than 40 collecting stations (PM10, NO2 and O3 in long-term exposure and NO2 and O3 in short-term exposure). The current coverage of the spatial predictions of these pollutants are in line with similar studies. Using the same model and the same data (PM10), but using two different methods for the inference, Camelleti *et al*. find coverage between 0.95 and 0.97 (using MCMC) (Camelleti *et al.,* 2011) and 0.897 (using INLA SPDE) (Camelleti *et al.,* 2013). Mukhopadhyay and Sahu (2018) find coverage between 0.91 and 0.92 for the spatial predictions for O3 (in our case, 0.89 for short-term exposure and 0.945 for long-term exposure), between 0.89 and 0.90 for PM10 (in our case, 0.917 for long-term exposure) and between 0.95 and 0.965 for NO2 (in our case, 0.905 for short-term exposure and 0.963 for long-term exposure). Note: we have preferred not to comment on the results in which we found poor predictive performance. Our coverage could also be comparable to those provided by Pirani et al. (2013) for the spatial predictions for PM10 (between 0.87 and 0.93), although it should be noted that these show the coverage at 90%. The correlation coefficients between the observed levels of air pollutants and the subsequent means of the spatial predictions were higher in our case than in Camelleti *et al*. (0.863 when the inferences were made with MCMC - Camelleti *et al.,* 2011- and 0.702 when they were made with INLA SPDE - Camelleti *et al.,* 2013-, compared to 0.917 in our case), and than in Pirani et al. (2013) (between 0.73 and 0.78), in both cases for PM10. However, they were somewhat lower than in Mukhopadhyay and Sahu (2018) (0.88-0.89 for PM10, 0.92-0.94 for NO2, and 0.93-0.94 for O3). It should be said, nonetheless, that the number of observations in Mukhopadhyay and Sahu range between 56,625 (for PM10) and 100,138 (for NO2), while in our case we had 11,157 observations. The reduction in the variability of the spatial prediction, can only be compared with Mukhopadhyay and Sahu (2018), since they are the only ones who show these standard deviations. In this sense, both Mukhopadhyay and Sahu and ourselves achieved a similar reduction in the variability of the spatial prediction. Although good, the results of the predictive performance were less so for the spatial prediction of long-term exposure to PM10 (although it was being observed in the largest number of collecting stations, see Table 1) and for the short-term exposure for gaseous pollutants (NO2 and O3). Regarding the spatial prediction of long-term exposure to PM10, we believe that it is a consequence of the location of the monitoring stations. The stations that measure PM10, although more abundant in urban areas, are also located in rural areas, while those that measure NO2 are located almost exclusively in urban areas. In the city of Barcelona, while 13% of NO2 is generated outside the municipality, it is 71% in the case of PM10 (Barcelona City Council, 2021; Saez *et al.,* 2020). It is not unreasonable to suppose that these figures can be extrapolated to the entire Barcelona Metropolitan Area, which comprises 41.75% of the total population of Catalonia and where the majority of PM10 and NO2 monitoring stations are located. In other words, while NO2 monitoring stations measured almost all NO2 pollution, PM10 monitoring stations did not collect all PM10 pollution data. This could also explain why the posterior means of the PM10 predictions exhibited less spatial variability than the NO2 predictions (Figures 3a and 3b). With regard to the spatial predictions of short-term exposure, the reduction in the number of monitoring stations during 2020 could have led to a deterioration in the predictive performance. However, we believe it could also be due to the data behaviour during 2020. As a consequence of the lockdown to flatten the COVID- 19 pandemic curve, mobility was greatly reduced in 2020. Specifically, mobility was reduced by 40% on average, compared to pre-COVID-19 levels, during the lockdown and did not fully recover in the September-November 2020 period (being 5 to 15% lower, depending on the area of Catalonia [26]). We are sure that this anomalous behaviour would have influenced the predictive performance of the spatial predictions of short-term exposure. The predictive performance of our model depends on the number of stations where pollutants are measured. We have found that with less than 40 stations, probably spread throughout the territory (although not necessarily homogeneously), the predictive performance deteriorates considerably. Our method is quite similar to that of Camelleti et al. (2013). However, as we show with the sensitivity analysis, our method, in which we perform the inference year by year (or week by week) and then merge the subsequent ones, has a much shorter computation time, in addition to somewhat better results, even for such an atypical year as 2020. Nevertheless, we are convinced that our results might not be as good if the spatial and temporal dimensions were dependent and not separable, that is, if the spatial dependence varied over time. Fortunately, the spatial dependence of air pollutants does not vary over time. Even during 2020, although air pollution levels decreased as a consequence of the reduction in mobility, the spatial dependence was more or less similar to previous years. For spatial dependence to vary over time, major changes in infrastructures or, likewise, lasting limitations in mobility that were not homogeneous throughout the territory, have to be produced. Of course, other types of spatiotemporal data could imply other results. ## 5. Conclusion In this work, we have shown a hierarchical Bayesian spatiotemporal model that has allowed us to make fairly accurate spatial predictions with a low computational cost. Our model provides predictions of both long-term and short- term exposure. The only requirements of the method that we propose lie in a minimum number of stations being distributed throughout the territory where the prediction is to be made and that the spatial and temporal dimensions are either independent or separable. ## Funding This work was partially financed by the SUPERA COVID19 Fund, from SAUN: Santander Universidades, CRUE and CSIC, and by the COVID-19 Competitive Grant Program from Pfizer Global Medical Grants. It also received funding, in the form of a free transfer of data, from the AEMET. The funding sources did not participate in the design or conduct of the study, the collection, management, analysis, or interpretation of the data, or the preparation, review, or approval of the manuscript. ## Data Availability We used open data with free access using these sources. Air pollutants Departament de Territori i Sostenibilitat, Generalitat de Catalunya [Available at: https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Qualitat-de-l-aire-als-punts-de-mesurament-autom-t/tasf-thgu, last accessed on March 14, 2021]. Meteorological variables METEOCAT, Generalitat de Catalunya. Meteorological data from XEMA [Available at: https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Dades-meteorol-giques-de-la-XEMA/nzvn-apee, last accessed on March 14, 2021]. AEMET. AEMET Open Data [in Spanish] [Available at: http://www.aemet.es/es/datos\_abiertos/AEMET\_OpenData, last accessed on March 14, 2021]. Digitized cartography of the ABS Departament de Salut. Cartography [Available at: https://salutweb.gencat.cat/ca/el\_departament/estadistiques\_sanitaries/cartografia/, accessed on March 14, 2021]. Code will be available at www.researchprojects.es [http://www.researchprojects.es](http://www.researchprojects.es) ## Ethics Not applicable. ## Acknowledgements This study was carried out within the ‘Cohort-Real World Data’ subprogram of CIBER of Epidemiology and Public Health (CIBERESP). ## Footnotes * **Authorship:** MS had the original idea for the paper and designed the study. The bibliographic search and the writing of the introduction were carried out by MS and MAB. The methods and statistical analysis were chosen and performed by MS. MAB created the tables and figures. MS and MAB wrote the results and the discussion. The writing and final editing was done by all authors. MS and MAB reviewed and approved the manuscript. * Received June 6, 2021. * Revision received June 6, 2021. * Accepted June 7, 2021. * © 2021, Posted by Cold Spring Harbor Laboratory The copyright holder for this pre-print is the author. All rights reserved. The material may not be redistributed, re-used or adapted without the author's permission. ## References 1. AEMET. AEMET Open Data [in Spanish] [Available at: [http://www.aemet.es/es/datos\_abiertos/AEMET\_OpenData](http://www.aemet.es/es/datos_abiertos/AEMET_OpenData), last accessed on February 27, 2021]. 2. Atenció Primària Girona. Institut Català de la Salut. Basic Health Areas (ABS) [In Catalan] [Available at: [http://www.icsgirona.cat/ca/contingut/primaria/370](http://www.icsgirona.cat/ca/contingut/primaria/370), last accessed on February 27, 2021]. 3. Bakar KS, Sahu SK. spTimer: Spatio-temporal Bayesian modelling using R. Journal of Statistical Software 2015; 63(15):1–32. doi: 10.18637/jss.v063.i15. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.18637/jss.v063.i15&link_type=DOI) 4. Barcelona City Council. Barcelona Air Quality Improvement Plan (2015-2018) [in Catalan] [Available at: [https://ajuntament.barcelona.cat/ecologiaurbana/ca/que-fem-i-per-que/ciutat-productiva-i-resilient/pla-de-qualitat-de-l-aire-de-bcn](https://ajuntament.barcelona.cat/ecologiaurbana/ca/que-fem-i-per-que/ciutat-productiva-i-resilient/pla-de-qualitat-de-l-aire-de-bcn), last accessed on March 14, 2021]. 5. Bivand R, Keitt T, Rowlingson B. rgdal: Bindings for the ‘Geospatial’ Data Abstraction Library. R package version 1.5–8 [Available at [https://CRAN.R-project.org/package=rgdal](https://CRAN.R-project.org/package=rgdal), last accessed on February 27, 2021]. 6. Calculli C, Fassò A, Finazzi F, Pollice A, Turnone A. Maximum likelihood estimation of the multivariate hidden dynamic geostatistical model with application to air quality in Apulia, Italy. Environmetrics. 2015; 26(6):406–417. doi: 10.1002/env.2345. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/env.2345&link_type=DOI) 7. Cameletti M, Ignaccolo R, Bande S. Comparing spatio-temporal models for particulate matter in Piemonte. Environmetrics. 2011; 22(8): 985–996. doi: 10.1002/env.1139. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/env.1139&link_type=DOI) 8. Cameletti M, Lindgren F, Simpson D, Rue H. Spatio-temporal modeling of particulate matter concentration through the SPDE approach. AStA Adv Stat Anal. 2013; 97(2):109–131. doi: 10.1007/s10182-012-0196-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s10182-012-0196-3&link_type=DOI) 9. Cheam ASM, Marbac M, McNicholas PD. Model-based clustering for spatiotemporal data on air quality monitoring. Environmetrics. 2017; 28(3):e2437. doi: 10.1002/env.2437. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/env.2437&link_type=DOI) 10. Chen L, Guo B, Huang J, He J, Wang H, Zhang S, Chen SX, Assessing air-quality in Beijing-Tianjin-Hebei region: The method and mixed tales of PM2.5 and O3. Atmos Environ. 2018; 193:290–301. doi: 10.1016/j.atmosenv.2018.08.047. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.atmosenv.2018.08.047&link_type=DOI) 11. Cheng J, Karambelkar B, Xie Y. leaflet: Create Interactive Web Maps with the JavaScript ‘Leaflet’ Library. R package version 2.0.3. [https://CRAN.R-project.org/package=leaflet](https://CRAN.R-project.org/package=leaflet), 2019. 12. Clifford S, Low-Choy S, Mazaheri M, Salimi F. A Bayesian spatiotemporal model of panel design data: Airborne particle number concentration in Brisbane, Australia. Environmetrics. 2019; 30(7):e2597. doi: 10.1002/env.2597. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/env.2597&link_type=DOI) 13. Departament de Salut. Cartography [Available at: [https://salutweb.gencat.cat/ca/el\_departament/estadistiques\_sanitaries/cartografia/](https://salutweb.gencat.cat/ca/el_departament/estadistiques_sanitaries/cartografia/), accessed on February 27, 2021]. 14. Departament de Territori i Sostenibilitat, Generalitat de Catalunya [Available at: [https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Qualitat-de-l-aire-als-punts-de-mesurament-autom-t/tasf-thgu](https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Qualitat-de-l-aire-als-punts-de-mesurament-autom-t/tasf-thgu), last accessed on February 27, 2021]. 15. Idescat. Statistical Institute of Catalonia [Available at: [https://www.idescat.cat/?lang=en](https://www.idescat.cat/?lang=en), last accessed on February 27, 2021]. 16. Krainski E, Gómez-Rubio V, Bakka H, Lenzi A, Castro-Camilo D, Simpson D, Lindgren F, Rue H. Advanced Spatial Modelling with Stochastic Partial Differential Equations Using R and INLA. London: Chapman and Hall/CRC, 2020, chapter 2.5 [Available at: [https://becarioprecario.bitbucket.io/spde-gitbook/](https://becarioprecario.bitbucket.io/spde-gitbook/), last accessed on March 6, 2021]. 17. Liang X, Zou T, Guo B, Li S, Zhang H, Zhang S, Huang H, Chen SX. Assessing Beijing’s PM2.5 pollution: severity, weather impact, APEC and winter heating. Proc R Soc A. 2015; 471(2182):20150257. doi: 10.1098/rspa.2015.0257. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rspa.2015.0257&link_type=DOI) 18. Liang X, Li S, Zhang S, Huang H, Chen SX. PM2.5 data reliability, consistency, and air quality assessment in five Chinese cities. JCR Atmospheres 2016; 121(17):10220–10236. doi: 10.1002/2016JC024877. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/2016JC024877&link_type=DOI) 19. Lindgren FK, Rue H, Lindström J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J R Stat Soc Series B Stat Methodol. 2011; 73(4):423–498. doi: j.1467-9868.2011.00777.x. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=j.1467-9868.2011.00777.x&link_type=DOI) 20. Lindgren F, Rue H. Bayesian spatial modelling with R-INLA. Journal of Statistical Software 2015; 63(19). doi:10.18637/jss.v063.i19. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.18637/jss.v063.i19&link_type=DOI) 21. METEOCAT, Generalitat de Catalunya. Meteorological data from XEMA [Available at: [https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Dades-meteorol-giques-de-la-XEMA/nzvn-apee](https://analisi.transparenciacatalunya.cat/en/Medi-Ambient/Dades-meteorol-giques-de-la-XEMA/nzvn-apee), last accessed on February 27, 2021]. 22. Mukhopadhyay S, Sahu SK. A Bayesian spatiotemporal model to estimate long- term exposure to outdoor air pollution at coarser administrative geographies in England and Wales. J R Stat Soc Ser A-Stat Soc. 2018; 181(2): 465–486. doi: 10.1111/rssa.12299. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/rssa.12299&link_type=DOI) 23. Nicolis O, Díaz M, Sahu SK, Marín JC. Bayesian spatiotemporal modeling for estimating short-term exposure to air pollution in Santiago de Chile. Environmetrics. 2019; 30(7):e2574. doi: 10.1002/env.2574. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/env.2574&link_type=DOI) 24. Pirani M, Gulliver J, Fuller GW, Blangiardo M. Bayesian spatiotemporal modelling for the assessment of short-term exposure to particle pollution in urban areas. Journal of Exposure Science & Environmental Epidemiology 2013; 27:1–9. doi: 10.1038/jes.2013.85. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/jes.2013.85&link_type=DOI) 25. QGIS. Version 2.18 [Available at: [https://qgis.org/downloads/](https://qgis.org/downloads/), last accessed on February 27, 2021]. 26. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. [https://www.R-project.org/](https://www.R-project.org/), 2021. 27. R INLA project 2021a. Random walk of order 1 (RW1) [Available at: [https://inla.r-inla-download.org/r-inla.org/doc/latent/rw1.pdf](https://inla.r-inla-download.org/r-inla.org/doc/latent/rw1.pdf), last accessed on March 12, 2021]. 28. R INLA project 2021b. Model for seasonal variation [Available at: [https://inla.r-inla-download.org/r-inla.org/doc/latent/seasonal.pdf](https://inla.r-inla-download.org/r-inla.org/doc/latent/seasonal.pdf), last accessed on March 12, 2021]. 29. R INLA project, 2021c [Available at: [http://www.r-inla.org/home](http://www.r-inla.org/home), last accessed on March 12, 2021]. 30. Ribas V, Miralles F, Rey O, Rafael X, Subías P, Torrent M, Vicens JA, Saez M, Barceló MA, Ponce-de-León M, Valencia A, Arenas A, Saura P. Big Data i Intel·ligència Artificial per a la prevenció d’epidèmies. Monitoratge i predicció per a la detecció primerenca de brots epidemics [in Catalan] Barcelona: Generalitat de Catalunya, 2021. 31. Rue H, Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models using integrated nested Laplace approximations (with discussion). J R Stat Soc Series B Stat Methodol. 2009; 71:319–392. doi:j.1467-9868.2008.00700.x. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=j.1467-9868.2008.00700.x&link_type=DOI) 32. Rue H, Riebler A, Sørbye H, Illian JB, Simpson DP, Lindgren FK. Bayesian computing with INLA: A review. Annual Reviews of Statistics and its Applications 2017; 4(March):395–421. doi: annurev-statistics-060116-054045. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=annurev-statistics-060116-054045&link_type=DOI) 33. Saez M, Tobias A, Barceló MA. Effects of long-term exposure to air pollutants on the spatial spread of COVID-19 in Catalonia, Spain. Environ Res. 2020; 191:110177. doi: 10.1016/j.envres.2020.110177. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.envres.2020.110177&link_type=DOI) 34. Shaddick G, Yan H, Vienneau D. A Bayesian hierarchical model for assessing the impact of human activity on nitrogen dioxide concentrations in Europe. Environ Ecol Stat. 2013; 20:553–570. doi: 10.1007/s10651-012-0234-z. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s10651-012-0234-z&link_type=DOI) 35. Simpson DP, Rue H, Martins TG, Riebler A, Sørbye SH. Penalising model component complexity: A principled, practical approach to constructing priors (with discussion). Statistical Science 2017; 32(1):1–46. doi: 10.1214/16-STS576. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1214/16-STS576&link_type=DOI) 36. Wan Y, Xu M, Huang H, Chen SX. A spatio-temporal model for the analysis and prediction of fine particulate matter concentration in Beijing. Environmetrics 2021; 32(1):e2648. doi: 10.1102/env.2648. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1102/env.2648&link_type=DOI) 37. Wannemuehler K, Lyles R, Waller L, Hoekstra R, Klein M, Tolbert P. A conditional expectation approach for associating ambient air pollutant exposures with health outcomes. Environmetrics. 2009; 20(7):877–894. doi: 10.1002/env.978.36 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/env.978.36&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20161413&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F06%2F07%2F2021.06.06.21258419.atom) [1]: /embed/graphic-5.gif [2]: /embed/inline-graphic-1.gif [3]: /embed/graphic-6.gif [4]: /embed/graphic-7.gif [5]: /embed/inline-graphic-2.gif [6]: /embed/inline-graphic-3.gif [7]: /embed/inline-graphic-4.gif [8]: /embed/graphic-8.gif [9]: /embed/graphic-9.gif [10]: /embed/graphic-10.gif [11]: /embed/graphic-11.gif [12]: /embed/graphic-12.gif [13]: /embed/inline-graphic-5.gif [14]: /embed/graphic-13.gif [15]: /embed/graphic-14.gif [16]: /embed/graphic-15.gif [17]: /embed/graphic-16.gif [18]: /embed/graphic-17.gif [19]: /embed/graphic-18.gif