Deep learning: Bài toán nhị phân

1 Đặt vấn đề

Trong năm 2017, sự ra đời của package keras đã cho phép người dùng R có thể thực hành được một phương pháp mô hình thuộc loại tân tiến nhất, đó là mạng thần kinh nhân tạo sâu (Deep neural network hay Deep learning) với TensofFlow (trước đây bạn phải dùng Python). Package keras cho phép R kết nối với TensorFlow back-end của Python để dựng mô hình Deep ANN với nhiều ứng dụng hữu ích như hồi quy, xếp loại nhị phân hay đa giá trị, computer vision (convolutional Neural network), dự báo chuỗi thời gian, mô hình tự sinh mẫu, vv

Ở bài trước, Nhi đã trình bày 1 ứng dụng Hồi quy điển hình với keras. Trong bài này, Nhi sẽ tiếp tục giới thiệu một mô hình phân loại với kết quả nhị phân (Binary Classification).Trong bài, Nhi sẽ thực hiện tuần tự từng bước của một quy trình mẫu mực cho nghiên cứu ứng dụng Machine learning. Khi thực sự bắt tay vào làm nghiên cứu bạn cũng sẽ phải đi qua những bước này, bao gồm:

• Biện luận lựa chọn Algorithm

• Thăm dò dữ liệu

• Hoán chuyển, sơ chế dữ liệu

• Huấn luyện mô hình

• Kiểm định mô hình

• Biện luận, giải nghĩa mô hình

Một mục tiêu khác của bài này, đó là giúp cho các bạn làm quen với một phong cách viết R code cũng như quy trình phân tích dữ liệu hoàn toàn mới lạ bằng cách ứng dụng những package mới nhất của 2 tác giả Hadley Wickham và Max Kuhn, như recipes (sơ chế dữ liệu), rsample (tái chọn mẫu), dplyr (thao tác dữ liệu và tóm tắt), và các toán tử pipe trong tidyverse.

Trong bài này, Nhi sử dụng bộ số liệu về bệnh Suy giáp của tác giả J Ross Quinlan và viện Garvan (Úc) (1987). Dữ liệu này gồm hơn 3700 trường hợp với mục tiêu nghiên cứu (giả định) là xây dựng một mô hình chẩn đoán bệnh Suy Giáp (Hypothyroid) dựa vào thông tin gồm 6 biến số liên tục : Tuổi và giá trị các biomarker như hormone T3, TT4, T4U và FTI, và 12 biến nhị phân gồm Giới tính, điều trị thyroxine, có thai, phẫu thuật tuyến giáp, suy tuyến não thùy (hypopituitary ), triệu chứng bứu (goitre,tumor),tâm lý (psych), điều trị I131, lithium. Phân loại bệnh nhược giáp trong dữ liệu nguyên thủy có đến 5 nhãn kết quả là Negative, hypothyroid, primary hypothyroid, compensated hypothyroid và secondary hypothyroid. Để đơn giản hóa, ta sử dụng phiên bản giản lược của Quan Sun trên thư viện OpenML trong đó biến kết quả được giản lược chỉ còn 2 nhãn (nhị phân) là Positive và Negative.

2 Bước 1: Thăm dò dữ liệu

library(tidyverse)

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

## v ggplot2 2.2.1     v purrr   0.2.4
## v tibble  1.4.2     v dplyr   0.7.5
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0

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

df=read.csv("https://www.openml.org/data/get_csv/53534/hypothyroid.csv",na.strings = "?")

df=df%>%dplyr::select(age,sex,pregnant,
                      on.thyroxine,query.on.thyroxine,
                      on.antithyroid.medication,
                      thyroid.surgery,
                      I131.treatment,sick,
                      lithium,goitre,tumor,hypopituitary,psych,
                      TSH,T3,TT4,T4U,FTI,binaryClass)

df=df%>%filter(.,age<100)%>%na.omit()

colnames(df)=c("AGE","SEX","PREG","ON_THY","QO_THY","ON_AT","SURG",
               "I131","SICK","LIT","GOIT","TUMOR","HPI","PSY",
               "TSH","T3","TT4","T4U","FTI","CLASS")

str(df)

## 'data.frame':    2642 obs. of  20 variables:
##  $ AGE   : int  41 70 80 66 68 84 71 59 28 63 ...
##  $ SEX   : Factor w/ 2 levels "F","M": 1 1 1 1 2 1 1 1 2 1 ...
##  $ PREG  : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ON_THY: Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ QO_THY: Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ON_AT : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ SURG  : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ I131  : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ SICK  : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 2 1 1 1 ...
##  $ LIT   : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ GOIT  : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ TUMOR : Factor w/ 2 levels "f","t": 1 1 1 2 1 2 1 1 1 1 ...
##  $ HPI   : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ PSY   : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ TSH   : num  1.3 0.72 2.2 0.6 2.4 1.1 0.03 2.8 3.3 1.5 ...
##  $ T3    : num  2.5 1.2 0.6 2.2 1.6 2.2 3.8 1.7 1.8 1.2 ...
##  $ TT4   : num  125 61 80 123 83 115 171 97 109 117 ...
##  $ T4U   : num  1.14 0.87 0.7 0.93 0.89 0.95 1.13 0.91 0.91 0.96 ...
##  $ FTI   : num  109 70 115 132 93 121 151 107 119 121 ...
##  $ CLASS : Factor w/ 2 levels "N","P": 2 2 2 2 2 2 2 2 2 2 ...
##  - attr(*, "na.action")= 'omit' Named int  2 3 4 6 7 12 16 17 21 24 ...
##   ..- attr(*, "names")= chr  "2" "3" "4" "6" ...

