Bài toán phân loại đa giá trị

Keras v.s RandomForest

1 Giới thiệu

Trong bài thực hành này, chúng ta tiếp tục khám phá giải thuật Deep learning sử dụng TensorFlow qua giao thức Keras trong R. Nhi sẽ dựng một mô hình Neural network để giải quyết một bài toán phân lớp đa giá trị (đến 10 nhãn) bằng Keras, sau đó so sánh độ chính xác của nó với một mô hình Random Forest.

2 Dữ liệu

Bài toán của chúng ta là một nghiên cứu có thực của tác giả Ayres de Campos et al. (Khoa Y, Đại học Porto, Bồ Đào Nha) vào năm 2000. Các tác giả tạo ra một phần mềm SisPorto 2.0 cho phép phân tích Tim thai tự động. Bộ dữ liệu của họ gồm 2126 trường hợp với 20 features và 10 nhãn phân loại chẩn đoán (kiểu tim thai) khác nhau.

Dữ liệu này có thể tải từ database UCI:

https://archive.ics.uci.edu/ml/datasets/Cardiotocography

Như thường lệ, Nhi dùng caret để cắt ngẫu nhiên dữ liệu thành 2 phần,mỗi phần chứa 50% (n=1060). Sau đó Nhi sẽ không nhìn đến tập kiểm định (testset) nhưng thăm dò tập huấn luyện (trainset) một cách trực quan bằng violin plot.

library(tidyverse)

my_theme <- function(base_size = 5, base_family = "sans"){
  theme_minimal(base_size = base_size, base_family = base_family) +
    theme(
      axis.text = element_text(size = 5),
      axis.text.x = element_text(angle = 0, vjust = 0.5, hjust = 0.5),
      axis.title = element_text(size = 5),
      panel.grid.major = element_line(color = "grey"),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = "white"),
      strip.background = element_rect(fill = "#0f0016", color = "#0f0016", size =0.5),
      strip.text = element_text(face = "bold", size =5, color = "white"),
      legend.position = "bottom",
      legend.justification = "center",
      legend.background = element_blank(),
      panel.border = element_rect(color = "grey30", fill = NA, size = 0.5)
    )
}
theme_set(my_theme())

df=read.csv("Cardiotocography.csv",sep=";")%>%.[,-c(21:31)]%>%.[,-22]

df[c(1:20)] <- lapply(df[c(1:20)], as.numeric)

df$CLASS%<>%as.factor()%>%
  recode_factor(.,`1` = "CalmSleep", 
                `2` = "REM",
                `3` = "CalmVig",
                `4` = "ActVig", 
                `5` = "Shiftpattern",
                `6` = "AD",
                `7` = "DE", 
                `8` = "LD",
                `9` = "Flatsinusoidal",
                `10` = "Suspect")

mycol=c("#ff0048","#fce628","#0de2aa",
        "#2dd6e2","#c362f7","#ffae22",
        "#b1e827","#3495ea","#ac9af9","#ff51c5")

library(caret)
set.seed(1234)

idTrain=createDataPartition(y=df$CLASS,p=0.5,list=FALSE)
trainset=df[idTrain,]
testset=df[-idTrain,]

head(trainset)%>%knitr::kable()

	LB	AC	FM	UC	DL	DS	DP	ASTV	MSTV	ALTV	MLTV	Width	Min	Max	Nmax	Nzeros	Mode	Mean	Median	Variance	CLASS
2	132	0.006	0	0.006	0.003	0	0.000	17	2.1	0	10.4	130	68	198	6	1	141	136	140	12	AD
3	133	0.003	0	0.008	0.003	0	0.000	16	2.1	0	13.4	130	68	198	5	1	141	135	138	13	AD
4	134	0.003	0	0.008	0.003	0	0.000	16	2.4	0	23.0	117	53	170	11	0	137	134	137	13	AD
5	132	0.007	0	0.008	0.000	0	0.000	16	2.4	0	19.9	117	53	170	9	0	137	136	138	11	REM
7	134	0.001	0	0.013	0.008	0	0.003	29	6.3	0	0.0	150	50	200	6	3	71	107	106	215	LD
8	122	0.000	0	0.000	0.000	0	0.000	83	0.5	6	15.6	68	62	130	0	0	122	122	123	3	Flatsinusoidal

