Uncertainty Quantification of Spatiotemporal Travel Demand with
  Probabilistic Graph Neural Networks
Recent studies have significantly improved the prediction accuracy of travel
demand using graph neural networks. However, these studies largely ignored
uncertainty that inevitably exists in travel demand prediction. To fill this
gap, this study proposes a framework of probabilistic graph neural networks
(Prob-GNN) to quantify the spatiotemporal uncertainty of travel demand. This
Prob-GNN framework is substantiated by deterministic and probabilistic
assumptions, and empirically applied to the task of predicting the transit and
ridesharing demand in Chicago. We found that the probabilistic assumptions
(e.g. distribution tail, support) have a greater impact on uncertainty
prediction than the deterministic ones (e.g. deep modules, depth). Among the
family of Prob-GNNs, the GNNs with truncated Gaussian and Laplace distributions
achieve the highest performance in transit and ridesharing data. Even under
significant domain shifts, Prob-GNNs can predict the ridership uncertainty in a
stable manner, when the models are trained on pre-COVID data and tested across
multiple periods during and after the COVID-19 pandemic. Prob-GNNs also reveal
the spatiotemporal pattern of uncertainty, which is concentrated on the
afternoon peak hours and the areas with large travel volumes. Overall, our
findings highlight the importance of incorporating randomness into deep
learning for spatiotemporal ridership prediction. Future research should
continue to investigate versatile probabilistic assumptions to capture
behavioral randomness, and further develop methods to quantify uncertainty to
build resilient cities.
Qingyi Wang, Shenhao Wang, Dingyi Zhuang, Haris Koutsopoulos, Jinhua Zhao
# Uncertainty Quantification of Spatiotemporal Travel Demand with Probabilistic Graph Neural Networks

###### Abstract

Recent studies have significantly improved the prediction accuracy of travel demand using graph neural networks. However, these studies largely ignored uncertainty that inevitably exists in travel demand prediction. To fill this gap, this study proposes a framework of probabilistic graph neural networks (Prob-GNN) to quantify the spatiotemporal uncertainty of travel demand. This Prob-GNN framework is substantiated by deterministic and probabilistic assumptions, and empirically applied to the task of predicting the transit and ridesharing demand in Chicago. We found that the probabilistic assumptions (e.g. distribution tail, support) have a greater impact on uncertainty prediction than the deterministic ones (e.g. deep modules, depth). Among the family of Prob-GNNs, the GNNs with truncated Gaussian and Laplace distributions achieve the highest performance in transit and ridesharing data. Even under significant domain shifts, Prob-GNNs can predict the ridership uncertainty in a stable manner, when the models are trained on pre-COVID data and tested across multiple periods during and after the COVID-19 pandemic. Prob-GNNs also reveal the spatiotemporal pattern of uncertainty, which is concentrated on the afternoon peak hours and the areas with large travel volumes. Overall, our findings highlight the importance of incorporating randomness into deep learning for spatiotemporal ridership prediction. Future research should continue to investigate versatile probabilistic assumptions to capture behavioral randomness, and further develop methods to quantify uncertainty to build resilient cities.

 probabilistic graph neural networks, uncertainty quantification, travel demand prediction

## I Introduction

Uncertainty prevails in urban mobility systems. An enormous number of uncertainty sources range from the long-term natural disasters and climate change, to real-time system breakdown, and to the inherent randomness in travel demand. These uncertainty sources exhibit spatial and temporal patterns: travel demand could drastically change during the holiday seasons, or becomes spatially concentrated in the special events (e.g. football games). Traditionally, the spatial uncertainty is analyzed by spatial econometrics or discrete choice models, while the temporal one by time series models. Recently, deep learning models have been designed to capture the spatiotemporal correlations, with the Long Short Term Memory (LSTM) networks for temporal correlations and convolutional or graph networks (CNN, GNN) for spatial dependencies. Deep learning has been shown to significantly improve the prediction accuracy of travel demand. However, most studies focus on designing _deterministic_ deep learning models to predict the _average_ travel demand, but largely ignore the uncertainty inevitably associated with any prediction.

Overlooking uncertainty leads to theoretical and practical problems. In theory, travel demand is inherently a random quantity, which is characterized by a rich family of probability distributions, sometimes even "fat-tail" ones that render a simple average approximation highly inappropriate [1]. In practice, while the prediction of average demand provides a valuable basis for system design, the lack of uncertainty analysis precludes the opportunity of using deep learning to provide robust real-time service or resilient long-term planning [2, 3]. Since the recent work has demonstrated the outstanding capability of deep learning in modeling spatiotemporal dynamics, it seems timely and imperative to enrich the deterministic deep learning models to capture the spatiotemporal uncertainty of travel demand.

To address this research gap, we propose a framework of probabilistic graph neural networks (Prob-GNN) to quantify the spatiotemporal uncertainty in travel demand. This framework is substantiated by deterministic and probabilistic assumptions: the former refer to the various architectural designs in deterministic neural networks, while the latter refer to the probabilistic distributions with a variety of uncertainty characteristics. The framework is used to compare the two types of assumptions by both uncertainty and point prediction metrics. This study makes the following contributions:

1. We identify the research gap in quantifying uncertainty for travel demand using deep learning. Despite a vast number of deep learning techniques for predicting average travel demand, only a limited number of studies have sought to quantify spatiotemporal uncertainty.
2. To fill the research gap, the paper introduces a general framework of Prob-GNN, which comprises deterministic and probabilistic assumptions. This framework is further substantiated by twelve architectures, all of which can quantify the spatiotemporal uncertainty in travel demand.
3. We find that probabilistic assumptions significantly influence uncertainty predictions, while only slightly influence point predictions. Conversely, the deterministic assumptions, such as graph convolution, network depth, and dropout rate, lead to a similar performance in both point and uncertainty predictions when probabilistic assumptions are fixed.
4. We further demonstrate the generalizability of the Prob-GNNs by comparing the performance across multiple
shifting domains caused by COVID-19. Although the mean prediction encounters significant prediction errors, the interval prediction is accurate across the domains.

## II Literature Review

The literature on spatiotemporal travel demand prediction using deep learning for has grown significantly in recent years. Among them, very limited number of studies have sought to understand the predictive uncertainty. This literature review will first examine the latest developments in demand predictions using deep learning, and subsequently, delve into uncertainty quantification.

### _Travel Demand Prediction with Deep Learning_

