Bayesian Hierarchical Modeling for Predicting Spatially Correlated Curves in Irregular Domains: A Case Study on PM₁₀ Pollution

Polynomials are an attractive class of functions for various scientific and engineering computations. They are concisely represented by coefficients on a suitable basis and are amenable to efficient evaluation by simple algorithms. The set of polynomials is closed under the arithmetic operations of addition, subtraction, multiplication, differentiation, integration, and composition. The approximative capabilities of polynomials are also of great practical interest in applications. Perhaps the most fundamental result in this context is the theorem of Weierstrass, which is stated below (Davis 1975):

In other words, it is possible to uniformly approximate any continuous function $g$, defined on a polynomial’s closed interval $[a,b]$. An elegant constructive proof of this theorem was published in 1912, in which Bernstein’s polynomial basis was first introduced; for more details, see Bernstein (1912) and Lorentz (2012).

An illustration of the Bernstein basis can be seen in Figure 1. The vector of basis $\mathbf{b}^{p}(t)=\left(b_{\,0}^{\,p}(t),\ldots,b_{\,p}^{\,p}(t)\right)$ has a weight role, which varies with $t$, as shown in both panels of Figure 1. Thus, the approximation of the target function $g(\cdot)$ is weighted by $p+1$ values coming from the basis vector. It is clear that when $p=9$, more information is available to weigh the function values, resulting in a more accurate approximation.

A spatial functional process is described as $\left\{X_{s}:s\in\mathcal{D}\subset\mathbb{R}^{d}\right\}$, where $X_{s}$ are the functional random variables, located at the site $s$ in the $d$-dimensional Euclidean space, usually $d=2$ or $3$. Each $X_{s}$ is defined on the interval $T=[a,b]\subseteq\mathbb{R}$ and is assumed to belong to a Hilbert space of square-integrable functions, that is, in

with the inner product $\left\langle X_{s},X_{s^{\ast}}\right\rangle=\int_{T}X_{s}(t)X_{s^{\ast}}(t)dt$.

Let $\left\{X_{s_{1}}(t),\ldots,X_{s_{m}}(t)\right\}$ be a sample of $m$ non-independent curves that are indexed by a domain $t\in T$ at each spatial location $s_{j}$, $j=1,\ldots,m$. The functions are observed through a finite set of discrete points, denoted as $t_{ij},\,i=1,\ldots,n_{j}$. It’s essential to note that these measurements are frequently affected by noise. Considering a specific location represented by $s_{j}$, we assume that the observed functions are generated according to the following model

where $\delta_{ij}$ represents a random effect, and $\epsilon_{s_{j}}(t_{ij})$ are independent random errors associated with each fixed point $t_{ij}$. These errors have a mean of zero, i.e., $E(\epsilon_{s_{j}}(t_{ij}))=0$, and a common unknown variance, denoted as $Var(\epsilon_{s_{j}}(t_{ij}))=\tau^{2}$. Henceforth, we assume that each realization $X_{s_{j}}(t_{ij})$ of an underlying random function $X_{s}(t)$ can be adequately represented by a finite set of BP basis functions denoted as $b_{\,0}^{\,p}(t_{ij}),\ldots,b_{\,p}^{\,p}(t_{ij})$ to obtain a reasonable approximation. Consequently, each curve can be expanded into this basis as follows

for $p\in\mathbb{N}$. The basis coefficients, denoted as $\theta_{r,j}$, are crucial in constructing the correlation between curves. Drawing inspiration from the works of Liu et al. (2017) and Aguilera-Morillo et al. (2017), it is assumed that $\theta_{r,j}$ is associated across different locations $j$ for each $r$, but not across different basis function indices $r$. To fully define the spatial correlation structure among $\theta_{r,j}$’s, we introduce the possibility of a non-zero covariance. Specifically, we express this covariance as $Cov(\theta_{r,j},\theta_{r,j^{*}})=C(s_{j},s_{j^{*}})=C(h)$, where $j\neq j^{*}$. Notably, the covariance function is contingent on the locations $s_{j}$ and $s_{j}$ solely through the Euclidean distance $h=\parallel s_{j}-s_{j^{*}}\parallel\in\mathbb{R}$. Taking Equation (6) into account, we can express the model described in Equation (5) as follows

