# Hierarchical Bayesian Spatio-Temporal Modeling for PM10 Prediction

Authors: Esam Mahdi; Sana Alshamari; Maryam Khashabi; Alya Alkorbi
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

Over the past few years, hierarchical Bayesian models have been extensively used for modeling the joint spatial and temporal dependence of big spatio-temporal data which commonly involves a large number of missing observations. This article represented, assessed, and compared some recently proposed Bayesian and non-Bayesian models for predicting the daily average particulate matter with a diameter of less than 10 (PM10) measured in Qatar during the years 2016–2019. The disaggregating technique with a Markov chain Monte Carlo method with Gibbs sampler are used to handle the missing data. Based on the obtained results, we conclude that the Gaussian predictive processes with autoregressive terms of the latent underlying space-time process model is the best, compared with the Bayesian Gaussian processes and non-Bayesian generalized additive models.


---

## Body

## 1. Introduction

Many of the environmental data contain different scales of variability over space and time. For example, scientists from environmental and public health sciences are typically interested in modeling the evolving process of the air pollution during the time over specified locations. Such a stochastic process is often high-dimensional, large, and complicated with nonstationary structures, so the traditional statistical methods are hampered by the need of advanced statistical techniques to specify the spatio-temporal dependency. This can be quite practical with modern computers with high-level computational programming. The spatio-temporal modeling ofPM10 and PM2.5 (particulate matter with diameters of less than 10 and 2.5 micrometers, respectively) is rapidly becoming an important component of most air quality studies [1–3]. Particulate matter (also called particle pollution) is a mixture of solid particles and liquid droplets found in the air as a result of dust, soot, dirt, smoke caused by road transportation, and complex chemical reactions in the atmosphere such as sulfur dioxide and nitrogen oxides. Exposure to particle pollution is a public health hazard and can cause acute and chronic heart and lung diseases [4]. The larger the values of particulate matters (PM) are, the more harm on short and long terms of public health is. World Health Organization’s (WHO) air quality guidelines recommend that the annual and 24-hour mean concentrations should not, respectively, exceed 20 and 50 microgram per cubic metre (μg/m3) for PM10, and 10 and 25 μg/m3 for PM2.5. Countries with fast-developed infrastructure such as Qatar usually suffer from relatively high levels of PM air pollutant. In this regard, the WHO classified the air quality in Qatar as poor and unsafe. In this paper, we focus our attention on modeling the PM10 in Qatar because the PM2.5 data is inaccessible. The most recent data indicates that country’s annual average 24 hr PM10 concentration levels were ranged from 126.69 μg/m3 to 184.55 μg/m3, which exceeds the recommended maximum of 50 μg/m3 [4–6].Several authors have developed spatio-temporal models for analyzing the ambient of air pollution. Research in this field is back in dates to Cressie [7] and Goodall and Mardia [8]. Cressie et al. [9] compared the performance of Markov-random field with the familiar geestatistical approach in prediction the PM10 concentrations around the city of Pittsburgh, United States of America. The authors did not model the joint temporal and spatial structure of observations. They first modeled the purely spatial structure of observations for one particular day, then they incorporated the temporal component in their final statistical modeling. Sun et al. [10] developed a spatial forecasting distribution for unmeasured daily log PM10 average concentration given from ten locations in Vancouver, Canada. At each monitoring site, Sun et al. [10] showed that the autoregressive model of order one (AR(1)) described quite well the daily log PM10 average values. The authors did not consider the Bayesian models in their approach. Golam Kibria et al. [11] proposed a multivariate spatial prediction methodology in a Bayesian approach that is suited for spatio-temporal data observed at a small set of ambient monitoring stations at successive time points. They demonstrated the usefulness of their approach by mapping PM2.5 at monitoring sites with different start-up times in the city of Philadelphia, USA. Golam Kibria et al. [11] did not compare the performance of their model with other non-Bayesian models that are commonly used in literature. Shaddick and Wakefield [12] and Sahu and Mardia [13] used short-term spatio-temporal predictive analysis for modeling the PM2.5 and PM10 concentration levels. Zidek et al. [3] presented predictive distributions on nonmonitored PM10 values in Vancouver, Canada. Smith et al. [14] proposed a spatio-temporal model to predict the weekly averages of particulate matter concentrations within three southeastern states in the USA. Sahu et al. [15] modeled the PM2.5 by combining the rural background and the urban areas into one process. Cocchi et al. [1] developed a hierarchical Bayesian model for the daily average PM10 values. Pollice and Lasinio [2] developed a Bayesian-based kriging method for estimating the daily average PM10 concentration levels. Wikle et al. [16] provided an excellent review of classical and Bayesian approaches for analyzing spatio-temporal data. Although Taylor et al. [6], Ahmadi et al. [5], and others have studied the relationship between some environmental conditions and particulate matter levels assessing the air quality based on particulate matter levels in different locations in Qatar, they did not study the spatio-temporal variability of PM10. The main objective of this research is to develop space-time models for daily PM10 air pollution levels in Qatar for the four years, 2016-2019 comparing the hierarchical Bayesian approach with other spatio-temporal recent methods. We develop spatial interpolation and forecasting model using iterative Markov chain Monte Carlo (MCMC) computation setup which is an effective method for modeling a data with large number of missing values [17]. To the best of our knowledge, this will be the first study in Qatar, and we hope that this research will be helpful to protect the environment and public health in Qatar.The rest of this article is organized as follows: Section2 provides a brief review of two-stage hierarchical Bayesian models that have been used for modeling spatio-temporal data. A numerical example is given in Section 3 to demonstrate that the Bayesian approach accurately predicts the daily average PM10 values with a large number of missing values comparing with non-Bayesian models. Finally, a conclusion is given in Section 4.


## 2. Hierarchical Spatio-Temporal Model

