Ambulance Demand Prediction via Convolutional Neural Networks
Minimizing response times is crucial for emergency medical services to reduce
patients' waiting times and to increase their survival rates. Many models exist
to optimize operational tasks such as ambulance allocation and dispatching.
Including accurate demand forecasts in such models can improve operational
decision-making. Against this background, we present a novel convolutional
neural network (CNN) architecture that transforms time series data into
heatmaps to predict ambulance demand. Applying such predictions requires
incorporating external features that influence ambulance demands. We contribute
to the existing literature by providing a flexible, generic CNN architecture,
allowing for the inclusion of external features with varying dimensions.
Additionally, we provide a feature selection and hyperparameter optimization
framework utilizing Bayesian optimization. We integrate historical ambulance
demand and external information such as weather, events, holidays, and time. To
show the superiority of the developed CNN architecture over existing
approaches, we conduct a case study for Seattle's 911 call data and include
external information. We show that the developed CNN architecture outperforms
existing state-of-the-art methods and industry practice by more than 9%.
Maximiliane RautenstrauÃ, Maximilian Schiffer
# Ambulance Demand Prediction via Convolutional Neural Networks

###### Abstract

Minimizing response times is crucial for emergency medical services to reduce patients' waiting times and to increase their survival rates. Many models exist to optimize operational tasks such as ambulance allocation and dispatching. Including accurate demand forecasts in such models can improve operational decision-making. Against this background, we present a novel convolutional neural network (CNN) architecture that transforms time series data into heatmaps to predict ambulance demand. Applying such predictions requires incorporating external features that influence ambulance demands. We contribute to the existing literature by providing a flexible, generic CNN architecture, allowing for the inclusion of external features with varying dimensions. Additionally, we provide a feature selection and hyperparameter optimization framework utilizing Bayesian optimization. We integrate historical ambulance demand and external information such as weather, events, holidays, and time. To show the superiority of the developed CNN architecture over existing approaches, we conduct a case study for Seattle's 911 call data and include external information. We show that the developed CNN architecture outperforms existing state-of-the-art methods and industry practice by more than 9%.

**Keywords:** AI in health care, spatio-temporal forecasting, convolutional neural networks
## 1 Introduction

Reducing response times is paramount for emergency medical services (EMS) to provide first aid in a timely manner. However, tight budgets pressure EMS systems to minimize operational expenses, resulting in limited availability of resources. Demographic changes such as an aging population can further increase ambulance demands, intensifying the need for an efficient use of resources. Climate change causing more frequent extreme weather conditions, such as heat spells or heavy rainfalls, can additionally challenge EMS systems in the future. To tackle these challenges and to reduce response times, several models exist to optimize operational tasks such as ambulance allocation and dispatching, see, e.g., Brotcorne et al. (2003), Aboueljinane et al. (2013), Belanger et al. (2019), and Farahani et al. (2019) for literature overviews. Embedding accurate demand forecasts in such models can significantly improve operational decision-making.

In literature, various machine learning approaches such as artificial neural networks (ANNs), tree-based models, or support vector regressions (SVRs) exist to predict ambulance demand. For other prediction tasks, such as traffic data and mobility demand predictions, convolutional neural networks (CNNs) yield highly accurate predictions, see, e.g., Wang et al. (2018) and Guo et al. (2019). Although CNNs can outperform state-of-the-art approaches for such prediction tasks, CNNs have not yet been applied for ambulance demand prediction. Against this background, we present a novel CNN architecture that transforms time series information into heatmaps to forecast ambulance demand.

Our contribution to the existing literature is three-fold: First, we present a flexible, generic CNN architecture, allowing for the inclusion of external features with varying dimensions. This concept differs from most state-of-the-art approaches that divide the given region into subregions and derive local forecasts for each of them separately, neglecting possible spatial correlations. Contrarily, we enable the detection of correlations in space and time by including three-dimensional convolutional layers in a CNN architecture and obtain the predictions for all subregions simultaneously. Most CNN architectures neglect external information as they have been mainly applied to image classification and recognition tasks. Nonetheless, effectively predicting ambulance demand with CNN architectures requires the incorporation of external features. Thus, we integrate historical ambulance demand and external information such as weather, events, holidays, time, weekdays, and months. Second, we show how to jointly perform feature selection and hyperparameter tuning by applying Bayesian optimization (BO). We treat the decision of whether to include a feature or not as an additional hyperparameter of the prediction model and add these parameters to the hyperparameter tuning process. To tackle the high-dimensional search space, we introduce a novel hierarchical BO approach and apply BO with dimension dropout. Third, we analyze the features' importance by calculating the SHapley Additive exPlanation (SHAP) values representing the contribution of each feature to the prediction. Additionally, we explore the CNN's performance for different forecasting horizons by conducting a sensitivity analysis.

We base our results on a case study of Seattle's 911 call data and include external information. Results show that the developed CNN outperforms state-of-the art methods and industry practice
by more than 9%. Incorporating feature selection in the BO reduces the number of parameters by 40.4% while the performance decreases by less than 0.05%. The SHAP values show that the CNN is more robust to changes in the historical ambulance demand, compared to its benchmarks. This enables a better handling of demand outliers or data inconsistencies. We further show that the BO with dimension dropout and the hierarchical BO outperform basic BO approaches and random search.

The remainder of this paper is structured as follows. We first present related literature in Section 2. Section 3 describes our problem setting. Section 4 introduces the CNN architecture and our algorithmic framework to conduct feature selection and hyperparameter tuning. In Section 5, we introduce the numerical case study for Seattle's 911 call data and present our results in Section 6. Section 7 concludes this work.

## 2 Related work

