Benford law as a feature engineering technique for Machine Learning Models

R Machine Learning Series

Benford law

Nếu chúng ta gặp câu hỏi “Tần suất xuất hiện các chữ số từ 1 đến 9 tại vị trí chữ số đầu tiên của bộ dữ liệu về thu nhập cá nhân của 9 triệu người là bao nhiêu?” thì chúng ta dễ bị lôi kéo bởi ý tưởng rằng tần suất đó tuân theo Uniform Distribution và do vậy tần suất xuất hiện của mỗi chữ số sẽ là tương tự nhau và xấp xỉ 1/9 = 11.11% như chúng ta có thể thấy:

# Clear workspace: 

rm(list = ls())

# Create our theme: 

library(tidyverse)
library(extrafont)
my_font <- "Roboto Condensed" 

my_theme <- function(...) {
  
  theme(panel.grid.major.x = element_blank()) + 
    theme(panel.grid.minor.y = element_blank()) + 
    theme(axis.ticks.x = element_blank()) + 
    theme(text = element_text(family = my_font, size = 16)) + 
    theme(plot.caption = element_text(size = 10)) + 
    theme(plot.title = element_text(vjust = 4)) + 
    theme(plot.margin = unit(rep(0.8, 4), "cm"))
}

random_df <- tibble(digit = sample(1:9, size = 9000000, replace = TRUE))

random_df %>% 
  group_by(digit) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(digit = as.factor(digit)) %>% 
  ggplot(aes(digit, n / sum(n))) + 
  geom_col(width = 0.65, fill = "steelblue") + 
  scale_y_continuous(breaks = seq(0, 0.12, 0.025), expand = c(0, 0), limits = c(0, 0.12), label = scales::percent) + 
  my_theme() + 
  labs(x = NULL, y = NULL, title = "Figure 1: Frequency of occurrences", 
       caption = "Source: Data based on Random Generator")

Giả định rằng các chữ số từ 1 đến 9 xuất hiện theo Uniform Distribution là một giả định hợp lí vì rõ ràng là: không có lí do gì mà một chữ số cụ thể nào đó, như số 1 chẳng hạn, lại có cơ hội xuất hiện nhiều hơn số 9 hay một số bất kì nào khác. Nếu khác đi, có vẻ thượng đế sẽ ưu ái với số 1.

Tuy nhiên năm 1938 một nhà vật lí kiêm kĩ sư điện tại GE phủ nhận giả thuyết này và cho rằng các chứ số từ 1 đến 9 (kí hiệu là d) sẽ có xác suất xuất hiện \(P(d)\) được thể hiện theo công thức sau:

\[P(d) = log_{10}(1 + \frac{1}{d})\]

Công thức này gọi là Benford law thực ra đã được đề cập trước đó vào năm 1881 bởi Simon Newcomb và do vậy còn có tên gọi khác là Newcomb–Benford law. Chúng ta minh họa định luật này bằng công cụ hình ảnh:

# Data Frame based on Benford's law (from https://en.wikipedia.org/wiki/Benford%27s_law): 

df_benford <- tibble(digit = 1:9, prob = sapply(1:9, function(d) {log10(1 + 1 / d)}))

df_benford %>% 
  mutate(digit = as.character(digit)) %>% 
  mutate(text = round(100*prob, 1)) %>% 
  mutate(text = paste0(text, "%")) -> df_plot

# Graph: 

df_plot %>% 
  ggplot(aes(digit, prob)) + 
  geom_col(width = 0.65, fill = "steelblue") + 
  geom_line(aes(as.numeric(digit), prob), color = "orange", size = 1) + 
  geom_point(aes(as.numeric(digit), prob), color = "firebrick", size = 2) + 
  geom_text(aes(label = text), vjust = 1.3, color = "white", family = my_font, size = 4.5) + 
  scale_y_continuous(breaks = seq(0, 0.305, 0.05), expand = c(0, 0), limits = c(0, 0.31), label = scales::percent) + 
  my_theme() + 
  labs(x = NULL, y = NULL, title = "Figure 2: The Law of Anomalous Numbers", 
       caption = "Source: Data based on Benford's law")

