Positional Encoder Graph Quantile Neural Networks for Geographic Data

## Abstract

Abstract Positional Encoder Graph Neural Networks (PE-GNNs) are a leading approach for modeling continuous spatial data. However, they often fail to produce calibrated predictive distributions, limiting their effectiveness for uncertainty quantification.
We introduce the Positional Encoder Graph Quantile Neural Network ( PE-GQNN ), a novel method that integrates PE-GNNs, Quantile Neural Networks, and recalibration techniques in a fully nonparametric framework, requiring minimal assumptions about the predictive distributions. We propose a new network architecture that, when combined with a quantile-based loss function, yields accurate and reliable probabilistic models without increasing computational complexity. Our approach provides a flexible, robust framework for conditional density estimation, applicable beyond spatial data contexts. We further introduce a structured method for incorporating a KNN predictor into the model while avoiding data leakage through the GNN layer operation. Experiments on benchmark datasets demonstrate that PE-GQNN significantly outperforms existing state-of-the-art methods in both predictive accuracy and uncertainty quantification.

## 1 Introduction

Large spatial datasets are collected in a wide range of applications in economics (Anselin,, ), meteorology (Bi et al.,, ), urban transportation (Lv et al.,, ; Derrow-Pinion et al.,, ; Kashyap et al.,, ), social networks (Xu et al.,, ), e-commerce (Sreenivasa and Nirmala,, ) and other fields. Gaussian Processes (GPs) (Rasmussen and Williams,, ; Cressie and Wikle,, ) are a fundamental tool for modelling spatial data on continuous domains. They are flexible and interpretable models for unknown functions, both in spatial and more general regression settings. However, with time complexity $O(n^{3})$ and storage complexity $O(n^{2})$, naive GP methods quickly become intractable for large datasets. This has led to a large range of approximate inference methods, such as those based on sparse approximations to covariance or precision matrices (Reinhard Furrer and Nychka,, ; Lindgren et al.,, ), low rank approximations (Cressie et al.,, ) or nearest neighbour approximations (Vecchia,, ; Datta et al.,, ; Katzfuss and Guinness,, ).

Given the difficulty of GP computations, it’s of interest to explore scalable methods for large spatial datasets using neural networks (NNs) and to enhance their ability to quantify uncertainty. A state-of-the-art method for making spatial predictions using Graph Neural Networks (GNNs) is the Positional Encoder Graph Neural Network (PE-GNN) of Klemmer et al., (). Our contribution is to make three key modifications to the PE-GNN architecture to enhance its ability to make accurate spatial predictions and to quantify uncertainty. These modifications will be explained further below.

NNs are popular in data modeling and prediction tasks like computer vision and natural language processing (NLP). However, traditional NNs struggle to handle spatial dynamics or graph-based data effectively. GNNs (Kipf and Welling,, ; Veličković et al.,, ; Hamilton et al.,, ) offer a powerful and scalable method for applying NNs to graph-structured data. The idea is to share information through the edges of a graph, allowing nodes to exchange information during learning. GNNs are versatile and can uncover nonlinear relationships among inputs, hidden layers, and each node’s neighborhood information. The success of GNNs in spatial applications largely depends on the spatial graph construction, including choice of distance metric and the number of neighboring nodes, and traditional GNNs often struggle to model complex spatial relationships. To address this, Klemmer et al., () introduced the PE-GNN, which enhances predictive performance in spatial interpolation and regression. However, PE-GNN is not designed to provide a full probabilistic description of the target’s distribution, and assuming a Gaussian distribution for predictions can lead to poorly calibrated intervals, such as 80% intervals that fail to contain the true outcome 80% of the time. Recently, Bao et al., () proposed a new framework called Spatial Multi-Attention Conditional Neural Processes (SMACNPs) for spatial small sample prediction tasks. SMACNPs use GPs parameterized by NNs to predict the target variable distribution, which enables precise predictions while quantifying the uncertainty of these predictions.

Methods based on quantile regression are an alternative approach to probabilistic forecasting making rapid progress in recent years. Si et al., () introduced a novel architecture for estimating generic quantiles of a conditional distribution, proposing a set of objective functions that lead to enhancements in density estimation tasks. In one dimension, this method produces quantile function regression and cumulative distribution function regression. Kuleshov and Deshpande, () argue that the method of Si et al., () is inefficient with high-dimensional predictors. To address this, they modify the original formulation to incorporate a post hoc recalibration procedure whereby an auxiliary model recalibrates the predictions of a trained model. The first model outputs features, usually summary statistics like quantiles, representing a low-dimensional view of the conditional distribution. The auxiliary model, the recalibrator, uses these features as input to produce calibrated predictions using Si et al.’s quantile function regression framework. The main drawback is that it requires training two separate models, each needing its own training set.

Our work makes three contributions. (1) We propose a new architecture that merges the two-step procedure of Kuleshov and Deshpande, () into a single model by postponing the concatenation of the $\tau$ value proposed by Si et al., (). In this way, we enhance the network’s ability to model uncertainty and introduce a regularization mechanism. The model becomes robust to high-dimensional predictor spaces, even though few assumptions are made about the form of the target’s conditional distribution. This change allows a single model to fully describe the predictive conditional distribution and to generate quantile predictions and prediction intervals as byproducts. It can be applied to any context, not just spatial regression or GNNs. We show how to integrate this strategy into the PE-GNN framework to create an intrinsically calibrated model with no extra computational cost. (2) We introduce a structural change to PE-GNN. Instead of applying the GNN operator to the concatenation of the nodes’ features and the spatial embedding, we apply it only to the features. (3) In PE-GNN, the GNN operator uses neighbours’ features to create new node representations but does not include the target value of neighboring nodes. Our third contribution introduces the mean target value of a node’s neighbours as a feature after the GNN layers, closer to the output. This allows the model to use neighboring observations of the target variable when making predictions.

