Graph neural networks (GNNs) are routinely employed for spatiotemporal forecasting, yet their performance across widely used benchmark datasets is inconsistent. Here, we perform an audit of dataset properties and baseline models to assess the quality of the benchmarks, and the robustness of the conclusions drawn from them. Using classical statistical tools, we characterise spatiotemporal lagged dependencies in benchmarks, and examine how temporal differencing changes these relationships and affects model rankings. Motivated by this, we re-evaluate temporal linear baselines, significantly reducing the apparent gains from GNNs on several benchmarks, and surpassing GNNs on others. Suspecting that GNNs struggle to extract linear, node-wise signals, we find that supplying them with autoregressive residuals improves their performance particularly on non-traffic benchmarks. Finally, controlled synthetic experiments reveal that GNNs are sensitive to heterogeneity in temporal dynamics and spatial graph interactions. Together, our findings demonstrate that baseline specification, data pre-processing and system heterogeneity shape the interpretations drawn from benchmark rankings, informing the design and robust evaluation of GNNs.
Figures & tables
Figure 1 : Temporal (blue) and partial spatial (green) correlations for levels (solid) and first differences (dashed). Lines show nodewise medians; shaded regions span the 2.5th–97.5th percentiles across nodes, with hatching for first differences. These regions describe between-node variation, not confidence intervals. Lags extend to 3H=12 for H=4 . Partial spatial correlation controls for the same-lag own-node observation and is undefined at lag zero.
Model
Chickenpox
PedalMe
WikiMaths
Neural and graph-based baselines
GConvGRU [ Seo et al., 2018 ]
0.779 ± 0.009
1.372 ± 0.005
0.443 ± 0.002
Graph WaveNet [ Wu et al., 2019 ]
0.799 ± 0.017
1.424 ± 0.121
0.431 ± 0.005
TDE-GNN [ Eliasof et al., 2024 ]
0.782 ± 0.006
0.714 ± 0.051
0.565 ± 0.017
LSTM [ Hochreiter and Schmidhuber, 1997 ]
0.825 ± 0.005
1.405 ± 0.026
0.504 ± 0.002
Statistical baselines and residual models
Table 1 : Test MSE ( ↓ ) under a chronological 9:1 train–test split. All results were run and independently confirmed by us. Neural and residual models are reported as mean ± standard deviation over 10 seeds; ARIMA models are deterministic. For WikiMaths, ARIMA-family models use a period-7 SARIMA specification. Bold , underline , and italic denote the first-, second-, and third-lowest MSE, respectively, within each dataset column; ties are retained.
Model
METR-LA
PEMS-BAY
PEMS03
PEMS04
PEMS07
PEMS08
Literature: BasicTS+
DCRNN [ Li et al., 2018 ]
3.03
1.59
15.54
19.66
21.16
15.23
Graph WaveNet [ Wu et al., 2019 ]
3.03
1.59
14.59
18.80
20.44
14.67
MTGNN [ Wu et al., 2020 ]
3.05
1.60
14.85
19.13
21.01
15.25
DGCRN [ Li et al., 2023 ]
2.94
1.58
14.60
18.84
20.04
14.77
D 2 STGNN † [ Shao et al., 2022b ]
2.88
1.52
14.63
18.32
19.49
14.10
Table 2 : Traffic forecasting on six datasets of BasicTS+ benchmark [ Shao et al., 2024 ] . Entries are test MAE averaged over all 12 forecast steps (12 input steps; lower is better). METR-LA and PEMS-BAY use 7:1:2 train/val/test splits; PEMS03/04/07/08 use 6:2:2. † D 2 STGNN and STID use learned day-of-week embeddings; our neural models receive time of day but no explicit day-of-week feature. Bold identifies the lowest mean among our displayed models for each dataset.
Figure 2 : Recovered share of spatiotemporal GNNs versus ARIMA in closed form ceiling: 100⋅(MSEARIMA−MSEmodel)/(MSEownpast−Noisefloor) % across three synthetic datasets varying node heterogeneity, spatial dispersion and spatiotemporal memory.
Appendix figures & tables8 assets
Supplementary material from the paper’s appendix.
Appendix
Dataset
Target
Nodes
T
Freq.
H→K
Split
Chickenpox
Cases
20
522
Weekly
4→1
9:1∗
PedalMe
Deliveries
15
36
Weekly
4→1
9:1∗
WikiMaths
Page views
1,068
731
Daily
8→1
9:1∗
METR-LA
Speed
207
34,272
5 min
12→12
7:1:2
PEMS-BAY
Speed
325
52,116
5 min
12→12
7:1:2
PEMS03
Flow
358
26,208
5 min
12→12
6:2:2
Appendix
Table 3: Dataset statistics and forecasting protocols. T counts observations before constructing forecasting windows. H→K denotes input and forecast lengths. Three-part splits denote chronological training/validation/test proportions. ∗ Nontraffic experiments reserve the final 10% for testing.
Figure 3 : Association between dataset-level correlation structure and the relative MSE improvement of the mean of Graph WaveNet and direct GRU-GCN over the preselected statistical baseline. Strength is the mean absolute nodewise correlation over the input window; relative dispersion is the RMS across-node standard deviation normalized by RMS correlation magnitude. Partial spatial correlation conditions on the node’s own value at the same lag. Labels identify dataset variants; colours distinguish traffic and non-traffic benchmarks.
Axis
Factor
Levels
Fixed at
A
node heterogeneity va=Var(ai)
0.005,…,0.030 (10, linear)
c=0 , p=1
B
weight variance c ( Var(wij)=c2/9 )
0,…,1.2 (10, linear)
va=0.005 , p=1
C
spatial lag p
1,…,10
va=0.005 , c=0
Appendix
Table 5: The three axes of the synthetic benchmark.
Model
Hyper-parameter
Values
Base Values
Learning rate
{3⋅10−4,10−3,3⋅10−3,10−2,3⋅10−2}
Weight decay
{0,10−5,10−4,10−3}
Hidden size
{16,32,64}
Batch size
{32,64,128}
Dropout
{0,0.1,0.3}
Appendix
Table 6: Hyper-parameter search spaces. T-GCN, A3T-GCN, EGCN-O, EGCN-H and MPNN-LSTM use the shared space only. For every model and dataset, 24 configurations are sampled uniformly (the 24 ARIMA orders are enumerated) and the one with the lowest validation MSE is refit with three training seeds.
Model
Type
Var(ai)
Var(wij)
spatial lag p
min
mean
max
min
mean
max
min
mean
max
A3T-GCN
F
30.9
62.3
92.3
51.0
73.6
92.3
−202.1
−122.7
92.3
AGCRN
S
−12.8
−6.3
3.8
−11.3
−1.8
3.8
−29.2
−20.8
3.8
DCRNN
S
55.8
75.3
98.0
33.3
66.1
104.0
13.4
82.4
122.0
DyGrAE
S
65.2
84.5
108.0
66.6
86.9
108.0
22.5
72.7
118.8
DynGESN
S
11.6
24.9
37.2
19.2
28.9
36.3
25.6
56.8
118.9
Appendix
Table 7 : Recovered share of the closed-form graph ceiling (%) along the three axes of the synthetic benchmark 100⋅(MSEARIMA−MSEmodel)/(MSEownpast−σ2) . Type: F = fused into one normalized A+I operator, S = separate pathway. Values obtained from three seeds (0,1,2) and input-output horizon of 12 → 1. Bold: best per column.
Figure 4 : Nodewise temporal, spatial and partial spatial correlations. Lines show medians; shading spans the central 95% across nodes, not confidence intervals. Lags extend to 3H ; the dashed vertical line marks H . A common set of nodes with finite estimates for all three statistics at every plotted lag is used.
Figure 5 : Lagwise relative dispersion of temporal and partial spatial node-wise correlations, H=SD/RMS . Note that high relative dispersion can occur even when correlations are weak. Lags extend to 3H ; the dashed vertical line marks H .
Dataset
Temporal HT
Partial spatial HG
METR-LA
0.223
0.455
PEMS-BAY
0.089
0.558
PEMS03
0.037
0.531
PEMS04
0.051
0.516
PEMS07
0.038
0.705
PEMS08
0.050
0.642
Appendix
Table 8 : Temporal and partial spatial heterogeneity within each dataset’s input window.
Centre for Artificial Intelligence Research (CAIR) University of Cape Town, Cape Town, South Africa · Nonlinear Dynamics and Chaos Group, Department of Mathematics and Applied Mathematics University of Cape Town, Cape Town, South Africa
Karlsruhe Institute of Technology, Germany · The Hong Kong University of Science and Technology, Hong Kong SAR · East China Normal University, China +1