{'fraction_non_alphanumeric': 0.032430066451792604, 'fraction_numerical': 0.010584537025983845, 'mean_word_length': 4.741832305231894, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatiotemporal Prediction of Ambulance Demand using Gaussian Process Regression

Accurately predicting when and where ambulance call-outs occur can reduce response times and ensure the patient receives urgent care sooner. Here we present a novel method for ambulance demand prediction using Gaussian process regression (GPR) in time and geographic space. The method exhibits superior accuracy to MEDIC, a method which has been used in industry. The use of GPR has additional benefits such as the quantification of uncertainty with each prediction, the choice of kernel functions to encode prior knowledge and the ability to capture spatial correlation. Measures to increase the utility of GPR in the current context, with large training sets and a Poisson-distributed output, are outlined.


Introduction
Fast ambulance responses are often vital in achieving a good clinical outcome following a medical emergency. Indeed, ambulance response time has been shown to correlate with mortality [27,29,3]. Accurate prediction of the density of the calls in time and space can improve readiness and therefore reduce response times. In addition, longer-term prediction of demand can inform staffing levels, ensuring that the workforce available is always sufficient to meet the demand. Here we assert that the locations and times of previous emergencies will inform the anticipation of future call-outs. We seek to fit patterns of the demand in time and space, which can be extrapolated to generate predictions for future timesteps. This paper begins with an overview of GPR, its suitability to spatiotemporal applications and how we go beyond standard GPR to improve performance with a challenging problem. We review existing work before describing how we evaluate predictive performance. Next, we outline our method for call-out prediction -spatiotemporal GP (STGP) -and present results, which are compared against those of the MEDIC benchmark. We end with some suggestions as to how this work might be exploited or extended.

Gaussian Process Regression
By definition, a Gaussian process is a stochastic process in which any finite subset of the random variables are Gaussian distributed [22]. They are commonly used as prior distributions in Bayesian inference. Covariance functions, which capture the similarity in output for a given pair of inputs, can be chosen based on known characteristics of the system, seasonality and trends for example. GPR has exhibited promising performance in a wide variety of spatiotemporal applications. A nonexhaustive list includes: determining the distribution of available protein (in space and time) during post-transcriptional regulation [26]; mapping brain structure developmental profiles of patients suffering from neurodegenerative diseases [16]; monitoring the spread of gastroenteritis [9], influenza [31] and malaria [24]. A number of properties of GPR are thought to produce superior performance in spatiotemporal applications. Being a member of the Bayesian nonparametric class of methods, the output can be modelled without knowledge of all its dependencies or an assumed parametric form. Though, prior knowledge or intuition for types of spatial or temporal variation can be incorporated in expressive covariance functions. Further, large training sets are often required to capture spatiotemporal interactions and while traditional GPR is not scalable, recent research into faster alternatives makes GPR with large training sets accessible. In addition to benefits of GPR listed, we see the ability to forecast with uncertainty as a significant advantage in the current application. Articulation of uncertainty allows any action made in light of predicted demand to be properly qualified.

Challenges
There are a number of challenges in applying GPR here: 1. Non-Gaussianity. Standard GPR assumes a Gaussian distributed noise model. In predicting ambulance demand, the output is the number of counts in a given bin, and therefore a Poisson noise model is more appropriate. 2. Number of Training Instances. As ambulance call-outs must be binnned in time and space, the number of training samples, n, scales cubically with the three-dimensional resolution. Any resolution sufficient to capture the spatiotemporal dynamics results in a large number of training samples. Standard GPR involves matrix inversion, with undesirable O n 3 computational complexity. 3. Sparsity. As observed by [40], the number of call-outs in most spatial bins is zero for most time steps. This can create a poor signal to noise ratio which causes difficulties during model training.

GPR for Poisson Regression
To make use of the desirable properties of GPR while correctly capturing the discrete nature of the output, we use a log-Gaussian Cox process (LGCP) [23]. A Cox process is an inhomogeneous Poisson process where the underlying rate λ = λ(x) is a realisation of a stochastic process, denoted as ρ(x). A LGCP extends a Cox process to model ρ(x) as a transformed GP. The complete model description is f ∼ GP(µ(X), k(., .)) (1) where µ(.) and k(., .) are the mean and covariance functions of the GP, y = {y i } n i=1 are the observed outputs, f = {f i } n i=1 is the GP model fit, and X is the design matrix. Note that the exponential transformation of f (.) ensures non-negativity and thus a well-defined Poisson process.

Scalable GPR
Many schemes have been proposed to make inference using GPs scalable [35,15,13,33,7,25,36]. We use the stochastic variational inference (SVI) approach outlined in [14] which uses an inducing point framework for variational approximation of the posterior p(f |y, X) and marginal likelihood p(y|X) distributions. This method was chosen over alternatives because of its apparently high accuracy which can be balanced against computational complexity by tuning the number of inducing inputs. Note, all distributions below are implicitly conditioned on the observed inputs, X, which we drop from the notation for brevity. The approach begins at the inequality Combining (4) and (5), and exploiting the factorisation of the likelihood, p(y|f ) = n i=1 p(y i |f i ), can be shown to give the bound As f ∼ GP(µ, k(., .)) identities for the derivatives of Gaussian-expectations (see [28]) can be exploited to find a convenient form for the derivative of the bound. The derivatives with respect to the parameters of the variational distribution, inducing inputs Z and covariance hyperparameters can be found, and thus their values optimised by gradient descent. The predictive distribution can then be calculated by integrating the approximate posterior over u and f for a novel input x * .

Related Work
A large number of early studies focused on ambulance demand modelling investigated how quasistatic socio-economic attributes of a geographic region correlated with the ambulance demand using linear models [1,20,5,30,19,34,8]. While these methods are of interest in understanding the causes of ambulance demand, they are not directly relevant to short-term demand prediction. Nonparametric time-series methods have been implemented for forecasting in this context. A nonexhaustive list includes autoregressive moving average (ARIMA, [37]); vector autoregression (VAR) to correlate demand at multiple sites [18]; Singular Spectral Analysis (SSA, [38]) and Winters' triple exponential smoothing [2]. MEDIC is a common benchmark approach which is discussed in [32]. MEDIC exploits seasonality at different scales; the prediction for a given hour is the average of the call rates observed in the same hour of the previous weeks over the past month, and those same hours in previous years. MEDIC was said to be industry practice for Mecklenberg County, North Carolina [32]. Kernel density estimates (KDEs) are used for demand prediction over continuous geographic space and hourly timebins in [39] and [40]. Predictions are generated by first creating a spatial field made up of two-dimensional Gaussian kernels, one at each historic observation. They are then weighted for a given test input according to their geographic proximity, similarity in hour-of-day and similarity in hour-of-week. This method is extended for Melbourne in [40]. Melbourne has a highly irregular city boundary which motivates kernel warping: tweaking the kernel by projecting it onto the graph Laplacian of previous call-out locations. By observing the variation in different strengths of seasonality between locations, a time-varying Gaussian mixture model (GMM) is proposed in [41]. The weight of each Gaussian component is modulated by the phase of a prediction time within the weekly cycle. Multiple studies have implemented Artificial Neural Networks (ANNs) for callout prediction, most notably in Taipei [6,21] and Mecklenberg County, North Carolina [32]. While the studies for Taipei show promising accuracy at coarse test resolution, they struggle due to sparsity at the scale of resolution required here. The resolutions tested in [32] are closer to that required here. Using spatial grids with 4 × 4mile cells, the ANN is found to have slightly improved mean square error scores relative to MEDIC. Given the above review, GPR is thought to be an appropriate method for the following reasons: • Bayesian inference ensures robustness in a poor signal-to-noise-ratio regime.
• The seasonality in demand observed in other studies can be captured with periodic covariance functions. Further, composing the overall covariance as the product of those for the temporal and spatial components allows location-specific seasonal patterns to be captured.
• GPR uses principled model fitting, rather than direct averaging over historic call-out locations and times. This limits the effect of noise in the training data on predictions.
• Standard GPR also predicts a distribution over call-densities, without requiring empirical approximation of uncertainty.

Evaluation
We use two methods to evaluate the predictive performance of the implemented models. Our predictions, which are estimates of the inhomogeneous Poisson rate, are denoted as λ(s, t) where s is the geographic space and t the call time. The set of observed call times and locations is denoted S.
1. Inhomogeneous Poisson Log-likelihood Density. This measure reflects the probability of a set of call outs, given the predictions λ(s, t) are the true underlying rate of the Poisson process. Its mathematical form can be found by taking the product of probabilities of zero events over intervals with no events and probabilities of single events over small intervals δsδt around call-out times and locations. Taking the limit δ → 0 and noting that only a (a) The number of calls recorded daily over the training set.
(b) Diurnal seasonality of ambulance demand in the training set. Figure 1: The temporal variation of call-out rates.
single ordering of events is realised, the following likelihood density is found: Intuitively, we can interpret the first term as rewarding the prediction of high densities at event locations, the second term as penalising the overall prediction density and the final term as a normalising constant which reflects event ordering. We approximate the integral as a summation over a 50×50 spatial grid (covering the region shown in Fig. 2b) at 4hourly resolution and evaluate the log k! term using the relation Γ(k + 1) = k! and the log-gamma function in Python programming language's SciPy package [17]. In general, STGP generates small but non-zero predictions in regions with low call density. However, MEDIC averages previous observations, and will often predict zero events in a low density region. To avoid an undefined log 0 when an event occurs in a region with zero predicted counts, we clip our predictions λ(s i , t i ) at 10 −4 km −2 hour −1 when evaluating the first term in (7).

Mean Absolute Error.
We also use a more easily interpretable measure of accuracy. The mean absolute error (MAE) is given by where c i is the count of events in bin i, with centre (s i , t i ); A is the geographic area covered by a test bin and τ is the test timebin length. The Aτ term rescales the predicted rate to a prediction of the number of counts within the test bin.
While MAE is straightforward and readily communicable to a non-technical audience, its value is dependent on how the data are binned. This introduces subjectivity into the measurement. Conversely, the log-likelihood is independent of how the input space is descritised. As it is assumed that ambulances will not be called-out to locations at sea, so any bins which do not cover any land are removed from the MAE calculation. Similarly, the geographic element of the integral s t λ(s, t) ds dt is calculated over land only.

Materials
The dataset comprises around 4.

Method
Here we present details of the kernel function and algorithm we use in our implementation of SVI LGCP which make up STGP. All analysis was carried out in the Python programming language, and the GPy library [12] was used for Gaussian process implementation 1 .

Preprocessing
Given the concentration of demand in Cape Town, we restrict the geographic domain to its approximate limits 2 . As the demand is greatest in this region, we expect accurate predictions in the region to be most valuable. The dataset contains a small number of non-emergency events such as inter-hospital transfers. As these are likely to obscure the true distribution of emergency events, they are removed. Sparsity is a recognised problem in other studies [39,40,41,32]. Our dataset comprises an average of around 1.5 calls per hour over the whole Western Cape, around 1.1 calls per hour within Cape Town. This translates to average rates of around 1.2 × 10 −5 and 6.9 × 10 −4 km −2 hour −1 respectively. For training we divide the geographic domain (see Fig. 2b) into a 6×6 grid within which events are counted for each 4 hour timebin. The dimensions of a cell in such a grid are 6.7km×6.7km×4hours. While this is a coarser level of descritisation than might be desirable, a finer grid was found to result in a poor signal to noise ratio which resulted in non-convergent training.

Covariance Function
Choice of covariance function is both a benefit and a drawback of using Gaussian process regression: it allows us to encode some high-level prior knowledge of the problem within our, but also introduces a seemingly arbitrary decision. While a method to remove subjectivity from this process have been proposed [10], it involves multiple training cycles which is not practical when the training set is large. Instead we use what is already known about the problem to inform our choice as far as possible.

Name Symbol Value
Time Period T 24 hours Time Lengthscale l t 8 hours Spatial Lengthscale l s 10 km Table 1: Fixed covariance hyperparameter values.

Temporal Covariance
While the daily call rate seems not to exhibit any significant variation over the training set (Fig. 1a), there is considerable hour-in-day seasonality (Fig. 1b). Inspection of day-of-week and month-ofyear seasonality did not yield any notable observations. We thus choose a single periodic kernel [22] for the temporal component of the input: where l t is the lengthscale of the process; θ t reflects its variance and T is the period of oscillation, in our case 24 hours.

Spatial Covariance
The spatial distribution of events (Fig. 2b) shows there is significant variation over the domain and some challenges caused by boundaries in the landscape. We choose the commonly used radial basis function (RBF) kernel for its simplicity and robustness: where σ 2 s is the variance, m is the number of dimensions (two in our case) and l s is the spatial lengthscale. We choose a relatively small value of l s in the hope that the issues caused by geographic boundaries may be mitigated.

Hyperparameter Tuning
Given the sparsity in the data, and therefore the poor signal to noise ratio, we expect the surface of the marginal likelihood, p(y|x) to be highly non-convex. This makes hyperparameter tuning by gradient descent methods problematic. We thus fix covariance hyperparameters at values known to be sensible from exploratory analysis and prior understanding, for all except the variance. These values are given in Table 1.

Composition
As described in [10], the sum or product of two valid kernels is itself a valid kernel, however summation and multiplication give rise to different behaviours. Multiplication captures interactions between kernels whereas summation models each as an independent process. This is analogous to AND and OR operators enacted on the similarity measured by each kernel [10]. Deciding which is appropriate is dependent on the level of spatiotemporal interaction in the data. Fig. 3 shows the spatial distribution of calls in the training set. Fig. 3a shows that between midnight and 1am, while Fig. 3b shows that between 3 and 4pm, and both are averaged over all days in the training set. While the general shape of the distribution is similar in both cases, there is a greater focus of event density in the downtown region between 3 and 4pm. We thus conclude that there is some spatiotemporal interaction. We choose the following kernel structure: where x = (s, t), x = (s , t ), both k t,m (., .) and k t,i (., .) are instances of k t (., .) without tied parameters, the same is true for k s,m (., .) and k s,i (., .). This composition is motivated by the majority of variance being marginal in either time or space, which can be captured by the first two components. We hope to capture some of the residual variance arising from spatiotemporal interaction using the product kernel.

Inducing Inputs
While the inducing point framework provides an interpretable trade-off between runtime and accuracy, it requires selection of the number of inducing inputs and their initial values. With SVI GPs [14], it is common practice to optimise the inducing inputs by gradient descent of the marginal likelihood surface. However, we find our marginal likelihood surface to be largely non-convex and optimisation in this manner to result at poor values of inducing inputs. To mitigate this, we randomly sample, without replacement, 180 inputs from the training set for use as inducing inputs. We find this number of points sufficient to achieve a good model fit.

Gradient Descent
As discussed in Section 2, we train STGP using a variant of gradient descent. There are a number of possible algorithms which could be used. We use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [11] as implemented in in Python language's SciPy library [17]. We find this converges in less time than Newton's method and scaled conjugate gradients.
(a) (b) Figure 4: The distribution of weekly predictive accuracies.

Benchmark
As discussed in Section 3, MEDIC is a benchmark of call-out prediction accuracy which has been applied in industry and is used in [39,40,32]. The method predicts the counts to be the average of those for the same timebin in the previous four weeks over the past 5 years [32]. This is generalised to the spatiotemporal domain, by applying it within each geographic area. However, no consensus could be found as to the size of each bin, which appears to be context dependent. [32] implements both 2×2 mile and 4×4 mile descritisation as both 1 and 3 hour intervals. [39] and [40] use 1 km 2 and 1 hour bins. Here, we divide the spatial domain into a 6×6 grid and implement MEDIC within each, over the range of training data available (∼3 years). We use this spatial resolution in conjunction with a 4-hourly temporal resolution as a compromise between capturing the spatiotemporal dynamics, and reducing the affect of noise through averaging.

Algorithm
To generate predictions over the test set (17th March to 14th September 2015 inclusive), we use a weekly rolling train/test scheme as described in Algorithm 1. The training set size for STGP is clipped at 6 months to reduce runtime. To ensure against a poor STGP model fit resulting from gradient descent converging to a local minimum, we assess the model likelihood after training. If the training log-likelihood falls below a set threshold (−1.07 × 10 4 ), we retrain the model with a different random seed. P ← append(P, P te ) 12: t te ← min(t te + 1week, t end )

Results
The test-set accuracies are shown in Table 2 and the distributions of weekly accuracies are given in Fig. 4. STGP exhibits superior predictive power according to both MAE and log-likelihood, for all quantiles. The difference in performance is most marked when measured by log-likelihood. This is likely to be caused by two differences between STGP and MEDIC. First, STGP predicts a smooth spatial distribution at each timestep, whereas MEDIC is restricted to the spatial grid used for training. This allows STGP to predict in a more accurate manner, with less predictive weight being wasted in in locations with few call-outs adjacent to locations of high density. Such an effect would cause STGP to have a smaller integral of predicted rate over time and space, and/or greater predicted rates at probable call locations. Further, events within bins where MEDIC has predicted a zero rate still reduce its test set log-likelihood score. While we round such predictions up to 10 −4 km −2 hour −1 , they still create significantly negative contributions to the summation n i=1 log λ(s i , t i ).
Method MAE Log-likelihood STGP 0.174 −6.83 × 10 4 MEDIC 0.181 −7.96 × 10 4 Table 2: Predictive accuracies over the full test set.  Exploring further the dependence of the accuracy on the value of the output , the distribution of residuals for all values of events per bin (during the MAE calculation) is given in Fig. 5. There are a number of comments to make here. First, STGP is much more accurate when it comes to predicting for bins with no events. It is likely that MEDIC is less accurate in these cases due to its quantised predictions: it predicts averages of relatively few noisy count observations, so even if one of the historic observations is non-zero, the resulting inaccuracy is significant. Average accuracies for bins containing events are more similar, however STGP is more consistent, with less variance in accuracy. It should be noted that neither method is accurate for bins with large number of events. However, Fig. 5 is misleading as it does not give an indication of the relative prevalence of each number of events. The large majority of bins contain zero events, few contain a single event and very few contain more than one event. It is expected that such skew results in under-prediction for bins with many events. Continuity of the input and output result in smaller bias in STGP predictions, relative to MEDIC. The average spatial distribution of residuals over the test set is presented in Fig. 6. In general, MEDIC exhibits positive forecast bias in areas of low demand. This is thought to result not only from quantisation of the output, but no encoding for spatial similarity between adjacent bins. The spatial kernel used in STGP enforces smoothness which reduces the effect of noise. Note also that MEDIC shows greater negative bias in one of the downtown regions with greatest call density. This is concerning as it could result in under-provision of resources to high-demand areas.

Conclusion
We have presented STGP, a novel method for spatiotemporal ambulance demand using Gaussian process regression (GPR). GPR has a number of desirable features in this context such as the ability to choose classes of pattern which the model can fit, and make predictions with uncertainty. Our method shows promise, significantly outperforming MEDIC, a commonly used benchmark, by both MAE and log-likelihood. In addition, STGP is more consistent, with less variability in accuracy. This work could be extended in a number of directions. Exogenous datasets relating to factors which influence ambulance demand (historic road traffic accidents, population density etc) may be used to increase prediction accuracy within a multi-task GPR framework (such as that described by [4]). The geography in Cape Town is challenging due to an irregular coastline and mountainous park areas which are largely uninhabited, causing a jagged spatial distribution of demand. This issue is addressed for Melbourne in [40] by using kernel-warping -projecting the predicted spatial distribution onto the distribution of sampled historic call outs. This could be augmented with the road networks and/or samples from the population density distribution. Such kernel warping would likely increase the accuracy of STGP. Methods to exploit the predictions might also be considered. For example, the predictive distribution, marginalised over space, is a forecast for total demand at that timestep. Percentiles of this distribution could be translated into possible staffing levels, with different levels of provision for different levels of certainty. Conversely, the predictions could be used as an indication of environmental dynamics in a planning problem, in which the resting places of ambulances between calls are optimised at different times of day.