Positional Encoder Graph Quantile Neural Networks for Geographic Data
Authors: William E. R. de Amorim
Scott A Sisson
Thais Carvalho Valadares Rodrigues
David J Nott
Guilherme S. Rodrigues
Venue: Accepted by TMLR ()

Abstract: Positional Encoder Graph Neural Networks (PE-GNNs) are among the most effective models for learning from continuous spatial data. However, their predictive distributions are often poorly calibrated, limiting their utility in applications that require reliable uncertainty quantification. We propose the Positional Encoder Graph Quantile Neural Network (PE-GQNN), a novel framework that combines PE-GNNs with Quantile Neural Networks, partially monotonic neural blocks, and post-hoc recalibration techniques. The PE-GQNN enables flexible and robust conditional density estimation with minimal assumptions about the target distribution, and it extends naturally to tasks beyond spatial data. Empirical results on benchmark datasets show that the PE-GQNN outperforms existing methods in both predictive accuracy and uncertainty quantification, without incurring additional computational cost. We also identify important special cases arising from our formulation, including the PE-GNN.

Published in Transactions on Machine Learning Research (12/2025)
Positional Encoder Graph Quantile Neural Networks for Geographic Data
William E. R. de Amorim*
wilkam walkom_.appe@hhotmal.oom
Deportment of Statistics
University of Brasilia Brazil
Scott A. Sisson
scoll.sisson Gwnsw. eduu
School of Mathematics and Statistics University of New South Wales Austratia
Thais C. V. Rodrigues
thaisrodrigues Wumbb. br
Department of Statistics University of Brasilia Brazil
David J. Nott Department of Statistics and Data Science National University of Singapore
standj@russ wwt.rrss.suu...... edu.sg sg
Singapore
Guilherme S. Rodrigues
guilhernnerudrigucss Gunb. br
Departmend of Statistics University of Brasília Brazil
Reviewed OpenReview:
Abstract
Positional Encoder Graph Neural Networks (PE-GNNs) are among the most effective mod- els for learning from continuous spatial data. However, their predictive distributions are often poorly calibrated, limiting their utility in applications that require reliable uncer- tainty quantification. We propose the Positional Encoder Graph Quantile Neural Network (PE-GQNN), a novel framework that combines PE-GNNs with Quantile Noural Networks, partially monotonie neural blocks, and post-hoe recalibration techniques. The PE-GQNN enables flexible and robust conditional density estimation with minimal assumptions about the target distribution, and it extends naturally to tasks beyond spatial data. Empirical re- sults on benchmark datasets show that the PE-GQNN outperforms existing methods in both predictive accuracy and uncertainty quantification, without incurring additional computa- tional cost. We also identify important special cases arising from our formulation, including the PE-GNN.
*Corresponding author.

---

