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

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).


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 . In particular, following  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 .
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 spatiotemporal 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.

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 R 2 × R, where a point (u, t) ∈ X corresponds to an event at u ∈ R 2 occurring at time t ∈ R. A typical realisation of a spatio-temporal point process X on R 2 × R is a finite set {(u i , t i )} n i=1 of distinct points within a bounded spatio-temporal region W × T ⊂ R 2 × R, with area |W | > 0 and length |T | > 0, where n ≥ 0 is not fixed in advance. In this context, N (A × B) denotes the number of points of a set (A × B) ∩ X , where A ⊆ W and B ⊆ T . As usual (Daley and Vere-Jones 2007), when N (W × T ) < ∞ 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 (u, t), the events that are close to (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 (u, t), which can be expressed by the Cartesian product as b ((u, t) where || · || denotes the Euclidean distance in R 2 . Note that b ((u, t), r , h) is a cylinder with centre (u, t), radius r , and height 2h. Product densities λ (k) , k ∈ N and k ≥ 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 ( 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 λ and the (second-order) product density λ (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 realvalued process, respectively. Then, the first-order intensity function is defined as where du×dt defines a small region around the point (u, t) and |du×dt| is its volume. The second-order intensity function is given by Finally, the pair correlation function can be interpreted formally as the standardised probability density that an event occurs in each of two small volumes, du × dt and dv × ds, in the sense that for a Poisson process, g((u, t), (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 secondorder 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 = ||u − v|| and h = |t − s|.
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 The simplest expression of an estimator of the spatio-temporal K -function is given aŝ For a homogeneous Poisson process E[K (r , h)] = πr 2 h, regardless of the intensity λ.
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,K (r , h) is compared with the theoretical E[K (r , h)] = πr 2 h. ValuesK (r , h) > πr 2 h suggest clustering, whileK (r , h) < π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) aŝ where λ(·, ·) is the first-order intensity at an arbitrary point. We know that E[K I (r , h)] = πr 2 h, that is the same as the expectation ofK (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 estimatorK 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 λ(u, t), the expectation ofK I (r , h) should be close to E[K (r , h)] = πr 2 h for the Poisson process. For instance, valuesK I (r , h) greater than πr 2 h indicates that the fitted model is not appropriate, since the distances computed among points exceed the Poisson theoretical ones.

Local Indicators of Spatio-Temporal Association functions
Definition 3  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 λ (2)i (·, ·) denotes the local version of the spatio-temporal product density for the event (u i , t i ), then, for fixed r and h, it holds that whereλ , with r > > 0 and h > δ > 0, and κ a kernel function with spatial and temporal bandwidths and δ, 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 spatiotemporal local indicator functions is proportional to the global statistic, can be called a LISTA statistic.
In , 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λ = n/(|W ||T |), they propose the local version of (2) for the i-th event (u i , t i ) aŝ and the local version of (3) aŝ with (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 spatiotemporal 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 π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.
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 = ∪ n i=1 l i ⊂ R 2 is commonly taken as a finite union of line segments l i ⊂ R 2 of positive length. A line segment is defined as are the endpoints of l i . For any i = 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 × T , where L is a linear network and T is a subset (interval) of R.
We hereafter focus on a spatio-temporal point process X on a linear network L with no overlapping points (u, t), where u ∈ L is the location of an event and t ∈ T (T ⊆ R + ) is the corresponding time occurrence of 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 where | · | is a numerical distance, and d L (·, ·) 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 ⊆ L × T , N (X ∩ A) ∈ 0, 1, . . ., is the number of points of X restricted to A, whose expected value is denoted by where ν, the intensity measure of X , is a locally finite product measure on L × 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 λ (k) exist, for any non-negative measurable function f (·) on the product space L k , we have where = indicates that the sum is over distinct values. Assume that X has an intensity function λ(·, ·), hence Eq. (7) reduces to The second-order Campbell's theorem is obtained from (7) Assuming that X has a second-order product density function λ (2) (·, ·), we then obtain Finally, an important result concerns the conversion of the integration over L × T to that over R × R (Rakshit et al. 2017). For any measurable function f : is the number of points lying exactly at the shortest-path distance r ≥ 0 and the time distance h ≥ 0 away from (u, t).

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 secondorder 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 × T with constant intensity function λ ≥ 0. Let the homogeneous K -function be given by Then, X is called second-order pseudostationary and isotropic if K L ((u, t), r , h) does not depend on the point (u, t), and we then write K L ((u, t), r , h) = K L (r , h).

Definition 6 Moradi and Mateu (2020) For a homogeneous Poisson point process on
The authors proposed a nonparametric estimator of both the homogeneous K -functionK and pair correlation function where |L| > 0 and |T | > 0 are the total length of the network L and of the time interval T , respectively, and κ 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 secondorder characteristics, namely the inhomogeneous K -function and the inhomogeneous pair correlation function whereλ(·, ·) is an estimate of the intensity function, and M((u, t), r , h) is the number of points lying exactly at the shortest-path distance r ≥ 0 and the time distance h ≥ 0 away from the point (u, t).
Moradi and Mateu (2020) recommend using some normalisation factor of the form in such a way that the updated estimates of the inhomogeneous K -and pair correlation functions 1 D(x)K (r , h) and 1 D(x)ĝ (r , h) provide estimators with low bias and variance. They also prove that for Poisson processes, and for any r , h > 0, it holds thatK L,I (r , h) = rh andĝ L,I (r , h) = 1, which is useful in model selection and hypothesis testing.
Therefore, a relevant result is that the expectation ofK L,I (r , h), when the intensity used for the weighting is the true generator model, is the same as the expectation ofK L (r , h) for the Poisson process, that is E[K L (r , h)] = rh. Therefore, as in the Euclidean case, the weighted estimatorK 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 λ(u, t), the estimated valueK L, 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).

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 (u i , t i ) on a linear network aŝ and the corresponding local pair correlation function weighted by the reciprocal of the normalising factor obtaining 1 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λ = n/(|L||T |), obtaininĝ whereλ 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
(By theorem (8) and since λ (2) ≡ λ 2 for a Poisson process) (By the change of variables in (9)) (taking the conditional expectation with respect to its filtration) Theorem 2 If the process is weighted by the true intensity function, the expectation of K i L,I (r , h) is the same as the expectation ofK i L (r , h).
(taking the conditional expectation with respect to its filtration) We know that the expectation ofK L,I (r , h), when the intensity used for the weighting is the true generator model, is the same as the expectation ofK L (r , h) for the Poisson process, that is E[K L (r , h)] = rh.
We have proved that the same result holds for the local versionK i L,I (r , h), meaning that the our proposed local estimatorK i L,I (r , h) for a general point process, weighted by the true intensity function, has the same expectation of the local estimatorK i L (r , h) under the Poisson case.
Based on this result, we can use our proposed local estimatorK i L,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 , 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 λ(·, ·).
Basically, departures of the LISTA functionsK i L,I (r , h) from the Poisson expected value rh directly suggest the unsuitability of the intensity function λ(·) 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.

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  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.

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.

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

Inhomogeneous spatio-temporal point process
To generate inhomogeneous spatio-temporal point processes on a linear network driven by the spatio-temporal intensityλ(u, t), we proceed with the following steps: 1. Define the generating intensity function λ 0 (u, t) =λ(u, t). 2. Set an upper bound λ max for λ 0 (u, t). 3. Simulate a homogeneous Poisson process with intensity λ max and denote by N the number of the generated points, with coordinates (u , t ). 4. Compute p(u , t ) = λ (u ,t ) max(λ(u ,t )) for each point (u , t ) of a homogeneous Poisson process. 5. Generate a sample x of size N from the uniform distribution on (0, 1). 6. Thin the simulated homogeneous Poisson process x retaining the n ≤ N locations for which x ≤ p.
The algorithm provides a point pattern 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 (u, t) of the generating intensity functionλ(u, t) occurs only on the network line segments. The simulation is performed through the function rpoistlpp() in the R package stlnpp (Moradi and Mateu 2020). In Fig. 1, a simulated inhomogeneous spatiotemporal Poisson process is depicted, with intensity λ(x, y, t) = exp(8.25 − 4y − 2t), decreasing together with the y-axis and with time.

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 selfexciting 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 where H t is the past history of the process, μ is the background general intensity, and f (u) is the spatial density. Concerning the triggered component, η j = β Z j is a linear predictor, with Z j the external known covariate vector, including the magnitude, and θ = (μ, κ 0 , c, p, d, q, β) are the parameters to be estimated. In the usual ETAS model, with the only external covariate representing the magnitude, β = α, as The algorithm provides a point pattern 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 θ = (μ, κ 0 , c, p, α, 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 β = α.

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 where S is a Gaussian process with E(S(u, t)) = μ = −0.5σ 2 and so E(exp S(u, t)) = 1 and with variance and covariance matrix C(S(u, t), S(v, s)) = σ 2 γ (r , h), with γ (·) 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(u, t), with covariance function C((u, t), (v, s)) and mean function μ(u, t). 1.b Define the generating intensity function λ 0 (u, t) =λ(u, t) exp(S(u, t)).
3. Simulate a homogeneous Poisson process with intensity λ max and denote by N the number of the generated points, with coordinates (u , t ). 4. Compute p(u , t ) = λ (u ,t ) max(λ(u ,t )) for each point (u , t ) of an homogeneous Poisson process. 5. Generate a sample x of size N from the uniform distribution on (0, 1). 6. Thin the simulated homogeneous Poisson process retaining the n ≤ N locations for which x ≤ p.
This algorithm provides a point pattern 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, . . . , 4. We consider a separable spatio-temporal exponential covariance with parameters σ = 2, φ = 0.2 and θ = 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).

Results
For each of the 3 previously defined scenarios, we compute the chi-squared statistics between the LISTA functions and their theoretical valuesK i L,T heo (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 χ 2 i values are obtained following the expression one for each point in the point pattern. We weight the LISTA functionsK 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.

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 -  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 χ 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. Therefore we compute the LISTA functions, and again compare them to their theoretical values, by computing the χ 2 i 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 χ 2 i values larger than a fixed threshold. For instance, fixing the 75th percentile of the distributions of the χ 2 i values, and considering as influential points those with a χ 2 i larger than this percentile, we obtain the black points in Fig. 9. 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 χ 2 i 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 χ 2 i 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. 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.

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).  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 We evaluate and compare an inhomogeneous intensity fitted with a kernel density with a constant one. Figure 13 depicts the distribution of the χ 2 i 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.

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.
Funding Open access funding provided by Universitá degli Studi di Palermo within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.