Theo định luật này thì rõ ràng số 1 có tần suất xuất hiện lớn nhất và các số càng lớn thì tần suất xuất hiện càng ít đi. Chi tiết hơn về định luật này cũng như các ứng dụng của nó bạn đọc có thể tìm hiểu thêm ở đây. Chúng ta sẽ kiểm tra xem liệu các chữ số từ 1 đến 9 fit ở mức độ như thế nào với định luật Benford với một số bộ dữ liệu cụ thể.

Trước hết chúng ta viết một hàm nhận input là chuỗi các số từ 1 đến 9 - là chữ số ở vị trí thứ nhất từ một chuỗi số bất kì và hiển thị tần suất xuất hiện của chúng đồng thời so sánh tần suất thực tế với tần suất được dự báo bởi công thức Benford:

compare_with_benford <- function(x) {
  
  tibble(digit = x) %>% 
    na.omit() %>% 
    group_by(digit) %>% 
    count() %>% 
    ungroup() %>% 
    mutate(digit = as.factor(digit)) %>% 
    ggplot(aes(digit, n / sum(n))) + 
    geom_col(width = 0.65, fill = "steelblue") + 
    geom_line(data = df_benford, aes(as.numeric(digit), prob), color = "orange", size = 1) + 
    geom_point(data = df_benford, aes(as.numeric(digit), prob), color = "firebrick", size = 2) + 
    my_theme() + 
    theme(axis.title = element_blank())

}

Với chuỗi poptotal từ bộ số liệu midwest đi kèm với thư viện ggplot2:

midwest %>% 
  pull(poptotal) %>% 
  as.character() %>% 
  str_sub(start = 1, end = 1) -> totalpop

compare_with_benford(x = totalpop) + 
    labs(title = "Figure 3: Frequency of occurrences, poptotal (437 obs)", 
       caption = "Source: Demographic information of midwest counties")

Rõ ràng phân bố của các chữ số chưa fit lắm với định luật Benford. Điều này có thể nằm ở nguyên nhân là mới chỉ có 437 quan sát - một số lượng chưa đủ nhiều. Nếu quy mô dữ liệu tăng lên thì mức độ fit với định luật Benford có thể cao hơn:

midwest %>% 
  select(contains("pop")) %>% 
  select(3:7) %>% 
  gather(a, b) %>% 
  filter(b != 0) %>% 
  pull(b) %>% 
  as.character() %>% 
  str_sub(start = 1, end = 1) -> popden_all_types

compare_with_benford(x = popden_all_types) + 
    labs(title = "Figure 4: Frequency of occurrences, poptotal for all types (2180 obs)", 
       caption = "Source: Demographic information of midwest counties")

Vậy là khi số quan sát tăng lên 2180 thì tần suất xuất hiện của các số từ 1 đến 9 tiệm cận hơn với các tần suất dự báo bởi quy luật Benford. Chúng ta lại xét một bộ dữ liệu khác là txhousing trong đó volume là tổng giá trị doanh thu các thương vụ giao dịch bất động sản ở Texas:

txhousing %>% 
  pull(volume) %>% 
  as.character() %>% 
  str_sub(start = 1, end = 1) -> volume

compare_with_benford(x = volume) + 
    labs(title = "Figure 5: Frequency of occurrences, volume (8602 obs)", 
       caption = "Source: Information about the housing market in Texas provided by the TAMU real estate center.")

Các bằng chứng thực nghiệm ở trên về tần suất suất xuất hiện các chữ số từ 1 đến 9 trong thế giới thực có vẻ xác nhận định luật Benford. Và bởi vậy định luật này có nhiều ứng dụng. Một trong những ứng dụng đó là xác định các báo cáo tài chính gian lận của công ti Enron. Bạn đọc có thể tham khảo bài viết của The Wall Street Journal về áp dụng Benford law với tình huống của Enron ở đây.

Benford law as a feature engineering technique

Định luật Benford có thể được sử dụng như là một kĩ thuật feature engineering cho các mô hình Machine Learning. Dưới đây chúng ta sẽ so sánh chất lượng dự báo/phân loại của mô hình phân loại Logistic Classifier cho hai tình huống:

1. Tạo ra các features mới bằng cách sử dụng Benford law
2. Sử dụng chỉ các features nguyên bản.