When data is collected at different points in space and time, we should use a model that can, simultaneously, describe the dependency structure coming from the three sources of variations: time variation, space variation, and joint variability between time and space. Such a model is called a space-time model (or spatio-temporal model, where spatio refers to space and temporal refers to time).In this article, we develop hierarchical models to predict the daily PM10 concentration levels which vary over time and locations. The PM10 concentration levels do not often follow the normal distribution. Thus, we usually model these values on the square-root scale or we use the log-transformation to stabilize the variance and enforce normality and to stabilize the variance [18]. We consider the square-root scale to alleviate the departure from normality in our research data.Letℓ and t denote the two units of time where ℓ=1,⋯,r represents the longer unit (e.g., year), and t=1,⋯,Tℓ represents the shorter unit (e.g., day). Let Zℓs,t denote the observed value of the PM10 concentration, after any necessary transformation, at a given location s and over a given discrete time t. We assume that the spatial location s is a two-dimensional vector describing the latitude-longitude (or equivalently northing and easting coordinates), and the time unit is typically hour, day, month, or year. We also assume that the Zℓs,t is observed at n monitoring sites denoted by si,i=1,⋯,n and at time points denoted by two indices ℓ and t so that the total number of observations is denoted by N=n∑ℓ=1rTℓ. In this article, we denote all the missing data by z⋆, whereas all the observed data will be denoted by z.The first stage of the hierarchy assumes that the observed valuesZℓt, where Zℓt=Zℓs1,t,⋯,Zℓsn,t′, can be decomposed into a true (latent) spatio-temporal process Yℓt=Yℓs1,t,⋯,Yℓsn,t′ with an error term εℓt=εℓs1,t,⋯,εℓsn,t′. More specifically, the data (or measurement error) model in the first stage of the hierarchy is
(1)Zℓt=Yℓt+εℓt,ℓ=1,⋯,r,t=1,⋯,Tℓ.The error termεℓt is assumed to be a Gaussian white noise process with mean zero and constant variance σε2, which is often called the nugget effect absorbing microscale variability.The second stage assumes that the true processYℓt has a systematic mean μℓt and a spatio-temporal error term. The mean can be modeled based on the past values of the unobserved variable or/and based on some relevant covariates. Typically, Yℓt can be specified in the following formula:
(2)Yℓt=μℓt+ηℓt,where ηℓt=ηℓs1,t,⋯,ηℓsn,t′ is an spatio-temporal residual random intercept assumed to follow N0,C, where C=ση2Hη, ση2 is the site invariant spatial variance, and Hη is the spatial correlation matrix. In this article, we consider four modelings for ηℓt. The first one assumes that the random effects, ηℓsi,t=0, for all locations si and times t. This implies that the model in (1)–(2) will be the simple regression model. The other models for ηℓt are assumed to be Gaussian process, independent over the time, which is specified in Sections 2.1, 2.2, 2.3, and 2.4.


### 2.1. Matérn Spatio-Temporal Covariance Function

For spatio-temporal modeling, we usually assume that the random effects process is a weakly stationary Gaussian with a zero mean and a valid isotropic covariance function. A valid covariance function implied that the covariance matrix is positive definite, and isotropic means that the separation vector between the two locations only depends on the distance and not on the direction. The class of the spatio-temporal covariance functions can be separable, product-sum, metric, and sum-metric [19]. In this paper, we use the separable covariance model which is simply the product of the pure spatial covariance function, Csh, by the pure temporal one, Ctu, given by
(3)Ch,u=CshCtu,where h=s−s′ is the separating spatial distance, and u=t−t′ is the temporal distance for any pair of points s,t×s′,t′ in the spatial and temporal study domain.The Matérn family provides a general choice of covariance functions. For each timet, the Matérn covariance is given by:
(4)Cws,t,ws′,t=ση22κ−1Γκ2κs−s′ϕκKκ2κ2κs−s′ϕϕ,ϕ>0,κ≥1,where Kκ· is the modified Bessel function of second kind of order κ which is a parameter controlling the smoothness of the realized random field [20], Γκ is the standard gamma function, and ϕ is a parameter which controls the decay rate in the correlation as the distance h increases. Popular special cases of Matérn model are (i) κ=0.5 leads to exponential covariance function Ch=ση2exp−ϕh and (ii) Gaussian model, Ch=ση2exp−ϕ2h2 when κ⟶∞.


### 2.2. Gaussian Process Model Specification

Considerμℓt=Xℓtβ in (2), where ℓ=1,⋯,r,t=1,⋯,Tℓ, Xℓt is the n×k+1 design matrix of spatially and temporally varying k-covariates and β=β0,β1,⋯,βk′ is the k+1-dimensional vector of regression coefficients. Thus, the Gaussian process (GP) two-stage model can be written as
(5)Zℓt=Yℓt+εℓt,Yℓt=Xℓtβ+ηℓt.We assume that the random effect processη· is independent from the white noise process ε·. Note that some of the covariates may vary spatially and not temporally or vice versa.One advantage of using Bayesian model is that we can use it to handle any missing data. This can be done by using the Markov chain Monte Carlo, MCMC, computation where the missing data is simulated from theNYℓt,σε2 distribution defined by (5) at each MCMC iteration using Gibbs sampling. The Gibbs sampler requires that the full conditional distributions of the parameters θ=β,σε2,ση2,ϕ,κ are given in a closed form.The logarithm of the joint posterior distribution of the missing data and the parameters in this case is given by:(6)logpθ,Y,z⋆∣z∝−N2logσε2−12σε2∑ℓ=1r∑t=1TℓZℓt−Yℓt′Zℓt−Yℓt−12∑ℓ=1rTℓlogση2Hη−12ση2∑ℓ=1r∑t=1TℓYℓt−Xℓtβ′Hη−1Yℓt−Xℓtβ+logpθ,where pθ is the prior distribution [21] and we refer the readers to Bakar and Sahu ([21], Appendix A) to obtain the Gibbs sampling using the full conditional distributions of θ.


### 2.3. Autoregressive Model Specification

