Section 10.5 Scan statistics
Subsection 10.5.1 The DCluster package
As we mentioned at the outset of the chapter, a number of alternative approaches using a different philosophy to those that are based on the idea of LISA have been developed within spatial epidemiology. Rather than focusing on local dependence, these tests focus on assessing homogeneity ([Marshall, 1991]):
The null hypothesis when clustering is thought to be due to elevated local risk is usually that cases of disease follow a non-homogeneous Poisson process with intensity proportional to population density. Tests of the null Poisson model can be carried out by comparing observed and expected count (p. 423).
Here we will look at two well-known approaches with a focus on their implementation in the
DCluster package: Openshaw’s GAM and Kulldorff’s scan statistic. The package DCluster ([Gomez-Rubio et al., 2005]) was one of the first to introduce tools to assess spatial disease clusters to R. This package precedes sf and it likes the data to be stored in a data frame with at least four columns: observed number of cases, expected number of cases, and longitude and latitude. It is a bit opinionated in the way it is structured. DCluster, like many packages that are influenced by the epidemiological literature, takes as key inputs the observed cases (in criminology this will be criminal occurrences) and expected cases. In epidemiology, expected cases are often standardised to account for population heterogeneity. In crime analysis, this is not as common (although in some circumstances it should be!).
Straight away we see the problem we encounter here. The criminology of place emphasises small places, but for this level of aggregation (street segments or micro-grid cells) we won’t have census variables indexing the population that we could easily use for purposes of standardisation. This means we need to (1) compromise and look for clusters at a level of aggregation for which we have some population measure that we could use for standardising or (2) get creative and try to find other ways of solving this problem. Increasingly, environmental criminologists are trying to find ways of taking this second route (see [Andresen et al., 2021]). For this illustration and for convenience reasons, we will take the first one. Although it is important to warn we should try to stick whenever possible to the level of the spatial process we are interested in exploring.
For this example, the expected counts will be computed using the overall incidence ratio (i.e., total number of cases divided by the total population; in this case the number of residential dwellings in the area). We will use the unit of analysis of neighbourhoods operationalised as a commonly used census unit of geography in the UK, the Lower Super Output Area (LSOA).
Let’s compute the expected count by dividing the total number of burglaries by the total number of dwellings (if burglaries were evenly distributed across all dwellings in our study area) and then multiplying this with the number of dwellings in each LSOA.
# calculate rate per dwelling in whole study area
even_rate_per_dwelling <- sum(burglary_lsoa$burglary) / sum(burglary_lsoa$dwellings)
# expected count for each LSOA times rate per number of dwellings
burglary_lsoa <- burglary_lsoa %>%
mutate(b_expected = dwellings * even_rate_per_dwelling)
This package requires as well the centroids for the polygons we are using. To extract this, we can use
st_centroid() function from the sf package. If we want to extract the coordinates individually, we also need to use the st_coordinates() function. All the GEOS functions underlying sf need projected coordinates to work properly, so we need to ensure an adequate projection (which we did when we loaded the data at the outset of the chapter):
# Get coordinates for centroids
lsoa_centroid <- st_centroid(burglary_lsoa) %>%
st_coordinates(lsoa_centroid)
# Place the coordinates as vectors in our sf data frame
burglary_lsoa$x <- lsoa_centroid[,1]
burglary_lsoa$y <- lsoa_centroid[,2]
For convenience, we will just place these four columns in a new object. Although this is not clear from the documentation, in our experience some of the code will only work if you initialise the data using these column names:
lsoas <- data.frame(Observed = burglary_lsoa$burglary,
Expected = burglary_lsoa$b_expected,
Population = burglary_lsoa$dwellings,
x = burglary_lsoa$x,
y = burglary_lsoa$y)
DCluster includes several tests for assessing the general heterogeneity of the relative risks. A chi-squared test can be run to assess global differences between observed and expected cases.
chtest <- achisq.test(Observed~offset(log(Expected)),
as(lsoas, "data.frame"),
"multinom",
999)
chtest
## Chi-square test for overdispersion ## ## Type of boots.: parametric ## Model used when sampling: Multinomial ## Number of simulations: 999 ## Statistic: 4612 ## p-value : 0.001
It also includes the Pottoff and Withinghill test of homogeneity. The alternative hypothesis of this test is that the observed cases are distributed following a negative binomial distribution ([Bivand et al., 2013]).
pwtest <- pottwhitt.test(Observed~offset(log(Expected)),
as(lsoas, "data.frame"),
"multinom", 999)
oplus <- sum(lsoas$Observed)
1 - pnorm(pwtest$t0, oplus * (oplus - 1), sqrt(2 * 100 * oplus * (oplus - 1)))
## [1] 0
Detail in those two tests and the used code is available from [Bivand et al. (2013)]. In what follows we introduce two of the approaches implemented here for the detection of clusters.
Subsection 10.5.2 Openshaw’s GAM
Stan Openshaw is a retired British geography professor with an interest in geocomputation. We already encountered him earlier, for he wrote one of the seminal papers on the modifiable areal unit problem. In 1987, together with several colleagues, he published his second most highly cited paper introducing geographical analysis machine (GAM). This was a technique developed to work with point pattern data with the aim of finding clusters ([Openshaw et al., 1987]). The machine is essentially an automated algorithm that tries to assess whether "there is an excess of observed points within \(x\) km of a specific location" (p. 338). A test statistic is computed for this circular search area, and some form of inference can then be made.
The technique requires a grid and will draw circles centred in the intersection of the cells of this grid. From each point of the grid backdrop, a radial distance is calculated assuming a Poisson distribution. The ratio of events to candidates is counted up and then the test is done. What we are doing is to observe case counts in overlapping circles and trying to identify potential clusters among these, by virtue of running separate significance test for each of the circles individually. The statistically trained reader may already note a problem with this. In a way, this is similar to tests for quadrat counting, only here we have overlapping "quadrats" and we have many more of them ([Kulldorff and Nagarwalla, 1995]).
The
opgam() function from DCluster does this for you. You can either specify the grid passing the name of the object containing this grid using the thegrid argument or specify the step size for computing this grid. For this you simply introduce a value for the step argument. You also need to specify the search radius and can specify an explicit significance level of the test performed with the alpha argument. The radius is typically 5 to 10 times the lattice spacing. It is worth mentioning that although one could derive this test using asymptotic theory, DCluster uses a bootstrap approach.
gam_burglary <- opgam(data = as(lsoas, "data.frame"),
radius = 50, step = 10, alpha = 0.002)
Once
opgam has found a solution, we could plot either the circles or the centre of these circles. We can use the following code to do the map, which will be shown in the next section.
gam_burglary_sf <- st_as_sf(gam_burglary, coords = c("x", "y"), crs = 27700)
gam_burglary_sf <- st_transform(gam_burglary_sf, crs = 4326)
tmap_mode("plot")
map_gam <- tm_shape(burglary_lsoa) + tm_borders(alpha = 0.3) +
tm_shape(gam_burglary_sf) + tm_dots(size=0.1, col = "red") +
tm_layout(main.title = "Openshaw's GAM", main.title.size=1.2)
One of the problems with this method is that the tests are not independent; very similar clusters (i.e., most of their regions are the same) are tested ([Bivand et al., 2013]). Inference then is a bit compromised and the approach has received criticism over the years for this failure to solve the multiple testing problem: "any Bonferroni type of procedure to adjust for multiple testing is futile due to the extremely large number of dependent tests performed" ([Kulldorff and Nagarwalla, 1995]). For these reasons and others (see [Marshall, 1991]), it is considered hard to ask meaningful statistical questions from a GAM. Notwithstanding this, this technique, as noted by [Bivand et al. (2013)], can still be helpful as part of the exploration of your data, as long as you are aware of its limitations.
Subsection 10.5.3 Kulldorff’s Scan Statistic
More popular has been the approach developed by [Kulldorff and Nagarwalla (1995)]. This test could be used for aggregated or point pattern data. It aims to overcome the limitations from previous attempts to assess clusters of disease. Unlike previous methods, the scan statistic constructs circles of varying radii ([Kulldorff and Nagarwalla, 1995]):
Each of the infinite number of circles thus constructed defines a zone. The zone defined by a circle consists of all individuals in those cells whose centroids lie inside the circle and each zone is uniquely identified by these individuals (p. 802).
The test statistic then compares the relative risk within the circle to the relative risk outside the circles (for more details, see [Kulldorff and Nagarwalla, 1995] and [Kulldorff, 1997]).
The scan statistic was implemented by Kulldorff in the SatScan™ software, freely available from www.satscan.org. This software extends the basic idea to various problems, including spatio-temporal clustering. A number of R packages use this scan statistic (
smerc, spatstat, SpatialEpi are some of the best known), although the implementation is different across these packages. There is also a package (rsatscan) that can allow us to interface with SatScan™ via R. Here we will continue focusing on DCluster for parsimonious reasons. The vignettes for the other packages provide a useful platform for exploring those other solutions.
For this we return to the
opgam() function. There are some key differences in the arguments. The icluster argument, which we didn’t specify previously, needs to make clear now we are using the clustering function for Kulldorff and Nagarwalla scan statistic. R is the number of bootstrap permutations to compute. The model argument allows us to specify the model to generate the random observations. [Bivand et al. (2013)] suggest the use of a negative binomial model ("negbin") to account for over-dispersion which may result from unaccounted spatial autocorrelation. And mle specifies parameters that are needed for the bootstrap calculation. These can be computed with the DCluster::calculate.mle() function, which simply takes as arguments the data and the type of model. When illustrating opgam earlier, we used the step method to draw the grid; in this case we will explicitly provide a grid based on our data points.
mle <- calculate.mle(as(lsoas, "data.frame"), model = "negbin")
thegrid <- as(lsoas, "data.frame")[, c("x", "y")]
kns_results <- opgam(data = as(lsoas, "data.frame"),
thegrid = thegrid, alpha = 0.02, iscluster = kn.iscluster,
fractpop = 0.05, R = 99, model = "negbin",
mle = mle)
Once the computation is complete, we are ready to plot the data and to compare the result to those obtained with the GAM.
kulldorff_scan <- st_as_sf(kns_results,
coords = c("x", "y"),
crs = 27700)
kulldorff_scan <- st_transform(kulldorff_scan, crs = 4326)
map_scan <- tm_shape(burglary_lsoa) + tm_borders(alpha = 0.3) +
tm_shape(kulldorff_scan) + tm_dots(size=0.1, col = "red") +
tm_layout(main.title = "Kulldorff Scan", main.title.size=1.2)
tmap_arrange(map_gam, map_scan)
Whereas the GAM picks up a series of locations in East and North Manchester known by their high level of concentrated disadvantage, the Kulldorff scan only identifies the clusters around the city centred on areas of student accommodation south of the University of Manchester.