Bộ dữ liệu sử dụng là creditcard.csv từ cuộc thi Credit Card Fraud Detection và có thể download tại đây. Ở đây chúng ta sẽ tạo ra các features mới từ Amount - một feature sẵn có của bộ dữ liệu. Trước hết ta so sánh tần suất xuât hiện các chữ số từ 1 đến 9:

# Load data: 
data <- read_csv("creditcard.csv")

# Data manipulation: 

data %>% 
  mutate(Class = case_when(Class == 1 ~ "Fraud", TRUE ~ "NonFraud")) %>% 
  mutate(Amount_text = as.character(Amount)) %>% 
  mutate(digit_0 = case_when(!str_detect(Amount_text, "\\.") ~ "Yes", TRUE ~ "No")) %>% 
  mutate(digit = str_sub(Amount_text, start = 1, end = 1)) -> data

# Calculate frequencies and join datasets: 

data %>% 
  filter(digit != "0") %>% 
  group_by(digit) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(actual_rate = n / sum(n)) %>% 
  full_join(df_plot %>% select(digit, prob)) %>% 
  mutate(diff_freq = actual_rate - prob) -> digit_from_creditdata

# Compare actual Frequency and Benford's law: 

x <- data %>% 
  filter(digit != "0") %>% 
  pull(digit)

compare_with_benford(x) + 
  scale_y_continuous(breaks = seq(0, 0.42, 0.05), expand = c(0, 0), limits = c(0, 0.42), label = scales::percent) + 
  labs(x = NULL, y = NULL, title = "Figure 6: Actual frequency compared with Benford's law", 
       caption = "Source: Credit Card Fraud Detection, https://www.kaggle.com")

Từ Figure 6 rõ ràng tần suất xuất hiện của số 1 là tương đối cao (giống như tình huống về các số liệu tài chính gian lận của công ti Enron). Chúng ta thử kiểm tra sự khác biệt về tỉ lệ các giao dịch gian lận giữa nhóm có và không xuất hiện chữ số 1:

data %>% 
  mutate(digit_1 = case_when(str_detect(digit, "1") ~ "Yes", TRUE ~ "No")) %>% 
  group_by(Class, digit_1) %>% 
  count() %>% 
  ungroup() %>% 
  spread(value = "n", key = "Class") %>% 
  mutate(rate = Fraud / (Fraud + NonFraud)) %>% 
  mutate(diff = rate[2] / rate[1])

## # A tibble: 2 x 5
##   digit_1 Fraud NonFraud    rate  diff
##   <chr>   <int>    <int>   <dbl> <dbl>
## 1 No        290   184062 0.00157  1.28
## 2 Yes       202   100253 0.00201  1.28

Rõ ràng nhóm giao dịch có xuất hiện chữ số 1 có tỉ lệ giao dịch gian lận cao hơn 28% so với nhóm không có chữ số 1. Mặc dù Benford law không cover cho chữ số 0 nhưng chúng ta vẫn nên điều tra xem giữa nhóm có chữ số 0 và phần còn lại có khác biệt đáng kể hay không về tỉ lệ các giao dịch tài chính gian lận:

data %>% 
  mutate(digit_0 = case_when(str_detect(digit, "0") ~ "Yes", TRUE ~ "No")) %>% 
  group_by(Class, digit_0) %>% 
  count() %>% 
  ungroup() %>% 
  spread(value = "n", key = "Class") %>% 
  mutate(rate = Fraud / (Fraud + NonFraud)) %>% 
  mutate(diff = rate[2] / rate[1])

## # A tibble: 2 x 5
##   digit_0 Fraud NonFraud    rate  diff
##   <chr>   <int>    <int>   <dbl> <dbl>
## 1 No        424   267579 0.00158  2.56
## 2 Yes        68    16736 0.00405  2.56

Kết quả này cho thấy các giao dịch có chữ số 0 có tỉ lệ giao dịch gian lận cao hơn 156% so với nhóm còn lại. Chúng là các giao dịch quy mô nhỏ có giá trị giao dịch chưa đến 1 hoặc bằng 0 như chúng ta có thể thấy:

data %>% 
  mutate(digit_0 = case_when(str_detect(digit, "0") ~ "Yes", TRUE ~ "No")) %>% 
  filter(str_detect(digit, "0")) %>% 
  select(Amount, digit_0) %>% 
  sample_n(10)