The autoregressive process (AR) model in Equations (1) and (2) can be written in the hierarchical form (e.g., [18] model) as follows:
(7)Zℓt=Yℓt+εℓt,Yℓt=ρYℓt−1+Xℓtβ+ηℓt,where −1<ρ<1 is parameter of the first-order autoregressive model and μℓ=ρYℓt−1+Xℓtβ. Note that if ρ=0, we get the GP model given by (5). Also, note that the autoregressive model requires to specify an independent spatial model with initial values of Yℓ0 for each ℓ=1,⋯,r, with mean μℓ and covariance σℓ2H0 obtained from the Matérn covariance function in (4) with the same set of parameters. In this case, logarithm of the joint posterior distribution of the missing data and the parameters in this case is given by:
(8)logpθ,Y,z⋆∣z∝−N2logσε2−12σε2∑ℓ=1r∑t=1TℓZℓt−Yℓt′Zℓt−Yℓt−12∑ℓ=1rTℓlogση2Hη−12ση2∑ℓ=1r∑t=1TℓYℓt−ρYℓt−1−Xℓtβ′Hη−1Yℓt−ρYℓt−1−Xℓtβ−12∑ℓ=1rlogσℓ2H0−12∑ℓ=1r1σℓ2Yℓ0−μℓ′H0−1Yℓ0−μℓ+logpθ,where pθ is the prior distribution for the parameter θ=β,ρ,σε2,ση2,ϕ,κ,μℓ,σℓ2 ([21], see Appendix B).


### 2.4. Gaussian Predictive Processes with Autoregressive Model Specification

We now introduce the case that the processηℓt in (2) has a nonstationary covariance structure. We follow Sahu and Mukhopadhyay [22] by choosing the integer m out of the sample prediction performance (crossvalidation) to get a smaller number of locations (denoted by knot), Sm⋆=s1⋆,⋯,sm⋆′. Given Sm⋆, we assume that ηℓ⋆=ηℓs1⋆,t,⋯,ηℓsm⋆,t′ is a Gaussian process with zero mean vector and covariance function given in (4). Thus, the nonstationary model is defined after replacing ηℓsi,t in (2) by
(9)η~si,t=Eηsi,t∣ηt⋆,i=1,⋯,n,hence, ηℓt,ηℓt⋆ is a vector of n+m-independent realization at each time point ℓ and t with the same Gaussian process in (4). Define η~ℓ=η~ℓs1,t,⋯,η~ℓsn,t′, we have
(10)η~t=C⋆ϕ,κHη⋆−1ϕ,κηt⋆,where C⋆ϕ,κ is the n×m crosscorrelation matrix between ηt and ηt⋆ having elements C⋆ij=Csi−sj⋆,i=1,⋯,n and j=1,⋯,m, and Hη⋆ϕ,κ is the m×m correlation matrix of ηt⋆ so that Hη⋆ij=Csk⋆−sj⋆;ϕ,κ, for k,j=1,⋯,m. Clearly, the process η~ℓt shows nonstationary structure with variance function that is given by
(11)Varη~t=C⋆ϕ,κHη⋆−1ϕ,κC⋆′ϕ,κ.The advantage of using the nonstationary model in (10) is the flexibility in the η~t surface which is based on m≪n linear functions of ηt⋆. When m is very small compared with n, this will lead to reduce the computational burden, especially for big data which is usually the case of the spatio-temporal data. Moreover, using nonstationary models usually provides more accurate results in prediction the nonstationary PM10 process. We specify the hierarchical Gaussian predictive processes (GPP) as follows:
(12)Zℓt=Yℓt+εℓt,Yℓt=Xℓtβ+η~ℓt,where η~ℓt is given in (10). The process ηℓt⋆, at the Sm⋆ knots, can be modeled according to the autoregressive model
(13)ηℓt⋆=ρηℓt−1⋆+ηℓt,where ηℓt⋆~Nm0,Ση⋆ for ℓ=1,⋯,r,t=1,⋯,Tℓ and Ση⋆=ση2Hη⋆. We assume that the initial conditions ηℓ0⋆ has normal distribution with mean zero and covariance matrix Σ0=σℓ2H0, and both Ση⋆,Σ0 can be obtained from the Matérn covariance function defined in (4), where Ση⋆ is an m×m matrix much lower dimensional than Ση of dimension n×n.In this case, logarithm of the joint posterior distribution of the missing data and the parameters in this case is given by:(14)logpθ,Y,z⋆∣z∝−N2logσε2−12σε2∑ℓ=1r∑t=1TℓZℓt−Yℓt′Zℓt−Yℓt−mN2logση2−N2logHη⋆−m2∑ℓ=1rlogσℓ2−r2logH0−12ση2∑ℓ=1r∑t=1Tℓηℓt⋆−ρηℓt−1⋆′Hη⋆−1ηℓt⋆−ρηℓt−1⋆−12∑ℓ=1r1σℓ2ηℓ0⋆′H0−1ηℓ0⋆+logpθ,where pθ is the prior distribution for the parameter θ=β,ρ,σε2,ση⋆2,ϕ,κ,σℓ2 and the Gibbs sampling procedure of these parameters can be obtained from the full conditional distributions required provided in the Appendix.


## 2.1. Matérn Spatio-Temporal Covariance Function

For spatio-temporal modeling, we usually assume that the random effects process is a weakly stationary Gaussian with a zero mean and a valid isotropic covariance function. A valid covariance function implied that the covariance matrix is positive definite, and isotropic means that the separation vector between the two locations only depends on the distance and not on the direction. The class of the spatio-temporal covariance functions can be separable, product-sum, metric, and sum-metric [19]. In this paper, we use the separable covariance model which is simply the product of the pure spatial covariance function, Csh, by the pure temporal one, Ctu, given by
(3)Ch,u=CshCtu,where h=s−s′ is the separating spatial distance, and u=t−t′ is the temporal distance for any pair of points s,t×s′,t′ in the spatial and temporal study domain.The Matérn family provides a general choice of covariance functions. For each timet, the Matérn covariance is given by:
(4)Cws,t,ws′,t=ση22κ−1Γκ2κs−s′ϕκKκ2κ2κs−s′ϕϕ,ϕ>0,κ≥1,where Kκ· is the modified Bessel function of second kind of order κ which is a parameter controlling the smoothness of the realized random field [20], Γκ is the standard gamma function, and ϕ is a parameter which controls the decay rate in the correlation as the distance h increases. Popular special cases of Matérn model are (i) κ=0.5 leads to exponential covariance function Ch=ση2exp−ϕh and (ii) Gaussian model, Ch=ση2exp−ϕ2h2 when κ⟶∞.


