Head/tail breaks in the classInt package

There are far more ordinary people (say, 80 percent) than extraordinary people (say, 20 percent), and this is often characterized by the 80/20 principle, based on the observation made by the Italian economist Vilfredo Pareto in 1906 that 80% of land in Italy was owned by 20% of the population. A histogram of the data values for these phenomena would reveal a right-skewed or heavy-tailed distribution. How to map the data with the heavy-tailed distribution?

Abstract

This vignette discusses the implementation of the “Head/tail breaks” style (Jiang (2013)) in the classIntervals() function from the classInt package. A step-by-step example is presented to clarify the method. A case study using spData::afcon is also included, making use of additional packages such as sf.

Introduction

Head/tail breaks is a classification scheme introduced by Jiang (2013) for revealing hierarchies in data with many small values and a few large ones. The related ht-index measures the depth of this hierarchy, rather than being another name for the classification method.

The distributions considered here are strongly right-skewed, with a minority of large values in the head and a majority of small values in the tail. This imbalance can be expressed as “far more small things than large things”.

Power-law and lognormal distributions are examples of heavy-tailed distributions. The head/tail breaks literature also considers exponential distributions under its broader description of this imbalance. Nature, society and finance (Vasicek (2002)) provide examples of rare and extreme events. Taleb (2008) discusses the impact of unexpected extreme events, offering another reason to pay attention to low-frequency observations.

library(classInt)

# 1. Characterization of heavy-tail distributions----
set.seed(1234)
# Pareto distribution a=1 b=1.161 n=1000
sample_par <- 1 / (1 - runif(1000))^(1 / 1.161)
opar <- par(no.readonly = TRUE)
par(mar = c(2, 4, 3, 1), cex = 0.8)
plot(
  sort(sample_par, decreasing = TRUE),
  type = "l",
  ylab = "F(x)",
  xlab = "",
  main = "80/20 principle"
)
abline(
  h = quantile(sample_par, .8),
  lty = 2,
  col = "red3"
)
abline(
  v = 0.2 * length(sample_par),
  lty = 2,
  col = "darkblue"
)
legend(
  "topleft",
  legend = c("F(x): p80", "x: Top 20%"),
  col = c("red3", "darkblue"),
  lty = 2,
  cex = 0.8
)

hist(
  sample_par,
  n = 100,
  xlab = "",
  main = "Histogram",
  col = "grey50",
  border = NA,
  probability = TRUE
)
par(opar)

Breaking method

The method itself consists of a four-step process performed recursively until a stopping condition is satisfied. Given a vector of values \(v = (a_1, a_2, ..., a_n)\), the process can be described as follows:

1. Compute the arithmetic mean, \(\mu = \frac{1}{n}\sum_{i=1}^{n} a_i\).
2. Break \(v\) into the \(tail\) and the \(head\): \(tail = \{ a_x \in v | a_x \lt \mu \}\) \(head = \{ a_x \in v | a_x \gt \mu \}\).
3. Check whether the proportion of observations in the head is at most the chosen threshold: \(\frac{|head|}{|v|} \le threshold\)
4. If step 3 is TRUE and the head contains at least two observations, replace \(v\) with \(head\) and repeat. Otherwise, stop.

Only values strictly above the mean enter the next iteration. Values equal to the mean are not part of the head. Stopping does not require an exact 50/50 split.

A threshold of 40% provides a practical stopping rule. It is not a statistical test for whether a distribution is heavy-tailed.

The final breaks are the vector of consecutive \(\mu\):

Step-by-step example

We reproduce here the pseudo-code¹ as per Jiang (2019):

Recursive function Head/tail Breaks:
 Rank the input data from the largest to the smallest
 Break the data into the head and the tail around the mean;
 // the head for those above the mean
 // the tail for those below the mean
 While (head <= 40%):
 Head/tail Breaks (head);
End Function

The following R example uses prop < thr, whereas the pseudocode above uses head <= 40%. At exactly 40%, the R loop stops and the pseudocode continues. This distinction does not affect the example below.

