Naive Bayes ứng dụng cho y học

1 Đặt vấn đề

Như đã trình bày trong trước, Nhi đang đặc biệt quan tâm đến việc vận dụng những phương pháp phân tích dữ liệu đặc biệt từ phái Machine learning vào nghiên cứu y học truyền thống, với hy vọng thực hiện được hai mục tiêu: 1) Giúp bác sĩ và sinh viên Y khoa thoát khỏi cái hộp của những quy trình cổ điển/lạc hậu, tăng khả năng giải quyết vấn đề và sáng tạo; và 2) Xóa bỏ dần khoảng cách và những dị biệt giữa 2 lĩnh vực: Khoa học máy tính (Machine learning) và Thống kê y học.

Trong bài này, Nhi muốn giới thiệu đến các bạn đồng nghiệp về một giải thuật (algorithm) cổ điển bên phái Machine learning, có tên là Naive Bayes, những ưu thế và hạn chế của nó khi áp dụng giải quyết vấn đề trong chuyên khoa Y.

2 Giới thiệu về Naive Bayes

Naive Bayes là một trong những giải thuật Machine learning cổ điển được đề cập nhiều trong thập niên 1950-1960 để giải quyết các bài toán phân loại văn bản. Phương pháp Naive Bayes có mối liên hệ mật thiết với ngành Thống kê vì cơ chế của nó dựa vào định lý Bayes và khi vận dụng vào Y học thì Naive Bayes cũng tương đồng với quy trình biện luận lâm sàng của người bác sĩ. Tuy cổ xưa và quá đơn giản, nhưng Naive Bayes vẫn còn chỗ đứng ở thời đại ngày nay; nó vẫn được nhắc đến trong mọi giáo trình về Machine learning bên cạnh những giải thuật phức tạp khác, điều này cho thấy Naive Bayes có hiệu quả thực sự, và đáng để ta tìm hiểu.

Naive Bayes là một giải thuật dựa vào lý thuyết xác suất điều kiện.Giả sử ta có một bài toán phân loại với kết quả là C gồm k nhãn giá trị. Mục tiêu của chúng ta là xếp một cá thể đặc trưng bởi vector dữ liệu X vào một phân lớp C gồm k loại. Điều chúng ta sẽ thực sự làm, đó là ước tính xác suất cho mỗi nhãn giá trị Ci, từ C1 đến Ck trong điều kiện X hiện có. Nhãn giá trị nào có xác suất cao nhất sẽ được chọn làm quyết định sau cùng. Theo định lý Bayes, ta có :

\[P(C_{i}|X) = \frac{P(C_{i})P(X|C_{i})}{P(X)}\]

Trong công thức này, P(Ci) được gọi là xác suất tiền nghiệm mà ta biết về Ci, trước khi tiếp cận dữ liệu X. Một thí dụ thường gặp về P(Ci) trong y học là tỉ suất mắc bệnh trong quần thể. Trong bài toán phân loại, P(Ci) được cung cấp từ chính tập dữ liệu ta dùng để huấn luyện mô hình, vì ta biết tỉ lệ phân bố của mỗi nhãn Ci trên toàn bộ mẫu.

Ở mẫu số, P(X) được hiểu như xác suất quan sát được (những) giá trị của vector dữ liệu X, trong toàn bộ khả năng có thể của chúng trên thực tế. Việc tính giá trị chính xác của P(X) gần như bất khả thi, nhưng may mắn thay, điều này không thực sự cần thiết, vì P(X) là mẫu số chung cho tất cả nhãn Ci. Do đó, hướng đi của chúng ta thay đổi một chút, đó là tối ưu hóa cho riêng tử số để đạt giá trị cực đại :

\[\arg max \frac{P(C_{i})P(X|C_{i})}{P(X)} = \arg max P(C_{i})P(X|C_{i})\]

Thành phần còn lại : p(X|Ci) là xác suất quan sát được giá trị dữ liệu X khi biết nhãn phân loại Ci. Việc tính xác suất điều kiện này cũng gần như bất khả thi trên thực tế, vì X là một không gian dữ liệu đa chiều gồm giá trị các biến ngẫu nhiên và tất cả những khả năng tổ hợp có thể trong số này một khi có liên hệ, tương tác giữa chúng.

Để giải quyết khó khăn này, người ta áp dụng một giả định quan trọng, đó là các yếu tố (biến số) xj trong dữ liệu X hoàn toàn độc lập với nhau. Lý thuyết xác suất cho ta biết rằng khi A và B là hai sự kiện độc lập với nhau, ta có xác suất cho trường hợp A và B đồng thời xảy ra bằng tích của xác suất riêng của A và B :