## 2.2. Gaussian Process Model Specification

Considerμℓt=Xℓtβ in (2), where ℓ=1,⋯,r,t=1,⋯,Tℓ, Xℓt is the n×k+1 design matrix of spatially and temporally varying k-covariates and β=β0,β1,⋯,βk′ is the k+1-dimensional vector of regression coefficients. Thus, the Gaussian process (GP) two-stage model can be written as
(5)Zℓt=Yℓt+εℓt,Yℓt=Xℓtβ+ηℓt.We assume that the random effect processη· is independent from the white noise process ε·. Note that some of the covariates may vary spatially and not temporally or vice versa.One advantage of using Bayesian model is that we can use it to handle any missing data. This can be done by using the Markov chain Monte Carlo, MCMC, computation where the missing data is simulated from theNYℓt,σε2 distribution defined by (5) at each MCMC iteration using Gibbs sampling. The Gibbs sampler requires that the full conditional distributions of the parameters θ=β,σε2,ση2,ϕ,κ are given in a closed form.The logarithm of the joint posterior distribution of the missing data and the parameters in this case is given by:(6)logpθ,Y,z⋆∣z∝−N2logσε2−12σε2∑ℓ=1r∑t=1TℓZℓt−Yℓt′Zℓt−Yℓt−12∑ℓ=1rTℓlogση2Hη−12ση2∑ℓ=1r∑t=1TℓYℓt−Xℓtβ′Hη−1Yℓt−Xℓtβ+logpθ,where pθ is the prior distribution [21] and we refer the readers to Bakar and Sahu ([21], Appendix A) to obtain the Gibbs sampling using the full conditional distributions of θ.


## 2.3. Autoregressive Model Specification

The autoregressive process (AR) model in Equations (1) and (2) can be written in the hierarchical form (e.g., [18] model) as follows:
(7)Zℓt=Yℓt+εℓt,Yℓt=ρYℓt−1+Xℓtβ+ηℓt,where −1<ρ<1 is parameter of the first-order autoregressive model and μℓ=ρYℓt−1+Xℓtβ. Note that if ρ=0, we get the GP model given by (5). Also, note that the autoregressive model requires to specify an independent spatial model with initial values of Yℓ0 for each ℓ=1,⋯,r, with mean μℓ and covariance σℓ2H0 obtained from the Matérn covariance function in (4) with the same set of parameters. In this case, logarithm of the joint posterior distribution of the missing data and the parameters in this case is given by:
(8)logpθ,Y,z⋆∣z∝−N2logσε2−12σε2∑ℓ=1r∑t=1TℓZℓt−Yℓt′Zℓt−Yℓt−12∑ℓ=1rTℓlogση2Hη−12ση2∑ℓ=1r∑t=1TℓYℓt−ρYℓt−1−Xℓtβ′Hη−1Yℓt−ρYℓt−1−Xℓtβ−12∑ℓ=1rlogσℓ2H0−12∑ℓ=1r1σℓ2Yℓ0−μℓ′H0−1Yℓ0−μℓ+logpθ,where pθ is the prior distribution for the parameter θ=β,ρ,σε2,ση2,ϕ,κ,μℓ,σℓ2 ([21], see Appendix B).


## 2.4. Gaussian Predictive Processes with Autoregressive Model Specification

We now introduce the case that the processηℓt in (2) has a nonstationary covariance structure. We follow Sahu and Mukhopadhyay [22] by choosing the integer m out of the sample prediction performance (crossvalidation) to get a smaller number of locations (denoted by knot), Sm⋆=s1⋆,⋯,sm⋆′. Given Sm⋆, we assume that ηℓ⋆=ηℓs1⋆,t,⋯,ηℓsm⋆,t′ is a Gaussian process with zero mean vector and covariance function given in (4). Thus, the nonstationary model is defined after replacing ηℓsi,t in (2) by
(9)η~si,t=Eηsi,t∣ηt⋆,i=1,⋯,n,hence, ηℓt,ηℓt⋆ is a vector of n+m-independent realization at each time point ℓ and t with the same Gaussian process in (4). Define η~ℓ=η~ℓs1,t,⋯,η~ℓsn,t′, we have
(10)η~t=C⋆ϕ,κHη⋆−1ϕ,κηt⋆,where C⋆ϕ,κ is the n×m crosscorrelation matrix between ηt and ηt⋆ having elements C⋆ij=Csi−sj⋆,i=1,⋯,n and j=1,⋯,m, and Hη⋆ϕ,κ is the m×m correlation matrix of ηt⋆ so that Hη⋆ij=Csk⋆−sj⋆;ϕ,κ, for k,j=1,⋯,m. Clearly, the process η~ℓt shows nonstationary structure with variance function that is given by
(11)Varη~t=C⋆ϕ,κHη⋆−1ϕ,κC⋆′ϕ,κ.The advantage of using the nonstationary model in (10) is the flexibility in the η~t surface which is based on m≪n linear functions of ηt⋆. When m is very small compared with n, this will lead to reduce the computational burden, especially for big data which is usually the case of the spatio-temporal data. Moreover, using nonstationary models usually provides more accurate results in prediction the nonstationary PM10 process. We specify the hierarchical Gaussian predictive processes (GPP) as follows:
(12)Zℓt=Yℓt+εℓt,Yℓt=Xℓtβ+η~ℓt,where η~ℓt is given in (10). The process ηℓt⋆, at the Sm⋆ knots, can be modeled according to the autoregressive model
(13)ηℓt⋆=ρηℓt−1⋆+ηℓt,where ηℓt⋆~Nm0,Ση⋆ for ℓ=1,⋯,r,t=1,⋯,Tℓ and Ση⋆=ση2Hη⋆. We assume that the initial conditions ηℓ0⋆ has normal distribution with mean zero and covariance matrix Σ0=σℓ2H0, and both Ση⋆,Σ0 can be obtained from the Matérn covariance function defined in (4), where Ση⋆ is an m×m matrix much lower dimensional than Ση of dimension n×n.In this case, logarithm of the joint posterior distribution of the missing data and the parameters in this case is given by:(14)logpθ,Y,z⋆∣z∝−N2logσε2−12σε2∑ℓ=1r∑t=1TℓZℓt−Yℓt′Zℓt−Yℓt−mN2logση2−N2logHη⋆−m2∑ℓ=1rlogσℓ2−r2logH0−12ση2∑ℓ=1r∑t=1Tℓηℓt⋆−ρηℓt−1⋆′Hη⋆−1ηℓt⋆−ρηℓt−1⋆−12∑ℓ=1r1σℓ2ηℓ0⋆′H0−1ηℓ0⋆+logpθ,where pθ is the prior distribution for the parameter θ=β,ρ,σε2,ση⋆2,ϕ,κ,σℓ2 and the Gibbs sampling procedure of these parameters can be obtained from the full conditional distributions required provided in the Appendix.