Many early but also recent studies apply regression models to predict ambulance demand, investigating the influence of socio-economic variables, e.g., population, age, gender, race, land use, living conditions, education, employment, income, or marital status, see, e.g., Aldrich et al. (1971), Siler (1975), Lowthian et al. (2011), and Steins et al. (2019). Wong & Lin (2020) apply multivariate forward regression and show that EMS demand for specific patient groups, such as elderly or critical patients, is sensitive to weather conditions and changes throughout the day. Another stream of literature focuses on time series models to predict ambulance demand. Baker & Fitzpatrick (1986) apply Winter's exponential smoothing (Winters, 1960) and present a goal programming approach to optimize the weightings given to different forecasting statistics. Tandberg et al. (1998) implement different time series models to predict hourly emergency incidents in Albuquerque and show that simple time series models can outperform more expensive, complex models. Channouf et al. (2007) include different day-of-week, month-of-year, and special-day effects in their time series models. They show that considering historic call volumes of previous hours can improve the hourly forecast accuracy. Matteson et al. (2011) make predictions on an hourly level and introduce a factor model combined with time series models. They include time information such as the weekday and week. Vile et al. (2012) present a singular spectrum analysis to forecast emergency calls in Wales and show that their approach outperforms Holt-Winters and autoregressive integrated moving-average (ARIMA) models for long-term forecasts. Wong & Lai (2014) integrate weather data in an ARIMA model and show that the prediction error of the 7-day forecast can be reduced by 10% by incorporating weather data in the model. Gijo & Balakrishna (2016) introduce a Seasonal ARIMA model to predict the emergency call volume for a state in India. Nicoletta et al. (2016) introduce a Bayesian approach, modeling ambulance demand as a generalized linear mixed model. Results show that the time-of-day, holidays, each zone's population, area, and type heavily influence the prediction's performance. Zhou et al. (2015) present a Gaussian mixture model to predict ambulance demand in Toronto, Canada, tackling data sparsity caused by fine temporal resolutions. Zhou & Matte
son (2016) apply kernel density estimation and introduce a kernel warping approach, a form of Laplacian eigenmaps belonging to the field of manifold learning. Setzler et al. (2009) present a multilayer perceptron (MLP) including four categorical variables: time-of-day, day-of-week, month, and season. Chen et al. (2015) apply an ANN, moving average, sinu-soidal regression, and support vector regression model for ambulance demand prediction. For each subregion, they select the best-performing model. Wang et al. (2021) include a heterogeneous multi-graph convolutional layer in a neural network to predict ambulance demand. To form the graph, they make use of dispatching areas. Jin et al. (2021) introduce a bipartite graph convolutional network to predict ambulance demand for region-hospital pairs. They include regional features and hospital features, e.g., capacity information.

In summary, many forecasting approaches exist to predict ambulance demand. These approaches include regression analyses, time series models, mixture models, ANNs, graph-based models, SVRs, and tree-based models. We provide a summary including the features considered in each study in Table 1 and refer to Reuter-Oppermann et al. (2017) for an overview. In the field of mobility demand forecasting, few approaches exist that base forecasting on CNNs, which have so far mostly been used for image classification, to implicitly capture spatio-temporal correlations (Wang et al., 2018; Guo et al., 2019). Although CNNs yield highly accurate predictions for such applications, they have not yet been applied to predict ambulance demand. Moreover, a framework that allows generic feature integration, efficient feature selection, and hyperparameter tuning for such architectures has neither been developed for general mobility demand nor for ambulance demand forecasting. To close this research gap, we introduce a novel CNN architecture for predicting ambulance demand that is able to generically incorporate external features and embed it in an efficient hyperparameter tuning framework.

## 3 Problem Setting

We divide the region for which we are making the forecast into \(q\times p\) subregions and aim at predicting the ambulance demand for all subregions for a time period \(t\). For each subregion, the ambulance demand can be represented as discrete time series \(z\in\mathbb{R}^{T}\).

**Definition 3.1** (Time series).: _Let \(z\in\mathbb{R}^{T}\) be a discrete time series containing a sequence of \(T\) observations sampled from a random variable \(\mathcal{Z}\) at equidistant timesteps._

