# Spatiotemporal Graph Neural Networks with Uncertainty Quantification for
Traffic Incident Risk Prediction

Xiaowei Gao<EMAIL_ADDRESS>Xinke Jiang<EMAIL_ADDRESS>Dingyi Zhuang<EMAIL_ADDRESS>Huanfa Chen<EMAIL_ADDRESS>Shenhao Wang
<EMAIL_ADDRESS>James Haworth<EMAIL_ADDRESS>
###### Abstract

Predicting traffic incident risks at granular spatiotemporal levels is
challenging. The datasets predominantly feature zero values, indicating no
incidents, with sporadic high-risk values for severe incidents. Notably, a
majority of current models, especially deep learning methods, focus solely on
estimating risk values, overlooking the uncertainties arising from the
inherently unpredictable nature of incidents. To tackle this challenge, we
introduce the Spatiotemporal Zero-Inflated Tweedie Graph Neural Networks
(STZITD-GNNs). Our model merges the reliability of traditional statistical
models with the flexibility of graph neural networks, aiming to precisely
quantify uncertainties associated with road-level traffic incident risks. This
model strategically employs a compound model from the Tweedie family, as a
Poisson distribution to model risk frequency and a Gamma distribution to
account for incident severity. Furthermore, a zero-inflated component helps to
identify the non-incident risk scenarios. As a result, the STZITD-GNNs
effectively capture the dataset’s skewed distribution, placing emphasis on
infrequent but impactful severe incidents. Empirical tests using real-world
traffic data from London, UK, demonstrate that our model excels beyond current
benchmarks. The forte of STZITD-GNN resides not only in its accuracy but also
in its adeptness at curtailing uncertainties, delivering robust predictions
over short (7 days) and extended (14 days) timeframes.

###### keywords:

Spatiotemporal Sparse Data Mining , Urban Traffic Risk Prediction ,
Uncertainty Quantification , Zero-Inflated Tweedie Distribution

[inst1]organization=SpaceTimeLab, University College London (UCL),
state=London, country=UK

[inst2]organization=School of Computer Science, Peking University (PKU),
state=Beijing, country=China

[inst3]organization=Department of Civil and Environmental Engineering,
Massachusetts Institute of Technology (MIT), state=Cambridge, MA, country=USA
[inst4]organization=The Bartlett Centre for Advanced Spatial Analysis,
University College London (UCL) , state=London, country=UK
[inst5]organization=Department of Urban and Regional Planning, University of
Florida, state=Gainesville, FL, country=USA

##  1 Introduction

Urban areas, defined by their complex transportation networks and dynamic
traffic patterns, are further accentuated by diverse socio-economic activities
(Dai et al.,, 2018; Tu et al.,, 2020). While these features present a plethora
of opportunities, they also pose significant challenges for policymakers,
transportation authorities, and urban planners in managing human mobility-
resulted traffic incident risks (Barahimi et al.,, 2021; Rodrigue,, 2020).
Such risks negatively impact urban sustainability, leading to human
casualties, economic losses, and diminished livability. Addressing these
concerns necessitates the development of predictive models that can
effectively identify high-risk zones and reliably anticipate potential
incidents (Laureshyn et al.,, 2017; Zhou et al.,, 2021).

Traffic incident prediction is an interdisciplinary domain, intersecting urban
planning, transportation engineering, data science, and social sciences (Zhou
et al.,, 2020; Yu et al., 2020b, ; Bibri,, 2018). The overarching aim is to
discern patterns from historical spatiotemporal incident risk data to profile
or forecast traffic incidents (Chakraborty et al.,, 2019; Liu et al.,, 2022).
The current literature in this realm can be classified into two predominant
categories:

  * •

Incident injury severity identification and analysis: Researchers focus on
extracting patterns, causes, and contributing factors related to injury
severity from historical incident data. The objective is to understand the
interplay between various determinants such as road attributes, traffic
conditions, weather conditions, and human behaviour, and their impact on
individual injury severity levels (Hu et al.,, 2010; Benlagha and
Charfeddine,, 2020; Ma et al.,, 2021; Yang et al.,, 2021, 2022).

  * •

Incident risk prediction by spatiotemporal modelling: Given the significance
of spatiotemporal attributes, researchers have increasingly emphasized
modelling the dynamics and exposures in urban environments related to incident
risk. This involves capturing the spatial dependencies among different road
segments and considering the temporal variations in traffic patterns (Wang et
al., 2022a, ). By incorporating spatiotemporal modelling techniques,
researchers can effectively account for the evolving nature of urban traffic
and improve the accuracy of incident risk predictions.

Both methodologies offer advantages yet bear intrinsic constraints. As Liu et
al., (2023) underscores, injury severity causal analysis may not reliably
forecast future incidents. On the flip side, spatiotemporal modelling, while
presenting a comprehensive approach with the potential for enhanced accuracy,
faces challenges regarding imbalanced incident data and unreliable results.

Delving deeper into traffic incident risk prediction research, these
challenges can be elaborated upon as follows. The first is summarized as the
zero-inflation issue, a characteristic of risk data at fine-grained levels
(Trirat et al.,, 2023; Zhou et al.,, 2020). The issue describes the excessive
quantity of zeroes due to data sparsity. Conventional modelling approaches
often dilute this issue by aggregating risk values at more extensive regional
scales (Bao et al.,, 2019; Wang et al., 2021c, ). While this might yield
reasonable predictions from both broader and smaller regional perspectives, it
leads to unsatisfactory estimates at road levels plagued by overdispersion.
Hence, there is a pressing need for innovative modelling techniques capable of
accurately capturing and representing zero-inflation patterns in traffic
incident risk data.

A second pressing issue pertains to  uncertainty quantification (UQ). Broadly,
uncertainty can be divided into epistemic and aleatoric types, signifying
model and data uncertainty, respectively (Qian et al.,, 2022). Due to their
unpredictable and sporadic nature, traffic incidents serve as typical examples
of aleatoric uncertainty, which is irreducible compared to the epistemic one
(Bao et al.,, 2019; Zhou et al.,, 2021). While current deep learning models
excel at learning intricate patterns and accurately predicting the average
risk, they often implicitly assume homoskedasticity (uniform variance) and
then inadequately account for the inherent variability in the results of their
predictions (Wang et al., 2023a, ; Zhuang et al.,, 2022; Jiang et al.,, 2023).
Additionally, while some strides have been made in UQ for time-series data
prediction (Wu et al.,, 2023; Zhu et al.,, 2022; Li et al., 2022c, ), such
efforts are still scarce in sparse data, particularly in risk prediction.

To address the aforementioned challenges, our approach prioritizes UQ, with a
specific focus on the prediction interval of risk estimation. We introduce a
novel approach — the Spatiotemporal Zero-Inflated Tweedie Graph Neural Network
(STZITD-GNNs)— for the finest granularity at the road level among existing
spatiotemporal deep learning models that consider UQ. This method is designed
to boost the proficiency of spatiotemporal deep learning models in UQ by
interacting with the statistical methods, consequently improving the
robustness of road-level risk prediction.

Our concentrated focus on finer spatiotemporal granularity naturally results
in a more substantial graph structure compared to traditional grid-based
datasets. This poses considerable challenges in efficiently handling node
embeddings and producing outputs. Furthermore, road segments typically present
a low degree of observed daily risk, with a large majority displaying no
incidents occurrence throughout the year. This results in an extremely
unbalanced dataset, exhibiting significant zero-inflation for risk assessment
(Zhou et al.,, 2020; Wu et al.,, 2022). Meanwhile, the distribution of risk
severities’ frequency also follows a skewed and long-tailed pattern. The
majority of incidents are categorized as minor in severity, contributing to
slight risk (Shirazi and Lord,, 2019; Saha et al.,, 2020). Given these
circumstances, the common assumption of a Gaussian distribution for entities
in graph neural network prediction becomes unsuitable (Wang et al., 2023a, ).
While the state-of-the-art approach proposed by Zhuang et al., (2022) provides
encouraging results for UQ in sparse travel demand scenarios, their
methodology—grounded in the negative binomial (NB) distribution—limits its
focus to binary outputs, neglecting the real full distribution of incident
risk values in our problems.

To further overcome the above limitations, we propose the use of the Tweedie
(TD) distribution, a compound Poisson-gamma distribution, which simultaneously
models the exact zero values and the continuous positive values. The TD
model’s variance function is analogous to that in the NB model but
incorporates a flexible index parameter $\rho\in(1,2)$ (Lord and Mannering,,
2010; Shirazi and Lord,, 2019). The TD distribution has been extensively used
in various domains with overdispersion data, including crash frequency
estimation (Wang et al.,, 2019), insurance claims (Smyth and Jørgensen,,
2002), actuarial studies (Shi, (2016)), and meteorological precipitation
(Dunn,, 2004).

However, when confronted with the challenge of extremely unbalanced zero-
inflated data and skewed positive data in urban risk values, the TD
distribution alone may not provide satisfactory performance as it cannot
capture both aspects of the distribution (Khosravi et al.,, 2010; Abe and
Yadohisa,, 2017)). Therefore, we propose a mixture model strategy known as the
”Zero-Inflated Tweedie” (ZITD) model. Unlike the traditional TD model, the
ZITD model incorporates an additional parameter $\pi$ to estimate the
probability mass of zero values, in addition to the parameters of the TD model
component. This allows for a more comprehensive understanding of risk values,
accounting for both the zero mass and the skewed positive component of the
distribution.

Unlike conventional methodologies that estimate spatial and temporal outputs
separately before fusing them (Zhuang et al.,, 2022; Wang et al., 2023a, ),
our approach harnesses the power of joint spatiotemporal graph neural networks
to simultaneously estimate and fit the four parameters of the ZITD
probabilistic distributions. In detail, we employ a Gated Recurrent Unit (GRU)
as the temporal encoder and a Graph Attention Network (GAT) as the spatial
encoder. These components work in tandem to simultaneously capture the
intricate spatiotemporal dependencies in the data. Through extensive empirical
evaluation using real data from three boroughs of London, UK, we provide
compelling evidence of the superior performance of our model in road-level
urban risk prediction tasks. Our major contributions can be summarized as
follows:

  * •

We propose a novel framework, STZITD-GNNs, for estimating the spatiotemporal
uncertainty of traffic incident risk, specifically addressing the challenges
posed by extremely unbalanced zero-inflated and long-tail data at the road
level. By incorporating probabilistic GNNs, our framework successfully
quantifies the sparse and discrete uncertainty, particularly in high-
resolution risk datasets.

  * •

Through the integration of spatiotemporal graph neural networks (ST-GNNs), our
framework jointly captures the temporal dynamics and spatial relationships
within urban data, enabling a holistic insight into the inherent uncertainties
and intricate risk patterns and also ensuring a robust and nuanced approach to
urban risk assessment.

  * •

Extensive experiments are conducted on real-world datasets to evaluate the
performance of the STZITD-GNNs. The results demonstrate that our proposed
method surpasses state-of-the-art methods, validating the effectiveness and
superiority of the STZITD-GNNs framework in tackling sparse urban risk
prediction challenges.

To our knowledge, our research advances traffic incident risk prediction by
transitioning from deterministic deep learning to a probabilistic approach,
emphasizing road-level over grid-based models. This method adeptly tackles the
zero-inflation and long-tail distribution in granular risk data, enhancing
accuracy in traffic incident risk assessment. Our study underscores the
importance of UQ in road-level risk prediction, providing pivotal insights for
future research and applications.

This manuscript is structured as follows: In Section 2, we conduct a rigorous
review of spatiotemporal modelling methodologies in the domain of traffic risk
prediction and then delve into Uncertainty Quantification (UQ). Section 3
provides a detailed description of the proposed model. Section 4 compares our
empirical results against state-of-the-art benchmark models, further
augmenting our assertions with illustrative visualizations.In Section 5, we
delineate the salient attributes of our model and suggest avenues for future
enhancements.

##  2 Literature Review

###  2.1 Traffic Incident Risk Prediction under Spatiotemporal Modelling

Traffic incident risk prediction is a crucial aspect of urban planning and
mobility management. Early studies in this field have primarily relied on
traditional statistical methods or classical machine learning techniques to
analyze and forecast incidents on specific road segments or highways. These
approaches include regression models (Caliendo et al.,, 2007; Bergel-Hayat et
al.,, 2013), Bayesian Networks (Hossain and Muromachi,, 2012), tree-based
models (Chong et al.,, 2004; Lin et al.,, 2015), and k-nearest neighbour
method (Lv et al.,, 2009). However, these studies have initially focused on
small areas, limiting their applicability to broader risk management knowledge
and considering only a limited set of influential factors.

Recognizing the need to overcome these limitations, researchers have expanded
their investigations to encompass broader traffic zones, such as borough or
citywide levels, to provide a more comprehensive understanding of urban risk
(Chang and Chen,, 2005; Anderson,, 2009; Zhang et al.,, 2014). However, it is
crucial to acknowledge that risk data possess inherent complexities in their
spatial and temporal attributes, which pose challenges that these classical
approaches may not effectively address. These methods often handle the
individual incident data separately without considering the interdependencies
among different locations, which can result in inadequate accuracy levels in
risk prediction(Wang et al., 2021a, ). Moreover, when dealing with larger
areas characterized by finer spatial units, the limitations of these
traditional methods become more pronounced.

To overcome these challenges, recent advancements in deep learning techniques
have shown promise in urban traffic risk prediction by leveraging the power of
spatiotemporal modelling. Convolutional Neural Networks (CNNs), Recurrent
Neural Networks (RNNs), and Graph-based networks, like Graph Convolution
Networks (GCNs), have demonstrated their effectiveness in capturing intricate
spatial patterns, capturing temporal dependencies, and modelling complex
interactions among urban areas (Wang and Cao,, 2021; Liu et al.,, 2020). As
summarized in Table 1, Chen et al., (2016) were among the first to explore
national-wide risk prediction by employing a Stack Denoise Autoencoder. They
integrated human mobility GPS data and historical incident point data at the
grid level to map the real-time risk situation in Tokyo. However, while their
approach attempted to capture a larger spatial area, it did not consider
geographical content for precise and long-term risk prediction. Later, Chen et
al., (2018) introduced a Stack Denoise Convolutional Autoencoder that
incorporated spatial dependencies in incidents using stacked CNNs.
Nonetheless, both of these studies overlooked the influence of temporal
factors. Although RNNs are capable of handling temporal information, they are
more suited for short-term temporal learning (Ren et al.,, 2018). To better
address the temporal aspect, Ren et al., (2018) utilized multiple Long Short-
Term Memory (LSTM) to consider temporal influential factors among several
spatial locations and predict risk at the city level. Their work represented
the first attempt to incorporate temporal factors in deep learning models for
city-wide traffic risk prediction. Similarly, Moosavi et al., (2019) employed
LSTM networks for national-wide risk prediction, demonstrating their efficacy
in predicting risk at 15-minute intervals. Building upon the LSTM framework,
Yuan et al., (2018) proposed a more advanced model known as the Convolutional
Long Short-Term Memory (Hetero-ConvLSTM) neural network based on heterogeneous
data. This model aimed to jointly address spatial heterogeneity and capture
temporal features in state-wide incident prediction. However, their approach
relied on pre-defined moving windows based on the spatial features of Iowa,
United States, limiting its generalizability to other cities. Bao et al.,
(2019) integrated CNNs, LSTM as well as ConvLSTM into a spatiotemporal
convolutional long short-term memory network (STCL-Net). Their work presented
a more practical and feasible approach for city-wide traffic risk prediction.
By combining the strengths of these models, the STCL-Net effectively captured
the spatiotemporal dependencies in the risk data and influential factors,
leading to improved prediction performance.