trainset%>%gather(LB:Variance,key="Feature",value="Value",na.rm=T)%>%
  ggplot(aes(x=CLASS,y=Value,fill=CLASS))+
  geom_violin(color="black",alpha=0.7,show.legend = F)+
  facet_wrap(~Feature,drop = TRUE,ncol=3,scales="free")+
  my_theme(5)+scale_fill_manual(values=mycol)

3 Mô hình Deep neural network với Keras

Bây giờ, ta sẽ thực hiện một quy trình dựng mô hình feed forward neural network bằng Keras. Như đã giới thiệu trong bài trước, Keras cho phép thực hành Deep learning với TensorFlow một cách thoải mái với cú pháp có tính thứ bậc, rõ ràng dễ hiểu (ngay cả khi máy tính của bạn không có GPU xịn thì tốc độ training cho các mô hình nhỏ trên CPU cũng khá nhanh)

Trở ngại lớn nhất khi bắt đầu dùng Keras không phải ở cú pháp hay cấu trúc các Layer, mà ở khâu chuẩn bị dữ liệu: keras chỉ tiếp nhận một cấu trúc dữ liệu đặc biệt gọi là tensor. Vì Nhi chỉ là 1 bác sĩ nên phải mất mấy ngày Nhi mới hiểu được khái niệm tensor và các phép ncho nó. Trên youtube có nhiều thầy dạy về những thứ này rất dễ hiểu. Sau đó bạn sẽ khám phá là tensor được lưu trong R dưới dạng những arrays nên toàn bộ dữ liệu phải được chuyển thành array, việc này mất thời gian và cần sự cẩn thận, nhưng không khó khăn lắm. Cuối cùng, bạn cần làm quen với 3 thủ thuật khi tiếp cận bài toán phân lớp trong keras: thứ nhất : tensor kết quả (1D) phải được tách rời khỏi tensor features (2D hay 3D), thứ hai là các nhãn giá trị phải được chuyển thành dummy variables với giá trị 0 và 1; thứ ba đó là toàn bộ features phải được xử lý hoán chuyển thành con số và chuẩn hóa.

Quy trình chuẩn bị dữ liệu được làm từng bước như sau:

library(keras)
library(tensorflow)

# A function for normalising features

normalize <- function(traindf,testdf) {
  mu <- apply(traindf,2,mean)
  sigma <- apply(traindf,2,sd)
  normtrain <- scale(traindf, center = mu, scale = sigma)
  normtest<- scale(testdf, center = mu, scale = sigma)
  outlist<-list(normtrain,normtest)
  return (outlist)
}

# Normalisation
outlist<-normalize(traindf=trainset[,-21],testdf=testset[,-21])

normtrain<-outlist[[1]]%>%as.data.frame()
normtest<-outlist[[2]]%>%as.data.frame()

# Make 1D tensor with dummy labels
train.lab<-to_categorical(trainset$CLASS)%>%.[,-1]
test.lab<-to_categorical(testset$CLASS)%>%.[,-1]

# Convert dataframe to array

train_data<-normtrain%>%as.matrix()%>%as.array.default(dimnames=NULL)  
dimnames(train_data) <- NULL

test_data<-normtest%>%as.matrix()%>%as.array.default(dimnames=NULL) 
dimnames(test_data) <- NULL

Đây là 1D tensor chứa 10 nhãn dummies dùng huấn luyện

is.array(train.lab)

## [1] TRUE

head(train.lab)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    0    0    0    0    0    1    0    0    0     0
## [2,]    0    0    0    0    0    1    0    0    0     0
## [3,]    0    0    0    0    0    1    0    0    0     0
## [4,]    0    1    0    0    0    0    0    0    0     0
## [5,]    0    0    0    0    0    0    0    1    0     0
## [6,]    0    0    0    0    0    0    0    0    1     0

