I. Giới thiệu

Thiếu sót dữ liệu là một vấn đề rất khó chịu khi làm nghiên cứu y học lâm sàng. Trên thực tế việc thu thập đầy đủ toàn bộ dữ liệu lâm sàng và cận lâm sàng cho toàn bộ bệnh nhân là bất khả thi, nên cỡ mẫu càng lớn, càng nhiều biến số thì nguy cơ thiếu sót dữ liệu càng cao.

Trong bài hôm nay, chúng ta sẽ cùng khảo sát 5 giải pháp cho phép bổ túc dữ liệu trống. Hoàn cảnh giả định của chúng ta là một nghiên cứu có mục tiêu tiên lượng, dựa trên dataset Pima. Mục tiêu nghiên cứu giả định là xây dựng một mô hình tiên lượng cho phép chẩn đoán bệnh Tiểu đường dựa vào 8 biến số độc lập bao gồm số lần mang thai, nồng độ Glucose/Plasma, Skinfold, huyết áp, nồng độ Insuline, BMI, hàm pedigree và Tuổi.

II. Thăm dò số liệu

Khi mới nhìn vào dataset Pima, tình hình có vẻ rất khả quan; n=768 trường hợp, và có vẻ như không có kí hiệu NA nào cả.

pima=read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data", sep =",",na.strings="0.0",strip.white=TRUE, fill = TRUE)%>%as_tibble()

pima%>%select(.,2,3,4,5,8)%>%map(~as.numeric(.))->pima[,c(2,3,4,5,8)]

pima$V9<-recode_factor(pima$V9,`0`="Normal",`1`="Diabetic")

names(pima)=c("pregnant", "plasmaglu", "DBP", "skinfold", "seruminsulin", "BMI", "pedigreefunction", "age", "Diagnosis")

Tuy nhiên, chỉ với 1 phép kiểm đơn giản, ta bắt đầu nghi ngờ: Ngoại trừ biến pregnant, các biến plasmaGlu, DBP, Skinfold, serum insulin đều có giá trị =0; đây là điều phi lý, vì những đại lượng sinh lý này không thể có giá trị =0 và nếu có thì những trường hợp này phải được xem như Thiếu sót dữ liệu. Do đó ta sẽ thay thế toàn bộ giá trị 0 bằng kí tự NA.

summary(pima==0)
##   pregnant       plasmaglu          DBP           skinfold      
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:657       FALSE:763       FALSE:733       FALSE:541      
##  TRUE :111       TRUE :5         TRUE :35        TRUE :227      
##  NA's :0         NA's :0         NA's :0         NA's :0        
##  seruminsulin       BMI          pedigreefunction    age         
##  Mode :logical   Mode :logical   Mode :logical    Mode :logical  
##  FALSE:394       FALSE:757       FALSE:768        FALSE:768      
##  TRUE :374       NA's :11        NA's :0          NA's :0        
##  NA's :0                                                         
##  Diagnosis      
##  Mode :logical  
##  FALSE:768      
##  NA's :0        
## 
temp=pima[,-1]
temp[temp==0]=NA
pima[,c(2:9)]=temp
rm(temp)

Thiết lập tùy chỉnh cá nhân cho ggplot2

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

theme_set(my_theme())

myfillcolors=c("#bc14ff","#f90e69", "#f9160e" , "#f94d0e", "#f9ab0e","#bbf90e","#0ef988","#0ec3f9")

2 đồ thị sau đây sẽ giúp ta đánh giá tổng quát hơn vấn đề :

Hình 1&2: Thăm dò dữ liệu thiếu sót

bar_missing <- function(x){
  library(dplyr)
  library(reshape2)
  library(ggplot2)
  x %>%
    is.na %>%
    melt %>%
    ggplot(data = .,
           aes(x = Var2)) +
    geom_bar(aes(y=(..count..),fill=value),alpha=0.7)+scale_fill_manual(values=c("skyblue","red"),name = "",
                                                                        labels = c("Available","Missing"))+
    theme_minimal()+
    theme(axis.text.x = element_text(angle=45, vjust=0.5)) +
    labs(x = "Variables in Dataset",
         y = "Observations")+coord_flip()
}


