Local inhomogeneous second-order characteristics for spatio-temporal point processes occurring on linear networks

Abstract

Point processes on linear networks are increasingly being considered to analyse events occurring on particular network-based structures. In this paper, we extend Local Indicators of Spatio-Temporal Association (LISTA) functions to the non-Euclidean space of linear networks, allowing to obtain information on how events relate to nearby events. In particular, we propose the local version of two inhomogeneous second-order statistics for spatio-temporal point processes on linear networks, the K- and the pair correlation functions. We put particular emphasis on the local K-functions, deriving come theoretical results which enable us to show that these LISTA functions are useful for diagnostics of models specified on networks, and can be helpful to assess the goodness-of-fit of different spatio-temporal models fitted to point patterns occurring on linear networks. Our methods do not rely on any particular model assumption on the data, and thus they can be applied for whatever is the underlying model of the process. We finally present a real data analysis of traffic accidents in Medellin (Colombia).

1 Introduction

Point processes on linear networks are increasingly being considered to analyse events occurring on particular network structures, such as traffic accidents, street crimes or ambulance calls and interventions. Such locations inherently live on the corresponding network structure, and considering such geometry as the support of data results in defining a more realistic scenario. Nevertheless, geometrical complexities of linear networks give rise to different mathematical and computational challenges. Point processes on linear networks were first introduced in the spatial context and then extended to the spatio-temporal case, focusing on the analysis of first- and second-order summary statistics (Ang et al. 2012; McSwiggan et al. 2017; Rakshit et al. 2017; Moradi et al. 2019; Rakshit et al. 2019; Moradi and Mateu 2020; D’Angelo et al. 2021). However, none of these approaches have considered local properties of the introduced second-order statistics to measure local structure and to further provide local diagnostics to fitted point process models on network-based structures. In general, local extensions of spatio-temporal statistics on networks are very welcome in many scientific fields, such as epidemiology, criminology, or sociology, where one could be interested in identifying those events that most contributed to the fitting of the model while accounting for the geometry of the underlying network.

D’Angelo et al. (2021) have extend local indicators of spatio-temporal association (Siino et al. 2018), known as LISTA functions, for spatio-temporal point processes occurring on linear networks. As the proposed local second-order statistics can be used for obtaining further insight into the local structure of the analysed point pattern, and on the characteristics of individual points, D’Angelo et al. (2021) have used the LISTA homogeneous versions to build a local test allowing to assess local differences between the spatio-temporal second-order structure of two point patterns occurring on the same linear network. This paper aims at developing inhomogeneous versions of two local second-order statistics for spatio-temporal point processes occurring on networks, namely the K-function and the pair correlation function (Moradi and Mateu 2020). In particular, following Adelfio et al. (2020) for the Euclidean case, we use LISTA functions based on local inhomogeneous second-order statistics on networks to assess the goodness-of-fit of general spatio-temporal models. Indeed, the peculiar lack of homogeneity in a network structure discourages the use of traditional spatial and spatio-temporal methods based on stationary processes (Baddeley et al. 2021). Therefore, weighted second-order statistics are appropriate diagnostic tools since they directly apply to data without assuming homogeneity (Adelfio et al. 2020).

We provide several simulation studies by weighting the contribution of each observed point by the inverse of the intensity function coming from several types of inhomogeneous point processes on the network. In particular, we consider Poisson, self-exciting and log-Gaussian spatio-temporal processes, and detail the corresponding simulation procedures on a linear network.

We finally provide an application of our methodology to traffic data. All the analyses are carried out through the software R Core Team (2020) and the codes are currently available upon request to the authors, although they will be published soon in a dedicated package.

The paper is structured as follows. Section 2 is devoted to the introduction of basic definitions of spatio-temporal point processes, and reports the LISTA functions for the Euclidean case. Section 3 reviews some basics about the analysis of spatio-temporal point patterns occurring on linear networks, focusing on their second-order characteristics. These first two sections are also aimed at defining the notation used throughout the paper. In Sect. 4 the proposed LISTA functions are defined and, through some simulation studies, in Sect. 5 we show how these can be used as diagnostic tools for different fitted models. Further, an application to a case study is presented in Sect. 6. Conclusions are drawn in Sect. 7.

2 Second-order properties of spatio-temporal point processes on Euclidean spaces

We consider a spatio-temporal point process with no multiple points as a random countable subset X of $\mathbb {R}^2 \times \mathbb {R}$, where a point $({\textbf {u}}, t) \in X$ corresponds to an event at ${\textbf {u}} \in \mathbb {R}^2$ occurring at time $t \in \mathbb {R}$. A typical realisation of a spatio-temporal point process X on $\mathbb {R}^2 \times \mathbb {R}$ is a finite set $\{({\textbf {u}}_i, t_i)\}^n_{ i=1}$ of distinct points within a bounded spatio-temporal region $W \times T \subset \mathbb {R}^2 \times \mathbb {R}$, with area $\vert W\vert > 0$ and length $\vert T\vert > 0$, where $n \ge 0$ is not fixed in advance. In this context, $N(A \times B)$ denotes the number of points of a set $(A \times B) \cap X$, where $A \subseteq W$ and $B \subseteq T$. As usual (Daley and Vere-Jones 2007), when $N(W \times T) < \infty$ with probability 1, which holds e.g. if X is defined on a bounded set, we call X a finite spatio-temporal point process.

For a given event $({\textbf {u}}, t)$, the events that are close to $({\textbf {u}}, t)$ in both space and time, for each spatial distance r and time lag h, are given by the corresponding spatio-temporal cylindrical neighbourhood of the event $({\textbf {u}}, t)$, which can be expressed by the Cartesian product as

$$\begin{aligned} b(({\textbf {u}}, t), r, h) = \{({\textbf {v}}, s) : \vert \vert {\textbf {u}} - {\textbf {v}}\vert \vert \le r, \vert t - s \vert \le h\} , \quad \quad ({\textbf {u}}, t), ({\textbf {v}}, s) \in W \times T, \end{aligned}$$

where $\vert \vert \cdot \vert \vert$ denotes the Euclidean distance in $\mathbb {R}^2$. Note that $b(({\textbf {u}}, t), r, h)$ is a cylinder with centre (u, t), radius r, and height 2h.

Product densities $\lambda ^{(k)}, k \in \mathbb {N} \text { and } k \ge 1$, arguably the main tools in the statistical analysis of point processes, may be defined through the so-called Campbell Theorem (see Daley and Vere-Jones 2007), which states that, given a spatio-temporal point process X, for any non-negative function f on $( \mathbb {R}^2 \times \mathbb {R} )^k$

$$\begin{aligned} \mathbb {E} \Bigg [ \sum _{\zeta _1,\dots ,\zeta _k \in X}^{\ne } f( \zeta _1,\dots ,\zeta _k)\Bigg ]=\int _{\mathbb {R}^2 \times \mathbb {R}} \dots \int _{\mathbb {R}^2 \times \mathbb {R}} f(\zeta _1,\dots ,\zeta _k) \lambda ^{(k)} (\zeta _1,\dots ,\zeta _k) \prod _{i=1}^{k}\text {d}\zeta _i, \end{aligned}$$

that constitutes an essential result in spatio-temporal point process theory. In particular, for $k=1$ and $k=2$, these functions are respectively called the intensity function $\lambda$ and the (second-order) product density $\lambda ^{(2)}$. Broadly speaking, the intensity function describes the rate at which the events occur in the given spatio-temporal region, while the second-order product densities are used when the interest is in describing spatio-temporal variability and correlations between pair of points of a pattern. They represent the point process analogues of the mean function and the covariance function of a real-valued process, respectively. Then, the first-order intensity function is defined as

$$\begin{aligned} \lambda ({\textbf {u}},t)=\lim _{\vert \text {d}{} {\textbf {u}} \times \text {d}t\vert \rightarrow 0} \frac{\mathbb {E}[N(\text {d}{} {\textbf {u}} \times \text {d}t )]}{\vert \text {d}{} {\textbf {u}} \times \text {d}t\vert }, \end{aligned}$$