The structure of this work is as follows: Section [2](#S2) offers a brief background overview, Section [3](#S3) outlines the proposed method for geographic data prediction, Section [4](#S4) shows experimental results on three real-world datasets, and Section [5](#S5) concludes.

## 2 Background

Positional Encoder    Inspired by the Transformer architecture (Vaswani et al.,, ) for geographic data (Mai et al.,, ), PE-GNN (Klemmer et al.,, ) employs a PE with two components: a sinusoidal transformation and a fully-connected NN. The first is a deterministic transformation formed by the concatenation of sinusoidal functions, including variations in frequency and scale. The spatial dimensions (typically represented as latitude and longitude) are handled separately. The second component is a fully-connected NN, denoted $NN(\Theta_{PE})$, taking the output produced by the sinusoidal transformation as input and processing it through a fully-connected NN. Let $\bm{C}_{B}=[\bm{c}_{1},\ldots,\bm{c}_{n_{B}}]^{\top}$ be the matrix containing the spatial coordinates of a batch of datapoints, typically of dimension $n_{B}\times 2$, where each $\bm{c}_{i}$ corresponds to the pair $(\text{latitude}_{i},\text{ }\text{longitude}_{i})$. This transformation results in the desired vector space representation, thereby generating the coordinate embedding matrix $\bm{C}_{B}^{emb}=PE(\bm{C}_{B},\sigma_{min},\sigma_{max},\Theta_{PE})=NN(ST(\bm{C}_{B},\sigma_{min},\sigma_{max}),\Theta_{PE})$.

Graph Neural Network    GNNs are powerful and scalable solutions for representation learning and inference with graph-structured data. They leverage the topological structure of correlations between nearby graph nodes and represent each node in a latent space embedding suitable for the specific downstream task (Wu et al.,, ). Popular GNN architectures use this graph structure to update the embeddings of each node, considering both the features of each node and its neighbors, in an iterative process (Wu et al.,, ). The first step comprises aggregating features from each node’s neighbours. After aggregation, we combine each node’s prior representation with the output of the first step. The initial embedding of each node is its feature vector, so $\bm{H}^{(0)}_{B}=\bm{X}_{B}$. Then, for each GNN layer $k\in\{1,\ldots,K\}$, an iteration of the two step process described above is executed.

The most popular GNN architectures follow this backbone, but differ in the way they aggregate neighbours messages and update the embeddings. Graph Convolutional Networks (GCNs) (Kipf and Welling,, ) are inspired by the convolution operation from Convolutional Neural Networks (CNNs). For weighted graphs, GCN layer $k$ has the following update equation

$$ $\bm{H}^{(k)}_{B}=f^{(k)}\left(\bm{D}^{-1/2}_{B}\left[\bm{A}_{B}+\bm{I}_{B}\right]\bm{D}^{-1/2}_{B}\bm{H}^{(k-1)}_{B}\bm{W}^{(k)}\right),\text{ for }k\in\{1,\ldots,K\}.$ (1) $$

Here, $f^{(k)}$ is an activation function (e.g., ReLU) and $\bm{W}^{(k)}$ is a matrix of learnable parameters, while the adjacency matrix $\bm{A}_{B}$ describes the connectivity of the constructed graph, where unconnected nodes have a value of 0, and connected nodes have an edge weight computed from their distance. $\bm{D}_{B}$ is the so-called *degree* matrix and $\bm{I}_{B}$ denotes the identity matrix.

Positional Encoder Graph Neural Network    Klemmer et al., () proposed a novel approach for applying GNNs to spatial data: PE-GNN. A PE was introduced that takes the set of spatial coordinates for each datapoint as input and produces a vector representing the learned spatial embedding. This vector is then column concatenated with the node features before application of the GNN operator. Thus, for a given batch $B$ of randomly sampled datapoints, the input to the first GNN layer is $\bm{H}^{(0)}_{B}=concat\left(\bm{X}_{B},\bm{C}^{emb}\right)$. PE-GNN also predicts the Local Moran’s I (Anselin,, ) as an auxiliary task (Klemmer and Neill,, ). Another innovation lies in the training process, as PE-GNN uses a batch-based procedure. At each training step, a random batch $B$ of nodes is sampled, given by $\{p_{1},\ldots,p_{n_{B}}\}\in B$. Using only the nodes belonging to the batch, the entire process of constructing the training graph, generating the spatial embedding, column concatenating with the features, and applying the GNN operator is carried out. The loss function used by Klemmer et al., () is given by $\mathcal{L}_{B}=MSE\left(\hat{\bm{y}}_{B},\bm{y}_{B}\right)+\lambda MSE\left(I(\hat{\bm{y}}_{B}),I(\bm{y}_{B})\right)$, where $\lambda$ denotes the auxiliary task weight.

Quantile regression Koenker and Bassett Jr, () proposed a linear quantile regression model to estimate conditional distribution quantiles. It uses the pinball loss $\rho_{\tau}(r_{i})=\max{(\tau r_{i},(\tau-1)r_{i})}$, where $r_{i}=y_{i}-\hat{q}_{i}(\tau)$, $\hat{q}_{i}(\tau)=\bm{X}_{i}\hat{\bm{\beta}}$, and $\tau$ is the desired cumulative probability associated with the predicted quantile $\hat{q}_{i}(\tau)$. The pinball loss for the $i$-th observation is $\rho_{\tau}(r_{i})$. The loss over a dataset is the average $\rho_{\tau}(r_{i})$ value over all datapoints. A natural extension of quantile linear regression is quantile neural networks (QNNs). This approach is illustrated in Figure [1(a)](#S2.F1.sf1), which seeks to estimate the conditional quantiles for a pre-defined grid $(\tau^{1},\ldots,\tau^{d})$. Each quantile is estimated by an independent model (Figure [1(a)](#S2.F1.sf1)). This can lead to quantile predictions with quantile crossing (e.g., a median prediction lower than the first quartile prediction).

Rodrigues and Pereira, () proposed an approach that outputs multiple predictions: one for the expectation and one for each quantile of interest (Figure [1(b)](#S2.F1.sf2)). The loss function is:

$$ $\mathcal{L}=\frac{1}{d+1}\left[\text{MSE}\left(\hat{\bm{y}},\bm{y}\right)+\sum_{i=1}^{n}\sum_{j=1}^{d}\frac{\rho_{\tau^{j}}\left(y_{i}-\hat{q}_{i}(\tau^{j})\right)}{n}\right].$ (2) $$

Figure: (a) Non-linear quantile regression using NN.
Refer to caption: /html/2409.18865/assets/Figures/Chap-2/QNN.png

Si et al., () proposed a method to generate a model that is independent of quantile selection (Figure [1(c)](#S2.F1.sf3)). For each datapoint sampled during training, $d$ Monte Carlo samples $\tau\sim U(0,1)$ are drawn. Each sample is concatenated with the datapoint features to obtain a corresponding quantile estimate, so for each datapoint there are $d$ predicted quantiles. The loss function is similar to Eqn. [2](#S2.E2), but they predict random quantiles $\mathcal{L}=\frac{1}{n\cdot d}\sum_{i=1}^{n}\sum_{j=1}^{d}\rho_{{\tau}_{i}^{j}}(y_{i}-\hat{q}_{i}({\tau}_{i}^{j}))$. As the network learns, it becomes able to provide
a direct estimate to any quantile of interest. Hence, this procedure outputs an inherently calibrated model suitable for conditional density estimation. However, Kuleshov and Deshpande, () argue this method is inefficient in mid-to-high predictor space dimensions.

Kuleshov and Deshpande, () adapted the architecture from Si et al., () into a two-step process for larger predictor spaces (Figure [1(d)](#S2.F1.sf4)). First, a model is trained to take the original features as inputs and generate low-dimensional representations of the predicted distribution. Next, a recalibrator is trained using these *new* features by minimizing the estimated expected pinball loss over $\tau$. During inference, the recalibrator takes the new features and an arbitrary $\tau$ as inputs to produce the quantile prediction. This method is highly dependent on the choice of recalibrator features.

## 3 Method

In this work, we propose a novel approach to spatial data prediction tasks: the Positional Encoder Graph Quantile Neural Network (PE-GQNN). Algorithm [1](#alg1) shows the step-by-step procedure to train a PE-GQNN model, and Figure [2](#S3.F2) illustrates its complete pipeline. Here, each rectangle labeled "GNN" and "LINEAR" represents a set of one or more neural network layers, with the type of each layer defined by the title inside the rectangle. At each layer, a nonlinear transformation (e.g. ReLU) may be applied. Each datapoint $p_{i}$ comprises three components $p_{i}=\{y_{i},\bm{x}_{i},\bm{c}_{i}\}$. The component $y_{i}$ is the target variable, and as the focus here is regression, then $y_{i}$ is a continuous scalar. Additionally, $\bm{x}_{i}$ is the feature vector and $\bm{c}_{i}$ contains the geographical coordinates associated with observation $i$.

Figure: Algorithm 1 PE-GQNN training

Figure: Figure 2: PE-GQNN compared to PE-GNN and GNN.
Refer to caption: /html/2409.18865/assets/Figures/Chap-3/PEGQNN-v11.png

After initializing the model and hyperparameters, the first step of PE-GQNN is to randomly sample a batch $B$ of datapoints, $p_{1},\ldots,p_{n_{B}}$. The batch can be fully represented by the target ${\bm{y}_{B}}_{(n_{B}\times 1)}$, features ${\bm{X}_{B}}_{(n_{B}\times p)}$, and coordinates matrices ${\bm{C}_{B}}_{(n_{B}\times 2)}$, respectively. The next step uses the matrix of geographical coordinates $\bm{C}_{B}=\left[\bm{c}_{1},\ldots,\bm{c}_{n_{B}}\right]^{\top}$ to obtain spatial embeddings for each datapoint (Algorithm [1](#alg1), Step [14](#alg1.l14)). This process receives $\bm{C}_{B}$ as input, and after passing through deterministic sinusoidal transformations and a fully-connected NN, outputs the spatial embedding matrix of the batch ${\bm{C}_{B}^{emb}}_{(n_{B}\times u)}$, containing the spatial context of each pair of coordinates. $\bm{C}_{B}$ is also used to compute the distance between each pair of datapoints (Algorithm [1](#alg1), Step [15](#alg1.l15)). From these distances and a predefined number of nearest neighbors, a graph can be constructed, with each datapoint as a node and edge weights computed from the distances, leading to the batch adjacency matrix $\bm{A}_{B}$.

At Step [22](#alg1.l22) of Algorithm [1](#alg1), the first distinction between PE-GQNN and PE-GNN arises: instead of using the concatenation of the feature matrix and the spatial embedding as the input for the GNN operator, we apply the GNN operator only to the feature matrix $\bm{X}_{B}$. One or more fully-connected layers are then used to reduce the feature embedding dimensionality. This process receives the constructed graph and the batch feature matrix ${\bm{X}_{B}}_{(n_{B}\times p)}$ as input and yields an embedding matrix of features as output: ${\bm{X}_{B}^{emb}}_{(n_{B}\times g)}$. This modification applies the GNN operators exclusively to the features, without smoothing out the PEs. The GNN layers can be of any desired type. Step [23](#alg1.l23) of Algorithm [1](#alg1) performs a column concatenation between the feature embedding ${\bm{X}_{B}^{emb}}_{(n_{B}\times g)}$ and the output obtained from the PE: ${\bm{C}_{B}^{emb}}_{(n_{B}\times u)}$ (Figure [2](#S3.F2)). This concatenation results in the matrix ${\bm{L}_{B}}_{(n_{B}\times(g+u))}$. After that, the other innovations of PE-GQNN come into play.

First, we use one or more fully-connected layers (Algorithm [1](#alg1), Step [24](#alg1.l24)) to reduce the dimensionality of $\bm{L}_{B}$, making it suitable for two of the three innovations in PE-GQNN. This set of fully-connected layers outputs the matrix ${\bm{\phi}_{B}}_{(n_{B}\times s)}$, which is then combined with $\bar{\bm{y}}_{B}$ and $\bm{\tau}_{B}$. $\bar{\bm{y}}_{B}$ represents a vector with one scalar for each datapoint in the batch, containing the mean target variable among the training neighbours for each node. It is computed using the graph constructed in previous steps (Algorithm [1](#alg1), Step [19](#alg1.l19)), and has dimensions $n_{B}\times 1$. It is comparable to a vector of predictions generated by a KNN regression model, where neighbours are determined using the distance calculated from geographical coordinates. Here, we used the simple average due to its relationship with KNN prediction; however, one could use a weighted average via the adjacency matrix $\bm{A}_{B}$. We introduce this input at a later stage to avoid data leakage. If the GNN operator received $\bar{\bm{y}}_{B}$ as input, after completing the message passing process in each GNN layer, the true node target value would inadvertently be transmitted to its neighbours, creating potential data leakage (Appleby et al.,, ).

In the same layer where $\bar{\bm{y}}_{B}$ is introduced, we apply a similar approach to Si et al., () to make PE-GQNN an inherently calibrated model suitable for probabilistic and quantile predictions. For each batch $B$, we create a $n_{B}\times 1$
vector ${\bm{\tau}_{B}}_{(n_{B}\times 1)}=[\tau_{1},\ldots,\tau_{n_{B}}]^{\top}$ of random $U(0,1)$ draws (Algorithm [1](#alg1), Step [25](#alg1.l25)). Then, we column concatenate $\bm{\phi}_{B}$ with $f\left(\bm{\tau}_{B}\right)$ and $\bar{\bm{y}}_{B}$ to create $\bm{\widetilde{\phi}}_{B}{}_{(n_{B}\times(s+2))}$ (Algorithm [1](#alg1), Step [26](#alg1.l26)), where $f(\makebox[4.30554pt]{{$\cdot$}})$ is an activation function. Here we propose use of $f(\makebox[4.30554pt]{{$\cdot$}})=logit(\makebox[4.30554pt]{{$\cdot$}})$, to facilitate the network’s learning. Subsequently, forward propagation is computed (Algorithm [1](#alg1), Step [27](#alg1.l27)) in one or more fully-connected layers, outputting predicted quantiles for each datapoint in the batch. The batch loss is the one proposed by Si et al., (), but with $d=1$ for the $\tau$ values.

This procedure aims to improve the model’s ability to learn the conditional probability distribution of the target variable, enhancing uncertainty estimation and quantile predictions. Instead of introducing $\tau$ values alongside features at the network’s input, as suggested by Si et al., (), we delay their entry into a reduced latent dimension to boost learning. This adjustment makes PE-GQNN suitable for both low- and high-dimensional predictor spaces. It also improves on the Kuleshov and Deshpande, () approach by merging the two-network process into a single, intrinsically calibrated model.

Incorporating $\tau$ values into the model architecture improves its ability to model uncertainty and serves as a regularization mechanism (Rodrigues and Pereira,, ). The use of pinball loss for quantile regression acts as a natural regularizer, producing a detailed description of the predictive density beyond just mean and variance estimation. For predictions, the quantile of interest, $\tau$, must be given, along with the basic data components (e.g. $\tau=0.25$ gives the first quartile). If interest is in predicting multiple quantiles for the same observation, the input can be propagated up to the layer where $\bm{\tau}$ is introduced. For each quantile of interest, propagation can be limited to the final layers.

Target domain    The final layer should use an activation function coherent with the domain of the target variable, ensuring model outputs are valid for target distribution support. E.g., an exponential function could be appropriate if the target variable is continuous, unbounded and positive.

Quantile crossing    This phenomenon occurs when estimated quantile functions for different quantile levels ($\tau$) intersect, violating the requirement that higher quantiles be greater than or equal to lower quantiles. In PE-GQNN, by utilizing the same latent representation up to the layer where the quantile level ($\tau$) is introduced, the architecture adopts a hard-parameter sharing multi-task learning framework. This severely mitigates the problem of quantile crossings by constraining the flexibility of independent quantile regression neural network models. If $\tau$ is introduced at the prediction layer, it is guaranteed that quantile crossings will be absent, as the layer equation would be

$$ $\hat{q}_{i}(\tau)=f\left(bias+w_{\tau}\tau+w_{\bar{y}_{i}}\bar{y}_{i}+\sum_{j=1}^{neurons}w_{j}u_{j}\right),\text{ }\forall i\in{1,\ldots,n_{B}}.$ (3) $$

Here, $neurons$ denotes the number of neurons in the prediction layer, excluding $\tau$ and $\bar{y}_{i}$. $bias$, $w_{\tau}$, $w_{\bar{y}_{i}}$, and $\{w_{j}\}$ are the prediction layer parameters, and $\{u_{j}\}$ are the activation values from the previous layer. Commonly, $f$ is chosen to be monotonic, resulting in a monotonic relationship between $\tau$ and $\hat{q}_{i}(\tau)$. When $\tau$ is introduced at a layer proximal to, but preceding, the prediction layer, the results in Section [4](#S4) suggest our approach is not prone to suffer from quantile crossing.

Number of Monte Carlo samples    When applying the framework proposed by Si et al., (), we chose to use $d=1$ for the $\tau$ values. Let $\mathcal{L}(\theta,\tau,\bm{x},y)$ be the loss function for a given quantile $\tau\sim\text{U}(0,1)$ and an observed pair $(\bm{x},y)\sim D_{\text{data}}$, where $D_{\text{data}}$ denotes the full data generative process. On each training iteration, we minimize $\mathcal{L}_{B}$, which, by the Law of Large Numbers, converges to $\tilde{\mathcal{L}}(\theta)=\mathbb{E}_{\tau,\bm{x},y}\mathcal{L}(\theta,\tau,\bm{x},y)$, as the batch size, ${n_{B}}$, goes to infinity. Therefore, the gradients converge to the same value for any $d$, provided that ${n_{B}}\to\infty$. This choice $(d=1)$ simplifies the implementation without sacrificing performance, as shown in Section [4](#S4).

## 4 Experiments

### 4.1 Experimental setup

PE-GQNN was implemented using PyTorch (Paszke et al.,, ) and PyTorch Geometric (Fey and Lenssen,, ). The source code is available at: https://github.com/WilliamRappel98/PE-GQNN. We conducted comprehensive simulations to explore the prediction performance and other properties of the proposed model, comparing it with state-of-the-art methods. Computation was performed on an Intel i7-7500U processor with 16 GB of RAM, running Windows 10.

Candidate models    The experiment was designed to compare six primary approaches for addressing spatial regression problems across three distinct real-world datasets. Table [1](#S4.T1) lists each candidate model and their applicable datasets.
All models were trained using the Adam optimizer (Kingma and Ba,, ), early stopping and, for all GNN-based models, $k=5$ nearest-neighbours.

**Table 1: Summary of candidate models.**
| Approach | Model | Type | Components | Loss | Datasets |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| PE | Moran’s I | $\bm{\tau}$ | Structure | $\bm{\bar{y}}$ |  |  |  |  |  |
| I | GNN | GNN | No | No | No | No | No | MSE | All |
| II | PE-GNN $\lambda=\text{best}$ | GNN | Yes | Yes | No | No | No | $\text{MSE}_{y}+\lambda\text{MSE}_{I(y)}$ | All |
| III | PE-GQNN $\tau$ | GNN | Yes | No | Yes | No | No | Pinball | California |
| IV | PE-GQNN $\tau$, Structure | GNN | Yes | No | Yes | Yes | No | Pinball | California |
| V | PE-GQNN | GNN | Yes | No | Yes | Yes | Yes | Pinball | All |
| VI | SMACNP | GP | No | No | No | No | No | Log Likelihood | All but 3D road |

Approach I involves the traditional application of GNNs to geographic data. Three types of GNN layers were considered: GCNs (Kipf and Welling,, ), GATs (Veličković et al.,, ), and GSAGE (Hamilton et al.,, ). For each of these, the architecture remains consistent to facilitate performance comparisons: two GCN/GAT/GSAGE layers with ReLU activation and dropout, followed by a linear prediction layer.

Approach II involves the application of PE-GNN (Klemmer et al.,, ) with optimal weights for each dataset and layer type combination, as demonstrated by the experimental findings of Klemmer et al., (). The GNN architecture used is the same as for approach I. It was implemented using the code available at: https://github.com/konstantinklemmer/pe-gnn.

Approach III is similar to PE-GNN but augmented with the first innovation proposed in this study: the quantile regression framework described in Section [3](#S3) is applied. Approach IV is similar to III, but augmented with an additional innovation: the structural alteration in the model’s architecture, where the GNN operator is applied only to the features. Approach V, which is the primary focus of this research, explores the utilization of PE-GQNN. The PE and GNN layers’ architectures remain identical to the previous approaches, with any alterations limited to the proposed innovations.

Finally, a benchmark approach that does not use GNNs but was recently proposed for modelling spatial data will be considered as approach VI: SMACNPs. This approach, proposed by Bao et al., (), has demonstrated superior predictive performance, surpassing GPs models in the three real-world datasets considered. This model was implemented following the specifications of Bao et al., (), using the code available at: https://github.com/bll744958765/SMACNP.

Approaches I and II do not inherently provide predicted conditional distributions. However, as they optimize the MSE metric, they implicitly learn a Maximum Likelihood Estimate (MLE) of a Gaussian model. Thus, the predictive distribution considered for these approaches was a Gaussian distribution centered on the point prediction with variance equal to the MSE of the validation set. For computational simplicity in the experiments, instead of calculating $\bar{\bm{y}}_{B}$ for each batch, we pre-calculated $\bar{\bm{y}}$ using the entire training set.

Performance metrics    We evaluate predictive accuracy using Mean Squared Error (MSE) and Mean Absolute Error (MAE). To assess calibration of the predictive distributions, we use Mean Pinball Error, $\text{MPE}=\frac{1}{n}\sum_{i=1}^{n}\rho_{\tau_{i}}\left(y_{i}-\hat{q}_{i}(\tau_{i})\right)$, where $\tau_{i}\sim U(0,1)$, and the Mean Absolute Distance of the Empirical Cumulative Probability, $\text{MADECP}=\frac{1}{99}\sum_{j=1}^{99}\left|\tau^{j}-\frac{1}{n}\sum_{i=1}^{n}\mathds{1}\left[y_{i}\leq\hat{q}_{i}(\tau^{j})\right]\right|$. For quantile predictions of a calibrated model for a given $\tau$, the proportion of observed values less than or equal to the predicted quantile should approximate $\tau$. Evaluating the MADECP helps determine whether the predicted quantiles are accurate and consistent across the entire space.

Figure: (a) GCN.
Refer to caption: /html/2409.18865/assets/Figures/Chap-4/california-housing_GCN_validation_mse.png

### 4.2 California Housing

This dataset comprises pricing information for $>$20,000 residential properties in California, recorded during the 1990 U.S. census (Pace and Barry,, ). The main objective is a regression task: predict housing prices, $y$, through the incorporation of six predictive features, $\bm{x}$, and geographical coordinates, $\bm{c}$. The predictive features are neighborhood income, house age, number of rooms, number of bedrooms, occupancy and population. All models were trained and evaluated using 80% of the data for training, 10% for validation, and 10% for testing. In the case of SMACNP, to adhere to the specifications of Bao et al., (), a training subsample was extracted to represent 10% of the entire dataset. The validation MSE curves throughout training are shown in Figure [3](#S4.F3). The number of training epochs and final test dataset performance metrics are in Table [2](#S4.T2).

**Table 2: Performance metrics on the California Housing test set.**
| Model | Epochs | Parameters | MSE | MAE | MPE | MADECP |
| --- | --- | --- | --- | --- | --- | --- |
| GCN | 441 | 1,313 | 0.0222 | 0.1101 | 0.0403 | 0.0475 |
| PE-GCN $\lambda=\text{best}$ | 170 | 24,129 | 0.0179 | 0.0935 | 0.0354 | 0.0450 |
| PE-GQCN $\tau$ | 179 | 25,217 | 0.0179 | 0.0914 | 0.0351 | 0.0428 |
| PE-GQCN $\tau$, Structure | 264 | 26,169 | 0.0138 | 0.0800 | 0.0302 | 0.0287 |
| PE-GQCN | 76 | 26,201 | 0.0114 | 0.0686 | 0.0272 | 0.0262 |
| GAT | 398 | 1,441 | 0.0227 | 0.1099 | 0.0410 | 0.0586 |
| PE-GAT $\lambda=\text{best}$ | 120 | 24,290 | 0.0183 | 0.0930 | 0.0352 | 0.0476 |
| PE-GQAT $\tau$ | 136 | 25,345 | 0.0179 | 0.0926 | 0.0355 | 0.0413 |
| PE-GQAT $\tau$, Structure | 261 | 26,297 | 0.0140 | 0.0829 | 0.0312 | 0.0193 |
| PE-GQAT | 68 | 26,329 | 0.0114 | 0.0685 | 0.0268 | 0.0254 |
| GSAGE | 348 | 2,529 | 0.0170 | 0.0945 | 0.0349 | 0.0569 |
| PE-GSAGE $\lambda=\text{best}$ | 222 | 27,426 | 0.0114 | 0.0732 | 0.0280 | 0.0464 |
| PE-GQSAGE $\tau$ | 243 | 28,481 | 0.0113 | 0.0686 | 0.0266 | 0.0478 |
| PE-GQSAGE $\tau$, Structure | 224 | 27,385 | 0.0100 | 0.0632 | 0.0248 | 0.0314 |
| PE-GQSAGE | 160 | 27,417 | 0.0089 | 0.0596 | 0.0229 | 0.0288 |
| SMACNP | 70 | 748,482 | 0.0160 | 0.0881 | 0.0466 | 0.1481 |

As shown in Table [2](#S4.T2), PE-GQNN achieves state-of-the-art performance metrics, with major improvements over traditional GNN, PE-GNN and SMACNP. For the GSAGE layers, PE-GQSAGE achieved the lowest MSE, MAE and MPE. For this type of layer, which gave the best results overall, we still encounter considerable relative improvements from PE-GQSAGE in comparison with PE-GSAGE, with a reduction of 22% in MSE, 19% in MAE, 18% in MPE, and 38% in MADECP. We can also explore in-depth, the contribution of each specific innovation. The $\tau$ innovation, which corresponds to the application of the quantile regression framework proposed by Si et al., (), improved the calibration of quantile predictions, reducing MPE and MADECP. The structural innovation, which involves applying the GNN operator only to the features, is instrumental in enhancing prediction performance and improving the calibration, as evidenced by reduced MSE, MAE, MPE and MADECP. Finally, the use of training neighbours’ target mean as a feature introduced at one of the last network layers also further improved the model. This last innovation also accelerated convergence during model training, requiring fewer epochs.

Figure: (a) PE-GQSAGE densities.
Refer to caption: /html/2409.18865/assets/Figures/Chap-4/california-housing_density_v2.png

Figure [4](#S4.F4) presents plots that elucidate the behavior of the PE-GQSAGE quantile predictions. Figure [4(a)](#S4.F4.sf1) illustrates the predicted density of a subsample of 10 observations from the test set. For each observation of this sample, the cumulative distribution function was approximated via the quantile predictions using $\tau$ values in $\left[0.001,0.002,\ldots,0.999\right]^{\top}$. While parametric models presume a rigid structure for their outputs (such as a Gaussian distribution), which constrains their expressiveness, for PE-GQNN, no assumptions are made about the form of the predictive distribution. However, as shown in Figure [4(a)](#S4.F4.sf1), despite the absence of explicit model restrictions, the model produced symmetric distributional shape across predictions, similar to a Gaussian distribution, in this case.

For all test set observations, we verified that no quantile crossings were observed in any of the PE-GQNN models, i.e., all predicted quantiles are monotonically increasing with respect to $\tau$, aligning with the expectations described in Section [3](#S3).

Lastly, Figure [4(b)](#S4.F4.sf2) displays the empirical cumulative probability (ECP) for the test dataset quantile predictions using each of the 99 $\tau$ values in $[0.01,0.02,\ldots,0.99]^{\top}$. This type of plot was proposed by Kuleshov et al., (). The closer a model gets to the dashed diagonal line, the closer the $\tau$ values and the ECP. The Gold Standard represents one Monte Carlo draw from a perfectly specified model, where for each quantile level, the ECP is the observed success rate in $n$ Bernoulli trials with a success probability of $\tau$, where $n$ is the number of test set instances. It is evident that PE-GQSAGE has by far the best calibration performance. This is particularly notable when compared to SMACNP, which exhibits substantial calibration deficiencies due to its tendency to overestimate the variance component.

### 4.3 All datasets

Experiments were conducted on two other geographic datasets used by Klemmer et al., () and Bao et al., (). The Air Temperature dataset (Hooker et al., ()) contains geographical coordinates for $\sim$3,000 meteorological stations worldwide, with the goal of predicting mean temperatures ($y$) using mean precipitation levels ($x$). Models were trained with 80% of the data, with 10% for validation and testing each, while SMACNP used a 30% subsample for training, following the specifications of Bao et al., (). The 3D road dataset (Kaul et al., ()), includes $>$ 430,000 points with latitude, longitude, and altitude for the Jutland, Denmark road network. The task is to interpolate altitude ($y$) using latitude and longitude ($\bm{c}$). The data were split into 90% for training, 1% for validation, and 9% for testing. SMACNP metrics are not reported due to high computational costs.

**Table 3: Performance metrics from three different real-world datasets.**
| Model | California Housing | Air Temperature | 3D road |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| MSE | MAE | MPE | MADECP | MSE | MAE | MPE | MADECP | MSE | MAE | MPE | MADECP |  |
| GCN | 0.0222 | 0.1101 | 0.0403 | 0.0475 | 0.0224 | 0.1158 | 0.0427 | 0.0334 | 0.0170 | 0.1029 | 0.0358 | 0.0560 |
| PE-GCN $\lambda=\text{best}$ | 0.0179 | 0.0935 | 0.0354 | 0.0450 | 0.0045 | 0.0467 | 0.0189 | 0.0640 | 0.0032 | 0.0406 | 0.0151 | 0.0476 |
| PE-GQCN | 0.0114 | 0.0686 | 0.0272 | 0.0262 | 0.0025 | 0.0327 | 0.0119 | 0.0713 | 0.0001 | 0.0053 | 0.0022 | 0.0439 |
| GAT | 0.0227 | 0.1099 | 0.0410 | 0.0586 | 0.0233 | 0.1166 | 0.0434 | 0.0497 | 0.0170 | 0.1030 | 0.0359 | 0.0601 |
| PE-GAT $\lambda=\text{best}$ | 0.0183 | 0.0930 | 0.0352 | 0.0476 | 0.0058 | 0.0566 | 0.0209 | 0.0960 | 0.0035 | 0.0430 | 0.0163 | 0.0551 |
| PE-GQAT | 0.0114 | 0.0685 | 0.0268 | 0.0254 | 0.0025 | 0.0340 | 0.0143 | 0.0677 | 0.0001 | 0.0053 | 0.0022 | 0.0545 |
| GSAGE | 0.0170 | 0.0945 | 0.0349 | 0.0569 | 0.0223 | 0.1152 | 0.0431 | 0.0361 | 0.0170 | 0.1031 | 0.0358 | 0.0582 |
| PE-GSAGE $\lambda=\text{best}$ | 0.0114 | 0.0732 | 0.0280 | 0.0464 | 0.0037 | 0.0449 | 0.0169 | 0.0720 | 0.0032 | 0.0422 | 0.0146 | 0.0417 |
| PE-GQSAGE | 0.0089 | 0.0596 | 0.0229 | 0.0288 | 0.0023 | 0.0326 | 0.0130 | 0.0785 | 0.0001 | 0.0054 | 0.0022 | 0.0786 |
| SMACNP | 0.0160 | 0.0881 | 0.0466 | 0.1481 | 0.0018 | 0.0290 | 0.0391 | 0.2160 | — - | — - | — - | — - |

Table [3](#S4.T3) showcases the experimental results obtained from all three datasets: California Housing, Air Temperature, and 3D road. Each GNN layer’s performance is evaluated across three approaches: the traditional GNN, PE-GNN, and PE-GQNN. The PE-GQNN models incorporate all three innovations discussed in Section [3](#S3). Additionally, we include the SMACNP results as a benchmark model based on GPs. PE-GQNN consistently outperforms both traditional GNN and PE-GNN across all datasets and GNN backbones. In every dataset, the PE-GQNN innovations lead to significant reductions in MSE, MAE, and MPE. In the California Housing dataset, PE-GQNN consistently outperforms SMACNP in predictive accuracy and provides enhanced uncertainty quantification across all types of GNN layers. Conversely, for the Air Temperature dataset, SMACNP achieves the lowest MSE and MAE but suffers from significantly uncalibrated predictions, reflected by a much higher MPE and MADECP compared to PE-GQNN.

## 5 Conclusion

In this work, we have proposed the Positional Encoder Graph Quantile Neural Network (PE-GQNN) as an innovative framework to enhance predictive modeling for geographic data. Through a series of rigorous experiments on real-world datasets, we have demonstrated the significant advantages of PE-GQNN over competitive methods. The empirical results underscored the capability of PE-GQNN to achieve lower MSE, MAE, and MPE across all datasets and GNN backbones compared to traditional GNN and PE-GNN. Notably, PE-GQNN demonstrated substantial improvements in predictive accuracy and uncertainty quantification, as evidenced by its consistent performance in quantile calibration metrics such as MPE and MADECP. The PE-GQNN framework’s ability to provide a full description of the predictive conditional distribution, including quantile predictions and prediction intervals, provides a notable improvement in geospatial machine learning. PE-GQNN provides a solid foundation for future advancements in the field of geospatial machine learning.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - Please elaborate on the limitation of the proposed method (e.g., computational cost).
- Is the proposed method robust against noise?
soundness: 3
strengths: - PE-GQNN combines quantile prediction and distribution recalibration in a single model (not two-stage models), enhancing both predictive accuracy and calibration efficiency.
- By limiting GNN operations to specific features and introducing target values near the output layer, the model effectively prevents data leakage and improves computational efficiency.
- The authors use pinball loss for quantile regression, which allows one to provide a regularization effect, improving prediction reliability across diverse quantile levels.
- Extensive experiments on multiple real-world datasets demonstrate PE-GQNN’s superior performance in both predictive accuracy and uncertainty estimation compared to existing methods.
summary: This paper presents the Positional Encoder Graph Quantile Neural Network (PE-GQNN), a novel model for spatial data prediction. PE-GQNN integrates Positional Encoder Graph Neural Networks (PE-GNN) with Quantile Neural Networks (QNN) to improve predictive accuracy and quantify uncertainty. Experiments on real-world datasets show that PE-GQNN consistently outperforms existing methods in predictive accuracy and uncertainty calibration.
weaknesses: - The proposed method is practically useful, but the way PE-GNN and QNN are combined is somewhat straightforward.
- There are several innovations in the architecture, but they are all practical techniques, not theoretically sophisticated.
- There is no discussion of computational cost for high-dimensional data.
- There is no discussion of the shortcomings of the proposed method.

---

confidence: 4
contribution: 2
initial_rating: 1
presentation: 2
questions: 1. The Figure 4(a) does not have legend. What does the 9 curves represent? It is unclear to me how does these 9 curves represent 10 samples.
2. Instead of using MADECP, what is the performance when using popular adopt interval? For example, 95% confidence interval?
3. The real-world geographic data usually collect in a temporal manner, so how does the data is separated into different train/val/test segments? It leads to data leakage issue for sure if the separation does not only use historical data for training.
4. The sensor of the geographic data also will be removed/replaced/re-deployed/defected from time to time in real world applications. Does the proposed method also robust in these cases?
soundness: 2
strengths: 1. The framework is easy to understand and follow.
2. The final experimental results indicate some effectiveness of the proposed method.
summary: The paper introduced a distribution-free uncertainty quantification framework by integrating PE-GNNs, Quantile Neural Networks for geographic data.
weaknesses: 1. The term "fully nonparametric" could be misleading. It would be more accurate to use "distribution-free."

2. The novelty of the proposed method is limited. Quantile regression is already proposed and widely used in distribution-free uncertainty quantification. For example:
   - "Single-model uncertainties for deep learning." Advances in Neural Information Processing Systems 32 (2019).
   - "Image-to-image regression with distribution-free uncertainty quantification and applications in imaging." International Conference on Machine Learning. PMLR, 2022.

   Additionally, using KNN-Graph with geographic coordinates has been proposed and widely applied in different domains:
   - "Dynamic graph CNN for learning on point clouds." ACM Transactions on Graphics (TOG) 38.5 (2019): 1-12.
   -  "Spatiotemporal graph convolutional networks for earthquake source characterization." Journal of Geophysical Research: Solid Earth 127.11 (2022): e2022JB024401.

3. PE-GNN does not claim the method as SOTA. So which baseline are considered as SOTA in this submission? PE-GNN is better than GNN but not the SOTA. At least there are more available models in applying GNN for earthquake data and traffic data.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: What does ST stand for in the last sentence of the 1st paragraph of Section 2?

In Figure 4(a), why does the predicted density of PE-GQSAGE look like Gaussian distribution? The proposed method outputs predicted quantiles at the given tau. So, the prediction by the proposed method does not necessarily to be Gaussian shape.
soundness: 2
strengths: Proposed some techniques for quantifying uncertainty in spatial regression.
Experiments with three datasets.
summary: This paper proposes a graph neural network-based method for uncertainty quantification in spatial regression.
weaknesses: The technical contribution of this paper is limited since the proposed method is a combination of Positional Encoder Graph Neural Networks and quantile regression. The proposed method introduces some techniques, such as the use of response variables y and quantile parameters tau in neural networks. However, the novelty of these techniques is incremental.

The experimental results are not convincing. There have been proposed many quantile regression methods in neural networks; i.e., outputting variance in the last layer, the use of sparse Gaussian processes at the last layer. These methods can be easily combined with graph neural networks. The comparison with such existing methods can demonstrate the effectiveness of the proposed method.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1.  It is better to add some simple examples to illustrate the integration process and novelty in the process on how to design this new graph. 

2. The paper shall demonstrate the complexity change for this new PE-GQNN graph.
soundness: 3
strengths: There are several strengths demonstrated in the paper:
1. The paper introduces the Positional Encoder Graph Quantile Neural Network (PE-GQNN), a new approach that integrates PE-GNNs, Quantile Neural Networks, and recalibration techniques in a fully nonparametric framework, requiring minimal assumptions about the predictive distributions.
2. The paper has demonstrated the results on three datasets: California Housing, Air Temperature, and 3Droad with 6 different approaches including the proposed PE-GQNN.
summary: In this paper, the authors proposed the Positional Encoder Graph Quantile Neural Network (PE-GQNN) as a new framework to enhance predictive modeling for geographic data.  The major contributions of this paper are listed as the following: The empirical results showed the capability of PEGQNN to achieve lower MSE, MAE, and MPE compared to traditional GNN and PE-GNN. Also, PE-GQNN demonstrated substantial improvements in predictive accuracy and uncertainty quantification.
weaknesses: The weaknesses of this paper are listed as the following:

1.	The innovation of this paper seems incremental. Positional Encoder Graph Quantile Neural Network (PE-GQNN) is just a simple combination of PE-GNN with Quantile regression model.

First, the paper shall illustrate in detail the challenges in the integration process. Normally a good integration will include some short cuts to reduce the total cost while comparing with the cost of simple addition of several algorithms together directly. Please try to add some "novel points" or new ideas to demonstrate your merits in integration.

Second, it is better to add some simple examples to illustrate the integration process and novelty in the process on how to design this new graph. Please illustrate which cost you saved compared the cost that you simply integrate several different phases from literature.

Third, will the new integrated framework achieve higher performance compared simply by adding several phases together? What other advantages do you have for the new framework?


2.	The paper shall demonstrate the complexity change for this new PE-GQNN graph. 

Is this new graph a simple integration? 
How much is the increase of the total complexity? 
Is there any new approach to lower the total complexity?

Please provide a specific complexity analysis comparing PE-GQNN to PE-GNN, including time and space complexity.


3.	For the experiments, the following should be addressed. 

First, the paper presented experimental results obtained from three datasets: California Housing, Air Temperature, and 3Droad. It seems that the paper lacks discussion about whether other kinds of datasets are suitable for this new approach. Also, how much the total cost changes to implement this new approach.

Second, please discuss the generalizability of their approach beyond the three datasets used. Also, please address what characteristics of a dataset make it suitable for PE-GQNN.

Third, additionally, it would provide valuable practical insights if the authors can demonstrate a comparison of implementation costs between PE-GQNN and existing methods.