matrix_missing <- function(x){
  library(dplyr)
  library(reshape2)
  library(ggplot2)
  x %>%
    is.na %>%
    melt %>%
    ggplot(data = .,
           aes(x = Var1,
               y = Var2)) +
    geom_tile(aes(fill = value),alpha=0.6) +
    scale_fill_manual(values=c("skyblue","red"),name = "",
                      labels = c("Available","Missing")) +
    theme_minimal()+
    theme(axis.text.x = element_text(angle=45, vjust=0.5)) +
    labs(x = "Variables in Dataset",
         y = "Total observations")+coord_flip()
}

pima[,-9]%>%bar_missing()

pima[,-9]%>%matrix_missing()

Vấn đề thiếu sót dữ liệu có vẻ khá nghiêm trọng : Có khoảng 50% trường hợp bị thiếu dữ liệu (hình 1), tập trung nhiều nhất vào 2 biến là serum insulin (sót 374) và Skinfold (sót 227). Đáng lo nhất là sự thiếu sót phân tán ngẫu nhiên (hình 2).

Dùng hàm na.omit không phải là 1 giải pháp hay (thậm chí rất tệ), vì bạn sẽ mất đi 376 trường hợp một cách hoang phí.

nrow(pima)-nrow(na.omit(pima))
## [1] 376

Một giải pháp khác có thể đề nhị, đó là liệu ta có thể bỏ 2 biến Serum insulin và Skinfold ra khỏi mô hình dự kiến được không ? Muốn trả lời chắc chắn, ta sẽ thực hiện 1 thăm dò bằng cách đánh giá vai trò quan trọng của từng biến vào mục tiêu dự báo, dựa vào mô hình Random Forest:

Hình 3: Vai trò của 8 biến số trong mô hình tiên lượng

library(mlr)

task=makeClassifTask(id = "Original", data=pima, target = "Diagnosis",positive = "Diabetic")

gt=generateFilterValuesData(task,method="gain.ratio")%>%.$data%>%ggplot(aes(x=reorder(name,gain.ratio),y=gain.ratio,fill=reorder(name,-gain.ratio)))+geom_bar(stat="identity",color="black",show.legend=F)+scale_fill_manual(values=myfillcolors,name="Gain ratio")+scale_x_discrete("Features")+coord_flip()

ig=generateFilterValuesData(task,method="information.gain")%>%.$data%>%ggplot(aes(x=reorder(name,information.gain),y=information.gain,fill=reorder(name,-information.gain)))+geom_bar(stat="identity",color="black",show.legend=F)+scale_fill_manual(values=myfillcolors,name="information gain")+scale_x_discrete("Features")+coord_flip()

library(gridExtra)

grid.arrange(gt,ig,ncol=1)

Kết quả thăm dò cho thấy rằng cả 8 biến số độc lập đều CÓ đón góp ít nhiều vào việc chẩn đoán, nhất là biến Serum insulin, do đó giải pháp tối ưu phải làm sao bảo tồn được cả 8 biến trong mô hình.

III. Hướng giải quyết vấn đề:

Có 5 Giải pháp sẽ được đặt ra, cho phép Bổ túc toàn bộ dữ liệu thiếu sót:

  1. KNN Imputation: Bổ túc các trường hợp NA bằng giá trị dự báo của đại lượng liên quan, dựa vào cùng một thuật toán hồi quy KNN

  2. Bagging imputation= Bổ túc các trường hợp NA bằng giá trị dự báo của đại lượng liên quan, dựa vào cùng 1 thuật toán hồi quy Bootstrap aggregating (Bagging)

  3. Median imputation : Thay thế toàn bộ NA bằng giá trị Trung vị của biến số liên quan.

  4. Ensemble imputation: Bổ túc các trường hợp NA bằng giá trị dự báo của 1 thuật toán hồi quy chuyên biệt cho mỗi đại lượng (như vậy ta sẽ tạo ra 1 tập hợp 5 mô hình hồi quy, mỗi mô hình phụ trách bổ túc dữ liệu cho 1 biến số)

  5. Không cần làm gì cả, chỉ việc dùn 1 thuật toán phân loại có khả năng tự xử lý dữ liệu sót, thí dụ như Random Forest SRC hay CART

3 giải pháp đầu tiên được hỗ trợ trong package CARET và được xem như 1 quy trình Preprocess cho dữ liệu. Quy trình Preprocess trong caret được thi hành theo 2 bước: 1) Tạo ra 1 thuật toán (hay hàm) hoán chuyển dữ liệu, 2) Áp dụng hàm này cho đối tượng cần xử lý