## 3. Numerical Example

### 3.1. Data Description

The data used in this article is obtained from three different sources. The first one was collected by the air pointer and the meteorological station and is managed by the Environmental Science Center at Qatar University over the years 2016-2019. It has hourly pollutant including PM10 measured in μg/m3, temperature, Temp, measured in degree Celsius, and relative humidity, RH, with several missing values.First, we sort the data into an ascending order by date, and then, we impute the missing values of the PM10 by averaging the two nonmissing values before and after this missing value. When two or more successive missing values exist, we impute them by the corresponding monthly average. We do the same for the hourly missing values of temperature and relative humidity. After that, we aggregate the hourly data into daily data for which we fit our models. The amended data from this resource is irregular time series which started on the fourth of January 2016 and ended on the twenty-eighth of August 2019. Table 1 provides a summary statistics of this data. Results suggest that the PM10 levels increased during the time, and the majority of missing observations (with more than 52%) have been occurred during the year 2018.Table 1
Summary statistics of the daily PM10 pollutant measured in μg/m3, temperature (denoted by Temp) measured in degree Celsius, and relative humidity (RH) obtained by air pointer and the meteorological station at Qatar University site. n and SD stand for sample size and standard deviation, respectively. The % n.a. denotes the percentage of missing values in the data.
Yearn% n.a.MinMaxMedianMeanSDPM1020163484.9218.35659.4976.95102.6389.2320172572.2866.221379.13103.15154.60157.84201808052.6625.89551.5163.0380.6768.90201920115.9037.871542.5772.29117.46157.24Overall88610.9618.351542.5783.35119.09104.68Temp20163630.828.7639.0828.3227.856.89201722115.9710.8641.6733.1029.777.67201813420.7118.5737.6328.6628.624.32201912647.2822.2840.2833.7532.375.12Overall8448.708.7641.6731.4029.146.30RH20163630.8215.1392.7756.3453.6718.31201722115.9718.1686.9452.1250.5716.74201813420.7123.4089.0557.7455.1916.86201912647.2820.8380.2946.4048.3915.99Overall8448.7015.1392.7750.8152.3116.66The second source of the data was the regular monthly observations obtained from nine meteorological stations over the years 2016-2019. The total number of observations is432 where the stations are located in the following sites: Qatar University, Abu Samra, Al Khor, Al Ruwais, Al Wakrah, Doha Airport, Dukhan, Mukenis-Al Karanaah, and Umm Said. Figure 1 shows the map locations of 9 meteorological sites. Although this data provides the monthly average values of temperature in degree Celsius and relative humidity (RH), it does not include measures for any type of pollutant. Table 2 summarizes the descriptive statistics of this data. We use this data to interpolate the daily temperature and relative humidity by disaggregating technique at each location. We also use our own portable devices to collect daily PM10 values, as a third source of data, over the period of time from November 5th to December 31th of 2020 close to these sites. We use these devices to collect and simulate more daily PM10 concentrations over the aforementioned locations. In particular, we recorded the average of the PM10 taken randomly in the morning and evening of each day over this period. Then, we use this collected data to simulate a new data in each location. Finally, we merge the simulated data with the previous one and obtained the data that we use in this article.Figure 1
A map of the state of Qatar showing locations of the 9 monitoring sites.Table 2
Summary statistics of the available temperature in Celsius, and relative humidity measured at 9 monitoring sites in Qatar; SD stands for standard deviation.
YearMinMaxMedianMeanSDTemp201616.5038.0027.6527.486.57201716.3037.9029.7027.906.99201816.4037.8028.9028.216.44201918.8037.6028.5028.006.42Overall16.3038.0028.3027.906.59RH201625.0074.0057.5055.2811.74201723.0074.0057.0053.9211.27201824.0073.0054.0052.7211.70201927.0069.0056.0054.1610.33Overall23.0074.0056.0054.0211.27The final data has 1037 observations (see Figure2) which has three variables (PM10, temperature, and relative humidity) measured over irregular time over the years 2016-2020 in nine spatial locations.Figure 2
Daily PM10 concentrations observed in nine locations over the years 2016-2021.The top panel in Figure3 shows the distribution of the daily average PM10 concentrations over the years 2016-2020, whereas the bottom one in the same figure shows these averages by the nine sites. Clearly that the distribution of the PM10 concentrations, in all locations, are right skewed with very high extreme values. Thus, to stabilize the variance and reduce the departure from normality, we transform the original scale to square root scale.Figure 3
Boxplots of daily average of PM10 concentrations over the years 2016-2019 and by 9 locations in Qatar showing the quartiles; the whiskers show the farthest observation from the median.