where $\text {d}{} {\textbf {u}} \times \text {d}t$ defines a small region around the point $({\textbf {u}},t)$ and $\vert \text {d}{} {\textbf {u}} \times \text {d}t\vert$ is its volume. The second-order intensity function is given by

$$\begin{aligned} \lambda ^{(2)}(({\textbf {u}},t),({\textbf {v}},s))=\lim _{\vert \text {d}{} {\textbf {u}} \times \text {d}t\vert ,\vert \text {d}{} {\textbf {v}} \times \text {d}s\vert \rightarrow 0} \frac{\mathbb {E}[N(\text {d}{} {\textbf {u}} \times \text {d}t )N(\text {d}{} {\textbf {v}} \times \text {d}s )]}{\vert \text {d}{} {\textbf {u}} \times \text {d}t\vert \vert \text {d}{} {\textbf {v}} \times \text {d}s\vert }. \end{aligned}$$

Finally, the pair correlation function

$$\begin{aligned} g(({\textbf {u}},t),({\textbf {v}},s))=\frac{ \lambda ^{(2)}(({\textbf {u}},t),({\textbf {v}},s))}{\lambda ({\textbf {u}},t)\lambda ({\textbf {v}},s)} \end{aligned}$$

can be interpreted formally as the standardised probability density that an event occurs in each of two small volumes, $\text {d}{} {\textbf {u}} \times \text {d}t$ and $\text {d}{} {\textbf {v}} \times \text {d}s$, in the sense that for a Poisson process, $g(({\textbf {u}},t),({\textbf {v}},s))=1$.

In this paper, the focus is on second-order characteristics of spatio-temporal point patterns, with an emphasis on the K-function (Ripley 1976). This is a measure of the distribution of the inter-point distances and captures the spatio-temporal dependence of a point process. Gabriel and Diggle (2009) extend the second-order methods provided by Baddeley et al. (2000) to the spatio-temporal setting, defining the spatio-temporal inhomogeneous K-function and proposing a nonparametric estimator.

Definition 1

Gabriel and Diggle (2009) A spatio-temporal point process is second-order intensity reweighted stationary and isotropic if its intensity function is bounded away from zero and its pair correlation function depends only on the spatio-temporal difference vector (r, h), where $r= \vert \vert {\textbf {u}}-{\textbf {v}} \vert \vert$ and $h= \vert t-s \vert$.

Definition 2

Gabriel and Diggle (2009) For a second-order intensity reweighted stationary, isotropic spatio-temporal point process, the space-time inhomogeneous K-function takes the form

$$\begin{aligned} K(r,h)=2 \pi \int _{0}^{r} \int _0^{h} g(r',h')r'\text {d}r'\text {d}h' \end{aligned}$$

(1)

where $g(r,h)=\lambda ^{(2)}(r,h)/(\lambda ({\textbf {u}},t)\lambda ({\textbf {v}},s)), r=\vert \vert {\textbf {u}}-{\textbf {v}}\vert \vert ,h= \vert t-s \vert$.

The simplest expression of an estimator of the spatio-temporal K-function is given as

$$\begin{aligned} \hat{K}(r,h)=\frac{1}{ \vert W \vert \vert T \vert }\sum _{i=1}^n \sum _{j > i} I( \vert \vert {\textbf {u}}_i-{\textbf {u}}_j \vert \vert \le r, \vert t_i-t_j \vert \le h). \end{aligned}$$

(2)

For a homogeneous Poisson process $\mathbb {E}[\hat{K}(r,h)]=\pi r^2 h$, regardless of the intensity $\lambda$.

Both the K-function and the pair correlation function can be used as a measure of spatio-temporal clustering and interaction (Gabriel and Diggle 2009; Møller and Ghorbani 2012). Usually, $\hat{K}(r,h)$ is compared with the theoretical $\mathbb {E}[\hat{K}(r,h)]=\pi r^2 h$. Values $\hat{K}(r,h) > \pi r^{2} h$ suggest clustering, while $\hat{K}(r,h) < \pi r^2 h$ points to a regular pattern.

Inhomogeneous second-order statistics can be constructed and used for assessing goodness-of-fit of fitted first-order intensities. Nevertheless, it is widespread practice in the statistical analysis of spatial and spatio-temporal point pattern data to focus primarily on comparing the data with a homogeneous Poisson process, which is generally the null model in applications, rather than the fitted model. Indeed, when dealing with diagnostics in point processes, often two steps are needed: the transformation of data into residuals (thinning or rescaling (Schoenberg 2003)) and the use of tests to assess the consistency of the residuals with the homogeneous Poisson process (Adelfio and Schoenberg 2009). Usually, second-order statistics estimated for the residual process (i.e. the result of a thinning or rescaling procedure) are analysed. Essentially, to each observed point a weight inversely proportional to the conditional intensity at that point is given. This method was adopted by Veen and Schoenberg (2006) in constructing a weighted version of the K-function of Ripley and Kelly (1977); the resulting weighted statistic is in many cases more powerful than residual methods (Veen and Schoenberg 2006).

The spatio-temporal inhomogeneous version of the K-function in (2) is given by Gabriel and Diggle (2009) as

$$\begin{aligned} \hat{K}_{I}(r,h)=\frac{ \vert W \vert \vert T \vert }{n(n-1)}\sum _{i=1}^n \sum _{j > i} \frac{I( \vert \vert {\textbf {u}}_i-{\textbf {u}}_j \vert \vert \le r,\vert t_i-t_j\vert \le h)}{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {u}}_j,t_j)}, \end{aligned}$$

(3)

where $\lambda (\cdot ,\cdot )$ is the first-order intensity at an arbitrary point. We know that $\mathbb {E}[\hat{K}_{I}(r,h)]=\pi r^2 h$, that is the same as the expectation of $\hat{K}(r,h)$ in (2), when the intensity used for the weighting is the true generator model. This is a crucial result that allows to use the weighted estimator $\hat{K}_{I}(r,h)$ as a diagnostic tool, for assessing the goodness-of-fit of spatio-temporal point processes with generic first-order intensity functions. Indeed, if the weighting intensity function is close to the true one $\lambda ({\textbf {u}},t)$, the expectation of $\hat{K}_{I}(r,h)$ should be close to $\mathbb {E}[\hat{K}(r,h)]=\pi r^2 h$ for the Poisson process. For instance, values $\hat{K}_{I}(r,h)$ greater than $\pi r^{2} h$ indicates that the fitted model is not appropriate, since the distances computed among points exceed the Poisson theoretical ones.

2.1 Local Indicators of Spatio-Temporal Association functions

Definition 3

Siino et al. (2018) Local Indicators of Spatio-Temporal Association (LISTA) are a set of functions that are individually associated with each one of the points of the point pattern, and can provide information about the local behaviour of the pattern.

This operational definition of local indicators was introduced by Anselin (1995) for the spatial case, and extended by Siino et al. (2018) to the spatio-temporal context. If $\lambda ^{(2)i}(\cdot ,\cdot )$ denotes the local version of the spatio-temporal product density for the event $({\textbf {u}}_i,t_i)$, then, for fixed r and h, it holds that

$$\begin{aligned} \hat{\lambda }^{(2)}_{\epsilon ,\delta }(r,h)=\frac{1}{n-1}\sum _{i=1}^n\hat{\lambda }^{(2)i}_{\epsilon ,\delta }(r,h), \end{aligned}$$

(4)

where $\hat{\lambda }^{(2)i}_{\epsilon ,\delta }(r,h)=\frac{n-1}{4\pi r \vert W \times T \vert }\sum _{j\ne i}\kappa _{\epsilon ,\delta }( \vert \vert {\textbf {u}}_i-{\textbf {v}}_j \vert \vert -r, \vert t_i-s_j \vert -h),$ with $r>\epsilon >0$ and $h>\delta >0$, and $\kappa$ a kernel function with spatial and temporal bandwidths $\epsilon$ and $\delta$, respectively.

Definition 4

Siino et al. (2018) Any second-order spatio-temporal summary statistic that satisfies the operational definition in (4), which means that the sum of spatio-temporal local indicator functions is proportional to the global statistic, can be called a LISTA statistic.