Published in Transactions on Machine Learning Research (12/2025)
1 Introduction
Large spatial datasets are naturally generated in wide range of applications in economics (Anselin 2022
meteorology (Bi et al,, 2023). urban transportation (Lv et al. 2014 Derrow-Pinion et al. 2021 Kashyap
et al. 2022), social networks (Xu el al.) 2020), e-commerce (Sreenivasa & Nirmala 2019) and other fields. Gaussian Processes (GPs) [Rasmussen & Williams 2006) Cressie & Wikle, 2011) are 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²) ) and storage complexity O(n²), naive GP methods quickly become intractable for large datasets (Liu et al, 2019). This has led to a large range of approximate inference methods, such as those based on sparse approximations to covariance approximations (Cressie et. al. 2022) or nearest neighbour approximations (Vecchia 1998 Dalta et al. 2016
precision matrices (Furrer et al.) 2006 Lindgren et al., 2011], low rank
Katzfuss & Guinness p021p.
Given the difficulty of GP computations, it is of interest to explore scalable methods for large spatial datasets using neural networks (NNs) and to enhance their ability to quantify uncertainty. A prominent method for making spatial predictions using Graph Neural Networks (GNNs) the Positional Encoder Graph Neural Network (PE-GNN) Klemmner et al. (2023). Our contribution is to make key modifications to the PE-GNN architecture to enhance its ability to make accurate spatial predictions and to better quantify uncertainty. These modifications will be explained further below.
NNs are popular in data modeling and prediction tasks like computer vision and natural language process- ing (NLP). However, traditional NNs struggle to handle spatial dynamics or graph-based data effectively. GNNs KKipf & Welling) 2017 Veličković et al. 2018: Hamilton et al., 2017) offer a powerful and scalable method for applying NNs to graph-structured data. The idea to share information through the edges of graph, allowing nodes to exchange information during learning. GNNs are versatile and can uncover non- linear 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 Iraditional GNNs often struggle to model complex spatial relationships. To address this, Klemmer et al. (2023) introduced the PE-GNN, which enhances predictive performance in spatial interpolation and regression. However, the 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 (for more detailed discussion of calibration definitions, see Kock ct al.) 2024). Recently, Bao et al. (2024) proposed a new framework called Spatial Multi-Autention Conditional Neural Processes (SMACNPs) for spatial small sample prediction tasks. SMACNPss use GPs parameterized by NNs to predict the target variable distribution, which enables procisc predictions while quantifying the uncertainty of these predictions. However, these methods remain constrained to producing Gaussian predictive distributions, limiting their capacity to represent asymmetric or multimodal random variates.
Methods based on quantile regression are an alternative approach to probabilistic forecasting making rapid progress in recent years. Tagasovska & Lopez- Paz 2019 introduced the Simultaneous Quantile Regression (SQR) loss function that we use in our formulation. Si et al. (2022) proposed a novel architecture for estimating generic quantiles of a conditional distribution. In one dimension, this method produces a quantile function regression ((T)|z that estimates the r-th quantile of the predictive distribution of the target variable given a feature vector ac Kuleshov & Deshpande (2022) argue that the method of et al. (2022) 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. dow-dinensi The first model outputs features, usually summary statistics like quantiles, representing 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 two sets of contributions. (1) We propose a new architocture that merges the two-step procedure of Kuleshov & Deshpande (2022) into a single model by postponing the concatenation of the 7 value used by Si et al. (2022). Additionally, 7 is transformed beforehand to facilitate learning and provide
2

---

Published in Transactions on Machine Learning Research (12/2025)
better control over the form of the predictive distribution. To ensure the model outputs valid, non-crossing quantile functions, we employ the partially monotonic blocks introduced by Nolte et. al. 2023). This approach enhances the network's s ability to model uncertainty and the model becomes more robust to high-dimensional prodictor spaces. (2) We introduce structural changes to the 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. In addition, we introduce the mcan target value of node's neighbours as a feature after the GNN layers, closer to the output.
The structure of this work is as follows: Section 2 offers a brief background overview, Section 3 outlines the proposed method for geographic data prediction, Section shows experimental results on four real-world datasets, and Section concludes.
2 Background
Graph Neural Networks (GNNs) are powerful and scalable tools for represen- tation learning and inference on graph-structured data. They exploit the topological relationships between adjacent nodes to produce context-aware embeddings, which can be effectively used in downstream tasks [Wu et al. 2022). GNN layers iteralively refine each node's embedding by aggregating information from its own features as well as those of its neighboring nodes.
Graph Neural Networks:
Graph Convolutional Networks (GCNs) Kipf & Welling 2017 are a specific type of GNN layer inspired by the convolution operations used in Convolutional Neural Networks (CNNs). For weighted graphs, a GCN layer updates node embeddings according to the following equation:
H (()) = f(h)
[A 1=
for k E {I....,K} K}.
(1)
In this formulation, the input to the network is the feature matrix H = X. The function f() denotes a non-linear activation function (e.g., ReLU), and ((k) is a learnable weight matrix. The adjacency matrix A encodes the graph structure, with edge weights typically derived from node distances and zero entries for unconnected pairs. The identity matrix adds self-loops, and D is the corresponding degree matrix. In our experiments (Section 4), we also consider other widely used GNN architectures, namely, Graph Attention Networks (GATs) /Veličković et al. 2018) and GraphSAGE (Hamilton al. 2017).
Positional Encoder Graph Neural Network: In typical spatial regression setting, each datapoint is represented as = {y: xi c where y, is a continuous scalar target variable, æ; is a vector of input features, and C₁ denotes the geographical coordinates associated with observation i. given batch of datapoints B = (. | can be fully represented by three matrices: the larget vector yn E R=ax1, the feature matrix XB xp and the coordinate matrix CB E Rnax2 respectively.
Klemmer et al. (2023) introduced a novel approach for incorporating spatial information into GNNs: the Positional Encoder Graph Neural Network (PF-GNN). In this framework, a positional encoder (PF) processes the spatial coordinate matrix CB to produce a Iearned spatial embedding matrix C.mab, The embedding is computed as Cgub = PE(CB, Tmin, maax, OPF) = (N(SS(CR, Tmin mman Opp), where ST denotes a deterministic set of sinusoidal transformations with hyperparameters Jmin and maax: and NN is a fully connected neural network with trainable parameters Org- A complete description of the PE mechanism is provided in Appendix A.2
The matrix Clmmb is concatenated column-wise with the node features before the application of GNN layers. Thus, the input to the first GNN layer is given by B = concat(X,,, n. Cgmhh). At each training step, a new random batch of nodes B is sampled, and the full pipeline - graph construction, spatial embedding gencration, feature concatenation, and GNN propagation For each node pt E (Pi,...,,P') the PE-GNN predicts a target value in and, as an auxiliary task (Klemmer & Neill 2021). the corresponding Local Moran's I statistic (Anselin 1995), denoted by [(y)). The total loss used by Klemmer et al. (2023) combines both objectives:
exccuted using only the nodes in that batch.
LB MSE, yn) + A XMSE(f(U)), ((y)))
3

---

Published in Transactions on Machine Learning Research (12/2025)
where À controls the contribution of the auxiliary task.
Quantile regression: Koenker & Bassett Jr| (1978) proposed a linear quantile regression model to estimate conditional distribution quantiles.
uses the pinball loss pr(ri) = max (Tri, (T - 1)r), where Tt = - [T),,, (T( = x, and t is the desired cumulative probability associated with the predicted quantile (T).). The pinball loss for the i-th observation is pr(r)). The loss over a dataset is the average pr value over all datapoints. A natural extension of quantile linear regression is quantile neural networks (QNNs). This approach is illustrated in Figure [aa] which seeks estimate the conditional quantiles for a pre-defined grid (7' (r',...,, =4). Each quantile is estimated by an independent model (Figure 4a). This lead to quantile predictions with quantile crossing (e-g., a median prediction lower than the first quartile prediction). Rodrigues & Pereira (2020) proposed an approach that outputs multiple predictions: one for the expectation and one for each quantile of interest. The loss function is:
1 +
a
Pet (u: (((T)))
f
MSE(O.) y)
(2)
n
(0199) proposed a method to generate model that is independent of quantile selection. During training, for each datapoint in the batch, a Monte Carlo sample T ~ (0, 1) is drawn and concatenated with the corrosponding datapoint feature vector. The SQR loss function is similar to Eqn. 2 but they predict random quantiles L = 1/ Σx= (,, - (T:)), with {m}] ~ ((, 1). As the network learns, it becomes able to provide direct estimate to any quantile of interest. Hence, this procedure outputs an inherently calibrated model suitable for conditional density estimation. Si ct al. 2022 construct NNs with
Tagasovska & Lopez- Paz
similar loss function (Figure 4cp.
Kuleshov & Deshpande 2022 adapted the architecture from Si et al 2022) into a two-step process for larger prodictor spaces (Figure 4d). 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 т. During inference, the recalibrator takes the new features and an arbitrary т as imputs to produce the quantile prediction. This method is highly dependent on the choice of recalibrator features.
3 Method
In this work, we introduce the Positional Encoder Graph Quantile Neural Network (PE-GQNN), a novel framework for predictive modeling on spatial data. Algorithm shows the step-by-step procedure to train a PE-GQNN model.
Figure I illustrates its complete pipeline. Here, each rectangle labeled "GNN', 'LINEAR' and 'MONO- TONIC" represents a set of one more neural network layers, with the type of each layer defined by the title inside the rectangle. At each layer, nonlinear transformation (e.g. ReLU) may be applied.
After initializing the model and byperparameters, the first step of PE-GQNN is to randomly sample a batch B of datapoints. The next step projects the matrix of geographical coordinates CB into the positional embeddings, Cgmb (nzxu) (Algorithm [ Step 5). CB also used to compute the distance between each pair of datapoints (Step 6j. From these distances and 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 Ag-
At Step 13 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 XB. One more fully-connected layers are then used to reduce the feature embedding dimensionality. This process receives the constructed graph and the batch feature matrix
X (mgx))) as input and yiclds an embedding matrix of features as output: Xemb (naxg)' Step 14 performs a column concatenation between the feature embedding Xgmb (naxg) and the output obtained from the PE: Cemb (Mgxu) x This concatenation results in the matrix L((nx((+xx)-
4

---

Published in Transactions on Machine Learning Research (12/2025)
Algorithm PE-GQNN training
Require:
Training data target, features, and coordinates matrices: M(xx)), X (mxp), and (nx2)).
positive integer k defining the number of neighbors considered in the spatial graph. Positive integers tsteps and BB, the number of training steps and the batch size. Positive integers u, 9, and s, the embedding dimensions considered in, respectively, the PE, the GNN layers, and the layer where we introduce and y. An activation function f() for T.
Ensure:
A set of learned weights for the model initialized at Step
1: Initialize model with random weights and hyperparameters.
2: Set optimizer with hyperparameters.
for b- 1 to tsteps do
Batched training
4: Sample minibatch B of un datapoints: X X((nxx)) CB (B
Input x2) into PE, which outputs the batch's spatial embedding matrix C (ngXu))
Compute the great-circle distance between each pair of datapoints from Cp.
Construct a graph using k-nearest neighbors from the distances computed in Step
8: Set AB as the adjacency matrix of the graph constructed in Step 7
for i + to ng do Using An, compute Tt = Eÿ= Mj. where j = 1,...,k are the neighbors of i.
10:
11:
end for
12:
Sct YB = .y...... Dnn
13:
Apply GNN layers to the features X X xp) followed by fully-connected layers to reduce dimension-
ality. This step outputs a feature embedding matrix Xemb
(naxg)'
14:
Column concatenate Xg (# with Cgmb (mgxu)² which results in LB(nxX(--u))-
15:
Apply fully-connected layers to reduce L ((++a) to (B(n))))
16:
Create vector with values sampled from (0, 1): (B = [r....,, Tan)
17:
Column concatenate ΦB with f (TB) and ys to create (B(nnx(*+2))
18:
Predict the target quantile vector (TI)...... Ing (Tna))] using ΦB. Compute loss LB ag Σn Pr: (v,- (T:)).
19:
20:
Update the parameters of the model using stochastic gradient descent.
21:
end for
Subsequently, we use one or more fully-connected layers (Step 15 reduce the dimensionality Ln, making it suitable for two innovations in PE-GQNN. This set of fully-connected layers outputs the matrix ( xs) which is then combined with ÜB and тв- ÜB represents vector with 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 (Step 10), and has dimensions n. x 1. It. is comparable to a vector of predictions generated by a KNN regression model, where neighbours are determined using the distance calculated from goographical 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 Ag. We introduce this input at a later stage to avoid data leakage. If the GNN operator received ye 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 al. 2020
In the same layer where YB is introduced, we apply similar approach to Si et al. 2022 to make PE- GQNN an inherently calibrated model suitable for probabilistic and quantile predictions. For each batch B, we create a ng X 1 vector (mgx)) = _
random U(0, 1) draws (Step 16). Then, we column concatenate ΦB with (TB) and DB create (=(xx()=)))) (Step 17. where f() is an activation function. Next, forward propagation is computed (Step 18 in one or more fully-connected layers, outputting predicted
5

---

Published in Transactions on Machine Learning Research (12/2025)
Yp.
GNN LINEAR
GNN
X
O
bn
MSE yn)
Construct KNNN-Graph with Cg geographic coordinates
Apply GNN model
Compute MSC loss
Randomiy sample fly datapoints
PE-GNN
LINEAR LINEAR
y
Construct Graph with Cy geagraphic coordinates
( (
concan m
GNN
Concatenate embadding with festures
Randomly sample "ly datapoints
2 2
Generrle bedlding using PE applied to Cp
Applyy GNN mode
Compute weighted MSE loss
(Kg. Yp)
X
LLIEER
Xami
32
0(0.1)
PE-GQNN
a
concat (X
Xemb cF
concet go Foi (F,)))
Corotruct SNN-Gruphrwill
(e(s))
Ca geegraphic coordinates
layers
Apply GNN
Ln
Yo ((F.))
MONO
La
Pr. ((())
fectures
Randomly sample Ry datapoints
LINVEN
Concatenate embeddingg
Reduce
Concatenate 4.. f(T)
Apply leyers
Compute
Ganerate embadding PE appliad to Co.
Figure PE-GQNN compared to PE-GNN and GNN
quantiles for each datapoint in the batch. These predictions are then used compute the SQR loss function introduced by (Tagasovska & Lopez-Paz (2019).
Incorporating values into the model architecture improves its ability to model uncertainty and serves as regularization mechanism (Rodrigues & Percira 2020). The use of SQR loss acts a natural regularizer, producing & detailed description of the predictive density beyond just mean and variance estimation. For predictions, the quantile interest, т, must be given, along with the basic data components (e.g. == 0.25 gives the first quartile). If interest is predicting multiple quantiles for the same observation, the input can be propagated up to the layer where introduced. For each quantile of interest, propagation can be limited to the final layers.
Homoscedastie Gaussian model as particular case: If in Steps and 18 one sets ʃ(T) = br), where 4() denotes the PDF of a standard Gaussian distribution, and use single linear predictive layer, the modol's quantile function,
(.(T) (b Jg + Σ w Φj) w
(r) #4 oT - - Vi I....,n,,
(3)
matches the quantile function of a Gaussian random variable, therefore reducing PE-GQNN a ho- moscedastic Gaussian regression model with mean μ4 and a common, learnable standard deviation wr. Here, b, Wy, Wg. and {w,} are the prediction layer parameters, and { are the activation values from the previous layer. This result may appear counterintuitive, as the SQR loss can produce a model strue- turally identical to one trained with a fundamentally different objective, namely, the MSE loss. Similarly, if ((r) = logit(7) and a single linear predictive layer is employed, the model yields logistic predictive distribu- tions, a continuous distribution that should not be confused with logistic regression. Therefore, a well-chosen activation function f()) not only facilitates learning but also explicitly shapes the structure of the predictive distributions. The use of a set of monotonic layers enhances the flexibility of our approach, significantly reducing the reliance the choice of f(T).

---

Published in Transactions on Machine Learning Research (12/2025)
Quantile crossing and Monotonic Blocks:
Several related approaches (Tagasovska & Lopez- Paz,
2019 Rodrigues & Pereira 2020 Si et al., 2022 Kuleshov & Deshpande, 2022) are subject to quantile
crossing, a phenomenon that occurs when the requirement that higher quantiles be greater than or equal to lower quantiles is violated. To ensure the network outputs valid probability distributions (i.c. ((T, S i('', Vr VT < T'), we propose using the Lipschitz Monotonic Networks (LMN) introduced by Nolte et al. (023)) to approximate q(T, Φ, 3). Those are highly expressive blocks of layers that can approximate all monotonic Lipschitz bounded functions. That is, any function for which there exists constant À such that @(T, v) ~(T', Φ, f) - A/T r'I. vO<T,T' < 1.
By design, there a trade-off wetn model expressiveness and probabilistic coherence. The decision to include monotonic layers should be guided by the specific application. In practice, it is straightforward to assess on a validation set whether quantile crossing is significant. When it not (especially since SQR (Tagasovska & Lopez-Pazz 2019) already mitigates this issue), adding monotonic layers may be unneces- sary even counterproductive. As usual, we recommend exploring multiple network configurations during hyperparameter tuning to achieve good balance between expressiveness and coherence.
Number of Monte Carlo samples: When applying the framework proposed by Si et al. (022)) we followed Tagasovska & Lopez-Paz (2019) and set d = for the T values. Let L(0, T, x. y) denote the loss function for a given quantile T ~ U 1) and an observed pair (x, y) ~ Daatar where Daata represents the full data-generating process. At cach training iteration, we minimize the empirical loss LB, which, by the Law of Large Numbers, converges to Σ())- Er, æ, yC(0, T, x. y) the batch size n: - xc. Consequently: the gradient estimates converge to the same limit for any d. Since the batch size in our experiments is sufficiently large (sce Section 4.2 each gradient estimate effectively aggregates losses over more than a thousand distinct quantile levels. Thus, setting (d = 1) simplifics the implementation without compromising performance, as remonstrated in Section 4
The final layer should preferably 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 (Goodfellow et al. 2016).
Target domain:
4 Experiments
4.1
Experimental setup
The PE-GQNN was implemented using PyTorch Paszke et al., 2019) and PyTorch Geometric (Fey & Lenssen 2019). The source code is available at: https://gitii com http//111aaaaapp1e///ppp.r GQNN We conducted comprehensive simulations to explore the prediction performance and other properties of the Intel i7 (7th Generation) processor.
proposed model. Computation performed
The experiment was designed to compare five primary approaches for addressing spatial regression problems across four diverse real-world datasets (California Housing, Air Temperature, 3D Road, and Australian Census). See Table Table lists each candidate model and their applicable datasets. All models were trained using the Adam optimizer (Kingma Ba 2015), with carly stopping employod prevent overfitting. For GNN-based approaches, we used k = 5 nearest neighbors to construct the graphs. The learning rate set to 0.001 across all models. A batch size of 1,024 was used for the Air Temperature dataset, while a batch size of 2,048 was used for the remaining three datascts. The architectural details are given in Appendix Section A.3
Candidate models:
Approach I involves the traditional application of GNNs to geographic data. Three types GNN layers were considered: GCNs (Kipf & Welling 2017), GATs (Velickovi al., 2018), and GSAGE (Hamilton et al., 2017). For each of these, the architecture remains consistent facilitate performance comparisons: two CCN/GAT/GSAGE layers with ReLU activation and dropout, followed by linear prediction layer.
7

---

Published in Transactions on Machine Learning Research (12/2025)
Table 1: Summary of candidate models.
Approach
Model
Type
PE
Innovations
Loss
Datasets
I
GNN
GNN
No
No
MSE
All
II
PE-GNN = best
GNN
Yes
No
MSE (M)
All
III
PE-GNN (with SQR)
GNN
Yes
Partial
SQR
California
IV
PE-GQNN
GNN
Yes
Yes
SQR
All California and Air Temp.
V
SMACNP
GP
No
No
Log Likelihood
Approach II involves the application of PE-GNN Kleemeer et al., 2023) with optimal weights for each dataset and layer type combination, demonstrated by the experimental findings of Klemmer et al. (2023) The GNN architecture used is the same as for approach I. The code is available at githuty
Approach III represents a naive combination of the PE-GNN with the quantile regression framework described in Section 3 Spocifically, we trained the PE-GNN with the SQR loss function, concatenating f(r) = 7 immediately after the GNN layers. Approach IV, which the primary focus of this research, is the PE- GQNN. Compared to Approach III, it introduces the following innovations: T is incorporated into the network through probit(T), but only after reducing LB into ΦB. The GNN layers no longer process the positional encoders, and the final fully connected layers are replaced with monotonic blocks. Additionally, is used as a feature in the final part of the network. The architectures of the PE and GNN layers remain identical to those in the previous approaches.
Finally, a benchmark approach that does not use GNNs but was recently proposed for modelling spatial data will be considered as approach V: SMACNPs. This approach, proposed by [ao et al. (2024), has demon- strated superior predictive performance, surpassing GPs models in the three real-world datasets considered. This model was implemented following the specifications of [ao et al. 20024) with code available at githul
Approaches I and 1l do not inherently provide predicted conditional distributions. However, they optimize the MSE metric, they implicitly learn a Maximum Likelihood Estimate (MLE) of 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 ÿB for each batch, we pre-calculated ý 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, report the SQR loss and calibration metric introduced by Kuleshov et al. (2018): calibration = E For quantile predictions of a calibrated model for a given T, the proportion of observed values less than or equal to the predicted quantile should approximate т. While informative, the calibration metric should be interpreted with caution: (i) it represents coarse global average, meaning that regional over- and underestinations can cancel each her out; and (ii) it is highly sensitive to random variation, as evidenced by the large standard deviations reported in Table 2 Therefore, we regard SQR as the more reliable metric for assessing uncertainty quantification.
( - 1 [y: S (r)
4.2 California Housing
This dataset comprises pricing information for 20,640 census block groups in California, recorded during the 1990 U.S. census (Pace & Barry 1997)- The main objective a regression task: predict housing prices, y, through the incorporation of six predictive features, æ, and geographical coordinates, 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. 2024 a training subsample was extracted to represent 10% of the entire dataset.
The
httte//gghhaaacc////mtttttttttte gnn
8

---

Published in Transactions on Machine Learning Research (12/2025)
Table Performance metrics on the California Housing test set. Values are reported as mean + f0.975.9 x SD/VI0, yiclding approximate 95% confidence intervals from the t-distribution with nine degrees of freedom. All metrics except # Param. and Time are multiplied by 100. Boldfaced entries represent the best mean value. The training Time is measured in minutes and Cover. refers to the observed coverage for 95% confidence intervals.
Model
Param.
Time
MSE
MAE
SQR
Cover. (%)
Calibr.
GON
1,313
40.53
18:77
2.42
11.39
1.0.5
34.60
PE-GON best
24.129
17.48
2.81
1.75
+-0.6
9.34 +0.10
3.49
93.02
+.555
40.06
PE-GON (with SQR)
25,217
12.64 1.00
1.79
0.00
9.20 =0.18
3.38 0.07
85.70
0.88
18.08
6.02
PE-GQUN (ours)
26,201
13.86
0.82
1.11
+ 0.02
6.83 0.10
2.61
= 0.06
91.47
+ 1.31
19.02
= 17.94
GAT
1,441
31.93
10.51
2.40
0:1.5
11.32 0.34
4.17 0.1.
92.95
+ 0.53
33.72
= 6.85
PE-GAT = best
24,290
26.62
4.38
1.77
+ 0.08
9.36 0015
3.51 0.07
92.98
0.47
11.32
7.69
PE-GAT (with SQR)
25,345
15.03
2.72
1.84
85.31
19.85
5,06
PE-GQAT (ours)
26,329
14.39
1,15
1.11
0.08
6.77 0.30
2.60 0.06
95.24 1: 0.80
11.51
6.76
GSAGE
2,529
48.33
5.07
1.62 + 0.06
9.21 0.20
3.40 0.08
93.88
0.31
37.55 = 5.84
PE-GSAGE A= best
27,426
29.37 = 4.70
1.07 1.08
7.21 0.18
2.73 4:4.0
93.94
0:48
39.83
2: 9.23
PE-GS GSAGE (with SQR)
28,481
16.26
2.48
1.14
1.1.2
7:23 0.30
2.66 1.0.12
85.50
1: 1.87
28.99
13:97
PE-GQSAGE (nurs)
27.417
1.20
0.00
2.20 +-005
1.04
17.00
SMACNP
748,482
73.46
7.49
2.35 1:
0.58
1.077
03
As shown in Table 2 PE-GQNN delivers strong performance across all evaluation metries, significantly outperforming traditional GNNs, PE-GNN, and SMACNP. Among models with GSAGE layers the overall top performers - PE- GQSAGE achieved the lowest values for Time, MSE, MAE, SQR, coverage and calibration error. Specifically, compared to PE-GS GSAG,, it reduced the training time by 49%, MSE by 21%, MAE by 18,, SQR by 17%, and calibration error by 63%. In contrast, a naive combination of PE-GNN with SQR offered no meaningful improvement over PE-GNN alone.
Figure 2 presents plots that elucidate the behavior of the PE-GQSAGE predictions. The validation MSE curves for given are shown in Figure 2a Panel 2b illustrates the predicted density of subsample of three observations from the test set. Parametric models typically assunee fixed output structure such Gaussian distribution - which can restrict their expressiveness. In contrast, PE-GQNN imposes minimal assumptions on the form of the predictive distribution, offering greater flexibility. As illustrated in Figure 26 PE-GQNN capable of producing predictive distributions with varying shapes and scales, whereas PE-GNN can only differentiate samples based on their location (mean value). is worth noting that, this example, assuming Gaussian distributions is reasonable. However, in scenarios involving asymmetric or multimodal distributions, the flexibility of model-free approach becomes even more advantageous.
0.06
DSAGE
D&4
PE PE-ORAGELAA A bee DEB.RRRR (MAN)
1.0
GLAGE
0.05
Observed Confidence Level
PE-CORACE
Grid Stnwdert
Prodicted Denaltty
0.6
0.03
0.4
0.02
0.01
0.0
0,00
0.2
0.8
1.0
Epoch
y
Expested Confidence Level
(a) Validation MSE.
(b) Predictive densities.
(c) Calibration Plot.
Figure 2: Diagnostics for the California Housing dataset. All plots refer a single run of each model. (a)
Validation MSE curves. (b) PE-GQSAGE predicted densities of 3 observations of the test set. (c) Calibration Plot.
Figure 2e displays the observed confidence level (averaged across ten runs) for the test set quantile predictions using each of the 99 values in [0.01, 0.02. 0..9 This type of plot was proposed by Kuleshov al.
9

---

Published in Transactions on Machine Learning Research (12/2025)
2018). The closer model gets to the dashed diagonal line, the closer the expected and observed confidence levels are. The Gold Standard represents one Monte Carlo draw from perfectly specified model, where for each quantile level, the observed confidence level is the observed success rate in n Bernoulli trials with a success probability of r, where n is the number of test set instances. It is evident that PE-GQSAGE has the best calibration performance. This is particularly notable when compared to SMACNP, which exhibits substantial calibration deficiencies due its tendency to overestimate the variance component.
As summarized in Figure 1 both PE-GNN and PE-GQNN reconstruct the graph for each mini-batch. According to Klemmer et al. (2023). "This approach brings unique advantage: When training with (randomly shuffied) batches, points may have different neighbors in different training iterations (...) This forces PE to learn generalizable features." While this finding is generally supported by their experimental analysis, the additional stochasticity introduced by graph reconstruction can degrade perfor- mance when hani-b are too small, as shown in Table & Conversely, using very large mini tini-batches can lead to a substantial increase in computational time. As usual, it is good practice tune this and other hyperparameters using efficient optimization framework.
Mini-batch size:
Table 3: Performance metrics on the California Housing test set for different mini-batch sizes. Values are reported as mean + f0.975,9 x SD/V10. All metrics except # Param. and Time are multiplied by 100.
Batch Size
# Param.
Time
MSE
MAE
SQR
Cover. (%)
Calibr.
128
27.417
4.73 = 0.45
0.96 0.08
6.47 =
2.46
0.10
95.55
0.96
83.45 = 54.06
250
27,417
4.76 = 0.33
0.90 0.03
6.18 0.10
2.35
+ 0.05
95.83
0.70
36.71 = 16.88
512
27.417
6.45 = 0.54
0.89 0.04
6.13 2014
2.31 0.07
95.54
0.82
39.77 = 18.61
1024
27.117
9.43 = 1.02
0.86 0.02
6.02 20.1
2.29
0.03
95.17 0.70
20.37
= 22.54
2048
27,417
14.99 1.25
(1855 + 0.02
0.08
2.26
+ 0.0.5
94.46 1.04
14.43 17:10
5096
27,417
30.01
= 2.11
0.84::0.0
5.87 0.09
2.23
0.05
94.80
1.18
8.41 6.18
4.3 All datasets
Experiments were also conducted on three additional geographic datasets, two of which were previously used by Klemmer et al. (2023) and Bao et al. (2024). See Table 4 The Air Temperature dataset (Hooker et al. (0018)) contains geographical coordinates for 3,076 meteorological stations worldwide, with the goal of predicting mean temperatures (y) using mean precipitation levels (z). 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 [ao et al. (2024). The 3D Road dataset (Kaul et al. (2013)), includes 434,873 points with latitude, longitude, and altitude for the Jutland, Demmark road network. The task is to interpolate altitude (y) using latitude and longitude (c). The data were split into 90% for training, 1% for validation, and 9% for testing. The 2021 Australian Census data were used to predict median total personal income (AUD/ ( across Statistical Areas Level 1 (SA1) based other census variables. The "2021 General Community Profile for SA1" dataset is publicly available and spatial coordinates were obtained from the absmapsdata R package (Mackey) (025). Tables containing additional income information (numbers 17, 32, 33, 38, 40, 57, 58, and 59), as well as columns with more than 90% zero values, were excluded. The resulting design matrix includes 5,978 features, providing a suitable benchmark for assessing model performance in high-dimensional scttings. The dataset was split into 80/10/105 training, validation, and test sets.
Table 5 showcases the experimental results obtained from all four datasets. The PE-GQNN models incor- porate all innovations discussed in Section 3 PE-GQN FE-GONN outperforms both traditional GNN and PE-GNN. Across all datasets, the innovations introduced by PE-GQNN result in substantial reductions in MSE, MAE, and SQR. 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 datasct, SMACNP achieves the lowest MSE and MAE but suffers from significantly uncalibrated predictions, reflected by a much higher SQR and calibration error compared to PE-GQNN.
10

---

Published in Transactions on Machine Learning Research (12/2025)
Table 4: Summary of benchmark datasets. Column Peatures refers to the number columns of the Design
Matrix (not including the geographical features Latitude and Longitude). The remaining columns indicate the number of samples in the respective sets.
Dataset
Features
Training
Validation
Test
Total
California Housing
16,512
2,064
2,064
20,640
Air Temperature
2,460
308
308
3,076
Road
391,385
4,348
39,140
434,873
Australian Census
5,978
48,471
6,059
6,059
60,589
Table 5: Performance metrics from four real-world datasets. In each line, the reported values refer to the implementation with the best performing kind of layer from GCN, GAT and GSAGE. Values are reported mean + to $0.975,9 SD/VI0. SMACNP metrics are not reported on the 3D Road and the Australian Census datascts due to high computational costs.
Malaset/Moder
Param.
Time
MSE
MAE
SQR
Cover. (%)
Calibr.
California Housing
GNN
2,529
48.33
5.07
1.62 6.00
9.24 0.000
3.40 + 0.08
93.88
0031
37.35 + 5.84
PE-GNN
27,426
29.37
4.7.0
1.07 1.008
7.21
1.0.7
2.73
1. 0.088
93.94
L.44
39.83
9.28
SMACNP
748,482
73.46
2.35 +-165
8.57
+.0.8
4.64
4-0.0
+-1.51
PE-GQNN (ours)
27,417
14.99
0.85 0.00
5.93 +
2.26
94.46
0:50
14.43
14.81
Air Temperature
GNN
1,281
1.06
= 0.06
2.62 + 0.20
12:43 = 0.45
4.19 0.20
92.63 + 1.65
65.25
+ 16,09
PE-GNN
27.106
12.83
1:2211
0.46
1: 0.07
5.04 0.37
1.89
1: 0.16
94.35
1: 1.49
124.72
60.40
SMACNP
741,482
23.32
0.18
2.77
3.80
100.00
600.60
12.07
PE-GQNN (owss)
26,169
8.63
0:23
3,08
1.19 = 0.06
92.01
= 2,58
25.29
10.18
3D Road
GNN
1,153
7.20
1. 0.67
1.71 1.
10.22 1. 1.006
3.57 J. 0.01
94.36 1.021
36.11 11
4.57
27,106
$7.27
8.83
0.30
4.11
0.03
94.89
0.23
26.91
9.27
PE-GQNN (ours)
23,897
6.85
0.01 1.
0.59 1.00
0,23 1 0.01
93.53 1: 1.55
179.81
88.10
Australian Census
GNN
384,737
22.75 6.35
0.19 +.008
2,90 +-0.6
1.08 + 0.08
95.69 + 0.53
51.51
10.08
PE-GNN
409 634
18.89 4.63
0.24 + 0.00
3.43 0.48
1.26 + 0.18
95.73 + 0.78
52.18 18.14
PE-GQNN (ours)
109,625
33.26
0.16 0.00
2.41 0.52
0.90
104.61 61.56
Figure g illustrates the spatial distribution of test average absolute residuals, aggregated within hexagonal bins of selected methods. The predictions obtained with PE-GQNN appear lighter, especially for the 3D Road dataset. As expected, the highest residuals are concentrated low-density areas.
5 Discussion
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 experiments real-world datasets, we have demonstrated the significant advantages of the PE-GQNN over competitive methods. The empirical results underscored the capability of the PE-GQNN to achieve lower MSE, MAE, and SQR across all datasets and GNN backbones compared to traditional GNN and PE-GNN, Notably, the PE-GQNN demonstrated substantial improvements in predictive accuracy and uncertainty quantification. The PE-GQNN framework's ability to provide a full description the predictive conditional distribution, including quantile predictions and prediction intervals, provides notable improvement geospatial machine learning. The PE-GQNN provides a solid foundation for future advancements in the field of geospatial machine learning.
Limitations: Although the Lipschitz Monotonic Networks (LMN) (Noltc et al. 2023 can approximate all monotonic Lipschitz bounded functions, many common quantile functions do not satisfy this condition. This often occurs because the derivative 0g(r)/0r tends to oo as T approaches 1, and to -00 as т approaches 0. The practical implications of this limitation when using LMN blocks remain to be fully understood. However, based on our experiments, it does not appear to pose major concern. LMNs are unlikely to struggle when approximating truncated distributions, where the support is slightly restricted and the derivative 0q(r)/Or
11

---

Published in Transactions on Machine Learning Research (12/2025)
Residual
(a) PE-GNN (California Housing).
(b) PE-GQNN (California Housing).
Residual
0.5
0.4
0.3
0.2
0.0
(c) PE-GNN (Australian Census).
(d) PE-GQNN (Australian Census).
Residual 0.20
0.15
0.10
0.00
(e) PE-GNN (3D Road).
(I) PE-GQNN (3D Road).
Figure 3: Geographical maps depicting the mean absolute test residuals across ten runs, aggregated within each hexagonal bin.
remains bounded. Naturally, any advances in monotonic neural network architectures would directly enhance the effectiveness of our approach.
Acknowledgments
Scott Sisson's 's research was supported by the Australian Research Council. David Nott's research was supported by the Ministry of Education, Singapore, under the Academic Research Fund Tier 2 (MOE- T2EP20123-0009). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior Brasil (CAPES) Finance Code 001.
References
Luc Anselin. Local indicators of spatial association LISA. Geogruphicul analysis, 27(2):93 115, 1995.
Luc Anselin. Spatial econometrics. Handbook of Spatial Analysis in the Social Sciences, pp. 101-122 2022.
Gabriel Appleby, Linfeng Liu, and Li-Ping Liu. Kriging convolutional networks. In Proceedings of the AAAI Conference Artificial Intelligence, volume 34(04), pp. 3187-3194, 2020.
12

---

Published in Transactions on Machine Learning Research (12/2025)
Li-Li Bao, Jiangshe Zhang, and Chunxia Zhang. Spatial multi-attention conditional neural processes. Neural networks: : the official journal of the International Neural Network Sociely, 173: 173:106201, 2024. URT. https:
sapiisemmattcschooaar org/CorpusID 268189461
Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Accurate medium-range global weather forecasting with 3d neural networks. Nature, 619(7970):533--5) 2023.
N. Cressie and C.K. Wikle. Statistics for Temporall Data. CourseSmart Series. Wiley, 2011.
Noel Cressie, Matthew Sainsbury-Dale, and Andrew Zammit-Mangion. Basis-function models in spatial statistics. Annual Review of Statistics and Its Application, 9:373-400, 2022.
Abhirup Datta, Sudipto Banerjee, Andrew O Finley, and Alan E Gelfand. Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514):800- 812, 2016.
Austin Derrow-Pinion, Jennifer She, David Wong, Oliver Lange, Todd Hester, Luis Perez, Marc Nunkesser, Seongjae Lee, Xueying Guo, Brett Wiltshire, et al. Eta prediction with graph neural networks in google maps. In Proceedings of the 30th ACM international conference on information & knowledge management, pp. 3767-3776, 2021.
Matthias Fey and Jan Eric Lenssen. Fast graph representation learning with pytorch geometric, 2019.
Reinhard Furrer, Mare G. Genton, and Douglas Nychka. Covariance tapering for interpolation of large spatial datasets. Journal of Computational and Graphical Statistics, 15(3):502-523, 2006.
Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning.
MIT Press, 2016. http://wwx.
deeplearningbook.org org.
Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. T. Guyon, Von Luxburg, Bengio, H. Wallach, Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Newral Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Josh Hooker, Gregory Duveiller, and Alessandro Cescatti. global dataset of air temperature derived from satellite remote sensing and weather stations. Scientific Data, 5(1):1-11, 2018.
Anirudh Ameya Kashyap, Shravan Raviraj, Ananya Devarakonda, Shamanth R Nayak K, Santhosh KV, and Soumya J Bhat. Traffic flow prediction models-a review of deep learning techniques. Cogent Engineering, (1):201055 2022.
Matthias Katzfuss and Joseph Guinness. A general framework for Vecchia approximations of Gaussian processes. Statistical Science, 36(1):124 141, 2021.
Manohar Kaul, Bin Yang, and Christian Jensen. Building accurate 3d spatial networks to enable next generation intelligent transportation systems. In 2013 IEEE 14th International Conference on Mobile Dula Management, volune 1, pp. 1471146.77 IEEE, 2013.
Diederik Kingma and Jimmy Ba. ADAM: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), San Diega, CA, USA, 2015.
Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations (ICLR), 2017.
Konstantin Klemmer and Daniel B Neill. Auxiliary-task learning for geographic data with autoregressive embeddings. In Proceedings of the 29th International Conference on Advances in Geographic Information Systems, pp. 141 144, 2021.
Konstantin Klemmer, Nathan Safir, and Danicl B Neill Positional encoder graph neural networks for geographic data. In International Conference on Artificial Intelligence and Statistics, pp. 1379-1389. PMI.R, 2023.
13

---

Published in Transactions on Machine Learning Research (12/2025)
Lucas Kock, G. S. Rodrigues, Scott A. Sisson, Nadja Klein, and David J. Nott. Calibrated multivariate regression with localized pit mappings, 2024. URL httpp:://arxiv.orr/ab//// /00bs///249.. 10855
Roger Koenker and Gilbert Bassett Jr. Regression quantiles. Econometrica: journal of the Econometric
Society, pp. 33-50 33-50, 1978.
Volodymyr Kuleshov and Shachi Deshpande. Calibrated and sharp uncertainties in deep learning via density estimation. In ICML, pp. 11683-11693, 2022.
Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon. Accurate uncertainties for deep learning using calibrated regression. In International conference on machine learning, pp. 2796-2804. PMLR, 2018.
Finn Lindgren, Håvard Rue, and Johan Lindström. An explicit link between Gaussian fields and Gaus- sian Markov random fields: The stochastic partial differential equation approach. Journal of the Royal Statistical Sociely Series B: Statislical Melhodology, 73(4):423-44) 498, 2011.
Haitao Liu, Jianfei Cai, Yew-Soon Ong, and Yi Wang, Understanding and comparing scalable gaus- sian process regression for big data. Knowledge-Baase Based Systems, 164:324-338 2019. ISSN 0950-7051. doi: httpss///do..org/10000 http://ddooooor//.....ojoooooo.. atticleepp1//09660061188833
URL http:///www.cieeeddreeee.. cwt.eeeddreet..oooocccteeeeeeeee
Yisheng Lv, Yanjie Duan, Wenwen Kang, Zhengxi Li, and Fei- Wang. Traffic flow prediction with big data: A deep learning approach. IEEE Transactions on Intelligent Transportation Systems, 16(2):865-873, 2014.
Will Mackey.
absmapsdata: A cotalogue of ready-to-ur use wattaackee/absmassaaae.o
ASGS (and other) sf objects, R package version 1.4, commit
2025. URL https://github.. com/ www.///wwww/////e////eeeeetwoeee 534-1119990000000077770700000200
Gengchen Mai, Kraysztof Janowiez, Bo Yan, Rui Zhu, Ling Cai, and Ni Lao. Multi-scale representation learning for spatial feature distributions using grid cells. In International Conference on Learning Repre- sentalions, 2020.
Niklas Nolte, Ouail Kitouni, and Mike Williams. Expressive monotonic neural networks. In The Eleventh International Conference on Learning Representations, 2023. URL sttp:::/oppmmeriisr.etfffoooo
w2P7fMy RH
R Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3): 291 297, 1997.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 32, 2019.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006.
Filipe Rodrignes and Francisco C Percira. Beyond expectation: Deep joint mean and quantile regression for spatiotemporal problems. IEEE Transactions on Newral Networks and Learning Systems, 31((2): 5377-5389, 2020.
Phillip Si, Volodymyr Kuleshov, and Allan Bishop. Autoregressive quantile flows for predictive uncertainty ostimation. In International Conference on Learning Representations, 2022.
BR Sreenivasa and CR Nirmala. Hybrid location-centric e-commerce recommendation model using dynamic behavioral traits of customer. Iran Journal of Computer Science, 273):179-188 2019.
Natasa Tagasovska and David Lopez-Paz. Singlc-model uncertaintics for deep lcarning. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, Red Hook, NY, USA, 2019. Curran Associates Inc.
14

---

Published in Transactions on Machine Learning Research (12/2025)
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in Neural Information Processing Systems, 30, 2017.
V. Vecchia. Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society: Series B (Methodological), 50(2))297-312, 1998.
Petar Veličkovié, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. Internalional Conference on Learning Represenlalions, 2018.
Lingfei Wu, Peng Cui, Jian Pei, Liang Zhao, and Xiaojie Guo. Graph neural networks: Foundation, frontiers and applications. Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 4840-4841, 2022.
Shuai Xu, Xiaoming Fu, Jiuxin Cao, Bo Liu, and Zhixiao Wang. Survey user location prediction based geo-social networking data. World Wide Web, 23(3):1621-1664 2020.
A
Appendix
A.1
Quantile Neural Networks and recalibration
Figure provides visual overview of some of the most closely related approaches.
UINEAR
dtt
pr) 4:(r"))
X
p
9
X
p. (z()))
X
:
(ye (x1)))
X
que
PraOi (((r))
d+1
MSE(P.))
n
(a) Non-linear quantile regression using NN.
(b) Non-linear multiple quantile regression.
x
T
X=
r!
X
C
MSE(Y,)) y)
U(0.1)
7
U(0.1)
Φ
m
4r
X
d
Pr p((T/)) n.d d
p, a((x/)) p
12 d
q,4
(c) Non-linear quantile function regression.
(d) Two-step density estimation.
Figure 4: (a) For each quantile of interest, a separate NN is trained. (b) Rodrigues Pereira 2020): one NN
outputs + predictions: one for the expectation and d for the quantiles. (i) et al. (2022): single NN trained predict any generic quantile of the conditional distribution. (d) Kuleshov & Deshpande (2022):
two-step procedure: the first model outputs low-dimensional representation of the conditional distribution, which a recalibrator then uses to produce calibrated predictions.
15

---

Published in Transactions on Machine Learning Research (12/2025)
A.2 Positional Encoder
Inspired by the Translormer architecture (Vaswani et al. 2017) Mai et al. (2020) introduced the Positional Encoder (PE) for geographic data. The PE maps the geographic coordinate vector с = (c1, 02) E RR? of a single datapoint, typically representing latitude and longitude, to a high-dimensional embedding using a set of deterministic sinusoidal transformations followed by fully connected neural network.
The encoded spatial embedding c E Rd is computed as:
c emmb PF((c, min, max: PED NN(ST minn max)) OpE),
(4)
where (., OPF) denotes a fully connected neural network with parameters OPF, and ST(-) is sinusoidal transformation defined by:
(T, Omin Tmax) (ST (c, minn Tmax);.... sT, minn Tmax))
(5)
with min and max hyperparameters. Each component ST, (c, min: max), for E (1,...,,S), .S) is given by:
2T jj
2Tez jj
2TC2
(c, mmins max)) = sin
.00
sin
.oo8
(6)
jj
where the values d og form a logarithmically spaced grid between amin and maxx computed as:
0 max
jj Mmin
The result of ST(c,-) is a 4S-dimensional vector that encodes spatial patterns at multiple scales. This vector is then passed through the neural network NN to produce the final embedding cemb.
A.3 Architectural details
To facilitate clear and detailed comparison between PE- GENN and PE-GNN, present the architectures for the California Housing dataset in Tanles The same architectures are used for the Air Temperature dataset, differing only in the mumber of node fcatures (p 1). For the 3D Road dataset, which is an interpolation task without node features, the GNN block is omitted.
16