Như vậy ngoài dataset pima nguyên thủy, ta sẽ tạo ra thêm lần lượt 3 dataset khác là: knnimp, bagimp và medImp tương ứng với 3 giải pháp vừa nêu trên.

knnimpfunc=caret::preProcess(pima,method="knnImpute")
knnimp=predict(knnimpfunc,pima)

bagimpfunc=caret::preProcess(pima,method="bagImpute")
bagimp=predict(bagimpfunc,pima)

medimpfunc=caret::preProcess(pima,method="medianImpute")
medimp=predict(medimpfunc,pima)

4 dữ liệu trước và sau khi Imput bằng KNN, Bagging và Median

knitr::kable(
  head(pima,10),booktabs = TRUE,
  caption = 'Pima Original dataset'
)
Pima Original dataset
pregnant plasmaglu DBP skinfold seruminsulin BMI pedigreefunction age Diagnosis
6 148 72 35 NA 33.6 0.627 50 Diabetic
1 85 66 29 NA 26.6 0.351 31 Normal
8 183 64 NA NA 23.3 0.672 32 Diabetic
1 89 66 23 94 28.1 0.167 21 Normal
0 137 40 35 168 43.1 2.288 33 Diabetic
5 116 74 NA NA 25.6 0.201 30 Normal
3 78 50 32 88 31.0 0.248 26 Diabetic
10 115 NA NA NA 35.3 0.134 29 Normal
2 197 70 45 543 30.5 0.158 53 Diabetic
8 125 96 NA NA NA 0.232 54 Diabetic 
knitr::kable(
  head(knnimp,10),booktabs = TRUE,
  caption = 'KNN imputation dataset'
)
KNN imputation dataset
pregnant plasmaglu DBP skinfold seruminsulin BMI pedigreefunction age Diagnosis
0.6395305 0.8617221 -0.0327232 0.5580405 0.9198147 0.1649875 0.4681869 1.4250667 Diabetic
-0.8443348 -1.2014407 -0.5172914 -0.0146435 -0.7606615 -0.8458446 -0.3648230 -0.1905477 Normal
1.2330766 2.0079237 -0.6788141 0.1380722 0.3590947 -1.3223797 0.6040037 -0.1055154 Diabetic
-0.8443348 -1.0704463 -0.5172914 -0.5873275 -0.5181880 -0.6292377 -0.9201630 -1.0408711 Normal
-1.1411079 0.5014873 -2.6170869 0.5580405 0.1048342 1.5368309 5.4813370 -0.0204831 Diabetic
0.3427574 -0.1862336 0.1287995 -0.6064169 -0.2807660 -0.9902491 -0.8175458 -0.2755801 Normal
-0.2507887 -1.4306810 -1.8094733 0.2716985 -0.5687033 -0.2104644 -0.6756927 -0.6157094 Diabetic
1.8266227 -0.2189822 0.0964950 0.3289669 0.2193356 0.4104753 -1.0197620 -0.3606124 Normal
-0.5475618 2.4664043 -0.1942460 1.5125138 3.2620416 -0.2826667 -0.9473263 1.6801637 Diabetic
1.2330766 0.1085039 1.9055495 -0.3009855 0.6520835 0.0552400 -0.7239831 1.7651961 Diabetic 
knitr::kable(
  head(bagimp,10),booktabs = TRUE,
  caption = 'Bagging imputation dataset'
)
Bagging imputation dataset
pregnant plasmaglu DBP skinfold seruminsulin BMI pedigreefunction age Diagnosis
6 148 72.00000 35.00000 182.04792 33.60000 0.627 50 Diabetic
1 85 66.00000 29.00000 69.14615 26.60000 0.351 31 Normal
8 183 64.00000 20.55747 209.35756 23.30000 0.672 32 Diabetic
1 89 66.00000 23.00000 94.00000 28.10000 0.167 21 Normal
0 137 40.00000 35.00000 168.00000 43.10000 2.288 33 Diabetic
5 116 74.00000 21.56593 116.60055 25.60000 0.201 30 Normal
3 78 50.00000 32.00000 88.00000 31.00000 0.248 26 Diabetic
10 115 75.00872 32.80297 126.39013 35.30000 0.134 29 Normal
2 197 70.00000 45.00000 543.00000 30.50000 0.158 53 Diabetic
8 125 96.00000 33.39806 163.10080 33.22051 0.232 54 Diabetic 
knitr::kable(
  head(medimp,10),booktabs = TRUE,
  caption = 'Median imputation dataset'
)
Median imputation dataset
pregnant plasmaglu DBP skinfold seruminsulin BMI pedigreefunction age Diagnosis
6 148 72 35 125 33.6 0.627 50 Diabetic
1 85 66 29 125 26.6 0.351 31 Normal
8 183 64 29 125 23.3 0.672 32 Diabetic
1 89 66 23 94 28.1 0.167 21 Normal
0 137 40 35 168 43.1 2.288 33 Diabetic
5 116 74 29 125 25.6 0.201 30 Normal
3 78 50 32 88 31.0 0.248 26 Diabetic
10 115 72 29 125 35.3 0.134 29 Normal
2 197 70 45 543 30.5 0.158 53 Diabetic
8 125 96 29 125 32.3 0.232 54 Diabetic