In Adelfio et al. (2020), local versions of both the homogeneous and inhomogeneous spatio-temporal K-functions on the Euclidean space are introduced. Defining an estimator of the overall intensity by $\hat{\lambda }=n/(\vert W \vert \vert T \vert )$, they propose the local version of (2) for the i-th event $({\textbf {u}}_i,t_i)$ as

$$\begin{aligned} \hat{K}^i(r,h)=\frac{1}{\hat{\lambda }^2 \vert W \vert \vert T \vert }\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} I( \vert \vert {\textbf {u}}_i-{\textbf {v}} \vert \vert \le r,\vert t_i-s\vert \le h) \end{aligned}$$

(5)

and the local version of (3) as

$$\begin{aligned} \hat{K}^i_{I}(r,h)=\frac{1}{ \vert W \vert \vert T \vert }\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{I(||{\textbf {u}}_i-{\textbf {v}} \vert \vert \le r,\vert t_i-s\vert \le h)}{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {v}},s)}, \end{aligned}$$

(6)

with $({\textbf {v}},s)$ being the spatial and temporal coordinates of any other point. The authors extended the spatial weighting approach of Veen and Schoenberg (2006) to spatio-temporal local second-order statistics, proving that the inhomogeneous second-order statistics behave as the corresponding homogeneous ones, basically proving that the expectation of both (5) and (6) is equal to $\pi r^2 h$. In our paper, we follow the same reasoning, extending the results to the setting of spatio-temporal point processes occurring on particular network structures. Therefore, next section is devoted to the review of some basics about point processes occurring on linear networks.

3 Point processes on linear networks and their second-order characteristics

Point processes on linear networks are recently considered to analyse events occurring on particular network structures such as the traffic accidents on a road network that we analyse in this paper (see Fig. 13). Spatial patterns of points along a network of lines are indeed found in many applications. The network might reflect a map of railways, rivers, electrical wires, nerve fibres, airline routes, irrigation canals, geological faults or soil cracks (Baddeley et al. 2021). Observations of interest could be the locations of traffic accidents, bicycle incidents, vehicle thefts or street crimes, and many others.

A linear network $L=\cup _{i=1}^{n} l_{i} \subset \mathbb {R}^{2}$ is commonly taken as a finite union of line segments $l_i\subset \mathbb {R}^{2}$ of positive length. A line segment is defined as $l_i=[u_i,v_i]=\{ku_i+(1-k)v_i: 0 \le k \le 1\}$, where $u_i,v_i \in \mathbb {R}^2$ are the endpoints of $l_i$. For any $i \ne j$, the intersection of $l_i$ and $l_j$ is either empty or an endpoint of both segments. A spatio-temporal linear network point process is a point process on the product space $L \times T$, where L is a linear network and T is a subset (interval) of $\mathbb {R}$.

We hereafter focus on a spatio-temporal point process X on a linear network L with no overlapping points $({\textbf {u}},t)$, where ${\textbf {u}} \in L$ is the location of an event and $t \in T (T \subseteq \mathbb {R}^+)$ is the corresponding time occurrence of ${\textbf {u}}$. Note that the temporal state-space T might be either a continuous or a discrete set. A realisation of X with n points is represented by ${\textbf {x}} = {({\textbf {u}}_i ,t_{i} ),i = 1,\dots ,n}$ where $({\textbf {u}}_i ,t_{i} ) \in L \times T$.

A spatio-temporal disc with centre $({\textbf {u}},t) \in L \times T$, network radius $r > 0$ and temporal radius $h > 0$ is defined as

$$\begin{aligned} b(({\textbf {u}},t ),r,h) = \{({\textbf {v}},s ) : d_L ({\textbf {u}},{\textbf {v}}) \le r , \vert t - s \vert \le h\},\quad \quad ({\textbf {u}}, t), ({\textbf {v}}, s) \in L \times T \end{aligned}$$

where $\vert \cdot \vert$ is a numerical distance, and $d_L(\cdot ,\cdot )$ stands for the appropriate distance in the network, typically taken as the shortest-path distance between any two points. The cardinality of any subset $A \subseteq L \times T, N(X \cap A) \in {0,1,\dots }$, is the number of points of X restricted to A, whose expected value is denoted by

$$\begin{aligned} \nu (A) = \mathbb {E}[N(X \cap A)],\quad A \subseteq L \times T, \end{aligned}$$

where $\nu$, the intensity measure of X, is a locally finite product measure on $L\times T$ (Baddeley et al. 2006).

We now recall Campbell’s theorem for point processes on linear networks (Cronie et al. 2020). Assuming that the product densities/intensity functions $\lambda ^{(k)}$ exist, for any non-negative measurable function $f(\cdot )$ on the product space $L^k$, we have

$$\begin{aligned} \mathbb {E} \Bigg [ \sum _{\zeta _1,\dots ,\zeta _k \in X}^{\ne } f( \zeta _1,\dots ,\zeta _k)\Bigg ]=\int _{L^k} f(\zeta _1,\dots ,\zeta _k) \lambda ^{(k)} (\zeta _1,\dots ,\zeta _k) \prod _{i=1}^{k}\text {d}\zeta _i, \end{aligned}$$

(7)

where $\ne$ indicates that the sum is over distinct values. Assume that X has an intensity function $\lambda (\cdot ,\cdot )$, hence Eq. (7) reduces to

$$\begin{aligned} \mathbb {E}[N(X \cap A)] =\int _{A} \nu (d({\textbf {u}},t )) = \int _{A} \lambda ({\textbf {u}},t)d_2({\textbf {u}},t), A \subseteq L \times T, \end{aligned}$$

where $d_2 ({\textbf {u}},t)$ corresponds to integration over $L \times T$.

The second-order Campbell’s theorem is obtained from (7) with $k=2$

$$\begin{aligned} \begin{aligned}&\mathbb {E} \Bigg [ \sum _{({\textbf {u}},t),({\textbf {v}},s)\in X}^{\ne } f\big (({\textbf {u}},t),({\textbf {v}},s)\big ) \Bigg ] \\&\quad = \int _{L \times T} \int _{L \times T} f\big (({\textbf {u}},t),({\textbf {v}},s)\big ) \lambda ^{(2)}\big (({\textbf {u}},t),({\textbf {v}},s)\big )\text {d}_2({\textbf {u}},t)\text {d}_2({\textbf {v}},s). \end{aligned} \end{aligned}$$

(8)

Assuming that X has a second-order product density function $\lambda ^{(2)} (\cdot ,\cdot )$, we then obtain

$$\begin{aligned} \mathbb {E}[N(X \cap A)N(X \cap B)] = \int _{A} \int _{B} \lambda ^{(2)} (({\textbf {u}},t ),({\textbf {v}},s ))d_2 ({\textbf {u}},t)d_2 ({\textbf {v}},s ), \quad A,B \subseteq L \times T. \end{aligned}$$

Finally, an important result concerns the conversion of the integration over $L \times T$ to that over $\mathbb {R} \times \mathbb {R}$ (Rakshit et al. 2017). For any measurable function $f: L \times T \rightarrow \mathbb {R}$

$$\begin{aligned} \int _{L \times T} f({\textbf {u}},t)\text {d}_2({\textbf {u}},t)=\int _0^{\infty } \int _0^{\infty } \sum _{\begin{array}{c} ({\textbf {u}},t)\in L \times T:\\ d_L({\textbf {v}}, {\textbf {u}})=r,\\ \vert s-t \vert =h \end{array}} f({\textbf {u}},t) \text {d}r\text {d}h. \end{aligned}$$

(9)

Letting $f({\textbf {u}},t) = \eta (d_L({\textbf {u}},{\textbf {v}}), \vert t-s\vert )$ then

$$\begin{aligned} \int _{L \times T} \eta (d_L({\textbf {u}},{\textbf {v}}), \vert t-s\vert ) \text {d}_2({\textbf {u}},t)= \int _0^{\infty } \int _0^{\infty } \eta (r,h)M(({\textbf {u}},t),r,h)\text {d}r \text {d}h \end{aligned}$$