Based on each time series, conventional algorithms, e.g., exponential smoothing (Baker & Fitzpatrick, 1986) or ARIMA models (Tandberg et al., 1998), can predict the future ambulance demand for every subregion in an iterative fashion. However, such iterative approaches neglect spatial correlations between the time series. For this reason, we aim at generating an integrated forecast for all subregions simultaneously. To enable an integrated forecast, we represent the ambulance demand as a spatio-temporal time series assuming a spatial correlation between the subregions' ambulance demands.
\begin{table}
\begin{tabular}{c
**Definition 3.2** (Spatio-temporal time series).: _Assuming a spatial correlation between the subregions' time series, we can represent the ambulance demand as a \(q\times p\) matrix of spatially correlated time series, each of length T, which we denote by \(Z\in\mathbb{R}^{q\times p\times T}\)._

In addition, effectively predicting ambulance demand requires the incorporation of external features. In the existing literature, various external features have been considered such as socio-economic or weather data. These external features strongly differ in their properties. For example, socio-economic data mostly differs across city districts but only changes gradually over time. In contrast, weather data can change quickly but does not strongly differ across city districts. To generically include external features, we distinguish four types:

**Three-dimensional input data.**: Let \(X^{\rm 3D}=\{X^{\rm 3D}_{i}\}_{i=1}^{n}\) be the three-dimensional features of our data set, where each data instance \(X^{\rm 3D}_{i}\in\mathbb{R}^{q\times p\times L}\) is a spatio-temporal time series representing the historic ambulance demand of \(L\) periods for all \(q\times p\) subregions.
**Two-dimensional input data.**: Let \(X^{\rm 2D}=\{X^{\rm 2D}_{i}\}_{i=1}^{n}\) be the two-dimensional features of our data set, where each data instance \(X^{\rm 2D}_{i}=\{x^{\rm 2D}_{i_{j}}\}_{j=1}^{m}\) is a set of \(m\) two-dimensional vectors, i.e., \(x^{\rm 2D}_{i_{j}}\in\mathbb{R}^{q\times p}\). We include two-dimensional data in the case that only the spatial dimension is considered, e.g., for representing events taking place in the predicted period.
**One-dimensional input data.**: Let \(X^{\rm 1D}=\{X^{\rm 1D}_{i}\}_{i=1}^{n}\) be the one-dimensional features of our data set, where each data instance \(X^{\rm 1D}_{i}=\{X^{\rm 1D}_{i_{j}}\}_{j=1}^{r}\) is a set of \(r\) one-dimensional vectors, i.e., \(x^{\rm 1D}_{i_{j}}\in\mathbb{R}^{l_{j}}\), where \(l_{j}\) is the length of vector \(j\). One-dimensional features are either time series such that \(l_{j}=L\), where \(L\) is the number of past periods to be considered, or one-hot-encoded information, e.g., for including time information such as the month, weekday, or hour.
**Scalar input data.**: Let \(X^{\rm S}=\{X^{\rm S}_{i}\}_{i=1}^{n}\) be the scalar features of our data set, where each data instance \(X^{\rm S}_{i}=\{x^{S}_{i_{j}}\}_{j=1}^{s}\) is a set of \(s\) scalar features, i.e., \(x^{\rm S}_{i_{j}}\in\mathbb{R}\). We include scalar inputs for data without spatial or temporal distribution which is not one-hot-encoded, e.g., the temperature prediction for period \(t\) for which we neglect the spatial distribution.

Let \((\mathbf{X},\mathbf{y})=\{\mathbf{x}_{i},y_{i}\}_{i=1}^{n}\) be our dataset. \(\mathbf{x}_{i}=\{X^{\rm 1D}_{i},X^{\rm 2D}_{i},X^{\rm 3D}_{i},X^{\rm S}_{i}\} \in\mathbb{R}^{v}\) represents a \(v\)-dimensional feature vector and \(y_{i}\in\mathbb{R}\) its respective dependent variable. Each data instance is independently sampled from an unknown distribution \((\mathcal{X},\mathcal{Y})\). Our model \(\mathcal{M}\) aims at learning a function \(F_{\mathcal{M}}:\mathcal{X}\rightarrow\mathcal{Y}\) that maps each input vector \(x_{i}\) to its associated variable \(y_{i}\).

## 4 Methodology

We first introduce our CNN architecture in Section 4.1. Then, we describe our BO approaches to perform hyperparameter optimization and feature selection in Section 4.2.
### Convolutional Neural Network Architecture

To predict ambulance demand, we apply a novel CNN architecture. Leveraging our neural network with convolutional layers enables the generation of feature maps, i.e., transformed images, highlighting patterns identified in the input image. To form these feature maps, filters, i.e., weight matrices, traverse the input image with a defined stride, i.e., step size, and conduct convolutions to detect patterns that the filters have learned to identify. So far, CNNs have mainly been used to solve image classification tasks, in which convolutional layers detect spatial patterns in the images. We apply this concept to reveal spatio-temporal patterns within our historic ambulance demand. We further integrate external features through concatenation at different stages of the CNN. We visualize our CNN architecture in Figure 1 and provide a more detailed overview in Appendix A.1.

We apply three-dimensional convolutional layers to detect spatio-temporal patterns in the historic ambulance demand, represented as spatio-temporal time series \(X^{3D}\in\mathbb{R}^{q\times p\times L}\), where \(L\) is the number of past periods that we take into account. Each convolutional layer \(l\) contains \(M_{l}\) filters. The weight of filter \(h\) in layer \(l\) at position \((w_{1},w_{2},w_{3})\) is \(w_{h}^{l}(w_{1},w_{2},w_{3})\). We define the filter size by \(k^{l}=(k_{1}^{l},k_{2}^{l},k_{3}^{l})\) and apply a stride of \(s^{l}=(s_{1}^{l},s_{2}^{l},s_{3}^{l})\). We denote the activation function in layer \(l\) by \(\vartheta^{l}(\cdot)\) and the bias by \(b_{h}^{l}\). By applying zero padding, i.e., surrounding the input image with additional zero-valued pixels, we maintain the output's size throughout consecutive layers. The coordinates \((x_{1},x_{2},x_{3})\) uniquely define the values within each feature map, assuming zero-based indexing. We calculate the number of parameters \(\phi_{l}\) learned for layer \(l\) by

\[\phi_{l}=(k_{1}^{l}k_{2}^{l}k_{3}^{l}M_{l-1}+1)M_{l}, \tag{4.1}\]

and derive the output of filter \(h\) in layer \(l\) by

Figure 1: Visualization: CNN Architecture
\[\begin{split} o_{h}^{l}\left(x_{1},x_{2},x_{3}\right)=\vartheta^{l} \Bigg{(}\sum_{m=0}^{M_{l-1}-1}&\sum_{\begin{subarray}{c}i_{1}=0, \\ \iota(x_{1},i_{1},s_{1}^{l})\end{subarray}}^{\left\lfloor\frac{x_{1}}{s_{1}^{l} }\right\rfloor}&\sum_{\begin{subarray}{c}i_{2}=0,\\ \iota(x_{2},i_{2},s_{2}^{l})\end{subarray}}^{\left\lfloor\frac{x_{2}}{s_{2}^{ l}}\right\rfloor}&\sum_{\begin{subarray}{c}i_{3}=0,\\ \iota(x_{3},i_{3},s_{3}^{l})\end{subarray}}^{\left\lfloor\frac{x_{3}}{s_{3}^{l} }\right\rfloor}& o_{m}^{l-1}\left(\left\lfloor\frac{x_{1}}{s_{1}^{l} }\right\rfloor-i_{1},\left\lfloor\frac{x_{2}}{s_{2}^{l}}\right\rfloor-i_{2}, \left\lfloor\frac{x_{3}}{s_{3}^{l}}\right\rfloor-i_{3}\right)\\ &\ast w_{h}^{l}\left(\iota\left(x_{1},i_{1},s_{1}\right)\iota \left(x_{2},i_{2},s_{2}\right)\iota\left(x_{3},i_{3},s_{3}\right)\right)+b_{h }^{l}\Bigg{)},\end{split} \tag{4.3}\]

where

\[\iota\left(x,i,s\right)=x-s\left(\left\lfloor\frac{x}{s}\right\rfloor-i\right). \tag{4.4}\]

After conducting the convolutions (4.2) and (4.3), we concatenate the generated feature maps. To enable this concatenation, we set the stride, padding, and filter sizes such that all feature maps are of shape \(q\times p\times L\). We then fuse the concatenated feature maps along the temporal dimension. We concatenate the temporal fusion's outputs, i.e., two-dimensional feature maps of shape \(q\times p\), with the two-dimensional input data. In the following, we add a locally connected layer to combine the identified patterns with the included inputs. Here, we substitute \(w_{h}^{l}\) by \(w_{hx_{1}x_{2}}^{l}\) and \(b_{h}^{l}\) by \(b_{hx_{1}x_{2}}^{l}\) to represent the weights and bias in filter \(h\) at position \((x_{1},x_{2})\). We compute the outputs as follows

\[\begin{split} o_{h}^{l}\left(x_{1},x_{2}\right)=\vartheta^{l} \Bigg{(}\sum_{m=0}^{M_{l-1}-1}&\sum_{w_{1}=0}^{k_{1}-1}\sum_{w_ {2}=0}^{k_{2}-1}o_{m}^{l-1}\left(x_{1}s_{1}^{l}+w_{1},x_{2}s_{2}^{l}+w_{2} \right)\\ &\ast w_{hx_{1}x_{2}}^{l}\left(w_{1},w_{2}\right)+b_{hx_{1}x_{2} }^{l}\Bigg{)}.\end{split} \tag{4.5}\]
For the locally connected layer \(l\), we train \(\phi_{l}\) parameters calculated by

\[\phi_{l}=(M_{l-1}+1)qp, \tag{4.6}\]

where \(qp\) is the number of output neurons.

We then embed one-dimensional data and scalar inputs through concatenation. Finally, we add two dense layers to enable the network to learn from these inputs. Each dense layer \(l\) consists of \(N^{l}\) neurons, each denoted by \(v_{i}^{l}\) where \(i=1,...,N^{l}\). The activation function and bias for layer \(l\) are \(\vartheta^{l}\left(\cdot\right)\) and \(b^{l}\). We denote the weights between neurons \(v_{i}^{l}\) and \(v_{j}^{l-1}\) by \(w_{ij}^{l}\). Thus, we calculate the output \(a_{i}^{l}\) of neuron \(v_{i}^{l}\) as follows

\[a_{i}^{l}=\vartheta^{l}\Big{(}\sum_{j=1}^{N^{l-1}}w_{ij}^{l}a_{j}^{l-1}+b^{l} \Big{)}. \tag{4.7}\]

We calculate the number of learned parameters for dense layer \(l\) by

\[\phi_{l}=(N^{l-1}+1)N^{l}. \tag{4.8}\]

In the last layer, we set the number of neurons to \(qp\) such that we can reshape its outputs into a \(q\times p\) heatmap, constituting the model's prediction.

### Hyperparameter Optimization and Feature Selection

To provide a lean model and to improve the model's performance, we conduct feature selection and tune the model's hyperparameters. Adding or eliminating features may change the network architecture, as some input layers are only added when including certain external information. Thus, optimizing hyperparameters before conducting feature selection may lead to a suboptimal choice of parameters, e.g., the number of filters and their sizes may depend on the features selected. Similarly, conducting feature selection before hyperparameter tuning may lead to suboptimal results, since the selection process may be based on a network with insufficiently tuned hyperparameters. For this reason, we aim at utilizing an intrinsic feature selection method tuning the CNN's hyperparameters and conducting feature selection simultaneously.

Therefore, we treat the decision of whether to include a feature or not as an additional hyperparameter of the prediction model and include these parameters in the hyperparameter tuning process. Pursuing such an approach entails the following challenges: First, the model has numerical and categorical hyperparameters that may depend on each other, e.g., adding layers requires their hyperparameters to be tuned. Vice versa, when removing layers, we can neglect their hyperparameters in the tuning process. Second, conducting feature selection for \(N\) features increases the search space by \(2^{N}\), resulting in a high-dimensional search space. Third, evaluating the predictive model is computationally expensive. Accordingly, common hyperparameter optimization strategies, such as manual or grid search, are unsuitable for adequately covering the search space in an acceptable amount of computational time. Fourth, the stochastic nature of the training process entails noise, as similar hyperparameter values can result in different function values. Finally, derivatives of the
underlying target function are not given, and multiple local optima may exist due to the target function's non-convexity. To tackle these challenges and to learn from information gained by previously performed parameter combinations, we apply BO.

Bayesian OptimizationBO is a sequential model-based optimization (SMBO) approach that is used to optimize the hyperparameters of a black-box function \(f\left(\cdot\right)\), which is expensive to evaluate. We leverage a surrogate model to approximate the target function to bypass its expensive evaluation. In each iteration of the BO, we update the surrogate model based on our previous observations. Then, we optimize an acquisition function utilizing the surrogate model, which derives the next promising hyperparameter settings to be evaluated.

We denote \(J\) hyperparameters by \(\theta=\lambda_{1},\lambda_{2},\lambda_{3},...,\lambda_{J}\) and the search space by \(\Theta=\Lambda_{1}\times\Lambda_{2}\times\Lambda_{3}\times...\times\Lambda_{J}\). Then, we aim to solve the following optimization problem:

\[\theta^{*}=\arg\min_{\theta\in\Theta}f\left(\theta\right). \tag{4.9}\]

As the target function \(f\left(\theta\right)\) is computationally expensive to evaluate, we model this function via the surrogate model \(\mathcal{M}\). As surrogates, we apply Gaussian processes (Snoek et al., 2012), random forests (Hutter et al., 2011), and extremely randomized trees (Geurts et al., 2006). In each iteration \(\iota\) of the BO, we optimize the acquisition function \(\mathcal{A}\) to derive a new promising parameter setting as follows

\[\theta_{\iota}\leftarrow\arg\max_{\theta\in\Theta}\mathcal{A}( \theta). \tag{4.10}\]

Here, we apply an _Expected Improvement_ strategy enabling a trade-off between exploration and exploitation: We define the best (minimal) observed function value of \(n\) iterations by \(f^{*}=\min(f(\theta_{1}),...,f(\theta_{n}))\). When sampling from \(f\) at a new, unknown point, we treat its realization \(Y\) as a normally distributed variable \(Y\sim\mathcal{N}(\mu,\sigma^{2})\). Thus, the expected improvement at point \(\theta\) is

\[\mathbb{E}[\mathbf{I}(\theta)]=\mathbb{E}[\max(f^{*}-Y,0)]=(f^{*}- \mu)\Phi\left(\frac{f^{*}-\mu}{\sigma}\right)+\sigma\phi\left(\frac{f^{*}-\mu }{\sigma}\right), \tag{4.11}\]

where \(\phi(\cdot)\) and \(\Phi(\cdot)\) are the standard normal density and distribution functions, correspondingly (Jones et al., 1998).

We present a basic pseudo-code for BO in Algorithm 1. After initializing the set of observations \(R=\{(\theta_{1},f(\theta_{1})),(\theta_{2},f(\theta_{2})),...,(\theta_{m},f( \theta_{m}))\}\) via \(m\) iterations of random search (l. 3), we conduct \(n\) iterations of BO. In each iteration \(\iota\), we update the surrogate model and optimize the acquisition function \(\mathcal{A}\) (l. 6). We then evaluate the target function \(f\) at the next promising point \(\theta_{\iota}\) and update the set of results \(R\) with our new observation (l. 7). Finally, we return the incumbent parameter settings \(\theta^{*}\) (l. 9).

High-Dimensional Bayesian OptimizationHigh-dimensional search spaces challenge existing BO methods. Although classic approaches are still applied, e.g., random search (Bergstra & Bengio, 2012) or BO with random forests (Hutter et al., 2011), more advanced methods have been developed
to tackle high-dimensional search spaces. For example, Wang et al. (2013) make use of low effective dimensions, i.e., dimensions which do not significantly change the objective function, to reduce the search space. However, assuming that only a subset of dimensions is effective may not be possible, e.g., if all dimensions highly influence the objective. To handle problems with high effective dimensionality, Li et al. (2018) introduce the concept of random dimension dropout, which optimizes the acquisition function only over a subset of dimensions. In this study, we implement two different approaches. First, we build upon Li et al. (2018) and apply BO with dimension dropout. Second, we present a novel hierarchical BO to reduce the search space and optimize hyperparameters sequentially. For comparison, we refer to the previously introduced BO approaches without dimension dropout as _basic_ BO approaches.

A. Bayesian Optimization with Dimension DropoutBased on Li et al. (2018), we implement BO with dimension dropout as detailed in Algorithm 2. First, we derive \(m\) initial points via random search and add them to the set of observations \(R=\{(\theta_{1},f(\theta_{1})),(\theta_{2},f(\theta_{2})),...,(\theta_{m},f( \theta_{m}))\}\) (l. 3). Then, we conduct \(n\) iterations of BO with dimension dropout. In each iteration \(\iota\), we randomly draw \(d\) dimensions out of the problem's search space dimensions \(\mathcal{D}\) and refer to the subset of drawn dimensions as \(D^{\prime}\). We ensure that \(d=|D^{\prime}|<|\mathcal{D}|\). To dynamically adapt \(d\) in accordance with the total number of dimensions, we calculate \(d\) by \(d=\left\lfloor\tilde{d}*|\mathcal{D}|\right\rfloor\), where \(\tilde{d}\in[0,1]\) (l. 6). We denote the parameters of the drawn dimensions by \(\theta^{[\mathcal{D^{\prime}}]}\), and refer to the remaining parameters by \(\theta^{[\mathcal{D^{\prime}}]}\). Similarly, we denote the search space for the drawn dimensions by \(\Theta^{[\mathcal{D^{\prime}}]}\) and the best-observed values by \(\theta^{*}=(\theta^{*})^{[\mathcal{D^{\prime}}]}\cup(\theta^{*})^{[\mathcal{D ^{\prime}}]}\). We apply a Gaussian process surrogate model and optimize the acquisition function considering only the drawn dimensions (l. 7) as follows

\[\theta_{i}^{[\mathcal{D^{\prime}}]}\leftarrow\arg max_{\theta^{[\mathcal{D^{ \prime}}]}\in\Theta^{[\mathcal{D^{\prime}}]}}\mathcal{A}(\theta^{[\mathcal{D^ {\prime}}]}). \tag{4.12}\]

We select the values for the remaining parameters, \(\theta_{t}^{[\mathcal{D^{\prime}}]}\), by applying a _Dropout-Mix_ strategy (Li et al., 2018): With a probability of \(p\), we randomly draw the values for the remaining parameters from their respective domains. Otherwise, with a probability of \((1-p)\), the parameter values yielding the best observed function value are copied for the left-out dimensions, i.e., \(\theta_{t}^{[D\setminus D^{\prime}]}=(\theta^{*})^{[\mathcal{D^{\prime}} \setminus\mathcal{D^{\prime}}]}\)
(ll. 8-13). We evaluate the target function and update the set of observations \(R\) and incumbent parameters (ll. 14-16). Finally, we return the incumbent parameter setting \(\theta^{*}\) (l. 18).

B. Hierarchical Bayesian OptimizationWe further present a hierarchical approach to apply BO in high-dimensional problem settings by decomposing the search space. Here, we optimize sets of hyperparameters sequentially, enabling the application of distinct surrogate models for each set of hyperparameters. Further, we decompose high-dimensional search spaces which exceed the capabilities of basic BO approaches, such that the resulting sub-problems can be solved via BO methods for low-dimensional search spaces. We show the pseudo-code for the hierarchical approach in Algorithm 3.

```
1:Given: target function \(f(\cdot)\), domain \(\Theta\), number of iterations \(n\), number of initial points \(m\), acquisition function \(\mathcal{A}\), observed results \(R=\emptyset\), search dimensions \(\mathcal{D}\), percentage share of search dimensions to be drawn \(\tilde{d}\), probability for fill-up strategy \(p\)
2:Output: Incumbent parameter setting \(\theta^{*}\)
3:\(R=\{(\theta_{1},f(\theta_{1})),...,(\theta_{m},f(\theta_{m}))\}\leftarrow\text {RandomSearch}(f(\cdot),\Theta,m)\)
4:\(\iota\gets m\)
5:for\(\iota=\iota+1\) to \(\iota+n\)do
6: Randomly select a subset of dimensions \(\mathcal{D}^{\prime}\) such that \(|\mathcal{D}^{\prime}|=\left\lfloor\tilde{d}*|\mathcal{D}|\right\rfloor\)
7:\(\theta_{i}^{[\mathcal{D}^{\prime}]}\leftarrow\arg max_{\theta^{[\mathcal{D}^{ \prime}]}\in\Theta^{[\mathcal{D}^{\prime}]}}\mathcal{A}(\theta^{[\mathcal{D}^ {\prime}]})\)
8:\(q\leftarrow\) random number between 0 and 1
9:if\(q<p\)then
10:\(\theta_{\iota}^{[\mathcal{D}^{\prime}]}\leftarrow\) random values within domain \(\Theta^{[\mathcal{D}^{\prime}]}\)
11:else
12:\(\theta_{\iota}^{[\mathcal{D}^{\prime}]}=(\theta^{*})^{[\mathcal{D}^{\prime} ]}\)
13:endif
14:\(\theta_{\iota}\leftarrow\theta_{\iota}^{[\mathcal{D}^{\prime}]}\cup\theta_{ \iota}^{[\mathcal{D}^{\prime}]}\)
15:\(R=R\cup(\theta_{\iota},f(\theta_{\iota}))\)
16:\(\theta^{*}=\arg\min_{\iota^{\prime}\leq_{\iota}}f(\theta_{\iota^{\prime}})\)
17:endfor
18:return\(\theta^{*}\)
```

**Algorithm 2** Bayesian Optimization with Dimension Dropout

As a basis, we assign the search space dimensions to disjoint sets \(k=1,...,K\). We express the assignment of search dimension \(\Lambda_{j}\) to set \(k\) with a binary variable \(\gamma_{\Lambda_{j}k}\):

\[\gamma_{\Lambda_{j}k}=\begin{cases}1&\text{if $\Lambda_{j}$ is assigned to set k}\\ 0&\text{otherwise}\end{cases}\forall k=1,...,K,\Lambda_{j}\in\Theta,j=1,...,J. \tag{4.13}\]

We assign each search dimension \(\Lambda_{j}\) to exactly one set such that

\[\sum_{k=1}^{K}y_{\Lambda_{j}k}=1\ \ \forall\Lambda_{j}\in\Theta,j=1,...,J. \tag{4.14}\]

We initialize the hyperparameters by applying random search (ll. 3-5). Second, we apply BO for all hyperparameter sets \(k=1,...,K\) sequentially. Thus, we determine the search space
of set \(k\) and refer to the parameters optimized in the set \(k\) by \(\theta^{[k]}\) (l. 7). We fix the remaining parameters denoted as \((\theta^{*})^{[K\setminus k]}\) by copying the values of the incumbent parameter setting \(\theta^{*}\) (l. 8). Before conducting BO for set \(k\), we determine \(m^{[k]}\) initial points via random search (l. 9). Before each function evaluation, we merge the varied parameters' values \(\theta^{[k]}\) with the fixed parameters' values \((\theta^{*})^{[K\setminus k]}\). We conduct \(n^{[k]}\) iterations of BO for set \(k\) (ll. 11-15). After optimizing the hyperparameters of set \(k\), we continue with the next set until all sets of parameters have been optimized.

## 5 Case Study: Ambulance Demand Prediction for Seattle

We conduct a numerical case study for Seattle's 911 call data1 considering incidents of the category _Medic Response_, i.e., incidents requiring paramedical staff qualified for Advanced Life Support during 2 years (2020-2021). We divide this period into 8-hour intervals, corresponding to commonly applied 8-hour shifts in EMSs. By including ambulance demand of previous periods, we consider local, short-term dynamics such as local virus outbreaks during the Covid-19 pandemic. We divide Seattle into a grid of \(11\times 6\) subregions each with a dimension of approximately \(2.5\times 2.5\) km and apply a look back \(L\) of 6 periods, i.e., we take into account the ambulance demand of two preceding days.
FeaturesWe base our feature choice on the literature review in Section 2 and focus on features that can only hardly be learned with our plain CNN architecture. We neglect infrastructural and demographic information, e.g., the mobility infrastructure or the age distribution among the population, as we only conduct short-term forecasts and these variables only change gradually over time.

Since the weather changes dynamically, we include the following data from _Weather Underground2_: temperature, wind speed, humidity, dew point, sea level pressure, and precipitation. Except for the precipitation, we consider the daily minimum, maximum, and average for all weather features. As some of these inputs may be correlated, we conduct a correlation analysis and eliminate features with a correlation of \(>80\%\). We present the results of the correlation analysis in Appendix A.2. We assume that the spatial distribution of weather data is negligible for our case study since the area is sufficiently small to assume that the measurements of different weather stations within this area are highly correlated. For larger areas, our generic approach enables the inclusion of two- or three-dimensional weather data, e.g., for (spatial) weather predictions or historic (spatio-temporal) weather data, correspondingly.

Footnote 2: [https://www.wunderground.com](https://www.wunderground.com)

However, Wong & Lin (2020) show that temporal temperature changes influence people's health conditions. Also, age, income, sex, and health conditions influence people's sensitivity to weather conditions. Thus, the impact of temporal weather changes may vary between subregions. Therefore, we perform transposed convolution on the historic one-dimensional weather data to learn its spatial impacts while the distribution of relevant socio-economic factors can be implicitly learned by the model.

Many studies include special-day effects for periods with an anomalous ambulance demand, e.g., during school holidays (Vile et al., 2012) or New Year's Day (Channouf et al., 2007). To consider such effects, we take public holidays, school holidays, and events3 into account. We represent holidays as binary features. As events take place at different locations across the depicted area, we include their locations on the two-dimensional, spatial level. We additionally consider the number of expected participants at each event to account for the events' sizes.

Footnote 3: [https://data.seattle.gov/](https://data.seattle.gov/)

We further include time information, as Channouf et al. (2007) show that the day-of-week, month-of-year, and hour-of-day impact ambulance demand. We embed this information as one-hot-encoded vectors. To increase the models' performance, we scale the input data such that all values range within \([0,1]\).

BenchmarksThe developed CNN serves as a representative example of an integrated spatio-temporal ANN architecture. To evaluate the numerical benefits of incorporating spatio-temporal dependencies in such architectures, we compare this approach to an iterative ANN architecture, applying an MLP as an exemplary model. MLPs are fully connected feedforward ANNs commonly applied for forecasting tasks, also in the domain of ambulance demand prediction (Setzler et al., 2009; Chen et al., 2015). We refer to Equation 4.7 for a definition of a fully connected architecture.
We further aim at evaluating the performance of the proposed intrinsic feature selection algorithm which integrates feature selection in the hyperparameter tuning process of the corresponding ANN architectures. For this evaluation, we compare our approach to other commonly applied intrinsic feature selection methods that are able to solve prediction tasks. One approach widely applied for regression tasks and capable of providing valuable insights into features' importance which make outcomes interpretable for practitioners are decision trees. However, decision trees are prone to overfitting, especially for high dimensional search spaces as there is mostly an insufficient number of samples for each parameter value. To apply a model more robust to overfitting, we additionally apply random forests, i.e., ensembles of multiple decision trees based on randomly selected data samples. The average prediction generated by the individual decision trees constitutes the random forest's prediction. In line with Setzler et al. (2009), Zhou et al. (2015) and Zhou & Matteson (2016), we further compare the performance of the developed CNN to a common industry practice, i.e., the Medic method. This method takes the average ambulance demand of similar time periods of the past 4 weeks, i.e., the historic ambulance demand associated with the same weekday and time, of the current and preceding years as prediction. To enable a fair comparison to the ANNs trained on data from 2020-2021, we consider the same data range for the Medic method.

Hyperparameter TuningAs the hyperparameter choices for decision trees and random forests are limited, we apply grid search for tuning. However, for the CNN and MLP, we optimize the hyperparameters by applying i) random search ii) BO with Gaussian process, random forest, and extremely randomized tree surrogate models iii) BO with dimension dropout, and iv) the novel hierarchical BO with Gaussian process, random forest, and extremely randomized tree surrogate models. First, we perform these approaches without feature selection, and second, including feature selection. We present the hyperparameters and their domains in Appendix A.3. We further show the set assignments for the hierarchical BO.