Để chuẩn bị cho một show diễn rất hay tiếp theo về Machine learning, ta sẽ lần lượt cắt 4 dataset này thành 2 phần với tỉ lệ 75% dùng làm tập huấn luyện, 25% dùng làm tập kiểm định.

Hàm createDataPartition của caret hoạt động cực kì tinh tế, như bạn có thể thấy: Nó thậm chí bảo tồn cả tỉ lệ dữ liệu thiếu sót cho từng biến, tỉ lệ này xấp xỉ như nhau ở tập huấn luyện và kiểm định.

Chú ý: Khi dùng thuật toán KNN để bổ túc dữ liệu, toàn bộ dữ liệu của chúng ta sẽ bị caret chuẩn hóa về thang đo (scaling và centering) ngay cả khi ta không mong muốn điều này. Việc hoán chuyển này không có gì nghiêm trọng, vì mục tiêu của ta là mô hình tiên lượng chứ không phải nghiên cứu diễn dịch, tuy nhiên đây cũng là nhược điểm của phương pháp KNN imputation.

set.seed(1234)
idTrain=caret::createDataPartition(y=pima$Diagnosis, p=0.75,list=FALSE)

train=pima[idTrain,]
test=pima[-idTrain,]

trainKimp=knnimp[idTrain,]
testKimp=knnimp[-idTrain,]

trainbagimp=bagimp[idTrain,]
testbagimp=bagimp[-idTrain,]

trainmedimp=medimp[idTrain,]
testmedimp=medimp[-idTrain,]

Hình 4: Sự tinh tế của hàm creatPartition

mvtrain=bar_missing(train)+ ggtitle("Train subset (n=576)")
mvtest=bar_missing(test)+ ggtitle("Test subset (n=192)")
grid.arrange(mvtrain,mvtest)

Sau đây là phân phối của 8 biến số khi dùng hàm na.omit (pima) và sau khi bổ túc dữ liệu bằng package caret (knnimp, bagimp và medImp)

Hình 5A,B,C,D: Đặc tính phân phối của 4 dữ liệu bao gồm N/A omit và 3 giải pháp bổ túc dữ liệu

mycolors=c("#0e9bf9","#f90e65")

pima%>%gather(pregnant:age,key="Features",value="Value")%>%ggplot()+geom_density(aes(x=Value,fill=Diagnosis),alpha=0.5,show.legend = F)+facet_wrap(~Features,scales="free",ncol=2)+scale_fill_manual(values=mycolors)+ggtitle("NA.Omitted dataset")

knnimp%>%gather(pregnant:age,key="Features",value="Value")%>%ggplot()+geom_density(aes(x=Value,fill=Diagnosis),alpha=0.5,show.legend = F)+facet_wrap(~Features,scales="free",ncol=2)+scale_fill_manual(values=mycolors)+ggtitle("KNN imputed dataset")

bagimp%>%gather(pregnant:age,key="Features",value="Value")%>%ggplot()+geom_density(aes(x=Value,fill=Diagnosis),alpha=0.5,show.legend = F)+facet_wrap(~Features,scales="free",ncol=2)+scale_fill_manual(values=mycolors)+ggtitle("Bagging imputed dataset")

medimp%>%gather(pregnant:age,key="Features",value="Value")%>%ggplot()+geom_density(aes(x=Value,fill=Diagnosis),alpha=0.5,show.legend = F)+facet_wrap(~Features,scales="free",ncol=2)+scale_fill_manual(values=mycolors)+ggtitle("Median imputed dataset")

