Section 9.1 Exploring spatial dependence in point pattern data
Subsection 9.1.1 The empirical K-function
In Chapter 7 we assessed homogeneity, or lack of, in point pattern data. Another important aspect of point location is whether they show some kind of interpoint dependence. In observing a point pattern, one is interested in what are referred to as first order properties (varying intensity across the study surface vs. homogeneity) and second order properties (independence vs. correlation among the points). As [170] explain:
Spatial variation occurs when the risk is not homogeneous in the study region (i.e. all individuals do not have the same risk) but cases appear independently of each other according to this risk surface, whilst clustering occurs when the occurrence of cases is not at random and the presence of a case increases the probability of other cases appearing nearby. (p. 204)
In Chapter 7 we covered the Chi Square test and we very briefly discussed how it can only help us to reject the notion of spatial randomness in our point pattern, but it cannot help us to determine whether this is due to varying intensity or to interpoint dependence.
Both processes may lead to apparent observed clustering:
-
First order clustering near particular locations in the region of study takes place when we observe more points in certain neighbourhoods or street segments or whatever our unit of interest is;
-
Second order clustering is present when we observe more points around a given location in a way that is associated with the appearance of points around this given location.
A technique commonly used to assess whether points are independent is Ripley’s K-function. For exploring global clustering in a point pattern, we use the K-function. [156] define the K-function as the "cumulative average number of data points lying within a distance \(r\) of a typical data point, corrected for edge effects, and standardised by dividing by the intensity" (p. 204). This is a test that, as we will see, assumes that the process is spatially homogeneous. As [172] indicate:
"The basic idea of K-function is to examine whether the observed point distribution is different from random distributions. In a random distribution, each point event has the equal probability to occur at any location in space; the presence of one point event at a location does not impact the occurrence possibilities of other point events." (p. 613)
Imagine we want to know if the locations of crime events are independent. Ripley’s K-function evaluates distances across crime events to assess this. In this case, it would count the number of neighbouring crime events represented by the points found within a given distance of each individual crime location. The number of observed neighbouring crime events is then traditionally compared to the number of crime events one would expect to find based on a completely spatially random point pattern. If the number of crimes found within a given distance of each individual crime is greater than that for a random distribution, the distribution is then considered to be clustered. If the number is smaller, the distribution is considered not to be clustered.
Typically, what we do to study correlation in a point pattern data is to plot the empirical K-function estimated from the data and compare this with the theoretical K-function of the homogeneous Poisson process. So, for this to work, we need to assume that the process is spatially homogeneous. If we suspect this is not the case, then we cannot use the standard Ripley’s K-methods. Let’s illustrate how all this works with the example of some homicide data in Buenos Aires.
Subsection 9.1.2 Homicide in Buenos Aires
For this exercise we travel to Buenos Aires, the capital and largest city in Argentina. The government of Buenos Aires provides open geocoded crime data from 2016 onwards. The data is freely available from the official source. But we have added it already as supplementary material available with this book. Specifically, we are after the file
delitos.csv.
We will filter on "homicidios dolosos" (intentional homicides) and we are using all homicides from 2016 and 2017. After reading the .csv file, we turn it into a
sf object and plot it. We are also reading a GeoJSON file with the administrative boundary of Buenos Aires. The code below uses functions in ways we have demonstrated in earlier chapters.
delitos <- read_csv("data/delitos.csv") %>% # read csv
filter(tipo_delito == "Homicidio Doloso") %>% # filter international homicides
mutate(year = year(fecha)) %>% # create new variable for 'year'
filter(year == "2016" | year == "2017") # filter only 2016 and 2017
# create sf object of crimes as points
delitos_sf <- st_as_sf(delitos,
coords = c("longitud", "latitud"),
crs = 4326)
# transform to projected coords
delitos_sf <- st_transform(delitos_sf, 5343)
# read in boundary geojson
buenos_aires <- st_read("data/buenos_aires.geojson", quiet = TRUE)
After reading the files and turning them into
sf objects, we can plot them. Let’s use the tmap package introduced in Chapter 3.
tm_shape(buenos_aires) +
tm_polygons() +
tm_shape(delitos_sf) +
tm_dots(size = 0.08, col = "red") +
tm_layout(main.title = "Homicides in Buenos Aires, 2016",
main.title.size = 1)
We will be using the
spatstat functionality; thus, we need to turn this spatial point pattern data into a ppp object, using the boundaries of Buenos Aires as an owin object. We are using code already explained in Chapter 7. As we will see, there are duplicate points. Since the available documentation does not mention any masking of the locations for privacy purposes, we cannot assume that’s the reason for it. For simplicity we will do as in Chapter 7 and introduce a small degree of jittering with a rather narrow radius.
# create window from geometry of buenos_aires object
window <- as.owin(st_geometry(buenos_aires))
# name the unit in window
unitname(window) <- "Metre"
# create matrix of coordinates of homicides
delitos_coords <- matrix(unlist(delitos_sf$geometry), ncol = 2, byrow = T)
# convert to ppp object
delitos_ppp <- ppp(x = delitos_coords[,1], y = delitos_coords[,2],
window = window, check = T)
# set seed for reproducibility and jitter homicide points because of overlap
set.seed(200)
delitos_jit <- rjitter(delitos_ppp, retry=TRUE, nsim=1,
radius = 2, drop=TRUE)
Subsection 9.1.3 Estimating and plotting the K-function
For obtaining the observed K, we use the
Kest function from the spatstat package. Kest takes two key inputs, the ppp object and a character vector indicating the edge corrections we are using when estimating K. Edge corrections are necessary because we do not have data on homicide outside our window of observation and, thus, cannot consider dependence with homicide events around the city of Buenos Aires. If we don’t apply edge correction, we are ignoring the possibility of dependence between those events outside the window and those in the window close to the edge. There is also a convenient argument (nlarge) if your data-set is very large (not the case here) that uses a faster algorithm.
As noted by [156], "most users... will not need to know the theory behind edge correction, the details of the technique, or their relative merits" and "so long as some kind of edge correction is performed (which happens automatically in
spatstat), the particular choice of edge correction technique is usually not critical" (p. 212). As a rule of thumb, they suggest using translation or isotropic in smaller data-sets, the border method in medium-sized data-sets (1000 to 10000 points), and that there is no need for correction with larger data-sets. Since this is a smallish data-set, we show how we use two different corrections: isotropic and translation (for more details, see [156]).
k1 <- Kest(delitos_jit, correction = c("isotropic", "translate"))
plot(k1)
The graphic shows the K-function using each of the corrections and the theoretical expected curve under a Poisson process of independence. When there is clustering, we expect the empirical K-function estimated from the observed data to lie above the empirical K-function for a completely random pattern. This indicates that a point has more "neighbouring" points than it would be expected under a random process.
It is possible to use Monte Carlo simulations to draw the confidence interval around the theoretical expectations. This allow us to draw envelope plots that can be interpreted as a statistical significance test. If the homicide data were generated at random, we would expect the observed K-function to overlap with the envelope of the theoretical expectation. [156], in particular, suggest the use of global envelopes as a way to avoid data snooping. For this we use the
envelope function from spatstat. As arguments we use the data input, the function we are using for generating K (Kest), the number of simulations, and the parameters needed to obtain the global envelope.
set.seed(200)
E_1 <- envelope(delitos_jit, Kest, nsim=39, rank = 1, global = TRUE,
verbose = FALSE) # We set verbose to false for sparser reporting
plot(E_1)
This envelope suggests significant clustering above the 200 metre mark. The problem, though, is that (as noted above) this test only works if we have a homogeneous process. The envelope will also fail to overlap with the observed K-function if there are variations in intensity. We also saw previously the reverse, that it is hard to assess homogeneity without assuming independence. In Chapter 7 we saw how if we suspect the process is not homogeneous, we can estimate it with non-parametric techniques such as kernel density estimation. We can apply this here (see also Chapter 7: Spatial Point Patterns of Crime Events in this book).
set.seed(200)
plot(density.ppp(delitos_jit, sigma = bw.scott(delitos_jit),edge=T),
main=paste("Density estimate of homicides in Buenos Aires"))
[156] discuss how we can adjust our test for independence when we suspect inhomogeneity, as long as we are willing to make other assumptions (e.g., the correlation between any two given points only depend on their relative location). This adjusted K-function is implemented in the
Kinhom() function in spatstat. This function estimates the inhomogeneous K-function of a non-stationary point pattern. It requires as necessary inputs the ppp object with the point pattern we are analysing and a second argument lambda which should provide the values of the estimated intensity function. The function provides flexibility in how lambda can be entered. It could be a vector with the values, an object of class im like the ones we generate with density.ppp(), or it could be omitted. If lambda is not provided, we can provide the parameters as we did in the density.ppp() function, and the relevant values will be generated for us that way. We could, for example, input this parameter as below:
set.seed(200)
inhk_hom <- Kinhom(delitos_jit,
correction = c("isotropic", "translate"),
sigma = bw.scott)
plot(inhk_hom)
Once again we can build an envelope to use in the manner of confidence intervals to interpret whether there is any significant clustering.
set.seed(200)
E_2 <- envelope(delitos_jit, Kinhom, nsim=39, rank = 1, global = TRUE,
verbose = FALSE) # We use this for sparser reporting
plot(E_2)
In this example, once we allow for varying intensity, the spatial dependence no longer holds for a range of distances.
Subsection 9.1.4 Network constrained K-function
We already warned in Chapter 8 about the dangers of using planar point pattern analysis to data that lies along a network. Not recognising this network structure can lead to errors. [172], for example, have shown how the standard K-functions are problematic when studying crime in an urban environment. They discuss how the planar K-function to analyze the spatial autocorrelation patterns of crimes (that are typically distributed along street networks) can result in false alarm problems, which depending on the nature of the street network and the particular point pattern, could be positive or negative. As [156] also indicate when the point locations are constrained to only exist in a network we are likely to find "spurious evidence of short-range clustering (due to concentrations of points in a road) and long-range regularity (due to spatial separation of different roads)" (p. 736).
[172] recommend to use a network K-function that examines distance along the network, by street distance. They argued that this should be "more relevant for crime control policing and police man-power dispatching purpose, considering that police patrol and many other crime management activities are commonly conducted following street networks" (p. 627). This recommendation is being followed in recent studies (see, for example, [177]), although it is still much more common in the professional practice of crime analysis to continue using, incorrectly, the planar K-function.
The network K-function was developed by [174] and was first implemented in SANET. Chapter 7 of [175] also provides a thorough discussion of this adaptation and its computation. We can use functionality implemented in the
spatstat package to use the network constrained K. For that we need a lpp object, which we discussed in Chapter 8. For convenience we will use a lpp data-set about crime around University of Chicago that ships with the spatstat package.
set.seed(200) # set seed for reproducibility
data(chicago) # load chicago data (in spatstat package)
# Estimate Ripley's K with envelope
E_3 <- envelope(as.ppp(chicago), Kest, nsim=99, rank = 1, global = TRUE, verbose = FALSE)
plot(E_3) # plot result
How would we do this now taking the network into consideration? We can use the methods introduced by [174] with the function
linearK(). The function linearK() computes an estimate of the linear K-function for a point pattern on a linear network. For this we need to specify the parameter for correction= to be "none". If correction="none", the calculations do not include any correction for the geometry of the linear network. The result is the network K-function as defined by [174]. We can wrap this within the envelope() function for inference purposes. This will generate an envelope assuming a homogeneous Poisson process. We also limit the number of simulations (nsim =) to 99, and we suppress a verbose output.
set.seed(200) # set seed for reproducibility
E_4 <- envelope(chicago, linearK, correction = "none", nsim= 99, verbose = FALSE)
plot(E_4) # plot result
As noted by [156], there is a difference between what we get with
Kest and with linearK. Whereas the two planar "K-function is a theoretical expectation for the point process that is assumed to have generated our data... the empirical K-function of point pattern data" (that we get with linearK) "is an empirical estimator of the K-function of the point process" (p. 739, the emphasis is ours). A limitation of the Okabe-Yamada network K-function is that, unlike the general Ripley’s K which is normalised, you cannot use it for "comparison between different point processes, with different intensities, observed in different windows... because it depends on the network geometry" (p. 739). This is a problem that was approached by [176] who proposed a "geometrically corrected" K-function. You can see the mathematical details for this function in either [156] or in greater depth in [176]. This geometrically corrected K-function is implemented as well in the linearK() function. If we want to obtain this solution, we need to specify as an argument correction="Ang" or simply not include this argument, for this is the default in linearK():
set.seed(200)
E_5 <- envelope(chicago, linearK, correction = "Ang", nsim= 99, verbose = FALSE)
plot(E_5)
Finally, there is also an equivalent to the geometrically corrected K-function that allows for an inhomogeneous point process that can be obtained using
linearKinhom() (see [173] for details).
These methods all help discuss spatial autocorrelation in point pattern data. In the next section, we will consider how we might find out whether such clustering exists in data measured at polygon (rather than point) level.