In the context of SFD, a common practice involves linking discrete realizations of each function with finite points $t$ that are regularly spaced along the functional domain. However, a fundamental and distinct aspect of the methodology proposed in this paper is the use of non-equidistant distances between measurement points. In other words, given a sequence of increasing consecutive points $t_{ij}$, denoted by $d_{ij}=(t_{ij}-t_{(i-1)j})\neq(t_{i^{*}j}-t_{(i^{*}-1)j})=d_{i^{*}j}$ for some $i\neq i^{*}\in\left\{2,\ldots,n_{j}\right\}$ are not necessarily equal. In equation (7), we use the locations’ geographic coordinates $s_{j}$ to introduce spatial dependence between the observed functional trajectories. This approach takes into account both the spatial relationships and an autoregressive random effect $\delta_{ij}$, which employs the distances $d_{ij}$ to measure the similarity between discrete measurements closely located within each function $Y_{s_{j}}(t)$. We include this structure to account for the non-equidistant spacing of observations throughout the functional domain. The mathematical structure of $\delta_{j}=\left(\delta_{1j},\ldots,\delta_{n_{i}j}\right)^{\top}$ is defined below.

We introduce a sequence of uncorrelated and identically distributed Gaussian random variables in the given specification, denoted by $\epsilon_{\delta_{ij}}$. These variables have a mean of zero and a variance of $\nu^{2}$, for $i=1,\ldots,n_{j}$ and $j=1,\ldots,m$. Regarding $\phi_{ij}$, it has the structure of an exponential correlation function, where $d_{ij}$ represents the interval between two observations and $\eta$ is the parameter that controls the correlation’s decay rate. A higher value of $\eta$ results in a faster decay, implying a lower correlation between farther-apart observations. Note that an autoregressive structure is defined to associate the random effects in $\delta$. In other words, it ensures that the level of relationship between $\delta_{ij}$ and $\delta_{(i-1)j}$, concerning the positions $t_{ij}$ and $t_{(i-1)j}$ of the functional domain, is controlled by the coefficient $\phi$. Assume $d_{0j}=0$ (so $\phi_{1j}=1$) and $\delta_{0j}=0$; it implies that $\delta_{1j}\sim N(0,\nu^{2})$.

The choice of a Bayesian hierarchical model is crucial when dealing with complex data sets that exhibit both within-group and between-group variation. By structuring the model hierarchically, we can capture these different sources of variability and estimate parameters at each level, providing more robust and interpretable results.

Let’s consider the $i$th discrete observation for the $j$th curve, denoted as $Y_{s_{j}}(t_{ij})$, along with a set of BP basis functions $\left\{b_{\,0}^{\,p}(t_{ij}),\ldots,b_{\,p}^{\,p}(t_{ij})\right\}$. The model incorporates parameters $\theta_{rj}$, $\delta_{ij}$, and $\tau^{2}$ (see Subsection 2.2). We can construct the following Bayesian hierarchical model: sampling for $i=1,\ldots,n_{j}$ and $j=1,\ldots,m$,