Furthermore, in order to capture both local and global dynamics with
hierarchical spatial information, Zhu et al., (2019) proposed the Deep
Spatiotemporal Attention Learning Framework. Unlike previous approaches that
primarily focused on local risk information, this framework incorporated a
spatiotemporal attention mechanism to effectively capture the dynamic impact
of traffic risk at different levels. One notable aspect of their work was the
emphasis on using real traffic administrative areas instead of manually
divided grids. This approach enabled more realistic modelling of spatial
heterogeneities in traffic incidents, thereby reducing prediction deviations
and enhancing the accuracy of risk prediction. In a similar vein, Huang et
al., (2019) proposed the deep Dynamic Fusion Network (DFN) framework, which
dynamically aggregated heterogeneous external factors. This framework allows
for the effective fusion of diverse data sources, enabling a comprehensive
understanding of the spatiotemporal dynamics of risk factors.

The utilization of Graph-based networks has emerged as a prominent approach in
addressing city-wide traffic risk prediction tasks by leveraging the spatial
graph features inherent in urban networks or regions. This approach takes
advantage of the spatial graph features derived from real urban elements,
enabling more effective modelling and prediction of traffic risk throughout a
city (Xue et al.,, 2022; Cheng et al.,, 2022; Li et al., 2022a, ; Li et al.,
2022b, ).

For instance, Zhou et al., (2020) constructed an urban graph based on urban
regions and proposed the RiskOracle framework, which integrates the
Differential Time-varying Graph Neural Network (DTGN). They learned dynamic
subregion-wise correlations between historical traffic status and traffic
risk, achieving joint prediction for traffic flow and risk. Notably, to
address the zero-inflation issue in the risk data, they employed an a priori
knowledge-based data enhancement (PKDE) strategy to define the subregional
risk labels during the training process. This strategy utilized two
experiencing parameters to transform zero risk values within the range of 0
and I while preserving the order of the real risk situation. Building upon
this work, Zhou et al., 2022b  enhanced the DTGN framework by incorporating an
LSTM-based hierarchical sequence learning structure. This enhancement
considered multi-scale historical spatiotemporal information for both long-
term and short-term predictions.

In a similar vein, Wang et al., 2021a  devised three distinct graph
structures: risk, road, and Points of Interest (POI), to capture the semantic
relationships among regions in their graph. They proposed two modules for risk
forecasting: a Spatiotemporal Geographical Module that captured the
geographical spatiotemporal correlations, and a spatiotemporal Semantic Module
that described the semantic spatiotemporal correlations based on the three
similarity graph structures. To address the zero-inflation issue, they
introduced a weighted loss function that focused on the non-zero risk values,
thereby mitigating the risk of overfitting. However, it is important to note
that this work primarily focuses on static geographic semantic similarities
and is not capable to deal with graph unaligned problems. In a similar vein,
Wang et al., 2021b  proposed a Multi-View Multi-Task Spatio-Temporal Networks
model for jointly predicting fine- and coarse-regional traffic incident risks.
Their approach integrated attention-based LSTM, channel-wise CNN, and multi-
view GCN along with multiple loss functions to facilitate multi-view
prediction and contain a mutual constraint for the finer-level loss. While
their work considered the zero-inflation issues in a cross-level corresponding
prediction task, it may not be an optimal solution when the primary focus is
on fine-level risk prediction. In contrast,Trirat et al., (2023) proposed a
multi-view graph neural networks framework with a multi-attention module to
incorporate dynamic similarity information from traffic status and also static
geographic semantic similarities. By integrating this dynamic information,
their approach provides a more comprehensive understanding of dangerous
driving behaviour and traffic incidents and enables a more accurate prediction
of urban traffic risk. The Huber loss function serves as a compromise between
squared error and absolute error losses, allowing for a robust regression
process that avoids being overly influenced by many zero values.

Recently, there have been a few works attempting to utilize graph-based deep
learning methods to address road-level risk prediction tasks. Yu et al.,
(2021) conducted road-level risk prediction using a novel Deep Spatiotemporal
Graph Convolutional Network. They defined three network layers, treating the
road graph as an independent layer, and incorporating heterogeneous
spatiotemporal data in the second Spatiotemporal Convolution Layer and an
Embedding Layer. To address the zero-inflation issue, they employed an
Undersampling method, which balanced the number of risky road samples and non-
risk road samples and treated incident risk as a binary classification task.
Although their approach demonstrated promising results at the road-level
granularity, the output only indicates the presence or absence of incidents,
thus masking the true risk situation of individual roads with limited sample
information. Wu et al., (2022) proposed a fully connected and dynamic network
framework for predicting road-level risk by capturing dynamic spatiotemporal
correlations. They addressed the zero-inflation issues by incorporating a
cost-sensitive learning module, which improved the classification accuracy of
positive samples without compromising the overall accuracy. This module
considered a probability fault cost based on the ratio of zero to non-zero
values, allowing for more constraints on non-zero risk value outputs.

Table 1: Summary of Deep Learning Methods in Spatiotemporal Prediction Reference |  Methodology |  Spatial Unit |  Zero-Inflate Issue
---|---|---|---
Chen et al., (2016) |  Stack Denoise Autoencoder |  $0.5km*0.5km$ Grid |  /
Ren et al., (2018) |  Long Short-term Memory network |  $1km*1km$ Grid |  /
Chen et al., (2018) |  Stack Denoise Convolutional Autoencoder |  Grid |  /
Yuan et al., (2018) |  Convolutional Long Short-Term Memory |  $5km*5km$ Grid |  /
Moosavi et al., (2019) |  Long Short-term Memory Network |  $5km*5km$ Grid |  /
Zhu et al., (2019) |  Deep Spatial-Temporal Attention Learning Framework |  Urban administrative areas |  /
Huang et al., (2019) |  Dynamic Fusion Network |  Grid |  /
Bao et al., (2019) |  Spatiotemporal Convolutional Long Short-term Memory Network |  $8*3$ Matrix Cells $15*5$ Matrix Cells $30*10$ Matrix Cells |  /
Zhou et al., (2020) |  Differential Timevarying Graph Neural Network |  $1.5km*1.5km$ Grid |  A priori knowledge-based data enhancement strategy
Wang et al., 2021a  |  Geographical and Semantic Spatiotemporal Network |  $2km*2km$ Grid |  Weighted Loss Function
Wang et al., 2021b  |  Multi-View Multi-Task Spatio-Temporal Networks |  $20*20$ Matrix Cells |  Cross-view mutual beneficial prediction
Yu et al., (2021) |  Deep Spatio-Temporal Graph Convolutional Network |  Road Segment |  Undersampling method
Wu et al., (2022) |  Multi-attention Dynamic Graph Convolution Network |  Road Segment |  Cost-sensitive Learning Loss Function
Zhou et al., 2022b  |  Differential Timevarying Graph Neural Network |  $1.5km*1.5km$ Grid |  A priori knowledge-based data enhancement strategy
Trirat et al., (2023) |  Multi-view Graph Neural Networks Framework with a Multi-attention Module |  Urban administrative areas |  Huber Loss Function

###  2.2 Uncertainty Quantification (UQ) in Spatiotemporal Prediction

With an increasing focus on solving a multitude of spatiotemporal prediction
problems, UQ has found broad applications in areas ranging from air pollution
(Wu et al.,, 2021) and traffic demand (Wang et al., 2023a, ; Zhuang et al.,,
2022) to traffic speed (Wu and James,, 2021), epidemic transitions (Lin et
al.,, 2021), and climate change (Thornton et al.,, 2021), among others.
Current UQ methods are generally bifurcated into two schools of thought: the
Bayesian and Frequentist perspectives.

In the Bayesian framework, the aim is centred around parameter inference,
grounded in Bayes’ theorem. This approach models the posterior distribution of
the network parameters, $p(\theta|x,y)$, for a given input-target pair,
$(x,y)$, predicated on the assumption of a known prior probability
distribution over these parameters, $p(\theta)$ (Gawlikowski et al.,, 2021).
According to comprehensive surveys by Gawlikowski et al., (2021) and He and
Jiang, (2023), primary methods for inferring the posterior distribution are
threefold: Laplace approximation based on the second derivative of likelihood,
Markov Chain Monte Carlo (MCMC) -based Sampling methods, and Variational
Inference (VI) predicated on a parametric distribution. Each of these methods
has found application in addressing Spatiotemporal Traffic-related UQ
problems. For instance, Qian et al., (2022) employed a Gaussian distribution
to encapsulate the aleatoric uncertainty, while combining VI and MC sampling
to estimate the epistemic uncertainty in the 5-minute traffic flow prediction
on an hourly basis. Kong et al., (2023) and Wu and James, (2021) incorporated
VI into GNN to quantify traffic condition uncertainty, specifically concerning
speed and flow. However, they mainly considered temporal information due to
computational complexity constraints. In contrast, Wu et al., (2021) utilized
the MCMC-based method to quantify uncertainty in limited node-restricted
forecasting tasks regarding air pollution, sensor-based traffic flow, and
state-based COVID death prediction. Liu et al., (2022) designed Bayes by
Backpropagation and local reparameterization techniques to capture the weight
distribution of models, thus achieving freeway incident detection.
Nevertheless, their method only provided short-term predictions of up to 30
minutes, owing to the disregard of spatial features. To summarize, the primary
challenge associated with the Bayesian approach lies in its computationally
intensive processing step, which poses significant scaling challenges for
finer-level spatiotemporal forecasting (Wang et al., 2023a, ; Zhuang et al.,,
2022; He and Jiang,, 2023). Moreover, the common approximation process shared
among the discussed Bayesian works contributes to an irreducible inaccuracy in
UQ (Gawlikowski et al.,, 2021; Wang and Yeung,, 2020).

In stark contrast to the Bayesian approach, the Frequentist perspective
primarily relies on Mean-Variance Estimation (MVE), enjoying the advantages of
computational efficiency and flexible model structure tuning with fewer prior-
knowledge constraints (Wu et al.,, 2023). Traditional models that produce
single-point outputs are supplanted by MVE’s ability to furnish a predictive
distribution. In this framework, it is postulated that the outputs are
consistent with a specific parameterized family of probability distributions,
typically falling within the Gaussian distribution (Chen et al.,, 2020; Wang
et al., 2023a, ). Therefore, MVE’s objective is the minimization of the loss
function associated with the parameters anticipated by the model.
Consequently, this technique provides a detailed representation of potential
outcomes, encouraging a deeper understanding of inherent uncertainties within
the predictions.

Two primary methods can be discerned within this perspective: parametric and
non-parametric. Parametric methods predict values based on the parameterized
distributions derived from the parameters of the neural networks (Wang et al.,
2023a, ; Khosravi and Nahavandi,, 2014). In contrast, non-parametric methods
concentrate on establishing the upper and lower bounds or employ quantile
regression to obtain the prediction intervals (PIs) (Pearce et al.,, 2018;
Khosravi et al.,, 2010; Koenker and Hallock,, 2001; He and Jiang,, 2023;
Gawlikowski et al.,, 2021).

For instance, both Guo et al., (2014) and Lin et al., (2018) applied non-
parametric methods to ascertain the upper and lower bounds of short-term
traffic flow predictions. Guo et al., (2014) employed an adaptive Kalman
filter to dynamically map volatile traffic flow, thus achieving real-time UQ.
Lin et al., (2018) combined Particle Swarm Optimization (PSO) with the Extreme
Learning Machine (ELM) neural network to predict PIs, leading to a more
precise individual assessment of upper and lower bounds. In a different
approach, Li et al., 2022c  adapted the traditional ConvLSTM model,
integrating it with a differentiable PI-oriented loss function to attain long-
term highway traffic flow predictions and UQ encapsulation. Nevertheless,
these works, primarily concerned with traffic flow prediction, tend to
approach UQ from a temporal series perspective, neglecting spatial UQ.
Moreover, Zhou et al., (2021) introduced a novel spatiotemporal UQ framework
that accounted for differences in regional traffic situations and a
hierarchical data turbulence scheme. They utilized an out-of-distribution
noise parameter to encapsulate the predicted PIs. Despite promising results,
their model depends on dense data input. While non-parametric methods are not
as restricted by prior data knowledge, they are heavily reliant on abundant
data input and advanced method structures, which reduces their applicability
in complex urban environments.

On the other hand, parametric methodologies typically introduce auxiliary
probabilistic layers to seize pertinent parameters, thus achieving minimal
alterations to the existing training structures. UQ is calculated in this
approach via parameter-driven intervals, in tandem with mean values and
variance. Wang et al., 2023a  offered a thorough examination of diverse
Gaussian-based probabilistic suppositions incorporated with GNNs, thereby
formulating a generalized framework for quantifying spatiotemporal uncertainty
in travel demand. They proposed that the adoption of suitable probabilistic
suppositions is pivotal for achieving stable UQ while exerting minimal impact
on single-point estimation. As noted by Wang et al., 2023a , these strategies
are heavily dependent on predefined data distributions. While parametric
methodologies for UQ present computational efficiency and facile adaptability,
exploration concerning sparse and event-centric traffic data remains
conspicuously scarce. The most recent pioneering work in this field is Zhuang
et al., (2022), who were among the first to address the issue of sparse data
in the context of UQ. Instead of using the Gaussian assumption, they
introduced a three-parameter zero-inflated negative binomial distribution
combined with spatiotemporal neural networks to simultaneously manage the
predominant zero instances and non-normal distribution in sparse travel demand
data. Despite the commendable merits of their innovative method and impressive
results, they encountered practical challenges when zero values were overly
extreme and in non-time series traffic scenarios, particularly in event-based
traffic risk situations.

##  3 Methodology

In this section, we start by introducing the important concepts and defining
the framework for the real-world problem of road-level traffic risk prediction
under study in Section 3.1. We then illustrate the progression from the TD
distribution to ZITD within the context of UQ in Section 3.2. Subsequently, we
introduce our novel STZITD-GNNs model, a spatiotemporal probabilistic
modelling framework designed for ZITD-based UQ, that combines the strengths of
the GRU and GAT in Section 3.3. The section finishes with an in-depth
exploration of the loss function employed within our model in Section 3.4. The
major notations used in this section are listed in Table 2.

Table 2: Main notations in this paper Notation | Description
---|---
$Y\leavevmode\nobreak\ /\leavevmode\nobreak\ X$ | Road-level traffic risk / Road-level spatiotemporal traffic feature
$V\leavevmode\nobreak\ /\leavevmode\nobreak\ E\leavevmode\nobreak\ /\leavevmode\nobreak\ A$ | Road set / Edge set / Road adjacency matrix
$v_{i}\leavevmode\nobreak\ /\leavevmode\nobreak\ (v_{i},v_{j})$ |  The $i^{th}$ road / Edge between the $i^{th}$ road and the $j^{th}$.
$N\leavevmode\nobreak\ /\leavevmode\nobreak\ T$ | Numbers of roads / Time length
$C_{it}$ |  Number of incidents on road $v_{i}$ at time slot $t$
$l_{it}^{j}$ |  The $j^{th}$ cardinality values of incident severities at time slot $t$ for road $v_{i}$
$d\leavevmode\nobreak\ /\leavevmode\nobreak\ F\leavevmode\nobreak\ /\leavevmode\nobreak\ F^{\prime}$ | Feature dimension/Temporal feature dimension/Spatio-Temporal feature dimension
$\mathcal{ST}$ | Spatiotemporal deep learning Model
$f_{\text{TD}}(\cdot)$ | Probabilistic density function of Tweedie Distribution
$\rho/\leavevmode\nobreak\ \phi/\leavevmode\nobreak\ \mu/\leavevmode\nobreak\ \pi$ | Index parameter/ Dispersion parameter/ Mean parameter/ Zero-inflated parameter
$\theta\leavevmode\nobreak\ /\leavevmode\nobreak\ \lambda\leavevmode\nobreak\ /\leavevmode\nobreak\ \alpha,\gamma$ | Natural parameter / Poisson mean parameter / Gamma parameters
$\mathcal{Z}_{T}\leavevmode\nobreak\ /\leavevmode\nobreak\ \mathcal{Z}$ | Temporal embedding / Spatiotemporal embedding
$\text{ReLU}(\cdot)/\leavevmode\nobreak\ \sigma(\cdot)/\leavevmode\nobreak\ \text{LeakyReLU}(\cdot)/\leavevmode\nobreak\ [\cdot||\cdot]$ | ReLU activation / Sigmoid activation / LeakyReLU activation / Concatenation
$W_{r},W_{u},W_{c},W_{a},b_{r},b_{u},b_{c},\vec{a}$ | Learnable weight matrics of GRU and GAT
$W_{\pi},W_{\mu},W_{\phi},W_{\rho},b_{\pi},b_{\mu},b_{\phi},b_{\rho}$ | Learnable weight matrics of the four Paramater Encoders
$h_{t}\leavevmode\nobreak\ /\leavevmode\nobreak\ r_{t}\leavevmode\nobreak\ /\leavevmode\nobreak\ u_{t}$ | Hidden feature / Reset gate / Update gate
$\alpha_{i,j}$ |  Attention value between road $v_{i}$ and $v_{j}$
$M$ | Number of attention heads
$\epsilon$ | The minimum value
$\eta$ | L2 normalization weight-parameter