where $M(({\textbf {u}},t),r,h)$ is the number of points lying exactly at the shortest-path distance $r \ge 0$ and the time distance $h \ge 0$ away from $({\textbf {u}},t)$.

3.1 Global second-order characteristics for spatio-temporal point processes on networks

Following the work of Moradi and Mateu (2020), we first recall the concept of second-order pseudostationary point processes, and then we turn to the inhomogeneous case.

Definition 5

Moradi and Mateu (2020) Assume X is a spatio-temporal point process on $L \times T$ with constant intensity function $\lambda \ge 0$. Let the homogeneous K-function be given by

$$\begin{aligned} K_L(({\textbf {u}},t),r,h)=\frac{1}{\lambda }\mathbb {E} \Bigg [ \sum _{({\textbf {u}},t)\in X} \frac{I(d_L({\textbf {u}},{\textbf {v}})\le r,\vert t-s\vert \le h)}{M(({\textbf {u}},t),d_L({\textbf {u}},{\textbf {v}}),\vert t-s\vert )} \Bigg \vert ({\textbf {u}},t \in X) \Bigg ]. \end{aligned}$$

Then, X is called second-order pseudostationary and isotropic if $K_L(({\textbf {u}},t),r,h)$ does not depend on the point $({\textbf {u}},t)$, and we then write $K_L(({\textbf {u}},t),r,h)=K_L(r,h)$.

Definition 6

Moradi and Mateu (2020) For a homogeneous Poisson point process on $L \times T$ with constant intensity function $\lambda$, $K_L(r,h)=rh$.

The authors proposed a nonparametric estimator of both the homogeneous K-function

$$\begin{aligned} \hat{K}_L(r,h)=\frac{ \vert L \vert \vert T \vert }{n(n-1)}\sum _{i=1}^n\sum _{i \ne j} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {u}}_j)<r,\vert t_i-t_j\vert <h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {u}}_j),\vert t_i-t_j\vert )}, \end{aligned}$$

and pair correlation function

$$\begin{aligned} \hat{g}_L(r,h)=\frac{ \vert L \vert \vert T \vert }{n(n-1)}\sum _{i=1}^n\sum _{i \ne j} \frac{\kappa ( d_L({\textbf {u}}_i,{\textbf {u}}_j)-r)\kappa (\vert t_i-t_j\vert -h) }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {u}}_j),\vert t_i-t_j\vert )}, \end{aligned}$$

where $\vert L \vert >0$ and $\vert T \vert >0$ are the total length of the network L and of the time interval T, respectively, and $\kappa$ is a one-dimensional kernel function.

As in the Euclidean case, for a Poisson process on the linear network L, these second-order statistics can be used to indicate clustering and inhibition, if compared to their Poisson process values rh. Both the global homogeneous K-function and pair correlation functions are used in Moradi and Mateu (2020) for detecting departure from homogeneity of spatio-temporal point patterns occurring on different linear networks.

Motivated by practical situations where homogeneity is not a realistic assumption, Moradi and Mateu (2020) further proposed the corresponding inhomogeneous second-order characteristics, namely the inhomogeneous K-function

$$\begin{aligned} \hat{K}_{L,I}(r,h)=\frac{1}{ \vert L \vert \vert T \vert }\sum _{i=1}^n\sum _{i \ne j} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {u}}_j)<r,\vert t_i-t_j\vert <h\} }{\hat{\lambda }({\textbf {u}},t)\hat{\lambda }({\textbf {v}},s)M(({\textbf {u}},t),d_L({\textbf {u}},{\textbf {v}}),\vert t-s\vert )}, \end{aligned}$$

(10)

and the inhomogeneous pair correlation function

$$\begin{aligned} \hat{g}_{L,I}(r,h)=\frac{1}{ \vert L \vert \vert T \vert }\sum _{i=1}^n\sum _{i \ne j} \frac{\kappa ( d_L({\textbf {u}}_i,{\textbf {u}}_j)-r)\kappa (\vert t_i-t_j\vert -h) }{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {u}}_j,t_j)M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {u}}_j),\vert t_i-t_j\vert )}, \end{aligned}$$

(11)

where $\hat{\lambda }(\cdot ,\cdot )$ is an estimate of the intensity function, and $M(({\textbf {u}},t),r,h)$ is the number of points lying exactly at the shortest-path distance $r \ge 0$ and the time distance $h \ge 0$ away from the point $({\textbf {u}},t)$.

Moradi and Mateu (2020) recommend using some normalisation factor of the form $D(X)=\frac{1}{( \vert L \vert \vert T \vert )^2}\sum _{i=1}^n \sum _{i \ne j} \frac{1}{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {u}}_j,t_j)}$ in such a way that the updated estimates of the inhomogeneous K- and pair correlation functions $\frac{1}{D({\textbf {x}})}\hat{K}(r,h)$ and $\frac{1}{D({\textbf {x}})}\hat{g}(r,h)$ provide estimators with low bias and variance. They also prove that for Poisson processes, and for any $r,h>0$, it holds that $\hat{K}_{L,I}(r,h)=rh$ and $\hat{g}_{L,I}(r,h)=1$, which is useful in model selection and hypothesis testing.

Therefore, a relevant result is that the expectation of $\hat{K}_{L,I}(r,h)$, when the intensity used for the weighting is the true generator model, is the same as the expectation of $\hat{K}_{L}(r,h)$ for the Poisson process, that is $\mathbb {E}[\hat{K}_{L}(r,h)]=rh$. Therefore, as in the Euclidean case, the weighted estimator $\hat{K}_{L,I}(r,h)$ can be used as a diagnostic tool, for assessing the goodness-of-fit of spatio-temporal point processes occurring on a linear network with generic first-order intensity functions. Consequently, if the weighting intensity function is close to the true one $\lambda ({\textbf {u}},t)$, the estimated value $\hat{K}_{L,I}(r,h)$ should be close to $\mathbb {E}[\hat{K}_{L}(r,h)]= r h$.

In D’Angelo et al. (2022a), the inhomogeneous global K-function in (10) is used as a diagnostic tool for comparing the goodness-of-fit of different spatio-temporal parametric models fitted to a point pattern of visitors’ stops by touristic attractions in Palermo (Italy).

4 Local Indicators of Spatio-Temporal Association on linear networks

In this section, we propose the local versions of the previously reviewed global summary statistics for spatio-temporal point processes occurring on a network.

Therefore, we propose the local spatio-temporal inhomogeneous K-function for the i-th event $({\textbf {u}}_i,t_i)$ on a linear network as

$$\begin{aligned} \hat{K}^i_{L,I}(r,h)=\frac{1}{ \vert L \vert \vert T \vert }\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert <h\} }{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {v}},s)M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )}, \end{aligned}$$

(12)

and the corresponding local pair correlation function

$$\begin{aligned} \hat{g}^i_{L,I}(r,h)=\frac{1}{ \vert L \vert \vert T \vert }\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{\kappa ( d_L({\textbf {u}}_i,{\textbf {v}})-r)\kappa (\vert t_i-s\vert -h) }{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {v}},s)M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )}, \end{aligned}$$

(13)

weighted by the reciprocal of the normalising factor

$$\begin{aligned} D(X) = \frac{n-1}{ \vert L \vert \vert T \vert }\sum _{i=1}^n\sum _{i \ne j} \frac{1}{\hat{\lambda }({\textbf {u}}_i,t_i)\hat{\lambda }({\textbf {u}}_j,t_j)}, \end{aligned}$$

(14)

obtaining $\frac{1}{D(X)}\hat{K}^i_{L,I}(r,h)$ and $\frac{1}{D(X)}\hat{g}^i_{L,I}(r,h)$. Basically, we avoid summing up all the points as in the global statistics counterparts, and we denote the individual contribution to the global statistics with the index i.

The homogeneous versions are computed by weighting the second-order summary statistics by the constant intensity $\hat{\lambda }=n/( \vert L \vert \vert T \vert )$, obtaining