opar <- par(no.readonly = TRUE)
par(mar = c(2, 2, 3, 1), cex = 0.8)
var <- sample_par
thr <- .4
brks <- c(min(var), max(var)) # Initialise with min and max

sum_table <- data.frame(
  iter = 0,
  mu = NA,
  prop = NA,
  n_var = NA,
  n_head = NA
)
# Pars for chart
limchart <- brks
# Iteration
for (i in 1:10) {
  mu <- mean(var)
  brks <- sort(c(brks, mu))
  head <- var[var > mu]
  prop <- length(head) / length(var)
  stopit <- prop < thr & length(head) > 1
  sum_table <- rbind(
    sum_table,
    c(i, mu, prop, length(var), length(head))
  )
  hist(
    var,
    main = paste0("Iter ", i),
    breaks = 50,
    col = "grey50",
    border = NA,
    xlab = "",
    xlim = limchart
  )
  abline(v = mu, col = "red3", lty = 2)
  ylabel <- max(hist(var, breaks = 50, plot = FALSE)$counts)
  labelplot <- paste0("PropHead: ", round(prop * 100, 2), "%")
  text(
    x = mu,
    y = ylabel,
    labels = labelplot,
    cex = 0.8,
    pos = 4
  )
  legend(
    "right",
    legend = paste0("mu", i),
    col = c("red3"),
    lty = 2,
    cex = 0.8
  )
  if (isFALSE(stopit)) {
    break
  }
  var <- head
}
par(opar)

The head proportion varies across iterations and reaches 50% in the fourth, which stops the loop. It does not increase monotonically.

The break vector includes the original minimum and maximum, plus the means computed during the iterations. Here the original data are sample_par, since var is replaced by the head inside the loop.

Implementation in the classInt package

The implementation in classIntervals() reproduces these results:

ht_sample_par <- classIntervals(sample_par, style = "headtails")
brks == ht_sample_par$brks

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

As stated in Jiang (2013), the number of breaks is determined by the data. However, the thr parameter can help adjust the final number. A lower value of thr can yield fewer breaks, while a larger thr can increase the number if the underlying distribution follows the “far more small things than large things” principle.

opar <- par(no.readonly = TRUE)
par(mar = c(2, 2, 2, 1), cex = 0.8)

pal1 <- c("wheat1", "wheat2", "red3")
# Minimum: single break
print(paste("number of breaks", length(classIntervals(sample_par, style = "headtails", thr = 0)$brks - 1)))
plot(
  classIntervals(sample_par, style = "headtails", thr = 0),
  pal = pal1,
  main = "thr = 0"
)

# Two breaks
print(paste("number of breaks", length(classIntervals(sample_par, style = "headtails", thr = 0.2)$brks - 1)))
plot(
  classIntervals(sample_par, style = "headtails", thr = 0.2),
  pal = pal1,
  main = "thr = 0.2"
)

# Default breaks: 0.4
print(paste("number of breaks", length(classIntervals(sample_par, style = "headtails")$brks - 1)))
plot(classIntervals(sample_par, style = "headtails"),
  pal = pal1,
  main = "thr = Default"
)

# Maximum breaks
print(paste("number of breaks", length(classIntervals(sample_par, style = "headtails", thr = 1)$brks - 1)))
plot(
  classIntervals(sample_par, style = "headtails", thr = 1),
  pal = pal1,
  main = "thr = 1"
)
par(opar)

For this example, even thr = 0 retains the mean as an internal break. The returned vector also contains the minimum and maximum, so this gives two classes, not one.

Case study

Jiang (2013) states that "the new classification scheme is more natural than the natural breaks in finding the groupings or hierarchy for data with a heavy-tailed distribution." (p. 482), referring to Jenks’ natural breaks method. In this case study, we compare headtails vs. fisher, which is the alias for the Fisher-Jenks algorithm and is always preferred to the jenks style (see ?classIntervals). For this example, we will use the afcon dataset from the spData package, plus some additional spatial information to create the data visualization.