###  3.1 Preliminaries and Problem Definition

This section lays out the foundational concepts and terminology pertinent to
the study, followed by a precise statement of the problem at hand.

Definition 1. Road Traffic Incident Risk: A traffic incident refers to an
incident involving one or more road users, leading to physical injury, loss of
life, or damage to property. We denote the set of roads as $V$ of size
$N=|V|$, where $v_{i}$ signifies the $i^{th}$ road. Traffic risk on roads is
represented as $Y\in\mathbb{R}^{N\times T}$ for a duration of $T$ minutes, and
$y_{it}$ denotes the traffic risk at the $t^{th}$ time slot for road $v_{i}$.
Thus, $Y_{t}\in\mathbb{N}^{N}$ signifies the traffic risk for all roads at the
$t^{th}$ time slot, where $y_{it}$ is its constituent element. The computation
of traffic risk $y_{it}$ for road $v_{i}$ at time slot $t$ is given by:

$\displaystyle y_{it}=\sum_{j=1}^{C_{it}}l_{it}^{(j)},$ (1)

where $C_{it}$ denotes the number of incidents, and $l_{it}^{(j)}$ represents
the $j^{th}$ severity level of incidents at time slot $t$ for road $v_{i}$. We
allocate the incident point to its closest road. As per prior research (Wang
et al., 2021a, ; Trirat et al.,, 2023), $l$ is assigned the values 1, 2, and
3, representing minor, injury, and fatal incident severities, respectively.
$y_{it}$ equals 0 if there are no traffic incidents, i.e., $C_{it}=0$.

Definition 2. Spatiotemporal Features: Road features are represented as
$X\in\mathbb{R}^{N\times T\times d}$, where $d$ is the dimension of the
spatiotemporal road features. We denote $x_{it}$ as the spatiotemporal feature
of road $v_{i}$ at the $t^{th}$ time slot. Accordingly,
$X_{t}\in\mathbb{R}^{N\times d}$ signifies the spatiotemporal feature of all
roads at time slot $t$. Our study includes spatial features like geographical
information, general road conditions, and census features. Temporal features,
which vary daily, encompass weather information such as sunrise and sunset
times, humidity, visibility, rainfall, etc., along with the designation of the
day as a holiday or a working day.

Definition 3. Road Connection Graph: Keeping the connectivity of roads in
mind, we construct the road network graph $\mathcal{G}=(V,E,A)$. Here, $E$
denotes the edge set, and $A\in\mathbb{R}^{N\times N}$ is the adjacency matrix
that describes the connections between roads. An entry $A_{i,j}=1$ indicates
the existence of an edge between road $v_{i}$ and $v_{j}$, denoted as
$(v_{i},v_{j})\in{E}$, while $A_{i,j}=0$ implies no connection, denoted as
$(v_{i},v_{j})\notin{E}$.

In line with Tobler’s First Law of Geography (Tobler,, 1970), spatial
correlation bears a significant relationship with geographic distance (Yu et
al., 2020a, ; Zhou et al.,, 2021; Wang et al., 2023b, ). This implies that the
geographical proximity between two roads directly influences the correlation
of traffic risk effects: the closer the roads, the stronger the correlation.
As a result, traffic incidents often lead to spatial spillover effects,
suggesting that a risk point on one road can have a significant impact on its
neighbouring roads. To encapsulate this spatial attribute, we employ random
walks (Perozzi et al.,, 2014; Li et al.,, 2017; Hamilton,, 2020) as a
simulation process. In this procedure, we use the adjacency matrix $A$ to
disperse road traffic risks to their closest neighbours, adhering to three
predefined thresholds: a diffusion probability, a road-level traffic risk
reduction factor, and a maximum hop count.

Problem Definition: The aim of our model is to predict the future expected
traffic incident risk and the confidence interval of predicted results per
road, with a forecast horizon of $p$ time windows, based on $N$ roads and
their traffic risk across time windows of length $T$ minutes. Essentially,
this is a sequence-to-sequence prediction task. The objective is to utilize
historical records $X_{1:t},Y_{1:t}$ and the graph structure $A$ as input data
for training, to predict the probabilistic density function $f(Y_{t+1:t+p})$
of the distribution of $Y_{t+1:t+p}$ (i.e., the traffic risk distribution for
the next $p$ time windows). This facilitates the analysis of future road-level
traffic risk expectations and confidence intervals:

$\mathcal{ST}_{\Theta}([X_{1:t},Y_{1:t}],A)\to f(Y_{t+1:t+p}),$ (2)

where $\Theta$ denotes the learned parameters for the ST-GNN model (also
referred to as $\mathcal{ST}$). The predictions are made for both a 14-day
input to a 14-day output and a 7-day input to a 7-day output.

Figure 1: The overall framework of STZITD-GNNs. STZITD-GNNs utilize ST-GNN
Encoder $\mathcal{ST}$ (composed by GRU and GAT encoders) to encode the
history time window 1:t incident risk $Y_{1:t}$ and road spatiotemporal
features $X_{1:t}$ with use of road connection graph into road spatialtemporal
embedding $\mathcal{Z}$ at first. After encoding, the four parameter decoders
map $\mathcal{Z}$ into ZITD parameter space and obtain
$\pi_{t+1:t+p},\mu_{t+1:t+p},\phi_{t+1:t+p},\rho_{t+1:t+p}$ for predicted time
window $t+1:t+p$, which determine the predict feature road traffic risk
distribution $f(Y_{t+1:t+p})$.

###  3.2 From TD Distribution to ZITD Distribution

This subsection provides an in-depth exploration of the progression from TD to
ZITD distribution in the specific context of road risk prediction. Building on
the work of Lord and Mannering, (2010) and Smyth and Jørgensen, (2002), we
recognize the effectiveness of TD distribution in addressing over-dispersed
data, indicating its potential applicability in road-level traffic risk
prediction.

As illustrated in Section 3.1, when we aggregate risk data (incorporating
various severities of crashes) from different road segments or identical
segments over different time intervals, the resulting distribution of risk
counts tends to exhibit over-dispersion of varying degrees. In these
scenarios, the variance of risk counts outstrips the mean risk count, creating
a deviation from the Gaussian distribution, as underlined by the work of
Halder et al., (2019) and Wang et al., 2023a .

In detail, a defining feature of overdispersed risk data—applying to both
incident counts and the associated severity—is the prevalence of ’excess’
zeros, implying more zero entries than would be predicted under a
corresponding Poisson process. Conversely, we observe a skewness in non-zero
outcomes, specifically a left-skewed distribution representing a smaller
proportion of higher risk values. This unique pattern, combined with the long-
tail distribution of risk data, necessitates a robust modelling approach. In
transitioning from the TD to the ZITD distribution, our model accommodates
these intrinsic traits of overdispersed and long-tail data. It not only
addresses the issue of excess zero counts innovatively but also accurately
captures the left-skewness of non-zero outcomes, making it a more fitting tool
for road-level traffic risk prediction.

####  3.2.1 TD Distribution

In adherence to Equation 1, we designate $y_{k}$ as the risk value of road $i$
on day $t$. As a special case of Exponential Dispersion Models (EDMs), the TD
model encompasses a broad family of statistical distributions (Tweedie et
al.,, 1984; Jørgensen,, 1987; Jiang et al.,, 2023). A random variable
following the TD distribution possesses a probability density function,
$f_{TD}$, with its probability density function computed as follows
(Jørgensen,, 1987):

$f_{\text{TD}}(y_{k}|\theta,\phi)\equiv
a(y_{k},\phi)\exp{\bigl{[}\frac{y_{k}(\theta)-\kappa(\theta)}{\phi}\bigl{]}},$
(3)

here, $\theta\in\mathbb{R}$ represents the natural parameter, while
$\phi\in\mathbb{R}^{+}$ stands for the dispersion parameter. The normalizing
functions $a(\cdot)$ and $\kappa(\cdot)$ are parameter functions for $\phi$
and $\theta$, respectively, detailed further in the subsequent discussion. For
a TD distribution of EDMs, the mean and variance of variable $y_{k}$ are given
by:

$E(y_{k})=\mu=\kappa(\theta)^{{}^{\prime}},\quad
Var(y_{k})=\phi\kappa(\theta)^{{}^{\prime\prime}},$ (4)

In these equations, $\kappa(\theta)^{{}^{\prime}}$ and
$\kappa(\theta)^{{}^{\prime\prime}}$ represent the first and second
derivatives of $\kappa(\theta)$, respectively. Meanwhile, $\mu\geq 0$
functions work as the mean parameter. Since the mean $E(y_{k})$ and variance
$Var(y_{k})$ are dependent, $\kappa(\theta)^{{}^{\prime\prime}}$ can be
directly expressed as $V(\mu)$, thereby providing a direct mean-variance
relationship (Mallick et al.,, 2022). Additionally, the TD family encompasses
a number of vital distributions for distinct index parameters $\rho$ as
demonstrated by $V(\mu)=\mu^{\rho}$. These include the Normal ($\rho=0$),
Poisson ($\rho=1$), Gamma ($\rho=2$), Inverse Gaussian ($\rho=3$), and
Compound Poisson-Gamma distribution ($1<\rho<2$) and could be refer to A. The
Compound Poisson-Gamma distribution is particularly applicable, given its
practical prowess in parameterizing zero-inflated and long-tailed data, which
aligns with the observed characteristics of real traffic risk data scenarios.

The frequency of traffic risk $C_{k}$ adheres to a Poisson distribution
represented by $Pois(\lambda)$ with a mean of $\lambda$, expressed as
$C_{k}\sim Pois(\lambda)$. This is natural as the frequency often appears as
the main prediction target in the previous literature (Zhou et al.,, 2020;
Trirat et al.,, 2023). Concurrently, the incident severity $l_{k}^{(j)}$ are
modelled as independent and identically distributed $i.i.d$ Gamma random
variables, $C_{k}\perp l_{k}^{(j)}$ and $\forall j$, denoted as
$Gamma(\alpha,\gamma)$. These variables bear a mean of $\alpha\gamma$ and a
variance of $\alpha\gamma^{2}$, as $l_{k}^{(j)}\stackrel{{\scriptstyle
iid}}{{\sim}}\operatorname{Gamma}(\alpha,\gamma)$. $C_{k}$ and $l_{k}^{(j)}$
together form the modelling of the occurrence of a single incident.
Consequently, the traffic incident risk $y_{k}$ is formulated as a Poisson sum
of these iid Gamma random variables. Corresponding to Equation 1 and
considering $1<\rho<2$, we can finally reformulate the random variable $y_{k}$
as follows:

$y_{k}=\left\\{\begin{array}[]{cc}0&\text{if }C_{k}=0,\\\
\sum_{j=1}^{C_{k}}l_{k}^{(j)}=l_{k}^{(1)}+l_{k}^{(2)}+\cdots+l_{k}^{(C_{k})}&\text{if
}C_{k}>0.\end{array}\right..$ (5)

In the scenario where $C_{k}=0$, we correspondingly have $y_{k}=0$. The
probability density at zero for traffic incident risk, therefore, becomes
$\mathbb{P}(y_{k}=0)=\mathbb{P}(C_{k}=0)=\exp(-\lambda)$. Conversely, when
$C_{k}>0$, $y_{k}$ is deduced from the sum of $C_{k}$ independent Gamma random
variables. For clarity, following the work of Zhou et al., 2022a , we re-
parameterize the TD distribution, setting $\theta=\mu^{1-\rho}/(1-\rho)$ and
$\kappa(\theta)=\mu^{2-\rho}/(2-\rho)$. This results in the revised TD
distribution as:

$f_{\text{TD}}(y_{k}|\theta,\phi)=f_{\text{TD}}(y_{k}|\mu,\phi,\rho)=a(y_{k},\phi,\rho)\exp{\bigl{[}\frac{1}{\phi}({y_{k}\frac{\mu^{1-\rho}}{1-\rho}-\frac{\mu^{2-\rho}}{2-\rho}})\bigl{]}},$
(6)

where

$a(y_{k},\phi,\rho)=\left\\{\begin{array}[]{cc}1&\text{if }y_{k}=0,\\\
\frac{1}{y_{k}}\sum_{j=1}^{\infty}\frac{y_{k}^{-j\alpha}(\rho-1)^{\alpha
j}}{\phi^{j(1-\alpha)}(2-\rho)^{j}j!\Gamma(-j\alpha)}&\text{if
}y_{k}>0.\end{array}\right..$ (7)

In this case, $\mu,\phi$, and $\rho$ are the three parameters governing the
probability and expected value of traffic incident risk. The parameters in
Gamma and Poisson distributions ($\lambda,\alpha,\gamma$) can be obtained as
(Kurz,, 2017; Halder et al.,, 2019; Mallick et al.,, 2022), :

$\lambda=\frac{1}{\phi}\frac{\mu^{2-\rho}}{2-\rho},\quad\alpha=\frac{{2-\rho}}{\rho-1},\quad\gamma=\phi(\rho-1)\mu^{\rho-1}.$
(8)

Generally, the TD is favoured because of its congruity with the nature of road
risk, defined as an expectation derived from three categories of severity and
their associated frequencies. This echoes the empirical observation that
traffic incident frequency aligns with a Poisson distribution, while incident
severity follows a Gamma distribution characterized by zero inflation and long
tails. A significant departure from previous studies such as STZINB-GNN, our
approach to the TD distribution innovatively parameterizes $\phi$ to model the
left-skewed zero- or minor-risk situations, and $\rho$ to address the issues
related to the long tail. This nuanced approach enables a more realistic
modelling of the incident risk landscape, providing improved prediction
accuracy and a better understanding of the complex dynamics of road risk.

####  3.2.2 ZITD Distribution

Traffic risk values often exhibit a characteristic abundance of zero values,
as depicted in the numerical analysis. This surplus often surpasses the zero
input expected in conventional models, thus risking the oversight of many
hazardous road situations due to methodological constraints. Likewise, while
the standard TD model proposed in section 3.2.1 accommodates a positive
probability density at zero, as illustrated by
$\mathbb{P}(y_{k}=0)=\mathbb{P}(C_{k}=0)=\exp(-\lambda)$, it encounters
performance shortfalls when the empirical distribution of risk data at a road-
level is exceedingly imbalanced and dominated by zero risk values (Mallick et
al.,, 2022; Zhou et al., 2022a, ).

To address this issue, we introduce a refinement to the model, known as the
Zero-Inflated Four-parametric Tweedie Compound (ZITD) model. To our knowledge,
this marks the inaugural application of a ZITD distribution along with its
estimation framework for the traffic incident domain, specifically tailored
for the characterization of road-level traffic risk. In the construction of
the ZITD distribution, we integrate a sparsity parameter, denoted as $\pi$, to
account for the skewness inherent in risk data, often manifested as an excess
of zeros. This sparsity parameter $\pi$ is utilized to denote a zero-inflation
scenario, while the probability $1-\pi$ corresponds to the TD distribution:

$y_{k}=\left\\{\begin{array}[]{cl}0&\text{with probability }\pi\\\
Y&\text{with probability }1-\pi,Y\sim
f_{TD}(Y|\mu,\phi,\rho)\end{array}\right..$ (9)

Finally, by denoting the ZITD model as $y_{k}\sim
f_{ZITD}(\pi,\mu,\phi,\rho)$, we formulate the probability density function of
this model as follows:

$f_{ZITD}(y_{k}|\pi,\mu,\phi,\rho)=\left\\{\begin{array}[]{cc}\pi+(1-\pi)f_{TD}(y_{k}=0|\mu,\phi,\rho)&\text{if
}y_{k}=0\\\ (1-\pi)f_{TD}(y_{k}|\mu,\phi,\rho)&\text{if
}y_{k}>0\end{array}\right..$ (10)

In accordance with this specification, we tackle the problem of extreme zero
inflation by applying a linear weighting to the zero value, such that
$\mathbb{P}(y_{k}=0)=\pi+(1-\pi)\exp{(-\frac{1}{\phi}\frac{\mu^{2-\rho}}{2-\rho})}$.
It’s important to note that this probability is significantly higher than the
probability of observing a zero count from the TD distribution. Its increased
probability intuitively signifies zero inflation. This modification ensures
that the mean value of the ZITD distribution is
$\mathbb{E}(y_{k})=(1-\pi)\mu$.

In the proposed ZITD model, each of the four parameters plays a critical role
in managing various facets of road-level incident risk values. Importantly,
the sparsity parameter, $\pi$, modulates the weights for situations devoid of
risk occurrence. This modulation leads to a significant reformation of the
dispersion parameter, $\phi$, which is now interpreted as weighted, thereby
refining the model’s approach towards left-skewed scenarios. In this capacity,
larger $\phi$ values symbolize an enhanced model ability to simulate instances
predominantly characterized by none- or lower-risk situations. Complementing
these, the index parameter $\rho$ becomes instrumental for scenarios
characterized by a long-tail distribution of road-level risk values. Higher
$\rho$ values, particularly those closer to 2, represent an enhanced
capability to accommodate such instances, distinguishing this model from
others that predominantly focus on mean values, represented by $\mu$.

These parameters improve the model’s representational capacity and deepen
understanding of uncertainty quantification for predicted outcomes. By moving
beyond point estimation to interval estimation, the ZITD model marks a
significant advancement in road risk representation and prediction.

###  3.3 The framework of Spatiotemporal Zero-Inflated Tweedie Graph Neural
Networks (STZITD-GNNs)

In relation to the four parameters $\pi,\mu,\phi,\rho$ of the ZITD probability
distributions detailed above, we have designed a probabilistic spatiotemporal
graph learning module, named STZITD-GNNs. This module extracts the
spatiotemporal correlations
$\mathcal{Z}=\mathcal{ST}_{\Theta}(X_{1:t},Y_{1:t},A)$ under the ZITD
assumption for each road’s risk and predicts future outcomes for the
subsequent $p$ time windows through $f_{ZITD}(X_{t+1:t+p}|\mathcal{Z})$. As
described in Equation 11, the comprehensive probabilistic Spatiotemporal Graph
Neural Network is denoted as $\mathcal{ST}_{\Theta}$. It jointly takes
historical spatiotemporal features, corresponding incident risk values, and
the adjacency matrix as inputs, learning the parameter embedding $\mathcal{Z}$
for future risk predictions.

$\displaystyle
f_{ZITD}(Y_{t+1:t+p}|\pi_{t+1:t+p},\mu_{t+1:t+p},\phi_{t+1:t+p},\rho_{t+1:t+p})$
(11)
$\displaystyle=f_{ZITD}(Y_{t+1:t+p}|\mathcal{ST}_{\Theta}(X_{1:t},Y_{1:t},A))=f_{ZITD}(X_{t+1:t+p}|\mathcal{Z}).$

As illustrated in Figure 1, instead of conducting two independent estimations,
our methodology commences with the deployment of a GRU acting as the temporal
encoder. Subsequently, a GAT is implemented for the purpose of spatial
encoding. The combined output from the GRU and GAT is utilized to parameterize
the ZITD distribution. The proposed unified approach is formulated as
showcased in Equation 12. Our architecture’s encoding component applies a
spatiotemporal embedding scheme complemented with an extra sparsity parameter
and distribution parameters. On the other hand, the decoding component is
responsible for probabilistic estimations of future risk values, which can be
construed as the parameter set of the future risk distribution. In contrast
with the principles of variational autoencoders (VAEs) (Kipf and Welling,
2016b ; Hamilton, (2020)), which constraint the latent embedding with specific
distributions such as Gaussian and t- distributions for a more robust
representation of latent variables, we echo the statistical domain by allowing
the integration of deep learning methodologies into parameters learning:

$\pi_{t+1:t+p},\mu_{t+1:t+p},\phi_{t+1:t+p},\rho_{t+1:t+p}=\mathcal{ST}_{\Theta}(X_{1:t},Y_{1:t},A).$
(12)

In the following sections, we present our STZITD-GNNs learning framework in
greater detail. Firstly, we introduce the temporal encoder, GRU, and then
delve into the crucial aspects of the spatial encoder, GAT.

####  3.3.1 Temporal encoder – GRU

Traffic incidents demonstrate significant temporal correlation, underscoring
the idea that the dynamic temporal features of roads encapsulate a wealth of
information. This observation has been acknowledged and exploited in prior
methodologies (Wang et al., 2021a, ; Zhao et al.,, 2019). As shown in Figure
1, to seize these temporal proximity features, we implement a GRU (Chung et
al.,, 2014). The GRU, in comparison to traditional LSTM networks, offers the
advantage of learning relatively long-term dependencies without facing issues
of vanishing and exploding gradients. In the GRU framework, there isn’t a
segregation of internal and external states, as seen in LSTM networks.
Instead, it addresses the gradient vanishing and exploding problem by directly
integrating a linear dependency between the current network state $h_{t}$ and
the preceding state $h_{t-1}$. Consequently, while a GRU preserves the
capabilities of an LSTM network, its architecture remains more streamlined (Gu
et al., (2019); Mahmoud et al., (2021); Ma et al., (2022)).

Specifically, this component accommodates the feature sequence over past time
windows $1:t$, thereby learning temporal dependencies. We represent the
temporal embedding as $\mathcal{Z}_{T}$, taking into account both historical
spatiotemporal road features and traffic risks as input:

$\displaystyle\mathcal{Z}_{T}=\text{GRU}(X_{1:t},Y_{1:t}),$ (13)

Once the input features have been processed via the GRU, we obtain historical
temporal embeddings $\mathcal{Z}T$ for all roads. These are subsequently
directed to graph neural encoders. The GRU, in particular, consists of two
gates: the update gate and the reset gate. The encoder processes a single
input $X_{t},Y_{t}$ at time slot $t$ and prior hidden features $h_{t-1}$ in
sequence:

$\displaystyle r_{t}=\sigma(W_{r}\cdot[h_{t-1},X_{t},Y_{t}]+b_{r}),$ (14)
$\displaystyle u_{t}=\sigma(W_{u}\cdot[h_{t-1},X_{t},Y_{t}]+b_{u}),$
$\displaystyle\tilde{h}_{t}=\text{tanh}(W_{c}\cdot[r_{t}\cdot
h_{t-1},X_{t},Y_{t}]+b_{c}),$ $\displaystyle h_{t}=(1-u_{t})\cdot
h_{t-1}+u_{t}\cdot\tilde{h}_{t}$

In the equation above, $h_{t}$ and $h_{t-1}$ represent hidden features at the
current and preceding time steps, respectively, while $\tilde{h}_{t}$ denotes
the candidate state that temporarily stores information from reset and update
gates ($r_{t}$,$u_{t}$). The variables $W$ and $b$ (i.e.
$W_{r},W_{u},W_{c},b_{r}$, $b_{u},b_{c}$) represent learnable weight and bias
matrices, respectively. Importantly, $r_{t}$ acts as the reset gate,
determining the extent of integration of hidden features from previous time
steps into features at the current time step. Meanwhile, $u_{t}$ serves as the
update gate, controlling how much information from previous time steps will be
carried over to the next time step. The $\text{tanh}(\cdot)$ is the hyperbolic
tangent function. We perform the GRU gate operations step by step and select
the final layer’s hidden features $h_{t}$ as the temporal embedding
$\mathcal{Z}_{T}\in\mathbb{R}^{N\times F}$, where $F$ denotes the dimension of
road temporal embedding.

####  3.3.2 Spatial encoder – GAT

Within our research, road data is fundamentally perceived as having an
inherent graph structure. This perspective facilitates the interpretation of
node features as graph signals. To accommodate such graph-structured
information, we employ the GAT (Veličković et al.,, 2017), an effective
methodology renowned for enabling potential model deepening while concurrently
incorporating dynamic spatial correlation into our network, achieved through
real-time adjustment of attention weights in response to data changes (Wang et
al., 2022b, ; Zhang et al., 2019a, ).

Distinct from GCN of a spectral-domain convolution (Kipf and Welling, 2016a,
), GAT operates as a spatial-domain convolutional network. Instead, GAT
emphasizes the magnitude of the graph signal value at distinct points and
their inter-distance to execute the convolution process. This trait empowers
GAT to alleviate the over-smoothing issue associated with GCN deepening,
thereby augmenting prediction precision regarding peaks and valleys (Zhao et
al.,, 2019).

Initially, we compute the attention coefficient $\alpha_{ij}$ between roads
$v_{i}$ and $v_{j}$ as per the following equation:

$\displaystyle\alpha_{ij}=\frac{\exp({\text{LeakyReLU}\bigl{(}\vec{a}^{\mathsf{T}}[W_{a}\mathcal{Z}{T_{i}}||W_{a}{\mathcal{Z}{T_{j}}])\bigl{)}}}}{\sum_{A_{i,o}>0}\exp({\text{LeakyReLU}\bigl{(}\vec{a}^{\mathsf{T}}[W_{a}{\mathcal{Z}{T_{i}}}||W_{a}{\mathcal{Z}{T_{o}}}])\bigl{)}}},$
(15)

In this equation, $A_{i,o}>0$ represents the neighbourhood or connected roads
of the road $v_{i}$. $W_{a}\in\mathbb{R}^{F\times F^{\prime}}$ and
$\vec{a}\in\mathbb{R}^{2F^{\prime}}$ are learnable weights, where $F^{\prime}$
denotes the hidden dimension. We apply the concatenation function
$[\cdot||\cdot]$ to concatenate the embedding of $v_{i}$ and its neighbours,
utilizing learnable parameter $W_{a}$ and mapping function $\vec{a}$.
Subsequently, we employ an activation function, $\text{LeakyReLU}(\cdot)$, and
normalize the attention coefficient $\alpha_{ij}$ using the $\exp(\cdot)$
function. The embedding of road $v_{i}$ is consequently updated with a linear
combination of the neighbourhood of road $v_{i}$ and itself, wherein these
computed attention coefficients are employed.

In order to preemptively address potential inaccuracies in attention
assignments, we have incorporated multiple attention assignments, otherwise
known as multiple attention heads, into our approach. This tactic greatly
amplifies the expressive capacity of the model. The multi-head attention
mechanism facilitates concurrent consideration of information from a wide
range of representation subspaces at various positions. This enables each head
to learn and focus on different facets of the input data, thereby enhancing
the model’s ability to comprehend and accurately represent the data (Reza et
al.,, 2022; Tang and Zeng,, 2022; Ding et al.,, 2023). Herein, $M$ symbolizes
the number of attention mechanisms, with the output embeddings from these
mechanisms being concatenated in the following manner:

$\displaystyle\mathcal{Z}_{i}=||_{m=1}^{M}\sigma\left(\sum_{A_{i,j}>0}\alpha_{ij}^{m}{W_{a}}^{m}\mathcal{Z}_{T_{j}}\right),$
(16)

In this equation, superscript $m$ denotes a different head. The embedding of
the road $v_{i}$ is optimized by aggregating its neighbours’ embedding with
attention coefficients $\alpha_{ij}^{m}$, and we apply learnable weights
${W_{a}}^{m}$ for different heads. After this, we apply a nonlinear sigmoid
function and concatenate different heads of embedding. The multi-head
attention mechanisms hereby ensure stability in the attention mechanisms
learning process and encapsulate diverse perspectives of road-level graph
features. Notably, in Eq 16, the dimension of the output feature is $M\times
F^{\prime}$, rather than $F^{\prime}$. To ensure consistency in the output
feature’s dimension, we compute the average of each road’s embedding in the
last layer as follows:

$\displaystyle\mathcal{Z}_{i}=\sigma\left(\frac{1}{M}\sum_{m=1}^{M}\sum_{A_{i,j}>0}\alpha_{ij}^{m}{W_{a}}^{m}\mathcal{Z}_{T_{j}}\right).$
(17)

Finally, we input temporal embedding $\mathcal{Z}_{T}$ into GAT to capture
spatial dependency and obtain spatiotemporal
$\mathcal{Z}_{i}\in\mathbb{R}^{F^{\prime}}$ of road $v_{i}$. After that, we
use $\mathcal{Z}$ to infer the four-parameter STZITD-GNNs model.

####  3.3.3 Parameter Decoders

To transmute the acquired spatiotemporal embeddings into ZITD parameter
values, we have devised four parameter decoders to operate on $\mathcal{Z}$.
Consequently, four parameters, namely
$\pi_{t+1:t+p},\mu_{t+1:t+p},\phi_{t+1:t+p},$ and $\rho_{t+1:t+p}$ can be
computed from $\mathcal{Z}$, as detailed below:

$\displaystyle\pi_{t+1:t+p}$
$\displaystyle=\sigma(W_{\pi}\cdot\mathcal{Z}+b_{\pi})$ (18)
$\displaystyle\mu_{t+1:t+p}$
$\displaystyle=\text{ReLU}(W_{\mu}\cdot\mathcal{Z}+b_{\mu})$
$\displaystyle\phi_{t+1:t+p}$
$\displaystyle=\text{ReLU}(W_{\phi}\cdot\mathcal{Z}+b_{\phi})+\epsilon$
$\displaystyle\rho_{t+1:t+p}$
$\displaystyle=\sigma(W_{\rho}\cdot\mathcal{Z}+b_{\rho})+1+\epsilon$

Here, $\pi$ lies within the range of $[0,1]$, $\mu$ falls within the interval
$[0,+\infty)$, $\phi$ exists within $(0,+\infty)$, and $\rho$ spans the range
of $(1,2)$. The learnable weight matrices are represented as
$W_{\pi},W_{\mu},W_{\phi},W_{\rho},b_{\pi},b_{\mu},b_{\phi},b_{\rho}\in\mathbb{R}^{F^{\prime}\times
p}$. $\lim\epsilon\to 0$ represents the smallest value, and
$\text{ReLU}(\cdot)$ signifies the ReLU activation function.

###  3.4 Learning Framework of STZITD-GNNs

####  3.4.1 ZITD Loss Function

As previously mentioned, our encoding mechanism employs a spatiotemporal
embedding architecture complemented by a sparsity parameter. On the other
hand, the decoding component involves probabilistic estimation of future road
traffic risks. To accurately predict these risks, the overall learning
objective can be seen as maximizing the log-likelihood function: $\max\log
f_{ZITD}(y_{k}|\pi,\mu,\phi,\rho)$. We utilize the negative likelihood as our
loss function to better fit the distribution into the data. Here, $y$ denotes
the ground-truth values corresponding to one of the predicted risks with
parameters $\pi,\mu,\phi,\rho$. The log-likelihood of ZITD is composed of the
$y=0$ and $y>0$ parts:

For $y>0$, $NLL_{y>0}=-\log f_{ZITD}(y>0|\pi,\mu,\phi,\rho)$:

$\displaystyle\log f_{ZITD}(y>0|\pi,\mu,\phi,\rho)$
$\displaystyle=\log(1-\pi)+\log f_{\text{TD}}(y>0|\mu,\phi,\rho)$ (19)
$\displaystyle=\log(1-\pi)+\frac{1}{\phi}\bigl{(}y\frac{\mu^{1-\rho}}{1-\rho}-\frac{\mu^{2-\rho}}{2-\rho}\bigl{)}+\log
a(y>0,\phi,\rho)$
$\displaystyle=\log(1-\pi)+\frac{1}{\phi}\bigl{(}y\frac{\mu^{1-\rho}}{1-\rho}-\frac{\mu^{2-\rho}}{2-\rho}\bigl{)}-\log
y+\log\sum_{j=1}^{\infty}\frac{y^{-j\alpha}(\rho-1)^{\alpha
j}}{\phi^{j(1-\alpha)}(2-\rho)^{j}j!\Gamma(-j\alpha)}$
$\displaystyle\geq\log(1-\pi)+\frac{1}{\phi}\bigl{(}y\frac{\mu^{1-\rho}}{1-\rho}-\frac{\mu^{2-\rho}}{2-\rho}\bigl{)}-\log(j_{max}\sqrt{-\alpha}y)+j_{max}(\alpha-1).$

In this equation, $j_{max}=\frac{y^{2-\rho}}{(2-\rho)\phi}$ and
$\alpha=\frac{2-\rho}{1-\rho}<0$. We optimize the lower bound, hence
$NLL_{y>0}$ can be optimized in this manner.

For $y=0$, $NLL_{y=0}=-\log f_{ZITD}(y=0|\pi,\mu,\phi,\rho)$:

$\displaystyle\log f_{ZITD}(y=0|\pi,\mu,\phi,\rho)$
$\displaystyle=\log(\pi)+\log(1-\pi)f_{\text{TD}}(y=0|\mu,\phi,\rho)$ (20)
$\displaystyle=\log(\pi)+\log(1-\pi)+\frac{1}{\phi}\bigl{(}-\frac{\mu^{2-\rho}}{2-\rho}\bigl{)}$

Here, $\pi,\mu,\phi,\rho$ are chosen and calculated in accordance with the
index of $y=0$ or $y>0$. The ultimate negative log-likelihood loss function is
given by:

$NLL_{STZITD}=NLL_{y=0}+NLL_{y>0}+\eta{\Theta^{2}}.$ (21)

Model parameter $\Theta$ can be optimized by minimizing the negative log-
likelihood loss: $\hat{\Theta}=\text{argmin}NLL_{STZITD}$. $\eta$ stands for
the weight parameter for $L_{2}$ Normalization.

It’s worth noting that our model can be extended to other distributions by
modifying the probability layer to fit other distributions. The main reason is
that we learn the parameters that determine distributions. For instance, if we
choose the Gaussian distribution, we can parameterize the probability layer
using the spatiotemporal embedding of the mean and variance, thereby
quantifying the data uncertainty that follows the Gaussian distribution. By
leveraging the flexibility of the probability layer, we design a range of
benchmark models for comparison to the STZITD-GNNs model. Our code is
available on Github 111https://github.com/STTDAnonymous/STTD.

####  3.4.2 Time analysis & Training algorithms

We summarize the training procedure of STZITD-GNNs in Algorithm 1. Based on
that, the time complexity can be analyzed as follows (for the sake of
simplicity, we denote the hidden union of embedding $F^{\prime}\approx F$ and
$F$ as $F$):

  * •

For temporal encoders, the time complexity can be formulated as
$O\bigl{(}tN(dF+F^{2}+F)\bigl{)}$.

  * •

For spatial encoders, we assume the number of road edges to be $|E|$. Hence,
when implementing the graph attention encoder, the complexity of computing
attention becomes $O(M|E|F)$. Furthermore, when executing $W_{a}$
multiplication, the cost is $O(MNF^{2})$. This leads to a total time
complexity of $O(M|E|F+MNF^{2})$.

  * •

When generating the four parameters, the four one-layer MLPs require a time
complexity of $O\bigl{(}F^{2}p\bigl{)}$.

Algorithm 1 The training of STZITD-GNNs.

1:Spatiotemporal road feature $X\in\mathbb{R}^{N\times T\times d}$, road
network graph $G=(V,E,A)$, roads traffic risk $Y\in\mathbb{R}^{N\times T}$,
model hyper-parameters $M,F^{\prime},F$, and total time slot $T$, time window
$t$, $p$.

2:Calculate road traffic risk $Y$ via Eq. (1).

3:Initialize parameters $\Theta$ via Xavier Initializer.

4:while STZITD-GNNs is not converged do $\triangleright$ Train

5: Calculate temporal embedding $\mathcal{Z}_{T}$ via Eq. (14) by road feature
$X_{1:t}$; $\triangleright$ Temporal Encode

6: Calculate spatiotemporal embedding $\mathcal{Z}$ via Eq. (17) by graph $G$
and $\mathcal{Z}_{T}$; $\triangleright$ Spatial Encode

7: Obtain four parameters
$\pi_{t+1:t+p},\mu_{t+1:t+p},\phi_{t+1:t+p},\rho_{t+1:t+p}$ via output layer
refer to Eq. (18);

8: Calculate distribution
$f_{ZITD}(y_{t+1:t+p}|\pi_{t+1:t+p},\mu_{t+1:t+p},\phi_{t+1:t+p},\rho_{t+1:t+p})$
from four parameters via Eq. (10);

9: Minimizing negative log-likelihood loss ${NLL}_{STZITD}$ against
$Y_{t+1:t+p}$ via Eq. (21) using Adam Optimizer;

10:end while

11:End optimizing model parameters $\Theta$;

##  4 Experiments

In this section, we present a thorough evaluation of the proposed STZITD-GNNs
model, employing real traffic incident data from three of London’s most
perilous boroughs in 2019: Lambeth, Tower Hamlets, and Westminster. Our
evaluation is designed to examine the model’s generalizability using two
distinct types of input data: 14 days of historical data for 14-day (Long-
term) predictions, and 7 days of historical data for 7-day (Short-term)
predictions. Our evaluations focus on two aspects: firstly, to illustrate the
ability of our proposed model to effectively quantify the uncertainty of road-
level traffic incident data, and secondly, to showcase its accuracy in
comparison to several relevant baseline models.

We initiate our discussion by presenting our datasets (Section 4.1), followed
by an outline of the experimental setup (Section 4.2). Next, we introduce the
evaluation metrics (Section 4.3) and baseline models (Section 4.4) utilized in
our study. Subsequently, we delve into a comparison of our model’s performance
against alternative methods (Section 4.5), rounded off with a discussion on
supplementary experiments as sensitive analysis (Section 4.6).

###  4.1 Data Description

In this study, we utilize the real-world UK road incidents and safety
statistics, known as STATS19, to validate our model. The STATS19 dataset,
maintained by the Department for Transport in the UK, is comprised of three
interlinked tables providing information about incidents, their severity, and
related vehicle data 222https://www.gov.uk/guidance/road-incident-and-safety-
statistics-guidance.

Figure 2: The map of three boroughs, London

The case study of our research is in three boroughs within Greater London, UK:
Westminster, Lambeth, and Tower Hamlets, as shown in Figure 2. These boroughs
were selected due to their notable high frequency of incidents from January to
December 2019. The detailed statistics for these areas are summarized in Table
3. Upon examination of these three boroughs—which boast the highest number of
incidents on an average of 3 to 5 incidents per day—it’s evident that the zero
rates remain considerably high, even after factoring in the spatial spillover
effects. This pattern arises due to the substantial number of road nodes
present in these boroughs. Despite the relatively frequent occurrence of
incidents, a significant proportion of roads experience no incidents,
resulting in a high zero rate. Furthermore, an examination of the risk value
frequency distribution presented in Figure 3 reveals not only the presence of
a zero-inflated issue but also a clear long-tail effect. Remember that The
long-tail effect refers to a particular type of distribution of variables
where a larger share of occurrences are far away from the ”head” or central
part of the distribution and are spread out, forming a ”long tail”. The data
distributional pattern from Figure 3 indicates that while the majority of risk
values are quite low (equating to 1), the frequency of occurrence decreases as
the risk level increases. Consequently, these conditions make the dataset
particularly suited for the application of our proposed model.

To assess the generalization ability of our model, it is imperative to
consider the unique geospatial aspects of each borough within Greater London.
These 32 local authority districts, along with the City of London, constitute
the administrative area of Greater London, England. Each borough is governed
by a separate London borough council and exhibits distinct socio-demographic
characteristics, road network formats, and transport risk patterns, adding
layers of complexity to the urban risk prediction task (Dixit and Sivakumar,,
2020; Yuan et al.,, 2018). By including boroughs of differing characteristics
in our validation, we strive to demonstrate the broad applicability of the
STZITD-GNNs model.

Table 3: The statistics of Westminster, Lambeth, and Tower Hamlets datasets.

Dataset |  Roads/Nodes |  Edges |  Features |  incident Counts |  Zero-Inflated
---|---|---|---|---|---
Westminster${}^{\text{1}}$ | 4,822 | 20,128 | 110 | 1,745 | 95.72%
Lambeth${}^{\text{1}}$ | 5,659 | 21,574 | 110 | 1,335 | 96.71%
Tower Hamlets${}^{\text{1}}$ | 4,688 | 17,928 | 110 | 1,236 | 96.28%

  * •

${}^{\text{1}}$ https://roadtraffic.dft.gov.uk/downloads

Figure 3: The distribution of Risk values frequency in the three boroughs,
London. Both zero-inflated issue and long-tail effect are significant in the
data distributions.

###  4.2 Experiment Setup

The parameters of STZITD-GNNs are optimized with Adam optimizer Kingma and Ba,
(2014), where the learning rate is set to 0.01, and weight decay is 0.01. In
STZITD-GNNs, $N_{epoch}=20$, the hidden unit is set to 42, GAT head is set to
3. The parameters of baselines are optimized using the Adam Kingma and Ba,
(2014) with $L_{2}$ regularization and a dropout rate of 0.2. The GNNs in
STZITD-GNNs and baselines are all two-layered. We set the hidden units to 42
as well. Similar to previous work, we also employ the early-stopping strategy
with patience equals 10 to avoid over-fitting.

STZITD-GNNs is implemented in Pytorch 1.9.0 with Python 3.8. All the
experiments are conducted on 1 NVIDIA GeForce RTX 3090, 24 GB. All baselines
are implemented by ourselves.

###  4.3 Evaluation Metrics

We use eight metrics to evaluate the performance of our STZITD-GNNs from four
aspects, including the predictive regression effect of road incident risk, the
predictive effect of distributional characteristics, the degree of zero
inflation fit, and the hit rate of high-risk roads.

####  4.3.1 Accuracy Metrics

For the point estimate aspect, we compare the numerical value prediction. For
probabilistic distribution outputs like our model, we use the mean value of
the distribution as the numerical output. For specific metrics, we use Mean
Absolute Error (MAE), Mean Absolute Percentage error (MAPE), and Root Mean
Squared Error (RMSE) to evaluate the predicted risks of all roads. Specific
definitions are summarized:

$\displaystyle\text{MAE}=\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}\left|y_{ij}-\hat{y}_{ij}\right|,$
(22)
$\displaystyle\text{MAPE}=\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}\left|\frac{y_{ij}-\hat{y}_{ij}}{y_{ij}}\right|,$
$\displaystyle\text{RMSE}=\sqrt{\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}(y_{ij}-\hat{y}_{ij})^{2}},$

where $\hat{y}_{ij}$ and $y_{ij}$ are the predicted (i.e. distribution mean)
and ground-truth values of $i^{th}$ road at time slot $j^{th}$ respectively.
Literally, $p$ = 14 or 7 days, and time slot $j^{th}$ is inside this predicted
window. Lower values of these metrics mean better performance.

####  4.3.2 Uncertainty Quantification Metric

To assess the measured uncertainty, we utilize a selection of metrics:
Kullback-Leibler Divergence (KLD), Mean Prediction Interval Width (MPIW), and
Prediction Interval Coverage Probability (PICP) predicated on a 5%-95%
confidence interval with a confidence level set at 5% ($\alpha=5\%$). The KLD
is deployed to gauge the degree of similarity between the output distribution
generated by our model and the distribution of the test set. Concurrently, the
MPIW serves as a tool to determine the mean width of the prediction interval,
with its output heavily contingent upon the distinct traits of the
distribution. It’s optimal for the values of MPIW and KLD to be as small as
possible. The PICP, on the other hand, is implemented to evaluate whether the
ground-truth value, denoted as $y_{ij}$, lies within the predicted interval. A
higher PICP is regarded as more advantageous. It is important to recognize the
inherent trade-off between PICP and MPIW: an ideal prediction interval would
be both narrow and cover a substantial portion of the data points. However, a
wider MPIW naturally leads to a larger PICP. This interaction poses a
considerable challenge for any model that is designed to cater to more precise
assumptions, as elucidated in Wang et al., 2023a  and Zhu et al., (2022).

$\displaystyle\text{KLD}=\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}(\hat{y}_{ij}\log\frac{\hat{y}_{ij}+\epsilon}{{y}_{ij}+\epsilon}),$
(23)
$\displaystyle\text{MPIW}=\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}(U_{ij}-L_{ij}),$
$\displaystyle\text{PICP}=\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}\mathcal{I}(U_{ij}<\hat{y}_{ij}<L_{ij}),$