Còn đây là 2D tensor chứa 20 features đã được chuẩn hóa

str(train_data)

##  num [1:1066, 1:20] -0.1589 -0.0581 0.0426 -0.1589 0.0426 ...

head(train_data)

##             [,1]        [,2]       [,3]       [,4]       [,5]        [,6]
## [1,] -0.15893720  0.74351612 -0.2049871  0.5515012  0.3865457 -0.07520018
## [2,] -0.05814776 -0.03838607 -0.2049871  1.2319413  0.3865457 -0.07520018
## [3,]  0.04264169 -0.03838607 -0.2049871  1.2319413  0.3865457 -0.07520018
## [4,] -0.15893720  1.00415019 -0.2049871  1.2319413 -0.6367764 -0.07520018
## [5,]  0.04264169 -0.55965420 -0.2049871  2.9330415  2.0920824 -0.07520018
## [6,] -1.16683162 -0.82028826 -0.2049871 -1.4898191 -0.6367764 -0.07520018
##            [,7]      [,8]       [,9]      [,10]      [,11]       [,12]
## [1,] -0.2606678 -1.738683  0.8683459 -0.5425538  0.3928203  1.51960804
## [2,] -0.2606678 -1.796799  0.8683459 -0.5425538  0.9264819  1.51960804
## [3,] -0.2606678 -1.796799  1.2048667 -0.5425538  2.6341989  1.18892231
## [4,] -0.2606678 -1.796799  1.2048667 -0.5425538  2.0827486  1.18892231
## [5,]  4.5302266 -1.041291  5.5796377 -0.5425538 -1.4572064  2.02835530
## [6,] -0.2606678  2.096974 -0.9264320 -0.2199869  1.3178337 -0.05750849
##          [,13]      [,14]      [,15]      [,16]       [,17]        [,18]
## [1,] -0.870266  1.8598883  0.6689788  1.0108018  0.19524090  0.064834109
## [2,] -0.870266  1.8598883  0.3327547  1.0108018  0.19524090  0.001659374
## [3,] -1.371330  0.3139695  2.3500996 -0.4452992 -0.03999708 -0.061515361
## [4,] -1.371330  0.3139695  1.6776513 -0.4452992 -0.03999708  0.064834109
## [5,] -1.471543  1.9703111  0.6689788  3.9230039 -3.92142375 -1.767233196
## [6,] -1.070692 -1.8944861 -1.3483661 -0.4452992 -0.92213950 -0.819612176
##            [,19]      [,20]
## [1,]  0.11002225 -0.2356141
## [2,] -0.02463185 -0.2025965
## [3,] -0.09195889 -0.2025965
## [4,] -0.02463185 -0.2686316
## [5,] -2.17909731  6.4669513
## [6,] -1.03453753 -0.5327721

Xong phần dữ liệu, bây giờ ta phác thảo cấu trúc của mạng neuron: Do đây chỉ là thí dụ minh họa, Nhi lựa chọn các thông số kỹ thuật một cách ngẫu nhiên mà không thông qua giai đoạn kiểm định chéo. Nhi sẽ dựng một mạng neuron với 5 lớp, lớp thứ nhất chứa 256 neuron,dữ liệu đầu vào gồm 20 biến, lớp thứ 2,3,4 đều chứa 256 neuron và có áp dụng L2 regularisation , lớp cuối cùng xuất ra kết quả với 10 neuron, sử dụng hàm softmax.

Mô hình này được đóng gói với hàm loss là categorical_crossentropy, hàm optimiser là “adam”, và tiêu chí huấn luyện là accuracy.

build_model <- function() {
  model <- keras_model_sequential() %>% 
    layer_dense(units = 256, 
                activation = "relu", 
                input_shape = dim(train_data)[[2]])%>%
    layer_dense(units =256,kernel_regularizer = regularizer_l2(0.001),
                activation = "relu")%>%
    layer_dense(units =256,kernel_regularizer = regularizer_l2(0.001),
                activation = "relu")%>%
    layer_dense(units =256, kernel_regularizer = regularizer_l2(0.001),
                activation = "relu")%>% 
    layer_dense(units = 10, activation = "softmax")
  
  model %>%compile(
    loss = 'categorical_crossentropy',
    optimizer = 'adam',
    metrics = 'accuracy'
  )
}