Như vậy, cấu trúc dữ liệu gồm có: 1 biến kết quả nhị phân là CLASS, 13 features là biến định tính nhị phân, 6 biến định lượng liên tục. Dữ liệu gồm 2642 trường hợp

table(df$CLASS)

## 
##    N    P 
##  216 2426

Tổng cộng có 2426 trường hợp Positive và 216 trường hợp Negative. Như vậy có sự bất xứng về nhãn kết quả, gợi ý rằng trong quá trình tái chọn mẫu và phân chia, ta phải cố gắng bảo toàn tỉ lệ giữa 2 nhãn P và N.

Điều đầu tiên Nhi sẽ làm, đó là cắt ngẫu nhiên dữ liệu gốc thành 2 phần bằng nhau (do ta có quá nhiều dữ liệu, nên tỉ lệ 50/50 là lý tưởng). Có thể làm việc này chỉ bằng 3 dòng code với hàm initial_split của gói rsample.

Với tùy chỉnh strata cho biến CLASS, ta có thể bảo toàn tỉ lệ giữa 2 nhãn P/N như nhau cho 2 tập trainset và testset là 1213/108

dat=df

# Prior Spliting 

library(rsample)

## Loading required package: broom

## 
## Attaching package: 'rsample'

## The following object is masked from 'package:tidyr':
## 
##     fill

set.seed(2405)

data_split <- initial_split(dat, strata = "CLASS",prop=0.5)

trainset <- training(data_split)
testset <- testing(data_split)

table(trainset$CLASS)

## 
##    N    P 
##  108 1213

table(testset$CLASS)

## 
##    N    P 
##  108 1213

Ta sẽ không nhìn vào tập testset, vì đây là dữ liệu dùng để kiểm định mô hình, nhưng ta có thể khảo sát tập trainset. Nhi sẽ dùng phương pháp trực quan để đánh giá phân bố của các biến và so sánh giữa 2 phân nhóm P/N

# First Exploration

trainset%>%gather(AGE,TSH,T3,TT4,T4U,FTI,key="Feature",value="Score")%>%
  ggplot(aes(x=Score,fill=CLASS))+
  geom_density(alpha=0.5,col="black")+
  theme_bw()+
  facet_wrap(~Feature,scales="free",ncol=3)+
  scale_fill_manual(values=c("blue","red"))

trainset%>%gather(PREG:PSY,key="Feature",value="Status")%>%
  ggplot(aes(x=Status,fill=CLASS))+
  geom_bar(stat="count",position="fill",alpha=0.5,col="black")+
  scale_y_continuous(labels=NULL)+
  theme_bw()+
  facet_wrap(~Feature,scales="free",ncol=6)+
  scale_fill_manual(values=c("blue","red"))

Ở hình thứ nhất, ta có thể nhận định rằng hầu hết biến định lượng có phân bố bất thường, không đồng dạng và không cân đối. Điều này gợi ý rằng một phép hoán chuyển để hiệu chỉnh Skewness (thí dụ Box-Cox) cần được áp dụng. Ngoài ra, do ta sẽ dùng ANN, và giải thuật này chỉ hoạt động tối ưu khi dữ liệu đầu vào được chuẩn hóa về cùng thang đo, nên ta sẽ áp dụng Centering và Scale cho toàn bộ biến định lượng.

Hình thứ hai cho thấy sự tương phản của các biến định tính trong dữ liệu là rất yếu, thậm chí gợi ý rằng ta phải loại bỏ một số biến không cho thấy sự tương phản rõ như I131, SICK.

3 Bước 2: Chuẩn bị dữ liệu

Bước tiếp theo sẽ là hoán chuyển dữ liệu, cũng như quá trình chuẩn bị nguyên liệu phù hợp với món ăn mà ta muốn nấu. Để làm việc này, ta sẽ dùng package recipes để tạo ra 1 công thức nấu ăn có nội dung như sau:

1. Đầu tiên, lấy toàn bộ features ,

2. Sau đó, làm hoán chuyển Box-Cox cho toàn bộ những biến định lượng

3. Tiếp theo, chuẩn hóa toàn bộ biến định lượng

4. Sau đó loại bỏ những biến định lượng tương quan mạnh với nhau, và chỉ giữ lại những biến không có tương quan

5. Tiếp theo, chuyển những biến định tính thành dummy variable có giá trị 0/1

6. Cuối cùng, loại bỏ tất cả những biến nhị phân lẫn định lượng nào có variance rất thấp (=0 hay gần =0).

Quy trình liên hoàn trên đây có vẻ phức tạp, nhưng khi dùng recipes mọi chuyện hết sức dễ dàng, bạn chỉ việc đặt bước này sau bước kia và xâu chuỗi bằng pipes, tận cùng bằng hàm prep.

Khi áp dụng công thức này trên trainset, ta sẽ thu được tập dữ liệu dùng để dựng mô hình,

Tương tự, khi áp dụng cho testset, ta sẽ thu được dữ liệu dùng để kiểm định mô hình

# recipes

library(recipes)

## 
## Attaching package: 'recipes'

## The following object is masked from 'package:stringr':
## 
##     fixed

## The following object is masked from 'package:stats':
## 
##     step