library(spData)
data(afcon, package = "spData")

Let’s have a look at the top 10 values and the distribution of the variable totcon (index of total conflict 1966-78):

# Top10
knitr::kable(head(afcon[order(afcon$totcon, decreasing = TRUE), c("name", "totcon")], 10))

opar <- par(no.readonly = TRUE)
par(mar = c(4, 4, 3, 1), cex = 0.8)
hist(afcon$totcon,
  n = 20,
  main = "Histogram",
  xlab = "totcon",
  col = "grey50",
  border = NA,
)
plot(
  density(afcon$totcon),
  main = "Distribution",
  xlab = "totcon",
)
par(opar)

The values for Egypt (EG) and Sudan (SU) stand out from the rest. The histogram shows a strongly right-skewed pattern with many small values and a few large ones.

In addition to headtails and fisher, we use quantile to compare the classification styles. Quantile breaks are based on ranks rather than the size of the gaps between observations.

Applying the three aforementioned methods to break the data:

brks_ht <- classIntervals(afcon$totcon, style = "headtails")
print(brks_ht)
# Same number of classes for "fisher"
nclass <- length(brks_ht$brks) - 1
brks_fisher <- classIntervals(afcon$totcon,
  style = "fisher",
  n = nclass
)
print(brks_fisher)

brks_quantile <- classIntervals(afcon$totcon,
  style = "quantile",
  n = nclass
)
print(brks_quantile)

pal1 <- c("wheat1", "wheat2", "red3")
opar <- par(no.readonly = TRUE)
par(mar = c(2, 2, 2, 1), cex = 0.8)
plot(brks_ht, pal = pal1, main = "headtails")
plot(brks_fisher, pal = pal1, main = "fisher")
plot(brks_quantile, pal = pal1, main = "quantile")
par(opar)

The top three classes of headtails contain five observations, whereas those of fisher contain 13. In this example, headtails gives more detail at the upper end of the distribution.

The next plot compares density estimates of the rescaled totcon values and the class assignments from each method. The original values are rescaled to [1, nclass] and shifted by -0.5 for the visual comparison.

# Helper function to reescale values
help_reescale <- function(x, min = 1, max = 10) {
  r <- (x - min(x)) / (max(x) - min(x))
  r <- r * (max - min) + min
  return(r)
}
afcon$ecdf_class <- help_reescale(afcon$totcon,
  min = 1 - 0.5,
  max = nclass - 0.5
)
afcon$ht_breaks <- cut(afcon$totcon,
  brks_ht$brks,
  labels = FALSE,
  include.lowest = TRUE
)

afcon$fisher_breaks <- cut(afcon$totcon,
  brks_fisher$brks,
  labels = FALSE,
  include.lowest = TRUE
)

afcon$quantile_break <- cut(afcon$totcon,
  brks_quantile$brks,
  labels = FALSE,
  include.lowest = TRUE
)

opar <- par(no.readonly = TRUE)
par(mar = c(4, 4, 1, 1), cex = 0.8)
plot(
  density(afcon$ecdf_class),
  ylim = c(0, 0.8),
  lwd = 2,
  main = "",
  xlab = "class"
)
lines(density(afcon$ht_breaks), col = "darkblue", lty = 2)
lines(density(afcon$fisher_breaks), col = "limegreen", lty = 2)
lines(density(afcon$quantile_break),
  col = "red3",
  lty = 2
)
legend("topright",
  legend = c(
    "Continuous", "headtails",
    "fisher", "quantile"
  ),
  col = c("black", "darkblue", "limegreen", "red3"),
  lwd = c(2, 1, 1, 1),
  lty = c(1, 2, 2, 2),
  cex = 0.8
)
par(opar)

In this example, the class assignments from headtails retain more of the original imbalance. By contrast, quantile aims for similar numbers of observations per class, regardless of the gaps between values.