### 3.2. Parameter Estimates, Model Validations, Prediction, and Comparison Results

The main objective of this research is to develop a hierarchical Bayesian model that can be used to select and validate the best model which fits the daily average of the PM10 air pollutant levels over different locations in Qatar. We consider the two covariates, temperature and relative humidity, for spatial interpolation of the PM10 concentration at a new location and any time.First, we use the MCMC algorithm with Gibbs sampler, based on 5,000 iterations discarding the first 1,000 values, to impute the missing dataz⋆ as explained in Appendix, and then, we utilize this data to predict the values of Zs′,t at new location s′∉s1,⋯,sn. The posterior predictive distribution is
(15)pzs′,t∣z=∫pzs′,t∣Sm⋆,η⋆,θ,zpSm⋆,η⋆,θ∣zdSm⋆dη⋆dθ,where θ=ρ,σε2,ση2,ϕ′ are the model parameters and η⋆=η1⋆,⋯,ηn⋆′ (see [17]). We fit the Bayesian GP, AR, and GPP models described in Section 2 using the spTimer package [21, 23] and compare these models with the non-Bayesian generalized additive model (GAM) using the R package mgcv [24–26]. We use the crossvalidation method to evaluate the predictive performance of these models. Here, data at locations s1,⋯,sm where m<n are used to fit the model while data at other sites sm+1,⋯,sn are used to assess the model using the mean squared error (MSE), mean absolute error (MAE), mean absolute prediction error (MAPE), relative bias (rBIAS), and relative mean separation (rMSEP) defined, respectively, as follows:
(16)MSE=1N∑i=1Nzi−z^i2,MAE=1N∑i=1Nzi−z^i,MAPE=1N∑i=1Nzi−z^izi,rBIAS=1N∑i=1Nzi−z^iz¯,rMSEP=1N∑i=1Nzi−z^i2z^¯−zi,where N is the total number of nonmissing observations, z^i is the posterior predictive value of zi, and z¯,z^¯−zi are the arithmetic means of zi,z^i for i=1,⋯,N.Specifically, we use the data from Abu-Samra and Doha airport stations with2×1037=2074 observations for validation purposes. The data from the remaining six stations with 7×1037=7259 total number of observations are used for model fitting (see Figure 1). The parameter estimates (posterior) for the Bayesian spatio-temporal models are given in Table 3. The 95% credible intervals for the parameters suggest that most of regression parameters for the AR and GP models are statistically significant, whereas the GPP model suggest that the variables are significant. In all models, the nugget effect σε2 is small. On the other hand, the spatial decay parameter ϕ=0.0033 for GP model and ϕ=0.001 for GPP model suggesting that the effective ranges will approximately be 909 and 3000 kilometers, respectively. These two values are unusual associated with very large spatial variance values (ση2) suggesting that the PM10 concentrations in Qatar are not significantly different over the space.Table 3
Posterior estimates and 95% credible interval for the parameters of the AR, GP, and GPP models fitted for PM10 concentration levels. SD stands for standard deviation.
ModelsParametersMeanMedianSD95% credible intervalARβ01.35211.35190.0206(1.3115, 1.3934)Temp0.00010.00010.0004(-0.0007, 0.0009)RH0.00000.00000.0002(-0.0004, 0.0003)ρ0.01600.01450.0267(-0.0328, 0.0734)σε20.02090.02080.0007(0.0199, 0.0226)ση20.01240.01230.0009(0.0107, 0.0144)ϕ0.01680.01680.0000(0.0168, 0.0168)GPβ012.130312.43781.1060(9.7744, 13.5797)Temp-0.0173-0.02670.0334(-0.0590, 0.0538)RH-0.0153-0.01630.0041(-0.0209, -0.0062)σε20.00520.00520.0001(0.0050, 0.0054)ση215.226212.83937.9727(6.6616, 35.1414)ϕ0.00320.00270.0018(0.0011, 0.0069)GPPβ00.47920.47840.0692(0.3425, 0.6131)Temp-0.0028-0.00280.0013(-0.0053, -0.0003)RH-0.0013-0.00130.0003(-0.0018, -0.0007)ρ0.97040.97040.0025(0.9654, 0.9754)σε20.00520.00520.0001(0.0050, 0.0054)ση21.74501.74530.0256(1.6955, 1.7954)ϕ0.00100.00100.0000(0.0010, 0.0010)After fitting the three models, we perform the crossvalidation and select the best one. Table4 summarized the validation statistics for the GAM, GP, and AR models. Clearly, the MSE, MAE, MAPE, rBIAS, and rMSEP are smaller for the Bayesian spatio-temporal methods providing better predictive performance compared to the non-Bayesian additive model. For example, the Bayesian AP, GP, and GPP models reduced the MSE by about 63.7%, 66.4%, and 79.6%, respectively, compared with the non-Bayesian GAM model. We conclude from this table that the Bayesian spatio-temporal GPP model is the best model that can be used to predict the PM10 concentrations.Table 4
Validation statistics for the GP, AR, and GAM models.
ModelsMSEMAEMAPErBIASrMSEPNon-BayesianGAM991.0479.95122.430.5570.773AR360.1536.3746.530.3370.342Spatio-temporalGP332.8236.6148.3790.3390.316GPP202.0416.6128.3510.1030.096Figure4 shows the time series plot of the true and fitted values for Qatar University location using the GPP model. We clearly see that the fitted values are very close to the true values. To further demonstrate the usefulness of the spatio-temporal GPP model for predicting the linear trend surfaces of the PM10 values with their standard values, we illustrate a prediction map for the 29th of March 2018 and 15 May of 2019 over a one-kilometer square grid (see Figure 5). Obviously, the graphs show that the GPP model correctly represented the PM10 concentration level over the space and time.Figure 4
Time series plots of the true (blue) and fitted values (green) at one randomly chosen location (Qatar University).Figure 5
Spatially interpolated plots of the daily average of PM10 concentration levels (a) and their standard deviations (b) obtained from the AR model for the 29th of March 2018 (top graphs) and 15 May of 2019 (bottom graphs). Actual values and associated locations are superimposed in plots, respectively.
(a)(b)