## # A tibble: 10 x 2
##    Amount digit_0
##     <dbl> <chr>  
##  1   0.99 Yes    
##  2   0.77 Yes    
##  3   0.89 Yes    
##  4   0    Yes    
##  5   0.01 Yes    
##  6   0.89 Yes    
##  7   0.89 Yes    
##  8   0.76 Yes    
##  9   0.89 Yes    
## 10   0.89 Yes

Từ các phân tích trên chúng ta mạnh dạn đưa ra một giả định rằng feature mới cho biết thông tin về giá trị của chữ số xuất hiện ở the first digit có thể là một chỉ báo tốt, hoặc chí ít là có ý nghĩa trong việc phân biệt giao dịch gian lận. Do vậy chúng ta tạo ra hai bộ số liệu như sau:

1. Bộ dữ liệu nguyên bản.
2. Bộ dữ liệu có các features mới dựa trên Benford law.

Dưới đây là R codes:

# Transform data: 

data %>% 
  mutate(digit_0 = case_when(str_detect(digit, "0") ~ "Yes", TRUE ~ "No")) %>% 
  mutate(digit_0 = as.factor(digit_0), Class = as.factor(Class)) %>%
  select(-Time, - Amount_text) %>% 
  full_join(digit_from_creditdata %>% select(digit, diff_freq)) %>% 
  mutate(diff_freq = replace_na(diff_freq, 0), digit = as.factor(digit)) -> df_modelling

# Split data:

library(caret)
set.seed(1)
id <- createDataPartition(y = df_modelling$Class, p = 0.2, list = FALSE)

df_modelling %>% slice(id) -> train_benford
df_modelling %>% slice(-id) -> test_benford

train_benford %>% select(-digit_0, -digit, -diff_freq) -> train_none
test_benford %>% select(-digit_0, -digit, -diff_freq) -> test_none

Huấn luyện Logistic Classifier với hai bộ dữ liệu đã được chuẩn bị và sử dụng các model có được thực hiện dự báo xác suất cho sự kiện “giao dịch gian lận” trên bộ dữ liệu test:

# Set conditions for training model and cross-validation: 

set.seed(1)
number <- 3
repeats <- 1
control <- trainControl(method = "repeatedcv", 
                        number = number, 
                        repeats = repeats, 
                        classProbs = TRUE, 
                        summaryFunction = twoClassSummary, 
                        allowParallel = TRUE)

# Parallel Computing with multicores from our computer: 

library(parallel)
nCores <- detectCores(logical = TRUE)
nThreads <- detectCores(logical = TRUE)
cl <- makeCluster(nThreads)


# Logistic Classifier with feature engineering using Benford law: 

model <- "glm"

set.seed(1)
glm1 <- train(Class ~., 
              data = train_benford %>% select(-Amount), 
              trControl = control,
              method = model,
              metric = "ROC")

# Logistic Classifier with original data: 

set.seed(1)
glm2 <- train(Class ~., 
              data = train_none, 
              trControl = control,
              method = model, 
              metric = "ROC")


# Use models for predicting probability of fraud: 

predict(glm1, test_benford %>% select(-Amount), type = "prob") %>% pull(Fraud) -> pd1
predict(glm2, test_none, type = "prob") %>% pull(Fraud) -> pd2

Vì mức độ chính xác cho dự báo một giao dịch gian lận thực sự là gian lận (dương tính thật) phụ thuộc vào ngưỡng (Threshold hay Cutoff) nên chúng ta khảo sát chất lượng của hai mô hình đồng thời trên một loạt ngưỡng được lựa chọn:

# Function calculates model performances by cutoff for model 1: 

byCutoff_rf1 <- function(cutoff) {
  
  pred <- case_when(pd1 >= cutoff ~ "Fraud", TRUE ~ "NonFraud") %>% as.factor()
  
  confusionMatrix(pred, test_benford$Class, positive = "Fraud") -> cm
  
  cm$table %>% as.vector() -> bg
  cm$overall %>% as.vector() -> acc
  cm$byClass %>% as.vector() -> sen
  
  data.frame(BB = bg[1], 
             BG = bg[2], 
             GB = bg[3], 
             GG = bg[4], 
             Accuracy = acc[1],
             Kappa = acc[2],
             Recall = sen[1],
             Specificity = sen[2], 
             Cutoff = cutoff) -> model_perCutoff
  
  return(model_perCutoff)
  
}


