Density Estimation of Spatio-temporal Point Patterns Using Moran ’ s Statistic

Moran’s Index is a statistic that measures spatial autocorrelation, quantifying the degree of dispersion (or spread) of objects in space. When investigating data in an area, a single Moran statistic may not give a sufficient summary of the autocorrelation spread. However, by partitioning the area and taking the Moran statistic of each subarea, we discover patterns of the local neighbors not otherwise apparent. In this paper, we consider the model of the spread of an infectious disease, incorporate time factor, and simulate a multilevel Poisson process where the dependence among the levels is captured by the rate of increase of the disease spread over time, steered by a common factor in the scale. The main consequence of our results is that our Moran statistic is calculated from an explicit algorithm in a Monte Carlo simulation setting. Results are compared to Geary’s statistic and estimates of parameters under Poisson process are given.


Introduction
Moran's Index is a measure of spatial autocorrelation within a domain area; it has and continues to be applied in many fields.As described by Moran (1950), when given a set of spatial variates such as (x, y) (defined on a two-dimensional discrete area), we may want to investigate whether there is any evidence that spatial autocorrelation is present overall or in neighboring clusters based on selected features.For example, in epidemiology, we may be interested in how a disease spreads based on the number of people infected in a given area.By adding a time component to the model, Moran's indices can lead to better understanding of the spread.Indeed, modeling disease spatio-temporal spread using ordinary least squares or weighted least squares will lead to much bias in the estimates and violation of model assumptions (Pace & LeSage, 2009).A simple global Moran's value cannot explain differences, area spread, or the relationship of some natural ubiquitous phenomena.This has led many authors such as Anselin (1995) to partition the global Moran value into a sum of local indices of spatial association (LISA), because data show local trends that are not shared globally.However, the impact of time was not incorporated.Spatio-temporal autocorrelation was first introduced by Cliff and Ord (1975) and the concept has been explored by many others over the years (see for example, Martin & Oeppen, 1975;Wang & He, 2007;and Lee & Li, 2017).We propose to evaluate spatiotemporal association, extending Moran's index spatial autocorrelation models over iterative time in dynamic subareas under the Poisson process.A full description of the model is presented as a way to measure spatio-temporal non-stationary association.
In this paper, we design a multilevel Poisson process that is time dependent, and the dependence among the levels is captured by a metric of closeness, the rates of increase of the disease spread over time, steered by a common factor in the scale at each time step.Meddens and Hicke (2014) describe spatial and temporal patterns of tree mortality from beetle ecology and dynamics in parts of Colorado and Wyoming.The spatial point process experiment in Vaillant et al. (2011) is revisited the same number of time intervals.In our case, we partition the domain area into subareas and propose incremental time rates of increase under critical scales from each subarea.The scale parameter effectively creates a region of finer detail around disease occurrence.Using Monte Carlo techniques, simulated data are used to investigate model fit based on Moran's statistic under discrete time Markov chain.Diagnostics and comparison statistics are performed.Calibration of clusters for predictive locations is performed to enhance understanding of Moran's autocorrelation, trend in time, and functionality of model analysis.

Temporal Version of Moran Spatial Autocorrelation
For a measurable space S = R d , a domain D ⊂ S , and a collection of random points x 1 , x 2 , . . ., let N(D) represent the count or number of points x i that capture location occurrence of some disease/event in the domain D. Then N(D) can be seen as a measure and is called a point process.For nonoverlapping areas D 1 , D 2 , . . ., N(D 1 ), N(D 2 ), . . .are mutually independent and for each area D, N(D) follows a Poisson distribution with mean λ|D|, where |D| is the area (Lebesgue) measure of D and λ is the intensity of the process.As a process, a counting may be associated to it and is defined as: To every domain, we consider an associated intensity measure λ with λ(D) = E(N(D)), the expected number of points that are obtained in D. Such measure can be assumed to be finite if N(D) < ∞, hence the name finite measure.
Here, we consider the nonhomogeneous Poisson point process with non-constant intensity function λ in time say, λ(t) with t ∈ [0, T ] with the main assumption that points are independent of each other.The nonhomogeneous Poisson process has the following three important properties as described in Baddeley et al. (2016): • conditional property: given exactly n points in a region D, these points are mutually independent and each point has the same probability distribution over D, with probability density f (t) = λ(t)/µ, where µ = ∫ D λ(t)dt.
• superposition property: We consider the conditional distribution of the occurrence of disease in the partitioned subareas defined by Voronoi cells induced by a set of points (sites) as described in Kallenberg (2001).These cells are defined as: where N(R d ) is the class of locally finite measures on R + and D r denotes the ball of radius r around the point s.Doing so, we have a local precise measurement of the density and neighborhoods.We set a partition of D into subareas D 1 , D 2 , . . .where a finite measure or "locally finite" measure is such that µD k < ∞ for all k ≥ 1.
Note that λD = λ(•, D) is a random variable in [0, +∞] for every D ⊆ S = R d .Such an idea is also described in Thäle and Yukich (2016) with properties of Poisson-Voronoi generated sets.
Using Equation (1), Vaillant et al. (2011) modeled the propagation of sugar cane yellow leaf virus with a focus on the spatial spread of disease over six time periods (weeks 6, 10, 14, 19, 23, and 30 in the growing season).For each pair of consecutive observation dates (t i−1 , t i ), i = 1, 2, . . ., 6, they defined a Moran-type index based on a nearest neighbor scheme as follows: where D denotes the discrete set of plant locations, T x denotes the date (time variable) of virus detection for plant x, 1 [0,t i−1 ] (T x ) is an indicator whether time T x falls in the interval [0, t i−1 ], and w x,y denotes the weights, which are non-zero only if x and y are neighbors, located within the same subarea.There are several possible weight functions as we will describe later in this section.
We will utilize a similar space-time definition of Moran's index.Since the partitioned areas vary from one time to the next, we define D t k as subarea k ≥ 1 at time t and then we will define a Moran-type autocorrelation index on each measurable subset D t k as: where w u,u ′ denotes the spatial weight between points u and u ′ of disease/event occurrence, T u and T u ′ denote the time of detection of u and u ′ , and In the definition of Moran's autocorrelation index, w u,u ′ represents a spatial weight between any two distinct event locations generated within disc D t k , which could be a function, e.g., (i) the inverse distance between two points; (ii) the inverse distance squared between two points; (iii) an estimate of the autocorrelation/semivariance statistic; (iv) the "geographical" weights defined as w i j = { e (−d i j /r) j i, 0 otherwise, where r is the maximum distance in the minimum tree that spans all points and does not have any nodes that link back to itself as in Murakami et al. (2017).
As stated in Chen (2012), selection of the weight function objectively is an open question.Different wight functions have different spatial effects.Later in Section 4, we give reasons why we choose a particular weight function.
Also, while values of the Moran-type statistic may not always be between -1 and 1 (the full range depends on the weights), they are essentially in the same spirit as larger values indicate larger autocorrelation measures.
We will use a modified version of Geary's C autocorrelation index (Geary, 1954;Sokal, Oden, & Thomson, 1998) as a comparison to determine and evaluate the trend between the two measures.We will define it as C t k : with w u,u ′ as in Equation ( 2) and w u,u ′ , the sum of weights for all points generated across all discs at time t.

Spatio-temporal Process
For a time period of length T partitioned into m time subintervals The number of occurrences is defined as X t r and the location of event occurrences is defined as Xt r .The process follows a discrete time Markov chain as defined in Resnick (2002): for some function of propagation from location i to location j and subarea containing points i and j in time interval (t r−1 , t r ).The associated count within D t k defines a sequence {X t r } of number of occurrences within [t r−1 , t r ) in subarea D t k and follows the Markov property; i.e. given the current state of the system, we can make predictions about the future state without regard for previous states.This is a discrete time finite Markov chain.
The transition probability from location i to location j in interval time (t, t+s), denoted by P i j (s), n i and n j .The conditional probability is: To minimize complexity, we define the count process {X t , t ≥ 0} of occurrences between consecutive times t r−1 and t r to be where λ(t) denotes the local intensity for any t within [t r−1 , t r ).Such a process has stationary transition probabilities and transition matrix of the counts Q = (Q i j ) constructed from the transition probabilities, P i j (s), where Q i j is the (i, j) count from P i j (s) from consecutive time periods t 1 , t 2 , . . ., t m adjusted with subarea containing locations i and j, and We extend D t k such that they form a nested sequence D t r k(t r ) ⊂ D t r−1 k(t r−1 ) , the subarea within time (t r−1 , t r ) and count X t k(t) , 1 ≤ k(t) ≤ n(t).Following Resnick (2002), we can define the nested subareas D t k(t) and each fixed time t, a sequence {E t k(t) } k≥0 of independent and identically distributed exponential random variables with unit mean such that , where |D t k(t) | is the area (Lebesgue) measure within D t k(t) .Then for u > 0, 1 ≤ k(t) ≤ n(t), and at location i ∈ D t k(t) .
By discretizing time, we then define two finite sequences {(X t r ), (t r )}, where X t k(t) with locations Xt k(t) for t n < t < t n+1 and {t r }, r = 1, . . ., m is the sequence of times when the process is observed.The sequence {X t } is the cumulative count of occurrence within interval time [t r−1 , t r ) in a subset of space D such that t r − t r−1 is conditionally independent and exponential given X t r−1 .More precisely, the sequence {X t k(t) } is defined as follows: and . . . and . . . and where n(t r ) is the total number of points generated at time , and X(t) is the total count for t r ≤ t < t r+1 , 1 ≤ r ≤ m.
Setting t ∞ = lim r→∞ t r the process {X t } becomes explosive.However, The pair sequence {X t , t r } has the following two properties.
(i) The sequence {t r − t r−1 , 1 ≤ r ≤ m} are conditionally independent and exponentially distributed given {X n }; i.e. memoryless in time: for u m > 0, m = 1, . . ., n, (ii) The distributional structure of {X t , (t r )} is a function of the transition matrix and the limiting function: For u > 0, i ∈ S , The nested sequence is shown in Figure 1, where E t k(t) equals the count at time t for subarea D t k(t) , with 1 ≤ k(t) ≤ n(t), and n(t) represents the number of points generated at time t and Q is the (i, j) transition probability of count from locations i to j between times t r−1 and t r .
. We propose a comparable implementation and compute the successive Moran's spatio-temporal autocorrelations under such setting.

Simulation Study
Here we provide an illustration of the theory built in the previous section.We conduct the simulation using R with the spatstat package and represent the weights of the Moran's index as the inverse distance between two points.The weight function was chosen as described in Chen (2012) since we are modeling spread of disease which is a large scale complex system and we intend to capture transitions globally to locally.We begin with an observed unit area (1 × 1 unit) and generate a Poisson point processes within that area.The Poisson process is used as general setup since it is widely applicable and is applicable to many real world situations such as disease cases, radioactive decay, and plant propagation.Other processes, for example negative binomial and Cox, may be utilized to explain other space-time point phenomena.
Next, we define the Voronoi cells around each point (e.g. the locus of locations nearer to that point than any other).Each location site is then used as the center point to generate a disc with diameter equal to the minimum of the distance to the Voronoi edge and the distance to the observed area edge (edge of D).This effectively means that an infected site can only subsequently infect sites that are closer to it than any other infected site.These Voronoi discs define the subareas D t k(t) for the point process generated at the next time interval (step).In application, subareas can be tied to locations such as country, state, county, and zip code.The choice of subareas could be used to determine whether more affluent regions are able to deal more effectively with disease outbreak, for example.The Moran statistic is calculated at each step as the sum of the inverse distances between an offspring site and another offspring site.This format will continue through time t = 4 adhering to the Markov property (t = 4 chosen to follow results from a forestry example in Meddens and Hicke, 2014;or as in Vaillant et al., 2011).Each subarea will be magnified by a scale parameter, α > 0, such that the new radius is α times the previous radius with intensity function λ * t, following the spatial resolution of aggregated grid cells that experienced a 1% tree mortality as in Meddens and Hicke (2014).

Algorithm: Iterative local Moran statistics
Procedure: (i) Define the partition and generate a Poisson process in the area with intensity λ.
(ii) Calculate the Moran statistic.
(iii) Define new subareas with point from previous process as center and radius equal to the minimum of distance to the edge of the window or distance to the edge of the Voronoi cell.
(iv) Rescale subarea and iterate local points from a Poisson process with intensity λ * t.
(v) Compute the local Moran statistic for each subarea at that time step.

Repeat steps (iii) through (v).
Stop when all times are reached.

End
Figure 2 shows the original observed area with the seven points generated at t = 1 alongside the discs generated from these points and the points generated at t = 2.

Results
The Moran value for t = 1 equals 42.57.The number of points generated in each Voronoi subarea, the rescaled area, and Moran values generated at time t = 2 are displayed in Table 1 with λ = 4 and α = 5.For example, at point 1 in Figure 2(a), we generate a disc of maximum radius within the Voronoi cell centered at the point generated at the previous time with 20 points as shown in Figure 2(b).The number of points generated in other Voronoi cells is displayed in Table 1.The Geary's values are also computed.The ordering of Moran and Geary's values display resemblance as they "rank" the strength of correlations in a similar way.The process is described at time t = 3 with disc 1 having 20 subareas (see Figure 3), disc 2 having 14 subareas, etc.The non-zero Moran values for disc 1 at time t = 3 are listed in Table 2.After four time intervals, the clustering of points generated is displayed in Figure 4, with red shades representing the most dense areas and white shades representing areas with no instances of disease.Such a representation would not have been clearly observed if the area was not partitioned and if the rate of spread was not tabulated in space and time.This framework will allow us to detect clusters and estimate the density.

Model Analysis
In this section estimation of the parameter, λ, for the Poisson distribution given the data generated in the simulation is proposed using the maximum likelihood method as presented in Baddeley et al. (2016).A scale parameter is incorporated into the Poisson distribution through the use of the area of the Voronoi subarea The parameter estimate becomes λ Using the data generated in the simulation above, we tabulate the estimates of λ at t = 2, 3, and 4 with different scale parameters, α.Table 3 shows these parameter estimates of λ with standard error in parenthesis.For α = 5 at time t = 2, λ = 8.14 with associated SE=1.07 and the true parameter is λ = 8.The data appear to follow a Poisson distribution with intensity λ * t; however, the estimated parameter values for α = 2.5 seem less reliable due to the sparse number of points generated at each time step.We also investigated increasing α as function of time, α(t) = α t−1 .Using an initial value of α = 2.5, the parameter estimates at t = 2, 3, 4 were 10.26 (1.21), 13.44 (0.76), and 15.98 (0.31), respectively.By making the scale parameter a function of time, the parameter estimates appear more precise for small starting values, but under this construct, using α(t) = 10 t−1 , generates up to 30,000 points per subarea at t = 4.As a result, we will use α = 5 for the scale parameter in the following analysis.
Next we investigate any trend in our Moran index versus number of points generated and area of the disc.Figure 5 shows Moran's values plotted versus area and versus number of points.The plots of Moran values for time t = 2 and for disc 1 at t = 3 show an approximate linear trend between Moran and number of points generated, even though there is no apparent relationship between area and Moran.The 3-D plot in Figure 6 shows that at each successive time point the number of points generated versus the Moran gets progressively larger.This is due to the number of points generated increasing as time increases, but the area increase is moderated as time increases.This results in points that are closer together, thus increasing the sum of inverse distances and the Moran value.Focusing on the seven points generated at t = 1 and the associated discs at t = 2, we investigate the Markov property of the simulation results.Table 4 shows the number of points generated within each of the original seven discs at times t = 2, 3, 4 and also includes the area of each disc.The transition matrix for each disc Q t k is developed using the ratio of points generated at each time to the total number of points generated over all times.This matrix is upper triangular since we cannot transition backward in time or stay in the same time.We then conduct cluster analysis to determine similarities between discs.Using the Canberra distance measure, , to create a distance matrix using the counts in Table 4, the dendrogram tree in Figure 7 was generated to look at similarities in clusters.It can be inferred that disc association is related to area and abundance of occurrence of events.The separation of clusters two and three to govern the tree shows that even though they are spatial neighbors, they are perceived to reflect distinct evolutionary paths.Such finding could add to the conclusion given in Meddens and Hick (2014) when they state that spectral responses of lodgepole pine, ponderosa pine, and mountain pine beetles lead to different stages of tree maturity.

Summary
Our designed model has shown that global properties of spatio-temporal disease spread can be tuned in by modifying the space and adding time efficiently.The results show that long range, spread feature and variability yielded spatial clustering of occurrence events.To obtain a quantitative analysis of the performance of the Moran's values, we compared them with adjusted Geary's C statistics.We controlled for targeted subareas and such an approach provides a framework for understanding "disorganized" or "disordered" evolution of some natural organisms in the form of cluster analysis.The distribution of subareas can be found based on our density estimation or concentration under the modified Moran's Index in space and time.
Extension to different dynamics other than the Poisson will be considered in future studies and widen the scope of the spatio-temporal model incorporating covariates.

Figure 6 .
Figure 6.Plot of Moran values t r ) for disc k.

Table 2 .
Non-zero Moran and Geary's c values for disc 1 at t = 3