Z-score

Z-score, also called a standard score, of an observation is [broadly speaking] a distance from the population center measured in number of normalization units. The default choice for center is sample mean and for normalization unit is standard deviation.

⬛ Observation is not an outlier based on z-score if its absolute value of default z-score is lower then some threshold (popular choice is 3).

Here is the function for identifying non-outliers based on z-score:

isnt_out_z <- function(x, thres = 3, na.rm = TRUE) {
  abs(x - mean(x, na.rm = na.rm)) <= thres * sd(x, na.rm = na.rm)
}

It takes a numeric vector as input and returns logical vector of the same length indicating whether input value is a non-outlier.

Z-score with MAD

Median Absolute Deviation is a robust normalization unit based on median as a population center. In order to use MAD “as a consistent estimator for the estimation of the standard deviation” one takes its value multiplied by a factor. This way base R function mad is implemented.

⬛ Observation is not an outlier based on MAD if its absolute value of z-score with median as center and MAD as normalization unit is lower then some threshold (popular choice is 3).

isnt_out_mad <- function(x, thres = 3, na.rm = TRUE) {
  abs(x - median(x, na.rm = na.rm)) <= thres * mad(x, na.rm = na.rm)
}

Tukey’s fences

Tukey’s fences is a technique used in box plots. The non-outlier range is defined with \([Q_1 - k(Q_3 - Q_1),~ Q_3 + k(Q_3 - Q_1)]\), where \(Q_1\) and \(Q_3\) are the lower and upper quartiles respectively, \(k\) - some nonnegative constant (popular choice is 1.5).

⬛ Observation is not an outlier based on Tukey’s fences if its value lies in non-outlier range.

isnt_out_tukey <- function(x, k = 1.5, na.rm = TRUE) {
  quar <- quantile(x, probs = c(0.25, 0.75), na.rm = na.rm)
  iqr <- diff(quar)
  
  (quar[1] - k * iqr <= x) & (x <= quar[2] + k * iqr)
}

Mahalanobis distance

All previous approaches were created for univariate numerical data. To detect outliers in multivariate case one can use Mahalanobis distance to reduce to univariate case and then apply known techniques.

⬛ Observation is not an outlier based on Mahalanobis distance if its distance is not an outlier.

maha_dist <- . %>% select_if(is.numeric) %>%
    mahalanobis(center = colMeans(.), cov = cov(.))

isnt_out_maha <- function(tbl, isnt_out_f, ...) {
  tbl %>% maha_dist() %>% isnt_out_f(...)
}

This function takes as input a data frame of interest (with possible non-numeric columns which are ignored) and function performing univariate outlier detection. It returns a logical vector of the same length as number of rows in input data frame.

To read more about practical usefulness of Mahalanobis distance in detecting outliers go to Steffen’s very helpful post.

Definition of non-outlier row

Package ruler, based on dplyr grammar of data manipulation, offers tools for validating the following data units: data as a whole, group [of rows] as a whole, column as a whole, row as a whole, cell. Our primary interest is row as a whole. However, using this framework, we can construct several approaches for definition of the non-outlier row:

1. Row is not an outlier based on some column if it doesn’t contain outlier (computed based on the target column) on the intersection with that column. In other words, first a univariate outlier detection is performed based solely on data from target column and then all rows containing non-outliers are named non-outlier rows.
2. Row is not an outlier based on Mahalanobis distance if its distance (computed based on the selected numeric columns) is not an outlier.
3. Row is not an outlier based on grouping if it is a part of a non-outlier group [of rows]. A group [of rows] is not an outlier if its summary value is not an outlier among summary values of other groups.

Note that all listed approached depend on the choice of the univariate outlier detection method. We will use all three previously listed univariate techniques.

isnt_out_funs <- funs(
  z = isnt_out_z,
  mad = isnt_out_mad,
  tukey = isnt_out_tukey
)

Implementation

In ruler framework rules are defined in packs (to learn more go to ruler README and vignettes).

diamonds %>% transmute_if(is.numeric, isnt_out_funs)
## # A tibble: 53,940 x 21
##   carat_z depth_z table_z price_z   x_z   y_z   z_z carat_mad depth_mad
##     <lgl>   <lgl>   <lgl>   <lgl> <lgl> <lgl> <lgl>     <lgl>     <lgl>
## 1    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE      TRUE
## 2    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE      TRUE
## 3    TRUE   FALSE   FALSE    TRUE  TRUE  TRUE  TRUE      TRUE     FALSE
## 4    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE      TRUE
## 5    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE      TRUE
## # ... with 5.394e+04 more rows, and 12 more variables: table_mad <lgl>,
## #   price_mad <lgl>, x_mad <lgl>, y_mad <lgl>, z_mad <lgl>,
## #   carat_tukey <lgl>, depth_tukey <lgl>, table_tukey <lgl>,
## #   price_tukey <lgl>, x_tukey <lgl>, y_tukey <lgl>, z_tukey <lgl>

