key: cord-0273153-5romalzr
authors: Saez, M.; Barcelo, M. A.
title: Spatial prediction of air pollution levels using a hierarchical Bayesian spatiotemporal model in Catalonia, Spain
date: 2021-06-07
journal: nan
DOI: 10.1101/2021.06.06.21258419
sha: a8bb79e095ac17856aa5f419b363c93fa9dd812f
doc_id: 273153
cord_uid: 5romalzr

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.

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.

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 (Cameletti et al., , 2013 Pirani et al., 2013; Shaddick et al., 2013; Liang et al., 2015 Liang et al., , 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 (Cameletti et al., , 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 Cameletti et al., , 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 (Cameletti et al., , 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 (Cameletti et al., , 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 (Cameletti et al., , 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 (Cameletti et al., , 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) . All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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.

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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (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/m 3 ) (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.

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 shortterm 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) .

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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 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 km 2 ). 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)

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/m 2 ) (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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 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) .

We specified a hierarchical spatiotemporal model as follows:

At the top of the hierarchy:

where i denoted the air pollution monitoring station where the pollutant was observed; t was the time unit; # ! was the location of the station; ((. , . ) 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) (, " # was the nugget effect).

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint At the next level, we specified the following measurement equation:

where -(. , . ), 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:

where Κ ( is the modified Bessel function of the second type and order : > 0. : is a parameter controlling the smoothness of the GF, , # is the variance and = > 0, is a scaling parameter related to the range, >, the distance to which the spatial correlation becomes small. We used > = √8 : = ⁄ , where > corresponded to the distance where the spatial correlation is close to 0.1 for each : (Lindgren et al., 2011 ). = = 2C√:, where C is a parameter controlling the rate of decay of the spatial correlation as the distance

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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 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:

Short-term exposure:

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 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); . !,+ = T(-!+ ), -!+ denoted the air pollutants of interest, PM10, NO2, O3 and PM2.5; F2GGH%IJ% -,!+ 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); INMI ! was the area of the ABS i; #L_-!,. , ! + and # . denoted random effects.

In the models, we included #L_-!,3456 , #L_-!,:44; 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 , the random walk of order 1 for the Gaussian vector x is constructed assuming independent increments (R INLA project, 2021a):

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.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint ! + 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 # -./01 and # 234 , 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):

Inferences for GMRFs were made following a Bayesian perspective, using the INLA approach (Rue et al., 2009 .

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) .

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 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.,

We separately estimated each year (long-term exposure) and each week (shortterm 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 .

All analyses were carried out using the free software R (version 4.0.3), through the INLA package (Rue et al., 2009 R INLA project, 2021c) . The maps were represented using the leaflet package (Cheng et al., 2019) .

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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint for long-term exposure and weeks 37 -September 9 -to 48 -November 29, 2020for short-term exposure).

The prediction accuracy was assessed by:

-Mean absolute percentage error (MAPE)

where N was the total number of available observations in the test set; -(# ! , %)

were the pollutant measurements (at station i and time unit t) at the test set; and -_(# ! , %) were the posterior means.

-Root mean square error (RMSE)

-Actual coverage of the 95% prediction intervals All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint

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 the realization of the spatiotemporal process, /(. , . ), was specified as, All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint

where |C| < 1.

Here, it was h(# ! , %) what was assumed to be zero mean Gaussian and a Matérn covariance function:

In 

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).

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint

In contrast, the daily average of O3 increased by 3.97% in 2019 compared to [2011] [2012] [2013] [2014] [2015] [2016] [2017] [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%. 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). All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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 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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; 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) .

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 (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 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.

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). 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 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 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 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).

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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 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.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 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.

In this work, we have shown a hierarchical Bayesian spatiotemporal model that has allowed us to make fairly accurate spatial predictions with a low All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint computational cost. Our model provides predictions of both long-term and shortterm 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.

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.

Abstraction Library. R package version 1.5-8 [Available at https://CRAN.Rproject.org/package=rgdal, last accessed on February 27, 2021].

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Coverage: Actual coverage of the 95% prediction credibility intervals All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted June 7, 2021. ; https://doi.org/10.1101/2021.06.06.21258419 doi: medRxiv preprint 

Institut Català de la Salut

spTimer: Spatio-temporal Bayesian modelling using R

Barcelona Air Quality Improvement Plan

Maximum likelihood estimation of the multivariate hidden dynamic geostatistical model with application to air quality in Apulia

Comparing spatio-temporal models for particulate matter in Piemonte

Spatio-temporal modeling of particulate matter concentration through the SPDE approach

Model-based clustering for spatiotemporal data on air quality monitoring

Assessing air-quality in Beijing-Tianjin-Hebei region: The method and mixed tales of

leaflet: Create Interactive Web Maps with the

All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder

data reliability, consistency, and air quality assessment in five Chinese cities

An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach

Bayesian spatial modelling with R-INLA

Meteorological data from XEMA

A Bayesian spatiotemporal model to estimate longterm exposure to outdoor air pollution at coarser administrative geographies in England and Wales

Bayesian spatiotemporal modeling for estimating short-term exposure to air pollution in

All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

R: A language and environment for statistical computing. R Foundation for Statistical Computing

Random walk of order 1 (RW1)

Model for seasonal variation

No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

Approximate Bayesian inference for latent Gaussian models using integrated nested Laplace approximations (with discussion)

Bayesian computing with INLA: A review

Effects of long-term exposure to air pollutants on the spatial spread of COVID-19 in Catalonia

A Bayesian hierarchical model for assessing the impact of human activity on nitrogen dioxide concentrations in Europe

No reuse allowed without permission. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

A spatio-temporal model for the analysis and prediction of fine particulate matter concentration in Beijing

All rights reserved. No reuse allowed without permission. (which was not certified by peer review) is the author/funder

This study was carried out within the 'Cohort-Real World Data' subprogram of CIBER of Epidemiology and Public Health (CIBERESP).

Not applicable.