Nhận xét: Giải pháp Knn làm thay đổi thang đo, nhưng việc bổ túc dữ liệu khá suôn sẻ, vì hình dạng đường cong phân phối là tương đương với dữ liệu hoàn chỉnh (sau khi dùng na.omit).

Giải pháp bagimp bảo tồn được cả thang đo, và không làm thay đổi đáng kể đặc tính phân phối

Giải pháp median imput làm thay đổi khá nhiều về đặc tính phân phối.

Tuy nhiên hiệu quả của chúng sẽ được xét đến sau , mục tiêu của chúng ta là Độ chính xác của mô hình tiên lượng…

IV. Tập hợp mô hình

Giải pháp thứ 4 là bổ túc dữ liệu bằng mô hình tập hợp. Để dễ hình dung về ý tưởng này, bạn hãy liên tưởng đến những bộ phim Science Fiction, hoạt động của một mô hình tập hợp cũng giống như Phi hành đoàn trên tàu USS Enterprise trong phim Star Trek, Biệt đội siêu anh hùng Avengers, hay các mãnh sư Robot trong phim Voltron…

Chúng ta sẽ tạo ra một thứ tương tự như vậy, một tập hợp nhiều mô hình mà mỗi mô hình thành viên đảm nhận 1 nhiệm vụ nhỏ, thí dụ bổ túc dữ liệu cho từng biến số, sau đó sẽ có 1 team leader là 1 thuật toán phân loại, có nhiệm vụ tiên lượng, dự báo cho biến số đích là Chẩn đoán.

Trước tiên, ta cần tuyển chọn 5 thành viên cho mô hình với vai trò bổ túc dữ liệu cho 1 biến số.

Lúc này, mục tiêu bổ túc số liệu có thể xem như 1 bài toán Hồi quy (thí dụ: dự báo giá trị của V3 dựa vào giá trị của V1,V2,V4,V5,V6,V7,V8), ngoại trừ biến số Outcome = diagnosis không tham gia vào bài toán này. Điểm tế nhị, đó là thuật toán hồi quy này phải có khả năng xử lý dữ liệu thiếu sót ở toàn bộ những biến số khác (vì việc thiệu sót có thể xảy ra ngẫu nhiên và đồng thời trên nhiều biến).

Trong package mlr, có hỗ trợ 1 vài thuật toán hồi quy có khả năng tự xử lý dữ liệu sót, như:

listLearners("classif", properties = "missings")[c("class", "package")]
##                         class         package
## 1         classif.bartMachine     bartMachine
## 2          classif.blackboost    mboost,party
## 3            classif.boosting    adabag,rpart
## 4                 classif.C50             C50
## 5             classif.cforest           party
## 6               classif.ctree           party
## 7                 classif.gbm             gbm
## 8                 classif.J48           RWeka
## 9                classif.JRip           RWeka
## 10         classif.naiveBayes           e1071
## 11               classif.OneR           RWeka
## 12               classif.PART           RWeka
## 13    classif.randomForestSRC randomForestSRC
## 14 classif.randomForestSRCSyn randomForestSRC
## 15              classif.rpart           rpart

Ta sẽ so sánh hiệu năng của 4 thuật toán hồi quy là blackboost, cubist, rpart và randomForestSRC cho 5 Task hồi quy tương ứng 5 biến cần bổ túc dữ liệu.

Tiêu chí lựa chọn sẽ là R2 (độ chính xác) và rmse (sai số) , theo đó 1 thuật toán hồi quy sẽ được chọn đảm nhận nhiệm vụ bổ túc cho 1 biến số nếu nó hoạt động chính xác và ổn định hơn 3 thuật toán còn lại.

taskBMI=makeRegrTask(id = "BMI", data=subset(pima,BMI!="NA"),target = "BMI")
taskInsulin=makeRegrTask(id = "Insulin", data=subset(pima,seruminsulin!="NA"),target = "seruminsulin")
taskskinfold=makeRegrTask(id = "Skinfold", data=subset(pima,skinfold!="NA"),target = "skinfold")
taskDBP=makeRegrTask(id = "DBP", data=subset(pima,DBP!="NA"),target = "DBP")
taskGlucose=makeRegrTask(id = "Glucose", data=subset(pima,plasmaglu!="NA"),target = "plasmaglu")

taskIMP=list(taskBMI,taskInsulin,taskskinfold,taskDBP,taskGlucose)