In this research, we are utilizing the Gaussian covariance function as described by Banerjee et al. (2014), represented as $C(h)=\kappa^{2}\exp(-(\varphi\,h)^{2})$, to model an isotropic spatial process within the $m$-order covariance matrix, denoted as $\bm{\Sigma}_{m}$. In this context, $\kappa^{2}$ signifies spatial variation, while $\varphi$ is the spatial decay parameter. The symbol $h$ represents the Euclidean distance between two locations $s_{j}$ and $s_{j^{*}}$. It is crucial to emphasize that when the locations are close in $\mathbb{R}^{d}$, the basis coefficients tend to show similarity, resulting in curves with similar shapes, where $C(h)\approx\kappa^{2}$. Conversely, as the distance between these locations increases, the dissimilarity in the curves also increases, leading to $C(h)\approx 0$. For the spatial parameters $\kappa^{2}$ and $\varphi$, we are using prior Inverse Gamma distributions (IG). This choice is common in Bayesian hierarchical models due to the IG conjugacy property, simplifying subsequent calculations. We use Gaussian priors for the coefficients $\bm{\theta}_{r}$ towards a mean value $\mu_{\theta_{r}}$ based on our prior knowledge or assumptions about their distribution. By representing $\bm{\theta}_{r}$ with a multivariate normal distribution, we can account for potential correlations between the coefficients across different locations, capturing spatial dependencies in the model. In this context, $\bm{\theta}_{r}=\left(\theta_{r1},\ldots,\theta_{rm}\right)^{\top}$ represents a vector containing the values of the $r$-th coefficient at the $m$ locations, and $\mu_{\theta_{r}}$ is the mean of the $r$-th basis coefficient. Additionally, utilizing a Gaussian prior for $\mu_{\theta_{r}}$ introduces regularization to the mean of the coefficients, enabling control over their overall magnitude and preventing overfitting to the data. The hyperparameters $o$ and $v^{2}$ determine the center and spread of the prior distribution, respectively, and can be selected based on prior knowledge. To capture an autoregressive process that depends on the distance $d_{ij}$ between consecutive discrete points in the continuous domain, we use a Normal prior with mean $\phi_{ij}\,\delta_{(i-1)j}$. This formulation reflects that past effects influence future effects. The priors for $\eta$, $\nu^{2}$, and $\tau^{2}$ are IG, reflecting uncertainty about the parameters controlling the smoothness, the random effect variance, and the observations.

Consider a scenario with 15 specific geographic locations, depicted in Figure 2. At each location, 200 discrete points are selected along the continuous domain of the function $t$. These points are associated with a measure $Y$, forming ordered pairs $(t,Y)$ that define the corresponding functions. It’s important to note that the spacing between the points is irregular and varies for each function, meaning it doesn’t follow a consistent pattern. Let us assume, without loss of generality, that the distances between discrete points in the functional domain are within the interval $(0,1)$. We will examine two situations: In the first case, we will use a Uniform distribution $U(0,1)$, and in the second, from a $Beta(1,2)$. The first option will result in an equal number of small and large distances, while the second will produce more small distances than large ones.

The generation of the synthetic data $Y_{s_{j}}(t_{ij})$, with $j=1,\cdots,15$ and $i=1,\ldots,n_{j}$, follows the model structure proposed in Equation (2.3), and is denoted as $\mathcal{M}_{BP_{p};\,\delta}$, which indicates a BP model of degree $p$ (order $p+1$) that includes the random effect $\delta$. The selection of the appropriate number of BP basis functions to use in a given dataset depends on various factors, such as the complexity of the data and the desired modeling accuracy. Although there is no fixed rule to determine the exact number of basis functions, it is generally accepted that increasing the value of $p$ increases flexibility to obtain a reasonable approximation. However, increasing $p$ also increases the complexity of the resulting polynomials, which can lead to numerical problems and higher computational costs. This study utilized a BP with $p=3$, corresponding to four basis functions. It’s important to emphasize that the model used to create the data is identical to the one used to fit the datasets.

The coefficients of the basis expansions, denoted by ${\bm{\theta}_{r}}=(\theta_{r1},\cdots,\theta_{r15})^{\top}$, have been derived from a normal distribution with mean vector $\mu_{\theta_{r}}\cdot{\bm{1_{15}}}$ and covariance matrix ${\bm{\Sigma}_{15}}$, where ${\bm{1_{15}}}$ represents a vector of ones with dimensions $15\times 1$. The specific mean values are $\mu_{\theta_{0}}=3$, $\mu_{\theta_{1}}=29$, $\mu_{\theta_{2}}=15$, and $\mu_{\theta_{3}}=7$, with $p=3$ and $r$ ranging from 0 to 3. The spatial correlation within the data is governed by the covariance matrix $\bm{\Sigma}_{15}$, constructed using the Gaussian covariance function $C(h)=\kappa^{2}\exp({-(\varphi\,h)^{2}})$, where $h$ denotes the Euclidean distance between the spatial locations $s_{j}$ and $s_{j^{*}}$. To illustrate the scenario of spatially dependent sets of curves, we consider a spatial variation $\kappa^{2}=2$, along with a moderate spatial correlation level $\varphi=1$. In addition, for the $\delta_{ij}$ structure, we use $\eta=0.2$ and variance $\nu^{2}=0.5$. Finally, the variance of the observations is assumed to be $\tau^{2}=1$.

