Abstract

This vignette discusses the implementation of the “Head/tail breaks” style (Jiang (2013)) on the classIntervals function. A step-by-step example is presented in order to clarify the method. A case study using spData::afcon is also included, as well as a test suite checking the performance and validation of the implementation.

Introduction

The Head/tail breaks, sometimes referred as ht-index (Jiang and Yin (2013)), is a classification scheme introduced by Jiang (2013) in order to find groupings or hierarchy for data with a heavy-tailed distribution.

Heavy-tailed distributions are heavily right skewed, with a minority of large values in the head and a majority of small values in the tail. This imbalance between the head and tail, or between many small values and a few large values, can be expressed as “far more small things than large things”.

Heavy tailed distributions are commonly characterized by a power law, a lognormal or an exponential function. Nature, society, finance (Vasicek (2002)) and our daily lives are full of rare and extreme events, which are termed “black swan events” (Taleb (2008)). This line of thinking provides a good reason to reverse our thinking by focusing on low-frequency events.

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 on a four-step process performed recursively until a stopping condition is satisfied. Given a vector of values var the process can be described as follows:

1. Compute mu = mean(var).
2. Break var into the tail (as var < mu) and the head (as var > mu).
3. Assess if the proportion of head over var is lower or equal than a given threshold (i.e. length(head)/length(var) <= thr)
4. If 3 is TRUE, repeat 1 to 3 until the condition is FALSE or no more partitions are possible (i.e. head has less than two elements expressed as length(head) < 2).

It is important to note that, at the beginning of a new iteration, var is replaced by head. The underlying hypothesis is to create partitions until the head and the tail are balanced in terms of distribution.So the stopping criteria is satisfied when the last head and the last tail are evenly balanced.

In terms of threshold, Jiang, Liu, and Jia (2013) set 40% as a good approximation, meaning that if the head contains more than 40% of the observations the distribution is not considered 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

A step-by-step example in R (for illustrative purposes) has been developed:

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)

As it can be seen, in each iteration the resulting head gradually loses the high-tail property, until the stopping condition is met.

The resulting breaks are then defined as breaks = c(min(var), mu(iter=1), ..., mu(iter), max(var)).

Implementation on classInt package

The implementation in the classIntervals function should replicate the results:

ht_sample_par <- classIntervals(sample_par, style = "headtails")
brks == ht_sample_par$brks
#> [1] TRUE TRUE TRUE TRUE TRUE TRUE
print(ht_sample_par)
#> style: headtails
#> [1.000295,5.675463) [5.675463,27.23693) [27.23693,85.17664) [85.17664,264.7126) 
#>                 855                 114                  25                   3 
#> [264.7126,523.6254] 
#>                   3

As stated in Jiang (2013), the number of breaks is naturally determined, however the thr parameter could help to adjust the final number. A lower value on thr would provide less breaks while a larger thr would 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(classIntervals(sample_par, style = "headtails", thr = 0))
#> style: headtails
#> [1.000295,5.675463) [5.675463,523.6254] 
#>                 855                 145
plot(
  classIntervals(sample_par, style = "headtails", thr = 0),
  pal = pal1,
  main = "thr = 0"
)

# Two breaks
print(classIntervals(sample_par, style = "headtails", thr = 0.2))
#> style: headtails
#> [1.000295,5.675463) [5.675463,27.23693) [27.23693,523.6254] 
#>                 855                 114                  31
plot(
  classIntervals(sample_par, style = "headtails", thr = 0.2),
  pal = pal1,
  main = "thr = 0.2"
)

# Default breaks: 0.4
print(classIntervals(sample_par, style = "headtails"))
#> style: headtails
#> [1.000295,5.675463) [5.675463,27.23693) [27.23693,85.17664) [85.17664,264.7126) 
#>                 855                 114                  25                   3 
#> [264.7126,523.6254] 
#>                   3
plot(classIntervals(sample_par, style = "headtails"),