Section 7.5 Mapping intensity estimates with spatstat
Subsection 7.5.1 Basics of kernel density estimation
One of the most common methods to visualise hot spots of crime, when we are working with point pattern data, is to use kernel density estimation. Kernel density estimation involves applying a function (known as a "kernel") to each data point, which averages the location of that point with respect to the location of other data points. The surface that results from this model allows us to produce isarithmic maps. It is essentially a more sophisticated way of visualising intensity (points per area) than the tessellations we have covered so far.
Kernel density estimation maps are very popular among crime analysts. According to [159], 9 out of 10 intelligence professionals prefer them to other techniques for hot spot analysis. As compared to visualisations of crime that rely on point maps or thematic maps of geographic administrative units (such as LSOAs), isarithmic maps are considered best for location, size, shape and orientation of the hot spot. [160] have also suggested that this method produces some of the best prediction accuracy. The areas identified as hot spots by KDE (using historical data) tend to be the ones that better identify the areas that will have high levels of crime in the future. Yet, producing these maps (as with any map, really) requires you to take a number of decisions that will significantly affect the resulting product and the conveyed message. Like any other data visualisation technique, they can be powerful, but they have to be handled with great care.
Essentially this method uses a statistical technique (kernel density estimation) to generate a smooth continuous surface aiming to represent the intensity or volume of crimes across the target area. The technique, in one of its implementations (quartic kernel), is described in this way by [161]:
a fine grid is generated over the point distribution; a moving three-dimensional function of a specified radius visits each cell and calculates weights for each point within the kernel’s radius. Points closer to the centre will receive a higher weight, and therefore contribute more to the cell’s total density value; and final grid cell values are calculated by summing the values of all kernel estimates for each location
The values that we attribute to the cells in crime mapping will typically refer to the number of crimes within the area’s unit of measurement (intensity). Let’s produce one of these density maps, just with the default arguments of the
density.ppp() function:
ds <- density.ppp(jitter_bur)
class(ds)
## [1] "im"
The resulting object is of class
im, a pixel image, a grid of spatial locations with values attached to them. You can visualise this image as a map (using plot()), contour plots (using contour()), or as a perspective plot (using persp()).
par(mfrow=c(1,3))
plot(ds, main='Pixel image')
contour(ds, main = 'Contour plot')
persp(ds, main = 'Perspective plot')
The density function is estimating a kernel density estimate. The result of
density.ppp() is not a probability density. It is an estimate of the intensity function of the point process that generated the point pattern. The units of intensity are "points per unit area". Density, in this context, is nothing but the number of points per unit area. This method computes the intensity continuously across the study area, and the object returns a raster image. Depending what package you use to generate the estimation (e.g., spatstat, SpatialKDE, etc) you can set a different metric.
Subsection 7.5.2 Selecting the appropriate bandwidth for KDE
The defaults in
density.ppp() are not necessarily the optimal ones for your application. When doing spatial density analysis, there are a number of things you have to decide upon (e.g., bandwidth, kernel type, etc.).
To perform this analysis, we need to define the bandwidth of the density estimation, which basically determines the area of influence of the estimation. There is no general rule to determine the correct bandwidth; generally speaking, if the bandwidth is too small the estimate is too noisy, while if bandwidth is too high the estimate may miss crucial elements of the point pattern due to oversmoothing. The key argument to pass to the
density.ppp() is sigma=, which determines the bandwidth of the kernel. Sigma is the standard deviation of the kernel, and you can enter a numeric value or a function that computes an appropriate value for it.
A great deal of research has been developed to produce algorithms that help in selecting the bandwidth (see [156] and [162] for details). In the
spatstat package the functions bw.diggle(), bw.ppl(), and bw.scott() can be used to estimate the bandwidth according to difference methods (mean square error cross-validation, likelihood cross-validation, and Scott’s rule of thumb for bandwidth section in multi-dimensional smoothing, respectively). The help files recommend the use of the first two. These functions run algorithms that aim to select an appropriate bandwidth.
set.seed(200) # We set this seed to ensure you get this same results when running
# these algorithms
bw.diggle(jitter_bur)
bw.ppl(jitter_bur)
bw.scott(jitter_bur)
## sigma ## 2.912 ## sigma ## 73.06 ## sigma.x sigma.y ## 608.4 1447.6
Scott’s rule generally provides larger bandwidth than the other methods. Don’t be surprised if your results do not match exactly the ones we provide if you do not use the same seed, the computation means that they will vary every time we run these functions. Keep in mind these values will be expressed in the unit of the projected coordinate system you use (in this case metres). We can visualise how they influence the appearance of ’hot spots’:
set.seed(200)
par(mfrow=c(1,3))
plot(density.ppp(jitter_bur, sigma = bw.diggle(jitter_bur),edge=T),
main = "Diggle")
plot(density.ppp(jitter_bur, sigma = bw.ppl(jitter_bur),edge=T),
main="PPL")
plot(density.ppp(jitter_bur, sigma = bw.scott(jitter_bur) ,edge=T),
main="Scott")
When looking for clusters, [156] suggest the use of the
bw.ppl() algorithm because in their experience it tends to produce the more appropriate values when the pattern consists predominantly of tight clusters. But they also insist that if your purpose is to detect a single tight cluster in the midst of random noise, then the bw.diggle() method seems to work best.
All these methods use a fixed bandwidth: they use the same bandwidth for all locations. As noted by [156] a "fixed smoothing bandwidth is unsatisfactory if the true intensity varies greatly across the spatial domain, because it is likely to cause oversmoothing in the high-intensity areas and undersmoothing in the low-intensity areas" (p. 174). When using these methods in areas of high point density, the resulting bandwidth will be small, while in areas of low point density the resulting bandwidth will be large. This allows us to reduce smoothing in parts of the window where there are many points so that we capture more detail in the density estimates there, whilst at the same time increasing the smoothing in areas where there are less points.
With
spatstat we can use an adaptive bandwidth method, where (1) a fraction of the points are selected at random and used to construct a tessellation; (2) a quadrat counting estimate of the intensity is produced; and (3) the process is repeated a set number of times and an average is computed. The adaptive.density() function (in spatstat) takes as arguments the fraction of points to be sampled (f) and the number of samples to be taken (nrep). The logical argument verbose when set to FALSE (not the default) means we won’t see a message in the console informing of the steps in the simulation.
plot(adaptive.density(jitter_bur, f=0.1, nrep=30, verbose=FALSE),
main = "adaptive bandwidth method")
Which bandwidth to use? [163] is critical of the various algorithms we have introduced. He suggests that "they tend to produce large bandwidths, and are considered unsuitable for exploring local spatial patterns of the density distribution of crime." As we have seen with our example this is not necessarily the case, if anything both mean square error cross-validation and likelihood cross-validation produce small bandwidths in this example. When it comes to crime, [164] suggests that the "best choice would be to produce an interpolation that fits the experience of the department and officers who travel an area. Again, experimentation and discussions with beat officers will be necessary to establish which bandwidth choice should be used in future interpolations" (p. 10.22). This is generally accepted by other authors in the field. We see no damage in using these algorithms as a starting point and then, as suggested by most authors (see as well [163]), to engage in a discussion with the direct users (considering what the map will be used for) to set in a particular bandwidth that is not too spiky, nor oversmoothed.
Subsection 7.5.3 Selecting the smoothing kernel
Apart from selecting the bandwidth, we also need to specify the particular kernel we will use. In density estimation there are different types of kernel. This relates to the type of kernel drawn around each point in the process of counting points around each point. The use of these functions will result in slightly different estimations. They relate to the way we weigh points within the radius:
The normal distribution weighs all points in the study area, though near points are weighted more highly than distant points. The other four techniques use a circumscribed circle around the grid cell. The uniform distribution weighs all points within the circle equally. The quartic function weighs near points more than far points, but the fall off is gradual. The triangular function weighs near points more than far points within the circle, but the fall off is more rapid. Finally, the negative exponential weighs near points much more highly than far points within the circle and the decay is very rapid. [165] Chapter 10, p. 10.
Which one to use? [165] produces the following guidance: "The use of any of one of these depends on how much the user wants to weigh near points relative to far points. Using a kernel function which has a big difference in the weights of near versus far points (e.g., the negative exponential or the triangular) tends to produce finer variations within the surface than functions which weight more evenly (e.g., the normal distribution, the quartic, or the uniform); these latter ones tend to smooth the distribution more." However, [166] has argued that it does not make that much difference as long as the kernel is symmetrical. [161] suggests that in his experience most crime mappers prefer the quartic function, since it applies greater weight to crimes closer to the centre of the grid. The authors of the CrimeStat workbook [167], on the other hand, suggest that the choice of the kernel should be based in our theoretical understanding of the data generating mechanisms. By this they mean that the processes behind spatial dependence may be different according to various crime patterns and that this is something that we may want to take into account when selecting a particular function. They provide a table with some examples that may help you to understand what they mean. For example, for residential burglaries they recommend 1 mile with normally dispersed quartic or uniform interpolation, while for domestic violence, 0.1 mile and the tightly focused, negative exponential method [167].
The default kernel in
density.ppp() is gaussian. But there are other options. We can use the epanechnikov, quartic or disc. There are also further options for customisation. We can compare these kernels, whilst using the ppl algorithm for sigma (remember, we set seed at 200 for this particular value):
par(mfrow=c(2,2),mar = c(0.1, 0.1, 0.6, 0.1))
# Gaussian kernel function (the default)
plot(density.ppp(jitter_bur, sigma = bw.ppl(jitter_bur), edge=TRUE), main="Gaussian")
# Epanechnikov (parabolic) kernel function
plot(density.ppp(jitter_bur, kernel = "epanechnikov", sigma = bw.ppl(jitter_bur),
edge=TRUE), main="epanechnikov")
# Quartic kernel function
plot(density.ppp(jitter_bur, kernel = "quartic", sigma = bw.ppl(jitter_bur), edge=TRUE),
main="Quartic")
plot(density.ppp(jitter_bur, kernel = "disc", sigma = bw.ppl(jitter_bur), edge=TRUE),
main="Disc")
You can see how the choice of the kernel has less of an impact than the choice of the bandwidth. In practical terms, choice of the specific kernel is of secondary importance to the bandwidth.
Equally, you may have noticed an argument in
density.ppp() we have not discussed: edge = TRUE. What is this? When estimating the density we are not observing points outside our window. In the case of an island like Manhattan, perhaps this is not an issue. But often when computing estimates for a given jurisdiction, we need to think about what is happening in the jurisdiction next door to us. Our algorithm does not see the crimes in jurisdictions neighbouring ours, and so the estimations around the edges of our boundaries receive fewer contributions in the computation. By setting edge = TRUE we are effectively applying a correction that tries to compensate for this negative bias.
Subsection 7.5.4 Weighted kernel density estimates
In the past examples we used the un-weighted version of the estimation methods, but if the points are weighted (like for example by multiplicity) we can use weighted versions of these methods. We simply need to pass an argument that specifies the weights.
set.seed(200)
plot(density.ppp(bur_ppp_m, sigma = bw.ppl(bur_ppp_m), edge = T,
kernel = "quartic", weights = marks(bur_ppp_m)),
main = "Weighted Points")
Subsection 7.5.5 Better KDE maps using leaflet
The plots produced by
spatstat are good for initial exploration, but if you want to present this to an audience, you may want to use one of the various packages we have covered so far for visualising spatial data. Here we show how to use leaflet for plotting densities. Hot spot analysis often has an exploratory character, and the interactivity provided by leaflet is particularly helpful, since it allows us to zoom in and out and check the details of the areas with high intensities.
Often it is convenient to use a basemap to provide context. In order to do that, we first need to turn the image object generated by the
spatstat package into a raster object, a more generic format for raster image used in R. Remember rasters from the first chapter? Now we finally get to use them a bit! First, we store the density estimates produced by density.ppp() into an object (which will be of class im, a two-dimensional pixel image), and then using the raster() function from the raster package, we generate a raster object.
set.seed(200)
dmap1 <- density.ppp(jitter_bur,
kernel = "quartic",
sigma = bw.ppl(jitter_bur),
edge=T)
# create raster object from im object
r1 <- raster(dmap1)
# plot the raster object
plot(r1, main = "Raster")
We can remove the areas with very low density of burglaries for better reading of the map:
#remove very low density values
r1[r1 < 0.00010 ] <- NA
#multiply by 1000 to re-express intensity in terms of burglaries per km
values(r1) <- values(r1)*1000
Now that we have the raster, we can add it to a basemap. Two-dimensional
RasterLayer objects (from the raster package) can be turned into images and added to Leaflet maps using the addRasterImage() function. The addRasterImage() function works by projecting the RasterLayer object to EPSG:3857 and encoding each cell to an RGBA colour, to produce a PNG image. That image is then embedded in the map widget. This is only suitable for small to medium-sized rasters.
It’s important that the
RasterLayer object is tagged with a proper coordinate reference system, so that leaflet can then reproject to EPSG:3857. Many raster files contain this information, but some do not. Here is how you would tag a raster layer object r1 which contains the EPSG:32118 data:
#make sure we have right CRS, which in this case is the New York 32118
crs(r1) <- "EPSG:32118"
The way we assign the CRS above is using the now most accepted PROJ6 format. You can find how these different formats apply to each coordinate system in the online archives.
Now we are ready to plot, although first we will create a more suitable colour palette and, since the values in the scale are quite close to 0, will multiply by a constant, to make the labels more readable. For this we use the
values() function from the raster package that allows us to manipulate the values in a raster object:
pal <- colorNumeric("Reds", values(r1),
na.color = "transparent") #Create colour palette
#and then make map!
nyc_map <- leaflet() %>%
addTiles() %>%
addProviderTiles(providers$Stamen.TonerLite) %>%
addRasterImage(r1, colors = pal, opacity = 0.7) %>%
addLegend(pal = pal, values = values(r1),
title = "Burglary map")
# print map
nyc_map
And there you have it. Perhaps those familiar with Manhattan have some guesses as to what may be going on there?
Subsection 7.5.6 Problems with KDE
There are at least two problems with these maps. First, they are just a convenient visualisation to represent the location of crime incidents. But they are not adjusting by the population at risk. If all you care about is where most crimes are taking place, that may then be fine. But if you are interested in relative risk, intensity adjusting by population at risk, then you are bound to be misled by these methods, which usually will plot very well the distribution of the population at risk. Second, the Poisson process, as we saw, creates what may appear as areas with more points than others, even though it is a random process. In other words, you may just be capturing random noise with these maps.
There are a number of ways to adjust for background population at risk when producing density estimates. Many of these methods have been developed in spatial epidemiology. In criminology, the population at risk is generally available only at more aggregated (at some census geography level) than in many of those applications (where typically you have points representing individuals without a disease being contrasted with points representing individuals with a disease). So, although one could extract the centroids of these polygons representing census tracts or districts as points and mark them with a value (say, population count or number of housing units), there isn’t a great advantage to go to all this trouble. In fact by locating all the mass of the population at risk in the centroids of the polygons, you are introducing some distortion. So if you want to adjust for the population at risk, it is probably better to map out rates in the ways we described in previous chapters. We just need to be aware, then, that KDE surfaces in essence are simply mapping counts of crime. We therefore need to interpret density maps with this very significant caveat.
The other problem we had is that a Poisson process under CSR can still result in a map that one could read as identifying hot spots, whereas all we have is noise. Remember the
poisson_ppp object we created earlier, which demonstrated complete spatial randomness? Well, if we were that way inclined, we could easily compute a KDE for this:
set.seed(200)
plot(density.ppp(poisson_ppp, kernel = "quartic", sigma = bw.diggle(poisson_ppp),
edge=T), main = "KDE for random poisson point process")
And there you go, you could come to the view looking at this density map that the intensity is driven by some underlying external risk factor when in fact the process is purely random. In subsequent chapters we will learn about methods that can be used to ensure we don’t map out hot spots that may be the product of random noise.
Subsection 7.5.7 Other packages for KDE
In this chapter we have focused on
spatstat, for it provides a comprehensive suite of tools for working with spatial point patterns. The snippets of code and examples we have used here only provide a superficial glimpse into the full potential of this package. But spatstat is by no means the only package that can be used for generating spatial density estimates. We need to mention at least two alternatives: SpatialKDE and sparr.
SpatialKDE implements kernel density estimation for spatial data with all the required settings, including the selection of bandwidth, kernel type, and cell size (which was computed automatically by density.ppp()). It is based in an algorithm used in QGIS, the best free GIS system relying primarily on a point-and-click interface. Unlike spatstat, SpatialKDE can work with sf objects. But it also finds it more convenient to have the data on a projected coordinate systems, so that it can compute distances on metres rather than with decimal degrees.
The package
sparr focuses on methods to compare densities and thus estimate relative risk, even at a spatio-temporal scale. It is written in such a way as to allow compatibility with spatstat. The basic functionality allows to estimate KDE maps through the bivariate.density() function. This function can take a ppp object and produce a KDE map. In some ways it gives the analyst more flexibility than density.ppp(). It allows to modulate the resolution of the grid, two different methods to control for edge effects, the type of estimate produced (intensity or density, that integrates to 1), etc. This package also introduces various methods for estimating the bandwidth different to those existing in spatstat. Critically, it introduces a function that allows for spatio-temporal kernel density estimation: spattemp.density(). For more details, see the help files for sparr, and the tutorial by [162].
