Spatio-temporal Joint Modelling on Moderate-Extreme Air Pollution in Spain

In mainland Spain, the $\mathrm{PM}_{10}$ concentration levels data is accessible by the Air Quality e-Reporting which is EEA’s air quality database consisting of a multi-annual time series data of air quality measurement and calculated statistics for a number of air pollutants. To work with a more robust dataset, we only retain 342 air pollution stations that have at least 60% valid observations in a year, with the geographical distribution of the stations shown in red circles embedded in the mesh constructed for the SPDE approach (Figure 1). This five-year dataset (2017–2021) with 1470 observations is divided into the training set (2017–2020; 1215 observations) and the validation set (2021; 255 observations), which are used to evaluate and compare the model fitness and predictive ability.

Figure 2 shows the temporal and spatial variations of extreme $\mathrm{PM}_{10}$ following the EEA’s air quality categories with 0–20$\mu\mathrm{g}/\mathrm{m}^{3}$ (Good), 20–40$\mu\mathrm{g}/\mathrm{m}^{3}$ (Fair), 40–50$\mu\mathrm{g}/\mathrm{m}^{3}$ (Moderate), 50–100$\mu\mathrm{g}/\mathrm{m}^{3}$ (Poor), 100–150$\mu\mathrm{g}/\mathrm{m}^{3}$ (Very poor), more than 150$\mu\mathrm{g}/\mathrm{m}^{3}$ (Extremely poor). Temporally, severe pollution seems to occur in 2017 and 2020, as indicated by the numerous monitors coloured in red (very poor) and purple (extremely poor). Spatially, compared with relatively low annual maxima recorded in the east (Valencian Community) and north (Basque Country), high $\mathrm{PM}_{10}$ concentrations are most prevalent in the centre, northwest and southeast, which correspond to the autonomous communities of Madrid, Galicia, Andalusia, and the Region of Murcia, respectively. This spatio-temporal pattern inspires our further investigation of spatio-temporal modelling taking appropriate topography covariates into account, which are stated in the following section.

A number of potential predictors are available based on prior findings in the air quality literature (Cameletti et al.,, 2013; Fioravanti et al.,, 2021; Castro-Camilo et al.,, 2021), and we choose to include a set of five main spatial and spatio-temporal varying predictors with the complete description list reported in Table 1.

In the following, we describe the selected predictors in detail.

Meteorological variables. The meteorological variables (temperature, precipitation and vapour pressure) of the monthly mean are collected from the CRU TS (Climatic Research Unit gridded Time Series; Harris et al.,, 2020) dataset and aggregated to be annual mean. Accordingly, CRU TS was first published in 2000, using ADW (angular-distance weighting) to interpolate anomalies of monthly observations onto a 0.5° grid over land surfaces (excluding Antarctica) for observed and derived variables (mean, minimum and maximum temperatures, precipitation, vapour pressure, wet days and cloud cover) with no missing values in the defined domain.

Elevation. The altitude data, height over sea level, matched with locations of all air pollution monitors, is accessible in the annual aggregated air quality values dataset provided by EEA, available at https://discomap.eea.europa.eu/App/AirQualityStatistics/index.html.