# Function calculates model performances by cutoff for model 2: 

byCutoff_rf2 <- function(cutoff) {
  
  pred <- case_when(pd2 >= cutoff ~ "Fraud", TRUE ~ "NonFraud") %>% as.factor()
  
  confusionMatrix(pred, test_benford$Class, positive = "Fraud") -> cm
  
  cm$table %>% as.vector() -> bg
  cm$overall %>% as.vector() -> acc
  cm$byClass %>% as.vector() -> sen
  
  data.frame(BB = bg[1],
             BG = bg[2],
             GB = bg[3],
             GG = bg[4], 
             Accuracy = acc[1],
             Kappa = acc[2],
             Recall = sen[1], 
             Specificity = sen[2], 
             Cutoff = cutoff) -> model_perCutoff
  
  return(model_perCutoff)
  
}

Rồi sử dụng hai hàm trên:

# Use this fuction: 

lapply(seq(0.05, 0.75, 0.05), byCutoff_rf1) -> list1
lapply(seq(0.05, 0.75, 0.05), byCutoff_rf2) -> list2

# Convert to data frame and combine the two datasets: 

do.call("bind_rows", list1) -> df1
do.call("bind_rows", list2) -> df2

bind_rows(df1 %>% mutate(`Feature Engineering` = "Benford"), 
          df2 %>% mutate(`Feature Engineering` = "None")) -> df_comparision

So sánh chất lượng của hai mô hình phân loại tương ứng với một loạt ngưỡng được lựa chọn:

# Compare model performance by Feature Engineering Method: 

my_colors <- c("firebrick", "steelblue")

df_comparision %>% 
  select(5:10) %>% 
  gather(Metric, Value, -Cutoff, -`Feature Engineering`) %>% 
  ggplot(aes(Cutoff, Value, color = `Feature Engineering`)) + 
  geom_line() + 
  geom_point() + 
  scale_color_manual(values = my_colors) + 
  facet_wrap(~ Metric, scales = "free") + 
  theme(legend.position = "top") + 
  my_theme() + 
  scale_y_continuous(labels = scales::percent) + 
  scale_x_continuous(breaks = seq(0.05, 0.75, 0.1)) + 
  labs(x = NULL, y = NULL, 
       title = "Figure 7: Model Performance by Feature Engineering Method", 
       caption = "Source: Credit Card Fraud Detection (https://www.kaggle.com)")

Figure 7 cho thấy tại tất cả các ngưỡng chất lượng phân loại (đặc biệt là khả năng xác minh đúng giao dịch gian lận) của mô hình sử dụng features được extract từ Amount dựa trên Benford law là tốt hơn.

Conclusion

Các cơ quan chuyên về điều tra các gian lận về thuế ở châu Âu đã công khai thừa nhận việc áp dụng định luật Benford cho điều tra kĩ hơn các cases gian lận. Các công ti kiểm toán (Auditing Firms) cũng áp dụng định luật này trong các nghiệp vụ liên quan đến hoạt động kiểm toán các số liệu tài chính bị nghi ngờ là giả mạo hoặc bị bóp méo. Ở Mĩ thì IRS (cơ quan chuyên thuế) không phủ nhận nhưng cũng không thừa nhận việc sử dụng định luật này để xác minh các gian lận về thuế. Bài viết này trình bày một hướng ứng dụng mới của định luật Benford như là một kĩ thuật cho Feature Engineering/Extraction. Sử dụng Logistic Classifier cho bộ dữ liệu creditcard.csv cho thấy rằng các features được tạo ra bằng cách vận dụng định luật Benford sẽ nâng cao tỉ lệ xác minh các giao dịch tài chính gian lận. Tuy nhiên hiệu lực của kĩ thuật feature engineering này cho các mô hình học máy khác hay cách tinh vi hơn của việc vận dụng định luật Benford cho fearure engineering là chưa được đề cập trong bài viết này.

References

1. An Introduction to Benford’s Law.
2. Application of Statistical Methods for Tax Inspection of Enterprises: A Case Study in Vietnam.
3. Benford’s Law: Theory, The General Law Of Relative Quantities, And Forensic Fraud Detection Applications.
4. Forensic Analytics: Methods and Techniques for Forensic Accounting Investigations.
5. Benford’s Law: Applying to Existing Data.