We now compare the methods using proportional-symbol maps. Symbol size represents totcon, while color distinguishes the classes.

The first map uses proportional symbols without class-based colors, providing a reference for the three classified maps.

library(sf)
library(giscoR)
library(cartography)

opar <- par(no.readonly = TRUE)

par(
  mfrow = c(2, 2),
  mar = c(1, 1, 1, 1),
  bg = "white"
)
africa <- gisco_get_countries(resolution = 60, region = "Africa", epsg = 3857)

afcon.sf <- st_as_sf(afcon, crs = 4326, coords = c("x", "y"))
afcon.sf <- st_transform(afcon.sf, st_crs(africa))
# afcon.sf <- st_join(africa[, "admin"], afcon.sf)
afcon.sf <- afcon.sf[order(afcon.sf$totcon), ]



# High granularity map
plot(st_geometry(africa), col = "grey80", border = NA)
propSymbolsLayer(
  afcon.sf,
  var = "totcon",
  inches = 0.2,
  col = adjustcolor("grey10", alpha.f = 0.5),
  border = NA
)
title(main = "High granularity map")

# Quantile

pal <- hcl.colors(5, palette = "inferno", alpha = 0.6)
plot(st_geometry(africa), col = "grey80", border = NA)
propSymbolsTypoLayer(
  afcon.sf,
  var = "totcon",
  inches = 0.2,
  col = pal,
  border = NA,
  legend.var.pos = "n",
  legend.var2.pos = "bottomleft",
  var2 = "quantile_break"
)
title(main = "Quantile")


# Fisher
plot(st_geometry(africa), col = "grey80", border = NA)
propSymbolsTypoLayer(
  afcon.sf,
  var = "totcon",
  inches = 0.2,
  col = pal,
  border = NA,
  legend.var.pos = "n",
  legend.var2.pos = "bottomleft",
  var2 = "fisher_breaks"
)
title(main = "Fisher")

# Head Tails
plot(st_geometry(africa), col = "grey80", border = NA)
propSymbolsTypoLayer(
  afcon.sf,
  var = "totcon",
  inches = 0.2,
  col = pal,
  border = NA,
  legend.var.pos = "n",
  legend.var2.pos = "bottomleft",
  var2 = "ht_breaks"
)
title(main = "Head Tails")
par(opar)

Compared with the unclassified proportional-symbol map, headtails makes the most extreme values easier to distinguish. The quantile style groups observations by rank, while fisher provides an intermediate view in this example.

It is also important to note that headtails and fisher reveal different information that can be useful depending on the context. While headtails highlights the outliers, it fails to provide good clustering on the tail, while fisher seems to reflect these patterns better. This can be observed in the values of Western Africa and the Niger River Basin, where headtails does not highlight any special cluster of conflicts, while fisher suggests a potential cluster consistent with the unclassified map. This can be confirmed in the histogram generated previously, where a concentration of totcon around 1,000 is visible.

References

Jiang, Bin. 2013. "Head/Tail Breaks: A New Classification Scheme for Data with a Heavy-Tailed Distribution." The Professional Geographer 65 (3): 482–94. DOI.

———. 2019. "A Recursive Definition of Goodness of Space for Bridging the Concepts of Space and Place for Sustainability." Sustainability 11 (15): 4091. DOI.

Jiang, Bin, Xintao Liu, and Tao Jia. 2013. "Scaling of Geographic Space as a Universal Rule for Map Generalization." Annals of the Association of American Geographers 103 (4): 844–55. DOI.

Jiang, Bin, and Junjun Yin. 2013. "Ht-Index for Quantifying the Fractal or Scaling Structure of Geographic Features." Annals of the Association of American Geographers 104 (3): 530–40. DOI.

Taleb, Nassim Nicholas. 2008. The Black Swan: The Impact of the Highly Improbable. 1st ed. London: Random House.

Vasicek, Oldrich. 2002. "Loan Portfolio Value." Risk, December, 160–62.