We determine the order in which we optimize the hyperparameter sets based on two criteria: The degree to which preliminary experiments can yield adequate values, and, the size of the search space. The better we can determine adequate hyperparameter values, the later we optimize them. This ensures that we optimize the hyperparameters that are more complex to determine on an appropriate basis, i.e., the hyperparameters that are currently not optimized are set to well-performing values which were either determined by applying random search for initialization, or, optimized at a previous level. Later, we can fine-tune the hyperparameters we fixed in the first place.

For our application, we distinguish three levels. At each level, one parameter set is tuned. The first level focuses on the main architecture decisions, such as the number of layers and filters per layer. The second level tunes the activation functions. The third level controls the regularization mechanisms to avoid overfitting and decides on the batch size, learning rate, and optimizer. In the case of incorporating feature selection, we include an additional binary hyperparameter for each feature, stating whether the feature is included (1) or not (0). In the hierarchical BO, we therefore
include an additional level for these decisions. We optimize this level first, such that we base the choice of the remaining hyperparameters tuning the model on the derived features.

TrainingWe train our model on 60% of the data and split the remaining data into a test and a validation set, each containing 20% of the data. We apply a mean-squared-error (MSE) loss function and train the model for 300 epochs. In the case that the validation loss does not improve for 20 epochs, we abort the training process and measure the MSE on the test set. We initialize the hyperparameters of the hierarchical BO based on 500 iterations of random search. To enable a comparison between the hierarchical approach and the remaining methods, we initialize them equally. After initialization, we execute 1000 iterations. As the number of parameters optimized at each level of the hierarchical BO is limited, we reduce the number of iterations to 250, applying random search for the first 25 iterations per hyperparameter set.