listLearners("regr", properties = "missings")[c("class", "package")]
##                     class         package
## 1         regr.blackboost    mboost,party
## 2            regr.cforest           party
## 3              regr.ctree           party
## 4             regr.cubist          Cubist
## 5                regr.gbm             gbm
## 6    regr.randomForestSRC randomForestSRC
## 7 regr.randomForestSRCSyn randomForestSRC
## 8              regr.rpart           rpart
# 4 learners assigned
implearners = list(
  makeLearner("regr.blackboost",id="BlackBoost"),
  makeLearner("regr.cubist",id="Cubist"),
  makeLearner("regr.rpart",id="DecisionTree"),
  makeLearner("regr.randomForestSRC",id="RandomForestSRC")
)

# This is a Benchmark study
rdesc=makeResampleDesc("CV",iters=10L)
mesreg=list(rsq,rmse)

set.seed(123)
bmrkIMP=benchmark(implearners,taskIMP,rdesc,mesreg)

bmimpdat=getBMRPerformances(bmrkIMP,as.df=TRUE)%>%gather(rsq:rmse,key="Metric",value="Value")

bmimpdat%>%ggplot()+geom_boxplot(aes(x=learner.id,y=Value,fill=learner.id),show.legend = F)+coord_flip()+facet_grid(task.id~Metric,scales="free")+scale_fill_manual(values=myfillcolors)

Kết quả của thí nghiệm thăm dò cho ohe&p ta đưa ra đội hình cho Tập hợp mô hình như sau:

Team leader sẽ là Thuật toán Logistic, nhận nhiệm vụ tiên lượng cho diagnosis. Việc lựa chọn 1 thuật toán đơn giản như Logistic nhằm mục tiêu kép: Thứ nhất nó cho phép diễn giải được nội dung mô hình sau cùng, thứ 2, để chứng minh giả thuyết là sức mạnh tập thể trong một mô hình tập hợp có hiệu năng không thua kém những mô hình phức tạp khác, ngay cả khi team leader là 1 algorithm yếu.

3 thành viên:

Ensemble=makeImputeWrapper("classif.logreg",cols=list(
  plasmaglu=imputeLearner(makeLearner("regr.cubist")),
  BMI=imputeLearner(makeLearner("regr.cubist")),
  DBP=imputeLearner(makeLearner("regr.cubist")),
  skinfold=imputeLearner(makeLearner("regr.blackboost")),
  seruminsulin=imputeLearner(makeLearner("regr.randomForestSRC"))
),
dummy.cols=c("plasmaglu","DBP","skinfold","BMI","seruminsulin"))

Ensemble=setPredictType(Ensemble,"prob")

Mô hình Ensemble của chúng ta sẽ thi đấu với 5 mô hình khác, bao gồm: M1,m2,m3 = 1 thuật toán SVM áp dụng cho 3 dữ liệu đã được bổ túc bởi KNN, Bagging và Median

Như vậy bản chất của 3 mô hình m1, m2, m3 cũng là mô hình tập hợp, team leader của chúng là Support Vector Machine

Mô hình m5 là Random Forest SRC, đại diện cho giải pháp số 5, nó có khả năng chấp nhận và tự bản thân mình xử lý dữ liệu thiếu sót

Mô hình m6 là mô hình cây kiểu CART thông thường. Như ta biết, CART cũng có khả năng tự nó giải quyết dữ liệu thiếu sót, nhưng không tinh tế bằng các mô hình tập hợp Bagging, Boosting khác.

taskNA=makeClassifTask(id = "Original", data=train, target = "Diagnosis",positive = "Diabetic")
taskTESTNA=makeClassifTask(id = "Test", data=test, target = "Diagnosis",positive = "Diabetic")

taskKnn=makeClassifTask(id = "KNN", data=trainKimp, target = "Diagnosis",positive = "Diabetic")
taskTESTKnn=makeClassifTask(id = "TestKNN", data=testKimp, target = "Diagnosis",positive = "Diabetic")

taskBag=makeClassifTask(id = "KNN", data=trainbagimp, target = "Diagnosis",positive = "Diabetic")
taskTESTBag=makeClassifTask(id = "TestKNN", data=testbagimp, target = "Diagnosis",positive = "Diabetic")

taskMed=makeClassifTask(id = "KNN", data=trainmedimp, target = "Diagnosis",positive = "Diabetic")
taskTESTMed=makeClassifTask(id = "TestKNN", data=testmedimp, target = "Diagnosis",positive = "Diabetic")