rec_obj<-recipe(CLASS~., data = trainset)%>%
  step_BoxCox(all_numeric(), -all_outcomes())%>%
  step_center(all_numeric(),-all_outcomes())%>%
  step_scale(all_numeric(),-all_outcomes())%>%
  step_corr(all_numeric(),-all_outcomes())%>%
  step_dummy(all_nominal(),-all_outcomes())%>%
  step_zv(all_predictors(),-all_outcomes())%>%
  step_nzv(all_predictors(),-all_outcomes())%>%
  prep(training = trainset)

train_df <- bake(rec_obj, newdata = trainset)
test_df <- bake(rec_obj, newdata = testset)

str(train_df)

## Classes 'tbl_df', 'tbl' and 'data.frame':    1321 obs. of  10 variables:
##  $ AGE     : num  0.9242 1.4924 0.9806 0.3102 0.0902 ...
##  $ TSH     : num  -0.302 0.366 -1.938 0.518 -1.006 ...
##  $ T3      : num  -1.007 -2.065 2.071 -0.287 -0.154 ...
##  $ TT4     : num  -1.375 -0.766 1.745 -0.255 0.783 ...
##  $ T4U     : num  -0.633 -1.903 0.77 -0.383 0.236 ...
##  $ FTI     : num  -1.2407 0.2277 1.3602 -0.0284 0.7348 ...
##  $ CLASS   : Factor w/ 2 levels "N","P": 2 2 2 2 2 2 2 2 2 2 ...
##  $ SEX_M   : num  0 0 0 0 1 0 0 0 0 0 ...
##  $ ON_THY_t: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ PSY_t   : num  0 0 0 0 0 0 0 0 0 0 ...

Như vậy, sau khi sơ chế ta đã loại bỏ được hầu hết biến định tính có phương sai quá thấp và chuẩn hóa hoàn toàn biến định lượng

train_df%>%
  gather(AGE:FTI,key="Feature",value="Score")%>%
  ggplot(aes(x=Score,fill=CLASS))+
  geom_density(alpha=0.5,col="black")+
  theme_bw()+
  facet_wrap(~Feature,scales="free",ncol=3)+
  scale_fill_manual(values=c("blue","red"))

train_df%>%
  gather(SEX_M:PSY_t,key="Feature",value="Status")%>%
  ggplot(aes(x=factor(Status),fill=CLASS))+
  geom_bar(stat="count",position="fill",alpha=0.5,col="black")+
  scale_y_continuous(labels=NULL)+
  theme_bw()+
  facet_wrap(~Feature,scales="free",ncol=6)+
  scale_fill_manual(values=c("blue","red"))

Ta hy vọng với dữ liệu này, mô hình sẽ hoạt động tối ưu :

4 Bước 3: Xây dựng cấu trúc mạng neuron

Trước khi lựa chọn một giải thuật (algorithm) Machine learning cho nghiên cứu của mình, ta cần biện luận về cơ sở của lựa chọn này. Trong thí dụ này, Nhi phải giải thích tại sao lại dùng một mô hình blackbox thay vì hồi quy logistic ?

Mô hình Neural network phù hợp trong 2 trường hợp :

• Một bệnh lý mà triệu chứng lâm sàng, rối loạn sinh lý, sinh hóa đã được mô tả một cách kinh điển, đầy đủ trong y văn và không còn gì phải bàn cãi nữa, tuy nhiên mục tiêu của nghiên cứu không phải là diễn dịch/phân tích, nhưng nhằm thiết lập một công cụ tự động (software) hỗ trợ bác sĩ đưa ra quyết định chẩn đoán trong những trường hợp khó.

• Dữ liệu đầu vào có bản chất phức tạp, vượt ra khỏi quy tắc thông thường của môn Thống kê, thí dụ hình ảnh (ảnh chụp nội soi, giải phẫu bệnh, MRI, …), âm thanh (tiếng tim, âm hô hấp), tín hiệu điện (ECG, EEG), văn bản (chữ viết tay, bệnh án điện tử)…, vì Deep learning cho phép phân tích rất tốt những loại dữ liệu này.

Thí dụ hiện thời thuộc trường hợp thứ nhất: Mục tiêu giả định là thiết lập một quy luật chẩn đoán tự động, càng chính xác càng tốt nhưng không quan tâm đến liên hệ giữa các biến trong dữ liệu.

Trước hết, ta trích xuất 2 vectors của biến CLASS sau khi chuyển nó thành biến nhị phân (0/1, cho tập train_df và test_df)

train_label <- ifelse(pull(train_df, CLASS) == "P", 1, 0)

test_label <- ifelse(pull(test_df, CLASS) == "P", 1, 0)

Tiếp theo, ta dựng cấu trúc (architecture) của mạng neuron:

Đây là một mạng ANN gồm 3 lớp, 1 lớp đầu vào gồm 9 perceptron để tiếp nhận 9 features. Tín hiệu trọng số sẽ được kích hoạt bằng hàm relu , sau đó đi qua 3 lớp ẩn, gồm 1 lớp ẩn thứ nhất có 64 neuron,liên kết bão hoàn với đầu vào, sau đó tín hiệu được lọc lần thứ nhất bằng 1 lớp drop_out với tỉ lệ 0.01, tiếp tục đi qua lớp ẩn thứ 2 có 32 neuron, rồi 1 lớp dropout thứ hai trước khi kết thúc ở lớp cuối cùng chỉ có 1 neuron để xuất kết quả, kích hoạt bằng hàm sigmoid.