We implement the algorithms in Python 3.8 using TensorFlow 2.10.1. We perform BO using Scikit-Optimize 0.9.0, which we adapt to apply dimension dropout. We train the CNNs on an NVIDIA A100 GPU with 80 GB RAM. The MLPs are trained with an Intel(R) Xeon(R) processor E5-2697 v3 with 15.7 GB RAM.

## 6 Experimental Results

We first compare the performance of the developed CNN architecture against its benchmarks and show the results we obtain for the different hyperparameter tuning approaches. Second, we analyze the selected features when conducting feature selection and calculate the contribution of each feature to the prediction. Third, we investigate the CNN's performance for different forecasting time intervals.

### Model and Hyperparameter Tuning Performances

In Table 2 we compare the MSE of the developed CNN architecture and its benchmark models, i.e., the Medic method, MLPs, decision trees, and random forests. The learning-based approaches were tuned with the introduced hyperparameter optimization approaches and initialized equally with random search.

As can be seen, the incumbent CNN achieves a 9.83% lower MSE than the best MLP. The CNN further outperforms the Medic method, decision tree, and random forest by 9.98%, 14.84%, and 11.26%. As expected, we see that the decision tree badly generalizes to unseen data leading to an incumbent tree with a depth of five. The incumbent random forest has a depth of eight. We visualize their tuning results in Appendix A.3. Although the random forest performs better on new data, it still remains inferior to the ANNs. Both tree-based approaches are unable to identify relevant patterns due to their limited depth and the resulting negligence of relationships among external features. Although the Medic method bases its prediction solely on historic ambulance demand neglecting all external features, it still outperforms the tree-based approaches and yields
comparable results to the MLP. Contrarily, the CNN benefits from its ability to detect spatio-temporal patterns across subregions within historic ambulance demand and combines these patterns with external information. This ability also encourages the use of the presented feature selection approach which enables the application of the developed CNN instead of the intrinsic benchmarks, i.e., decision trees, and random forests. The advantage of the introduced intrinsic approach is the exchangeability of the model: We are not bound to an intrinsic tree-based approach to make the prediction. Instead, we can apply the developed CNN resulting in the lowest MSE. We further observe that retaining less than 50% of the features in the CNN reduces the number of trainable parameters by 40.4% and only slightly increases the MSE by less than 0.05%. For practitioners, this performance decrease remains negligible and a leaner model is easier to maintain, as less input data must be collected, validated, and pre-processed.

