Phát hiện bất thường bằng Cluster

1 Mở đầu câu chuyện

Cách đây vài ngày, GS. Nguyễn Văn Tuấn đã vận dụng lý thuyết xác suất và kiểm định thống kê để khẳng định có sự bất thường trong phân bố điểm thi ở tỉnh Hà Giang so với dữ liệu toàn quốc.

Trong bài thực hành hôm nay, Nhi sẽ thử áp dụng một cách “điều tra” khác, khách quan và dựa vào dữ liệu quan sát thực tế, với mục tiêu nhận diện những trường hợp và đặc điểm phân bố bất thường về điểm thi tại Hà Giang. Phương pháp thống kê sẽ được sử dụng là “K-mean clustering” (Phân tích cụm).

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 3.5.1

## -- Attaching packages ----------------------------- tidyverse 1.2.1 --

## v ggplot2 3.0.0     v purrr   0.2.5
## v tibble  1.4.2     v dplyr   0.7.6
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0

## Warning: package 'ggplot2' was built under R version 3.5.1

## Warning: package 'tibble' was built under R version 3.5.1

## Warning: package 'tidyr' was built under R version 3.5.1

## Warning: package 'readr' was built under R version 3.5.1

## Warning: package 'purrr' was built under R version 3.5.1

## Warning: package 'dplyr' was built under R version 3.5.1

## Warning: package 'stringr' was built under R version 3.5.1

## Warning: package 'forcats' was built under R version 3.5.1

## -- Conflicts -------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Đầu tiên, Nhi tải bộ dữ liệu điểm thi cho từng thí sinh ở Hà Giang (xin cảm ơn Gs. Tuấn đã công khai dữ liệu này).

hg_data=read.csv("THPT 2018 Ha Giang.csv")%>%as_tibble()
vn_data=read.csv("THPT 2018 Quoc gia.csv")%>%as_tibble()

Câu chuyện bắt đầu bằng những lời bàn tán trên mạng về phân bố điểm thi rất kì dị ở Hà Giang, như chúng ta có thể thấy trong biểu đồ sau:

library(ggridges)

## Warning: package 'ggridges' was built under R version 3.5.1

## 
## Attaching package: 'ggridges'

## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual

library(viridis)

## Warning: package 'viridis' was built under R version 3.5.1

## Loading required package: viridisLite

## Warning: package 'viridisLite' was built under R version 3.5.1

hg_data%>%
  gather(Math:Geography,key="Exam",value="Score")%>%
  ggplot()+
  geom_density_ridges_gradient(aes(x=Score,y=reorder(Exam,Score),fill=..x..),
                               scale=1,
                               show.legend = T)+
  scale_fill_viridis(option="D")+
  labs(y="Exams")+
  theme_bw()+ggtitle("Ha_Giang")

## Picking joint bandwidth of 0.305

Ở các môn: Toán, Vật lý, ngoại ngữ, Hóa học, mật độ điểm 9 cao một cách bất thường, so với phân bố chung của cả nước:

vn_data%>%
  gather(Math:Geography,key="Exam",value="Score")%>%
  ggplot()+
  geom_density_ridges_gradient(aes(x=Score,y=reorder(Exam,Score),fill=..x..),
                               scale=1,
                               show.legend = T)+
  scale_fill_viridis(option="C",direction = -1)+  
  labs(y="Exams")+
  theme_bw()+ggtitle("National")

## Picking joint bandwidth of 0.088

Với dữ liệu tóm tắt cho từng môn thi riêng lẻ như báo chí đã sử dụng, ta chỉ có thể quan sát đặc tính phân bố, đặt nghi vấn và dừng lại ở đó. Dữ liệu chi tiết toàn bộ môn thi, đến từng cá thể mà Gs. Tuấn công bố cho phép chúng ta đi xa hơn, thí dụ với phương pháp Clustering mà Nhi sắp thực hiện.

2 Phân tích cụm

Trước hết, Nhi đặt ra giả định như sau:

Tổ hợp giữa các môn Toán, Vật Lý, Hóa học, Sinh học có ý nghĩa quan trọng, quyết định khả năng vào Đại học của thí sinh, nên một học sinh nếu đã đặt ra kế hoạch nhắm vào một trường Đại học nào đó, chắc chắn họ sẽ phải đầu tư tốt nhất có thể và đồng đều vào CẢ 3 môn: Toán-Lý-Hóa hoặc Toán-Hóa-Sinh.

Các môn Toán, Lý, Hóa, Sinh đều thuộc ngành khoa học tự nhiên, chúng có chung yêu cầu về năng lực trí tuệ, phương pháp học tập và kỹ năng làm bài thi.

Vì 2 lý do này, nếu mọi sự diễn ra thuận theo tự nhiên, có nhiều khả năng điểm của 3/4 môn thi trong tổ hợp Toán/Lý/Hóa/Sinh phải tương quan chặt chẽ với nhau.

Do đó, Nhi sẽ trích xuất dữ liệu của riêng 4 môn thi nói trên để làm Clustering:

sub_hg_data=hg_data%>%select(Math,Physics,Chemistry,Biology)%>%
  na.omit()

head(sub_hg_data)

## # A tibble: 6 x 4
##    Math Physics Chemistry Biology
##   <dbl>   <dbl>     <dbl>   <dbl>
## 1   5.6    6.5       3.75    3.5 
## 2   3.8    2.75      2.75    2.75
## 3   2.4    2         3.5     3.25
## 4   9.4    9         3.25    3   
## 5   2.8    3         3.75    3.5 
## 6   3.6    2.5       4       4.75

Có 630 thí sinh ở Hà Giang đã thi đồng thời 4 môn Toán-Lý-Hóa-Sinh, dữ liệu này khá đủ cho một phân tích cụm:

Sau khi thăm dò số cluster cho phương pháp K-means, Nhi quyết định chọn số lượng cluster = 5 (đây là tham số chủ quan duy nhất trong quy trình)

fviz_nbclust(sub_hg_data, kmeans, method = "wss") +
  theme_bw()+
  geom_vline(xintercept =5, linetype = 2)

Việc phân cụm rất dễ dàng chỉ với 1 dòng code, kết quả được lưu lại trong object km.res

km.res <- kmeans(sub_hg_data, 5, nstart = 25)

Kết quả phân cụm có thể được biểu diễn trực quan như sau:

fviz_cluster(km.res,
             data=sub_hg_data,
             ellipse.type = "t", 
             ggtheme = theme_classic()
)

Để hiểu rõ bản chất của mỗi cluster là gì, ta thử tính Mode của điểm 4 môn thi cho từng cluster:

getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

aggregate(sub_hg_data, by=list(cluster=km.res$cluster), getmode)

##   cluster Math Physics Chemistry Biology
## 1       1  9.4    9.50      3.25    2.00
## 2       2  2.6    2.50      2.75    3.25
## 3       3  9.0    9.25      9.00    1.75
## 4       4  6.8    6.00      5.50    4.75
## 5       5  4.8    3.75      4.25    4.75

Ta cũng có thể dùng Trung bình để suy diễn:

aggregate(sub_hg_data, by=list(cluster=km.res$cluster), mean)

##   cluster     Math  Physics Chemistry  Biology
## 1       1 8.973913 9.184783  2.945652 2.554348
## 2       2 2.974194 2.620968  2.707373 3.443548
## 3       3 8.782353 9.022059  9.066176 3.007353
## 4       4 6.540260 5.943182  5.905844 4.951299
## 5       5 4.826733 3.883663  4.134901 4.518564

Kết quả này cho thấy:

1. Cluster thứ 1 là những trường hợp đáng ngờ nhất vì tổ hợp điểm thi rất kì lạ : môn Toán và Vật lý rất cao, nhưng môn Hóa và Sinh học đều có điểm cực thấp.

2. Cluster thứ 2 bao gồm những học sinh loại yếu kém, với điểm cả 4 môn thi đều thấp như nhau, những trường hợp này cũng có thể xem là outliers trong dữ liệu, tuy nhiên lại không bị nghi ngờ về sự gian lận.

3. Cluster thứ 3 có vẻ đại diện cho những thí sinh có sở trường Khối A, với tổ hợp Toán-Lý-Hóa điểm rất cao, nhưng điểm môn Sinh học lại cực kì thấp, đây có thể, nhưng chưa hẳn là điều đáng ngờ, vì có nhiều học sinh không có lập trường xác định, họ thi cả 2 khối A,B một cách cầu may.

4. Cluster thứ 4 tương ứng với những thí sinh có học lực trung bình, khá với điểm thi cả 2 môn Toán, Lý đều nhau và ở mức khá,môn Hóa và Sinh thì ở mức trung bình. Có thể loại nhóm này khỏi diện nghi ngờ.

5. Cluster thứ 5 đại diện cho các thí sinh có học lực trung bình, yếu cho cả 4 môn, ta cũng có thể loại nhóm này khỏi sự nghi ngờ.

Ta cùng so sánh phân bố điểm cho từng môn thi giữa 5 Clusters:

library(lvplot)

## Warning: package 'lvplot' was built under R version 3.5.1

sub_hg_data%>%mutate(cluster=factor(km.res$cluster))%>%
  gather(Math:Biology,key="Exam",value="Score")%>%
  ggplot()+
  geom_lv(aes(x=Exam,y=Score,fill=..LV..),col="black",show.legend = F,)+
  facet_wrap(~cluster,ncol=2)+
  coord_flip()+
  scale_fill_brewer(palette="Reds",direction = -1)+
  theme_bw()

## Warning: package 'bindrcpp' was built under R version 3.5.1

3 Nhận diện cá thể bất thường

Để nhận diện những trường hợp bất thường, ta có thể: 1) dựa hoàn toàn vào kết quả K-mean cluster, ở đây là toàn bộ cluster 1 và cluster 3; hoặc 2) qua trung gian một chỉ số thống kê là “khoảng cách”trong không gian dữ liệu, hoặc kết hợp cả 2 tiêu chí.