Table 1 presents the prior along with the arguments assigned to each parameter and unknown quantity in the model. The pair of values (2,1) chosen for the first four Inverse Gamma distributions indicates high uncertainty about $\varphi$, $\eta$, $\tau^{2}$, and $\kappa^{2}$. This configuration results in distributions with a mean of 1 but an undefined variance. The fifth Inverse Gamma is defined by $a_{\nu^{2}}=3$ and $b_{\nu^{2}}=2$, which results in a prior distribution with a mean and variance of 1. This specification reflects an informative expectation about the magnitude of $\nu^{2}$; we expect the variability of the random effects $\delta_{ij}|\delta_{(i-1)j}$ to be moderate, i.e., neither too small nor excessively large. It’s essential because if $\nu^{2}$ were very large, the variability between consecutive observations would be too high, reducing the influence of the mean $\phi_{ij}\delta_{(i-1)j}$ on $\delta_{ij}$; see Equation (2.3). In addition, the hyperparameter $\mu_{\theta_{r}}$, a Normal prior with a mean of 0 and variance of $50^{2}$, is proposed to represent low certainty about the real values of the coefficient means.

To summarize the inference results, we use the posterior means as estimators. We obtain posterior samples using the MCMC algorithm with HMC dynamics implemented in Stan. We select every tenth iteration from 25,000 runs to reduce autocorrelation while discarding the initial 5,000 runs as the burn-in period. Each parameter is estimated using a single chain, and we visually confirm the convergence of all inspected chains during the analysis.

The objective of this initial simulation is to evaluate the effectiveness of the proposed model (as outlined in Equation 2.3) in recovering the dependent curves produced by the model, along with the corresponding configurations described in Subsection 3.1.

In Figure 3, boxplots display the results of 100 measurements of the Integrated Squared Error (ISE) metric (see Appendix A) obtained using the Monte Carlo (MC) scheme. These plots illustrate the approximation performance between smoothed and the target curves across two configurations of irregular domains. In Panel $(a)$, the observations are associated with discrete points $t_{ij}$, where the distances follow a Uniform distribution $U(0,1)$. In this case, there are fewer small distances compared to the $Beta(1,2)$ distribution, indicating a lower level of dependence between the observations of each curve. Conversely, Panel $(b)$ shows observations at discrete points where distances are generated from a $Beta(1,2)$ distribution, resulting in a higher concentration of smaller distances and thus indicating greater dependence between observations within each series.

The plots suggest that the model $\mathcal{M}_{BP_{3};\,\delta}$ accurately recovers the spatially correlated curves, regardless of the irregularity in the spacing of the domain. Notably, approximately 50% of the ISE values are below 0.15, as marked by the red line. Additionally, the boxplots in Panel $(a)$ exhibit slightly higher variability than those in Panel $(b)$. In summary, comparing the two panels reveals that the variation in ISE values is more stable under the $Beta(1,2)$ configuration, as indicated by the narrower range of extreme ISE values. On the other hand, the $U(0,1)$ displays a weaker correlation between observations within each series, resulting in a broader range of outliers. This implies that, in some simulations, the approximation performs significantly worse under the uniform distribution.

This study uses the same grid of geographical locations and the functional domain pattern for each curve described in the first simulation study (refer to Figure 2). The dependent curves are generated using a linear combination of 9 Fourier bases.

where the constant $\omega$ is related to the period $\mathcal{T}$ by the relation $\omega=2\pi/\mathcal{T}$. We assume tha $\epsilon_{s_{j}}(t_{ij})\sim N(0,\sigma^{2}=2)$ independently for all $i=1,\ldots,200$ and $j=1,\ldots,15$. In addition, $\bm{\beta_{l}}=(\beta_{l\,s_{1}},\ldots,\beta_{l\,s_{15}})^{\top}\sim N_{15}(10\cdot\bm{1}_{15},\bm{\Sigma_{15}})$, where $l=0,1,\ldots,8$ and $\bm{\Sigma}_{15}$ is defined according to the Exponential covariance function $C(h)=2\exp(-(1\cdot h))$. This study also explores the MC scheme with 100 replicas generated under these conditions. The simulated data sets, with and without noise, can be seen in Figure 4.