**Result 1**.: _The CNN architecture outperforms the incumbent MLP, Medic method, random forest, and decision tree by 9.83%, 9.98%, 11.26%, and 14.84%, correspondingly._

**Result 2**.: _The tree-based approaches badly generalize to unseen data, resulting in shallow trees which are unable to identify underlying patterns._

**Result 3**.: _Incorporating feature selection in the hyperparameter tuning process reduces the number of trainable parameters by 40.4% while the performance decrease of less than 0.05% remains negligible._

When including feature selection decisions in the hyperparameter tuning process, we integrate a binary variable for each feature selection decision. This exponentially increases the search space. However, it enables the model to simultaneously decide on features and make hyperparameter decisions. This is beneficial, as the CNN's architecture, e.g., the number of layers in the model, depends on the selected features. Therefore, we evaluate the influence of including feature selection in the hyperparameter tuning process. We compare the convergence of the different tuning approaches with and without feature selection for the CNN in Figures 2 and 3. Based on the hyperparameters

\begin{table}
\begin{tabular}{l||c|c|c c|c c c|c c c|c||c c} \cline{2-13}  & No tuning & Grid Search & \multicolumn{2}{c|}{Random Search} & \multicolumn{2}{c|}{BO with} & \multicolumn{2}{c||}{} & \multicolumn{2}{c||}{} & \multicolumn{2}{c||}{} & \multicolumn{2}{c||}{} & \multicolumn{2}{c||}{} \\ \cline{2-13} Model & & & & & & & & & & & & & & & & & \\ \hline CNN & - & & & & & & & & & & & & & & & \\ CNN & - & - & 14.74 & 14.74 & 14.74 & 14.74 & 14.72 & 14.73 & **14.67** & 14.74 & 14.70 & 14.70 & 14.70 & \(\checkmark\) \\ CNN & - & - & 14.72 & 14.72 & 14.71 & 14.72 & 14.72 & 14.72 & 14.70 & - & 14.67 & 14.66 & **14.66** & \\ MLP & - & - & 16.34 & 16.31 & 16.34 & 16.30 & 16.28 & 16.27 & **16.26** & 16.27 & 16.27 & 16.27 & \(\checkmark\) \\ MLP & - & - & 16.32 & 16.32 & 16.32 & 16.29 & 16.28 & 16.28 & 16.29 & - & 16.29 & 16.29 & **16.28** & \\ RF & - & **16.52** & - & - & - & - & - & - & - & - & - & - & - & \(\checkmark\) \\ DT & - & **17.22** & - & - & - & - & - & - & - & - & - & - & - & \(\checkmark\) \\ Medic & **16.29** & - & - & - & - & - & - & - & - & - & - & - & - & \(\checkmark\) \\ \hline \# iter. & & & & 500 & 1500* & & & & & 750* & 1000* & 1250* & \\ \hline \# iter. & & & & 500 & 1500* & & & & & & 750* & 1000* & 1250* & 1500* & \(\checkmark\) \\ \hline \end{tabular}
\end{table}
Table 2: **Comparison of MSE (*incl. 500 iterations random search for initialization) [\(\times\)100].** We compare the MSE of the convolutional neural networks (CNNs), multilayer perceptrons (MLPs), decision trees (DTs), random forests (RFs) and the Medic method. For hyperparameter tuning, we apply Bayesian optimization (BO) with Gaussian process (GP), RF and extremely randomized tree (ET) surrogate models. Hierarchical BO is applied for sets \(k\in\{1,2,3\}\). We add set k=0 when conducting feature selection. We train the model on 60% of the data and apply a validation and test set, each containing 20% of the data.
determined by the corresponding tuning approach in each iteration, we show the MSE achieved on the test set after training the model. In Figure 2 we neglect feature selection in the hyperparameter tuning process. Here, the basic BO can only slightly decrease the MSE in 1000 iterations when utilizing a Gaussian process surrogate model. In contrast, BO with dimension dropout and the hierarchical BO can better improve the MSE. In Figure 3 we include feature selection decisions. Here, the basic BO can only slightly decrease the MSE in 1000 iterations when utilizing a random forest surrogate model. For including feature selection, we add a binary variable for each selection decision, which a tree-based surrogate can better optimize than Gaussian processes. However, for both cases, i.e., with and without feature selection, the BO with dimension dropout and the hierarchical BO are superior to tackle the high dimensional search space.

