Spatial point process via regularisation modelling of ambulance call risk

Turning to the dataset at hand, $W$ will represent the Swedish municipality of Skellefteå, and the time period under consideration is given by the years 2014–2018. Our dataset consists of 14,919 events where the mark structure is such that for each event, and thereby the person associated with the ambulance call in question, there is a recording of the event’s priority label (1 indicates the highest severity, implying turned on sirens, and 2 is a lower severity level) and a recording of the sex/gender of the person (female or male). Among the 14,919 events, 7,204 and 7,715 of them have priority labels 1 and 2, respectively. Unfortunately, a missing data issue is present here: 5,236 and 5,238 of the events are recorded as male and female, respectively, and the remaining events do not have any sex recording for the person related to the event, so they are labelled as "missing". Since our interest in the sex label lies mostly in highlighting possible structural differences related to the different covariates at hand, we have decided to proceed by studying the data as two separate marked point patterns: one where the marks are given by the priority labels and one where the marks are given by the sex labels. It should be acknowledged that there are other ways of dealing with this issue.

In order to properly model the ambulance call risk, we also need a range of spatial covariates. To begin with, a closer look at Figure 1 and the relation between the road network and the call locations justifies the inclusion of road network related covariates. In addition, demographic spatial covariates should also play a role here, given that different demographic zones have different behavioural patterns. The demographic spatial covariates considered have been supplied by Statistics Sweden (SCB), and the road network related covariates considered have been provided by the Swedish Transport Administration (Trafikverket).

Here, we will distinguish between two categories of covariates. The first is the collection of ’original’ covariates, which are raw covariates retrieved from SOS alarm, Statistics Sweden (SCB), and Trafikverket. These covariates will, in turn, be used to generate a collection of new covariates, which we will refer to as ’created’ covariates. In Section 4, we outline the construction of the created covariates in detail.

All covariates considered have been based on the year 2018 and in their final forms they are given as functions of the form $z_{i}(s)$, $s\in W$. Moreover, in addition to including each individual covariate $z_{i}$ in the analysis, we also include each interaction term $z_{ij}(s)=z_{i}(s)z_{j}(s)$, $s\in W$, which makes sense because we suspect that many of the covariates interact in one way or another. We have a total of 43 individual spatial covariates and, all in all, we consider a total of 989 $(={43\choose 2}+2\times 43)$ covariates as candidates to be included in the modelling of the ambulance call risk.

In application areas such as environmental science, epidemiology, ecology, etc., aside from the spatial locations of the events, additional information about the events may be available, which can be associated to the locations of the events. Such information pieces are referred to as marks, and by including them in the analysis, we can often obtain more realistic spatial point process models for the events – note that, in contrast, a covariate reflects information which is known throughout the spatial domain before realisation of the events. For instance, as we saw above, it is both common and practically appropriate in emergency medical services to document the priority levels, the gender, the incident time, etc. for a call/patient; recall that by spatial points of the events we mean the set of spatial locations of the patients to which ambulances have been dispatched or the set of spatial locations of patients from which calls have been made to the dispatcher/emergency alarm center.

As previously indicated, marked point patterns are modelled by marked point processes (baddeley2015spatial; moller2003statistical; diggle2013statistical; van2000markov). Given a spatial domain $W\subseteq\mathbb{R}^{2}$ and a mark space¹¹1Formally, $\mathcal{M}$ is assumed to be a Polish space. $\mathcal{M}$, a point process:

is a random subset such that, with probability 1, $Y\cap A\times B$ has finite cardinality $Y(A\times B)$ for any²²2Throughout, any set under consideration is a Borel set. $A\times B\subseteq W\times\mathcal{M}$; this is referred to as local finiteness. If we additionally have that $Y_{g}=\{x_{i}\}_{i=1}^{N}\subseteq W$ is a well-defined point process on $W$ (locally finite) in its own right, then we say that $Y$ is a marked point process. Note that each $(x_{i},m_{i})$ is a random variable and note in particular that if $W$ is bounded then $Y$ is automatically a marked point process.