SVM=makeLearner(id="SVM","classif.svm",predict.type ="prob",fix.factors.prediction = TRUE)

RFSRC=makeLearner(id="RFSRC","classif.randomForestSRC",predict.type ="prob",fix.factors.prediction = TRUE)
cart=makeLearner(id="CART","classif.rpart",predict.type ="prob",fix.factors.prediction = TRUE)

m1=train(SVM,taskKnn)
m2=train(SVM,taskBag)
m3=train(SVM,taskMed)
m4=train(Ensemble,taskNA)
m5=train(RFSRC,taskNA)
m6=train(cart,taskNA)

# Now, we have the verdict:

p1=predict(m1,taskTESTKnn)
p2=predict(m2,taskTESTBag)
p3=predict(m3,taskTESTMed)
p4=predict(m4,taskTESTNA)
p5=predict(m5,taskTESTNA)
p6=predict(m6,taskTESTNA)

mes=list(auc,acc,bac,tpr,tnr,mmce,ber,fpr,fnr,logloss)

pf1=performance(p1,mes)
pf2=performance(p2,mes)
pf3=performance(p3,mes)
pf4=performance(p4,mes)
pf5=performance(p5,mes)
pf6=performance(p6,mes)

pdf=as.data.frame(rbind(pf1,pf2,pf3,pf4,pf5,pf6))
pdf$model=c("KNNIMP","BAGIMP","MEDIMP","ENSEMBLE","RFSRC","CART")

Khi kiểm định độc lập, hiệu năng của 6 mô hình được so sánh như sau:

Hình 7: So sánh hiệu năng của 6 giải pháp bổ túc dữ liệu

library(viridis)

pdflong=gather(pdf,Metrics,Value,auc:logloss,factor_key=TRUE)
ggplot(pdflong)+geom_tile(aes(x=reorder(model,Value),y=reorder(Metrics,-Value),fill=Value),color="black")+geom_text(aes(x=model,y=Metrics,label=round(Value,4)),color="white",fontface = "bold")+scale_fill_viridis(option="A",begin=0.75,end=0.1)+scale_x_discrete("Exams")+ggtitle("Performance of 6 solutions")

Về tiêu chí cân bằng giữa độ chính xác và sai số (dựa vào AUC và BAC): CART là mô hình kém hiệu quả nhất, RFSRC là mô hình hiệu quả nhất, tiếp theo là Tập hợp mô hình có hiệu năng khá cao, gần như tương đương với RFSRC. Cùng một thuật toán SVM nhưng khi áp dụng cho 3 cách bổ túc số liệu khác nhau sẽ cho ra hiệu năng cao thấp khác nhau: Bổ túc bằng Median có hiệu năng tốt hơn so với KNN và Bagging.

Tiếp theo, ta sẽ dùng Kiểm chứng chéo ngẫu nhiên với tỉ lệ 50/50 v à10 lượt để so sánh 3 mô hình: RFSRC, CART và mô hình Ensemble:

Learners = list(RFSRC,Ensemble,cart)

taskRes=makeClassifTask(id = "Original", data=pima, target = "Diagnosis",positive = "Diabetic")

CV=makeResampleDesc("RepCV",folds=2,reps=10L)

bmrk2=benchmark(Learners,taskRes,CV,mes)

bmrkdata=getBMRPerformances(bmrk2,as.df=TRUE)

bmrkdatalong=bmrkdata[,-c(1,5,13)]%>%gather(auc:fnr,key="Metric",value="Value")

Kết quả cho thấy: Mô hình Ensemble có hiệu năng tương đương với mô hình RFSRC, và cả 2 đều tốt hơn so với CART. Thậm chí sai số (mmce) của mô hình Ensemble còn thấp hơn RFSRC.

Hình 8: So sánh hiệu năng giữa 3 mô hình: Ensemble, RFSRC và CART

bmrkdatalong%>%ggplot(aes(x=reorder(Metric,Value),y=Value,fill=Metric))+geom_boxplot(show.legend = F)+facet_grid(learner.id~.,scales="free")+geom_hline(yintercept =0.8,color="blue")+geom_hline(yintercept =0.25,color="red")+coord_flip()+scale_y_continuous(breaks = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1))+scale_fill_manual(values=myfillcolors)