**Result 4**.: _The BO with dimension dropout and the hierarchical BO outperform basic BO approaches and random search._

### Features' Importance

To further investigate the reasons for the CNN's superiority, we analyze the importance of the external features for the different models. We base this analysis on the incumbent models derived in Section 6.1 for 8-hour time-intervals. Here, we distinguish two cases: first, when including feature selection in the hyperparameter tuning approach, we present the selected features in Table 3. Second, when excluding feature selection in the hyperparameter tuning approach, we keep all features and present their importance in Figure 4. Since the decision tree always conducts feature

Figure 2: Convergence of the hyperparameter tuning approaches of the CNN without Feature Selection (after applying 500 iterations random search for initialization)
selection, we only present its selected features. As the random forest is an ensemble of decision trees with randomly selected features, we only present its features' importance.

Selected FeaturesTable 3 shows that the amount and type of features selected by the ANNs and decision tree strongly vary. The tuning process of the decision tree reveals that the tree overfits for high depths, i.e., we observe that the model's performance decreases on unseen data when increasing the tree's depth beyond five, visualized in Appendix A.3. We further see that the decision nodes' conditions made at depths below four are exclusively based on historic ambulance demand. Only in the final splits, weather features, e.g., the maximum wind speed of past periods, or the time-of-day can be decisive. As these external features are only used by decision nodes preceding the leaf nodes, the model cannot identify relevant patterns among these features.