Mạng ANN này sau đó được compile với tiêu chí huấn luyện là binary_accuracy, hàm loss là binary_crossentropy và optimizer là hàm “adam”

5 Bước 5: Huấn luyện mô hình

Ta chạy thử cấu trúc này với 50 lượt lấy mẫu với kích thước 100 trường hợp và kiểm định trên 30%

# ANN

library(keras)

build_model <- function(){
  model_keras <- keras_model_sequential()%>% 
    # 1st hidden layer
    layer_dense(
      units              = 64, 
      kernel_initializer = "uniform", 
      activation         = "relu", 
      input_shape        = ncol(train_df)-1) %>% 
    # 1st Dropout 
    layer_dropout(rate = 0.01) %>%
    # 2nd hidden layer
    layer_dense(
      units              = 32, 
      kernel_initializer = "uniform", 
      activation         = "relu") %>% 
    # Dropout 
    layer_dropout(rate = 0.01) %>%
    #  Output layer
    layer_dense(
      units              = 1, 
      kernel_initializer = "uniform", 
      activation         = "sigmoid") %>% 
    # Compile ANN
    compile(
      optimizer = 'adam',
      loss      = 'binary_crossentropy',
      metrics   = c('binary_accuracy')
    )
}

model<-build_model()

fit_keras <- fit(
  object=model,
  x=as.matrix(train_df[,-7]), 
  y=train_label,
  batch_size       = 100,
  epochs           = 50, 
  validation_split = 0.3,
  verbose = 0
)

fit_keras

## Trained on 924 samples, validated on 397 samples (batch_size=100, epochs=50)
## Final epoch (plot to see history):
##            val_loss: 0.04597
## val_binary_accuracy: 0.9849
##                loss: 0.04735
##     binary_accuracy: 0.9848

Kết quả huấn luyện rất khả quan, accuracy đã được tối ưu đến 98.59% và loss hạ xuống còn 0.04

plot(fit_keras)+theme_bw()

6 Bước 4: Kiểm chứng chéo

Tuy nhiên, phẩm chất của mô hình còn tùy thuộc vào nhiều tham số trong quá trình huấn luyện, thí dụ đơn giản nhất ở đây là số lượt (epoch), Nhi muốn biết liệu tăng epoch cao hơn có làm thay đổi hiệu quả của mô hình ANN hay không ?

Để có câu trả lời, Nhi sẽ thực hiện một quy trình kiểm chứng chéo 10 blocks, mỗi 9 block như vậysẽ dùng để huấn luyện, Nhi sẽ thử lần lượt giá trị epoch tăng dần từ 20,40,60,80,100,120 và kiểm định hiệu năng mô hình trên 1 block còn lại

Nhi viết 1 hàm tương đối “khủng” với nội dung như sau:

1. Khai báo giá trị epoch, danh sách tái chọn mẫu split, và seed number.

2. Sau mỗi lần chạy, xóa bộ nhớ của keras để chuẩn bị cho lần chạy tiếp theo

3. Dựng 1 mạng neuron có cấu trúc như trên

4. Từ danh sách chọn mẫu, trích xuất 2 matrices tập train, test và vector label

5. Huấn luyện ANN sử dụng 9 blocks

6. Phân loại 1 block (test) dưới dạng probability (0:1) và dán nhãn P/N

7. Tính tất cả (xin nhấn mạnh: tất cả) tiêu chí phẩm chất của model vừa dựng, bao gồm 4 nhóm: tính hữu dụng, tính chính xác, nguy cơ sai lầm, sự tương hợp