---

Published in Transactions on Machine Learning Research (12/2025)
Table 6: Architecture of PE-GQSAGE applied to the California Housing dataset, illustrating input/ dimensions, parameter counts, layer-level annotations, and activation functions.
Layer (type & shape)
# Param.
Notes
Activation
GNN block: 6 + 32
Input: node features (dim
Node features input
SAGEConv (hidden): 6 32
GraphCony layer
Aggregation
224
GSAGE internal layer
Update
192
GSAGE internal layer
ReI.U
SAGEConv (out)::
GraphCony layer
Aggregation
1,056
GSAGE internal layer
Update
1,024
GSAGE internal layer
ReLUU
Positional encoder (PE) block: 64
Input: coordinates (dim 2)
Spatial coordinates input Fourier feature mapping
Sinusoidal Trausformation (ST): 2 - 64
Dropout regularization
Dropout rate: p 0.5
Linear (hidden): +
8,320
Feedforward layer
ReLU
Linear (hidden): 128 + 64
8,256
Feedforward layer
Tanh
Linear (hidden): 64 - 32
2,080
Feedforward layer
Tanh
Linear (out): 32 - 64
2,112
Feedforward layer
Identity
Concatenation: GNN (32) + (64)
Combined vectors
Fully-comnected (FC) block: 96 +
Dimension reduction block
Linear (hidden): 96 + 32
3,104
Feedforward layer
Tanh
Linear (hidden): ->
528
Feedforward layer
Tanh
Linear (out): 16 4 8
136
Feedforward layer
Identity
Concatenation: FC (8) ÿ (1) (1)
Combined vectors
Monotonic Layers: +
Final quantile regressor
Lipschitz Linear (hidden): 10 + 32
352
Monotonic layer
GroupSort(2)
Lipschitz Linear (out): 32
33
Monotonic layer 2
Identity
Total
27,417
trainable parameters
Table 7: Architccture of PE-GSAGE. Positional encodings derived from spatial coordinates are concate- nated with node features prior to the GNN layers. This architecture follows the specification introduced by Klemmer et al. )023).
Layer (type & shape)
# Param.
Notes
Activation
Positional Encoder (PE) block: 2
Input: coordinates (dim 2)
Spatial coordinates input
Sinusoidal Transformation (ST): 2 + 61
0
Fourier feature mapping
Dropout regularization
Dropout rate: p 0.5
Linear (hidden): 64 + 128
8,320
Feedforward layer
ReLU
Linear (hidden): 128 + 64
8,256
Feedforward layer
Tanh
Lincar (hidden): - 32
2,080
Foodforward layer
Tanh
Linear (out): 32 + 61
2,112
Feedforward layer
Identity
Concatenation: PE (64) Features (6)
Combined input vector
GNN block: 70 - 32
Input: node features (dim 70)
PE Node attributes
SAGEConv (hiddeu): 70 - 32
Graph layer
Aggregation Linear
2,272
GSAGE internal layer
Update Linear
2,240
GSAGE internal layer
ReLU
SAGEConv (out): 32
Graph layer
Aggregation Linear
1,056
GSAGE internal layer
Update Linear
1,024
GSAGE internal layer
ReLU
Fully-connected layers: 32 -
Linear (to target): 32
33
Final output layer
Idetiity Identity
Linear (to Moran's I): 32 -
Auxiliary output layer
Total
27,426
All trainable parameters
17


Reviews:
[{"review_id": "xr4JZ5ZAb5", "reviewer": "Reviewer_kFfa", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces the Positional Encoder Graph Quantile Neural Network (PE-GQNN), which integrates positional-embedding-aware graph neural networks with quantile regression to enable uncertainty-aware spatial prediction. Building on prior work, the approach removes the need for post-hoc calibration by combining prediction and uncertainty estimation within a single end-to-end framework. The method appears robust in high-dimensional predictor spaces and produces well-calibrated predictive distributions. Empirical results across several real-world datasets show that PE-GQNN achieves strong and often superior performance compared to established baselines across multiple evaluation metrics.", "claims_and_evidence_yn": "No", "claims_explanation": "The paper is well-written and presents an interesting contribution, but there are several concerns regarding the evidence provided and some methodological choices.\n\nFirst, while the paper emphasizes the scalability of the proposed components, there is no theoretical analysis of the model’s computational complexity. Table 2 provides runtime estimates for a single dataset, but without explicit details on dataset sizes, it is difficult to assess the actual computational cost or scalability limits of the approach.\n\nSecond, the claims regarding uncertainty quantification are largely theoretical. The experimental evaluation provides limited evidence that the model delivers reliable or calibrated uncertainty estimates in practice. \n\nFrom a methodological standpoint, there are no theoretical or empirical guarantees of model stability or convergence. The graph is recomputed for every batch, which means the adjacency structure changes dynamically. This introduces additional stochasticity that may not always improve generalization and could affect training stability. Moreover, the impact of batch size on both performance and computational cost is not analyzed, although it is likely to play an important role.\n\nThe decision to use a single Monte Carlo sample for τ (quantile level) is also questionable. While the law of large numbers offers theoretical justification, it would be important to examine in practice how this choice affects the variance of gradient estimates and the stability of training.\n\nFinally, although Lipschitz Monotonic Networks are introduced to prevent quantile crossing, these constraints may restrict model flexibility and still yield only approximately monotonic quantiles. A deeper empirical investigation of this trade-off would strengthen the work.\n\nOverall, the robustness analysis is insufficient to fully validate the proposed method. Performance results should be reported as averages over multiple runs with associated error estimates to provide a more reliable assessment of model consistency and generalization.", "audience_yn": "Yes", "audience_explanation": "Despite some limitations in the experimental setup, this paper represents a strong contribution to the development of GNN-based methods for geographical data. The integration of quantile neural networks with positional encoders is novel and adds clear value to the graph-ML community. By producing full predictive distributions rather than single point estimates, the proposed model enables uncertainty-aware spatial prediction, an aspect still uncommon in graph-based frameworks. Although primarily motivated by spatial applications, the framework’s quantile-based uncertainty modeling and monotonic design principles could, in principle, extend to non-spatial or relational graph problems, underscoring its potential broader impact.", "requested_changes": "Please see above in the claims and evidence section.", "broader_impact_concerns": "Not applicable", "additional_comments": "Some references should be added for specific arguments in the paper: for example, the introduction states “...poorly calibrated intervals, such as 80% intervals that fail to contain the true outcome 80% of the time,” but no citation is provided.", "extra_scores": {}, "extra_text": {}}, {"review_id": "hQY4wG8xM5", "reviewer": "Reviewer_aaBq", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a few modifications to the PE-GNN method to enhance its ability to make accurate spatial predictions and to quantify uncertainty. The following modifications are introduced:\n(1) a new method fro quantile regression task that merges the two-step procedure of Kuleshov & Deshpande (2022) into a single model by postponing the concatenation of the τ value. τ is transformed beforehand to facilitate learning and provide better\ncontrol over the form of the predicted distribution. This approach aims to \nenhances the network’s ability to model uncertainty and to make it more robust to high-dimensional\npredictor spaces. \n(2) the proposed approach modifies the way the nodes’ features and the spatial embeddings are concatenated in PE-GNN.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The proposed approach teases on 3 geospatial datasets and the performance is compared to current GNN approaches (GNNs, PE-GNN, SMACNP and GSAGE). The approach demonstrates some improvement on a variety of metrics, including running time, MSE, MAE etc.", "audience_yn": "Yes", "audience_explanation": "Making geospatial predictions more robust is of interest to geospatial and, potentially, machine learning audience. The performance  improvement is not significant, though.", "requested_changes": "As with many geospatial datasets, the visualisation of the results might be useful. Figure 4 in appendix might be brought up in the main paper (or even better, it can be dome for one of the other datasets). However, the captions should clearly describe what does the colour represent. Additionally, it'll be good to demonstrate not the predictions themselves, but the difference between real and predicted values: it will make the improvement in the model performance more obvious.", "broader_impact_concerns": "I do not see major ethical implications of this work", "additional_comments": "In the introduction, the contributions are quite wordy and not easy to follow", "extra_scores": {}, "extra_text": {}}, {"review_id": "5KsJp8oEtu", "reviewer": "Reviewer_q7n9", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors introduced the Positional Encoder Graph Quantile Neural Network (PE-GQNN), a novel framework designed to improve predictive modeling for geographic data. The main contributions of the study are as follows: Experimental results demonstrated that PE-GQNN consistently outperformed traditional GNN and PE-GNN models, achieving lower values in MSE and MAE metrics. Furthermore, the model showed significant gains in both predictive accuracy and the quality of uncertainty quantification.", "claims_and_evidence_yn": "Yes", "claims_explanation": "**Strengths** :\n\n- the proposed method combines both quantile prediction and distribution calibration in a single model for more efficiency\n- the proposed method allows for good predictions accross any quantile level\n- the experiments accross multiple dataset are convincing\n\n**Weaknesses** :\n\n- The technical contribution of the paper is limited. Indeed, the proposed method is essentially a combination of PE-GNN and quantile regression. Techniques such as response variables and quantile parameters into neural networks are included. However, these contributions are incremental rather than fundamentally novel.\n\n- The paper does not address the computational cost associated with applying the method to high-dimensional data, which could be a significant limitation in real-world scenarios.", "audience_yn": "Yes", "audience_explanation": ".", "requested_changes": "See weaknesses section.", "broader_impact_concerns": "None", "additional_comments": "None", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)