Figure 3: Convergence of the hyperparameter tuning approaches of the CNN with Feature Selection (after applying 500 iterations random search for initialization)

\begin{table}
\begin{tabular}{l||c c c|c c c c c c c c c c c c c c c c c c c} Features & \multicolumn{6}{c}{Historic data from periods \([t-L,t-1]\)} & \multicolumn{6}{c}{Data for predicted period \(t\)} \\ \hline
**Result 5**.: _The decision tree only uses few external features in decision nodes preceding the leaf nodes, preventing it from identifying patterns among these features._

In contrast to the decision tree, the ANNs combine information from multiple external features. Table 3 shows that although the features selected by the CNN and MLP differ, they mostly include similar types of information; while one model selects a feature with predicted data, the other model includes the same type of information but includes the feature with the corresponding historic data, and vice versa. For example, both models include information about the sea level pressure and humidity. This is in line with many studies showing the effect of humidity and air pressure on people's health conditions influencing ambulance demand, e.g., the risk for cardiovascular diseases (Ou et al. 2014, Borghei et al. 2020). In addition, both models include the hour and weekday which confirms their relevance observed in earlier studies (Setzler et al. 2009, Matteson et al. 2011). Simultaneously, some features are neglected by all models. As we include historic ambulance demand of previous periods, some features whose values remain unchanged for several periods, can be indirectly learned by the model. For this reason, the month and school holidays, mostly lasting several days, are neglected.

**Result 6**.: _The CNN and MLP select common types of information, such as weather data that has already been identified in medical studies to harm people's health conditions. In addition, in line with earlier studies, our results confirm the relevance of considering time information, i.e., the time-of-day and weekday._

Features' ImportanceIn Figure 4 we present the features' importance derived for the case of retaining all features during the tuning process. To express the importance, we calculate the SHAP values using shap 0.41.0 (Lundberg & Lee 2017). We determine these SHAP values for 50 randomly selected time periods covering all subregions from the test set. Each value allows us to quantify the feature's contribution to the corresponding prediction compared to the base value. We calculate the average prediction output over 300 random time periods over all subregions from the training set as a base value.

While historical data highly contributes to the prediction of the MLP and random forest, the CNN's output is mainly influenced by data characterizing the period to be predicted. This indicates that the CNN is more robust to changes in the upsampled and three-dimensional input data. Contrarily, the MLP and random forest are more sensitive to such changes since they conduct local forecasts for each subregion separately. Thus, changes in the historical data of the corresponding subregion strongly influence the forecast. In contrast, the CNN considers the historic data of all subregions simultaneously. Thus, it is more robust to changes in the historic data of a single subregion, which is advantageous, e.g., for handling demand outliers or data inconsistencies.

**Result 7**.: _The features' importance strongly differ between the MLP, random forest, and CNN. The CNN is more robust against changes in the upsampled data and historic ambulance demand._
### Sensitivity Analysis: Time Granularity

Earlier studies address the challenge arising from sparse ambulance demand which we face in the case of fine spatio-temporal granularities (Zhou et al. 2015, Zhou & Matteson 2016). For an extensive study of an iterative ANN architecture, we refer to Setzler et al. (2009) who investigate the influence of fine spatio-temporal granularities on the performance of an MLP and compare it to the Medic method. To close this gap for the domain of integrated spatio-temporal ANN architectures, we conduct a sensitivity analysis and investigate the CNN's performance for five different time intervals: 2 hours, 4 hours, 8 hours, 12 hours, and 24 hours. Similar to Setzler et al. (2009), we compare the performance of the developed CNN against the Medic method. Previous results show that the hierarchical BO and BO with dimension dropout yield superior results for high-dimensional search spaces. For this reason, we apply these approaches in the following experiments when tuning the

Figure 4: Features’ Importance: SHAP values
CNN. For both models, we calculate the MSE obtained on the test set for each time interval. In addition, we compute the MSE for instances with zero and non-zero demands individually. Moreover, to take the amplitude of each time interval into account, we include the normalized root-mean-square error (NRMSE) in the analysis, calculated as follows

\[NRMSE=\frac{\sqrt{MSE}}{y_{max}-y_{min}}, \tag{6.1}\]

where \(y_{max}\) and \(y_{min}\) are the maximum and minimum observed ambulance demand across all time intervals and subregions. We show the obtained MSEs and NRMSEs for the incumbent CNN and Medic method in Figure 5.

For all analyzed performance measures and time horizons, the CNN outperforms the Medic method. We see that both demand types, i.e., instances with zero and non-zero ambulance demand, can be better predicted by the CNN. The lower the granularity, the better the performance of the CNN for non-zero demands compared to the Medic method. This indicates, that in the case of higher demands, the CNN can better identify underlying patterns within the data. Considering the time interval's amplitude in this analysis, we further see that the NRMSE remains stable across all time intervals for both models. This shows that the prediction performance remains stable if we take into account the time intervals' emergency call volumes. Disregarding these amplitudes, the MSEs and the NRMSEs improve for finer time granularities due to the reduction in ambulance demand and the increase in zero-demands which can be well predicted.

**Result 8**.: _The CNN outperforms the Medic method for all analyzed time intervals. Both demand types, i.e., instances with zero and non-zero ambulance demand, can be better predicted by the CNN._

Figure 5: Sensitivity Analysis: Time Granularity
## 7 Conclusion

We presented a novel CNN architecture that transforms time series information into heatmaps to predict ambulance demand. It can incorporate historical ambulance demand and external features of varying dimensions and uses three-dimensional convolutional layers to detect correlations in space and time. We further introduced a hyperparameter tuning framework utilizing BO that enables to intrinsically tune the hyperparameters while selecting features. To tackle the high-dimensional search space, we applied BO with dimension dropout and introduced a novel hierarchical BO. We further calculated SHAP values to investigate the contribution of each feature to the prediction. We further analyzed the CNN's performance for different forecasting horizons. To apply the developed CNN to real data, we conducted a numerical case study for Seattle's 911 call data.

Results show that the CNN outperforms all benchmarks by more than 9%. Although the ANNs select different features, they include similar types of information confirming their relevance observed also in earlier studies in the domain of ambulance demand prediction and medical studies investigating the influence of weather data on people's health conditions. We further showed that the CNN is more robust against changes in the three-dimensional and upsampled input data. When including feature selection in the hyperparameter tuning process, we observed that excluding more than 50% of the features reduces the trainable parameters by 40.4% while the performance loss of \(<0.05\%\) is negligible. Moreover, the CNN is superior to the current industry practice for all analyzed time horizons, ranging between 2-hour and 24-hour time-intervals.

In future work, the CNN architecture could be further optimized for fine space and time granularities, e.g., by giving non-zero demands more weight. Moreover, the derived features' importance could be used to improve the feature selection process. In addition, further investigating the convergence behavior of the BO when including feature selection decisions remains an interesting subject for future research.