$$\begin{aligned} \hat{K}_L^i(r,h)=\frac{1}{\hat{\lambda }^{2} \vert L \vert \vert T \vert }\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert <h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )}, \end{aligned}$$

(15)

and

$$\begin{aligned} \hat{g}_L^i(r,h)=\frac{1}{\hat{\lambda }^{2} \vert L \vert \vert T \vert }\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{\kappa ( d_L({\textbf {u}}_i,{\textbf {v}})-r)\kappa (\vert t_i-s\vert -h) }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )}. \end{aligned}$$

(16)

Following the Euclidean case, D’Angelo et al. (2021) provide the operational definition based on the local second-order spatio-temporal summary statistics on a linear network.

Definition 7

D’Angelo et al. (2021) Any local second-order spatio-temporal summary statistic $\lambda ^{i}_{2,L}(({\textbf {u}},t),({\textbf {v}},s))$ computed on an observed point pattern ${\textbf {x}} = {({\textbf {u}}_i ,t_{i} ),i = 1,\dots ,n}$ on a linear network L, where $({\textbf {u}}_i ,t_{i} ) \in L \times T$, that satisfies the following operational definition

$$\begin{aligned} \hat{\lambda }_{2,L}(({\textbf {u}},t),({\textbf {v}},s))=\frac{1}{n-1}\sum _{({\textbf {u}}_i,t_i)\in {\textbf {x}}}\hat{\lambda }^{i}_{2,L}(({\textbf {u}}_i,t_i),({\textbf {v}},s)), \end{aligned}$$

(17)

where $\hat{\lambda }_{2,L}$ is an estimator of the global second-order intensity function, can be called a LISTA function on a linear network.

Knowing that the sum of the individual contributions given by the local estimator of the K-function in Eq. (15) are actually equal to the global statistics, Eq. (17) is satisfied, and therefore it can be called LISTA function on a linear network.

Throughout the paper we will only focus on the K-function, but of course also the pair correlation function can serve as an alternative.

Theorem 1

For a Poisson process, we have $\mathbb {E}[\hat{K}_L^i(r,h) ]= rh$.

Proof

$$\begin{aligned} \mathbb {E}[\hat{K}_L^i(r,h) ]&= \frac{1}{{\lambda }^{2}\vert L\vert \vert T\vert }\mathbb {E}\Bigg [\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \Bigg ] \\&\qquad \mathrm{(By~theorem~8~and~since}\, \lambda ^{(2)} \equiv \lambda ^2\, \mathrm{for~a~Poisson~process)} \\&= \frac{1}{{\lambda }^{2}\vert L\vert \vert T\vert } \int _{L \times T} \int _{L \times T} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} {\lambda }^{2} \text {d}_2({\textbf {v}},s) \\ \\&= \frac{1}{{\lambda }^{2}\vert L\vert \vert T\vert } {\lambda }^{2} \int _{L \times T} \int _{L \times T} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \text {d}_2({\textbf {v}},s) \\ \\&= \frac{1}{\vert L\vert \vert T\vert } \int _{L \times T} \int _{L \times T} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \text {d}_2({\textbf {v}},s) \\&\qquad \mathrm{(By~the~change~of~variables~in~9)} \\&=\frac{1}{\vert L\vert \vert T\vert } \int _0^{\infty } \int _0^{\infty } \sum _{\begin{array}{c} ({\textbf {v}},s)\in L \times T:\\ d_L({\textbf {u}}_i,{\textbf {v}})=r,\\ \vert t_i-s\vert =h \end{array}} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \text {d}r\text {d}h \\&=\frac{1}{\vert L\vert \vert T\vert } \int _0^{\infty } \int _0^{\infty }\sum _{\begin{array}{c} ({\textbf {v}},s)\in L \times T:\\ d_L({\textbf {u}}_i,{\textbf {v}})=r,\\ \vert t_i-s\vert =h \end{array}} \frac{I\{ (\vert \vert {\textbf {u}}_i-{\textbf {v}}\vert \vert )<r,\vert t_i-s\vert<h\} }{M({\textbf {u}}_i,t_i)} \text {d}r\text {d}h \\&=\frac{1}{\vert L\vert \vert T\vert } \int _0^{\infty } \int _0^{\infty } \frac{I\{ (\vert \vert {\textbf {u}}_i-{\textbf {v}}\vert \vert )<r,\vert t_i-s\vert<h\}M({\textbf {u}}_i,t_i) }{M({\textbf {u}}_i,t_i)} \text {d}r\text {d}h \\&=\frac{1}{\vert L\vert \vert T\vert } \int _0^{\infty } \int _0^{\infty } I\{ (\vert \vert {\textbf {u}}_i-{\textbf {v}}\vert \vert )<r,\vert t_i-s\vert <h\} \text {d}r\text {d}h\\&\qquad \mathrm{(taking~the~conditional~expectation~with~respect~to~its~filtration)} \\&=\frac{1}{\vert L\vert \vert T\vert } \vert L\vert \vert T\vert rh \\&=rh \end{aligned}$$

$\square$

Theorem 2

If the process is weighted by the true intensity function, the expectation of $\hat{K}_{L,I}^i(r,h)$ is the same as the expectation of $\hat{K}_L^i(r,h)$.

Proof