since many $y_{ij}$ and $\hat{y}_{ij}$ values are likely to be zeros, a
minimum value $\epsilon=10^{-5}$ is used to avoid numerical issues because of
division by 0. Since the KLD measures the difference between two
distributions, smaller values are desirable. Moreover, $L_{ij}$ and $U_{ij}$
correspond to the lower and upper bound of the confidence interval for
observation road $i$ at time slot $j^{th}$. Indicate function
$\mathcal{I}(U_{ij}<\hat{y}_{ij}<L_{ij})$ equals 1 if the condition is true,
else 0.

####  4.3.3 Zero Inflated Metrics

Given the notable data sparsity issue, as evident in Figure 3, it becomes
crucial to assess the proficiency of models in capturing these zero entries.
To evaluate the model’s performance on the discrete road risk entries, we
leverage the true-zero rate (ZR) measurements. The ZR acts as a quantifier of
how well the model mirrors the sparsity observed in the ground-truth data. A
ZR score value closer to the ground-truth zero-inflated ratios signifies
superior model performance. The mathematical expression for the ZR is given
as:

$\displaystyle\text{ZR}={\frac{1}{pN}\sum_{j=1}^{p}\sum_{i=1}^{N}\mathcal{I}(y_{ij}=0\cap\hat{y}_{ij}<\epsilon)},$
(24)

where the indicator function, $\mathcal{I}(y{ij}=0\cap\hat{y}{ij}<\epsilon)$,
equals 1 if the condition holds true, and 0 otherwise. We set a limit
approaching zero ($\epsilon\to 0$) as the threshold and consider
$\hat{y}{ij}<\epsilon$ as equivalent to zero.

####  4.3.4 Hit Rate Metric

In order to evaluate the proficiency of our model in identifying roads with
increased incident risk, we utilize the Average Hit Rate (HR) metric. This
metric assesses the model’s performance by determining the fraction of risk
scenarios accurately flagged within the top $a$ percent of locations predicted
to bear the highest incident probability. In particular, we scrutinize the top
20% ($a=20\%$) of high-risk road length coverage and ascertain whether these
roads indeed exhibit a high risk with actual incidents transpiring at these
specific locales and times (Zhang and Cheng,, 2020; Zhang et al., 2019b, ). By
correlating high-risk forecasts with actual incidents, $HR(a)$ serves as a
robust gauge of our model’s predictive precision for areas demonstrating a
heightened incident risk. The mathematical expression for HR is given as:

$\displaystyle
HR(a)=\frac{1}{p}\sum_{j=1}^{p}\frac{\sum_{i=1}^{N}\mathcal{I}(y_{ij}>0\cap\hat{y}_{ij}\geq\text{Top}(\hat{y}_{:j},a))}{\sum_{i=1}^{N}\mathcal{I}(y_{ij}>0)},$
(25)

