Section 3.6 Smoothing rates: adjusting for small sample noise
In previous sections we discussed how to map rates. It seems a fairly straightforward issue: you calculate a rate by dividing your numerator (e.g., number of crimes) by an appropriately selected denominator (e.g., daytime population). You get your variable with the relevant rate and you map it using a choropleth map. However, things are not always that simple.
Rates are funny animals. [Gelman and Price (1999)] go so far as to suggest that all maps of rates are misleading. The problem at hand is well known in spatial epidemiology: "plotting observed rates can have serious drawbacks when sample sizes vary by area, since very high (and low) observed rates are found disproportionately in poorly-sampled areas" ([Gelman and Price, 1999], p. 3221). There is associated noise for those areas where the denominators give us a small sample size. And it is hard to solve this problem.
Let’s illustrate with an example. We are going to use historic data on homicide across US counties. The dataset was used as the basis for the study by [Baller et al. (2001)]. It contains data on homicide counts and rates for various decades across the US, as well as information on structural factors often thought to be associated with violence. The data is freely available through the webpage of Geoda, a clever point-and-click interface developed by Luc Anselin (a spatial econometrician and coauthor of the above paper) and his team, to make spatial analysis accessible. It is also available as one of the datasets in the
geodaData package.
To read data available in a library we have loaded, we can use the
data() function. If we check the class() of this object, we will see it was already stored in geodaData as a sf object.
data("ncovr")
class(ncovr)
## [1] "sf" "data.frame"
Let’s look at the "ncovr" data. We can start by looking at the homicide rate for 1960.
summary(ncovr$HR60)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.00 0.00 2.78 4.50 6.89 92.94
We can see that the county with the highest homicide rate in the 1960s had a very high rate of homicides per 100,000 individuals. Just to put it into context, in the UK the homicide rate is about 0.92 per 100,000 individuals. Where is that place? I can tell you it is a place called Borden. Check it out:
borden <- filter(ncovr, NAME == "Borden")
borden$HR60
## [1] 92.94
Borden county is in Texas. You may be thinking "Texas Chainsaw Massacre" perhaps? No, not really. Ed Gein, who inspired the film, was based and operated in Wisconsin. Borden’s claim to fame is rather more prosaic: it was named after Gail Borden, the inventor of condensed milk. So, what’s going on here? Why do we have a homicide rate in Borden that makes it look like a war zone? Is it that it is only one of the six counties where alcohol is banned in Texas?
To get to the bottom of this, we can look at the variable
HC60 which tells us about the homicide count in Borden:
borden$HC60
## [1] 1
What? A total homicide count of 1. How can a county with just one homicide have a rate that makes it look like the most dangerous place in the US? To answer this, let’s look at the population of Borden county in the 1960s, contained in the
PO60 variable.
borden$PO60
## [1] 1076
Well, there were about 1000 people living there. It is among some of the least populous counties in our data:
summary(ncovr$PO60)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 208 9417 18408 57845 39165 7781984
If you contrast that population count with the population of the average county in the US, that’s tiny. One homicide in such a small place can end up producing a big rate. Remember that the rate is simply dividing the number of relevant events by the exposure variable (in this case, population) and multiplying by a constant (in this case, 100,000 since we expressed crime rates in those terms). Most times Borden looks like a very peaceful place:
borden$HR70
borden$HR80
borden$HR90
## [1] 0 ## [1] 0 ## [1] 0
It has a homicide rate of 0 in most decades. But it only takes one homicide and, bang, it goes top of the league. So a standard map of rates is bound to be noisy. There is the instability that is introduced by virtue of having areas that may be sparsely populated and in which one single event, like in this case, will produce a very noticeable change in the rate.
In fact, if you look at the counties with the highest homicide rate in the "ncovr" dataset, you will notice all of them are places like Borden, areas that are sparsely populated, not because they are that dangerous, but because of the instability of rates. Conversely, the same happens with those places with the lowest rate. They tend to be areas with a very small sample size.
This is a problem that was first noted by epidemiologists doing disease mapping. But a number of other disciplines have now noted this and used some of the approaches developed by public health researchers that confronted this problem when producing maps of disease (aside: techniques and approaches used by spatial epidemiologists are very similar to those used by criminologists, in case you ever think of changing careers or need inspiration for how to solve a crime analysis problem).
One way of dealing with this is by smoothing or shrinking the rates. This basically, as the word implies, aims for a smoother representation that avoids hard spikes associated with random noise. There are different ways of doing that. Some ways use a non-spatial approach to smoothing, using something called a empirical Bayesian smoother. How does this work? This approach takes the raw rates and tries to "shrink" them towards the overall average. What does this mean? Essentially, we compute a weighted average between the raw rate for each area and the global average across all areas, with weights proportional to the underlying population at risk. What this procedure does is adjusting considerably (brought closer to the global average) the rates of smaller areas (those with a small population at risk), whereas the rates for the larger areas will barely change.
Here we are going to introduce the approach implemented in
DCluster, a package developed for epidemiological research and detection of clusters of disease. Specifically, we can implement the empbaysmooth() function which creates a smooth relative risks from a set of expected and observed number of cases using a Poisson-Gamma model as proposed by [Clayton and Kaldor (1987)]. The function empbaysmooth() expects two parameters, the expected value and the observed value. Let’s define them first.
#First we define the observed number of cases
ncovr$observed <- ncovr$HC60
#To compute the expected number of cases through indirect standardisation
#we need the overall incidence ratio
overall_incidence_ratio <- sum(ncovr$HC60)/sum(ncovr$PO60)
#The expected number of cases can then be obtained by multiplying the overall
#incidence rate by the population
ncovr$expected <- ncovr$PO60 * overall_incidence_ratio
With this parameters we can obtain the raw relative risk:
ncovr$raw_risk <- ncovr$observed / ncovr$expected
summary(ncovr$raw_risk)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.000 0.000 0.614 0.994 1.520 20.510
And then estimate the smoothed relative risk:
res <- empbaysmooth(ncovr$observed, ncovr$expected)
In the new object we generated, which is a list, you have an element which contains the computed rates. We can add those to our dataset:
ncovr$H60EBS <- res$smthrr
summary(ncovr$H60EBS)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.234 0.870 0.936 0.968 1.035 2.787
We can observe that the dispersion narrows significantly and that there are fewer observations with extreme values once we use this smoother.
Instead of shrinking to the global relative risk, we can shrink to a relative rate based on the neighbours of each county. Shrinking to the global risk ignores the spatial dimension of the phenomenon being mapped out and may mask existing heterogeneity. If instead of shrinking to a global risk, we shrink to a local rate, we may be able to take unobserved heterogeneity into account. [Marshall (1991)] proposed a local smoother estimator in which the crude rate is shrunk towards a local, "neighbourhood", rate. To compute this, we need the list of neighbours that surround each county (we will discuss this code in Chapter 11, so for now just trust we are computing the rate of the areas that surround each country):
ncovr_sp <- as(ncovr, "Spatial")
w_nb <- poly2nb(ncovr_sp, row.names=ncovr_sp$FIPSNO)
eb2 <- EBlocal(ncovr$HC60, ncovr$PO60, w_nb)
ncovr$HR60EBSL <- eb2$est * 100000
We can now plot the maps and compare them:
tmap_mode("plot")
current_style <- tmap_style("col_blind")
map1 <- tm_shape(ncovr) +
tm_fill("HR60", style="quantile",
title = "Raw rate",
palette = "Reds") +
tm_layout(legend.position = c("left", "bottom"),
legend.title.size = 0.8, legend.text.size = 0.5)
map2 <- tm_shape(ncovr) +
tm_fill("HR60EBSL", style="quantile",
title = "Local Smooth",
palette = "Reds") +
tm_layout(legend.position = c("left", "bottom"),
legend.title.size = 0.8, legend.text.size = 0.5)
tmap_arrange(map1, map2)
Notice that the quantiles are not the same, so that will make your comparison difficult. Let’s look at a boxplot of these variables. In the map of raw rates, we have the most variation.
#Boxplots with base R graphics
boxplot(ncovr$HR60, ncovr$HR60EBSL,
main = "Homicide per 100,000 people",
names = c("Raw rate",
"Local base smoother"),
ylab = "Rates")
The range for the raw rates is nearly 93. Much of the variation in observed homicide rates by county is attributable to statistical noise due to the small number of (observed and expected) homicides in low-population counties. Because of this noise, a disproportionate fraction of low-population counties are observed to have extremely high (or low) homicide rates when compared to typical counties in the United States. With smoothing, we reduce this problem, and if you contrast the maps you will see how this results in a clearer and smoother spatial pattern for the rate that is estimated borrowing information from their neighbours.
So to smooth or not to smooth? Clearly we can see how smoothing stabilises the rates and removes noise. But as [Gelman and Price (1999)] suggests this introduces other artifacts and autocorrelation into our estimates. Some people are also not too keen on maps of statistically adjusted estimates. Yet, the conclusions one can derive from mapping raw rates (when the denominator varies significantly and we have areas with small sample size) means that smoothing is often a preferable alternative ([Waller and Gotway, 2004]). The problem we have with maps of estimates is that we need information about the variability and it is hard to map this out in a convenient way ([Gelman and Price, 1999]). [Lawson (2006)], in relation to the similar problem of disease mapping, suggests that "at the minimum any map of relative risk for a disease should be accompanied with information pertaining to estimates of rates within each region as well as estimates of variability within each region" (p. 38), whereas "at the other extreme it could be recommended that such maps be only used as a presentational aid, and not as a fundamental decision-making tool" (p. 38).