Mô hình được huấn luyện với 500 lượt, mỗi lần sử dụng 100 đơn vị dữ liệu, và kiểm định trên 20%

model<-build_model()

tocomod<-model%>% fit(
  train_data, 
  train.lab, 
  epochs = 500, 
  batch_size = 100, 
  validation_split = 0.2
)

Bên phải màn hình Rstudio, bạn sẽ thấy biểu đồ diễn tiến của ACCuracy và hàm Loss theo thời gian thực (đây là tính năng thú vị mà hiếm package machine learning nào cung cấp).

Sau khi mô hình được huấn luyện xong, ta có thể vẽ lại biểu đồ này:

plot(tocomod)+my_theme()

Kết quả có vẻ khả quan, nhưng thứ mà ta đang nhìn thấy có thể là “overfitting”, do đó ta phải kiểm chứng lại trên test data:

Ta làm một confussion matrix:

PredLab<-predict_classes(model,test_data)

PredLab%<>%as.factor()%>%
  recode_factor(.,`0` = "CalmSleep", 
                `1` = "REM",
                `2` = "CalmVig",
                `3` = "ActVig", 
                `4` = "Shiftpattern",
                `5` = "AD",
                `6` = "DE", 
                `7` = "LD",
                `8` = "Flatsinusoidal",
                `9` = "Suspect")

confusionMatrix(PredLab,reference=testset$CLASS,
                mode="everything")

## Confusion Matrix and Statistics
## 
##                 Reference
## Prediction       CalmSleep REM CalmVig ActVig Shiftpattern  AD  DE  LD
##   CalmSleep            137  24       7      0            8   2   4   0
##   REM                   11 238       0      2            0  10   0   0
##   CalmVig               10   0      19      0            0   2   1   0
##   ActVig                 0  13       0     38            0  13   0   0
##   Shiftpattern           2   3       0      0           22   0   0   0
##   AD                     1   9       0      0            0 112  10   5
##   DE                     3   0       0      0            0  24 111  31
##   LD                     0   0       0      0            0   3   0  17
##   Flatsinusoidal         3   0       0      0            0   0   0   0
##   Suspect               25   2       0      0            6   0   0   0
##                 Reference
## Prediction       Flatsinusoidal Suspect
##   CalmSleep                   0      25
##   REM                         0       0
##   CalmVig                     0       0
##   ActVig                      0       0
##   Shiftpattern                0       3
##   AD                          0       0
##   DE                          0       1
##   LD                          0       0
##   Flatsinusoidal             31       4
##   Suspect                     3      65
## 
## Overall Statistics
##                                           
##                Accuracy : 0.7453          
##                  95% CI : (0.7179, 0.7713)
##     No Information Rate : 0.2726          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.698           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: CalmSleep Class: REM Class: CalmVig
## Sensitivity                    0.7135     0.8235        0.73077
## Specificity                    0.9194     0.9702        0.98743
## Pos Pred Value                 0.6618     0.9119        0.59375
## Neg Pred Value                 0.9355     0.9362        0.99319
## Precision                      0.6618     0.9119        0.59375
## Recall                         0.7135     0.8235        0.73077
## F1                             0.6867     0.8655        0.65517
## Prevalence                     0.1811     0.2726        0.02453
## Detection Rate                 0.1292     0.2245        0.01792
## Detection Prevalence           0.1953     0.2462        0.03019
## Balanced Accuracy              0.8164     0.8968        0.85910
##                      Class: ActVig Class: Shiftpattern Class: AD Class: DE
## Sensitivity                0.95000             0.61111    0.6747    0.8810
## Specificity                0.97451             0.99219    0.9720    0.9368
## Pos Pred Value             0.59375             0.73333    0.8175    0.6529
## Neg Pred Value             0.99799             0.98641    0.9415    0.9831
## Precision                  0.59375             0.73333    0.8175    0.6529
## Recall                     0.95000             0.61111    0.6747    0.8810
## F1                         0.73077             0.66667    0.7393    0.7500
## Prevalence                 0.03774             0.03396    0.1566    0.1189
## Detection Rate             0.03585             0.02075    0.1057    0.1047
## Detection Prevalence       0.06038             0.02830    0.1292    0.1604
## Balanced Accuracy          0.96225             0.80165    0.8234    0.9089
##                      Class: LD Class: Flatsinusoidal Class: Suspect
## Sensitivity            0.32075               0.91176        0.66327
## Specificity            0.99702               0.99318        0.96258
## Pos Pred Value         0.85000               0.81579        0.64356
## Neg Pred Value         0.96538               0.99706        0.96559
## Precision              0.85000               0.81579        0.64356
## Recall                 0.32075               0.91176        0.66327
## F1                     0.46575               0.86111        0.65327
## Prevalence             0.05000               0.03208        0.09245
## Detection Rate         0.01604               0.02925        0.06132
## Detection Prevalence   0.01887               0.03585        0.09528
## Balanced Accuracy      0.65889               0.95247        0.81292