When the mark space is discrete, e.g. $\mathcal{M}=\left\{j:j=1,2,\ldots,k,\;k>1\right\}$, we say that $Y$ is multitype and we note that we may split $Y$ into the marginal (purely spatial) point processes:

This collection may formally be represented by the vector $(Y_{1},Y_{2},\ldots,Y_{k})$, which is referred to be a multivariate point process, and one commonly uses the two notions interchangeably.

Our main interest here is to create a set of "heat maps" which describe the risk of a call occurring at a given spatial location $x\in W$, given an associated mark. This is accomplished by assuming that our data are generated by a multivariate point process $Y\subseteq W\times\{1,2,\ldots,k\}$, $k>1$ and then modelling the spatial intensity function of each of the component $Y_{j},j=1,2,\ldots,k.$ Formally, the spatial intensity function of $Y$ is defined as the function $\rho_{Y}$ satisfying³³3We assume that $\rho$ is the Radon-Nikodym derivative of the Borel measure $(A\times B)\mapsto E[Y(A\times B)]$ w.r.t. the product measure given by the product of a Lebesgue measure and a counting measure.

where $\mathbf{1}\{\cdot\}$ denotes the indicator function. Heuristically, $\rho_{Y}(x,j)dx$ may be interpreted as the probability that $Y$ has an event with mark $j$ in an infinitesimal neighbourhood of $x$ with size $dx$. Since $E[Y(A\times\{j\})]=E[Y_{j}(A)]$, it follows that the spatial intensity function of $Y_{j}$ satisfies

Hence, by modelling the marginal processes $\rho_{j}$ separately, we obtain a model for $\rho_{Y}$. Letting $X$ denote an arbitrary $Y_{j}$, $j=1,2,\ldots,k$, we see that its spatial intensity function $\rho$ satisfies

or equivalently, $\rho(x)=\lim_{\lvert dx\rvert\to 0}E\left[X\left(dx\right)\right]/\lvert dx\rvert.$ If a point process has constant intensity function then we say that it is homogeneous, otherwise we refer to it as inhomogeneous. In a similar fashion, we may define higher-order intensity functions $\rho_{n}(x_{1},x_{2},\ldots,x_{n})$, $x_{1},x_{2}\ldots,x_{n}\in W$, $n\geq 2$ as $\rho_{n}(x_{1},x_{2},\ldots,x_{n})=E[X(dx_{1})X(dx_{2})\cdots X(dx_{n})]=P(X(dx_{1})=1,X(dx_{2})=1,\ldots,X(dx_{n})=1)$ for disjoint infinitesimal neighbourhoods $dx_{1},dx_{2},\ldots,dx_{n}\subseteq W$ of $x_{1},x_{1},\ldots,x_{n}$. If $x_{i}=x_{j}$ for any $i\neq j$, then $\rho_{n}(x_{1},x_{2},\ldots,x_{n})=0$.

Poisson processes, for which $\rho_{n}(x_{1},x_{2},\ldots,x_{n})=\rho(x_{1})\rho(x_{2})\cdots\rho(x_{n})$, $n\geq 1$ and $X(A)$ is Poisson distributed with mean $\int_{A}\rho(x)dx$ for any $A\subseteq W$, are used as baseline models for the case of complete spatial randomness, i.e.,  the case where there is no spatial interaction/dependence present. Poisson processes are completely governed by their spatially varying intensity functions and may be viewed as generalisations of random samples to the case where the size of the sample is random.

We have observed a high degree of inhomogeneity in the ambulance call locations and we believe that a large portion of this inhomogeneity can be attributed to (some of) the spatial covariates which we have access to. Hence, as a starting point, for a given mark, we will consider parametric modelling of the call locations.