Thí dụ nếu Nhi dùng 2 chỉ số khoảng cách Euclide và Manhattan, có thể thấy tỉ lệ thí sinh có giá trị khoảng cách cao bất thường trong dữ liệu hiện thời :

sub_hg_data%>%get_dist(method = "euclidean")%>%
  fviz_dist(gradient = list(low = "white", mid = "gold", high = "red"))+
  scale_x_discrete(labels = NULL,breaks=NULL)+
  scale_y_discrete(labels = NULL,breaks=NULL)+
  coord_equal()+
  ggtitle("Euclidian distance")

sub_hg_data%>%get_dist(method = "manhattan")%>%
  fviz_dist(gradient = list(low = "white", mid = "pink", high = "purple"))+
  scale_x_discrete(labels = NULL,breaks=NULL)+
  scale_y_discrete(labels = NULL,breaks=NULL)+
  coord_equal()+
  ggtitle("Manhattan distance")

Những vùng màu đậm nhất tương ứng với những thí sinh có thể xem là outliers (ghi chú: điểm thi của họ có thể rất thấp, hoặc rất cao so với những quan sát xung quanh).

Nếu ấn định một tỉ lệ nhất định các trường hợp có khoảng cách cao nhất, thí dụ 5% ta có thể phân lập những cá thể thuộc diện tình nghi này khỏi quần thể một cách trực quan như thí dụ minh họa sau

centers <- km.res$centers[km.res$cluster, ]

distances <- sqrt(rowSums((sub_hg_data - centers)^2)) 

outliers <- order(distances, decreasing=T)[1:30]

plot_df<-mutate(sub_hg_data,cluster=factor(km.res$cluster))

out_df<-plot_df%>%.[outliers,]

ggplot()+
  geom_point(data=plot_df,aes(x=Math,y=Physics,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Math,y=Physics,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Math,y=Physics),shape=21,size=6,stroke=1,col="red",alpha=0.9)+
  theme_bw()

ggplot()+
  geom_point(data=plot_df,aes(x=Math,y=Chemistry,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Math,y=Chemistry,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Math,y=Chemistry),shape=21,size=6,stroke=1,col="red",alpha=0.9)+
  theme_bw()

ggplot()+
  geom_point(data=plot_df,aes(x=Math,y=Biology,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Math,y=Biology,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Math,y=Biology),shape=21,size=6,stroke=1,col="red",alpha=0.9)+
  theme_bw()

ggplot()+
  geom_point(data=plot_df,aes(x=Chemistry,y=Biology,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Chemistry,y=Biology,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Chemistry,y=Biology),shape=21,size=6,stroke=1,col="red",alpha=0.9)+
  theme_bw()

Ta cũng có thể phân lập toàn bộ clusters 1 và 3 , được xem như “tình nghi cao”:

outliers <- order(distances, decreasing=T)[1:500]

out_df<-plot_df%>%.[outliers,]%>%filter(cluster== "1" | cluster =="3")

ggplot()+
  geom_point(data=plot_df,aes(x=Math,y=Physics,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Math,y=Physics,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Math,y=Physics),shape=23,size=6,stroke=1,col="blue",alpha=0.9)+
  theme_bw()

ggplot()+
  geom_point(data=plot_df,aes(x=Math,y=Chemistry,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Math,y=Chemistry,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Math,y=Chemistry),shape=23,size=6,stroke=1,col="blue",alpha=0.9)+
  theme_bw()

ggplot()+
  geom_point(data=plot_df,aes(x=Math,y=Biology,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Math,y=Biology,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Math,y=Biology),shape=23,size=6,stroke=1,col="blue",alpha=0.9)+
  theme_bw()

ggplot()+
  geom_point(data=plot_df,aes(x=Chemistry,y=Biology,col=cluster),alpha=0.5)+
  geom_point(data=as.data.frame(km.res$centers),aes(x=Chemistry,y=Biology,col=factor(c(1:5))),shape=15,size=3)+
  geom_point(data=out_df,aes(x=Chemistry,y=Biology),shape=23,size=6,stroke=1,col="blue",alpha=0.9)+
  theme_bw()

4 Kết luận

Qua bài thực hành này, các bạn đã làm quen với phương pháp nhận diện giá trị bất thường trong dữ liệu là K-means clustering analysis và khoảng cách Euclide.

Điểm thú vị của phương pháp K-means Cluster, đó là tính khách quan, không phụ thuộc vào bất cứ giả thuyết chủ quan nào của người điều tra, khác với cách làm kiểm định giả thuyết vô hiệu hoặc suy diễn dựa vào biểu đồ. Việc suy diễn chỉ được làm sau khi có kết quả phân cụm, và kết quả này dựa vào dữ liệu.

P.S: Chú ý một điều: Cả 2 phương pháp không thể khẳng định điều gì về sự can thiệp hay gian lận, nhưng có thể nhận diện những đối tượng tình nghi.