$$\begin{aligned}&\mathbb {E}[\hat{K}^i_{L,I}(r,h) ] \\&\quad = \frac{1}{\vert L\vert \vert T\vert }\mathbb {E}\Bigg [\sum _{({\textbf {u}}_i,t_i)\ne ({\textbf {v}},s)} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{{\lambda }({\textbf {u}}_i,t_i){\lambda }({\textbf {v}},s)M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \Bigg ] \\&\qquad \qquad \mathrm{(Second-order~Campbell's~theorem~in~8)}\\&\quad = \frac{1}{\vert L\vert \vert T\vert } \int _{L \times T} \int _{L \times T} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{{\lambda }({\textbf {u}}_i,t_i){\lambda }({\textbf {v}},s)M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} {\lambda }({\textbf {v}},s)\text {d}_2({\textbf {v}},s) \\&\quad = \frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} \int _{L \times T} \int _{L \times T} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )}\frac{1}{{\lambda }({\textbf {v}},s)} {\lambda }({\textbf {v}},s)\text {d}_2({\textbf {v}},s)\\&\quad = \frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} \int _{L \times T} \int _{L \times T} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \text {d}_2({\textbf {v}},s)\\&\qquad \qquad \mathrm{(Change~of~variables~in~9)}\\&\quad =\frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} \int _0^{\infty } \int _0^{\infty } \sum _{\begin{array}{c} ({\textbf {v}},s)\in L \times T:\\ d_L({\textbf {u}}_i,{\textbf {v}})=r,\\ \vert t_i-s\vert =h \end{array}} \frac{I\{ d_L({\textbf {u}}_i,{\textbf {v}})<r,\vert t_i-s\vert<h\} }{M(({\textbf {u}}_i,t_i),d_L({\textbf {u}}_i,{\textbf {v}}),\vert t_i-s\vert )} \text {d}r\text {d}h \\&\quad =\frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} \int _0^{\infty } \int _0^{\infty }\sum _{\begin{array}{c} ({\textbf {v}},s)\in L \times T:\\ d_L({\textbf {u}}_i,{\textbf {v}})=r,\\ \vert t_i-s\vert =h \end{array}} \frac{I\{ (\vert \vert {\textbf {u}}_i-{\textbf {v}}\vert \vert )<r,\vert t_i-s\vert<h\} }{M({\textbf {u}}_i,t_i)} \text {d}r\text {d}h\\&\quad =\frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} \int _0^{\infty } \int _0^{\infty } \frac{I\{ (\vert \vert {\textbf {u}}_i-{\textbf {v}}\vert \vert )<r,\vert t_i-s\vert<h\}M({\textbf {u}}_i,t_i) }{M({\textbf {u}}_i,t_i)} \text {d}r\text {d}h\\&\quad =\frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} \int _0^{\infty } \int _0^{\infty } I\{ (\vert \vert {\textbf {u}}_i-{\textbf {v}}\vert \vert )<r,\vert t_i-s\vert <h\} \text {d}r\text {d}h\\&\qquad \qquad \mathrm{(taking~the~conditional~expectation~with~respect~to~its~filtration)}\\&\quad =\frac{1}{\vert L\vert \vert T\vert } \frac{1}{{\lambda }({\textbf {u}}_i,t_i)} {\lambda }({\textbf {u}}_i,t_i) \vert L\vert \vert T\vert rh\\&\quad =rh \square \end{aligned}$$

We know that the expectation of $\hat{K}_{L,I}(r,h)$, when the intensity used for the weighting is the true generator model, is the same as the expectation of $\hat{K}_{L}(r,h)$ for the Poisson process, that is $\mathbb {E}[\hat{K}_{L}(r,h)]=rh$.

We have proved that the same result holds for the local version $\hat{K}_{L,I}^i(r,h)$, meaning that the our proposed local estimator $\hat{K}^i_{L,I}(r,h)$ for a general point process, weighted by the true intensity function, has the same expectation of the local estimator $\hat{K}_{L}^i(r,h)$ under the Poisson case.

Based on this result, we can use our proposed local estimator $\hat{K}_{L,I}^i(r,h)$ as a diagnostic tool for general spatio-temporal point processes occurring on a linear network. Since the local inhomogeneous estimator behaves as the corresponding homogeneous one of a Poisson process, we can follow the approach for diagnostics in a local scale in Adelfio et al. (2020), and use our proposed LISTA functions based on the K-function for assessing the goodness-of-fit of spatio-temporal point processes on linear networks with any generic first-order intensity function $\lambda (\cdot ,\cdot )$.

Basically, departures of the LISTA functions $\hat{K}_{L,I}^i(r,h)$ from the Poisson expected value rh directly suggest the unsuitability of the intensity function $\lambda (\cdot )$ used in the weighting of the LISTA functions. This means that, if the estimated intensity function used for weighting in our proposed LISTA functions (15) is the true one, then the LISTA functions should behave as the corresponding ones of a homogeneous Poisson process (12), that corresponds to the reference model.

5 Local space-time diagnostics on networks

This section is devoted to the use of the proposed LISTA functions to assess the goodness-of-fit of different spatio-temporal models fitted to point patterns occurring on linear networks.

Following Adelfio et al. (2020) for the Euclidean case, we show some simulation studies by generating different spatio-temporal point patterns, and then performing diagnostics on different fitted intensities. By comparing the values of the LISTA functions and their theoretical values, we evaluate whether the LISTA functions can correctly identify the true intensity, when this is constrained on a network.

Since in simulations the weights are obtained by considering the real intensity function, the inhomogeneous statistics are expected to behave as the ones of a homogeneous Poisson process. If departures from such behaviour were observed, that would be an indication that the data comes from a model identified by a first-order intensity function different from the one used in the weighting procedure.

5.1 Simulation set up

We simulate 100 space-time point patterns from 3 different inhomogeneous processes with 150 number of points on average: (a) a spatio-temporal inhomogeneous Poisson process on a linear network; (b) a spatio-temporal ETAS model; (c) a log-Gaussian Cox process.

For assessing the usefulness of the weighted LISTA functions also in network domains, we carry out a simulation study considering three different scenarios, in which the K-function is weighted by:

• The real intensity function, known in simulations

• The estimated intensity by a nonparametric kernel density function

• The intensity function estimated by an homogeneous Poisson intensity (wrong model).

We expect that the expected value of the spatio-temporal LISTA K-function, when weighting by the true or estimated intensity, resembles to that of a Poisson process, and therefore, when a wrong model is used for the intensity, the expected value would be far from the Poisson one.

5.2 Simulation schemes

As we simulate from different patterns on networks, we outline the simulation schemes used for simulating such processes.

5.2.1 Inhomogeneous spatio-temporal point process

To generate inhomogeneous spatio-temporal point processes on a linear network driven by the spatio-temporal intensity $\hat{\lambda }({\textbf {u}},t)$, we proceed with the following steps:

• 1. Define the generating intensity function $\lambda _0({\textbf {u}},t)=\hat{\lambda }({\textbf {u}},t)$.

• 2. Set an upper bound $\lambda _{max}$ for $\lambda _0({\textbf {u}},t)$.

• 3. Simulate a homogeneous Poisson process with intensity $\lambda _{max}$ and denote by N the number of the generated points, with coordinates $({\textbf {u}}', t')$.

• 4. Compute $p({\textbf {u}}',t')=\frac{\lambda ({\textbf {u}}',t')}{max(\lambda ({\textbf {u}}',t'))}$ for each point $({\textbf {u}}',t')$ of a homogeneous Poisson process.

• 5. Generate a sample ${\textbf {x}}'$ of size N from the uniform distribution on (0, 1).

• 6. Thin the simulated homogeneous Poisson process ${\textbf {x}}$ retaining the $n \le N$ locations for which ${\textbf {x}}' \le p$.

The algorithm provides a point pattern ${\textbf {x}}$ with n points. This scheme follows the procedure for simulating inhomogeneous point patterns outlined in Gabriel et al. (2013). In the context of point patterns occurring on linear networks, the events are guaranteed to be constrained on the network since each point $({\textbf {u}},t)$ of the generating intensity function $\hat{\lambda }({\textbf {u}},t)$ occurs only on the network line segments.

Fig. 1

Spatio-temporal inhomogeneous Poisson process with intensity function as in (18). Left panel: The projection of data onto the network. Right panel: The cumulative number of data points versus occurrence time

The simulation is performed through the function rpoistlpp() in the R package stlnpp (Moradi and Mateu 2020). In Fig. 1, a simulated inhomogeneous spatio-temporal Poisson process is depicted, with intensity

$$\begin{aligned} \lambda (x,y,t)=\exp (8.25 -4y -2t), \end{aligned}$$

(18)

decreasing together with the y-axis and with time.

5.2.2 Spatio-temporal Epidemic Type Aftershock Sequence point process

We now outline the generating scheme for simulating a pattern from an Epidemic Type Aftershocks-Sequence (ETAS) process (Ogata and Katsura 1988) with conditional intensity function as in Adelfio and Chiodi (2020) but adapting it on the network. We do this by defining the space-time region in a fixed temporal range and a linear network.

According to a branching structure, the conditional intensity function of the self-exciting model is defined as the sum of a term describing the large-time scale variation (spontaneous activity or background, generally assumed homogeneous in time but not in space) and one relative to the small-time scale variation due to the interaction with the events in the past (induced or triggered activity). The ETAS model can be written as

$$\begin{aligned} \lambda _{\theta }(t,{\textbf {u}}\vert \mathcal {H}_t)=\mu f({\textbf {u}})+\sum _{t_j<t} \frac{\kappa _0 \exp (\eta _j)}{(t-t_j+c)^p} \{ ({\textbf {u}}-{\textbf {u}}_j)^2+d \}^{-q} , \end{aligned}$$

(19)

where $\mathcal {H}_t$ is the past history of the process, $\mu$ is the background general intensity, and $f({\textbf {u}})$ is the spatial density. Concerning the triggered component, $\eta _j=\varvec{\beta }' {\textbf {Z}}_j$ is a linear predictor, with ${\textbf {Z}}_j$ the external known covariate vector, including the magnitude, and $\varvec{\theta }= (\mu , \kappa _0, c, p, d, q, \varvec{\beta })$ are the parameters to be estimated. In the usual ETAS model, with the only external covariate representing the magnitude, $\varvec{\beta }=\alpha$, as $\eta _j = \varvec{\beta }' {\textbf {Z}}_j = \alpha (m_j-m_0)$, where $m_j$ is the magnitude of the jth event and $m_0$ the threshold magnitude, that is, the lower bound for which earthquakes with higher values of magnitude are surely recorded in the catalogue.

Fig. 2

Spatio-temporal ETAS model with conditional intensity function as in (19). Left panel: The projection of data onto the network. Right panel: The cumulative number of data points versus occurrence time

The steps for the simulations are the following:

• 1. Input the true parameters values $\mu$, $\kappa _0$, c, p, d, q, the $\beta$ parameters related to the covariates. Also provide other control parameters such as the boundaries of the space-time region, and parameter b of the Gutenberg-Richter distribution.

• 2. Generate a random number $n_0$ from a Poisson distribution with parameter $E(N)=\mu \Delta t$.

• 3. Generate $n_0$ triples of space-time uniform coordinates in the space-time region.

• 4. For each point ${\textbf {x}}_j$, $j = 1, \dots , n_0$:

  • (a) Generate a random number $n_j$ from a Poisson distribution with parameter proportional to $exp(\eta _j)$.

  • (b) Generate $n_j$ triples of space-time coordinates of aftershocks in the space-time region.

  • (c) Add the $n_j$ new points to the set of events.

• 5. Proceed from step (2) until all the events inside the region are involved in the simulation process as possible generators of further events.

The algorithm provides a point pattern ${\textbf {x}}$ with n points. It also generates n magnitudes and other covariates’ values, if specified.

The simulation is performed through the function sim_ETASnet of the package LISTAnet. The simulation on the network is guaranteed by the homogeneous spatial Poisson process being generated on the network. Figure 2 shows a simulated ETAS process, with conditional intensity function as in (19) with parameters $\varvec{\theta }=(\mu , \kappa _0,c,p,\alpha ,d,q)=(0.079,0.004,0.013,1.2,0.5,0.424,1.165)$, displaying the typical epidemic structure, with very dense spatial and temporal clusters. We considered an ETAS process simulated with the only magnitude covariate, and therefore $\varvec{\beta }= \alpha$.

5.2.3 Spatio-temporal log-Gaussian Cox process

Following the inhomogeneous specification in Diggle (2013), a log-Gaussian Cox process for a generic point in space and time has the random intensity

$$\begin{aligned} \Lambda ({\textbf {u}},t)=\lambda ({\textbf {u}},t)\exp (S({\textbf {u}},t)) \end{aligned}$$

where S is a Gaussian process with $\mathbb {E}(S({\textbf {u}},t))=\mu =-0.5\sigma ^2$ and so $\mathbb {E}(\exp {S({\textbf {u}},t)})=1$ and with variance and covariance matrix $\mathbb {C}(S({\textbf {u}},t),S({\textbf {v}},s))=\sigma ^2 \gamma (r,h)$, with $\gamma (\cdot )$ the correlation function of the corresponding Gaussian random field (GRF).

For the generation of a spatio-temporal log-Gaussian Cox process on a linear network we proceed with the following algorithm:

• 1.a Generate a realisation from a GRF $S({\textbf {u}},t)$, with covariance function $\mathbb {C}(({\textbf {u}},t),({\textbf {v}},s))$ and mean function $\mu ({\textbf {u}},t)$.

• 1.b Define the generating intensity function $\lambda _0({\textbf {u}},t)=\hat{\lambda }({\textbf {u}},t)\exp (S({\textbf {u}},t))$.

• 2. Set an upper bound $\lambda _{max}$ for $\lambda _0({\textbf {u}},t)$.

• 3. Simulate a homogeneous Poisson process with intensity $\lambda _{max}$ and denote by N the number of the generated points, with coordinates $({\textbf {u}}', t')$.

• 4. Compute $p({\textbf {u}}',t')=\frac{\lambda ({\textbf {u}}',t')}{max(\lambda ({\textbf {u}}',t'))}$ for each point $({\textbf {u}}',t')$ of an homogeneous Poisson process.

• 5. Generate a sample ${\textbf {x}}'$ of size N from the uniform distribution on (0, 1).

• 6. Thin the simulated homogeneous Poisson process retaining the $n \le N$ locations for which ${\textbf {x}}' \le p$.

This algorithm provides a point pattern ${\textbf {x}}$ with n points. The Gaussian Random Field is generated using the function RFsim of the CompRandFld (Padoan and Bevilacqua 2015) package of R, and the overall simulation is performed thorugh the function sim_LGCPnet of the package LISTAnet. A simulated log-Gaussian Cox process is displayed in Fig. 3, with intensity function as in (18), but with discrete time $t=0,\dots , 4$. We consider a separable spatio-temporal exponential covariance with parameters $\sigma = 2$, $\phi = 0.2$ and $\theta = 1$. This process exhibits a lower but still present clustered structure, if compared to the ETAS process. The GRF is actually generated on the whole plane and then the results are constrained on the network since there is not a clear procedure for generating Gaussian Random Fields on networks. Therefore for the log-Gaussian Cox process we use Euclidean distances (Baddeley et al. 2021; D’Angelo et al. 2022a).

Fig. 3

Spatio-temporal log-Gaussian Cox process with parameters as in Sect. 5.2.3. Left panel: The projection of data onto the network. Right panel: The cumulative number of data points versus occurrence time

5.3 Results

For each of the 3 previously defined scenarios, we compute the chi-squared statistics between the LISTA functions and their theoretical values $\hat{K}^i_{L,Theo}(r,h)=rh$. In particular, we use the LISTA functions based on the inhomogeneous K-function, weighted each time with a different fitted intensity. The $\chi _i^2$ values are obtained following the expression

$$\begin{aligned} \chi _i^2=\int _L \int _T \Bigg ( \frac{\big (\hat{K}^i_{L,I}(r,h)-rh\big )^2}{rh} \Bigg ) \text {d}h \text {d}r , \end{aligned}$$

one for each point in the point pattern. We weight the LISTA functions $\hat{K}^i_{L,I}(r,h)$ with different intensities: the true one, a kernel density function, and a constant one. The latter of course is not a good fit, as we have simulated from inhomogeneous processes. In Fig. 4, as expected, the lowest values are found in correspondence of the true intensity and the non-parametric kernel-based estimation, while the largest values are encountered for the constant intensity. The same holds for the ETAS process, and for the log-Gaussian Cox process, in Figs. 5 and 6, respectively.

In summary, LISTA functions are able to correctly identify the true intensity and to spot the best fit among candidate models.

We now report an example of application on a ETAS process, simulated through the same procedure outlined in Sect. 5.2.2, to show further advantages of using LISTA functions for carrying out diagnostics.

Fig. 4

Distributions of the $\chi _i^2$ statistics computed for assessing the differences between the local LISTA functions on the network and the theoretical expected value, calculated for each individual point from a simulation of an inhomogeneous Poisson process with intensity as in (18)

Fig. 5

Distributions of the $\chi _i^2$ statistics computed for assessing the differences between the local LISTA functions on the network and the theoretical expected value, calculated for each individual point from a simulation of an ETAS process with intensity as in (19)

Fig. 6

Distributions of the $\chi _i^2$ statistics computed for assessing the differences between the local LISTA functions on the network and the theoretical expected value, calculated for each individual point from a simulation of an log-Gaussian Cox process process following the parameters in Sect. 5.2.3

5.4 Motivating example: a simulated ETAS process

To show some advantages of using LISTA functions, instead of their global counterparts, we simulate an ETAS process on the network (under the same conditions as in Fig. 5).

We then compute the LISTA functions weighting them with three different intensities: the true one, a fitted one, and a constant one. For assessing their goodness-of-fit, the standard approach (Baddeley et al. 2015) would require the computation of the global inhomogeneous second-order summary statistics in Sect. 3, weighted with their respective intensities. In Fig. 7 we represent the difference between the global K-functions and their theoretical values, weighted by the three considered intensities. The surfaces whose values are closest to zero are the ones whose intensity best resembles the true one. According to these results, the global second-order characteristic is able to correctly identify the true intensity. This result in corroborated by the values of the $\chi ^2$, equal to 7.05, 8.32, and 19.26, for the true, the ETAS, and the constant intensities, respectively. However, the main reason for using LISTA functions is that the global summary statistics do not provide information on the individual points, and thus cannot identify which points influence most the goodness-of-fit of the model.

Fig. 7

Difference between the global K-functions and their theoretical values, weighted by: (a) the true intensity; (b) the intensity estimated with the etasFLP (Adelfio and Chiodi 2015) package of R; (c) the constant intensity given by $\hat{\lambda }=n/\vert L \vert \vert T \vert$

Therefore we compute the LISTA functions, and again compare them to their theoretical values, by computing the $\chi _i^2$ values (see Fig. 8). Clearly, LISTA functions help to identify the true intensity as well as the best fit, in a similar manner to their global counterparts. Nonetheless, LISTA functions present further advantages. The main usefulness of carrying out diagnostics with LISTA functions is to get an insight into the local behaviour of the points of the analysed pattern, and therefore to identify the influential points. We indeed identify influential points by retaining only those points that exhibit $\chi _i^2$ values larger than a fixed threshold. For instance, fixing the 75th percentile of the distributions of the $\chi _i^2$ values, and considering as influential points those with a $\chi _i^2$ larger than this percentile, we obtain the black points in Fig. 9.

Fig. 8

Distributions of the $\chi _i^2$ statistics computed for assessing the differences between the local LISTA functions on the network and the theoretical expected value, calculated for each individual point of the simulated ETAS process

Fig. 9

In black: Influential points for the three intensities, which exhibit the $\chi _i^2$ values above the 75th percentile. The triangles indicate the points with $\chi _i^2$ values above the 95th percentile

We note that the influential points vary with the fitted intensity. In detail, we notice that the influential points obtained for the constant intensity are mainly placed in the clusters, and therefore indicate that we might need a better model in those areas, such as an inhomogeneous one. Furthermore, the resemblance between the true intensity and the kernel one (in terms of actual points, and not in the $\chi _i^2$ values only), indicates that the kernel achieves a better fit, if compared to the constant one. In Fig. 10, we report the distributions of the $\chi _i^2$ values of the influential points detected and shown in Fig. 9. It is clear that the influential points of the wrong (constant) intensity are the ones contributing most to the bad fit of the model, if compared to those of the true (or the most suitable) intensity.

Fig. 10

Distributions of the $\chi _i^2$ values above the 75th percentile of the detected influential points detected, and shown in Fig. 9 in black

Figure 11 shows an additional advantage of using LISTA functions. Indeed, their surfaces can be displayed and the different clustered structures of each individual point can be elicited. Of course, depending on the application, this can provide useful information. It can be interesting to plot the LISTA functions of the influential points, providing information on the clustered structure of individual points. Also, the degree of the identified clustered structure can be elicited and contextualised in the application. In particular, Fig. 11 depicts the LISTA functions of the influential points, above the 95th percentile, detected and shown in Fig. 9 as triangles. Under the Poisson case, we expect the K-function to take rh values. In graphical terms, we expect the LISTA functions in Fig. 11 to increase with those two distances, namely, to have constantly increasing values, from bottom-left to top-right. Points with different behaviour, such as the $41^{st}$, indicate that the local structure around that point is somewhat deviating from the Poissonian case. In particular, for that point, we expect neighbouring points to start grouping together abruptly after a specific spatial distance. Overall, the temporal distance does not seem to influence the local structure.

Fig. 11

LISTA functions of the influential points, above the 95th percentile, detected and shown in Fig. 9 as triangles. From top to bottom: influential points of the true intensity, the ETAS one, and the constant one

6 Application to traffic data

We analyse traffic accidents in the city of Medellin (Colombia), a dataset containing the locations of traffic accidents consisting of 1215 traffic accidents in 2019 in downtown Medellin (see Fig. 12).

We evaluate and compare an inhomogeneous intensity fitted with a kernel density with a constant one. Figure 13 depicts the distribution of the $\chi _i^2$ statistics computed for assessing the differences between the local LISTA functions on the network and the theoretical expected value, weighted by an inhomogeneous intensity, and by a constant one. Clearly, the inhomogeneous one seems to report a much better fit.

We therefore separate the most influential points from the rest, retaining only those who exhibit chi-squared values larger that the 99th percentile. We note in Fig. 14 that the two fitted intensities generally identify different influential points. It is also possible to represent the individual LISTA functions, and look at the structures that are more clustered (Fig. 15). We indeed can notice different behaviours. For instance points 49 and 777 display a quite regular and common K-function with the lowest values at the lowest spatial and temporal ranges; however, other points display peculiar behaviours, with the largest values for some specific spatial and temporal ranges. The same considerations can be drawn for the LISTA functions weighted by the constant intensity (Fig. 16).

So, on the one hand, the constant intensity is not appropriate and this means that an inhomogeneous intensity should be fitted to the data. Nevertheless, this information could have been elicited from the inspection of the global statistics. On the other hand, using LISTA functions has allowed to identify the most influential points and this result could lead to a further study on those particular points. This could be done by treating the point patterns as marked, considering some available characteristics of the events. Then, visual inspection of the areas where most of the influential points occur may guide further analysis on these spatio-temporal areas. For instance, one could think of fitting a parametric intensity that accounts for some characteristics of the network, by means of some available spatial covariates.

Fig. 12

The dataset represents the spatio-temporal locations of traffic accidents in the downtown of Medellin (Colombia) in the period of 2019. Left panel: The projection of data onto the network. Right panel: The cumulative number of data points versus occurrence time

Fig. 13

Distributions of the $\chi _i^2$ statistics computed for assessing the differences between the local LISTA functions on the network and the theoretical expected value, for the Medellin dataset. Left panel: The LISTA functions are weighted by an inhomogeneous intensity. Left panel: The LISTA functions are weighted by an inhomogeneous intensity. Right panel: The LISTA functions are weighted by a homogeneous intensity

Fig. 14

In black: Influential points for the Medellin dataset, identified as those that exhibit chi-squared values larger than the 99th percentile. Left panel: The LISTA functions are weighted by an inhomogeneous intensity. Left panel: The LISTA functions are weighted by an inhomogeneous intensity. Right panel: The LISTA functions are weighted by an homogeneous intensity

Fig. 15

LISTA functions of the influential points associated to the inhomogeneous intensity, for the Medellin dataset

Fig. 16

LISTA functions of the influential points associated to the homogeneous intensity, for the Medellin dataset

7 Conclusions

In this work, we have extended Local Indicators of Spatio-Temporal Association (LISTA) functions on linear network to their inhomogeneous versions proposing the inhomogeneous version of the local spatio-temporal K and pair correlation functions on the networks, introduced by D’Angelo et al. (2021). Through simulations, we have shown that the LISTA functions are useful for diagnostics of models specified on the networks. We have proved, and also shown by simulations, that the proposed methods do not rely on any particular model assumption on the data, and thus they can be applied whatever is the generator model of the process. Therefore, an appealing feature of the method is that it can be used for assessing the goodness-of-fit of spatio-temporal point processes occurring on linear networks characterised by any first-order intensity function.

Given the growing availability of proposed models on the networks, we believe that the proposed diagnostic approach could be of great interest. Just to cite some recent papers using this methodology, D’Angelo et al. (2022a) dealt with parametric intensity specification of inhomogeneous first-order intensities on networks to analyse the spatio-temporal distribution of visitors’ stops by touristic attractions in Palermo (Italy). The authors fitted a Gibbs point process model with mixed effects for the purely spatial component, as well as a spatio-temporal log-Gaussian Cox process, adapting them to the underlying road network. Then, Gilardi et al. (2021) proposed a spatio-temporal model to analyse the ambulance interventions that occurred in the road network of Milan (Italy) from 2015 to 2017, adopting a non-separable first-order intensity function, based on a Poisson regression model for the temporal component, and a network kernel function for the spatial dimension. Finally, D’Angelo et al. (2022c) has taken into account the self-exciting behaviour of points, proposing a spatio-temporal Hawkes point process model adapted to linear networks to analyse crime data in Bucaramanga (Colombia).

We note that the LISTA functions can be used to fit the first-order intensity function. Namely, summary statistics such as the K-function or the pair correlation function are commonly employed for some fitting procedures such as the minimum contrast for the estimation of the correlation parameters of Cox processes. Considering the individual contributions of such statistics could provide local estimates: see D’Angelo et al. (2022b) for a recent application to local log-Gaussian Cox Processes. Consequently, the LISTA functions proposed in our paper could be used to fit intensities on linear networks, even though, to the authors knowledge, no attempt has been made in this direction so far.