Population density. The population densities are calculated in each autonomous community. The original data is collected from the statistics report (available at https://stats.oecd.org/) of the Organisation for Economic Co-operation and Development (OECD). The OECD statistics contain data and metadata for economic and education indexes of OECD countries and some selected non-member economies.

We see in Table 2 that the correlations of potential predictors to extreme and average PM₁₀ display a similar or different direction, which also vary year by year. This will be further investigated by our spatio-temporal generalized extreme model and joint model with sharing effects, see details in both Sections 3.1 and 3.2. Furthermore, considering the correlation between location variables and meteorological variables (e.g., latitude and temperature), we adjust the location variables (longitude and latitude) as covariates to investigate the impact of meteorological factors.

The generalized extreme value (GEV) distribution is widely employed for modelling extremes in environmental science (Reiss and Thomas,, 2007), such as temperature (Cheng et al.,, 2014), precipitation (Panagoulia et al.,, 2014), air pollution (Deng and Zhang,, 2018; Martins et al.,, 2017) and sea level (Lobeto et al.,, 2018). The GEV distribution has three parameters, location parameter ($\mu$, $-\infty<\mu<\infty$), scale parameter ($\sigma$, $\sigma>0$) and tail parameter ($\xi$, $-\infty<\xi<\infty$) with the cumulative distribution function

The case $\xi=0$ is interpreted as the limit case of $\xi\rightarrow 0$, leading to the Gumbel family with distribution function

Considering the potential spatial and temporal dependence among the PM₁₀ concentrations, we use two typical approaches for measurement: Matérn covariance function (Matérn,, 1986; Guttorp and Gneiting,, 2006) with the Stochastic Partial Differential Equations approach (SPDE; Lindgren et al.,, 2011) approximation for spatial correlation and the auto-regressive dynamic model (AR; Shumway and Stoffer,, 2011) for temporal dependence. We establish two groups of generalised extreme value models with similar fixed effects and varying random effects to model the annual maxima $\mathrm{PM}_{10}$ concentrations.

Model 1: Gumbel model with fixed effect and spatio-temporal random effect. Let $y_{\text{max}}(\boldsymbol{s},t)$ denote the logarithm transform of annual maxima $\mathrm{PM}_{10}$ concentrations at location $\boldsymbol{s}\in\mathcal{S}$ and year $t\in\mathcal{T}$, where $\mathcal{S}$ is the study area and $\mathcal{T}$ is the time period in focus. Under the assumption of constant scale ($\sigma$) and tail ($\xi$) parameters, we use a linear combination of fixed effects with explanatory variables and spatio-temporal varying random effect to model the location parameter ($\mu(\boldsymbol{s},t)$) in Gumbel model below. Suppose that

and

Here, the fixed effects are associated with the vector $\mathbf{x}(s,t)$ including an intercept and the explanatory variables of location variables, meteorological variables and human-effect variables listed in Table 1, and the vector $\boldsymbol{\beta}$ corresponds to the regression coefficients associated with the fixed effects. The term $u(\boldsymbol{s},t)$ represents a spatio-temporal varying random effect that incorporates spatio-temporal interaction (Cameletti et al.,, 2013). It temporally changes according to AR(1) dynamics with autocorrelation parameter $a$ and spatial correlated and serially independent innovations $w\left(\boldsymbol{s},t\right)$.

Given two locations $s_{i}$ and $s_{j}$ separated by $h=d(s_{i},s_{j})$ (normally Euclidean) units, the Gaussian process $w(s,t)$ with mean 0 and Matérn covariance function is in the form of

where $\Gamma$ is the gamma function, $K_{\nu}$ is the modified Bessel function of the second kind, $\rho>0$ is the range parameter, $\nu>0$ is the smoothness parameter, and $\sigma^{2}>0$ for the marginal variance.

Model 2: Gumbel model with fixed effect and separated spatial and temporal random effects. Our second model is a modification of Model 1 specified in Eq.(1) with spatial and temporal random effects and separable interaction effects specified in Eq.(2), namely, we keep the Gumbel distribution assumption on $y_{\text{max}}(\boldsymbol{s},t)$ as below.

Note that $u(\boldsymbol{s},t)$ is defined the same as in Eq.(2), indicating the spatio-temporal interaction term with the Kronecker product, see e.g., Cameletti et al., (2011, 2013); Fioravanti et al., (2021). The $f(t)$ denotes the non-linear random effect in the temporal structure of AR($1$), the $w(\boldsymbol{s},t)$ is the spatially dependent only random effect with SPDE structure. Specifically, the implementation of the AR($1$) model in INLA generally assumes the Gaussian white noise with mean $0$ and precision $\tau_{AR}$. For $f(t)$ defined over the naturally binned covariate (Year), let $\boldsymbol{t}=\left(t_{1},\ldots,t_{5}\right)^{\top}$ denotes the time from the first year ($t_{1}$) to the last year ($t_{5}$),

where $-1<a<1$ is a numeric constant, the so-called autocorrelation, by which we multiply the lagged variable $f(t_{i-1})$, and $\epsilon_{i}$ denotes the unpredictable error in the form of Gaussian white noise.

Model 3: GEV model with fixed effect and spatio-temporal random effect. The generalised extreme value (GEV) models are basically following the same structure as Gumbel models except the generalized extreme value distribution for the response. To be specific, we suppose that

where the location parameter $\mu(\boldsymbol{s},t)$ is of the same form of mixed effects as in Eq.(1) and random effects in Eq. (2).

Model 4: GEV model with fixed effect and separated spatial and temporal random effect. Similar consideration of Model 2 modified from Model 1, we consider the following model in parallel with Model 3, i.e., we take spatial and temporal random effects with an interaction effect into the GEV model. Suppose that

with the location parameter $\mu(\boldsymbol{s},t)$ is of spatio-temporal structure of the form in Eq.(3).

It is worth pointing out that the well-known limit type theorem in Coles, (2001) motivates our GEV distribution assumption of annual maxima of PM₁₀, and its reduced case (i.e., Gumbel corresponds to GEV with $\xi=0$). The latter one becomes more parsimonious, and its exponential decay tail might be appropriate to fit extreme air quality (Deng and Zhang,, 2018). The Gumbel model is generally needed when the uncertainty of $\xi$ being non-zero arises in the GEV model with an estimate of $\xi$ close to zero.

In order to identify the potential varied effect levels of main explanatory variables, with inspiration from applications of the joint model with sharing effects on wildfire (Koh et al.,, 2021), we model both moderate and extreme PM₁₀ pollution simultaneously in two respective sub-models linked by the sharing effects and the sharing coefficients (scaling factors).

Let $y_{\text{mean}}(\boldsymbol{s},t)$ denote the logarithm transform of annual mean $\mathrm{PM}_{10}$ at location $\boldsymbol{s}\in\mathcal{S}$ and year $t\in\mathcal{T}$. We perform the Gaussian sub-model and Gumbel sub-model on annual mean and annual maxima simultaneously, with the structure of the best extreme value model as Model 1 according to the model fitness and prediction analysed in Section 4.1.

where the terms $\boldsymbol{\beta}^{S}$ and $u^{S}(\boldsymbol{s},t)$ with superscript $S$ denote the sharing effects, and $\mathbf{x}^{S}(\boldsymbol{s},t)$ are corresponding variables (geographical and meteorological). The term $\boldsymbol{\beta}^{N\!S}$ with superscript $NS$ denotes the non-sharing effects with corresponding covariates $\mathbf{x}^{N\!S}(\boldsymbol{s},t)$ which includes the intercept, location variables and human-effect variables (longitude, latitude, population density). For the selection of sharing and non-sharing variables, to avoid the potential issue of uncertainty of significant sharing effects (${\beta}^{S}=0$), we treat two significant predictors (altitude and precipitation) as sharing terms to investigate the potential similar and reverse effects by the ratios ($\boldsymbol{\beta}_{1}^{\text{Gaussian-Gumbel}}$) while considering all other variables, including three non-significant ones (temperature, vapour pressure, and population density), and location variables as non-sharing terms.

Additionally, the spatio-temporal random effect is also treated as sharing effects with respect to the simplicity of computation. The sharing effects $\boldsymbol{\beta}_{1}^{\text{Gaussian-Gumbel}}$ and $\boldsymbol{\beta}_{2}^{\text{Gaussian-Gumbel}}$ scale the common components of covariate vector $\boldsymbol{x}^{S}(\boldsymbol{s},t)$ and random effect $u^{S}(\boldsymbol{s},t)$, and control how much information is shared from the average predictor towards the annual max predictor, and determine the strength of interaction between the two processes. Precisely, it allows for capturing the positive or negative correlations.

For further explanation, the sharing effects (including fixed effects and random effects) are the same in two sub-models ($\boldsymbol{\beta}^{S}$, $u^{S}(\boldsymbol{s},t)$). Meanwhile, they are linked by the sharing coefficients (scaling factors) $\boldsymbol{\beta}_{1}^{\text{Gaussian-Gumbel}}$ and ${\beta}_{2}^{\text{Gaussian-Gumbel}}$. On the one hand, these sharing coefficients relax the strictly equal relation between the sharing effects in two sub-models. On the other hand, more importantly, their posterior distributions can also measure the similarities and differences in the effects of predictors in the sub-models. For instance, a significantly negative ${\beta}_{1}^{\text{Gaussian-Gumbel}}$ implies that the corresponding predictors oppositely influence the moderate and extreme air pollution cases.

In a Bayesian context, in order to finalize the model, we need to define prior distributions for the remaining parameters in Gumbel ($\sigma_{\text{Gumbel}}$) and GEV distributions ($\sigma_{\text{GEV}}$, $\xi_{\text{GEV}}$), the regression coefficients ($\boldsymbol{\beta}$), the sharing coefficients in the joint model ($\boldsymbol{\beta}^{\text{Gaussian-Gumbel}}$), parameters in the Matérn covariance function ($\sigma_{M},\rho_{M}$, $\nu_{M}$) and the parameters in AR(1) dynamic model ($a$, $\tau_{AR}$).

We use vague Gaussian priors for the tail parameter ($\xi_{\text{GEV}}$) in GEV distribution and the elements of coefficients ($\boldsymbol{\beta}$, $\boldsymbol{\beta}^{\text{Gaussian-Gumbel}}$). The smooth parameter $\nu_{M}$ is treated here as a fixed value with $\nu_{M}=1$, as in most spatial analyses. The parameters $\sigma_{M}$ and $\rho_{M}$ in the Matérn function and autocorrelation parameter $a$ in AR(1) model are defined by penalized complexity (PC) priors (Simpson et al.,, 2017) with knowledge from Moraga, (2019) and Fuglstad et al., (2019). PC prior for the range parameter ($\rho_{M}$) is defined with $\operatorname{Prob}\left(\rho_{M}<10^{4}\right)=0.01$, which means the probability that the range is less than $10\mathrm{km}$ is very small, and the PC prior for variance parameter as $\operatorname{Prob}\left(\sigma_{M}>3\right)=0.01$, indicating the probability for variance greater than $3$ is low. Similarly, we apply the auto-correlation ($a$) PC prior following the recommendation of $\operatorname{Prob}\left(a>0\right)=0.9$.

Note that INLA often uses the precision parameter ($\tau$) to replace the scale parameter ($\sigma$) by $\sigma={1}/{\sqrt{\tau}}$. The PC priors for all precision parameters are given by $\operatorname{Prob}\left({1}/{\sqrt{\tau}}>3\right)=0.01$ for Gumbel and GEV likelihood and $\operatorname{Prob}\left({1}/{\sqrt{\tau_{AR}}}>5\right)=0.01$ for the AR model.

Traditional Bayesian model performance is evaluated by two popular criteria, the deviance information criterion (DIC) and the Watanabe-Akaike information criterion (WAIC). The deviance information criterion (DIC) proposed by Spiegelhalter et al., (2002), is a popular criterion for model choice similar to the Akaike information criterion (AIC).

where $D(\widehat{\boldsymbol{\theta}})$ is the deviance function with Bayes estimate $\widehat{\boldsymbol{\theta}}$, and $p_{D}$ is the effective number of parameters. The Watanabe-Akaike information criterion, also known as the widely applicable Bayesian information criterion, is similar to the DIC, but the effective number of parameters is computed in a different way (Watanabe,, 2013).

However, DIC may under-penalize complex models with many random effects. Alternatively, for prediction performance, INLA suggests applying the leave-one-out cross-validation criteria, conditional predictive ordinates (CPO; Pettit,, 1990) with its summative version, the Logarithmic Score (LS; Gneiting and Raftery,, 2007) and predictive integral transform (PIT; Marshall and Spiegelhalter,, 2003), which facilitates the computation of the cross-validated log-score for model choice, and enables the calibration assessment of out-of-sample predictions, respectively.

where $y_{i}^{\text{obs }}$ denotes the $i$-th observation and $y_{-i}$ denotes the observations $y$ with the $i$-th component omitted. A model with a small value of LS and a standard uniform distributed PIT is preferable (Gómez Rubio,, 2020).

Note that numerical problems may occur when CPO and PIT values are computed with the complicated Gaussian process with INLA algorithm (Held et al.,, 2010). Hence, we take a few other criteria introduced in Bayesian literature into account: The coverage probability of 95% CI is computed as the proportion of the validation observations that the observed value lies between the 2.5% quantile to the 97.5% quantile of the predicted value (posterior distribution). The correlation coefficient denotes the correlation between observed values and predicted values in the validation set. The root mean square error (RMSE) is defined as the square root of the second sample moment of the differences between predicted values and observed values.

As we separate the original five-year dataset (2017–2021) into training (2017–2020) and validation (2021) sets, we decide to apply DIC, WAIC, LS, PIT and RMSE to compare the model fitness on the training set and use the coverage probability, correlation coefficient, and RMSE for predictive ability evaluation.

In order to examine the fitting and prediction ability of the extreme value models, we apply some evaluation criteria (DIC, WAIC, LS, PIT and RMSE) on the training set (Table 3), and use other criteria (coverage probability, correlation and RMSE) on the validation set (Table 4). We see from Table 3 with model performance on the training set, Models 2 and 4 outperform in DIC (36.27 and 73.43, respectively), WAIC (193.67 and 223.06, respectively) and RMSE (0.24 and 0.18, respectively). However, Model 1 is preferable in the leave-one-out cross-validation criteria with LS (8019.07) and PIT plots (Figure 3). Table 4 shows that all four models are comparable in the coverage probability, correlation and RMSE.

Figure 4 shows the simultaneous visualization of model performance on both training and validation sets, where the scatters’ distribution along the line with intercept 0 and slope 1 indicates the estimated (predicted) values are best suited to the observations. All four models exhibit a strong fit for the trend.

To summarise, no model excels in all criteria. Considering the similar performance of all four models on the validation set, we choose the parsimonious model, Model 1 (the Gumbel model with fixed effects and spatio-temporal random effects defined in Eq.(1)), as the best model and bring its structure into further analysis.

The summary statistics for fixed effects are shown in Table 5. We see that altitude and population density are significantly negatively associated with annual maxima $\mathrm{PM}_{10}$ concentrations, whereas the effects of other covariates are not statistically significant. High temperature, low precipitation and low vapour pressure are likely to associate with extreme $\mathrm{PM}_{10}$ concentrations, consistent with Kalisa et al., (2018) and Li et al., (2014).

The negative association between population density and annual maxima $\mathrm{PM}_{10}$ concentrations is different from the positive correlation found in Table 2 and the positive association demonstrated by Borck and Schrauth, (2021). Together with the facts of relatively low correlation for population density (Table 2) and the phenomena that severe $\mathrm{PM}_{10}$ pollution in 2017, 2020 and 2021 (Figure 2) is usually associated with high temperature and low precipitation, our results probably imply that the generation and spread of extreme $\mathrm{PM}_{10}$ concentrations depend heavily on the climate conditions. This evidence coincides with the findings of the overwhelming role of adverse meteorology in severe air pollution events (Wang et al.,, 2020; Morawska et al.,, 2021).

The spatial and temporal dependence is accessible by spatial heat plot (Figure 5) and posterior estimates of autocorrelation coefficient ($a$ in Table 6). Spatially, similar values of mean random effects occur in groups, especially, a large cluster of high values happens in the centre of the mainland. Temporally, the estimated mean correlation coefficient ($a$) is 0.80 with 95% credible interval (0.74, 0.85), providing evidence of relatively strong dependence between two consecutive years.

Combining the best extreme value model (Model 1) with the Gaussian model, we establish the joint model in Section 3.2 with sharing effects to estimate both maxima and mean concentration levels simultaneously. We keep the two significant effects of predictors (altitude and precipitation) as sharing effects, and all other effects (temperature, vapour pressure, population density and location) are treated as non-sharing ones. Table 7 summarizes the posterior estimates of these effects from the joint model. Precipitation shows significant negative associations with both annual mean and maxima concentrations, which is consistent with the notion that $\mathrm{PM}_{10}$ concentrations are generally low in wet regions (Li et al.,, 2014). In contrast, we find that altitude is negatively connected with the mean but positively connected with the maxima.

This is also supported by the detailed sharing coefficients analysed below. We see from Figure 6 that the posterior distribution of the sharing coefficients of precipitation (see detailed definition of $\boldsymbol{\beta}_{1}^{\text{Gaussian-Gumbel}}$ in Section 3.2) almost lies between 0 and 1, showing that the influences of precipitation are similar on extreme and moderate air pollution. The coefficient for altitude is negative, suggesting certain evidence that altitude can impact reversely in different scaled pollution.

Furthermore, the effect of temperature is not significant with respect to 95% credible interval for mean but becomes significant for maxima. This result implies that high temperatures weakly impact the general air quality, while playing an essential role in the dispersion of extremely poor pollution. Nevertheless, considering the possibility of unmeasured confounders potentially distorting the coefficients, the strength of the aforementioned inference might be compromised.

In terms of model performance evaluation, the joint model demonstrated promising results in the validation set. In particular, the Gumbel sub-model (max sub-model) exhibited improvements in the coverage probability ($88.63\%$) and correlation ($46.71\%$), compared with the univariate Gumbel model (Model 1), as shown in Tables 4 and 8. The insights of the joint model performance are depicted in Figure 7, where both the performance of the Gumbel model (maxima model) and the Gaussian model (mean model) are satisfied. Because of this, despite the less parsimonious nature of the joint model, which may introduce more volatility in the prediction of the Gumbel sub-model, the overall satisfying evaluation results bolster the credibility of the joint model’s findings. With this confidence, we proceed to utilize the joint model for subsequent predictions using excursion functions.

Hot-spot region identifications and predictive concentration level plots are important tools in air pollution studies, because they intuitively imply the air quality in specific regions and are easily interpreted to environmental agencies and the general public. Bolin and Lindgren, (2015) proposed the positive ($\mathrm{E}_{u,\alpha}^{+}(X)$) and negative excursion sets ($\mathrm{E}_{u,\alpha}^{-}(X)$) that determine the largest set that simultaneously exceeds or below the risk level ($u$) with a small error probability ($\alpha$), employing a parametric family and sequential importance sampling method for estimating joint probabilities. To visualize the excursion sets simultaneously, we apply the positive and negative excursion functions, $F_{u}^{+}(s)=1-\inf\left\{\alpha\mid s\in\mathrm{E}_{u,\alpha}^{+}\right\}$ and $F_{u}^{-}(s)=1-\inf\left\{\alpha\mid s\in\mathrm{E}_{u,\alpha}^{-}\right\}$. To explain, the term $\inf\left\{\alpha\mid s_{0}\in\mathrm{E}_{u,\alpha}^{+}\right\}$ denotes the ”smallest” $\alpha$ required that the location ($s_{0}$) can be included into the positive excursion set $\mathrm{E}_{u,\alpha}^{+}$ at the first time, while the higher $1-\inf\left\{\alpha\mid s_{0}\in\mathrm{E}_{u,\alpha}^{+}\right\}$ reported by positive excursion function generally indicates higher probabilities for the location ($s_{0}$) to exceed the risk threshold $u$ simultaneously.

In our case, to discover the areas that are most likely and most unlikely to suffer from severe $\mathrm{PM}_{10}$ pollution simultaneously, we utilize predicted $\mathrm{PM}_{10}$ concentration levels from the max sub-model of the joint model to generate both positive and negative excursion functions with the thresholds 50$\mu\mathrm{g}/\mathrm{m}^{3}$ (poor) and 100$\mu\mathrm{g}/\mathrm{m}^{3}$ (very poor) at 548 locations distributed throughout the mainland of Spain, including a set of 0.5° $\times$ 0.5° (50km $\times$ 50km) grids (206 locations) and locations of all $\mathrm{PM}_{10}$ stations (342 monitors).

In the case of exceeding 50$\mu\mathrm{g}/\mathrm{m}^{3}$ (Figure 8), the probability for simultaneously exceeding is high in the northwest, middle, and south, meaning that poor $\mathrm{PM}_{10}$ pollution probably hazard these locations during the year. In contrast, the negative excursion function with threshold 50$\mu\mathrm{g}/\mathrm{m}^{3}$ indicates that the regions in the north and east enjoy good or moderate air quality throughout the whole year. In the case of 100$\mu\mathrm{g}/\mathrm{m}^{3}$ (Figure 9), the probabilities for very poor $\mathrm{PM}_{10}$ pollution occurrence are low (in white) in most regions, but still likely to appear in certain areas in the community of Madrid, meanwhile, most areas in the north, northeast and east are expected to be below this threshold.