diamonds %>%
  transmute(maha = maha_dist(.)) %>%
  transmute_at(vars(maha = maha), isnt_out_funs)
## # A tibble: 53,940 x 3
##   maha_z maha_mad maha_tukey
##    <lgl>    <lgl>      <lgl>
## 1   TRUE     TRUE       TRUE
## 2   TRUE    FALSE      FALSE
## 3   TRUE    FALSE      FALSE
## 4   TRUE     TRUE       TRUE
## 5   TRUE     TRUE       TRUE
## # ... with 5.394e+04 more rows

data_tbl <- diamonds %>%
  unite(col = "group", cut, color, clarity)

compute_group_non_outliers <- . %>%
  # Compute per group mean values of columns
  group_by(group) %>%
  summarise_if(is.numeric, mean) %>%
  ungroup() %>%
  # Detect outliers among groups
  mutate_if(is.numeric, isnt_out_funs) %>%
  # Remove unnecessary columns
  select_if(Negate(is.numeric))

data_tbl %>% compute_group_non_outliers()
## # A tibble: 276 x 22
##        group carat_z depth_z table_z price_z   x_z   y_z   z_z carat_mad
##        <chr>   <lgl>   <lgl>   <lgl>   <lgl> <lgl> <lgl> <lgl>     <lgl>
## 1  Fair_D_I1   FALSE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE     FALSE
## 2  Fair_D_IF    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE
## 3 Fair_D_SI1    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE
## 4 Fair_D_SI2    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE
## 5 Fair_D_VS1    TRUE    TRUE    TRUE    TRUE  TRUE  TRUE  TRUE      TRUE
## # ... with 271 more rows, and 13 more variables: depth_mad <lgl>,
## #   table_mad <lgl>, price_mad <lgl>, x_mad <lgl>, y_mad <lgl>,
## #   z_mad <lgl>, carat_tukey <lgl>, depth_tukey <lgl>, table_tukey <lgl>,
## #   price_tukey <lgl>, x_tukey <lgl>, y_tukey <lgl>, z_tukey <lgl>

Exposure

Column and Mahalanobis based definition of non-outlier rows can be expressed with row packs and group based - as group packs.

row_packs_isnt_out <- row_packs(
  # Non-outliers based on some column
  column = . %>% transmute_if(is.numeric, isnt_out_funs),
  # Non-outliers based on Mahalanobis distance
  maha = . %>% transmute(maha = maha_dist(.)) %>%
    transmute_at(vars(maha = maha), isnt_out_funs)
)

group_packs_isnt_out <- group_packs(
  # Non-outliers based on grouping
  group = compute_group_non_outliers,
  .group_vars = "group"
)

Application of all those packs is called exposing process. The result is an exposure from which we can extract tidy data validation report using get_report.

# Don't remove obeyers to compute total number of applied rules
full_report <- data_tbl %>%
  expose(row_packs_isnt_out, group_packs_isnt_out,
         .remove_obeyers = FALSE) %>%
  get_report()

used_rules <- full_report %>%
  distinct(pack, rule)

breaker_report <- full_report %>%
  filter(!(value %in% TRUE))

used_rules contains data about all definitions of non-outlier rows applied to data. They are encoded with combination of columns pack and rule.

breaker_report contains data about data units that break certain rules. Packs column and maha has actual row numbers of data_tbl listed in id column of report (for rows which should be considered as outliers).

On the other hand, pack group defines group pack and is represented in breaker_report with id 0. To obtain row outliers based on grouping we need to expand those rows with information about rows in the data that belong to those groups. This can be done using dplyr::left_join():

group_breakers <- breaker_report %>%
  # Filter group packs
  filter(pack == "group") %>%
  # Expand rows by matching group with its rows
  select(-id) %>%
  left_join(
    y = data_tbl %>% transmute(var = group, id = 1:n()),
    by = "var"
  ) %>%
  select(pack, rule, var, id, value)

outliers <- bind_rows(
  breaker_report %>% filter(pack != "group"),
  group_breakers
) %>%
  select(pack, rule, id)

# Not all group based definitions resulted with outliers
outliers %>%
  count(pack, rule) %>%
  filter(pack == "group") %>%
  print(n = Inf)
## # A tibble: 13 x 3
##     pack        rule     n
##    <chr>       <chr> <int>
##  1 group   carat_mad    37
##  2 group carat_tukey    37
##  3 group     carat_z    29
##  4 group   depth_mad  1093
##  5 group depth_tukey  1016
##  6 group     depth_z   156
##  7 group   price_mad   209
##  8 group price_tukey  1146
##  9 group     price_z    44
## 10 group   table_mad   920
## 11 group table_tukey     8
## 12 group     table_z     7
## 13 group         z_z    23