## 3.1. Data Description

The data used in this article is obtained from three different sources. The first one was collected by the air pointer and the meteorological station and is managed by the Environmental Science Center at Qatar University over the years 2016-2019. It has hourly pollutant including PM10 measured in μg/m3, temperature, Temp, measured in degree Celsius, and relative humidity, RH, with several missing values.First, we sort the data into an ascending order by date, and then, we impute the missing values of the PM10 by averaging the two nonmissing values before and after this missing value. When two or more successive missing values exist, we impute them by the corresponding monthly average. We do the same for the hourly missing values of temperature and relative humidity. After that, we aggregate the hourly data into daily data for which we fit our models. The amended data from this resource is irregular time series which started on the fourth of January 2016 and ended on the twenty-eighth of August 2019. Table 1 provides a summary statistics of this data. Results suggest that the PM10 levels increased during the time, and the majority of missing observations (with more than 52%) have been occurred during the year 2018.Table 1
Summary statistics of the daily PM10 pollutant measured in μg/m3, temperature (denoted by Temp) measured in degree Celsius, and relative humidity (RH) obtained by air pointer and the meteorological station at Qatar University site. n and SD stand for sample size and standard deviation, respectively. The % n.a. denotes the percentage of missing values in the data.
Yearn% n.a.MinMaxMedianMeanSDPM1020163484.9218.35659.4976.95102.6389.2320172572.2866.221379.13103.15154.60157.84201808052.6625.89551.5163.0380.6768.90201920115.9037.871542.5772.29117.46157.24Overall88610.9618.351542.5783.35119.09104.68Temp20163630.828.7639.0828.3227.856.89201722115.9710.8641.6733.1029.777.67201813420.7118.5737.6328.6628.624.32201912647.2822.2840.2833.7532.375.12Overall8448.708.7641.6731.4029.146.30RH20163630.8215.1392.7756.3453.6718.31201722115.9718.1686.9452.1250.5716.74201813420.7123.4089.0557.7455.1916.86201912647.2820.8380.2946.4048.3915.99Overall8448.7015.1392.7750.8152.3116.66The second source of the data was the regular monthly observations obtained from nine meteorological stations over the years 2016-2019. The total number of observations is432 where the stations are located in the following sites: Qatar University, Abu Samra, Al Khor, Al Ruwais, Al Wakrah, Doha Airport, Dukhan, Mukenis-Al Karanaah, and Umm Said. Figure 1 shows the map locations of 9 meteorological sites. Although this data provides the monthly average values of temperature in degree Celsius and relative humidity (RH), it does not include measures for any type of pollutant. Table 2 summarizes the descriptive statistics of this data. We use this data to interpolate the daily temperature and relative humidity by disaggregating technique at each location. We also use our own portable devices to collect daily PM10 values, as a third source of data, over the period of time from November 5th to December 31th of 2020 close to these sites. We use these devices to collect and simulate more daily PM10 concentrations over the aforementioned locations. In particular, we recorded the average of the PM10 taken randomly in the morning and evening of each day over this period. Then, we use this collected data to simulate a new data in each location. Finally, we merge the simulated data with the previous one and obtained the data that we use in this article.Figure 1
A map of the state of Qatar showing locations of the 9 monitoring sites.Table 2
Summary statistics of the available temperature in Celsius, and relative humidity measured at 9 monitoring sites in Qatar; SD stands for standard deviation.
YearMinMaxMedianMeanSDTemp201616.5038.0027.6527.486.57201716.3037.9029.7027.906.99201816.4037.8028.9028.216.44201918.8037.6028.5028.006.42Overall16.3038.0028.3027.906.59RH201625.0074.0057.5055.2811.74201723.0074.0057.0053.9211.27201824.0073.0054.0052.7211.70201927.0069.0056.0054.1610.33Overall23.0074.0056.0054.0211.27The final data has 1037 observations (see Figure2) which has three variables (PM10, temperature, and relative humidity) measured over irregular time over the years 2016-2020 in nine spatial locations.Figure 2
Daily PM10 concentrations observed in nine locations over the years 2016-2021.The top panel in Figure3 shows the distribution of the daily average PM10 concentrations over the years 2016-2020, whereas the bottom one in the same figure shows these averages by the nine sites. Clearly that the distribution of the PM10 concentrations, in all locations, are right skewed with very high extreme values. Thus, to stabilize the variance and reduce the departure from normality, we transform the original scale to square root scale.Figure 3
Boxplots of daily average of PM10 concentrations over the years 2016-2019 and by 9 locations in Qatar showing the quartiles; the whiskers show the farthest observation from the median.


## 3.2. Parameter Estimates, Model Validations, Prediction, and Comparison Results