Kết quả kiểm định trên test_data không tệ lắm đối với một bài toán 10 nhãn !

Do confussion matrix này tương đối khó đọc, ta sẽ chuyển nó thành 1 hình vẽ bằng một hàm như sau:

Cuối cùng, ta dựng một mô hình chính thức với số epochs = 50 và batch_size=16

multiclasscmplot=function(Labs,Truth){
  cm=caret::confusionMatrix(Labs,reference=Truth)
  cmt=cm$table%>%as_tibble()
  freq=Hmisc::describe(Truth)%>%.$values%>%.$frequency
  labels=row.names(cm$table)
  
  k=1
  cmt$rate=cmt$n/1
  for(i in 1:10){
    for(k in 1:100){
      if(cmt$Reference[k]==labels[i]){
        cmt$rate[k]<-100*(cmt$n[k]/freq[i])
        k=k+1
      }
      else{k=k}
    }
  }
  
  cmt%>%ggplot(aes(x=Prediction,y=Reference,fill=rate))+
    geom_tile(color="black")+my_theme(10)+
    scale_fill_gradient2(low="white",mid="#008cff",high="#ff0059",midpoint=50)+
    theme(axis.text.x = element_text(angle =45, hjust = 0.5))+
    geom_text(aes(y=Reference,x=Prediction,label=round(rate,1)),color="black",size=3)+coord_fixed()
}

p1=multiclasscmplot(PredLab,testset$CLASS)+ggtitle("Keras")

p1

4 Random Forest

Tiếp theo, ta sẽ cho mô hình Deep neural net đấu với một mô hình Random Forest. Việc này vô cùng dễ dàng và nhanh chóng như ta từng làm nhiều lần trước đây:

library(randomForest)

rfmod=randomForest(CLASS~.,
                   data=trainset,
                   ntree=500,
                   mtry=5,
                   replace=TRUE
)

confusionMatrix(reference=testset$CLASS,
                data=predict(rfmod,newdata=testset),
                mode="everything")