Consider a family of spatial intensity functions $\rho_{\boldsymbol{\theta}}(x)$, $x\in W$, which depends on spatial covariates through a parameter vector $\boldsymbol{\theta}\in\mathbb{R}^{K+1}$. A common and convenient approach when modelling the intensity function of an arbitrary point process is to proceed as if we are considering a Poisson process, which is commonly referred to as composite likelihood estimation and is often motivated by good large sample properties (coeurjolly2019understanding). Accordingly, suppose that a point pattern $\mathbf{x}=\left\{x_{1},x_{2},\cdots,x_{n}\right\}\subseteq W$ represents a realisation of a spatial point process $X$ which is observed within a bounded study region $W$. If $X$ is an inhomogeneous Poisson point process, then the associated log-likelihood function is given by

Using any quadrature rule, the integral in (1) can be approximated by a finite sum

where the positive numbers $w_{j}$, $j=1,2,\ldots,m$, are quadrature weights summing to the area $\lvert W\rvert$ and $s_{j}\in W$, $j=1,2,\ldots,m,$ are quadrature points. Following this approximation, the approximated log-likelihood function may be expressed as

where $\textbf{s}=\{s_{1},\ldots,s_{m}\}\subseteq W$ represents the union of the observed spatial locations of events $\mathbf{x}=\{x_{1},\ldots,x_{n}\}$ and the set of dummy points $\textbf{s}\setminus\mathbf{x}=\{v_{1},\ldots,v_{q}\}$, $q=m-n$. Here we assume that $m$ is much larger than $n$ for a better approximation of the log-likelihood function $\log\mathcal{L}\left(\boldsymbol{\theta}\mid\mathbf{x}\right)$. The approximated log-likelihood function in expression (3) can then be rewritten as follows:

where

Note that $a_{j}=0$ means that $s_{j}$ is a dummy point. Exploiting the approximation in equation (2), a large number of dummy points are required to obtain accurate parameter estimation using equation (4). waagepetersen2008estimating proposed two ways of obtaining dummy points and quadrature weights. With regard to the quadrature weights, the first method is a grid approach in which the observation window $W$ is partitioned into a collection of rectangular tiles. The quadrature weight for a quadrature point $s\in\textbf{s}$ falling in a tile $R$ is the area of $R$ divided by the total number of quadrature points falling in $R$. This approach is advantageous since the computation of the quadrature weights is easy. The second approach is the Dirichlet approach (okabe2009spatial) in which the quadrature weights are the areas of the tiles of the Dirichlet/Voronoi tessellation generated by the quadrature points in s. With regard to the dummy points, waagepetersen2008estimating proposed two ways of generating dummy points. The first approach is to use stratified dummy points combined with grid-type weights while the second approach is to exploit binomial dummy points with the Dirichlet-type weights.

According to baddeley2000practical and thurman2015regularized, a computationally cheaper approach to generate dummy points and compute quadrature weights is to partition the study region $W$ into tiles $R$ of equal area. To generate dummy points, we place one dummy point exactly in each tile either systematically or randomly. It follows that the quadrature weights for quadrature points $s_{j}$ can be set to $w_{j}=\Delta/E_{j}$, where $\Delta$ is the area of each tile and $E_{j}$ is the number of events and dummy points in the same tile as point $s_{j}$.

Modelling spatial intensity functions parametrically, in particular modelling based on spatial covariates, we often assume that $\rho_{\boldsymbol{\theta}}$ has a log-linear form. More specifically, letting $\boldsymbol{\beta}=\left(\beta_{1},\beta_{2},\ldots,\beta_{K}\right)$ and $\boldsymbol{\theta}=\left(\beta_{0},\boldsymbol{\beta}\right)$, we assume that

where $\mathbf{z}(x)=(z_{1}(x),z_{2}(x),\cdots,z_{K}(x))$ is a vector of spatial covariates at location $x\in W$. Combining equations (4) and (6), we thus find that the log-likelihood function $\log\mathcal{L}\left(\boldsymbol{\theta}\mid\mathbf{x}\right)$ can be approximated by