Researchers have achieved a lot of success in improving the prediction accuracy of spatiotemporal travel demand using advanced neural network architectural designs [4, 5, 6, 7]. Temporally, LSTM networks and Gated Recurrent Units (GRU) layers are used to analyze the seasonal, weekly, and time-of-day trends [8, 9, 10]. Spatially, the analysis unit is often station/stop, urban grid, individuals, or census tract. The urban grid is often analyzed with convolutional neural networks (CNN) [4]. Individual travel demand is analyzed by neural networks with behavioral insights for innovative architectural design [11, 12, 13]. More complex urban structures, such as transit lines that connect upstream and downstream stops, can be represented by graphs. Graph convolutions (GCN) are then often used to model such spatial propagation of traffic/demand flow into adjacent nodes [14, 15, 16, 17]. The baseline GCNs have been expanded by defining multiple graphs to describe the spatial similarity in points of interest (POIs) and transportation connectivity [18, 19]. An alternative to GCN is Graph Attention Network (GAT), which automatically learn the attentional weightings of each node, thus detecting the long-range dependencies without handcrafted weighting mechanisms through adjacency matrices [20, 21].

Despite the abundance of deep learning models, few studies sought to quantify the travel demand uncertainty, typically caused by narrowly defined prediction objectives and methods. The prediction objective - the average ridership - can hardly represent the randomness in travel demand, particularly for distributions with fat tails, in which the outliers deviating from the average value could predominate [1]. The prediction method - the family of deterministic neural networks - is only a specific example of probabilistic neural networks with homogeneous assumptions on the variance. Demand uncertainty, if unaccounted for, can negatively influence downstream planning tasks [3, 22]. Uncertainty can also propagate and be significantly enhanced at each stage in multi-stage models [23]. Despite its importance, in the travel prediction field, only three studies tackled uncertainty with deep learning in transportation. One analyzed heteroskedastic noises with quantile regressions [24], one studied the uncertainty in sparse travel demand prediction [25], and the other quantified model uncertainty with Monte-Carlo dropouts [26].

### _Uncertainty Quantification Methods_

Research on uncertainty quantification has been developed in other deep learning fields. There are two major categories of uncertainty: data uncertainty and model uncertainty. Data uncertainty refers to the irreducible uncertainty inherent in the data generation process, while model uncertainty captures the uncertainties in the model parameters. For example, in linear regression, data uncertainty refers to the residuals, and model uncertainty refers to the standard errors of the estimated coefficients.

To characterize data uncertainty, either a parametric or non-parametric method can be used. Parametric methods refer to the models that parameterize a probabilistic distribution. Parametric models are often estimated by either Bayesian methods or mean-variance estimation (MVE). Although conceptually appealing, Bayesian methods often involves unnecessarily intense computation, which relies on sampling methods and variational inference [27, 28]. MVE minimizes the negative log-likelihood (NLL) loss based on a pre-specified distribution of the dependent variable [29, 30]. The MVE is computationally efficient, but it could lead to misleading results if probabilistic distributions are misspecified. Non-parametric methods quantify uncertainty without explicitly imposing any parametric form. Typical examples include quantile regression [31] and Lower Upper Bound Estimation (LUBE) [27, 28]. Non-parametric methods do not have misspecification problems, but it is hard to be optimized in complex neural networks. The pros and cons of both methods are systematically compared in review articles [32, 33].

Besides data uncertainty, uncertainty also emerges in model construction and training. Bayesian deep neural networks refer to the models that characterize all parameters of neural networks by probability distributions, typically Gaussian distribution [25, 34]. A popular, low-cost way to do approximate Bayesian inference for neural networks with Gaussian parameters is to use dropouts [26, 35, 36]. Another source of model uncertainty comes from model training. Neural networks tend to converge to different local minima, resulting in different predictions. This type of uncertainty is usually addressed by bootstrapping or model ensembles [37, 38]. Bootstrapping refers to the process of repeated training with re-sampling to capture the distributional uncertainty in sampling. Model ensembling refers to an average of multiple trainings procedures with different parameter initializations.

Although the literature on travel demand uncertainty within the deep learning framework is scarce, classical statistics have decades of effort in quantifying travel demand uncertainty, with a focus on evaluating a rich family of probabilistic assumptions. Using high-resolution temporal ridership data, researchers adopted the ARIMA-GARCH model to quantify the volatility of subway demands during special events [39] and adaptive Kalman filters to analyze real-time transit ridership [40]. Using cross-sectional data, researchers used bootstrapping [41] to analyze parameter uncertainty, ensembling [42] to analyze activity uncertainty, and heteroskedastic errors to describe the association between social demographics and ridership uncertainty in discrete choice models [41].
The common debates regarding the probabilistic assumptions include homoskedastic vs. heteroskedastic errors, Gaussian vs. exponential tails, real-line vs. positive support, and many others [39, 43, 44], through which researchers could learn the detailed characteristics of travel demand uncertainty. The primary focus from the classical statistics on the probabilistic assumptions presents an interesting difference from the primary focus of the deep learning field in refining the deterministic assumptions. Indeed, it is quite ambiguous which set of assumptions is more crucial in quantifying the uncertainty of travel demand.

## III Theoretical Framework

The authors propose a probabilistic graph convolutional neural network to compare the deterministic and probabilistic assumptions. The probabilistic GNN is represented as:

\[y\sim\mathcal{G}(\theta)=\mathcal{G}(\mathcal{F}(\mathbf{X},w)) \tag{1}\]

The first statement is a probability assumption about \(y\), specified by the model family \(\mathcal{G}(\theta)\), and the second uses the model family \(\mathcal{F}(\mathbf{X},w)\) to parameterize the probability distributions of \(\mathcal{G}(\theta)\), with \(\mathbf{X}\) being the inputs and \(w\) being the model weights. While the recent studies focus on enriching \(\mathcal{F}(\mathbf{X},w)\) through spatiotemporal deep learning architectures, the potentially more critical probability assumption \(\mathcal{G}(\theta)\) is largely neglected. The authors seek to compare the effects of the probabilistic assumptions \(\mathcal{G}\), and the deterministic assumptions \(\mathcal{F}\) in determining the model performance.

Table I summarizes the specifics of probabilistic and deterministic assumptions that substantiate the probabilistic GCN framework. A cross product of six probabilistic assumptions in \(\mathcal{G}\) and two deterministic assumptions in \(\mathcal{F}\) are tested, leading to twelve base models. Six probabilistic assumptions are: Homoskedastic Gaussian (HomoG), Poisson (Pois), Heteroskedastic Gaussian (HetG), Truncated Gaussian (TG), Gaussian Ensemble (GEns), and Laplace (LAP). Two deterministic architectures, Graph Convolutional Networks (GCN) and Graph Attention Networks (GAT) are used to specify the probabilistic parameters.

### _Probabilistic Assumptions \(\mathcal{G}\)_

The Gaussian distribution with a deterministic and homoskedastic variance term is chosen as the benchmark in specifying \(\mathcal{G}\) because it is stable and has simple ensembling properties. This homoskedastic Gaussian assumption represents the vast number of deterministic deep learning models that use mean squared error as the training objective. The Gaussian benchmark facilitates the comparison across the probabilistic assumptions, as listed below.