\[P(A\cap B) = p(A)p(B)\]

Như vậy với giả định các biến xj trong dữ liệu X độc lập với nhau, ta có thể ước tính P(X|Ci) như sau:

\[P(X|Ci)=\prod_{j=1}^{n}p(xj|Ci)\]

Giả định này rất phi lý và chắc chắn không thể tồn tại trong thế giới thực, do đó phương pháp này mới có tên gọi là “Naive Bayes”, tạm dịch “suy luận Bayes ngây thơ”.

Như vậy, quy trình phân loại vector dữ liệu đầu vào X chính là sự tối ưu hóa giá trị cực đại của tích các xác suất điều kiện cho từng biến riêng lẻ. Cuối cùng, ta chuyển từ tích sang tổng bằng hàm logarit :

\[\arg max P(C_{i}) \prod_{j=1}^{n} p(x_{j}|C_{i}) \propto \arg max \left ( log(p(C_{i}) +\sum_{j=1}^{n} log(p(x_{j}|C_{i}) \right )\]

Bây giờ chỉ còn việc tính xác suất điều kiện đơn biến độc lập p(xj|Ci). Tùy vào bản chất của yếu tố/đại lượng được khảo sát, ta có thể áp dụng quy luật phân bố Gaussian cho biến số liên tục, Bernoulli cho biến nhị phân hay ước tính tần suất đơn giản (biến nhị phân hoặc rời rạc nhiều giá trị).

Nếu xj là một biến định danh (rời rạc), thí dụ Giới tính Nam/nữ, trạng thái Có/không của triệu chứng, khi đó xác suất cần tìm chính là tỉ lệ trường hợp có giá trị Cj trên tổng số trường hợp.

\[p(x_{j}|C)=\frac{N_{Cj}}{N_{C}}\]

Trong trường hợp phân bố của biến rời rạc xj bị mất cân bằng (một level của x không hiện diện trong một phân lớp Ci), ta có thể hiệu chỉnh bằng Laplace smoothing:

\[\frac{N_{Cj} +\alpha}{N_{C} + d\alpha}\]

Với alpha = 1 cho tử số và da cho mẫu số để đảm bảo tổng xác suất = 1.

Một số package còn áp dụng kỹ thuật Kernel cho biến liên tục và kết quả thường tốt hơn so với giả định máy móc về phân bố Gaussian.

Nếu xj là biến liên tục, hoặc ta dùng Gaussian NB, hoặc ta có thể hoán chuyển nó thành một biến rời rạc bằng cách cắt tại nhiều ngưỡng, thí dụ tứ phân vị, hoặc nhiều khoảng nhỏ hơn nữa.

\[p(x_{j}|C) = p(x_{j}|\mu_{cj},\sigma_{Cj}^{2})=\frac{1}{\sqrt{2\pi}\sigma_{Cj}^{2}}exp\left ( -\frac{(x_{j}-\mu_{Cj})^{2}}{2\sigma_{Cj}^{2}} \right)\]

3 Thí dụ minh họa:

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.

Lưu ý: Trong dữ liệu gốc có 1 trường hợp mà biến Age bị nhập liệu sai, với giá trị 445 tuổi,ta loại trường hợp này khỏi dữ liệu phân tích.

library(tidyverse)

df=read.csv("https://www.openml.org/data/get_csv/53534/hypothyroid.arff",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)

str(df)

## 'data.frame':    3770 obs. of  20 variables:
##  $ age                      : int  41 23 46 70 70 18 59 80 66 68 ...
##  $ sex                      : Factor w/ 2 levels "F","M": 1 1 2 1 1 1 1 1 1 2 ...
##  $ pregnant                 : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ on.thyroxine             : Factor w/ 2 levels "f","t": 1 1 1 2 1 2 1 1 1 1 ...
##  $ query.on.thyroxine       : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ on.antithyroid.medication: Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ thyroid.surgery          : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ I131.treatment           : 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 1 1 1 1 ...
##  $ lithium                  : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ goitre                   : 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 1 1 1 1 1 2 1 ...
##  $ hypopituitary            : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ psych                    : Factor w/ 2 levels "f","t": 1 1 1 1 1 1 1 1 1 1 ...
##  $ TSH                      : num  1.3 4.1 0.98 0.16 0.72 0.03 NA 2.2 0.6 2.4 ...
##  $ T3                       : num  2.5 2 NA 1.9 1.2 NA NA 0.6 2.2 1.6 ...
##  $ TT4                      : num  125 102 109 175 61 183 72 80 123 83 ...
##  $ T4U                      : num  1.14 NA 0.91 NA 0.87 1.3 0.92 0.7 0.93 0.89 ...
##  $ FTI                      : num  109 NA 120 NA 70 141 78 115 132 93 ...
##  $ binaryClass              : Factor w/ 2 levels "N","P": 2 2 2 2 2 2 2 2 2 2 ...

Đầu tiên, ta phân chia ngẫu nhiên dữ liệu gốc thành 2 phần, một dùng cho việc huấn luyện mô hình, một dùng để kiểm định mô hình này. Chúng ta không đụng đến testset cho đến khi kiểm định.

library(caret)

set.seed(1234)
idTrain=caret::createDataPartition(y=df$binaryClass, p=0.5,list=FALSE)

trainset=df[idTrain,]
testset=df[-idTrain,]

Tiếp theo ta sẽ thăm dò đặc tính phân bố của dữ liệu bằng biểu đồ. Đây là công đoạn đầu tiên trong mỗi nghiên cứu, nhưng đối với Naïve Bayes thì hình ảnh trực quan còn có ý nghĩa đặc biệt, vì những gì bạn nhìn thấy và cảm nhận lúc này cũng chính là những gì mà bạn sẽ thu được trong nội dung mô hình.

Do NB không xét liên hệ và tương tác giữa các biến, chúng ta chỉ quan tâm đến 2 điều: tính tương phản giữa các phân lớp và đặc tính phân bố của biến liên tục.

trainset[,-c(2:14)]%>%gather(age:FTI,key="Features",value="Score")%>%
  ggplot()+
  geom_density(aes(x=Score,fill=binaryClass),alpha=0.5)+
  facet_wrap(~Features,ncol=2,scales = "free")+
  theme_bw()+scale_fill_manual(values=c("blue","red"))

trainset[,-c(2:14)]%>%gather(age:FTI,key="Features",value="Score")%>%
  ggplot()+
  geom_violin(aes(x=binaryClass,y=Score,fill=binaryClass),alpha=0.5)+
  coord_flip()+
  facet_wrap(~Features,ncol=2,scales = "free")+
  theme_bw()+scale_fill_manual(values=c("blue","red"))

trainset[,c(2:14,20)]%>%gather(sex:psych,key="Features",value="Level")%>%na.omit()%>%
  ggplot()+
  geom_bar(aes(x=Level,y=..count..,fill=binaryClass),alpha=0.7)+
  coord_flip()+
  facet_wrap(~Features,ncol=3,scales = "free")+
  theme_bw(8)+scale_fill_manual(values=c("blue","red"))

plotfuncLow <- function(data,mapping){
  p <- ggplot(data = data,mapping=mapping)+
    stat_density2d(geom="polygon",aes(fill=trainset$binaryClass,alpha = ..level..))+
    geom_point(aes(fill=trainset$binaryClass,col=trainset$binaryClass),shape=".")+
    scale_fill_manual(values=c("blue","red"))+
    scale_color_manual(values=c("blue4","red4"))+
    theme_bw()
  p
}

plotdensfunc <- function(data,mapping){
  p <- ggplot(data = data,mapping=mapping)+
    stat_density(aes(fill=trainset$binaryClass,alpha = 0.7))+
    geom_rug(aes(col=trainset$binaryClass))+
    scale_fill_manual(values=c("blue","red"))+
    scale_color_manual(values=c("blue4","red4"))+
    theme_bw()
  p
}

library(GGally)

ggpairs(trainset[,-c(2:14,20)],lower = list(continuous=plotfuncLow),
        diag=list(continuous=plotdensfunc))

Dựa vào các biểu đồ, ta có thể phán đoán là: sự tương phản lớn nhất giữa 2 phân lớp P/N được biểu hiện tại các biến T4U, TT4,FTI, gợi ý đây là những biến số hữu hiệu trong mô hình.Ngoại trừ biến Age và TSH, các trị số hormone còn lại có hình dạng phân bố tương đối bình thường. Các biến nhị phân có giá trị False hầu hết đều có liên hệ với phân loại Positive của bệnh Nhược giáp. Do có sự mất cân bằng lớn trong các biến định tính, điển hình là biến hypopituitary thậm chí không có trường hợp nào thuộc loại T, do đó ta có thể phải dùng đến Laplace smoothing.

4 Package e1071

Mô hình Naive Bayes trong R có thể được dựng bằng 2 packages, đầu tiên là e1071.

hàm naiveBayes của package e1071 dựng mô hình Naive Bayes với tốc độ rất nhanh.Trong tùy chỉnh, ta khai báo na.omit để xử lý trường hợp thiếu sót dữ liệu riêng cho từng biến, và laplace=1.

Nội dung của mô hình cho thấy xác suất tiền nghiệm cho nhãn giá trị P là 0.915 và N là 0.085. Sau đó là một loạt các bảng chéo trình bày xác suất điều kiện từng giá trị của biến đầu vào trong điều kiện mỗi nhãn giá trị của biến đích nếu đầu vào là một biến rời rạc hay nhị phân (thí dụ Sex), hoặc Trung bình cho biến đầu vào liên tục.

Mô hình NBmod1 này rất nhạy, nhưng chưa đặc hiệu, nó vẫn còn chẩn đoán sai trong một số trường hợp Negative.

library(e1071)

NBmod1=e1071::naiveBayes(formula=binaryClass~.,
                         data=trainset,
                         na.action = na.omit,laplace=1)

NBmod1

## 
## Naive Bayes Classifier for Discrete Predictors
## 
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
## 
## A-priori probabilities:
## Y
##          N          P 
## 0.08530084 0.91469916 
## 
## Conditional probabilities:
##    age
## Y       [,1]     [,2]
##   N 52.66071 17.07564
##   P 53.60783 19.40314
## 
##    sex
## Y           F         M
##   N 0.7982456 0.2017544
##   P 0.6608479 0.3391521
## 
##    pregnant
## Y            f          t
##   N 0.99122807 0.00877193
##   P 0.98254364 0.01745636
## 
##    on.thyroxine
## Y            f          t
##   N 0.97368421 0.02631579
##   P 0.89359933 0.10640067
## 
##    query.on.thyroxine
## Y             f           t
##   N 0.991228070 0.008771930
##   P 0.991687448 0.008312552
## 
##    on.antithyroid.medication
## Y            f          t
##   N 0.99122807 0.00877193
##   P 0.98171239 0.01828761
## 
##    thyroid.surgery
## Y            f          t
##   N 0.99122807 0.00877193
##   P 0.98254364 0.01745636
## 
##    I131.treatment
## Y            f          t
##   N 0.98245614 0.01754386
##   P 0.98420615 0.01579385
## 
##    sick
## Y            f          t
##   N 0.92105263 0.07894737
##   P 0.95261845 0.04738155
## 
##    lithium
## Y             f           t
##   N 0.991228070 0.008771930
##   P 0.995012469 0.004987531
## 
##    goitre
## Y             f           t
##   N 0.991228070 0.008771930
##   P 0.992518703 0.007481297
## 
##    tumor
## Y            f          t
##   N 0.96491228 0.03508772
##   P 0.97838736 0.02161264
## 
##    hypopituitary
## Y              f            t
##   N 0.9912280702 0.0087719298
##   P 0.9991687448 0.0008312552
## 
##    psych
## Y            f          t
##   N 0.97368421 0.02631579
##   P 0.93682461 0.06317539
## 
##    TSH
## Y        [,1]      [,2]
##   N 37.877679 70.383568
##   P  1.888293  6.009277
## 
##    T3
## Y       [,1]      [,2]
##   N 1.425893 0.7839316
##   P 2.025221 0.7287387
## 
##    TT4
## Y        [,1]     [,2]
##   N  70.93482 35.34950
##   P 111.21232 32.15326
## 
##    T4U
## Y        [,1]      [,2]
##   N 1.0155357 0.1973831
##   P 0.9920017 0.1890732
## 
##    FTI
## Y        [,1]     [,2]
##   N  71.26875 34.58914
##   P 113.08243 28.93413

predNB1=predict(NBmod1,testset)

caret::confusionMatrix(reference=testset$binaryClass,
                       data=predNB1,
                       positive="P",
                       mode="everything")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    N    P
##          N   60   13
##          P   85 1726
##                                          
##                Accuracy : 0.948          
##                  95% CI : (0.937, 0.9576)
##     No Information Rate : 0.923          
##     P-Value [Acc > NIR] : 1.136e-05      
##                                          
##                   Kappa : 0.526          
##  Mcnemar's Test P-Value : 7.387e-13      
##                                          
##             Sensitivity : 0.9925         
##             Specificity : 0.4138         
##          Pos Pred Value : 0.9531         
##          Neg Pred Value : 0.8219         
##               Precision : 0.9531         
##                  Recall : 0.9925         
##                      F1 : 0.9724         
##              Prevalence : 0.9230         
##          Detection Rate : 0.9161         
##    Detection Prevalence : 0.9613         
##       Balanced Accuracy : 0.7032         
##                                          
##        'Positive' Class : P              
## 

5 Package klaR

Một package khác cho phép dựng mô hình Naive Bayes đó là klaR.

Mô hình NBmod2 của klaR có nội dung và hiệu năng tương tự như mô hình NBmod1 của e1071.

library(klaR)

NBmod2=NaiveBayes(formula=binaryClass~.,
                  data=trainset,na.action=na.omit,fL=1,
                  usekernel=F)

NBmod2$tables

## $age
##       [,1]     [,2]
## N 52.66071 17.07564
## P 53.60783 19.40314
## 
## $sex
##         var
## grouping         F         M
##        N 0.7982456 0.2017544
##        P 0.6608479 0.3391521
## 
## $pregnant
##         var
## grouping         f          t
##        N 0.9912281 0.00877193
##        P 0.9825436 0.01745636
## 
## $on.thyroxine
##         var
## grouping         f          t
##        N 0.9736842 0.02631579
##        P 0.8935993 0.10640067
## 
## $query.on.thyroxine
##         var
## grouping         f           t
##        N 0.9912281 0.008771930
##        P 0.9916874 0.008312552
## 
## $on.antithyroid.medication
##         var
## grouping         f          t
##        N 0.9912281 0.00877193
##        P 0.9817124 0.01828761
## 
## $thyroid.surgery
##         var
## grouping         f          t
##        N 0.9912281 0.00877193
##        P 0.9825436 0.01745636
## 
## $I131.treatment
##         var
## grouping         f          t
##        N 0.9824561 0.01754386
##        P 0.9842062 0.01579385
## 
## $sick
##         var
## grouping         f          t
##        N 0.9210526 0.07894737
##        P 0.9526185 0.04738155
## 
## $lithium
##         var
## grouping         f           t
##        N 0.9912281 0.008771930
##        P 0.9950125 0.004987531
## 
## $goitre
##         var
## grouping         f           t
##        N 0.9912281 0.008771930
##        P 0.9925187 0.007481297
## 
## $tumor
##         var
## grouping         f          t
##        N 0.9649123 0.03508772
##        P 0.9783874 0.02161264
## 
## $hypopituitary
##         var
## grouping         f            t
##        N 0.9912281 0.0087719298
##        P 0.9991687 0.0008312552
## 
## $psych
##         var
## grouping         f          t
##        N 0.9736842 0.02631579
##        P 0.9368246 0.06317539
## 
## $TSH
##        [,1]      [,2]
## N 37.877679 70.383568
## P  1.888293  6.009277
## 
## $T3
##       [,1]      [,2]
## N 1.425893 0.7839316
## P 2.025221 0.7287387
## 
## $TT4
##        [,1]     [,2]
## N  70.93482 35.34950
## P 111.21232 32.15326
## 
## $T4U
##        [,1]      [,2]
## N 1.0155357 0.1973831
## P 0.9920017 0.1890732
## 
## $FTI
##        [,1]     [,2]
## N  71.26875 34.58914
## P 113.08243 28.93413

predNB2=predict(NBmod2,testset)

caret::confusionMatrix(reference=testset$binaryClass,
                       data=predNB2$class,
                       positive="P",
                       mode="everything")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    N    P
##          N   60   13
##          P   85 1726
##                                          
##                Accuracy : 0.948          
##                  95% CI : (0.937, 0.9576)
##     No Information Rate : 0.923          
##     P-Value [Acc > NIR] : 1.136e-05      
##                                          
##                   Kappa : 0.526          
##  Mcnemar's Test P-Value : 7.387e-13      
##                                          
##             Sensitivity : 0.9925         
##             Specificity : 0.4138         
##          Pos Pred Value : 0.9531         
##          Neg Pred Value : 0.8219         
##               Precision : 0.9531         
##                  Recall : 0.9925         
##                      F1 : 0.9724         
##              Prevalence : 0.9230         
##          Detection Rate : 0.9161         
##    Detection Prevalence : 0.9613         
##       Balanced Accuracy : 0.7032         
##                                          
##        'Positive' Class : P              
## 

Lưu ý là cả 2 mô hình này đều dùng giả định phân bố Gaussian cho tất cả biến liên tục. Ta có thể vẽ biểu đồ để quan sát hiện tượng các biến liên tục đã bị áp đặt phân bố Gaussian như thế nào: đây là hình ảnh hoàn toàn khác so với thực tế :

par(mfrow=c(2,3))

plot(0:0.1, xlim=c(0,98),ylim=c(0,0.03),
     ylab="density",type="n", 
     xlab="age")
curve(dnorm(x, NBmod2$tables$age[1,1], NBmod2$tables$age[1,2]), add=TRUE, col="blue")
curve(dnorm(x, NBmod2$tables$age[2,1], NBmod2$tables$age[2,2]), add=TRUE, col="red")

plot(0:300, xlim=c(0,300),ylim=c(0,0.09),
     ylab="density",type="n", 
     xlab="TSH")
curve(dnorm(x, NBmod2$tables$TSH[1,1], NBmod2$tables$TSH[1,2]), add=TRUE, col="blue")
curve(dnorm(x, NBmod2$tables$TSH[2,1], NBmod2$tables$TSH[2,2]), add=TRUE, col="red")

plot(0:6, xlim=c(0,6),ylim=c(0:1.5),
     ylab="density",type="n", 
     xlab="T3")
curve(dnorm(x, NBmod2$tables$T3[1,1], NBmod2$tables$T3[1,2]), add=TRUE, col="blue")
curve(dnorm(x, NBmod2$tables$T3[2,1], NBmod2$tables$T3[2,2]), add=TRUE, col="red")

plot(0:6, xlim=c(0,6),ylim=c(0,3),
     ylab="density",type="n", 
     xlab="T4U")
curve(dnorm(x, NBmod2$tables$T4U[1,1], NBmod2$tables$T4U[1,2]), add=TRUE, col="blue")
curve(dnorm(x, NBmod2$tables$T4U[2,1], NBmod2$tables$T4U[2,2]), add=TRUE, col="red")

plot(0:400,xlim=c(0,400),ylim=c(0,0.015),
     ylab="density",type="n", 
     xlab="TT4")
curve(dnorm(x, NBmod2$tables$TT4[1,1], NBmod2$tables$TT4[1,2]), add=TRUE, col="blue")
curve(dnorm(x, NBmod2$tables$TT4[2,1], NBmod2$tables$TT4[2,2]), add=TRUE, col="red")

plot(0:300,xlim=c(0,300),ylim=c(0,0.015),
     ylab="density",type="n", 
     xlab="FTI")
curve(dnorm(x, NBmod2$tables$FTI[1,1], NBmod2$tables$FTI[1,2]), add=TRUE, col="blue")
curve(dnorm(x, NBmod2$tables$FTI[2,1], NBmod2$tables$FTI[2,2]), add=TRUE, col="red")

Bây giờ chúng ta thực hiện một mô hình NB khác trong klaR, nhưng lần này có áp dụng phương pháp kernel để khảo sát chính xác hơn phân bố các biến liên tục.

NBmod2k=klaR::NaiveBayes(formula=binaryClass~.,
                        data=trainset,
                        na.action = na.omit,fL=1,
                        usekernel=T)

predNB2k=predict(NBmod2k,testset)

caret::confusionMatrix(reference=testset$binaryClass,
                       data=predNB2k$class,
                       positive="P",
                       mode="everything")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    N    P
##          N  107   43
##          P   38 1696
##                                           
##                Accuracy : 0.957           
##                  95% CI : (0.9468, 0.9657)
##     No Information Rate : 0.923           
##     P-Value [Acc > NIR] : 1.472e-09       
##                                           
##                   Kappa : 0.7021          
##  Mcnemar's Test P-Value : 0.6567          
##                                           
##             Sensitivity : 0.9753          
##             Specificity : 0.7379          
##          Pos Pred Value : 0.9781          
##          Neg Pred Value : 0.7133          
##               Precision : 0.9781          
##                  Recall : 0.9753          
##                      F1 : 0.9767          
##              Prevalence : 0.9230          
##          Detection Rate : 0.9002          
##    Detection Prevalence : 0.9204          
##       Balanced Accuracy : 0.8566          
##                                           
##        'Positive' Class : P               
## 

Mô hình NBmod2k có sử dụng kernel tốt hơn so với mô hình NBmod2 và 1; khi kiểm định trên testset, ta thấy mô hình NBmod2k có độ chính xác sau cân bằng cao hơn rõ rệt: BAC=0.86 so với 0.70, mức độ tương hợp thực tế và chẩn đoán rất cao với Kappa=0.70. Mô hình này rất nhạy và khá đặc hiệu.

Khi dùng phương pháp kernel, phân bố các biến liên tục có dạng như sau:

NBmod3=NaiveBayes(formula=binaryClass~.,
                 data=trainset[,-c(2:14)],
                 na.action = na.omit,fL=1,
                 usekernel=T)

par(mfrow=c(2,3))

plot(NBmod3,legendplot=T,col=c("blue","red"),lty=c(1,1))

Liên hệ giữa mỗi biến liên tục đầu vào và xác suất điều kiện cho nhãn P, N được khảo sát bằng biểu đồ sau:

testdf<-testset%>%mutate(Negative=predNB2k$posterior[,1],
                       Positive=predNB2k$posterior[,2],
                       Label=predNB2k$class)

testdf[,-c(2:14)]%>%gather(age:FTI,key="Features",value="Score")%>%na.omit()%>%
  gather(Negative,Positive,key="Type",value="Posterior")%>%
  ggplot()+
  geom_smooth(aes(x=Score,y=Posterior,
                  col=Type,fill=Type),
              method="loess",alpha=0.2)+
  facet_wrap(~Features,ncol=2,scales = "free")+
  theme_bw()+
  scale_color_manual(values=c("blue","red"))+
  scale_fill_manual(values=c("blue","red"))

Còn đây là hiệu năng của mô hình NBmod2k

testdf[,-c(2:14)]%>%gather(age:FTI,key="Features",value="Score")%>%na.omit()%>%
  mutate(.,Accuracy=if_else(binaryClass==Label,"Correct","Failed"))%>%
  ggplot()+
  geom_bar(aes(x=binaryClass,y=..count..,fill=Accuracy),position="fill",alpha=0.8)+
  theme_bw()+coord_flip()+
scale_fill_manual(values=c("blue","red"))

Ta làm thêm 2 biểu đồ khác để đánh giá về mức quan trọng của từng biến liên tục và nhị phân:

numdfy=data.frame(age=NBmod2k$tables$age$P$y,
                  TSH=NBmod2k$tables$TSH$P$y,
                  T3=NBmod2k$tables$T3$P$y,
                  TT4=NBmod2k$tables$TT4$P$y,
                  T4U=NBmod2k$tables$T4U$P$y,
                  FTI=NBmod2k$tables$FTI$P$y)%>%gather(age:FTI,key="Features",value="Density")

catdf=data.frame(PregF=NBmod2k$tables$pregnant[2,1],
                 ThyroxF=NBmod2k$tables$on.thyroxine[2,1],
                 QueryF=NBmod2k$tables$query.on.thyroxine[2,1],
                 MedF=NBmod2k$tables$on.antithyroid.medication[2,1],
                 SurgF=NBmod2k$tables$thyroid.surgery[2,1],
                 I131F=NBmod2k$tables$I131.treatment[2,1],
                 SickF=NBmod2k$tables$sick[2,1],
                 GoitreF=NBmod2k$tables$goitre[2,1],
                 LithiumF=NBmod2k$tables$lithium[2,1],
                 TumorF=NBmod2k$tables$tumor[2,1],
                 HypoPitF=NBmod2k$tables$hypopituitary[2,1],
                 PsychF=NBmod2k$tables$psych[2,1])

numdfy%>%
  ggplot(aes(x=reorder(Features,Density),y=Density,fill=reorder(Features,Density)))+
  geom_boxplot(show.legend = F)+
  scale_y_continuous(labels=NULL,"Contribution level")+
  scale_x_discrete("Features")+
  coord_flip()+theme_bw()

catdf%>%gather(PregF:PsychF,key="Feature",value="Posterior")%>%
  ggplot(aes(x=reorder(Feature,Posterior),y=Posterior,fill=reorder(Feature,Posterior)))+
  geom_text(aes(label=Feature),nudge_x = 0.1,nudge_y = 0.01)+
  geom_point(shape=21,col="black",size=5,show.legend = F)+
  theme_bw()+
  scale_y_continuous(labels=NULL,"Contribution level")+
  scale_x_discrete(labels=NULL,"Features")

6 Ưu thế và Giới hạn

Sau đây là những ưu thế và nhược điểm của giải thuật Naïve Bayes khi áp dụng cho nghiên cứu y học :

6.1 Ưu thế

Cơ chế hoạt động và kết quả mô hình dựng bằng Naive Bayes tương đồng với quy trình suy luận trong thực hành lâm sàng. Thật vậy, quy trình chẩn đoán mà bạn đang làm hằng ngày trên mỗi bệnh nhân là sự tổng hợp của hàng loạt suy luận Bayes theo kiểu : Với dữ liệu lâm sàng trong bệnh án (tiền sử, bệnh sử, triệu chứng chức năng/thực thể, kết quả xét nghiệm…) thì xác suất mắc bệnh của bệnh nhân là bao nhiêu ?

Về mặt kỹ thuật NB có nhiều ưu thế: nó khai phá dữ liệu với tốc độ cực kì nhanh, vì 2 lý do – thứ nhất NB xét riêng lẻ từng biến nhưng không cần biết mối liên hệ, tổ hợp giữa chúng, thứ hai vì tại mỗi biến, xác suất điều kiện riêng phần được ước lượng đơn giản bằng phép đếm tần suất hoặc giả định phân phối chuẩn. Tương tự, khi thi hành nhiệm vụ trên dữ liệu mới thì Naive Baye cũng cực kì nhanh. Quả thực, hiếm giải thuật nào có tốc độ nhanh trong cả 2 quá trình: học từ dữ liệu và thi hành nhiệm vụ như Naive Bayes. Do đó, NB rất thích hợp cho những bộ dữ liệu kích thước lớn, cả về số lượng biến và số trường hợp.

Một ưu thế hiển nhiên khác là tính phổ quát: Naive Bayes chấp nhận tất cả các loại biến trong dữ liệu đầu vào, từ liên tục, rời rạc cho đến nhị phân.

Mặt khác, NB chỉ cần cỡ mẫu vừa đủ cho mỗi biến, vì thực chất nó không dùng hết toàn bộ từng trường hợp mà chỉ quan tâm đến tỉ lệ phân bố cho mỗi bậc giá trị (biến rời rạc) hoặc đặc tính phân phối (biến liên tục). Cũng vì lý do này NB ít nhạy cảm với nhiễu và chấp nhận dữ liệu bị thiếu sót rải rác cho từng biến. Trong thực hành lâm sàng không phải lúc nào ta cũng thu thập được đầy đủ thông tin, một mô hình chấp nhận thiếu sót dữ liệu như Naive Bayes có thể sẽ có ích trong trường hợp này.

Việc dựng mô hình NB vô cùng đơn giản. Một mặt, do giả định về tính độc lập, những biến số vô hiệu, ảnh hưởng rất ít đến kết quả của mô hình, trong khi các biến số quan trọng vẫn thực hiện độc lập vai trò của mình. Điều này có nghĩa là khi dùng NB, bạn không cần quan tâm đến việc chọn lọc biến số và mất thời gian xét các giả thuyết về tương tác đa chiều, mà chỉ cần chuyển toàn bộ dữ liệu đầu vào cho NB tự lo. Bạn cũng có thể dùng NB vô tư cho những bộ dữ liệu có rất nhiều biến và chưa có bất cứ ý tưởng nào rõ ràng về vai trò từng biến cũng như liên hệ giữa chúng.

Nội dung mô hình Naive Bayes cũng rất đơn giản và dễ hiểu, nó gần giống như một phân tích thống kê đơn biến hàng loạt.

Nếu mục tiêu của nghiên cứu là phân loại chính xác, Naive Bayes có hiệu quả một cách đáng ngạc nhiên. Tuy đơn giản (và phi lý), nhưng hiệu năng của mô hình NB không thua kém bất cứ giải thuật nào khác, kể cả những phương pháp phức tạp hơn rất nhiều. Độ chính xác của Naive Bayes từng được chứng minh trên thực tế, thậm chí cho những bài toán phân loại tới 10-20 nhãn giá trị.

6.2 Hạn chế

Bất lợi lớn nhất của Naive Bayes chính là sự đơn giản quá mức của nó. Giả định về tính độc lập tuyệt đối giữa các biến đầu vào là rất vô lý và hoàn toàn mâu thuẫn với cơ chế sinh lý bệnh. Do đó mô hình Naive Bayes thường không cho phép diễn giải về tương tác đa chiều hoặc khai phá những cơ chế sinh lý bệnh học mới.

Một nhược điểm khác của Naive Bayes đó là nó nhạy cảm với vấn đề mất cân bằng giữa các nhãn phân loại trong dữ liệu. Tuy nhiên đây chỉ là vấn đề kỹ thuật và chúng ta có thể khắc phục nó.

7 Kết luận

Naive Bayes là một giải thuật Machine learning tuy đơn giản nhưng rất hiệu quả khi áp dụng cho bài toán phân loại. Cơ chế hoạt động của Naive Bayes gần như tương đồng với biện luận chẩn đoán trên lâm sàng, ngoại trừ giả định phi lý về tính độc lập giữa các triệu chứng và biomarker. Tuy nhiên so với những phương pháp thống kê cổ điển khác như kiểm định Ki bình phương, t test, kết quả của Naive Baye cung cấp nhiều thông tin hơn, như xác suất điều kiện cho từng nhãn giá trị, có sử dụng xác suất tiền định. Naive Bayes còn cho phép tổng hợp thông tin từ phân tích đơn biến thành quy luật chẩn đoán, và quy luật này rất hiệu quả như ta thấy.