It follows that the expression on the right-hand side of the approximation sign in equation (7) is a weighted log-likelihood function of independent Poisson random variables $Y_{j},j=1,2,\cdots,m$. That is, for $j=1,2,\cdots,m$, $y_{j}$ are the observations, $\rho_{\boldsymbol{\theta}}(s_{j})=\exp\{\beta_{0}+\mathbf{z}(s_{j})\boldsymbol{\beta}^{\prime}\}$ are the intensities of the Poisson distributions, and $w_{j}$ are the weights. Thus, the weighted log-likelihood function in equation (7) can be maximised using standard software for fitting generalised linear models (McCullagh1989generalised).

Incorporating regularisation into the log-likelihood function in equation (7) can help to simultaneously select variables and estimate the parameters of the model. A penalised log-likelihood function based on equation (7) may be given by

where $R(\boldsymbol{\beta})$ is a regularisation method or penalty function, and $\lambda\geq 0$ is a tuning or smoothing parameter determining the strength of the penalty, or the amount of shrinkage. Several penalties such as garrote (breiman1995better), least absolute shrinkage and selection operator (lasso), elastic-net and fused lasso (tibshirani2005sparsity), group lasso (yuan2006model), Berhu penalty (owen2007robust), adaptive lasso (zou2006adaptive), and LAD-lasso (wang2007robust) have been developed for penalised regression modelling. The most common regularisation methods are lasso, ridge regression, elastic-net, and adaptive lasso. In short, lasso has a tendency to shrink several coefficients to zero, leaving only the most influential ones in the model, while ridge regression shrinks the coefficients of correlated covariates towards each other to borrow strength from each other (friedman2010regularization). The elastic-net penalty provides a mix between the ridge and the lasso penalty, and it is useful in cases where there are many correlated covariates or when the number of covariates exceeds the size of observations. The elastic-net regularisation penalty has the form

where the elastic-net parameter $\alpha\in[0,1]$ turns (9) into a ridge penalty if $\alpha=0$ and a lasso penalty if $\alpha=1$. If $\alpha=1-\epsilon$ for small $\epsilon>0$, then elastic-net performs like lasso but it avoids unstable behaviour due to extreme correlation (yue2015variable); empirical studies have indicated that elastic-net technique tends to outperform lasso on data with highly correlated features (comber2018geographically; friedman2010regularization). zou2006adaptive proposed an adaptive lasso to address the shortcomings of lasso, such as biased estimates of large coefficients and conflict between optimal prediction and consistent variable selection. According to kramer2009regularized, however, the performance of adaptive lasso is poor in the presence of highly correlated variables. Hence, our way forward here is to consider elastic-net penalisation and substituting equation (9) into equation (8), the optimization problem of the elastic-net penalization of the log-likelihood function can be summarized as

We propose the approach outlined below for the setting where i) the main goal is a prediction/predictive model, i.e., one wants to predict a collection of further/future events as precisely as possible, ii) one believes that the observed covariates can only describe a part of the spatial intensity variation, and iii) added spatial flexibility is warranted in the modelling. In the case of our ambulance data, as our end goal is to build optimal dispatching strategies, we mainly want a predictive model. We want to mention again that the demographic covariates we have at hand only reflect where different demographic groups live but not how they move about. Explicitly, we do not have access to explicit mobility covariates such as aggregated movement patterns in the population, and, as one may guess, people do not only need access to ambulances when they are at home.

