Spatial Modelling of Emergency Service Response Times

This article concerns the statistical modelling of emergency service response times. We apply advanced methods from spatial survival analysis to deliver inference for data collected by the London Fire Brigade on response times to reported dwelling fires. Existing approaches to the analysis of these data have been mainly descriptive; we describe and demonstrate the advantages of a more sophisticated approach. Our final parametric proportional hazards model includes harmonic regression terms to describe how response time varies with time-of-day and shared spatially correlated frailties on an auxiliary grid for computational efficiency. We investigate the short-term impact of fire station closures in 2014. Whilst the London Fire Brigade are working hard to keep response times down, our findings suggest there is a limit to what can be achieved logistically: the present article identifies areas around the now closed Belsize, Downham, Kingsland, Knightsbridge, Silvertown, Southwark, Wesminster and Woolwich fire stations in which there should perhaps be some concern as to the provision of fire services.


Introduction
The thought of a fire in the home terrifies most people. In 2013-2014 UK fire and rescue services attended over half a million calls; there were a total of 322 fire-related deaths and 9748 non-fatal casualties due to fires of which 80% occurred in dwellings. There were 39600 dwelling fires in the UK in 2013-2014, with most occurring between the hours of 8 and 9pm at night and with misuse of equipment/appliances being the leading cause of around 1/3 1 arXiv:1503.07709v1 [stat.AP] 26 Mar 2015 these incidents (Department for Communities and Local Government, 2014). In London, 70% of fire-related deaths have been attributed to reporting delays, but crucial in saving lives is the efficient response of emergency services to calls to 999/112 (London Fire Brigade, 2013b).
The choice of where to locate emergency service stations (police, fire, ambulance) in cities has a direct impact on possible response times. From an academic perspective, this can be treated as an optimisation problem: how do we balance the need to respond quickly to emergency situations given a finite resource allocation (Toregas et al., 1971;Kolesar and Blum, 1973)?
Whilst such approaches are very useful in helping to decide where might be best to build emergency service stations in the first place, the urban environment is constantly changing and there is therefore a need to continually monitor and improve information on response times in order to ensure safety standards are maintained. One aspect of this is the profiling and mapping of high risk groups, which has been undertaken in a limited way in the UK, one example being a pair of studies in Merseyside (Higgins et al., 2013(Higgins et al., , 2014. Another aspect is the study of response times to emergency calls, the subject of the present article. The analysis of emergency response times has received a modest amount of attention in the literature. Scott et al. (1978) is one exception to this, these authors sought to form a mathematical model for ambulance response times in Houston. Other statistical approaches have focussed on predicting demand for emergency services such as Matteson et al. (2011), Vile et al. (2012) and Zhou et al. (2014). Recent concern over the UK Government's cuts to public services and their potential impact on the ability of fire services to maintain safety standards has resulted in a resurgence of interest, albeit primarily from the media and opposition parties (London Fire Brigade, 2013a,d;Open Data Institute, 2013;Bannister, 2014;Caven, 2014;Foley, 2014;Johnson, 2014;Westminster's Labour Councillors, 2014;Read, 2014;McCartney, 2015). The year 2014 saw 10 fire stations close: Belsize, Bow, Clerkenwell, Downham, Kingsland, Knightbridge, Silvertown, Southwark, Westminster and Woolwich.
In the present article our goal is to form a model for response times with the aim of providing emergency services with probabilistic information on where in space response times could be improved; clearly, faster response to fires in the home means more saved lives. Whilst in this article we focus on urban fires, it is worth noting the marked difference in response times to fires in urban compared with rural areas (Claridge and Spearpoint, 2013;Torney, 2013). In rural areas response times are usually longer and the results presented in the present article corroborate these findings: response times in the outskirts of the city, where there are fewer fire stations and there is generally more open space, are typically longer than near the centre of the city. The London Fire Brigade (LFB) is one of the largest fire and rescue services in the world; they collect and analyse substantial amounts of data on incident response times (London Fire Brigade, 2014a). In addition to presenting tables of average response time by ward, there has also been some investment in the mapping of response times at this level of aggregation (London Fire Brigade, 2013c, 2014b and in preparation for the proposed closures of 2014, an assessment of the potential impact of the closures was also carried out (Open Data Institute, 2013). Whilst these analytical efforts are to be highly commended, they could be further improved by the use of formal statistical models.
The main purpose of this article is to demonstrate ways in which the modelling-oriented approach could help to improve in the description and presentation emergency response time data and consequently inform city planners in their decision-making processes. To achieve this aim, we apply recently-developed techniques in survival analysis specifically designed for the modelling of large spatially referenced time-to-event data like the LFB response times (Taylor, 2015).
The remainder of this article is organised as follows. In Section 2 we introduce some basic concepts in spatial survival analysis and give details of the modelling approach proposed; in Section 3 we present results from the analysis of the LFB data; and in Section 4, the article concludes with a discussion.

Data and Model
The statistical analysis of time-to-event data is the realm of survival analysis (Cox and Oakes, 1984;Klein et al., 2013). Most often, survival methods are applied in clinical studies assessing the potential effect of treatments or exposures on the survival time of patients. Due to patients dropping out and the fact that studies are finite in duration, survival data are typically 'censored', which means the event of interest was not necessarily observed for all individuals. Survival methods handle the time to observed and censored events in a formal way.
Let T be a random variable, denoting the time after the call to 999/112 that the first fire engine arrives at the scene of a dwelling fire. We will shortly introduce a model for the hazard function defined formally by, ∆t .
For any time t, the interpretation of h(t) in this case is as the instantaneous arrival rate of the first engine at the scene of the fire conditional on the engine having not arrived before time t. The survival function, which is of particular interest in our current context can be derived from the hazard function, and represents the probability that the engine will arrive on the scene after time t.
Our ideal proportional hazards spatial survival model for the response times postulates the following form for the hazard function for the ith call: where X i are covariates associated with the ith call, β are parameters, h 0 is the baseline hazard function (see below) and Y i is the value a spatially continuous Gaussian process, Y , at the location of the ith call. We assume that the Gaussian process Y has associated covariance function, where the parameter σ 2 is the unconditional variance of the process at any point and ρ describes the interdependence of the points i and j at distance d ij apart.
In this model h 0 describes the part of the hazard function that is common to all individuals and the remaining term, exp{X i β + Y i }, describes the relative risk for the ith call. The relative risk splits into two parts: the first, exp{X i β}, is the part of the risk that we can explain by the available covariates; and the second, exp{Y i } is the unexplained risk. With regards the latter, we choose E[Y ] = −σ 2 /2, so that E[exp{Y i }] = 1. When survival models are used to measure time to death, a high hazard is regarded as bad, since it means there is a high chance of an individual experiencing the event. In the present context however, we need to adopt the opposite meaning in which a low hazard is interpreted as bad: at time t if no engine has yet arrived and h(t) is low, the chance of one arriving in the immediate future is low.
The two main options for modelling the baseline hazard are: (i) to assume a parametric form for h 0 ; or (ii) to leave h 0 unspecified, which results in a semiparametric model. Whilst the main advantage of the semiparametric approach is flexibility, in this article we opt for the former of these modelling paradigms because we are interested in probabilistic prediction.
We considered two different parametric models for the hazard function for these data. The first is a simple Weibull model where, so that the baseline cumulative hazard takes the form: In the second parametric model we mimic the flexibility of a semiparametric approach by modelling the baseline hazard using B-splines as in Rosenburg (1995), setting: where ω 1 , . . . , ω d are parameters to be estimated and B i (t) is a B-spline basis function, a piecewise positive polynomial of degree d, see Younes and Lachin (1997) for details on how to construct these. Being piecewise polynomial, the baseline cumulative hazard function H 0 (t) = t 0 h 0 (s)ds, required in likelihood computation, is trivial to compute provided we store the (piecewise) coefficients of the integrated basis functions: Our final model is a slight modification of (1) introduced in Taylor (2015); this modification concerns the frailties Y i . Rather than assume these are individual-specific we adopt a shared frailty approach, introducing an auxiliary grid of cells on which we wish to predict the response times. Our model for the hazard takes the form: where Y G[i] denotes the value of the process Y in the cell containing observation i: we approximate the spatially continuous process Y by a piecewise constant process on a fine grid. The reason for doing this is primarily computational efficiency: model (1)  For a discussion of other potential models for these data, see Section 4.

An Analysis of the 2012-2014 Data
The data on LFB response times analysed in this article is available from the London Datastore (London Fire and Emergency Planning Authority, 2015). All analyses were carried out in the R statistical software (R Core Team, 2014). In this section, we present an analysis of the LFB data from 2012 to 2014 inclusive. Table 1 shows the number of call-outs to fires, false alarms and special services by property type for the years 2012-2014 inclusive. We restrict our attention to analysing data from call-outs to fire events in dwellings because as mentioned above, dwellings are where most fire-related deaths occur. Figure   It can also be seen from Figure 1 that the intensity of locations of fire stations is also more concentrated towards the centre of the city, a bivariate isotropic Gaussian smoothing of these points is shown in Figure 2 (the intensity was scaled by a factor of 10 8 ). The bandwidth used to compute the smooth intensity was chosen using the rule of thumb method of Baddeley and

Preliminary Analyses and Model Choice
Turner's density.ppp function from the spatstat package, resulting in a kernel standard deviation of 5737; the result is quite a smooth approximation to the intensity of the fire stations which is desirable for the use we put it to. The relationship between response time to proximity of the nearest fire station is not clear cut: it is not necessarily the closest station that will respond to a call and for this reason, we use the smoothed intensity as a covariate in our model. We expect the coefficient of this covariate to be positive, since in places where there are a greater concentration of fire stations, fires are likely to be attended more swiftly.
Whilst the plot of response time by time of day, shown in Figure 3 suggests no obvious trends, it is reasonable to assume that time of day does in fact influence response time due to traffic congestion and local demand on the fire service. The demand by time of day is illustrated in Figure 4, this shows that it is lowest at around 5:30am, highest around 7pm and generally higher in the day-time than at night. There was no obvious difference in response times comparing weekdays with weekends: the median response time in 2014 for weekdays was 299 seconds and for weekends it was 296.5; again, though there is no obvious difference between these values, it is reasonable to assume that we might expect a difference in response times due in part to differing traffic patterns on weekdays compared with weekends.
Our overall model for the hazard function therefore includes the intensity of fire stations and a time trend through harmonic regression terms: sin(2πkt/24) and cos(2πkt/24) with k = 1, 2, 3, 4 i.e. with periodicity 1 day, 1/2 day, 1/3 day and 1/4 day. The number of harmonic terms was chosen by fitting a non-spatial version of the model to the 2010 data using forward selection, terminating the process when both the sin and cos contributions were not significant at the 5% level. We initially fitted the spatial Weibull model which included a weekday/weekend indicator variable as an additional covariate, but since this was not found to be significant, we fitted the B-spline model without this term. For the Bspline model, the baseline hazard function, h 0 was modelled using a piecewise cubic B-spline function of the form (3) with 5 internal knots at the minimum, 0.2, 0.4 and 0.6 quantiles and at the maximum response time; with a repeated knot at each end-point the spline has a total of seven parameters. We used an exponential covariance function for the spatial random effects, that is where d ij is the distance between the centroid of the cell representing Y i and that representing Y j .

Description of Inferential Procedure
Model (4) was fit separately to the 2012, 2013 and 2014 datasets using the R package spatsurv (Taylor and Rowlingson, 2014). This package implements a fully Bayesian adaptive MCMC algorithm, which delivers inference for the model parameters β and ω, the shared frailties Y G[i] and the parameters of the process Y (σ and φ). We used Gaussian priors for all parameters on an appropriately transformed scale: π(β) ∼ N(0, 100 2 ), log ω ∼ N(0, 10 2 ), log α ∼ N(0, 10 2 ), log λ ∼ N(0, 10 2 ), log σ ∼ N(0, 0.5 2 ) and log φ ∼ N(log 1000, 0.5 2 ) . As per the method described in Taylor (2015), we used a N(0, 1) prior for a whitened version of the spatial process, Γ where Y = −σ 2 /2 + Σ 1/2 σ,φ Γ, where Σ σ,φ is the covariance matrix of Y on the auxiliary grid; we retained the transformed samples We chose the size of cells in the computational grid to be 500m×500m.  and was run for 600,000 iterations with a 110,000 iteration burnin, again retaining every 490th sample. Plots of autocorrelation in the Y chains (for our final Weibull model) at lags 1, 5 and 10 are shown in Figure 11, these confirm that the the chain was mixing well: the autocorrelation in all cases had dropped to a negligible amount on or before the 5th lag. The B-spline chain mixed more slowly but by lag 10 autocorrelation in the 2013 and 2014 chains was low; the 2012 chain was mixing a little more slowly. For model comparison purposes we considered these chains to be sufficiently well mixing to decide between the Weibull and Bspline models. Diagnostic plots for the other chains in our final Weibull model are available by following the links appearing in Section 5. This is a challenging sampling problem: with around 6,000 observations and for technical reasons 256 × 256 = 65, 536 prediction points (reducing to an output grid of size 128 × 128) each chain takes around 4 days to run.
We chose between the Weibull and B-spline models using the DIC, shown in Table 3. It is interesting to note that the DIC from the simpler Weibull model is lower than that for the

Summary of Results
The result of fitting this model using the is a sample from the joint posterior density of all model parameters, π(β, ω, Y, η|data), from which we can compute expectations of quantities of interest; Figure 5 shows several such quantities.
Beginning with the baseline hazard, we report here the shape of the baseline hazard and Whilst the baseline hazard and cumulative hazard describe global properties of the process generating these data, the arrival rate of fire engines does depend on time of day and space as will now be illustrated. The bottom left plot in Figure 5 shows the relative risk by time of day: we do not present the coefficients of the harmonic regression terms here, as a plot is much simpler to interpret. This plot shows that there are two main times of day when services take longer than usual to arrive on the scene of a fire: between 3am and 7am and between 11am and 6pm the relative risk is significantly below 1 and reaches its lowest value of around 0.7 at 5am. It is interesting to note that this pattern was not observed in Figure 3, in which arrival times appear to be independent of the time of day; adjusting for the spatial variation in risk space has therefore bought out this trend.
The bottom right plot in Figure 5 shows the posterior covariance function with 95% confidence interval. This shows that spatial correlation is over quite a short range, around 0-1000 metres. Plots comparing the prior to the posterior for the parameters σ and φ showed that these were well identified by the data (the identifiability of φ is a common problem in spatial analyses). Table 2 gives the estimated coefficients from the Weibull model for the three years under consideration, it can be seen that the coefficients are quite similar for each year.
Using the spatial survival modelling framework, we can also illustrate answers to questions of substantive interest including (i) where in space is the London Fire Brigade's target response time of 6 minutes not being met; and (ii) what have been the effects on target response times of the 2014 fire station closures?
We answer the first question using the expected probability that the response time to a call will take longer than 6 minutes. For a response time in cell i of the computational grid, this is evaluated as, is the survival function in cell i evaluated for the jth retained sample of each parameter in the model. In computing S in the above, we assumed that the fire occurred on a weekday at 8:30pm, since as stated above, most dwelling fires occur between 8 and 9pm. The resulting plot is shown in Figure 6, where we have masked the computational grid cells appearing outside of the London boroughs. This plot shows that whilst there has been an improvement over the last three years in responding to reported fires in the outskirts of the city, there are some areas in the inner part of the city in 2014 where the expected probability the response will take longer than 6 minutes is slightly elevated compared with the other years -the colour has shifted from dark blue to light blue/neutral and red in some small areas, the largest stretching from below Westminster to around Camden Town and also of note the area around Gallions Point Mariana. Some caution must be maintained in not over-interpreting these maps, as these are point estimates, they are subject to uncertainty.
To account for this uncertainty and to make a comparison between years and thus address the second substantive question, in Figure 7, we plotted P[S y 1 (360) > S y 2 (360) + c] for y 1 , y 2 ∈ {2014, 2013, 2012} with y 1 > y 2 and c ∈ {0.1, 0.25}. Since we have samples from π(β, ω, η, Y |data) in each year. In the case that y 1 = 2014, y 2 = 2012 and c = 0.1 these probabilities were computed as: 1 1000 1000 j=1 I [S 2014 (360; β (j) , ω (j) , η (j) , Y (j) S 2012 (360; β (p(j)) , ω (p(j)) , η (p(j)) , Y (p(j)) where p(j) denotes a permuted index of the sample; in computing S we again assumed that  We used these probabilities to identify small 500×500m squares near the closed fire stations that are of potential concern in terms of response time. We identified regions as those in these small areas. The right hand plot is the same but illustrates times for all small regions inside the area of interest; there were are total of 8276 of these. Calls within the red regions of potential concern therefore accounted for 9% of all calls in the area of interest surrounding the closed stations. It can be seen from the two box and whisker plots that whilst in the region as a whole, the fire brigade appears to mainly be meeting their 6 minute target (right plot), in these small areas of potential concern near to the closed fire stations, there is a definite increase in response times in 2014; the median in these areas is above the 6 minute target.

Discussion
In this article we have shown how advanced methods from spatial survival analysis can be used to model emergency service response times. We have applied these methods to the London Fire Brigade data and have illustrated the impact of the 2014 closures on response times. We have identified areas of potential concern surrounding the recently closed stations; in these areas the median response time exceeds the London Fire Brigade's target of six minutes.
Whilst there may be simpler ways to model these data such as kriging the log-response times, the fitting of a spatial survival model is advantageous, being a 'natural' model for these time-to-event data. Secondly, the hazard, cumulative hazard and survival functions are useful for describing properties of the data that are of genuine interest in this context.
The proposed inferential framework is advantageous as it provides a sample from the joint posterior of all model parameters including the spatial process Y on all prediction cells.
These samples can be used to deliver posterior expectations of functionals of interest, some of which have been illustrated in this article. The main drawback with the proposed method is computational cost. We estimate it would take over 5 months to run the full model in Equation (1), so the 4 days the sampler takes represents a substantial reduction in cost.
Other techniques to speed up the MCMC such as using a Gaussian Markov random field to represent the spatial process on the auxiliary grid are also not as fast as the proposed  Figure 10: Plot of log posterior evaluated at the initial value and over all retained iterations of the chain. This shows that initially the chain started far away from a mode, but had found one before the burnin had finished.