## Confusion Matrix and Statistics
## 
##                 Reference
## Prediction       CalmSleep REM CalmVig ActVig Shiftpattern  AD  DE  LD
##   CalmSleep            169   8       7      0            8   1   2   0
##   REM                   14 270       3      8            4  15   0   0
##   CalmVig                4   0      16      0            0   0   0   0
##   ActVig                 0   4       0     32            0   0   0   0
##   Shiftpattern           1   0       0      0           16   0   0   0
##   AD                     0   4       0      0            0 141   9   0
##   DE                     3   1       0      0            0   7 113   3
##   LD                     0   0       0      0            0   2   2  50
##   Flatsinusoidal         1   0       0      0            0   0   0   0
##   Suspect                0   2       0      0            8   0   0   0
##                 Reference
## Prediction       Flatsinusoidal Suspect
##   CalmSleep                   2       8
##   REM                         0       1
##   CalmVig                     0       0
##   ActVig                      0       0
##   Shiftpattern                0       0
##   AD                          0       0
##   DE                          0       1
##   LD                          0       0
##   Flatsinusoidal             29       2
##   Suspect                     3      86
## 
## Overall Statistics
##                                           
##                Accuracy : 0.8698          
##                  95% CI : (0.8481, 0.8895)
##     No Information Rate : 0.2726          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8436          
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: CalmSleep Class: REM Class: CalmVig
## Sensitivity                    0.8802     0.9343        0.61538
## Specificity                    0.9585     0.9416        0.99613
## Pos Pred Value                 0.8244     0.8571        0.80000
## Neg Pred Value                 0.9731     0.9745        0.99038
## Precision                      0.8244     0.8571        0.80000
## Recall                         0.8802     0.9343        0.61538
## F1                             0.8514     0.8940        0.69565
## Prevalence                     0.1811     0.2726        0.02453
## Detection Rate                 0.1594     0.2547        0.01509
## Detection Prevalence           0.1934     0.2972        0.01887
## Balanced Accuracy              0.9194     0.9379        0.80576
##                      Class: ActVig Class: Shiftpattern Class: AD Class: DE
## Sensitivity                0.80000             0.44444    0.8494    0.8968
## Specificity                0.99608             0.99902    0.9855    0.9839
## Pos Pred Value             0.88889             0.94118    0.9156    0.8828
## Neg Pred Value             0.99219             0.98082    0.9724    0.9861
## Precision                  0.88889             0.94118    0.9156    0.8828
## Recall                     0.80000             0.44444    0.8494    0.8968
## F1                         0.84211             0.60377    0.8812    0.8898
## Prevalence                 0.03774             0.03396    0.1566    0.1189
## Detection Rate             0.03019             0.01509    0.1330    0.1066
## Detection Prevalence       0.03396             0.01604    0.1453    0.1208
## Balanced Accuracy          0.89804             0.72173    0.9174    0.9404
##                      Class: LD Class: Flatsinusoidal Class: Suspect
## Sensitivity            0.94340               0.85294        0.87755
## Specificity            0.99603               0.99708        0.98649
## Pos Pred Value         0.92593               0.90625        0.86869
## Neg Pred Value         0.99702               0.99514        0.98751
## Precision              0.92593               0.90625        0.86869
## Recall                 0.94340               0.85294        0.87755
## F1                     0.93458               0.87879        0.87310
## Prevalence             0.05000               0.03208        0.09245
## Detection Rate         0.04717               0.02736        0.08113
## Detection Prevalence   0.05094               0.03019        0.09340
## Balanced Accuracy      0.96971               0.92501        0.93202

plot(rfmod, 
     main = "RF Learning curve")

p2=multiclasscmplot(predict(rfmod,newdata=testset),testset$CLASS)+ggtitle("RF")

p2

Kết quả cho thấy mô hình Random Forest đã chiến thắng trong cuộc thi đấu này. Tuy nhiên, nguyên nhân của sự khác biệt này có thể vì cấu trúc của mạng neuronfeed forward này còn quá đơn giản và chưa thực sự tối ưu. Hy vọng các bạn có thể cải tiến nó trong quá trình tự học keras của mình.

gridExtra::grid.arrange(p1,p2,ncol=2)

5 Nhận xét

Nếu bạn muốn dùng Deep learning để giải quyết vấn đề của mình,nên sử dụng package kerasvì nó có rất nhiều ưu điểm so với h2o, như quy trình có tính thứ bậc và cấu trúc mạng neuron được phác thảo dưới bàn tay của bạn một cách tường minh, tiến trình huấn luyện được báo cáo ở thời gian thực và kết quả có thể khai thác bằng những thủ thuật quen thuộc trong R, thí dụ confussion matrix của caret. Tuy nhiên giai đoạn chuẩn bị dữ liệu thực sự là một thử thách với các bạn chưa quen.

Chúc các bạn thành công (Nếu Nhi làm được thì các bạn bác sĩ khác cũng hoàn toàn có thể :))