1. Homoskedasticity vs. heteroskedasticity We compare the heteroskedastic and homoskedastic Gaussian assumptions to examine how the data variance influences the model performance. The homoskedastic assumption assumes the same variance for all observations, whereas the heteroskedastic assumption estimates variance for every observation. It is highly likely that the travel demand variances has spatiotemporal patterns.
2. Continuous vs. discrete support We compare the Poisson distribution to the homoskedastic Gaussian benchmark to examine the effectiveness of continuous vs. discrete supports in determining the model performance. Since ridership takes integer values, the Poisson distribution could be more appropriate. However, the Poisson distribution uses only one parameter to represent both mean and variance, which could be overly restrictive.
3. Real-line vs. non-negative support We compare the truncated heteroskedastic Gaussian distribution to the heteroskedastic Gaussian distribution to examine the effectiveness of the real-line vs. non-negative distribution support. Travel demand is non-negative, but the support of the Gaussian distribution covers the entire real line \((-\infty,\infty)\). Therefore, a normal distribution left-truncated at zero is implemented to test whether non-negativity should be strictly imposed.
4. Gaussian vs. exponential tails We compare the Laplace distribution to the heteroskedastic Gaussian distribution to examine whether the tail behavior of the distribution matters. The probability density function of the Gaussian distribution decays at a fast rate of \(e^{x^{2}}\), while the Laplace distribution has a heavier tail with the decay rate of \(e^{|x|}\).
5. Single vs. ensembled models We compare the ensembled heteroskedastic Gaussian distributions to a single heteroskedastic Gaussian model to test whether an ensemble of distributions can outperform a single distributional assumption. The ensemble model is created by uniformly mixing \(K\) estimates, which are trained independently with different parameter initializations. Suppose \(\mathcal{Y}_{k},\sigma_{k}^{2}\) are the estimated mean and variance for model \(k\), and \(\mathcal{Y}_{*}\) is the ensembled random variable. For the Gaussian distribution, the mixture can be further approximated as a Gaussian distribution, whose mean and variance are respectively the mean and variance of the mixture. The ensembled distribution is given by [38]: \(\mathcal{Y}_{*}\sim\mathcal{N}(\frac{1}{K}\sum_{k}\mathcal{Y}_{k},\frac{1}{K} \sum_{k}\left(\sigma_{k}^{2}+\mathcal{Y}_{k}^{2}\right)-\mathcal{Y}_{*}^{2})\).

With the distributional assumption, maximum likelihood estimation (MLE) can be used to learn the parameters, which translates to minimizing the negative log-likelihood (NLL) loss in implementation. The NLL loss function based on the joint density of all observations is simply the negative sum of the log of the probability density of all observations:

\[\text{NLL}=-\sum_{s,t}\text{log }\mathcal{G}(y_{st}|\mathbf{X_{t}};w) \tag{2}\]

In the uncertainty quantification literature, researchers typically differentiate between data (aleator) and model (epistemic) uncertainty [47], and their sum is referred to as prediction uncertainty. The model uncertainty refers to the uncertainty resulting from the difference between \(\mathcal{F}(\mathbf{X},w)\) and the estimated \(\hat{\mathcal{F}}(\mathbf{X},w)\). The data uncertainty refers to the randomness in the data generation process and is represented by probabilistic assumptions \(\mathcal{G}\), of which \(y\) has the distribution. The prediction uncertainty combines model and
data uncertainty, and describes the overall difference between the actual \(y\) and the predicted \(\hat{y}\). The relationship between the three quantities is directly given by \(\sigma_{y}^{2}=\sigma_{model}^{2}+\sigma_{data}^{2}\). In our framework, assumptions 1-4 deal with the data uncertainty alone; while assumption 5 quantifies the prediction uncertainty.

### _Deterministic Assumptions in \(\mathcal{F}\)_

The deterministic assumptions in \(\mathcal{F}\) consist of the spatial encoding, temporal encoding, and associated hyperparameter specifications. The following two subsections introduce the formulation of common spatial and temporal encodings.

#### Iii-B1 GCN and GAT for Spatial Encoding

Two predominant spatial encoding methods - GCN and GAT - are adopted and compared to examine the effect of the deterministic architectures on model performance. The GCN layers need to access the global structure of the graph in the form of adjacency matrice(s), while the GAT layers aim to learn the spatial correlations from data. The propagation formula of both layers is introduced in Table I.

To construct the multi-graph proximity in GCNs and GATs, four types of adjacency matrices are computed, including direct connectivity \([A_{Con}]_{ij}=1\) if two stations are adjacent, \(=0\) otherwise, network distance \([A_{Net}]_{ij}=\text{Network Distance}(i,j)^{-1}\), Euclidean distance \([A_{Euc}]_{ij}=\text{Euclidean Distance}(i,j)^{-1}\), and functional similarity \([A_{Func}]_{ij}=\sqrt{(F_{i}-F_{j})^{T}(F_{i}-F_{j})}^{-1}\), where \(F_{i}\) is the vector of functionalities represented by population, jobs, percent low income, percent minority, percent adults, number of schools, shops, restaurants, etc.

#### Iii-B2 LSTM for Temporal Encoding

The LSTM network is chosen for temporal encoding. In short, the LSTM layers take the spatially encoded tensors for \(l\) previous time periods as inputs and propagate them down the layers through a series of input, forget, cell state (memory), and output gates. To make the prediction, the encoded output from the last time step \(h_{t-1}\) is decoded by fully connected layers to get the contribution (weights) of the time series on the final prediction. The LSTM structure is well-documented in a lot of papers [24, 39], and will not be discussed in detail here.

### _Specific Examples: Gaussian Distributions and Mean-Variance Estimation_

As an important case, the Gaussian distribution with homoskedastic variance (HomoG-GCN and HomoG-GAT) represents the vast number of deterministic deep learning models that use mean squared errors as the training objective. With the homoskedastic Gaussian assumption, the dependent variable \(y\) follows Gaussian distribution with mean \(F_{1}(\mathbf{X},w_{1})\) and a constant variance \(c\). The MLE with this homoskedastic Gaussian distribution is the same as the deterministic deep learning with the mean squared errors as the training objective because

\[\log P(x;\mu,\sigma=c)=\log\frac{1}{c\sqrt{2\pi}}-\frac{1}{2}\frac{(x-\mu)^{2} }{c^{2}} \tag{3}\]

In other words, our benchmark example represents the predominant deterministic modeling technique in this field.

The homoskedastic Gaussian can be extended to the heterskedastic Gaussian distribution, which resembles the mean-variance estimation (MVE), the most dominant method in uncertainty quantification literature. With the heteroskedastic Gaussian distribution, \(y\sim\mathcal{N}(\theta)=\mathcal{N}(F_{1}(\mathbf{X},w_{1}),F_{2}(\mathbf{X},w_{2}))=F_{1}(\mathbf{X},w_{1})+\mathcal{N}(0,F_{2}(\mathbf{X},w_{2}))\). Essentially, the MVE uses two graph neural networks to capture the mean and variance separately, which is the same as the HetG-GCN and HetG-GAT models. Therefore, the two Gaussian examples in our probabilistic GNN framework can represent the two most common research methods, namely the deterministic spatiotemporal models and the MVE for uncertainty quantification.