The main objective of this research is to develop a hierarchical Bayesian model that can be used to select and validate the best model which fits the daily average of the PM10 air pollutant levels over different locations in Qatar. We consider the two covariates, temperature and relative humidity, for spatial interpolation of the PM10 concentration at a new location and any time.First, we use the MCMC algorithm with Gibbs sampler, based on 5,000 iterations discarding the first 1,000 values, to impute the missing dataz⋆ as explained in Appendix, and then, we utilize this data to predict the values of Zs′,t at new location s′∉s1,⋯,sn. The posterior predictive distribution is
(15)pzs′,t∣z=∫pzs′,t∣Sm⋆,η⋆,θ,zpSm⋆,η⋆,θ∣zdSm⋆dη⋆dθ,where θ=ρ,σε2,ση2,ϕ′ are the model parameters and η⋆=η1⋆,⋯,ηn⋆′ (see [17]). We fit the Bayesian GP, AR, and GPP models described in Section 2 using the spTimer package [21, 23] and compare these models with the non-Bayesian generalized additive model (GAM) using the R package mgcv [24–26]. We use the crossvalidation method to evaluate the predictive performance of these models. Here, data at locations s1,⋯,sm where m<n are used to fit the model while data at other sites sm+1,⋯,sn are used to assess the model using the mean squared error (MSE), mean absolute error (MAE), mean absolute prediction error (MAPE), relative bias (rBIAS), and relative mean separation (rMSEP) defined, respectively, as follows:
(16)MSE=1N∑i=1Nzi−z^i2,MAE=1N∑i=1Nzi−z^i,MAPE=1N∑i=1Nzi−z^izi,rBIAS=1N∑i=1Nzi−z^iz¯,rMSEP=1N∑i=1Nzi−z^i2z^¯−zi,where N is the total number of nonmissing observations, z^i is the posterior predictive value of zi, and z¯,z^¯−zi are the arithmetic means of zi,z^i for i=1,⋯,N.Specifically, we use the data from Abu-Samra and Doha airport stations with2×1037=2074 observations for validation purposes. The data from the remaining six stations with 7×1037=7259 total number of observations are used for model fitting (see Figure 1). The parameter estimates (posterior) for the Bayesian spatio-temporal models are given in Table 3. The 95% credible intervals for the parameters suggest that most of regression parameters for the AR and GP models are statistically significant, whereas the GPP model suggest that the variables are significant. In all models, the nugget effect σε2 is small. On the other hand, the spatial decay parameter ϕ=0.0033 for GP model and ϕ=0.001 for GPP model suggesting that the effective ranges will approximately be 909 and 3000 kilometers, respectively. These two values are unusual associated with very large spatial variance values (ση2) suggesting that the PM10 concentrations in Qatar are not significantly different over the space.Table 3
Posterior estimates and 95% credible interval for the parameters of the AR, GP, and GPP models fitted for PM10 concentration levels. SD stands for standard deviation.
ModelsParametersMeanMedianSD95% credible intervalARβ01.35211.35190.0206(1.3115, 1.3934)Temp0.00010.00010.0004(-0.0007, 0.0009)RH0.00000.00000.0002(-0.0004, 0.0003)ρ0.01600.01450.0267(-0.0328, 0.0734)σε20.02090.02080.0007(0.0199, 0.0226)ση20.01240.01230.0009(0.0107, 0.0144)ϕ0.01680.01680.0000(0.0168, 0.0168)GPβ012.130312.43781.1060(9.7744, 13.5797)Temp-0.0173-0.02670.0334(-0.0590, 0.0538)RH-0.0153-0.01630.0041(-0.0209, -0.0062)σε20.00520.00520.0001(0.0050, 0.0054)ση215.226212.83937.9727(6.6616, 35.1414)ϕ0.00320.00270.0018(0.0011, 0.0069)GPPβ00.47920.47840.0692(0.3425, 0.6131)Temp-0.0028-0.00280.0013(-0.0053, -0.0003)RH-0.0013-0.00130.0003(-0.0018, -0.0007)ρ0.97040.97040.0025(0.9654, 0.9754)σε20.00520.00520.0001(0.0050, 0.0054)ση21.74501.74530.0256(1.6955, 1.7954)ϕ0.00100.00100.0000(0.0010, 0.0010)After fitting the three models, we perform the crossvalidation and select the best one. Table4 summarized the validation statistics for the GAM, GP, and AR models. Clearly, the MSE, MAE, MAPE, rBIAS, and rMSEP are smaller for the Bayesian spatio-temporal methods providing better predictive performance compared to the non-Bayesian additive model. For example, the Bayesian AP, GP, and GPP models reduced the MSE by about 63.7%, 66.4%, and 79.6%, respectively, compared with the non-Bayesian GAM model. We conclude from this table that the Bayesian spatio-temporal GPP model is the best model that can be used to predict the PM10 concentrations.Table 4
Validation statistics for the GP, AR, and GAM models.
ModelsMSEMAEMAPErBIASrMSEPNon-BayesianGAM991.0479.95122.430.5570.773AR360.1536.3746.530.3370.342Spatio-temporalGP332.8236.6148.3790.3390.316GPP202.0416.6128.3510.1030.096Figure4 shows the time series plot of the true and fitted values for Qatar University location using the GPP model. We clearly see that the fitted values are very close to the true values. To further demonstrate the usefulness of the spatio-temporal GPP model for predicting the linear trend surfaces of the PM10 values with their standard values, we illustrate a prediction map for the 29th of March 2018 and 15 May of 2019 over a one-kilometer square grid (see Figure 5). Obviously, the graphs show that the GPP model correctly represented the PM10 concentration level over the space and time.Figure 4
Time series plots of the true (blue) and fitted values (green) at one randomly chosen location (Qatar University).Figure 5
Spatially interpolated plots of the daily average of PM10 concentration levels (a) and their standard deviations (b) obtained from the AR model for the 29th of March 2018 (top graphs) and 15 May of 2019 (bottom graphs). Actual values and associated locations are superimposed in plots, respectively.
(a)(b)


## 4. Conclusion

The potential of applying nonstationary hierarchical Bayesian spatio-temporal models in PM10 prediction with a large number of missing values is presented in this paper. The predicting model is developed by comparing the Gaussian predictive processes (GPP) with Gaussian processes (GP), autoregressive (AR) Bayesian models, and non-Bayesian generalized additive model (GAM) models using the datasets from the state of Qatar. The numerical results show that the GPP model outperforms other alternative models providing forecasting with good accuracy and interpretability. We applied the disaggregating technique and simulated a daily spatio-temporal PM10 data using the available and collected data. Then, we used the Markov chain Monte Carlo with Gibbs sampler to impute the missing data in real collected data. We believe that our statistical data analysis approach will give similar results for future available real data. In many applications, the support vector machine (SVM) algorithm has shown a superior forecasting performance compared with several evolutionary algorithms. This could be an interesting possible extension to this article where further research can be done by comparing the performance of the SVM algorithms with Bayesian models in predicting the daily average PM10 concentration levels.



---

*Source: 8003952-2021-09-13.xml*