The indicator function,
$\mathcal{I}(y_{ij}>0\cap\hat{y}_{ij}\geq\text{Top}(\hat{y}{:j},a))$, equals 1
if the condition is satisfied, and 0 otherwise. $\mathcal{I}(y_{ij}>0)=1$ if
an incident occurs, and it equals 0 otherwise. In this context, $y_{ij}>0$
denotes that an incident indeed transpired on the $i^{th}$ road during the
$j^{th}$ time slot, indicated by a non-zero risk value. Conversely,
$\text{Top}(\hat{y}_{:j},a)$ represents the boundary value corresponding to
the top $a$ percentage of predicted risk levels for all road length coverage
at the $j^{th}$ time slot.

###  4.4 Baseline Models

The following baselines are considered in our comparative evaluation of our
proposed model, including the mean-predict method (HA), spatial-temporal graph
deep learning methods such as STGCN and STGAT, parameterization models such as
STG, STNB, STTD and STZINB:

  * •

HA(Liu and Guan, (2004)): Historical Average serves as the statistic baseline.
It is calculated by averaging the road traffic risks in the same time
intervals from the historical data to predict one step ahead of future road
traffic risk.

  * •

STGCN(Yu et al., (2018)): Spatiotemporal Graph Convolution Network, which
utilizes graph convolution and 1D convolution to capture spatial dependencies
and temporal correlations, respectively.

  * •

STGAT(Veličković et al., (2017)): Graph Attention Network is able to attend
over their neighbourhoods’ features, specifying different weights to different
nodes in a neighbourhood. We replace GCN with GAT in STGCN and call this
STGAT.

  * •

STG-GNNs(Wang et al., 2023a ): Same as STZITD-GNNs, we replace the statistical
assumption of road traffic risks with Gaussian Distribution instead.

  * •

STNB-GNNs(Zhuang et al., (2022)): Same as STZITD-GNNs, we replace the
statistical assumption of road traffic risks with Negative Binomial
Distribution instead.

  * •

STTD-GNNs: Same as STZITD-GNNs, we replace the statistical assumption of road
traffic risks with TD Distribution instead.

  * •

STZINB-GNNs(Zhuang et al., (2022)): Same as STZITD, we replace the statistical
assumption of road traffic risks with Zero-Inflated Negative Binomial
Distribution instead.

###  4.5 Performance Comparison

A series of experiments were executed to accomplish the task of quantifying
uncertainty in road-level traffic incident risk across three diverse datasets
under two scenarios: long-term and short-term. The experimental results
encompassing the performance evaluation of our proposed model, STZITD-GNNs,
and the existing baselines are displayed in Tables 4 and 5. The tabulated
results use a colour-coded scheme to aid interpretation, with the red and blue
zones highlighting the best and second-best performance scores, respectively.
The direction of the vectors serves to visually represent the comparative
performance of the STZITD-GNNs, depicting whether it surpasses (upward vector)
or lags behind (downward vector) the best scores of the baseline models.

In terms of predictive accuracy, it is evident that our STZITD-GNNs model
consistently outperforms the baselines in the majority of scenarios. This
superiority is especially pronounced in the context of long-term prediction.
When juxtaposed with the second-best baseline model, the accuracy improvements
offered by our model span from $11.15\%$ to $38.34\%$ in terms of MAE,
$10.16\%$ to $49.06\%$ for MAPE, and a substantial $9.08\%$ to $34.60\%$ for
RMSE. This applies to both long-term and short-term prediction in node
estimation tasks. It is worth noting that STZITD-GNNs shares similar
spatiotemporal graph neural network structures as STG-GNNs, STNB-GNNs, STZINB-
GNNs, and STTD-GNNs. The consistently superior results demonstrate the
efficacy of adopting the ZITD distribution, which, due to its distinct
characteristics, is able to capture the nuances of the ground truth data more
comprehensively.

In terms of the uncertainty quantification task, the STZITD-GNNs model
outstrips all baseline models for long-term prediction. This provides robust
evidence that our model exhibits superior capacity in capturing long-term
spatiotemporal dependencies compared to short-term prediction. The STTD-GNNs
model closely follows STZITD-GNNs in the PICP metric for the Lambeth and
Westminster datasets in short-term prediction scenarios. Moreover, the STZITD-
GNNs model delivers a substantial performance uplift in the KLD metric.
Particularly in zero-inflated situations, the STZITD-GNNs model excels at
modelling zero-risk roads, presenting improvements ranging from $12.60\%$ to
$24.11\%$ across all tested scenarios. STZITD-GNNs also outperform the
baselines in HR (20%), underscoring its capability to accurately identify
high-incident risk roads and offer valuable preventive recommendations.
Consequently, integrating the sparsity parameter $\pi$ bolsters the model’s
ability to handle the excessive zeroes in the data, enabling it to form more
accurate prediction confidence intervals.

An intriguing observation from the experimental results is the relatively
comparable performance of STTD-GNNs and STZINB-GNNs to the STZITD-GNNs model.
These models incorporate both temporal dependencies and spatial relations and
account for the zero-inflated parameters. While STZINB-GNNs showcase superior
performance in uncertainty quantification owing to its explicit modelling of
the zero-inflated parameter, STTD manages to outperform STZINB in other
metrics. This can be attributed to STTD-GNNs’s unique ability to address both
the long-tail problem as well as the zero-inflated issue. However, it tends to
struggle with extremely zero-inflated data. Contrarily, the HA model, which
neglects spatial dependencies, yields the poorest performance when compared to
STGCN and STGAT. This emphasizes the indispensability of considering spatial
and temporal dependencies for spatial-temporal prediction tasks.

We must acknowledge the well-documented trade-off between PICP and MPIW, as
discussed in Section 4.3.2. However, our STZITD-GNNs model, by integrating the
zero-inflation parameter $\pi$ into the traditional TD, and transforming it
into a ZITD, remarkably manages to alleviate this trade-off. The zero-
inflation component of the ZITD plays a pivotal role in this context as it
effectively captures the substantial proportion of zero-valued (non-risk)
roads in the dataset. This enhances the model’s ability to represent the true
distribution of road risks more accurately, which, in turn, facilitates the
formation of narrower prediction intervals for a majority of cases, thus
lowering the MPIW. Concurrently, the model retains its capacity to encompass a
substantial proportion of the data points, thereby maintaining a high PICP.
This dual ability of the STZITD-GNNs model significantly mitigates the
traditionally challenging trade-off between MPIW and PICP, which considerably
boosts the model’s overall predictive performance. Additionally, it enhances
the model’s capacity to effectively capture and represent the complexity and
characteristics of the urban road-level risk data.

In summary, the four-parameter STZITD-GNNs model emerges as the best-
performing model, outperforming the three-parameter models such as STTD-GNNs
and STZINB-GNNs, as well as the two-parameter models and single-parameter
methods like STGCN and STGAT. These plots provide insights into how the
STZITD-GNNs model effectively learns and adjusts the TD distribution
parameters to cater to the characteristics of urban traffic risk data.
Experiment results underscore the robustness of our proposed STZITD-GNNs model
and its ability to effectively quantify uncertainty in road-level traffic
incident risk.

Table 4: Comparison of Long-Term (14-day) Prediction Performance Across Three Boroughs. Red shading indicates the best-performing model, while blue signifies the second-best. To maintain clarity, detailed names of all GNNs in the final five statistical-integrated deep learning models are omitted. Instead, abbreviations such as STG, STNB, STTD, STZINB, and STZITD are used to underscore the most significant modifications within the GNN models. Metric | Model | HA | STGCN | STGAT | STG | STNB | STTD | STZINB | STZITD
---|---|---|---|---|---|---|---|---|---
Dataset | Lambeth
ACC | MAE | 0.1348 | 0.0613 | 0.0302 | 0.1182 | 0.0802 | 0.0315 | 0.0544 |  ${0.0238}_{\uparrow 21.19\%}$
MAPE | 0.4142 | 0.4023 | 0.3203 | 0.4046 | 0.0363 | 0.0426 | 0.0265 |  ${0.0135}_{\uparrow 49.06\%}$
RMSE | 0.2107 | 0.1436 | 0.1122 | 0.1827 | 0.1390 | 0.1098 | 0.1191 |  $0.0947_{\uparrow 13.75\%}$
Uncertainty | KLD | 0.2691 | 0.1068 | 0.0512 | 0.5039 | 1.5580 | 0.0849 | 0.2596 |  $0.0248_{\uparrow 51.56\%}$
MPIW(a=5%) | / | / | / | 0.4871 | 0.1946 | 0.1349 | 0.0454 |  $0.0204_{\uparrow 55.07\%}$
PICP(a=5%) | / | / | / | 0.7328 | 0.8903 | 0.9873 | 0.9796 |  $0.9899_{\uparrow 00.26\%}$
Zero Inflated | ZR | 0.5203 | 0.5305 | 0.5986 | 0.1992 | 0.5621 | 0.6838 | 0.6408 |  $0.7870_{\uparrow 15.09\%}$
Hit Rate | HR(20%) | 0.4520 | 0.6113 | 0.6422 | 0.2666 | 0.4471 | 0.7123 | 0.6184 |  $0.7659_{\uparrow 07.39\%}$
Dataset | Tower Hamlets
ACC | MAE | 0.0939 | 0.2385 | 0.0305 | 0.1741 | 0.2475 | 0.0433 | 0.2135 |  $0.0271_{\uparrow 11.15\%}$
MAPE | 0.2896 | 0.2017 | 0.1205 | 0.2662 | 0.0399 | 0.0246 | 0.0594 |  $0.0221_{\uparrow 10.16\%}$
RMSE | 0.1987 | 0.1325 | 0.1355 | 0.3896 | 0.1548 | 0.1186 | 0.1368 |  $0.0918_{\uparrow 22.60\%}$
Uncertainty | KLD | 0.2888 | 0.0612 | 0.0319 | 0.8090 | 1.9791 | 0.1096 | 0.6278 |  $0.0132_{\uparrow 58.62\%}$
MPIW(a=5%) | / | / | / | 0.5321 | 0.2263 | 0.0810 | 0.0296 |  $0.0156_{\uparrow 47.30\%}$
PICP(a=5%) | / | / | / | 0.6439 | 0.9442 | 0.9845 | 0.9845 |  $0.9871_{\uparrow 00.26\%}$
Zero Inflated | ZR | 0.4998 | 0.4896 | 0.5234 | 0.1351 | 0.5158 | 0.6420 | 0.7526 |  $0.8876_{\uparrow 17.94\%}$
Hit Rate | HR(20%) | $0.4752$ | 0.5869 | 0.6950 | 0.2647 | 0.5022 | 0.6368 | 0.5827 |  $0.7224_{\uparrow 03.94\%}$
Dataset | Westminster
ACC | MAE | 0.1467 | 0.0784 | 0.0847 | 0.1055 | 0.3266 | 0.0579 | 0.3240 |  $0.0357_{\uparrow 38.34\%}$
MAPE | 0.1884 | 0.8976 | 0.6012 | 2.9295 | 0.0330 | 0.0973 | 0.0415 |  $0.0259_{\uparrow 21.52\%}$
RMSE | 0.5999 | 0.2384 | 0.2433 | 0.6681 | 0.2909 | 0.1552 | 0.2751 |  $0.1015_{\uparrow 34.60\%}$
Uncertainty | KLD | 0.4149 | 0.2218 | 0.2626 | 2.6582 | 2.3624 | 0.1452 | 0.3889 |  $0.0331_{\uparrow 77.20\%}$
MPIW(a=5%) | / | / | / | 0.8542 | 0.6272 | 0.2021 | 0.0493 |  $0.0259_{\uparrow 47.46\%}$
PICP(a=5%) | / | / | / | 0.6036 | 0.9363 | 0.9801 | 0.9796 |  $0.9893_{\uparrow 00.94\%}$
Zero Inflated | ZR | 0.5981 | 0.1431 | 0.2922 | 0.2071 | 0.5393 | 0.5079 | 0.6288 |  $0.7328_{\uparrow 16.54\%}$
Hit Rate | HR(20%) | 0.4217 | 0.5020 | 0.4808 | 0.2933 | 0.4503 | 0.6075 | 0.5139 |  $0.6898_{\uparrow 13.55\%}$
Table 5: Comparison of Short-Term (7-day) Prediction Performance Across Studied Boroughs. Red shading indicates the best-performing model, while blue signifies the second-best. To maintain clarity, detailed names of all GNNs in the final five statistical-integrated deep learning models are omitted. Instead, abbreviations such as STG, STNB, STTD, STZINB, and STZITD are used to underscore the most significant modifications within the GNN models.  Metric | Model | HA | STGCN | STGAT | STG | STNB | STTD | STZINB | STZITD
---|---|---|---|---|---|---|---|---|---
Dataset | Lambeth
ACC | MAE | 0.1339 | 0.0815 | 0.0316 | 0.0482 | 0.1051 | 0.0378 | 0.0774 |  $0.0249_{\uparrow 21.20\%}$
MAPE | 0.4848 | 0.3553 | 0.2025 | 0.4427 | 0.0781 | 0.0495 | 0.0245 |  $0.0213_{\uparrow 13.06\%}$
RMSE | 0.2377 | 0.1712 | 0.1084 | 0.1847 | 0.1715 | 0.1151 | 0.1398 |  ${0.0978}_{\uparrow 09.78\%}$
Uncertainty | KLD | 0.2641 | 0.1561 | 0.0628 | 0.5218 | 0.7585 | 0.1031 | 0.4734 |  ${0.0291}_{\uparrow 54.66\%}$
MPIW(a=5%) | / | / | / | 0.4920 | 0.2296 | 0.1170 | 0.0070 |  ${0.0098}_{\downarrow 40.00\%}$
PICP(a=5%) | / | / | / | 0.6944 | 0.8897 | 0.9807 | 0.9874 |  ${0.9908}_{\uparrow 00.34\%}$
Zero Inflated | ZR | 0.5031 | 0.3610 | 0.6359 | 0.1023 | 0.5177 | 0.6506 | 0.5797 |  ${0.7611}_{\uparrow 16.98\%}$
Hit Rate | HR(20%) | 0.4432 | 0.5462 | 0.6592 | 0.3006 | 0.4417 | 0.6867 | 0.5749 |  ${0.7364}_{\uparrow 07.24\%}$
Dataset | Tower Hamlets
ACC | MAE | 0.1957 | 0.0375 | 0.0398 | 0.0997 | 0.1159 | 0.0547 | 0.0827 |  ${0.0311}_{\uparrow 17.07\%}$
MAPE | 0.4997 | 0.4860 | 0.2405 | 0.4542 | 0.0509 | 0.0574 | 0.0562 |  ${0.0423}_{\uparrow 16.90\%}$
RMSE | 0.2258 | 0.1273 | 0.1319 | 0.2811 | 0.1013 | 0.1261 | 0.1219 |  ${0.0921}_{\uparrow 09.08\%}$
Uncertainty | KLD | 0.2680 | 0.0596 | 0.0704 | 0.3696 | 0.7678 | 0.1322 | 0.4842 |  ${0.0155}_{\uparrow 73.99\%}$
MPIW(a=5%) | / | / | / | 0.5688 | 0.1223 | 0.0918 | 0.0052 |  ${0.0043}_{\uparrow 17.31\%}$
PICP(a=5%) | / | / | / | 0.6218 | 0.9250 | 0.9450 | 0.9748 |  ${0.9887}_{\uparrow 01.43\%}$
Zero Inflated | ZR | 0.4943 | 0.5417 | 0.5132 | 0.0890 | 0.5219 | 0.6052 | 0.6886 |  ${0.8546}_{\uparrow 24.11\%}$
Hit Rate | HR(20%) | 0.4693 | 0.6123 | 0.5661 | 0.2919 | 0.4882 | 0.6150 | 0.6021 |  ${0.6885}_{\uparrow 11.95\%}$
Dataset | Westminster
ACC | MAE | 0.1448 | 0.0676 | 0.0673 | 0.4521 | 0.1874 | 0.0594 | 0.1295 |  ${0.0427}_{\uparrow 28.11\%}$
MAPE | 0.4473 | 1.0104 | 0.7151 | 2.3351 | 0.0388 | 0.1183 | 0.0569 |  ${0.0310}_{\uparrow 20.10\%}$
RMSE | 0.1504 | 0.2044 | 0.1985 | 0.6541 | 0.2175 | 0.1548 | 0.1573 |  ${0.1357}_{\uparrow 12.34\%}$
Uncertainty | KLD | 0.4485 | 0.1789 | 0.1795 | 2.8741 | 0.9136 | 0.1595 | 0.7628 |  ${0.0628}_{\uparrow 60.63\%}$
MPIW(a=5%) | / | / | / | 0.8789 | 0.3813 | 0.1594 | 0.0032 |  ${0.0095}_{\downarrow 196.9\%}$
PICP(a=5%) | / | / | / | 0.5766 | 0.9151 | 0.9804 | 0.9762 |  ${0.9857}_{\uparrow 00.54\%}$
Zero Inflated | ZR | 0.5340 | 0.2030 | 0.2439 | 0.2032 | 0.5091 | 0.5258 | 0.6112 |  ${0.6882}_{\uparrow 12.60\%}$
Hit Rate | HR(20%) | 0.3827 | 0.5595 | 0.5274 | 0.2852 | 0.4456 | 0.6332 | 0.4859 |  ${0.6167}_{\downarrow 02.60\%}$

####  4.5.1 Visualization of Predicted Road-Level Risk

To illustrate the practical implications of our STZITD-GNNs model, we have
generated choropleth maps that visually represent road-level risks across
three distinct boroughs over two separate time windows. These maps are
presented in Figure 4, which depicts risks over a long-term period of 14 days,
and Figure 5, which captures risks over a short-term period of 7 days. Each
borough encompasses unique risk patterns, which are juxtaposed against our
model’s predicted traffic risk values to provide a comprehensive perspective
of its predictive capacity.

Upon examination of these visualizations, our model’s high predictive accuracy
is notably consistent across both the 7-day and 14-day time windows. It
accurately identifies roads with a range of risk levels, as well as those with
negligible risk, across all three boroughs. A particularly striking instance
of this precision is the model’s predictions for the Westminster borough, an
area characterized by a lower zero-inflation rate and a more dense graph
structure.

Long-term risk representations, as seen in Figure 4, offer a robust overview
of the risk landscape. They enable the identification of a larger set of risk-
prone roads with remarkable accuracy. This reiterates the model’s capacity to
leverage long-term data effectively, enhancing its predictive strength.
However, it is also evident that there are opportunities for improvement,
particularly in predicting zero-risk roads. In the cases of Lambeth and Tower
Hamlets, the model occasionally identifies zero-risk roads as low-risk. While
these errors predominantly fall within the lower-risk category and do not
greatly impact the identification of the highest-risk roads, addressing this
can help further refine the model’s performance.

Figure 5 showcases the model’s capacity to effectively highlight high-risk and
zero-risk roads, despite the reduced 7-day input data window. The importance
of this capacity cannot be understated, as it offers invaluable insights for
transport officers and road users alike. These insights can be utilized to
plan and prioritize mobility options to prevent unanticipated losses, even
when only a week’s worth of data is available. Moreover, our model
consistently and accurately identifies the spatial compactness (or
connectedness) of high-risk roads across all three boroughs, irrespective of
the prediction window. Although this was not a primary focus of our model, the
detection of compactness can prove highly beneficial in practical contexts.
Compact and connected high-risk roads can be patrolled more efficiently with
limited resources, a fact that is echoed in previous studies Adepeju et al.,
(2016); Zhang and Cheng, (2020).

In light of the impressive PICP score exceeding 98%, as per Table 4 and Table
5, we confidently assert that the STZITD-GNNs model excels in predicting road-
level risks through learning from historical spatiotemporal patterns. This
proficiency persists despite the varying spatial patterns and daily road risk
counts across the three boroughs, highlighting the model’s adaptability and
robust performance across diverse contexts. This consistent accuracy
underscores the model’s potential as a powerful instrument in urban risk
management, capable of providing critical insights to guide both policy making
and daily navigation decisions.

Figure 4: Predicted long-term (14-day) risk on Lambeth, Tower Hamlets, and
Westminster. Figure 5: Predicted short-term (7-day) risk on Lambeth, Tower
Hamlets, and Westminster.

####  4.5.2 Hit Rate study under different coverage of the roads

Figure 6 presents a mean hit rate line chart, spanning the three boroughs
under consideration, averaged over the entire testing period. In this
depiction, the coverage of the road length gradually increases, extending up
to 20%. On an overarching level, the STZITD-GNNs model maintains a high mean
hit rate across varying incident scenarios, with a minor exception noted in
the short-term predictive performance for Westminster. More specifically, the
accuracy exceeds 80% until the coverage broadens to 20%, thereby underscoring
the model’s consistent and dependable performance.

When drawing comparisons between the three boroughs, the parallels between
short and long-term predictive hit rates are most prominent within the Tower
Hamlets dataset. This resemblance could potentially be ascribed to the less
sparse distribution of traffic incidents within this region over the testing
period, thereby resulting in less pronounced variability in the predictions.
Alternatively, the winter season, which overlaps with our testing period,
could potentially bring about more sparse traffic incident scenarios in
Westminster. Consequently, due to the scarcity of input information, the
short-term predictions within this borough might be marginally affected, which
is reflected in a slightly lower hit rate.

It is worth noting that the ability of the STZITD-GNNs model to maintain high
predictive accuracy, even with the increasing coverage of the road length,
indicates its capacity to efficiently handle larger spatial scales. This
characteristic is particularly beneficial in practical applications, where the
model may need to predict risks for extensive road networks. Moreover, the
model’s resilience against changes in the spatiotemporal distribution of
incidents demonstrates its adaptability to varying urban dynamics, further
enhancing its applicability in real-world scenarios.

Figure 6: Hit rate increase figure under different road network coverage rates
for long and short prediction intervals.

####  4.5.3 Comparative Analysis of UQ for Baseline Models and STZITD-GNNs
Through a Time-wise MPIW Study

To evaluate the UQ performance of our STZITD-GNNs model in contrast to the
baseline models, we conduct a time-wise study of the MPIW. This assessment is
depicted in Figure 7 and Figure 8 and extends to each prediction step for both
short and long-term predictions across the three boroughs. The objective is to
establish an MPIW that is narrow yet encompasses a significant proportion of
the test samples, consequently ensuring a robust and dependable prediction
interval for traffic risk.

Figure 7: Time-wise MPIW changing lines for long intervals predictions for
baselines and STZITD-GNNs. Figure 8: Time-wise MPIW changing lines for short
intervals predictions for baselines and STZITD-GNNs.

Our STZITD-GNNs model consistently narrows down the prediction intervals by
dynamically fine-tuning the prediction interval in accordance with the risk
particular to each road segment. This dynamism effectively allows the STZITD-
GNNs model to outperform the baseline models in the majority of tasks. This
success can be attributed to the enhanced proficiency of the ZITD model in
handling overdispersed, zero-inflated urban road-level risk data, surpassing
the capabilities of other distribution models.

The Gaussian and NB distributions have limitations in this context: Gaussian
distributions often overlook the preponderance of zero values, while NB
distributions concentrate only on the binary risk outcomes, thereby neglecting
long-tail risk values. Both these shortcomings lead to a larger MPIW. While
ZINB and ZITD both consider zero-value skewness, the ZITD distribution
utilizes sparsity parameter $\pi$ weighted dispersion parameter $\phi$ to map
the substantial proportion of zero-valued (non-risk) roads and TD parameters
to map low-risk roads as well as long-tail middle and high-risk roads. Thanks
to its ability to capture the distinct characteristics of this complex road-
level risk data distribution, it helps to mitigate overfitting to the data,
leading to a narrower MPIW with a higher level of confidence.

Our analysis shows that the smallest MPIW typically occurs on the first
testing day, followed by a rise over time. In figure 7, this temporal increase
is less erratic in the long-term (14-day input) predictions compared to the
baseline models. However, figure 8 exhibits noticeable fluctuations in both
baseline models and STZITD-GNNs. This trend may reflect the inherent
uncertainty and randomness of the risk situation with only a 7-day input
period, resulting in less stability due to the limited spatiotemporal
information.

Moreover, a fluctuating MPIW appears in the short-term predictions for
Westminster and Lambeth, potentially stemming from a higher prevalence of
zero-risk values. This could introduce greater uncertainties for the
prediction intervals (Zhuang et al.,, 2022), given the infrequent occurrence
of traffic risks within the 7-day testing period, thereby augmenting
prediction uncertainty. Nevertheless, the MPIW in these two regions is
commensurate with that of the state-of-the-art STZINB-GNN model and is notably
lower than other baseline models.

Though a larger MPIW can ensure a good coverage rate for the true risk values,
it may also frequently underestimate or overestimate risk information, even
when the model has relative confidence in the prediction (Wang et al., 2023a,
; Zhou et al.,, 2021). Thus, we recognize the merits of our proposed model in
maintaining stability and reducing variance in predicted risk values at the
road level. Notably, our model’s superior performance remains consistent
across the three boroughs, demonstrating its adaptability to different road
conditions, even under diverse spatiotemporal embeddings.

###  4.6 Sensitivity Analysis

####  4.6.1 Analysis of long-tail risk data characteristics

The long-tail effect is an integral attribute of data distribution, playing a
key role in the selection of statistical components for traffic risk
evaluation. To delve deeper into the influence of this long-tail effect on our
model, we produce a four-dimensional 3D surface plot, as depicted in Figure 9.
However, due to the inherent limitations of visualization techniques,
representing all four parameters in a single five-dimensional plot is not
feasible. Moreover, as the sparsity parameter $\pi$ acts as a weight for the
original TD parameters, we focus on the TD parameters ($\phi$, $\rho$, $\mu$)
in our graphic representation. The points plotted to represent predicted risk
values, while the actual risk values $y$ are expressed through the colour bar
in the plot.

Upon visualizing data across three boroughs under two temporal scenarios, it
becomes clear that the learned $\rho$ values exhibit a consistent trend
towards larger values, often approaching 2. Furthermore, This pattern
underscores our model’s capacity to capture the long-tail characteristic
inherent in the risk distribution across diverse geographical and temporal
contexts. This observed behaviour also stems from the inherent properties of
the TD. Herein, $\rho$ significantly influences the distribution’s dispersion,
with larger $\rho$ values denoting a heavier-tailed distribution, a
characteristic that aligns well with our risk data’s long-tail nature. The
increase in $\rho$ can be considered an adaptive response to our loss
function’s component, $y\frac{\mu^{1-\rho}}{1-\rho}(y>0)$. This component
penalizes underestimations of true risk values situated within the long tail,
thereby steering the model towards more accurate risk predictions within these
zones.

Concurrently, we observe a marked increase in the $\phi$ values linked with
zero-valued data, particularly pronounced in short-term scenarios. This can be
attributed to the notably zero-inflated nature of urban risk data,
characterized by significant segments of the road network without recorded
incidents, a trait that becomes more apparent in short-term periods due to
lower incident occurrence. As stated in the section 3.2.1. in the original TD
distribution, the $\phi$ parameter, associated with the frequency of zero
occurrences, is essential in addressing this zero-inflation. However, with the
exceedingly imbalanced nature of our data, this parameter may not adequately
capture the prevalence of zeros. Hence, in our ZITD model, the amplified
$\phi$, under the influence of the sparsity parameter $\pi$, skillfully
accommodates this feature, leading to a more accurate depiction of the
inherent sparsity and discrete uncertainty in road-level risk prediction.

Further scrutiny reveals that the model parameter $\mu$, which serves as the
mean of the predicted risk values, exhibits a colour distribution that is
well-aligned with the real risk values $y$. This direct correspondence further
underlines the model’s robust performance in mean point estimation, adding
another layer of validation to the effectiveness of the STZITD-GNNs model
within the sphere of urban risk prediction tasks. The discernible correlation
between $\mu$ and $y$ also testifies to the model’s capacity to seamlessly
capture and adapt to the subtleties and nuances of urban risk data. This
observation further consolidates the argument for the adoption of the STZITD-
GNNs model for urban risk prediction and underscores its efficacy in handling
the complexity of the task.

Figure 9: Three Tweedie parameters visualization on Lambeth, Tower Hamlet and
Westminster with long (14 days) and short (7 days) predict. The colour bar
shows the predicted risk values from highest to zero.

####  4.6.2 Analyzing Risk Data Through Lower-Dimensional Representation

The Multidimensional Scaling (MDS) algorithm, which is employed for the
purpose of data dimensionality reduction and visualization, affords us the
ability to transform high-dimensional data into a lower-dimensional space
(Wickelmaier,, 2003). It achieves this while preserving the relative distances
between individual data points. This particular algorithm permits the
extraction of inherent structural information and latent relationships present
within the dataset, thereby enabling an intuitive visual observation and
subsequent analysis of the data. With this approach, we utilize the MDS
algorithm to validate the efficacy of our model in capturing long-tail risk
efficiently, by visualizing the learned node embeddings in a 2D space across
six distinct scenarios.

For the purposes of our analysis, risks greater than 1 are regarded as a long-
tail or high-risk road. We then visualize the learned spatiotemporal embedding
$\mathcal{Z}$ on the test set. A method’s representational ability is judged
on its capacity to learn representative embeddings in which similar nodes map
to proximate objects. As such, nodes that belong to the same class are
expected to cluster together in the MDS representation. The outcomes of these
visualizations are displayed in Figure 10.

Analysis of these six figures indicates that the long-tail class tends to be
naturally clustered and embedded into the low-risk class. This phenomenon can
be attributed to the risk spatial spillover effects, wherein adjacent roads
may share similar embeddings, despite having differing risk levels. While the
clusters of low-risk and high-risk spatiotemporal embeddings are distinctly
separate, high-risk classifications are still discernible, particularly within
the Westminster dataset. Both Lambeth and Tower Hamlets datasets continue to
display spatial aggregations, signifying that our model effectively
incorporates the severity of risk in its embeddings. These observations, taken
together, suggest that our proposed methodology is capable of adequately
capturing the long-tail risk within the data, thereby validating its
effectiveness.

This ability of the STZITD-GNNs model to identify and distinctly represent
various risk categories in the low-dimensional embedding space enhances its
interpretability. It provides an intuitive way to explore and comprehend the
high-dimensional spatiotemporal risk patterns, contributing towards its
potential applicability in real-world urban risk prediction tasks.
Furthermore, the distinct clustering of low-risk and high-risk road segments,
despite the inherent spatial spillover effects, underlines the model’s
resilience in capturing the nuances of the risk landscape.

(a) Lambeth, Long

(b) Tower Hamlets, Long

(c) Westminster, Long

(d) Lambeth, Short

(e) Tower Hamlets, Short

(f) Westminster, Short

Figure 10: Low risk (green color) and high risk (blue color) spatiotemporal
embedding $\mathcal{Z}$ MDS-visualization on Lambeth, Tower Hamlets and
Westminster with long (14 days) and short (7 days) predict.

##  5 Conclusion and Discussion

Road traffic incidents remain a significant obstacle to the sustainability of
urban environments and the dependability of transportation systems (Yang et
al.,, 2021; Zhou et al.,, 2020). These incidents not only endanger the safety
of road users but also introduce complexities in resource allocation and the
efficacy of transportation infrastructures, underscoring the necessity for
marked enhancements in predictive methodologies. The urgency of this
advancement extends beyond academic or economic boundaries; it is intimately
tied to the preservation of human lives (Ibrahim et al.,, 2021). Given the
comprehensive implications of these incidents, the development of models
capable of delivering both precise and reliable predictions for practical
applications has emerged as a pressing priority.