Our solution is quite pragmatic and simple. We simply add a further spatial covariate to the existing collection of covariates, which is given by a non-parametric spatial intensity estimate $\widetilde{\rho}(x)$, $x\in W$ and will be referred to as the benchmark spatial intensity. Hence, we include the spatial intensity estimate $\widetilde{\rho}$ as a covariate in the approximated log-likelihood expression in (7) and the inclusion has the effect that the modelling steps away from a purely parametric setting to a semi-parametric approach. This added covariate should pick up on regions where there is an increased intensity due to human mobility, which cannot be explained by the existing list of covariates. Note that this is similar in nature to the semi-parametric (spatio-temporal) log-Gaussian Cox process modelling approach advocated for in e.g. diggle2013statistical. To be able to discern whether this added covariate is in fact necessary/useful in the presence of the other covariates, we carry out elastic-net regularisation-based variable selection (see Section 3.5) to indicate whether the benchmark spatial intensity has any added value in terms of describing the true intensity function.

A natural question here is what kind of non-parametric intensity estimator one should use to generate the benchmark spatial intensity. There are different candidates for this, and the main distinction one usually makes is between global and adaptive/local smoothers. Adaptive smoothing techniques include adaptive kernel intensity estimation (davies2018tutorial) and (resample-smoothed) Voronoi intensity estimation (ogata2003modelling; ogata2004space; moradi2019resample). These have some clear benefits (in particular the latter, moradi2019resample), but here we do not want to put too much weight on the local features since we may run the risk of overfitting. Instead, we here consider (global) kernel intensity estimation (diggle1985kernel; baddeley2015spatial), which is arguably the most prominent approach to global smoothing and is defined as

where $\kappa_{h}(\cdot)=h^{-1}\kappa(\cdot/h)$, $\kappa$ is a symmetric density function and the smoothing parameter $h>0$ is the bandwidth. The function $w_{h}(x,y)$ is a suitable edge correction factor which adjusts the effect of unobserved events outside $W$ on the intensity of the observed events (baddeley2015spatial); we here use the local corrector $w_{h}(x,y)=\int_{W}\kappa_{h}(u-y)du$ which ensures mass preservation, i.e. that $\int_{W}\widetilde{\rho}(x)dx=n$, the number of observed events. Practically, to carry out kernel intensity estimation we make use of the function density.ppp in the R package spatstat (baddeley2015spatial).

Although the choice of kernel may play a certain role, the choice of bandwidth is absolutely the main determinant for the quality of the intensity estimate; recall that the bandwidth governs how much we smooth the data. However, optimal bandwidth selection is a well-studied and challenging problem. Concerning the state of the art, the bandwidth selection criterion of cronie2018non is generally the most stable with respect to accounting for spatial interaction; observed clusters of points in a point pattern may be the effect of aggregation/clustering (dependence) or intensity peaks, or a combination of the two. However, there is one scenario where it tends to not perform too well, and that is when there are large regions in $W$ where there are no points present, which is the case for our ambulance dataset. Other standard methods for bandwidth selection include the state estimation approach of diggle1985kernel (called bw.diggle in spatstat), the Poisson process likelihood leave-one-out cross-validation approach in baddeley2015spatial and loader1999local (called bw.ppl in spatstat), and the recent machine learning-based approach of bayisa2020large (see Algorithm 2).

The role of the road network-related covariates is to control that the fitted model generates events close to the underlying road network. But, as previously indicated, the original road network-related covariates are not of the form $z_{i}(s)$, $s\in W$.

• •

  Complete road networks line density. It is a spatial pattern of line segments, which is converted to a pixel image. The value of each pixel in the image is measured as the total length of intersection between the pixel and the line segments [1].

• •

  Main road networks line density. It is a spatial pattern of line segments, which is converted to a pixel image. The value of each pixel in the image is measured as the total length of intersection between the pixel and the line segments [1].

• •

  Densely populated line density. It is a spatial pattern of line segments, which is converted to a pixel image. The value of each pixel in the image is measured as the total length of intersection between the pixel and the line segments [1].

• •

  Bus stops density. It is a spatial point pattern, which is converted to a pixel image. The value of each pixel is an intensity, which is measured as "points per unit area" [1].