\begin{table}
\begin{tabular}{l|c|c} \hline \hline
**Panel 1: Probabilistic Assumptions in \(\mathcal{G}\)** & \multicolumn{1}{c}{} \\ \hline Probabilistic assumptions & Probability density function & Distribution parameters \(\theta\) \\ \hline Homoskedastic Gaussian (HomoG) & \(f_{HomoG}(x;\mu,\sigma)=\frac{1}{c\sqrt{2\pi}}\text{exp}(-\frac{1}{2}\frac{(x- \mu)^{2}}{\sigma^{2}})\), \(\sigma=c\) & \(\mu\) \\ Poisson (Pois) & \(f(x;\lambda)=\frac{\lambda_{\sigma}^{2}-c^{2}}{a^{2}}\) & \(\mu\) \\ Heteroskedastic Gaussian (HetG) & \(f_{HetG}(x;\mu,\sigma)=\frac{1}{c\sqrt{2\pi}}(-\frac{1}{2}\frac{(x-\mu)^{2}}{ \sigma^{2}})\) & \(\mu,\sigma\) \\ Truncated Gaussian (TG) & \(f_{TC}(x;\mu,\sigma)=\frac{1}{c\sqrt{2\pi}}\frac{1-f_{(0;\mu,\sigma)}}{\sum_{ \lambda\in\mathcal{N}}\mu_{i}}\frac{1-f_{(0;\mu,\sigma)}}{\sum_{\lambda\in \mathcal{N}}\mu_{i}}\) & \(\mu,\sigma\) \\ Gaussian Ensemble (GEns) & \(y_{\star}\sim\mathcal{N}(\frac{1}{N}\sum_{k}\mu_{k},\frac{1}{N}\sum_{k}(\alpha _{k}^{2}+\frac{1}{k}\mu_{k}^{2})-\mu_{k}^{2})\) & \(\mu_{k},\sigma_{k}\) for \(k=1..K\) models \\ Laplace (Lap) & \(f(x;\mu,b)=\frac{1}{c}\text{exp}(-\frac{|x-\mu|}{b})\) & \(\mu,b\) \\ \hline \hline
**Panel 2: Deterministic Assumptions in \(\mathcal{F}\)** & \multicolumn{1}{c}{} \\ \hline Deterministic assumptions & Graph convolutional iteration function & Deterministic parameters \(w\) \\ \hline Graph Convolutional Networks (GCN) [45] & \(h^{(l+1)}=\sigma(\sum_{\gamma=1}^{R}\hat{D}_{\gamma}^{-\frac{1}{2}}\hat{A}_{ \gamma}\hat{D}_{\gamma}^{-\frac{1}{2}}h^{(l)}(W_{r}^{(l)})\) & \(W_{r}^{(l)}\) \\ Graph Attention Networks (GAT) [46] & \(h^{(l+1)}=\sigma(\sum_{j\in\mathcal{N}(0)}\alpha_{ij}W_{j}^{(l)})\), \(\alpha_{ij}=\frac{c\sigma_{W}(W_{j}^{(l)})}{\sum_{\lambda\in\mathcal{N}, exp}(W_{j}^{(l)})}\) & \(\alpha_{ij},W\) \\ \hline Other hyperparameters in deterministic assumptions \(\mathcal{F}\): Number of GCN/GAT/LSTM layers, Number of hidden layer neurons, Dropouts, Weight decay. \\ \hline \hline \end{tabular}
\end{table} TABLE I: Model Design. Upper panel: Probabilistic assumptions and their probability density functions. Bottom panel: Deterministic architecture and hyperparameters. A cross product of the two panels are used to substantiate the probabilistic GCNs, leading to twelve base models: HomoG-GCN, HomoG-GAT, HetG-GCN, HetG-GAT, TG-GCN, TG-GAT, Gems-GCN, GENs-GAT, Pois-GCN, Pois-GAT, Lap-GCN, and Lap-GAT.
### _Evaluation_

Performance metrics can be grouped into three categories: composite measures, point prediction quality, and uncertainty prediction quality. Table II summarizes the evaluation metrics. The NLL loss is a composite metric to evaluate the joint quality of point and uncertainty estimates. The standard mean absolute error (MAE) and mean absolute percent error (MAPE) are used to evaluate point predictions.

The quality of the uncertainty estimates is less straightforward to evaluate, because the ground truth distributions are unknown. Therefore, we design the calibration error (CE) metric and also visualize quantile-quantile plots to measure the distributional fit. Let \(F_{i}(y)\) represent the cumulative distribution function of observation \(i\). Define \(q(p)=\mathbb{P}(F_{i}(y_{i})\leq p)\) as the proportion of observations that actually fall into the \(p\)-th quantile of the estimated distribution. A well-calibrated model should generate distributions that align with the empirical distribution so that \(q(p)=p\). Therefore, the calibration error associated with the quantile-quantile plots is defined as the sum of the deviation between the empirical quantiles and the predicted quantiles, approximated by a number of discrete bins from \([0,1]\): \(CE=\sum_{p=0}^{1}|q(p)-p|\). An alternative metric is the simultaneous use of the Prediction Interval Coverage Probability (PICP) and the Mean Prediction Interval Width (MPIW) [27, 28, 32]. Formally, \(PICP=\frac{1}{n}\sum_{i=1}^{n}c_{i}\), where \(c_{i}=\mathds{1}\{L_{i}\leq y_{i}\leq U_{i}\}\), that is an indicator variable of whether observation \(i\) falls within the prediction interval bounded by \(L_{i}\) and \(U_{i}\). \(MPIW=\frac{1}{n}\sum_{i=1}^{n}\left(U_{i}-L_{i}\right)\) measures the average width of the intervals. With significance level \(1-\alpha\), the lower bound of the prediction interval is given by \(L=F_{i}^{-1}(\frac{\alpha}{2})\), and the upper bound of the prediction interval by \(U=F_{i}^{-1}(1-\frac{\alpha}{2})\), where \(F_{i}^{-1}\) is the predicted inverse cumulative function of observation \(i\). Using this approach, the evaluation of uncertainty quantification involves a tradeoff between PICP and MPIW. A high-quality prediction interval should be narrow while covering a large proportion of data points; however, wider MPIW naturally leads to larger PICP. Therefore, this tradeoff poses a challenge in model selection since it is difficult for one model to dominate in both metrics.

## IV Case Study

Two case studies were conducted to examine the effect of deterministic and probabilistic assumptions on predicting travel demand, with a focus on estimating uncertainty. The case studies use data from the Chicago Transit Authority's (CTA) rail system and ridesharing in Chicago. This section describes data sources, the spatial and temporal characteristics of the two datasets, and the experiment setup. Our experiments are implemented in PyTorch and the source code is available at [https://github.com/sunnyqwang/uncertainty](https://github.com/sunnyqwang/uncertainty).

### _Data_

The CTA rail and bus data was obtained through collaboration with the CTA, while the ridesharing data was sourced from the City of Chicago Open Data Portal1. The CTA rail system has tap-in records for each trip, which were aggregated into 15-minute intervals for temporal analysis. The system comprises 141 stations, arranged in a graph structure for spatial analysis. The spatial relationships can be identified by adjacency matrices constructed from Euclidean, connectivity, network and functional similarity between stations. Ridesharing trips are available at the census tract level, in 15-minute intervals. The graph structure is determined by the relationships between census tracts. As there is no explicit network, the ridesharing adjacency matrices are defined by Euclidean, connectivity (neighboring census tracts are considered connected), functional similarity only. Most census tracts have close to none ridesharing trips in most 15-minute intervals. Since learning sparse graphs are a topic of its own [48], the ridesharing dataset is filtered to include only those census tracts that, on average, had more than 30 trips per hour pre-COVID. This resulted in a total of 59 census tracts being used in the analysis.

Footnote 1: [https://data.cityofchicago.org/Transportation/Transportation-Network-Providers-Tripps/m6dm-c72p](https://data.cityofchicago.org/Transportation/Transportation-Network-Providers-Tripps/m6dm-c72p)

Moran's I was calculated for each 15min ridership snapshot using each of the weight matrices to demonstrate that spatial autocorrelation exists through the defined adjacency matrices. Figure 1 shows the histograms of Moran's I for both data sources. Most statistics have averages well above 0, indicating the existence of spatial autocorrelations defined by the adjacency matrices. Ridesharing trips appear to be less correlated with Euclidean distance and functional similarity, but the correlation between bordering tracts is very strong.

Five different explanatory variables are obtained for each data source. First, the system's own history, both recent (in the past 1.5 hours) and historical (last week same time period) is used to account for the temporal correlation. Second, the history of the other travel modes is included to mine the correlation between different modes. In the case of CTA rail, ridesharing and bus counts are included; and in the case of ridesharing, CTA rail and bus counts are included.

\begin{table}
\begin{tabular}{c|c|c} \hline \hline
**Category** & **Metric** & **Formula** \\ \hline Composite & Negative Log Likelihood & \(NLL=-\sum_{i}\text{log}\ P_{W}(y_{i}|\mathbf{X}_{i})\) \\ \hline \multirow{2}{*}{Point} & Mean Absolute Error & \(MAE=\frac{1}{4}\sum_{i}\text{log}\ \hat{y}_{i}(y_{i}-\hat{y}_{i})\) \\  & Mean Absolute Percent Error & \(MAE=\frac{1}{4}MAE/\hat{y}\) \\ \hline \multirow{3}{*}{Uncertainty} & Calibration Error & \(CE=\sum_{p=0}^{1}|q(p)-p|\) \\  & Mean PI Width & \(MPIW=\sum_{i=1}^{n}\left(U_{i}-L_{i}\right)\) \\ \cline{1-1}  & PI Coverage Probability & \(PICP=\frac{1}{N}\sum_{i=1}^{n}\mathds{1}\{L_{i}\leq y_{i}\leq U_{i}\}\) \\ \hline \hline \end{tabular}
\end{table} TABLE II: Three Categories of Evaluation Metrics

Fig. 1: Spatial Autocorrelation (Moran’s I)
Third, demographics2 and POIs 3 are spatial covariates used in the calculation of functional similarity. Fourth, weather4 is a temporal covariate that is assumed to be the same in the whole study area. Lastly, the frequency of services in each spatial unit during each time period is used to indicate supply levels. In the case of ridesharing, we do not have supply information. Instead, we use the bus frequency as a proxy. Bus schedules are slower to respond during normal times, but after March 2020 the CTA is actively adjusting service to account for labor shortages and changing demand.

Footnote 2: [https://www.census.gov/programs-surveys/acs/data.html](https://www.census.gov/programs-surveys/acs/data.html)

Footnote 3: [https://planet.openstreetmap.org/](https://planet.openstreetmap.org/)

Footnote 4: [https://www.ncdc.noaa.gov/cdo-web/](https://www.ncdc.noaa.gov/cdo-web/)

### _Experiment Setup_

The datasets were split into three subsets along the temporal dimension: train, validation, and test. The training set consists of data from August 1, 2019 to Feb 16, 2020, the validation set from Feb 17, 2020 to March 1, 2020, and the testing set consists of four different two-week periods during the course of COVID-19 pandemic in 2020: immediately before (March 2 to March 15), stay-at-home (March 16 to March 29), initial recovery (June 22 to July 5), and steady recovery (Oct 18 to Oct 31). Significant changes emerged starting March 2020 when the confirmed COVID-19 cases were growing rapidly and the city issued stay-at-home orders. Since then, the ridership has seen different stages of recovery. The changes in ridership of a few stations/census tracts are shown in Figure 2. Meanwhile, changes induced by COVID-19 provide an opportunity to test the temporal generalizability from pre- to post-COVID periods.

Several other architectural consideration have been taken into account besides the GCN/GAT spatial convolutions. Figure 3 summarizes the deterministic architecture, consisting of three components: spatiotemporal layers (GCN/GAT+LSTM) for recent observed demand, and linear layer to connect last week's observation (history), and weather (temperature and precipitation). The three components are indexed by \(r,h,p\), respectively, and summed to the final prediction. First, due to the cyclic nature of travel demand, the travel demand time series can be decomposed into a weekly component (\(\hat{\mathcal{Y}}_{t}^{h}\), \(\hat{\sigma}_{t}^{h}\)), and a time-of-day component (\(\hat{\mathcal{Y}}_{t}^{r}\), \(\hat{\sigma}_{t}^{r}\)). The weekly component is calculated from a reference demand. A good reference demand is the observed value for the same region and time in the previous week. The weights \(w_{h}\) are obtained from the LSTM encodings. Additionally, the recent demand used for LSTM network inputs are the deviations from last week's demand, and the decoded outputs are used to produce both the weekly and the time-of-day components. Intuitively, if the recent residual demand is very different from the reference, the reference should have a smaller weight on the final prediction, and the final prediction is more uncertain. The time-of-day component is directly obtained from LSTM encodings of the recent residual demand. Next, weather is another source of temporal variation because extreme weather could have an evident impact on transit ridership. The model takes daily deviations from average of precipitation \(P_{T}\) and temperature \(T_{E}\), and multiply them with spatiotemporal weights \(W^{p},W^{t}\in\mathbb{R}^{2\times T\times S}\) to get \([\hat{\mathcal{Y}}_{t}^{p},\hat{\sigma}_{t}^{p}]\).

For each dataset and each deterministic architecture, six probabilistic assumptions are tested: homoskedastic Gaussian (homoG), Poisson (Pois), heteroskedastic Gaussian (HetG), truncated Gaussian (TG), Laplace (Lap), ensembled Gaussian (GEns). In the homoG model, a search for the best standard deviation was done. Values are searched in multiples of \(\bar{y}\): \(\frac{1}{4},\frac{1}{2},\frac{3}{4},1\), and the best value for the variance is \(\frac{1}{2}\bar{y}\) for both datasets. In Gems, five top HetG models by validation set loss are selected to create an ensemble model. The ensemble distribution is given by [38]: \(\mathcal{Y}_{*}\sim\mathcal{N}(\frac{1}{K}\sum_{k}\mathcal{Y}_{k},\frac{1}{K} \sum_{k}(\sigma_{k}^{2}+\mathcal{Y}_{k}^{2})-\mathcal{Y}_{*}^{2})\).

## V Results and Discussion

We present and compare the model performances of the Prob-GNNs under six probabilistic and two deterministic assumptions on two separate datasets. First, we analyze the model performance for periods immediately after the training period, that is, from March 2 to March 15, 2020. We demonstrate that probabilistic assumptions can influence the model performance more significantly than deterministic architectures. Next, we show that uncertainty predictions are important especially when point predictions become unreliable during system disruptions, by applying models trained with the pre-COVID training set to the post-COVID test sets. Lastly,

Fig. 3: Proposed model architecture. The model processes past demand information using spatial and temporal layers, with the addition of auxiliary information to produce estimates of distributions parameters

Fig. 2: Average Daily Ridership (T1: Immediately Before; T2: Stay-at-home; T3: Initial Recovery; T4: Steady Recovery)
we discuss the spatiotemporal patterns of uncertainty revealed by the models.

First, to illustrate the models' proper convergence, Figure 4 presents the GCN model training curves for all probabilistic assumptions trained on CTA Rail. All models demonstrated successful convergence under the same learning rate. However, the number of epochs required for convergence varied among the models. The two-parameter distributions (HetG, TG, Laplace) converged faster, and achieved a lower NLL compared to the one-parameter distributions (HomoG, Pois).

Model performances of CTA Rail and ridesharing data, on the test set between March 2 and March 15, 2020, are tabulated in Tables III. The table has two panels for GCN and GAT models, respectively. Each panel presents the six different probabilistic assumptions presented in Section IV-B. Model performances are measured on three categories of metrics: composite (NLL), uncertainty (calibration error, MPIW, and PICP), and point (MAE, MAPE). In the tables, the intended level of coverage was set to 95%. The next two sections discuss the findings from Table III in detail.

### _Significance of Probabilistic and Deterministic Assumptions_

The quality of uncertainty quantification varies remarkably with probabilistic assumptions. The performance is significantly higher in the probabilistic models with two parameters (HetG and Lap) compared to the single-parameter models (HomoG and Pois). Among all GCN models, the average NLL and calibration error for two-parameter distributions in the CTA Rail data is \(465.9\) and \(0.033\), and for one-parameter distributions \(570.8(+22.5\%)\) and \(0.11(3.3\times)\). Similar observations can be made in the ridesharing dataset, where the average NLL is \(218.8\) vs. \(189.3(+15.4\%)\) and the calibration error is \(0.017\) vs. \(0.099(5.8\times)\). Similar observations can be made for GAT models.

Truncated Gaussian and Laplace distribution have the overall best distributional fit. Figure 5 illustrates the distributional fit by plotting the empirical quantiles \(q(p)\) against the theoretical quantiles \(p\) for both datasets across all the GCN models. The line \(y=x\) represents a perfectly-calibrated model. The calibration error in Table III corresponds to the areas between the \(y=x\) line and the empirical quantiles. Poisson and homoskedastic Gaussian are very far away from the \(y=x\) line, while truncated Gaussian and Laplace trace the line most closely. Although the Heteroskedastic Gaussian works the best with GCN models on ridesharing data, the truncated Gaussian and Laplace distribution had more stable performances across the datasets and models, and the curves are closer to the \(y=x\) line.

In contrast to the strong influence probabilistic assumptions have on the quality of uncertainty predictions, little influence is exerted on the point prediction quality from the probabilistic assumptions. For both datasets and both GCN and GAT, between different probabilistic assumptions, the performance gap between the best and the worst point estimate is around 4%, significantly less than that of the composite and the uncertainty metrics.

The variations in predictive performance caused by the deterministic assumptions are also small compared to the probabilistic assumptions. Comparing the GCNs and GATs, the error metrics are similar across the models, and the performance patterns across probabilistic assumptions remain the same for GAT models. The one-parameter probabilistic assumptions are likely to be affected more by deterministic architectures. For all probabilistic assumptions except for Poisson, less than a 3% difference in NLL loss is observed between GCN and GAT models for the same dataset. Although the GAT allows the model to learn spatial relationships and is hence more flexible, the predefined adjacency matrices serve as domain knowledge to reduce the learning complexity. In theory, the GAT setup should perform better with more complex relationships or with larger sample sizes, and the GCN setup is better for more efficient learning under limited sample sizes. In this case, there is not a large difference between the two architectures.

### _Implications of Probabilistic Assumptions_

The probabilistic assumptions not only suggest the distributional fit, but also have practical implications on the ridership patterns. We then discuss these implications in detail below.

#### Iv-B1 Heteroskedasticity vs. Homoskedasticity (HetG vs. HomoG)

HetG forces a constant variance across all observations, resulting in an inaccurate representation of the data. The HomoG models achieve around \(20\%\) higher NLL loss and \(1.5\) to \(2\) times higher calibration error than the HetG models. For example, the NLL loss of HetG-GCN is 462.9, 24% lower than

Fig. 4: Training Curves. Top: GCN; Bottom: GAT

Fig. 5: Calibration Plot of GCN Models (a) CTA rail and (b) Ridesharing. The line \(y=x\) represents a perfectly-calibrated model. The calibration error (CE) corresponds to the areas between the \(y=x\) line and the empirical quantiles.
the NLL loss (612.5) of HomoG-GCN for CTA Rail. Although predicting the same variance is inaccurate distributionally, the average MPIW and PICP are not much worse compared to their heteroskedastic counterparts since the variance scale is searched and fixed. Additionally, the inaccurate distributional assumption does not affect the quality of the point estimate, and point estimate errors from HomoG models are only 1-2% higher than the best model, as the likelihood improvement can only come from a more accurate prediction of the mean. Regardless, having one fixed uncertainty parameter not only is an inaccurate representation, but also limits the model's flexibility to adapt to sharp changes in ridership magnitude.

#### Iv-B2 Continuous vs. Discrete (HomoG vs. Pois)

Both HomoG and Pois have the worst NLL loss among all probabilistic assumptions tested. Despite its discreteness, the Poisson assumption yields similar or worse NLL and calibration error than HomoG since it forces equal mean and variance. However, for both datasets, the prediction intervals significantly under-predict the demand, meaning that the variance should be larger than the mean. The PICP of the Poisson prediction interval indicates the magnitude of variance compared to the mean. The smaller the PICP, the larger the variance.

#### Iv-B3 Real-line vs. non-negative support (HetG vs. TG)

Since the time resolution is 15min, left-truncation at 0 significantly improves the model performance as many predictions have means close to 0. TG is the best-performing model for the CTA Rail and the second-best for ridesharing. Even though the average ridership is around 40 per station and 23 per census tract, significant heterogeneity exists among different stations/census tracts. Sparsity (zero ridership) is an issue to be considered in short-term or real-time demand predictions.

#### Iv-B4 Gaussian vs. Exponential tails (HetG vs. Lap)

Both distributions are characterized by two parameters and can be trained to describe the data accurately. The weights of the tails are different, and the model performances on the two probabilistic assumptions suggest behavioral differences. Compared to HetG, Lap has a heavier tail; hence the prediction intervals tend to be larger, covering more observations in more extreme cases. In all cases, Lap covers more points than intended.

#### Iv-B5 Single vs. ensembled models (HetG vs. GEns)

Ensembling only makes marginal improvements to the model results. The NLL loss, calibration error, and MAE improve by less than 1% in all cases. For example, the biggest improvement occurs in the GAT models for CTA Rail data, with NLL loss for HetG-GAT being 462.9 and for GEns-GAT 459.6. Since ensembling aims at reducing model uncertainty, its ineffectiveness suggests that different training instances of the neural networks produce similar results, and the model uncertainty is low.

We also constructed prediction intervals with Monte Carlo dropouts to further illustrate that model uncertainty is much smaller than the data uncertainty inherent in the data generation process. Monte Carlo dropout approximates Bayesian neural networks and measures the model uncertainty by applying dropout at test time, treating each as an instance from the space of all available models. However, since the different training instances produce similar results, Monte Carlo dropout fails to capture the full picture. Its prediction intervals were very narrow with PICPs between 30% - 40% for both datasets. Additionally, since we applied dropout at test time, the point prediction loses the benefit of dropout regularization and performs worse than other models.

### _Model Performance under System Disruption_

The drastic change in pre- and post-COVID periods presents a unique opportunity to test the generalizability of Prob-GNNs. The previous sections use the periods immediately following the training set for validation and testing. This section compares the predictive performance at different stages of the COVID-19 pandemic, by applying the models trained with the pre-COVID training set to three post-COVID time periods. Table IV presents the results for CTA Rail and ridesharing. Since HomoG and Pois have poor performance in the previous test set, they are excluded from this comparison.

All the error metrics increase under the significant domain shifts from pre- to post-COVID. In the stay-at-home period, the average ridership for both systems dropped to less than 10% of pre-COVID levels. The NLL, MPIW, and MAE are not comparable to pre-COVID levels because they are influenced

\begin{table}
\begin{tabular}{l|c|c|c|c|c||c|c|c|c|c|c} \hline \hline  & \multicolumn{6}{c||}{CTA Rail} & \multicolumn{6}{c}{Ridesharing} \\ \hline  & Comp & \multicolumn{3}{c|}{Uncertainty Prediction} & \multicolumn{3}{c|}{Point Prediction} & \multicolumn{3}{c|}{Comp} & \multicolumn{3}{c|}{Uncertainty Prediction} & \multicolumn{3}{c}{Point Prediction} \\ \hline Model & NLL & Cal. Err & MPIW & PICP & MAE & MAPE & NLL & Cal. Err & MPIW & PICP & MAE & MAPE \\ \hline HomoG-GCN & 612.5 & 0.146 & 64.53 & 97.9\% & 7.64 & 18.5\% & 225.6 & 0.114 & 42.63 & 97.7\% & 5.81 & 25.2\% \\ Pois-GCN & 529.2 & 0.075 & 18.04 & 84.9\% & **7.18** & **17.4\%** & 211.9 & 0.083 & 15.53 & 82.7\% & 5.81 & 25.2\% \\ HetG-GCN & 462.9 & 0.049 & 39.04 & 96.5\% & 7.67 & 18.6\% & 187.6 & 0.019 & 27.32 & 93.3\% & 5.80 & 25.1\% \\ TG-GCN & 480.7 & **0.021** & 39.37 & 95.6\% & 7.61 & 18.5\% & 196.7 & 0.029 & **36.63** & **95.2\%** & 6.06 & 26.3\% \\ Lap-GCN & 460.2 & 0.026 & 43.48 & 97.0\% & 7.52 & 18.2\% & 187.1 & 0.022 & 34.12 & 96.5\% & 5.81 & 25.2\% \\ GENs-GCN & **459.6** & 0.036 & **37.48** & **95.7\%** & 7.42 & 18.0\% & **185.6** & **0.017** & 28.01 & 93.9\% & **5.70** & **24.7\%** \\ \hline HomoG-GAT & 620.4 & 0.141 & 64.54 & 97.8\% & 8.13 & 19.7\% & 230.7 & 0.112 & 43.26 & 96.9\% & 6.52 & 28.2\% \\ Pois-GAT & 572.2 & 0.100 & 18.11 & 81.8\% & 7.98 & 19.4\% & 214.9 & 0.093 & 15.62 & 82.0\% & **6.06** & **26.3\%** \\ HetG-GAT & 472.3 & 0.058 & **39.58** & **96.3\%** & 7.84 & 19.0\% & 188.3 & 0.041 & **31.43** & **95.75** & 6.07 & 26.3\% \\ TG-GAT & 481.6 & **0.026** & 40.96 & 96.6\% & 7.69 & 18.7\% & 199.9 & 0.024 & 32.54 & 95.7\% & 6.55 & 28.4\% \\ Lap-GAT & **464.3** & 0.037 & 42.77 & 97.7\% & 7.63 & 18.5\% & 189.5 & **0.021** & 36.08 & 96.5\% & 6.16 & 26.7\% \\ GENs-GAT & 470.5 & 0.037 & 42.49 & 96.8\% & **7.52** & **18.3\%** & **188.1** & 0.044 & 32.28 & 96.0\% & 6.18 & 26.8\% \\ \hline \hline \end{tabular}
\end{table} TABLE III: Model Performance (Test Period: Immediately Before - March 2 to March 15, 2020)
by the magnitude of ridership, while calibration error, PICP, and MAPE are unit-free and can be compared to values in Table III. Unsurprisingly, the performance during the stay-at-home period is relatively poor but slowly rebounds with the recovery of the ridership.

In the case where significant disruptions happen in the system, the point predictions fail miserably, but the uncertainty predictions stay accurate and indicative of the changing situation. The MAPE for the test set in Table III was 18% for CTA and 25% for ridesharing. For the three additional periods, the MAPEs are 93%, 43%, 36% for CTA and 133%, 54%, and 47% for ridesharing, respectively. The uncertainty predictions recovered a lot faster than MAPE. The calibration error returned to pre-COVID levels at the steady recovery stage, although the point prediction error is still 10% higher than before. If only a 95% prediction interval is considered, even at the stay-at-home home stage, we can achieve pretty precise prediction intervals (95.6% and 96.7% for CTA and ridesharing).

The generalizability between different two-parameter distributions is similar, although each has distinct characteristics. TG-GCN excels at low ridership; Lap-GCN is heavy-tailed and more conservative; while HetG-GCN relies on the similarity between traning and testing domains. Since the situation is evolving constantly, there is no probabilistic assumption that dominates all scenarios. But some general conclusions can be drawn. In most cases, TG-GCN has the best calibration error due to its left-truncation. Enforcing non-negativity is beneficial since the ridership haven't recovered to pre-COVID levels for either CTA or ridesharing. Lap-GCN typically produces the best NLL and point prediction due to its more distributed shape, reducing overfitting to pre-COVID data. As ridesharing trips become more spontaneous in post-COVID times, Lap-GCN's conservative prediction intervals outperform others.

### _Spatiotemporal Uncertainty_

Travel demand uncertainty has both spatial and temporal patterns: spatially, uncertainty is higher at stations with higher volumes, and temporally, uncertainty is higher in the afternoons. Figure 6 shows the spatial distribution of predicted uncertainty for CTA Rail at different times of the day during the steady recovery test period (Oct 12 - 25, 2020). The left bottom corner of each figure zooms in on the "loop", the downtown of Chicago. Uncertainty is proportional to both the size and the color of the circles, with darker and larger circles indicating higher uncertainty.

Spatially, uncertainty is higher at busier stations. During the same time period, darker and larger circles appear near downtown, transfer stations, and airports. Statistically, if the occurrence of each potential trip is a random Bernoulli process with probability \(p\), and we do not have further information on the values of \(p\) for each trip, having more potential trips \(n\) will yield higher uncertainty in the sum. The number of observed trips has a binomial distribution with variance \(np(1-p)\), which is proportional to \(n\). Practically, the trips from downtown, transfer stations, and airports are usually more diverse and complex in nature, which promotes spontaneity in trip-making, resulting in higher uncertainty.

Temporally, uncertainty is higher during the afternoon peak. The uncertainty during the morning peak is even generally smaller than that during midday and in the evening, although the morning peak has significantly higher ridership. Statistically, this observation could be attributed to having knowledge about some of the \(p\)'s. Since the morning peak primarily consists of commuting trips, the probability of the trips happening is higher than recreational trips, which tend to happen during midday and afternoons.

Spatiotemporal uncertainty predictions can inform strategic decision-making regarding capacity buffers. First, the framework could be used to identify bottlenecks and outliers in the system. For example, the station Division/Milwaukee on the
blue line has unusually high uncertainty during the morning peak. Further investigation can be done to identify the reasons for the abnormal behavior and perform demand management. Moreover, different strategies are needed for different types of systems. In systems with a fixed, relatively large capacity, such as stations along the same subway line, uncertainty at the busiest stations at peak times is the most critical, as the rest of the line will have quite a lot of excess capacity. Therefore, understanding uncertainty in low-demand regions is important for behavioral analysis, but less critical for service planning. However, in systems that are built to meet demand, such as ride-hailing, uncertainty in lower-demand regions is as important as higher-demand regions, as the supply in those regions will be proportionally low. Re-balancing actions will be needed across the system.

## VI Conclusion

Despite the importance of uncertainty in travel demand prediction, past studies use deep learning to predict only the average travel demand but not quantify its uncertainty. To address this gap, this study proposes a framework of Prob-GNNs to quantify the spatiotemporal uncertainty of travel demand. The framework is concretized by six probabilistic assumptions (HomoG, Pois, HetG, TG, GEns, Lap) and two deterministic ones (GCN and GAT), which are applied to transit and ridesharing data, yielding the following conclusions.

First, the Prob-GNN framework can successfully quantify spatiotemporal uncertainty in travel demand with empirical evidence demonstrated on the transit and ridesharing datasets. In both cases, the Prob-GNNs can accurately characterize the probabilistic distributions of travel demand, while retain the point predictions of a similar quality to the deterministic counterpart (HomoG). Second, the probabilistic assumptions have a much more substantial impact on model performance than the deterministic ones. The two deterministic architectures lead to less than a 3% difference in NLL loss, while a wise choice of probabilistic assumptions could drastically improve model performance. Specifically, the two-parameter distributions (e.g., truncated Gaussian) achieve the highest predictive performance, which is 20% higher in NLL and 3-5 times lower in calibration errors in comparison to the one-parameter baseline distributions. Third, the Prob-GNNs enhance model generalizability under significant system disruptions. By applying the models trained on pre-COVID data to three post-COVID periods, we show that the point predictions fail to generalize but the uncertainty predictions remain accurate and successfully reflect the evolving situations. Even under significant domain shifts, the difference of predictive performance among the two-parameter distributions is minor. Lastly, Prob-GNNs can reveal spatiotemporal uncertainty patterns. Uncertainty is spatially concentrated on the stations with higher travel volume, and temporally concentrated on the afternoon peak hours. In addition, the Prob-GNNs can identify the stations with abnormally large uncertainty, which can inform real-time traffic controls to proactively address potential system disruptions.

Future studies could advance this work by taking further theoretical or empirical efforts. Theoretically, Prob-GNN is a parametric uncertainty quantification method, which should be compared to Bayesian and non-parametric methods in terms of prediction quality. Empirically, our framework can be generally applied to the prediction of origin-destination flows, travel safety, or even climate risks. Since uncertainty in urban systems has broad policy implications, future studies could integrate the Prob-GNNs with robust optimization methods to enhance urban resilience.

## Acknowledgment

This material is based upon work supported by the U.S. Department of Energy's Office of Energy Efficiency and Renewable Energy (EERE) under the Vehicle Technology Program Award Number DE-EE0009211. The views expressed herein do not necessarily represent the views of the U.S. Department of Energy or the United States Government. The

Fig. 6: Spatiotemporal Uncertainty. Standard deviations of estimated CTA Rail station tap-ins in the 15min periods starting at 8 A.M. (morning peak), 1 P.M. (middle), 5 P.M. (afternoon peak), and 9 P.M. (evening).
authors would also like to thank the Chicago Transit Authority for providing access to public transit ridership data.