The current landscape of incident risk prediction methodologies, while
providing insightful contributions, also exposes areas for improvement. Low-
resolution models, often preferred for their computational efficiency, tend to
compromise on precision, leading to a high bias that may inhibit effective
risk mitigation efforts. Additionally, frequency-based estimation techniques,
while competent at pinpointing potential incident risks, often fall short in
providing a complete representation of confidence in their predicted risk
values, thus introducing a layer of uncertainty. Further, the distinct
characteristic of traffic risk data, namely zero-inflation and overdispersion,
is often inadequately addressed by conventional models when assimilating
spatiotemporal risk features.

To address these gaps, we proposed the STZITD-GNNs model - a methodology
purposefully designed to bridge the demand for high-resolution, road-level
prediction and the requirement for rigorous uncertainty quantification. It
takes into account the unique distribution attributes of risk data and,
consequently, delivers predictions that are not only more precise and
efficient, but also contributive to a nuanced understanding of uncertainties
in traffic risks and the effectiveness of mitigation strategies. The primary
contributions of this study are encapsulated in the following points:

  * •

Guided by the overdispersed and zero-inflated nature of incident risk data, we
extended the conventional TD Distribution to a ZITD. This extension enhanced
the model’s ability to accurately capture the unique features of road-level
incident risk data. By introducing a sparsity parameter, $\pi$, and especially
integrating its effects into the dispersion parameter, $\phi$, our model
significantly bolsters its capacity to accurately depict the substantial
proportion of zero-valued (non-incident) roads. Simultaneously, it assures a
firm estimation of various risk levels, ranging from low-risk to long-tail
middle and high-risk roads.

  * •

Further, by integrating the ZITD encoder into the ST-GNNs, our model’s
capabilities are significantly elevated from traditional prediction to
pioneering UQ. This inventive amalgamation not only allows for effective
capture of the complex spatiotemporal dynamics inherent to urban road-level
risk but also adeptly transmits the distribution parameters to the output,
thereby achieving superior mapping of output variance and intervals compared
to traditional point estimation methodologies.

  * •

The proposed model showcased remarkable predictive accuracy across three
boroughs of London, underscoring its effectiveness across different urban
environments, each exhibiting unique risk patterns. Notably, our model
achieved high hit rates and ensured the greatest possible coverage of true
risk values within its prediction intervals across diverse time windows and
risk levels. This also underscores the model’s inherent resilience and
adaptability across a range of risk scenarios.

  * •

Lastly, our model adeptly handles the traditionally challenging trade-off
between MPIW and PICP when mapping the UQ, providing both precise and
confident risk predictions. By employing the zero-inflation component, it
achieves narrower prediction intervals (reduced MPIW) while ensuring
comprehensive coverage of actual risk data points (high PICP). This attribute
makes our model particularly valuable for practical applications, where both
precision and confidence in prediction intervals are critical.

The proposed STZITD-GNNs model, with its unique combination of spatial,
temporal, and statistical modelling capabilities, provides a comprehensive and
innovative solution for predicting and understanding urban road-level risks
and the prediction of UQ.

While we are encouraged by the results of our current research, we also
acknowledge that there remains scope for further refinement of our work. In
our pursuit to map risk values more accurately, one of our future endeavours
involves expanding the model’s applicability to different countries. Rather
than adhering to the current individual training processes, we plan to develop
a more versatile model that can be fine-tuned in response to geographic
heterogeneities, thereby extending its transferable applicability to diverse
areas.

Furthermore, we are keen to explore the integration of additional data sources
into our model. Specifically, the incorporation of GPS data could provide
valuable traffic exposure information, thereby further refining the risk
predictions.

In our continued commitment to improving the accuracy and utility of our
model, we also plan to evaluate our methodology using more comprehensive
datasets. For example, incorporating urban characteristics extracted from
street images or satellite imagery could enrich our model’s contextual
understanding, potentially highlighting factors that exacerbate the impacts of
traffic incidents in both spatial and temporal dimensions.

In conclusion, while we have made promising strides with the STZITD-GNNs
model, we view this as the beginning of an ongoing journey towards perfecting
urban road-level risk prediction. By continually refining our model and
integrating new sources of data, we aim to provide a tool that can greatly
assist in the mitigation of traffic risks and the enhancement of urban safety
strategies worldwide.

##  6 Declaration of Generative AI and AI-assisted technologies in the writing
process

During the preparation of this work the author(s) utilized GPT-4 for
preliminary grammar checks and proofreading. After using this tool/service,
the author(s) reviewed and edited the content as needed and take(s) full
responsibility for the content of the publication.

## References

  * Abe and Yadohisa, (2017) Abe, H. and Yadohisa, H. (2017).  A non-negative matrix factorization model based on the zero-inflated tweedie distribution.  Computational Statistics, 32(2):475–499.
  * Adepeju et al., (2016) Adepeju, M., Rosser, G., and Cheng, T. (2016).  Novel evaluation metrics for sparse spatio-temporal point process hotspot predictions-a crime case study.  International Journal of Geographical Information Science, 30(11):2133–2154.
  * Anderson, (2009) Anderson, T. K. (2009).  Kernel density estimation and k-means clustering to profile road accident hotspots.  Accident Analysis & Prevention, 41(3):359–364.
  * Bao et al., (2019) Bao, J., Liu, P., and Ukkusuri, S. V. (2019).  A spatiotemporal deep learning approach for citywide short-term crash risk prediction with multi-source data.  Accident Analysis & Prevention, 122:239–254.
  * Barahimi et al., (2021) Barahimi, A. H., Eydi, A., and Aghaie, A. (2021).  Multi-modal urban transit network design considering reliability: multi-objective bi-level optimization.  Reliability Engineering & System Safety, 216:107922.
  * Benlagha and Charfeddine, (2020) Benlagha, N. and Charfeddine, L. (2020).  Risk factors of road accident severity and the development of a new system for prevention: New insights from china.  Accident Analysis & Prevention, 136:105411.
  * Bergel-Hayat et al., (2013) Bergel-Hayat, R., Debbarh, M., Antoniou, C., and Yannis, G. (2013).  Explaining the road accident risk: Weather effects.  Accident Analysis & Prevention, 60:456–465.
  * Bibri, (2018) Bibri, S. E. (2018).  A foundational framework for smart sustainable city development: Theoretical, disciplinary, and discursive dimensions and their synergies.  Sustainable Cities and Society, 38:758–794.
  * Caliendo et al., (2007) Caliendo, C., Guida, M., and Parisi, A. (2007).  A crash-prediction model for multilane roads.  Accident Analysis & Prevention, 39(4):657–670.
  * Chakraborty et al., (2019) Chakraborty, P., Hegde, C., and Sharma, A. (2019).  Data-driven parallelizable traffic incident detection using spatio-temporally denoised robust thresholds.  Transportation research part C: emerging technologies, 105:81–99.
  * Chang and Chen, (2005) Chang, L.-Y. and Chen, W.-C. (2005).  Data mining of tree-based models to analyze freeway accident frequency.  Journal of safety research, 36(4):365–375.
  * Chen et al., (2018) Chen, C., Fan, X., Zheng, C., Xiao, L., Cheng, M., and Wang, C. (2018).  Sdcae: Stack denoising convolutional autoencoder model for accident risk prediction via traffic big data.  In 2018 sixth international conference on advanced cloud and big data (CBD), pages 328–333. IEEE.
  * Chen et al., (2016) Chen, Q., Song, X., Yamada, H., and Shibasaki, R. (2016).  Learning deep representation from big and heterogeneous data for traffic accident inference.  Proceedings of the AAAI Conference on Artificial Intelligence, 30(1).
  * Chen et al., (2020) Chen, Y., Kang, Y., Chen, Y., and Wang, Z. (2020).  Probabilistic forecasting with temporal convolutional neural network.  Neurocomputing, 399:491–501.
  * Cheng et al., (2022) Cheng, T., Zhang, Y., and Haworth, J. (2022).  Network spacetime ai: Concepts, methods and applications.  Journal of Geodesy and Geoinformation Science, 5(3):78.
  * Chong et al., (2004) Chong, M. M., Abraham, A., and Paprzycki, M. (2004).  Traffic accident analysis using decision trees and neural networks.  arXiv preprint cs/0405050.
  * Chung et al., (2014) Chung, J., Gulcehre, C., Cho, K., and Bengio, Y. (2014).  Empirical evaluation of gated recurrent neural networks on sequence modeling.  arXiv preprint arXiv:1412.3555.
  * Dai et al., (2018) Dai, L., Derudder, B., and Liu, X. (2018).  The evolving structure of the southeast asian air transport network through the lens of complex networks, 1979–2012.  Journal of Transport Geography, 68:67–77.
  * Ding et al., (2023) Ding, C., Sun, S., and Zhao, J. (2023).  Mst-gat: A multimodal spatial–temporal graph attention network for time series anomaly detection.  Information Fusion, 89:527–536.
  * Dixit and Sivakumar, (2020) Dixit, M. and Sivakumar, A. (2020).  Capturing the impact of individual characteristics on transport accessibility and equity analysis.  Transportation research part D: transport and environment, 87:102473.
  * Dunn, (2004) Dunn, P. K. (2004).  Occurrence and quantity of precipitation can be modelled simultaneously.  International Journal of Climatology: A Journal of the Royal Meteorological Society, 24(10):1231–1239.
  * Gawlikowski et al., (2021) Gawlikowski, J., Tassi, C. R. N., Ali, M., Lee, J., Humt, M., Feng, J., Kruspe, A., Triebel, R., Jung, P., Roscher, R., et al. (2021).  A survey of uncertainty in deep neural networks.  arXiv preprint arXiv:2107.03342.
  * Gu et al., (2019) Gu, Y., Lu, W., Qin, L., Li, M., and Shao, Z. (2019).  Short-term prediction of lane-level traffic speeds: A fusion deep learning model.  Transportation research part C: emerging technologies, 106:1–16.
  * Guo et al., (2014) Guo, J., Huang, W., and Williams, B. M. (2014).  Adaptive kalman filter approach for stochastic short-term traffic flow rate prediction and uncertainty quantification.  Transportation Research Part C: Emerging Technologies, 43:50–64.
  * Halder et al., (2019) Halder, A., Mohammed, S., Chen, K., and Dey, D. (2019).  Spatial risk estimation in tweedie compound poisson double generalized linear models.  arXiv preprint arXiv:1912.12356.
  * Hamilton, (2020) Hamilton, W. L. (2020).  Graph representation learning.  Synthesis Lectures on Artifical Intelligence and Machine Learning, 14(3):1–159.
  * He and Jiang, (2023) He, W. and Jiang, Z. (2023).  A survey on uncertainty quantification methods for deep neural networks: An uncertainty source perspective.  arXiv preprint arXiv:2302.13425.
  * Hossain and Muromachi, (2012) Hossain, M. and Muromachi, Y. (2012).  A bayesian network based framework for real-time crash prediction on the basic freeway segments of urban expressways.  Accident Analysis & Prevention, 45:373–381.
  * Hu et al., (2010) Hu, S.-R., Li, C.-S., and Lee, C.-K. (2010).  Investigation of key factors for accident severity at railroad grade crossings by using a logit model.  Safety science, 48(2):186–194.
  * Huang et al., (2019) Huang, C., Zhang, C., Dai, P., and Bo, L. (2019).  Deep dynamic fusion network for traffic accident forecasting.  In Proceedings of the 28th ACM international conference on information and knowledge management, pages 2673–2681.
  * Ibrahim et al., (2021) Ibrahim, M. R., Haworth, J., Christie, N., Cheng, T., and Hailes, S. (2021).  Cycling near misses: a review of the current methods, challenges and the potential of an ai-embedded system.  Transport reviews, 41(3):304–328.
  * Jiang et al., (2023) Jiang, X., Zhuang, D., Zhang, X., Chen, H., Luo, J., and Gao, X. (2023).  Uncertainty quantification via spatial-temporal tweedie model for zero-inflated and long-tail travel demand prediction.  arXiv preprint arXiv:2306.09882.
  * Jørgensen, (1987) Jørgensen, B. (1987).  Exponential dispersion models.  Journal of the Royal Statistical Society: Series B (Methodological), 49(2):127–145.
  * Khosravi and Nahavandi, (2014) Khosravi, A. and Nahavandi, S. (2014).  An optimized mean variance estimation method for uncertainty quantification of wind power forecasts.  International Journal of Electrical Power & Energy Systems, 61:446–454.
  * Khosravi et al., (2010) Khosravi, A., Nahavandi, S., Creighton, D., and Atiya, A. F. (2010).  Lower upper bound estimation method for construction of neural network-based prediction intervals.  IEEE transactions on neural networks, 22(3):337–346.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * (37) Kipf, T. N. and Welling, M. (2016a).  Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907.
  * (38) Kipf, T. N. and Welling, M. (2016b).  Variational graph auto-encoders.  arXiv preprint arXiv:1611.07308.
  * Koenker and Hallock, (2001) Koenker, R. and Hallock, K. F. (2001).  Quantile regression.  Journal of economic perspectives, 15(4):143–156.
  * Kong et al., (2023) Kong, J., Fan, X., Jin, X., Lin, S., and Zuo, M. (2023).  A variational bayesian inference-based en-decoder framework for traffic flow prediction.  IEEE Transactions on Intelligent Transportation Systems.
  * Kurz, (2017) Kurz, C. F. (2017).  Tweedie distributions for fitting semicontinuous health care utilization cost data.  BMC Medical Research Methodology, 17:1–8.
  * Laureshyn et al., (2017) Laureshyn, A., de Goede, M., Saunier, N., and Fyhri, A. (2017).  Cross-comparison of three surrogate safety methods to diagnose cyclist safety problems at intersections in norway.  Accident Analysis & Prevention, 105:11–20.
  * (43) Li, R., Zhong, T., Jiang, X., Trajcevski, G., Wu, J., and Zhou, F. (2022a).  Mining spatio-temporal relations via self-paced graph contrastive learning.  In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pages 936–944.
  * (44) Li, R., Zhong, T., Jiang, X., Trajcevski, G., Wu, J., and Zhou, F. (2022b).  Mining spatio-temporal relations via self-paced graph contrastive learning.  In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD ’22, page 936–944.
  * (45) Li, Y., Chai, S., Wang, G., Zhang, X., and Qiu, J. (2022c).  Quantifying the uncertainty in long-term traffic prediction based on pi-convlstm network.  IEEE Transactions on Intelligent Transportation Systems, 23(11):20429–20441.
  * Li et al., (2017) Li, Y., Yu, R., Shahabi, C., and Liu, Y. (2017).  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  arXiv preprint arXiv:1707.01926.
  * Lin et al., (2018) Lin, L., Handley, J. C., Gu, Y., Zhu, L., Wen, X., and Sadek, A. W. (2018).  Quantifying uncertainty in short-term traffic prediction and its application to optimal staffing plan development.  Transportation Research Part C: Emerging Technologies, 92:323–348.
  * Lin et al., (2015) Lin, L., Wang, Q., and Sadek, A. W. (2015).  A novel variable selection method based on frequent pattern tree for real-time traffic accident risk prediction.  Transportation Research Part C: Emerging Technologies, 55:444–459.
  * Lin et al., (2021) Lin, Y. T., Neumann, J., Miller, E. F., Posner, R. G., Mallela, A., Safta, C., Ray, J., Thakur, G., Chinthavali, S., and Hlavacek, W. S. (2021).  Daily forecasting of regional epidemics of coronavirus disease with bayesian uncertainty quantification, united states.  Emerging infectious diseases, 27(3):767.
  * Liu et al., (2022) Liu, G., Jin, H., Li, J., Hu, X., and Li, J. (2022).  A bayesian deep learning method for freeway incident detection with uncertainty quantification.  Accident Analysis & Prevention, 176:106796.
  * Liu and Guan, (2004) Liu, J. and Guan, W. (2004).  A summary of traffic flow forecasting methods.  Journal of highway and transportation research and development, 21(3):82–85.
  * Liu et al., (2020) Liu, J., Li, T., Xie, P., Du, S., Teng, F., and Yang, X. (2020).  Urban big data fusion based on deep learning: An overview.  Information Fusion, 53:123–133.