bmrkdatalong%>%ggplot()+geom_path(aes(x=iter,y=Value,color=Metric),size=1,alpha=0.8,show.legend = F)+facet_grid(Metric~learner.id,scales="free")+scale_color_manual(values=myfillcolors)

bmrkdata%>%gather(auc,bac,ber,mmce,logloss,key="Metric",value="Value")%>%ggplot()+geom_histogram(color="black",aes(x=Value,fill=Metric),alpha=0.7,show.legend = F)+facet_grid(learner.id~Metric,scales="free")+scale_fill_manual(values=myfillcolors)

PMCMR::posthoc.kruskal.nemenyi.test(x=bmrkdata$auc,g=bmrkdata$learner.id,method="Tukey")
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  bmrkdata$auc and bmrkdata$learner.id 
## 
##                        RFSRC   classif.logreg.preproc
## classif.logreg.preproc 0.86    -                     
## CART                   5.3e-05 4.6e-06               
## 
## P value adjustment method: none
PMCMR::posthoc.kruskal.nemenyi.test(x=bmrkdata$mmce,g=bmrkdata$learner.id,method="Tukey")
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  bmrkdata$mmce and bmrkdata$learner.id 
## 
##                        RFSRC classif.logreg.preproc
## classif.logreg.preproc 0.290 -                     
## CART                   0.016 6.2e-05               
## 
## P value adjustment method: none

Điều này có nghĩa gì ? Có vẻ như việc kết hợp nhiều mô hình trong 1 tập thể và hoạt động hiệp đồng giữa chúng cho ra hiệu năng tương đương thậm chí tốt hơn so với một thuật toán phức tạp như SVM hay RFSRC.

Nội dung mô hình Logistic Ensemble

summary(m4$learner.model$next.model$learner.model)
## 
## Call:
## stats::glm(formula = f, family = "binomial", data = getTaskData(.task, 
##     .subset), weights = .weights, model = FALSE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6373  -0.7230  -0.3910   0.7055   2.4621  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            -9.530243   0.995779  -9.571  < 2e-16 ***
## pregnant                0.116256   0.037707   3.083  0.00205 ** 
## plasmaglu               0.039393   0.005196   7.582 3.41e-14 ***
## DBP                    -0.013567   0.010279  -1.320  0.18688    
## skinfold               -0.001950   0.016202  -0.120  0.90418    
## seruminsulin           -0.001068   0.001411  -0.757  0.44901    
## BMI                     0.109411   0.023989   4.561 5.09e-06 ***
## pedigreefunction        0.922745   0.358707   2.572  0.01010 *  
## age                     0.013142   0.011437   1.149  0.25054    
## plasmaglu.dummyTRUE     1.056417   1.522373   0.694  0.48773    
## DBP.dummyTRUE           1.111087   0.560874   1.981  0.04759 *  
## skinfold.dummyTRUE     -0.274500   0.325517  -0.843  0.39907    
## BMI.dummyTRUE          -1.683837   1.154381  -1.459  0.14466    
## seruminsulin.dummyTRUE  0.614933   0.301442   2.040  0.04135 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 745.11  on 575  degrees of freedom
## Residual deviance: 531.31  on 562  degrees of freedom
## AIC: 559.31
## 
## Number of Fisher Scoring iterations: 5

Hình 9: Kết quả phân loại của mô hình Ensemble cho toàn bộ dataset pima

p4all=predict(m4,taskRes)

p4df=pima%>%mutate(.,Predicted=p4all$data$response)

plotfuncLow <- function(data,mapping){
  p <- ggplot(data = data,mapping=mapping)+geom_point(aes(fill=p4df$Predicted),shape=21,color="black")+stat_density2d(geom="polygon",aes(fill=p4df$Predicted,alpha = ..level..))+scale_fill_manual(values=mycolors)+geom_rug(aes(color=p4df$Predicted))+scale_color_manual(values=mycolors)
  p
}

plotfuncmid <- function(data,mapping){
  p <- ggplot(data = data,mapping=mapping)+geom_density(aes(fill=p4df$Predicted),alpha=0.5,color="black")+scale_fill_manual(values=mycolors)+geom_rug(aes(color=p4df$Predicted))+scale_color_manual(values=mycolors)
  p
}

library(GGally)

ggpairs(p4df,columns=c(1:8),lower=list(continuous=plotfuncLow),diag=list(continuous=plotfuncmid))

V.Kết luận

Bài hôm nay đưa ra một số thông điệp quan trọng như sau: