Crime Mapping and Spatial Data Analysis using R

The R package spatstat was created for spatial statistics with a strong focus on analysing spatial point patterns. It is the result of 15 years of development by leading researchers in spatial statistics. We highly recommend the book "Spatial Point Patterns: Methodology and Applications with R" to give a comprehensive overview of this package [156].

When working within the spatstat environment, we will be working with planar point pattern objects. This means that the first thing we need to do is to transform our sf object into a ppp (planar point pattern) object. This is how spatstat likes to store its point patterns. The components of a ppp object include the coordinates of the points and the observation window (defining the study area).

It is important to realise that a point pattern is defined as a series of events in a given area, or window, of observation. Therefore, we need to precisely define this window. In spatstat the function owin() is used to set the observation window. However, the standard function takes the coordinates of a rectangle or of a polygon from a matrix, and therefore it may be a bit tricky to use with a sf object. Luckily, the package maptools provides a way to transform a SpatialPolygons into an object of class owin, using the function as.owin(). Here are the steps.

First, we transform the CRS of our Manhattan polygon into projected coordinates as opposed to the original geographic coordinates (WGS84), since only projected coordinates may be converted to spatstat class objects. As we saw in Chapter 2, measuring distance in metres is a bit easier, so since with points a good deal of the analysis we subject them to involves measuring their distance from each other, moving to a projected system makes sense. For NYC, the recommended projection is the New York State Plane Long Island Zone (EPSG 2263), which provides a high degree of accuracy and balances size and shape well. This projection, however, uses US feet (which was deprecated in 2022), so instead we are going to use EPSG 32118, as an alternative that uses metres as unit of measurement. So first, transform our object using st_transform():

manhattan <- st_transform(manhattan, 32118)

Then we use the as.owin() function to define the window based on our manhattan object.

window <- as.owin(st_geometry(manhattan))

We can plot this window to see if it looks similar in shape to our polygon of Manhattan.

plot(window) # owin similar to polygon

We can see this worked and created an owin object.

class(window)
window

## [1] "owin"
## window: polygonal boundary
## enclosing rectangle: [295965, 307869] x [57329, 79111] 
## units

You can see it prints the units of measurement. But it does not identify them:

unitname(window)

## unit / units

As noted above, EPSG:32118 uses metres. So we can pass this information to our owin object:

unitname(window) <- "Meter"

Now that we have created the window as an owin object, let’s get the points. In the first place, we need to reproject to the same projected coordinate systems we have for our boundary data using the familiar st_transform() function from sf. Then we will extract the coordinates from our sf point data into a matrix using unlist() from base R, so that the functions from statspat can read the coordinates. If you remember sf objects have a column called geometry that stores these coordinates as a list. By unlisting we are flattening this list, taking the elements out of the list.

# first transform to projected CRS too to match our window
nyc_burg <- st_transform(nyc_burg, 32118)

# create matrix of coordinates
sf_nyc_burg_coords <- matrix(unlist(nyc_burg$geometry), ncol = 2, byrow = T)

Then we use the ppp() function from the spatstat package to create the object using the information from our matrix and the window that we created. The check argument is a logical that specifies whether we want to check that all points fall within the window of observation. It is recommended to set it to TRUE, unless you are absolutely sure that this is unnecessary.

bur_ppp <- ppp(x = sf_nyc_burg_coords[,1], y = sf_nyc_burg_coords[,2],
                   window = window, check = T)

## Warning: data contain duplicated points

Notice the warning message about duplicated points? In spatial point pattern analysis an issue of significance is the presence of duplicates. The statistical methodology used for spatial point pattern processes is based largely on the assumption that processes are simple, that is, that the points cannot be coincident. That assumption may be unreasonable in many contexts (for example, the literature on repeat victimisation indeed suggests that we should expect the same addresses to be at a higher risk of being hit again). Even so, the point (no pun intended) is that "when the data has coincidence points, some statistical procedures will be severely affected. So it is always strongly advisable to check for duplicate points and to decide on a strategy for dealing with them if they are present" [156].

So how to perform this check? When creating the object above, we have already been warned this is an issue. But more generally we can check the duplication in a ppp object with the following syntax: it will return a logical indicating whether we have any duplication:

any(duplicated.ppp(bur_ppp))

## [1] TRUE

To count the number of coincidence points, we use the multiplicity() function from the spatstat package. This will return a vector of integers, with one entry for each observation in our dataset, giving the number of points that are identical to the point in question (including itself).

multiplicity(bur_ppp)

If you suspect many duplicates, it may be better to start asking how many there are. For this you can use:

sum(multiplicity(bur_ppp) > 1)

## [1] 915

That’s quite something! 915 points out of 1727 burglaries share the same coordinates. And we can plot this vector if we want to see the distribution of duplicates:

hist(multiplicity(bur_ppp))

In the case of crime, as we have hinted some, of this may be linked to the nature of crime itself. Hint: repeat victimisation. It is estimated that between 7 and 33% of all burglaries in several western and developed countries are repeat offenses [158]. Our 52% is quite high by those standards. In this particular situation, it likely is the result of geomasking (the introduction of small inaccuracies into their geocoded location) due to privacy reasons. Most police data stored in open data portals do not include the exact address of crime events. In the case of US police data, as noted by the documentation of CODE, it relies on addresses geocoded to the nearest hundred block. So say "370 Gorse Avenue" could be rounded down and geocoded to "300 Gorse Avenue". The same would happen to "307 Gorse Avenue", in the public release version of the data it would feature as "300 Gorse Avenue". In the UK there exist a pre-determined list of points so that each crime event gets "snapped" to its nearest one. So, the coordinates provided in the open data are not the exact locations of crimes. This process is likely inflating the amount of duplication we may observe, because each snap point or hundred block address might have many crimes near it, resulting in those crimes being geo-coded to the same exact location. So keep it in mind when analysing and working with this kind of data set that it may not be the same as working with the real locations. (If you are interested in the effects of this, read the paper by [119]).

What to do about duplicates in spatial point pattern analysis is not always clear. You could simply delete the duplicates, but of course that may ignore issues such as repeat victimisation. You could also use jittering, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space. This, again, may ignore things like repeat victimisation. Another alternative is to make each point "unique" and then attach the multiplicities of the points to the patterns as marks, as attributes of the points. Then you would need analytic techniques that take into account these marks.

If you were to be doing this for real, you would want access to the original source data, not this public version of the data, as well as access to the data producers to check that duplicates are not artificial, and then go for the latter solution suggested above (using multiplicities as marks once we adjust for any artificial inflation). We don’t have access to the source data, so for the sake of simplicity and convenience at this point and so that we can illustrate how spatstat works, we will first add some jittering to the data. For this we use rjitter().

The first argument for the function is the object, retry asks whether we want the algorithm to have another go if the jittering places a point outside the window (we want this so that we don’t have to lose points), and the drop argument is used to ensure we get a ppp object as a result of running this function (which we do). The argument nsim sets the number of simulated realisations to be generated. If we just want one ppp object, we need to set this to 1. You can also specify the scale of the perturbation; the points will be randomly relocated in a uniform way in a circle within this radius. The help file notes that there is a "sensible default", but it does not say much more about it, so you may want to set the argument yourself if you want to have more control over the results. Here we are specifying 1.5 metres to minimise the effects of this random displacement.

set.seed(200)
jitter_bur <- rjitter(bur_ppp, retry=TRUE, nsim=1, radius = 1.5, drop=TRUE)

And you can check duplicates with the resulting object:

sum(multiplicity(jitter_bur) > 1)

## [1] 0

Later we will demonstrate how to use the number of incidents in the same address as a mark or weight, so let’s create a new ppp object that includes this mark. We need to proceed in steps. First, we create a new ppp, just a duplicate of the one we have been using as a placeholder.

bur_ppp_m <- bur_ppp

Now we can assign to this ppp object the numerical vector created by the multiplicity function as the mark using the marks function.

marks(bur_ppp_m) <- multiplicity(bur_ppp)

If you examine your environment, you will be able to see the new object and that it has a 6th element, which is the mark we just assigned.