(Xin xem bài này để hiểu rõ hơn: http://rpubs.com/lengockhanhi/347941 )

8. Xuất kết quả kiểm định

library(yardstick)

## 
## Attaching package: 'yardstick'

## The following object is masked from 'package:readr':
## 
##     spec

library(scoring)
library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(rsample)

# A large Cross validation function


large_cv_func<- function(epoch, split, seed, ...) {
  # Set the seed to get reproducible starting values and dropouts
  set.seed(seed)
  # clears memory used by the last trial  
  on.exit(keras::backend()$clear_session())

  # Define ANN
  
  model<- keras_model_sequential()%>% 
    # 1st hidden layer
    layer_dense(
      units              = 64, 
      kernel_initializer = "uniform", 
      activation         = "relu", 
      input_shape        = ncol(train_df)-1) %>% 
    # 1st Dropout 
    layer_dropout(rate = 0.01) %>%
    # 2nd hidden layer
    layer_dense(
      units              = 32, 
      kernel_initializer = "uniform", 
      activation         = "relu") %>% 
    # Dropout 
    layer_dropout(rate = 0.01) %>%
    #  Output layer
    layer_dense(
      units              = 1, 
      kernel_initializer = "uniform", 
      activation         = "sigmoid") %>% 
    # Compile ANN
    compile(
      optimizer = 'adam',
      loss      = 'binary_crossentropy',
      metrics   = c('binary_accuracy')
    )  
  # data used for training ("analysis" set)
  block_df<- analysis(split)

  block_lab<-ifelse(pull(block_df,CLASS)=="P",1,0)
  
  model%>%fit(
    x=as.matrix(block_df[,-7]),
    y=block_lab,
    epochs = epoch, 
    ...
  )
  
  # Now obtain the holdout set for prediction
  blocktest_df<- assessment(split)
  
  blocktest_lab<-ifelse(pull(blocktest_df,CLASS)=="P",1,0)
  
  
  Pred_class_vector<- 
    predict_classes(object = model, 
                    x = as.matrix(blocktest_df[,-7]))%>%as.vector()
  
  Pred_Prob_vector<- 
    predict_proba(object = model, 
                  x = as.matrix(blocktest_df[,-7]))%>%as.vector()
  
  valid_df <- tibble(
    truth      = as.factor(blocktest_lab) %>% 
      fct_recode (Positive = "1", Negative = "0"),
    estimate   = as.factor(Pred_class_vector) %>% 
      fct_recode (Positive = "1", Negative = "0"),
    class_prob = Pred_Prob_vector )
  
  dfpred<-valid_df%>%mutate(Truth=.$truth,
                    BinTruth=if_else(.$truth=="P",1,0),
                    ClassNeg=1-.$class_prob,
                    ClassPos=.$class_prob,
                    ClassLab=.$estimate)
  
  TP=with(dfpred,sum(ClassLab==Truth & ClassLab=="Positive"))
  TN=with(dfpred,sum(ClassLab==Truth & ClassLab!="Positive"))
  FN=with(dfpred,sum(ClassLab!=Truth & Truth=="Positive"))
  FP=with(dfpred,sum(ClassLab!=Truth & Truth!="Positive"))
  
  TPR=TP/(TP+FN)
  TNR=TN/(TN+FP)
  SEN=TPR
  SPEC=TNR
  Recall=SEN
  
  ACC=(TP+TN)/(TP+FP+TN+FN)
  BAC=(TPR+TNR)/2
  PPV=TP/(TP+FP)
  NPV=TN/(TN+FN)
  Precision=PPV
  G=sqrt(Recall*Precision)
  F2=(1+2^2)*(Precision*Recall)/((2^2)*Precision+Recall)
  F0.5=(1+(0.5)^2)*(Precision*Recall)/((0.5^2)*Precision+Recall)
  F1=2*TP/(TP+TP+FP+FN)
  Fscore=(1+1)*(Precision*Recall)/(1*Precision+Recall)
  
  FNR=FN/(FN+TP)
  FPR=FP/(FP+TN)
  FDR=FP/(TP+FP)
  FOR=1-NPV
  BER=with(dfpred,mean((ClassPos>0.5 & BinTruth==0) | (ClassPos<0.5 & BinTruth==1)))
  logloss=with(dfpred,-mean(log10(ClassPos)))
  LSR=with(dfpred,mean(log10(ClassPos)))
  mmce=with(dfpred,mean(Truth!=ClassLab))
  
  brier= dfpred%>%scoring::brierscore(data=.,BinTruth~ClassPos)%>%mean()
  
  SSR=dfpred%>%
    scoring::sphscore(data=.,BinTruth~ClassPos)%>%mean()
  

  BM=TPR+TNR-1
  MK=PPV+NPV-1
  MCC=(TP*TN-FP*FN)/(sqrt((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)))
  Kappa=psych::cohen.kappa(as.matrix(dfpred[,c(1,2)]))%>%.$weighted.kappa
  
  rbind(BM,MK,MCC,Kappa)
  
  LRpos=SEN/(1-SPEC)
  LRneg=(1-SEN)/SPEC
  
  
  myroc <- roc(Truth~ClassPos, data=dfpred)
  
  AUC=myroc$auc[1]
  
  return(cbind(
    TP,TN,FN,FP,
    TPR,TNR,SEN,SPEC,Recall,
    ACC,BAC,PPV,NPV,Precision,G,F1,F2,F0.5,Fscore,
    FNR,FPR,FDR,FOR,BER,logloss,LSR,mmce,brier,SSR,
    BM,MK,MCC,Kappa,LRpos,LRneg,AUC)
    )
}

Ta test thử hàm này với 1 block CV và 30 lượt (epoch=30)

# V-folds Cross-validation

set.seed(2405)
cv_splits <- vfold_cv(data=train_df,v = 10,strata ="CLASS")

cvdf=large_cv_func(epoch = 30,
              cv_splits$splits[[1]],
              seed=2405,
              batch_size = 100, 
              verbose = 0)

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

Hàm chạy tốt,

cvdf

##       TP TN FN FP TPR       TNR SEN      SPEC Recall       ACC       BAC
## [1,] 122  9  0  2   1 0.8181818   1 0.8181818      1 0.9849624 0.9090909
##           PPV NPV Precision         G        F1       F2     F0.5
## [1,] 0.983871   1  0.983871 0.9919027 0.9918699 0.996732 0.987055
##         Fscore FNR       FPR        FDR FOR       BER    logloss
## [1,] 0.9918699   0 0.1818182 0.01612903   0 0.9323308 0.07647281
##              LSR       mmce     brier       SSR        BM       MK
## [1,] -0.07647281 0.01503759 0.9079204 0.9147508 0.8181818 0.983871
##            MCC     Kappa LRpos LRneg AUC
## [1,] 0.8972098 0.8919578   5.5     0   1

Ta bắt đầu chạy thử quy trình kiểm chứng chéo 10 blocks này (sẽ hơi mất thời gian, bạn có thể đi ăn/uống cái gì đó)

epoch_values <-seq(20,120,by=20)

out=tibble(epoch=NA,fold=NA)

out=cbind(out,cvdf)

for(k in c(1:6)){
  for(i in c(1:10)){
    temp_id=tibble(epoch=epoch_values[k],fold=i)
    temp_out=large_cv_func(epoch=epoch_values[k],
                      cv_splits$splits[[i]],
                      seed=2405,
                      batch_size = 100, 
                      verbose = 0)
    temp=cbind(temp_id,temp_out)
    out=rbind(out,temp)
  }
}

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

mean_out<-out%>%.[-1,]%>%
  group_by(epoch)%>%
  summarize_if(is.numeric,mean)

mean_out%>%gather(TPR,TNR,SEN,SPEC,Recall,AUC,LRpos,LRneg,key="Utility",value="Score")%>%
ggplot(.,aes(x = epoch, y = Score,col=Utility)) + 
  geom_point()+ 
  geom_path()+ 
  theme_bw(8)+
  facet_wrap(~Utility,scales="free")+
  ggtitle("Utility")

mean_out%>%gather(ACC,BAC,PPV,NPV,Precision,G,F1,F2,F0.5,Fscore,
                  key="Accuracy",value="Score")%>%
  ggplot(.,aes(x = epoch, y = Score,col=Accuracy)) + 
  geom_point()+ 
  geom_path()+ 
  theme_bw(8)+
  facet_wrap(~Accuracy,scales="free",ncol=3)+
  ggtitle("Accuracy")

mean_out%>%gather(FNR,FPR,FDR,FOR,BER,logloss,LSR,mmce,brier,SSR,
                  key="Error",value="Score")%>%
  ggplot(.,aes(x = epoch, y = Score,col=Error)) + 
  geom_point()+ 
  geom_path()+ 
  theme_bw(8)+
  facet_wrap(~Error,scales="free",ncol=3)+
  ggtitle("Error")

mean_out%>%gather(BM,MK,MCC,Kappa,
                  key="Agreement",value="Score")%>%
  ggplot(.,aes(x = epoch, y = Score,col=Agreement)) + 
  geom_point()+ 
  geom_path()+ 
  theme_bw(8)+
  facet_wrap(~Agreement,scales="free")+
  ggtitle("Agreement")

Kết quả của quy trình kiểm chứng chéo 10 block cho thấy giá trị epoch tối ưu nằm giữa 80 và 100, khi đó các tiêu chí phẩm chất mô hình bắt đầu ổn định.

Ta huấn luyện lại mô hình lần cuối với epoch=100

model<-build_model()

fit_keras <- fit(
  object=model,
  x=as.matrix(train_df[,-7]), 
  y=train_label,
  batch_size       = 100,
  epochs           = 100,
  verbose=0
)

fit_keras

## Trained on 1,321 samples (batch_size=100, epochs=100)
## Final epoch (plot to see history):
##            loss: 0.01747
## binary_accuracy: 0.9955

plot(fit_keras)+theme_bw()

7 Bước 6: Kiểm định độc lập

Bây giờ ta sẽ kiểm định mô hình cuối cùng trên tập test_df

Pred_class_vector<- 
  predict_classes(object = model, 
                  x = as.matrix(test_df[,-7]))%>%as.vector()

Pred_Prob_vector<- 
  predict_proba(object = model, 
                x = as.matrix(test_df[,-7]))%>%as.vector()

Đầu tiên, một confusion matrix cho thấy quy luật chẩn đoán dựa vào Deep neural network có độ chính xác rất cao:

valid_df <- tibble(
  truth      = as.factor(test_label) %>% 
    fct_recode (Positive = "1", Negative = "0"),
  estimate   = as.factor(Pred_class_vector) %>% 
    fct_recode (Positive = "1", Negative = "0"),
  class_prob = Pred_Prob_vector )

valid_df %>% conf_mat (truth, estimate)

##           Truth
## Prediction Negative Positive
##   Negative       96        8
##   Positive       12     1205

Tính hữu dụng của mô hình được đo bằng 6 tiêu chí, bao gồm AUC= 0.997, cả độ nhạy và độ đặc hiệu đều cao

valid_df%>%mutate(Truth=.$truth,
                  BinTruth=if_else(testset$CLASS=="P",1,0),
                  ClassNeg=1-.$class_prob,
                  ClassPos=.$class_prob,
                  ClassLab=.$estimate)->dfpred

TP=with(dfpred,sum(ClassLab==Truth & ClassLab=="Positive"))
TN=with(dfpred,sum(ClassLab==Truth & ClassLab!="Positive"))
FN=with(dfpred,sum(ClassLab!=Truth & Truth=="Positive"))
FP=with(dfpred,sum(ClassLab!=Truth & Truth!="Positive"))

TPR=TP/(TP+FN)
TNR=TN/(TN+FP)
SEN=TPR
SPEC=TNR
Recall=SEN

myroc <- roc(Truth~ClassPos, data=dfpred)

AUC=myroc$auc

data_frame(fpr=1-myroc$specificities,
           tpr=myroc$sensitivities)%>%
  ggplot(aes(x=fpr,ymin=0,ymax=tpr))+
  geom_polygon(aes(y=tpr),fill="red",alpha=0.3)+
  geom_path(aes(y=tpr),col="red3",size=1)+
  geom_abline(linetype='dashed')+
  theme_bw()+
  coord_equal()+
  labs(x="False positive rate",y="True Positive rate")+
  ggtitle(paste0("AUC=",round(AUC,3)))

rbind(TPR,TNR,SEN,SPEC,Recall)

##             [,1]
## TPR    0.9934048
## TNR    0.8888889
## SEN    0.9934048
## SPEC   0.8888889
## Recall 0.9934048

Tiếp theo, 10 tiêu chí đánh giá độ chính xác của mô hình.2 tiêu chí quan trọng nhất bao gồm BAC =0.99, F1 score rất cao ~ 0.99

ACC=(TP+TN)/(TP+FP+TN+FN)
BAC=(TPR+TNR)/2
PPV=TP/(TP+FP)
NPV=TN/(TN+FN)
Precision=PPV
G=sqrt(Recall*Precision)
F2=(1+2^2)*(Precision*Recall)/((2^2)*Precision+Recall)
F0.5=(1+(0.5)^2)*(Precision*Recall)/((0.5^2)*Precision+Recall)
F1=2*TP/(TP+TP+FP+FN)
Fscore=(1+1)*(Precision*Recall)/(1*Precision+Recall)

rbind(ACC,BAC,PPV,NPV,Precision,G,F1,F2,F0.5,Fscore)

##                [,1]
## ACC       0.9848600
## BAC       0.9411468
## PPV       0.9901397
## NPV       0.9230769
## Precision 0.9901397
## G         0.9917709
## F1        0.9917695
## F2        0.9927500
## F0.5      0.9907910
## Fscore    0.9917695

Nguy cơ sai lầm của mô hình rất thấp: logloss=0.155, brier score=0.01 và BER=0.01

FNR=FN/(FN+TP)
FPR=FP/(FP+TN)
FDR=FP/(TP+FP)
FOR=1-NPV
BER=with(dfpred,mean((ClassPos>0.5 & BinTruth==0) | (ClassPos<0.5 & BinTruth==1)))
logloss=with(dfpred,-mean(log10(ClassPos)))
LSR=with(dfpred,mean(log10(ClassPos)))
mmce=with(dfpred,mean(Truth!=ClassLab))

library(scoring)
brier= dfpred%>%scoring::brierscore(data=.,BinTruth~ClassPos)%>%mean()

SSR=dfpred%>%
  scoring::sphscore(data=.,BinTruth~ClassPos)%>%mean()

rbind(FNR,FPR,FDR,FOR,BER,logloss,LSR,mmce,brier,SSR)

##                 [,1]
## FNR      0.006595218
## FPR      0.111111111
## FDR      0.009860312
## FOR      0.076923077
## BER      0.015140045
## logloss  0.165678366
## LSR     -0.165678366
## mmce     0.015140045
## brier    0.011227514
## SSR      0.012143483

Cuối cùng, mô hình có tính tương hợp cao giữa kết quả phân loại và quan sát thực tế: Kappa=0.908

BM=TPR+TNR-1
MK=PPV+NPV-1
MCC=(TP*TN-FP*FN)/(sqrt((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)))

## Warning in (TP + FP) * (TP + FN) * (TN + FP) * (TN + FN): NA produit par
## débordement d'entier par le haut

Kappa=psych::cohen.kappa(as.matrix(dfpred[,c(1,2)]))%>%.$weighted.kappa

rbind(BM,MK,MCC,Kappa)

##            [,1]
## BM    0.8822937
## MK    0.9132166
## MCC          NA
## Kappa 0.8974331

8 Bước 7: Giải nghĩa mô hình bằng LIME

Tuy quy tắc chẩn đoán mà ta vừa tạo ra có độ chính xác rất cao, nhưng nó cũng là một dạng mô hình bất khả tri (blackbox). Như ta biết, có một package tên là lime cho phép giải thích cách thức hoạt động của mô hình blackbox, thông qua cơ chế có tên là : Local Interpretable Model-Agnostic Explanations, tạm dịch : « Phép diễn giải cục bộ cho mô hình bất khả tri ».

Phương pháp LIME được giới thiệu lần đầu tiên năm 2016 trong một bài báo của 3 data scientists tại ĐH Washinton là Marco Tulio Ribeiro, Sameer Singh và Carlos Guestrin. https://arxiv.org/abs/1602.04938

Sau khi ứng dụng trong Python năm 2016, một phiên bản ứng dụng cho R do Thomas Lin Pedersen và Michaël Benesty đồng tác giả vừa mới được công bố trên CRAN tháng 9 năm 2017. https://cran.r-project.org/web/packages/lime/index.html https://github.com/thomasp85/lime

Một tutorial ngắn được giới thiệu kèm theo package https://cran.r-project.org/web/packages/lime/vignettes/Understanding_lime.html

Cơ chế hoạt động của LIME

Lime diễn giải bất cứ mô hình classifier nào, bao gồm các mô hình « bất khả tri » hay hộp đen (blackbox) theo quy trình như sau :

Phương pháp LIME dựa trên một giả định cơ bản, đó là bất kể mô hình phức tạp đến đâu, thì tại một miền cục bộ trong không gian dữ liệu (không gian nhỏ xíu xung quanh một điểm), mô hình có thể được ước lượng xấp xỉ bằng quy luật tuyến tính. Bạn có thể hình dung trên một vài bệnh nhân có đặc tính input data xấp xỉ như nhau, thì kết quả output của model sẽ như nhau, đó là bằng chứng về sự tuyến tính cục bộ.

Đầu tiên, LIME sẽ lấy thông tin về đặc tính phân phối của feature (dữ liệu đầu vào) dựa vào training dataset và nội dung (feature nào được dùng) trong model. Thông tin này được lưu trữ trong một object gọi là explainer.

Việc diễn giải sẽ được áp dụng cho 1 trường hợp cá thể, mới (unseen case), thí dụ một bệnh nhân bất kì trích từ tập kiểm định (testing subset).

LIME sẽ mô phỏng một lượng lớn các trường hợp giả định (một đám mây nhiễu của features) nằm kề cận chung quanh trường hợp (điểm) đang được xét, dựa vào quy luật phân phối của features mà nó đã ghi nhận từ trước.

Sau đó LIME áp dụng mô hình cho toàn bộ những điểm trong không gian nhiễu này, đồng thời tính khoảng cách giữa các điểm mô phỏng đến điểm trung tâm là trường hợp được xét. Khoảng cách này sẽ được chuyển thành thang điểm (score).

Tiếp theo, LIME chọn một số lượng M features tiêu biểu nhất cho phép mô tả tốt nhất khoảng cách nói trên.

Cuối cùng, LIME dựng một mô hình rất đơn giản cho các điểm mô phỏng, sử dụng M features được chọn làm predictor, để giải nghĩa cho outcome của model. Mô hình này có dạng Tuyến tính, hoặc mô hình Decision tree. Tham số hồi quy cho mỗi Features được điều chỉnh bằng một trong số (Weight) tỉ lệ với khoảng cách sai biệt với giá trị feature có thực của cá thể.

Việc diễn giải tính hợp lý của kết quả được thực hiện dựa vào Weight coefficient và danh sách M features được chọn. Nếu Weight coefficient > 0, thì giá trị quan sát của feature Mi đang ủng hộ cho kết quả tiên lượng (outcome) P, ngược lại, Weight Coefficient <0 thì giá trị feature Mi chống lại kết quả P.

Thí dụ, ta thử giải thích kết quả của mô hình trên 3 bệnh nhân (nhãn P) và 3 người nhãn N:

library(lime)

## 
## Attaching package: 'lime'

## The following object is masked from 'package:dplyr':
## 
##     explain

model_type.keras.models.Sequential <- function(x, ...) {
  "classification"}

predict_model.keras.models.Sequential <- function (x, newdata, type, ...) {
  pred <- predict_proba (object = x, x = as.matrix(newdata))
  data.frame (Positive = pred, Negative = 1 - pred) }

predict_model(x= model, 
               newdata = test_df[,-7], 
               type    = 'raw')%>%as_tibble()

## # A tibble: 1,321 x 2
##    Positive    Negative
##       <dbl>       <dbl>
##  1   1.000  0.00000238 
##  2   1      0          
##  3   1.000  0.0000186  
##  4   1.000  0.000000119
##  5   1.000  0.000114   
##  6   1.000  0.00000620 
##  7   0.0250 0.975      
##  8   1.000  0.000000834
##  9   1      0          
## 10   1.000  0.00000989 
## # ... with 1,311 more rows

explainer <- lime (
  x              = train_df[,-7], 
  model          = model, 
  bin_continuous = FALSE)

test_pos=test_df%>%dplyr::filter(CLASS=="P")%>%.[1:3,-7]
test_neg=test_df%>%dplyr::filter(CLASS=="N")%>%.[1:3,-7]

testlime=rbind(test_pos,test_neg)

explanation <- lime::explain (
  testlime, # Just to show first 10 cases
  explainer = explainer, 
  n_labels = 1, # explaining a `single class`(Polarity)
  n_features = 5, # returns top four features critical to each case
  kernel_width = 0.5) 

plot_features (explanation)+theme_minimal(5)

Ta cũng có thể huy động hàng loạt trường hợp để tổng kết về vai trò của từng features một cách tổng quát trong mô hình ANN:

test_pos=test_df%>%dplyr::filter(CLASS=="P")%>%.[1:50,-7]
test_neg=test_df%>%dplyr::filter(CLASS=="N")%>%.[1:50,-7]

testlime=rbind(test_pos,test_neg)

explanation2 <- lime::explain (
  testlime,
  explainer = explainer, 
  n_labels = 1, # explaining a `single class`(Polarity)
  n_features = 5, # returns top four features critical to each case
  kernel_width = 0.5) 

plot_explanations(explanation2)+theme_minimal(8)

Kết quả cho thấy, 3 features có vai trò quan trọng nhất lần lượt là TSH, On_Thyroxine treatment và FTI, cho phép phân định nhãn P hoặc N trong đa số trường hợp. Cụ thể, TSH thấp và FTI cao sẽ làm tăng xác suất của nhãn P và ngược lại.

9 Kết luận

Bài thực hành đến đây là hết, Nhi đã chuyển đến các bạn một quy trình tương đối hoàn chỉnh và chuẩn mực cho một nghiên cứu với mục tiêu xây dựng quy luật chẩn đoán có sự hỗ trợ của máy tính.

Như đã bàn ở trên, mô hình Deep neural network giải quyết khá hiệu quả bài toán chẩn đoán bệnh lý theo kiểu nhị phân. Việc huấn luyện mô hình tương đối nhẹ nhàng với điều kiện dữ liệu được chuẩn bị phù hợp và một chút kinh nghiệm về cấu trúc layer và các hàm activation, hàm loss cũng như tiêu chí huấn luyện.

Ta hoàn toàn có thể tinh chỉnh các tham số trong cấu trúc ANN cũng như tham số huấn luyện như bacht size, epoch,… package rsample cho phép tạo ra các quy trình tái chọn mẫu như ý thích, tứ boostrap cho đến kiểm chứng chéo ngẫu nhiên Montecarlo, kiểm chứng chéo K block …

Chúc các bạn thành công.