Tibble outliers contains data about outlier rows. Combination of columns pack and rule defines non-outlier/outlier definition approach and column id defines row number of input data frame that should be considered an outlier based on the definition.

Definitions with most outliers are as follows:

outliers %>%
  count(pack, rule, sort = TRUE)
## # A tibble: 37 x 3
##     pack        rule     n
##    <chr>       <chr> <int>
## 1   maha    maha_mad  6329
## 2   maha  maha_tukey  5511
## 3 column   price_mad  5386
## 4 column price_tukey  3540
## 5 column   table_mad  2560
## # ... with 32 more rows

Two out of three Mahalanobis based definition yielded the most row outliers.

Combination

Given outliers data frame, one can do whatever he/she wants to identify outliers. Here we will use the basic combination approach based on average score.

Combined outlier detection score for certain row can be defined as share of applied methods that tagged it as outlier. Alternatively one can define it just as number of those methods as it will only change absolute value of the result and not the order.

outlier_score <- outliers %>%
  group_by(id) %>%
  # nrow(used_rules) equals total number of applied methods
  summarise(score = n() / nrow(used_rules))

# Top 10 outliers
outlier_score %>% arrange(desc(score)) %>% slice(1:10)
## # A tibble: 10 x 2
##       id     score
##    <int>     <dbl>
##  1 26432 0.5777778
##  2 27416 0.5777778
##  3 27631 0.5777778
##  4 27131 0.4666667
##  5 23645 0.4222222
##  6 26445 0.4222222
##  7 26745 0.4000000
##  8 27430 0.4000000
##  9 15952 0.3777778
## 10 17197 0.3777778

Finally we will tag those rows as strong outliers which has score more than 0.2 (subjective threshold which should be researched more).

diam_tbl <- diamonds %>%
  mutate(id = 1:n()) %>%
  left_join(y = outlier_score, by = "id") %>%
  mutate(
    score = coalesce(score, 0),
    is_out = if_else(score > 0.2, "Outlier", "Not outlier")
  )

# Total number of outliers
sum(diam_tbl$score > 0.2)
## [1] 161

Tibble diam_tbl is basically the diamonds but with three more columns: id for row number, score for combined outlier score and is_out for non-outlier/outlier tag.

Plots illustrating strong outliers:

theme_set(theme_bw())

plot_outliers <- function(tbl, x, y, facet_var) {
  tbl %>%
    arrange(is_out) %>%
    ggplot(aes_string(x, y, colour = "is_out")) +
      geom_point() +
      facet_wrap(facets = facet_var) +
      scale_colour_manual(values = c("#AAAAAA", "#004080")) +
      guides(colour = guide_legend(title = NULL,
                                   override.aes = list(size = 4))) +
      labs(title = paste0("Strong outliers illustration by ", facet_var)) +
      theme(legend.position = "bottom",
            legend.text = element_text(size = 14))
}

diam_tbl %>% plot_outliers("carat", "price", facet_var = "cut")

diam_tbl %>% plot_outliers("x", "depth", facet_var = "color")

diam_tbl %>% plot_outliers("price", "table", facet_var = "clarity")

Based on those plots we see the complicated nature of “strong outliers”. They are not necessary located “on the edge” of two-dimensional scatter plots, but most extreme cases are tagged as outliers.

Also one interesting observation: most outliers are concentrated in the combination of “Fair” cut, “J” colour and “I1” clarity which are worst options among their features. The reason of this effect is group-based definitions of non-outliers which tagged certain groups more than others:

breaker_report %>%
  filter(pack == "group") %>%
  count(var, sort = TRUE) %>%
  print(n = 10)
## # A tibble: 47 x 2
##            var     n
##          <chr> <int>
##  1   Fair_D_I1     7
##  2   Fair_J_I1     7
##  3 Fair_H_VVS1     6
##  4  Ideal_J_I1     6
##  5 Fair_J_VVS1     5
##  6 Fair_G_VVS1     4
##  7 Fair_D_VVS1     3
##  8   Fair_E_I1     3
##  9   Fair_F_I1     3
## 10   Fair_H_I1     3
## # ... with 37 more rows

Here we see that “Fair” cut is among majority of top breaker groups. There are also some interesting combinations: Fair_D_I1 (“worst”-“best”-“worst”), Fair_J_I1 (“worst”-“worst”-“worst”), Ideal_J_I1 (“best”-“worst”-“worst”).

This fact might be interpreted as suggested combined outlier detection approach discovered notion of diamond quality without prior knowledge about it.