• •

  Shortest distance to bus stops. It is the shortest distance in meters from the ambulance location data to the bus stops [1].

• •

  Shortest distance to densely populated areas. It is the shortest distance in meters from the ambulance location data to the densely populated areas [1].

• •

  Shortest distance to main road networks. It is the shortest distance in meters from the ambulance location data to the main road networks [1].

• •

  Shortest distance to complete road networks. It is the shortest distance in meters from the ambulance location data to the complete road networks [1].

We here propose to treat the road networks under consideration as line segment patterns (baddeley2015spatial), which essentially means that each road network considered is a realisation of a point process in the space of line segments in $\mathbb{R}^{2}$. The spatial covariate corresponding to a given line segment pattern is then given by the estimated line segment intensity, which is obtained as the convolution of an isotropic Gaussian kernel with the line segments of the pattern in question. Practically, such an estimate may be obtained through the function density.psp in the R package spatstat, and the standard deviation of the Gaussian kernel, the bandwidth, determines the degree of smoothing. The default bandwidth choice in density.psp is given by the diameter of the observation window multiplied by 0.1. As an alternative, we propose to use our new heuristic algorithm for bandwidth selection, which is achieved by letting $W=[0,\infty)$, and letting the observations $\mathbf{x}=\left\{x_{1},x_{2},\cdots,x_{n}\right\}\subseteq W$ considered in Algorithm 2 represent the lengths of the line segments. Figure 4 compares the default bandwidth choice of density.psp to our heuristic algorithm, and it clearly suggests that the line segment intensities generated using the heuristic algorithm bandwidth selection have captured the spatial pattern of the line segments in the observed data better than the default bandwidth choice. Note that what interests us here are the relative scales rather than the raw scales; the values in the plots have been multiplied by 1000 for ease of visualisation.

Figure 5 shows the estimated coefficients of the spatial covariates, using both lasso regularisation and lasso-like elastic-net regularisation, for the spatial point pattern constituting the events with priority label 1. Note that the numbers at the top of each panel in Figure 5 indicate the number of spatial covariates with non-zero coefficients, i.e. the number of covariates that are associated with the fitted spatial intensity functions for the indicated value of $\lambda$. We clearly see how a large number of covariates are quickly excluded as we increase $\log\lambda$.

The two vertical lines in the figure have been drawn to show the location of the logarithm of the optimal estimate of the tuning parameter (blue) and the location of the logarithm of the estimate of the tuning parameter that is one standard error away from the optimal estimate (red); the one-standard-error-rule (hastie2017elements) says that one should go with the simplest model, which is no more than one standard error worse than the best model.

The aforementioned analyses have utilised the spatial covariates $z_{i}(s)$, $s\in W$, and the interaction terms of the form $z_{i}(s)z_{j}(s)$, which holds for the cases $i=j$ and $i\neq j$. However, with this approach, it is not easy to interpret the results and present the estimated models. Therefore, we are interested in further exploring the estimation of the intensity of the emergency alarm call events using only the original spatial covariates $z_{i}(s)$, $s\in W$, i.e., no first-order and higher-order interactions of spatial covariates in the modelling of the spatial intensity. Modelling the spatial intensity of the emergency alarm call events using only the original spatial covariates can help to identify the spatial covariates that play a key role in determining the spatial distribution of the emergency alarm call events, and it also ease the interpretation and presentation of the results. Moreover, we can also compare the estimated intensities that are obtained by the two approaches (i.e., modelling with only the original spatial covariates and modelling with the original spatial covariates and the first-order interaction terms) to discern the approach that can be practicable. Using only the original spatial covariates, a lasso-like elastic-net ($\alpha=0.95$) with an optimal tuning parameter estimate does not provide sparse solutions, i.e., most of the coefficients of the spatial covariates are non-zero. That is to say, the lasso-like elastic-net provides dense solutions, i.e., it has a low rate of variable exclusion, see Table LABEL:tab:EstimatedModels1.