Each synthetic dataset is fitted using the proposed $\mathcal{M}_{BP_{3};\,\delta}$ model and the model without the random effect component $\mathcal{M}_{BP_{3}\,;\,\bullet}$. Both models employ four bases of BP functions of degree 3. The objective is to compare their fitting performances and highlight the importance of considering both spatial dependence and the associations between the discretely observed measurements of each curve. This consideration is motivated by the irregular spacing of the points in the functional domain.

In Figure 5, Panel $(a)$ shows that most of the ISE values presented in each curve’s boxplot are below 0.4. As the number of bases increases to 18, the median of the majority of the boxplots declines to values below 0.15, as illustrated in panel $(b)$. The variability of the values remains consistently similar across all 15 plots in both panels. This pattern, along with the values obtained from the ISE measurement, indicates the effectiveness of the proposed methodology in successfully recovering curves that exhibit spatial dependencies differing from the model structure $\mathcal{M}_{BP_{p};\,\delta}$. In other words, the smoothed trajectories align closely with the target curves. In contrast, panels $(b)$ and $(c)$ demonstrate that the performance of the $\mathcal{M}_{BP_{p};\,\bullet}$ model is less effective, mainly when using only 4 BP functional bases. This is due to the model’s failure to account for potential dependencies between observations within each series, as the $\mathcal{M}_{BP_{p};\,\bullet}$ model does not include a random effect component.

This second part of the study aims to assess how well the model can predict curves in geographical areas that have yet to be observed. To achieve this, we will use sites 1 and 11 from the grid shown in Figure 2 as the basis for our predictions. These locations are not part of the sample used to train the model, but we know their actual values, which will help us evaluate the accuracy of the predictions.

Figure 6 shows the average curve of 100 replicas of the a posteriori means at unobserved locations 1 and 11. These locations are known for having few close neighbors. The observed curves were created under a moderate spatial dependence ($\varphi=1$). The main goal of the analysis is to assess how well the model predicts these curves using different numbers of bases: 4, 12, and 18. First, it was observed that using a reduced number of basis functions is not enough to accurately capture the spatial structure, which leads to decreased prediction accuracy. However, when 12 BP bases of degree 11 were used, the model showed greater flexibility, resulting in better capture of the trends of the target trajectories. Finally, with 18 functional bases, the smoothed curves closely fit the target curves, demonstrating a remarkable approximation. Additionally, it was observed that, as the number of bases increased, the average (out of 100 replicates) of the 95% highest posterior density (HPD) intervals (indicated by the gray area) became progressively more compressed, indicating a decrease in uncertainty and increased accuracy in the model predictions.

It should be noted that there is no fixed methodology to determine the optimal number of basis functions. However, increasing their number provides greater flexibility in capturing trends in locations with few neighbors. However, this increase also entails higher computational costs, which can lead to problems related to numerical stability.

Figure 7 presents the boxplots for 100 measurements of the ISE metric, with one measurement for each replicate, which allows us to evaluate the proximity of the smoothed curves concerning the target functions. Firstly, when a reduced number of basis functions is used (Panel a), the prediction curve at location $s_{11}$, which is more isolated from neighboring locations, demonstrates better performance compared to $s_{1}$. However, increasing the number of basis functions enhances the model’s flexibility, allowing it to capture spatial dependence better and improve estimates at both locations, with $s_{1}$ exhibiting superior performance. However, this increased flexibility also results in more significant variability in the ISE values, especially in the $s_{11}$ boxplot (Panel b and c). Finally, as the number of basis functions increases, the ISE values decrease, indicating an improvement in overall performance. This effect is especially noticeable in $s_{1}$, whose proximity to other locations strengthens the inference process and improves smoothing. The Appendix B offers supplementary details on other scenarios explored in this study.