DALEX: giải pháp toàn diện giải thích mô hình

1 Đặt vấn đề

Trong bài hôm nay, một lần nữa Nhi quay lại với chủ đề diễn giải mô hình. Thực ra không phải chờ đến thời đại hiện nay người ta mới quan tâm đến vai trò quan trọng của mô hình thống kê trong nghiên cứu Y học. Đã từ rắt lâu, khi phân tích dữ liệu bằng thống kê cổ điển, nghiên cứu sinh đã luôn vận dụng mô hình (ngay cả khi họ không ý thức về điều này), mặc dù chỉ có 1 loại mô hình duy nhất là hồi quy tuyến tính.

Như nhận xét của Georges Box : « tất cả mô hình đều sai, nhưng đôi khi chúng hữu dụng », tuy mô hình tuyến tính có nội dung rất đơn giản, nhưng nó vừa đủ để thỏa mãn hầu hết mục tiêu nghiên cứu. Thật vậy, trong hàng chục năm qua, rất nhiều kiến thức y học – bao gồm những quy luật sinh lý-bệnh học, quyết định điều trị và chẩn đoán đã được xác định dựa vào những mô hình tuyến tính đơn giản này.

Sơ đồ sau đây tóm tắt quy trình suy diễn thống kê dựa vào mô hình heo cách làm việc cũ:

Toàn bộ dữ liệu thực nghiệm được dùng để tạo ra mô hình thống kê. Mô hình này có bản chất là sự phản ánh một giả thuyết chủ quan, đơn giản của nhà nghiên cứu. Đây là sự giản lược, thu nhỏ thế giới phức tạp trên thực tế thành một quy luật đơn giản và tương đối chính xác, cho phép trả lời một số câu hỏi nhất định. Phẩm chất của mô hình thường được kiểm chứng trên chính dữ liệu gốc bằng một vài công cụ thống kê đơn giản như LR test, R2, BIC hay AIC… Sau đó việc suy diễn thống kê (thí dụ : so sánh, tương quan, liên hệ, nguy cơ…) được làm ở cấp độ quần thể, trực tiếp dựa vào các hệ số hồi quy trong mô hình. Kết quả suy diễn thống kê được kiểm chứng bằng phản nghiệm giả thuyết vô hiệu, và nếu có ý nghĩa, chúng sẽ được xem như tri thức mới, và quy trình dừng lại ở đó. Hiếm khi mô hình được áp dụng cho mục tiêu thực hành, thí dụ tiên lượng nguy cơ, phân loại chẩn đoán hoặc ước lượng trị số tham chiếu.

Gần đây, với sự phát triển thần tốc của trường phái Machine learning (thực ra: nó đã âm thầm song hành với thống kê quy ước từ thập niên 80), có nhiều ý kiến ủng hộ cho một sự hợp nhất giữa Thống kê cổ điển, Machine learning và khoa học máy tính, để cả ba có thể cùng tham gia hỗ trợ một cách hiệu quả nhất cho các nhà khoa học trong việc khai thác thông tin, kiến thức từ dữ liệu, cũng như những ứng dụng thiết thực và ở cấp độ cá thể.

Sự hợp nhất này không chỉ liên quan đến yếu tố công cụ (du nhập những algorithm mới lạ như Random Forest, SVM, Neural network… vào nghiên cứu khoa học), nhưng còn có thể dẫn đến sự thay đổi mang tính cách mạng về cách suy nghĩ, quy trình làm việc. Một số ý tưởng mới là hệ quả của sự giao thoa giữa 3 chuyên ngành khoa học dữ liệu được tóm tắt trong sơ đồ sau :

Nội dung những ý tưởng này bao gồm :

1. Không chỉ có mô hình tuyến tính, nhưng bất kì một mô hình nào, bao gồm những mô hình hộp đen (blackboxes) đều có thể sử dụng trong nghiên cứu khoa học, cho cả bài toán hồi quy (ước lượng giá trị biến liên tục) hoặc phân loại.

2. Những mô hình này có thể được ứng dụng cho nhiều mục tiêu khác nhau , bao gồm : Tiên lượng (phân loại, ước tính và đưa ra quyết định), kiểm tra phẩm chất của mô hình trên thực tế, suy diễn thống kê ở mức độ quần thể hay cá thể nhằm cung cấp kiến thức mới, kiểm tra giả thuyết nghiên cứu và giải thích cơ chế của mô hình.

3. Quy trình xây dựng mô hình sẽ áp dụng một số quy tắc mới, thí dụ : không sử dụng toàn bộ dữ liệu gốc cho việc tạo ra mô hình, nhưng chuẩn bị riêng các bộ dữ liệu kiểm định độc lập. Một số phương pháp kiểm chứng chéo, tái chọn mẫu (bootstrap) phải được áp dụng. Phẩm chất mô hình sẽ được kiểm tra chặt chẽ hơn bằng nhiều tiêu chí khác nhau trên dữ liệu độc lập trước khi suy diễn thống kê hoặc áp dụng vào thực hành.

4. Việc sử dụng các giải thuật (algorithm) phức tạp khác ngoài mô hình tuyến tính có thể mang lại nhiều lợi ích, như giảm nhẹ công sức của con người trong việc đặt giả thuyết, chọn lọc biến thủ công …, mô hình không phụ thuộc vào các giả định về phân bố nhưng có khả năng phát hiện được những hiệu ứng tương tác, quy luật phi tuyến tính/cục bộ, có thể tự chọn lọc biến, chấp nhận dữ liệu thiếu sót… và nhất là tính chính xác cao hơn.

Tuy nhiên, để có thể thực hiện 4 ý tưởng kể trên cho những nghiên cứu thuộc loại Suy diễn / Diễn dịch, nhà nghiên cứu bắt buộc phải hiểu được mô hình mình đang sử dụng. Việc « Hiểu » này đặt ra nhiều câu hỏi thú vị như :

1. Mô hình này chính xác đến đâu ? Có đáng tin cậy hay không ? Câu hỏi này dẫn đến việc kiểm định phẩm chất mô hình, không chỉ 1 lần trong khi làm nghiên cứu mà trong suốt quá trình lâu dài khi mà mô hình đó vẫn còn được ứng dụng trong thực hành lâm sàng (đối với mục tiêu tiên lượng) hoặc tồn tại trong thư viện y văn (đối với mô hình lý thuyết) .

2. Mô hình hoạt động như thế nào ? Nó đã học được gì từ dữ liệu ? Đây là câu hỏi rất quan trọng làm nền tảng cho việc diễn giải nội dung/cơ chế của mô hình. Giải đáp được câu hỏi này cho phép rút ra hàng loạt thông tin, suy diễn quan trọng, hữu ích bao gồm: Vai trò của mỗi biến ? Mối liên hệ bộ phận giữa mỗi biến và kết quả là gì ? (Ghi chú: nếu kết quả là một chẩn đoán / tiên lượng – câu hỏi này cho phép xác định yếu tố nguy cơ, nếu kết quả là một đại lượng hay con số - câu hỏi trên cho phép so sánh, phân cụm hoặc phân tích tương quan bộ phận). Kết quả của việc giải thích cơ chế này còn có thể được sử dụng để quay ngược lại cải thiện mô hình bằng cách bỏ bớt những biến không quan trọng, tinh chỉnh tham số của algorithm, thêm dữ liệu mới

3. Giải thích được cơ chế của mô hình ở cấp độ cá thể: Mô hình hoạt động có chính xác không cho trường hợp này ? Tại sao kết quả lại như vậy ? Biến nào có vai trò/ ảnh hưởng quan trọng nhất ở cá thể này ? v.v Trong bài trước Nhi đã bàn luận về những lợi ích của việc diễn giải mô hình cho từng cá thể, các bạn có thể đọc thêm ở đây:

2 Giới thiệu DALEX: một giải pháp mới

Trong bài thực hành hôm nay, Nhi sẽ giới thiệu với các bạn một công cụ mới – mà theo đánh giá của Nhi – là một giải pháp hữu hiệu cho phần lớn (nếu không muốn nói là tất cả ) những vấn đề nêu trên. Đó là package DALEX của tác giả Przemyslaw Biecek, mới công bố trên CRAN vào giữa tháng 6 năm 2018 vừa qua.

DALEX là tên rút gọn của « Descriptive mAchine Learning EXplanations ». Tác giả Biecek đã đi xa hơn bất cứ người nào khác trong việc diễn giải nội dung mô hình, với 3 ý tưởng độc đáo:

1. Đưa ra một quy trình giải nghĩa phổ quát cho mọi mô hình, bất kể bản chất của algorithm và mục tiêu nghiên cứu (hồi quy/tiên lượng, phân loại hay suy diễn) . DALEX tương thích với cả 3 giao thức Machine learning phổ biến hiện nay gồm caret, mlr và h2o, cũng như từng algorithm đơn lẻ khác như e1071, rf

2. Cho phép trình bày trực quan kết quả diễn giải của hàng loạt mô hình. Như vậy một công dụng khác của DALEX là so sánh và chọn lọc mô hình/algorithm.

3. Giải đáp hầu hết câu hỏi quan trọng để “hiểu” mô hình, bao gồm: Nội dung và cơ chế hoạt động : Tầm quan trọng của các biến, Quan hệ riêng phần của từng biến đối với kết quả (đặc biệt hữu ích cho bài toán hồi quy), độ chính xác của mô hình và diễn giải cho từng cá thể (theo phương pháp breakdown).

3 Cấu trúc của DALEX

Khi sử dụng DALEX, chỉ có ba bước cần nhớ trong một quy trình chung cho mọi mô hình :

1. Đầu tiên, ta sẽ dùng hàm explain trên mô hình để tạo ra một metadata object cần thiết cho những bước diễn giải tiếp theo.

Hàm explain có cú pháp : explain(model, data, y, predict_function, link, …, label)

Các tùy chỉnh bao gồm:

• model là object kết quả mô hình (DALEX tương thích với object kết quả của caret, mlr, h2o).

• data là một bộ dữ liệu có cấu trúc giống như dữ liệu dùng để dựng mô hình nhưng không gồm biến kết quả. Có thể dùng dữ liệu gốc (trainset) cho mục tiêu giải nghĩa đơn thuần, hoặc dữ liệu mới (testset) cho mục tiêu kiểm định. Phần data này cần cho các quy trình: kiểm định (model performance), khảo sát vai trò các biến (variable importance)

• y là vector biến kết quả từ tập dữ liệu data nói trên, cần để đánh giá variable importance

• predict_function : là 1 hàm do người dùng quy định, cho phép trả kết quả là con số (outcome định lượng cho bài toán hồi quy, hay probability của 1 nhãn trong bài toán phân loại). Nếu không được xác định, hàm predict mặc định trong R sẽ được dùng, việc tùy chỉnh 1 hàm predict riêng là bắt buộc khi mô hình dựng bằng package h2o, cũng như các mô hình binary , multilabel classification…

• link function cho phép tùy chỉnh 1 hàm hoán chuyển kết quả , để chuyển kết quả prediction của mô hình về thang đo và đơn vị nguyên thủy của outcome. Mặc định dùng hàm identity

• label : dùng để đặt cho mỗi mô hình một tên gọi ngắn gọn, tên gọi này sẽ được dùng khi vẽ biểu đồ

2. Từ object này, ta có thể áp dụng tiếp một số hàm chuyên biệt nhằm trả lời cho từng câu hỏi nhất định, thí dụ đánh giá vai trò của các biến, quan hệ phụ thuộc riêng phần, giải thích cho cá thể, …

3. Sau cùng, kết quả từ bước 2 có thể được trình bày trực quan bằng biểu đồ với hàm plot( ). Đặc biệt kết quả diễn giải của nhiều mô hình có thể kết hợp với nhau trong hàm plot cho mục tiêu so sánh.

4 Diễn giải quan hệ bộ phận trong mô hình (1)

Đầu tiên, Ta sẽ dùng một dữ liệu mô phỏng để kiểm tra tính năng khảo sát quan hệ phụ thuộc riêng phần giữa kết quả Y là biến liên tục và những biến số Xi khác trong mô hình.

library(tidyverse)

fsim=function(X1,X2,X3,X4,X5){
  (0.6*X1-0.2)^2+
    sin(8*X2)+cos(5*X2)+
    X3^0.5+
    exp(X4)/(exp(X4)+1)+
    log2(X5)+2
}


print(fsim)

## function(X1,X2,X3,X4,X5){
##   (0.6*X1-0.2)^2+
##     sin(8*X2)+cos(5*X2)+
##     X3^0.5+
##     exp(X4)/(exp(X4)+1)+
##     log2(X5)+2
## }

N=300
set.seed(12345)

X1=runif(N)
X2=runif(N)
X3=runif(N)
X4=runif(N)
X5=runif(N)
Y=fsim(X1,X2,X3,X4,X5)

sim_df=data_frame(X1,X2,X3,X4,X5,Y)

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

X1	X2	X3	X4	X5	Y
0.7209039	0.0526883	0.6992631	0.7621013	0.0508540	0.6492452
0.8757732	0.5986573	0.0104070	0.1196780	0.8032334	0.4356465
0.7609823	0.7082906	0.5039071	0.8480882	0.5837758	1.1997625
0.8861246	0.6534407	0.7310592	0.3288868	0.5854832	0.9118694
0.4564810	0.5774140	0.3297780	0.2841243	0.3210260	-0.4524134
0.1663718	0.6782855	0.8204991	0.0524973	0.9584395	1.6429598

Do dữ liệu là mô phỏng nên ta có thể chủ động xác định trước bản chất quan hệ phụ thuộc riêng phần này bằng các hàm khác nhau

sim_df%>%gather(X1:X5,key="Predictors",value="Value")%>%
  ggplot(aes(x=Value,y=Y,col=Predictors,fill=Predictors))+
  geom_smooth(alpha=0.2)+
  facet_wrap(~Predictors,ncol=3,scales="free")+
  theme_bw()

## `geom_smooth()` using method = 'loess'

Bây giờ giả định mục tiêu của chúng ta là xây dựng một mô hình hồi quy cho phép ước lượng Y từ X1,X2,X3,X4 và X5. Trong thống kê cổ điển, ta chỉ có thể dùng 1 công cụ là mô hình tuyến tính. Nhưng với Machine learning, ta có nhiều, thậm chí rất nhiều công cụ thay thế.

Trong thí dụ này, Nhi sẽ dựng 4 mô hình hồi quy khác nhau: 1 mô hình tuyến tính glm, một mô hình Random Forest, một mô hình Support vector machine (e1071) và 1 mô hình Gradient boosting GBM, sử dụng caret

library(randomForest)
library(e1071)
library(DALEX)

model_rf <- caret::train(Y~.,data= sim_df,method="rf",verbose=F)
model_svm <- svm(Y ~ ., sim_df)
model_gbm <- caret::train(Y~. , data = sim_df, method="gbm",verbose=F)
model_lm <- glm(Y ~ ., sim_df,family="gaussian")

Tiếp theo, giả định câu hỏi nghiên cứu của chúng ta là xác định mối liên hệ bộ phận giữa kết quả Y và mỗi predictors trong 4 mô hình trên: Trước hết ta dùng hàm explain của DALEX để tạo ra 4 objects metadata.

Ngoài ra, ta tạo thêm object thứ 5 chứa mô hình tuyến tính đã được dùng để mô phỏng dữ liệu ban đầu (hàm fsim). Mô hình này chắc chắn sẽ cho ra kết quả gần với sự thực nhất nên được xem như mốc tham chiếu

ex_rf <- explain(model_rf,data=sim_df[,-6],y=sim_df$Y,label="RandomForest")
ex_svm <- explain(model_svm,data=sim_df[,-6],y=sim_df$Y,label="SVM")
ex_lm <- explain(model_lm,data=sim_df[,-6],y=sim_df$Y,label="Linear")
ex_gbm <- explain(model_gbm,data=sim_df[,-6],y=sim_df$Y,label="GBM")

ex_tr <- explain(model_lm,data = sim_df[,-6], y=sim_df$Y,
                 predict_function = function(m, x) fsim(x[,1], x[,2], x[,3], x[,4], x[,5]), label = "True Model")

Để khảo sát liên hệ phụ thuộc riêng phần giữa mỗi predictor và kết quả, ta có thể sử dụng một trong 2 hàm : single_variable( ) hoặc variable_response().

Hai phương pháp có thể được chọn:

1. tùy chỉnh type=”pdp” sẽ áp dụng « Partial Dependence Plot» của Greenwell (2017 : The R Journal 9(1):421–36.https://journal.r-project.org/archive/2017/RJ-2017-016/index.html) ,

2. tùy chỉnh type=”ale” sẽ áp dụng phương pháp “Accumulated Local Effects” của Apley (2017: ALEPlot: Accumulated Local Effects (Ale) Plots and Partial Dependence (Pd) Plots. https://CRAN.R-project.org/package=ALEPlot)

Phương pháp ALE chính xác hơn trong trường hợp có hiện tượng đa cộng tuyến. Các biểu đồ này có ý nghĩa tương tự như “Marginalized effect plot” dành cho mô hình hồi quy, khảo sát quan hệ riêng phần giữa predictor và outcome đều là biến liên tục.

Nhi áp dụng hàng loạt hàm single_variable và kết nối chúng với nhau trong hàm plot cho từng predictor để so sánh kết quả phân tích giữa 5 mô hình

plot(single_variable(ex_rf, "X1"),
     single_variable(ex_svm, "X1"),
     single_variable(ex_lm, "X1"),
     single_variable(ex_gbm, "X1"),
     single_variable(ex_tr, "X1")
     ) +theme_bw() + ggtitle("X1")

plot(single_variable(ex_rf, "X2"),
     single_variable(ex_svm, "X2"),
     single_variable(ex_lm, "X2"),
     single_variable(ex_gbm, "X2"),
     single_variable(ex_tr, "X2")
) +theme_bw() + ggtitle("X2")

plot(single_variable(ex_rf, "X3"),
     single_variable(ex_svm, "X3"),
     single_variable(ex_lm, "X3"),
     single_variable(ex_gbm, "X3"),
     single_variable(ex_tr, "X3")
) +theme_bw() + ggtitle("X3")

plot(single_variable(ex_rf, "X4"),
     single_variable(ex_svm, "X4"),
     single_variable(ex_lm, "X4"),
     single_variable(ex_gbm, "X4"),
     single_variable(ex_tr, "X4")
) +theme_bw() + ggtitle("X4")

plot(single_variable(ex_rf, "X5"),
     single_variable(ex_svm, "X5"),
     single_variable(ex_lm, "X5"),
     single_variable(ex_gbm, "X5"),
     single_variable(ex_tr, "X5")
) +theme_bw() + ggtitle("X5")

Kết quả phân tích mối liên hệ phụ thuộc riêng phần không chỉ gợi ý về quy luật tương quan giữa mỗi predictor và kết quả Y, nhưng đồ thị của các hàm bộ phận này còn cho phép hình dung về cơ chế hoạt động của mỗi giải thuật (algorithm) của từng loại mô hình, giúp ta hiểu chúng đã học được thông tin gì từ dữ liệu và làm cách nào để có thể tái hiện lại mối tương quan thực tế giữaY và predictor.

Kết quả này rất thú vị, vì ta thấy rằng mô hình hồi quy tuyến tính (glm) chỉ có khả năng nhìn mọi thứ như một đường thẳng, dù cho quy luật tương quan trên thực tế là phi tuyến tính. Nó chỉ có thể thích nghi với quy luật này bằng cách gia giảm Intercept và Slope của đồ thị. Trái lại, những algorithm khác như SVM, GBM và RF có khả năng linh hoạt hơn nhiều, cho phép chúng “học” được chính xác hơn các quy luật phi tuyến tính và cố gắng nắm bắt, tái hiện lại những quy luật phức tạp này (bao gồm hàm sin và đa thức bậc cao) theo cách tốt nhất có thể.

5 So sánh phẩm chất giữa nhiều mô hình

DALEX cho phép so sánh hiệu năng giữa nhiều mô hình bằng hàm model_performance, tuy nhiên tính năng này chỉ có ý nghĩa với mô hình hồi quy, vì mặc định tiêu chí được khảo sát là root_mean_square. Tính năng này cũng có vẻ thừa vì ta hoàn toàn có thể tính thủ công hàng chục tiêu chí khác nhau, hoặc dựng 1 confusion matrix cho bài toán phân loại mà không cần dùng đến package nào.

Với thí dụ mô phỏng trên, không ngạc nhiên khi biết rằng mô hình tuyến tính là yếu kém nhất, trong khi Random Forest là chính xác tối ưu với rmse thấp nhất.

mp_rf <- model_performance(ex_rf)
mp_gbm <- model_performance(ex_gbm)
mp_lm <- model_performance(ex_lm)
mp_svm <- model_performance(ex_svm)

plot(mp_rf,mp_gbm,mp_lm,mp_svm,geom="boxplot")+theme_bw()

6 Đánh giá vai trò của từng biến trong mô hình

Khảo sát vai trò của mỗi predictor (feature) trong mô hình cũng là một cách để hiểu cơ chế hoạt động của mô hình. Hầu hết những công cụ Machine learning như mlr, h2o đều có hỗ trợ tính năng khảo sát variable importance được hỗ trợ, tuy nhiên khái niệm “mức độ quan trọng” có thể được định nghĩa theo nhiều cách khác nhau , thí dụ trong package randomForestexplainer có đến 7 tiêu chí khác nhau cho phép xác định một feature có quan trọng hay không.

Package DALEX chỉ sử dụng 1 phương pháp duy nhất khảo sát Variable Importance dựa vào kết quả của hàm loss qua nhiều lượt permutation. Để thử nghiệm tính năng này, Nhi sẽ đưa ra một thí dụ minh họa có thực: Đây là một nghiên cứu của Nierenberg DW et al. năm 1989 nhằm khảo sát quan hệ giữa nồng độ retinol trong máu (biến kết quả) với một số yếu tố khác.

betadf<-read.csv("https://www.openml.org/data/get_csv/52623/plasma_retinol.csv")

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

AGE	SEX	SMOKSTAT	QUETELET	VITUSE	CALORIES	FAT	FIBER	ALCOHOL	CHOLESTEROL	BETADIET	RETDIET	BETAPLASMA	RETPLASMA
64	Female	Former	21.48380	Yes_fairly_often	1298.8	57.0	6.3	0.0	170.3	1945	890	200	915
76	Female	Never	23.87631	Yes_fairly_often	1032.5	50.1	15.8	0.0	75.8	2653	451	124	727
38	Female	Former	20.01080	Yes_not_often	2372.3	83.6	19.1	14.1	257.9	6321	660	328	721
40	Female	Former	25.14062	No	2449.5	97.5	26.5	0.5	332.6	1061	864	153	615
72	Female	Never	20.98504	Yes_fairly_often	1952.1	82.6	16.2	0.0	170.8	2863	1209	92	799
40	Female	Former	27.52136	No	1366.9	56.0	9.6	1.3	154.6	1729	1439	148	654

Đây là một mạng lưới tương quan giữa các biến định lượng trong dữ liệu, các mối liên kết liên tục tương ứng với p_value<0.05 của 1 phân tích tương quan Pearson , ngược lại những đoạn đứt khúc có p_value > 0.05 (không có tương quan).

cormat=betadf%>%
  dplyr::select(-c(SEX,SMOKSTAT,VITUSE))%>%
  as.matrix()%>%
  scale()

m=as.matrix(cor(cormat,
                method="pearson",
                use="pairwise.complete.obs"))

cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

# matrix of the p-value of the correlation
p.mat <- cor.mtest(cormat)

library(igraph)

diag(p.mat)<-0

library(ggraph)

cg1=data.frame(row=rownames(p.mat)[row(p.mat)[upper.tri(p.mat)]], 
               col=colnames(p.mat)[col(p.mat)[upper.tri(p.mat)]], 
               corr=p.mat[upper.tri(p.mat)])

cg2=data.frame(row=rownames(m)[row(m)[upper.tri(m)]], 
               col=colnames(m)[col(m)[upper.tri(m)]], 
               corr=m[upper.tri(m)])

names(cg1)=c("from","to","pval")
names(cg2)=c("from","to","corr")

cg2$pval=cg1$pval
cg2=cg2%>%mutate(pVal=cg1$pval,Sig=if_else(cg1$pval<0.05,1,0))

gdf=filter(cg2,pval<0.05)

graph<-graph_from_data_frame(cg2)

ggraph(graph,circular=F,layout="kk")+
  geom_edge_fan(aes(linetype=factor(Sig),col=factor(Sig)),
                width=1,
                show.legend = F,
                alpha=0.5)+
  geom_node_label(aes(label = name))+
  coord_fixed()+
  theme_graph()+
  scale_edge_color_manual(values=c("red","violet"))

Nhi dựng 3 mô hình hồi quy khác nhau để ước lượng biến RetinolPlasma theo các biến còn lại

model_rf <- caret::train(RETPLASMA~.,data=betadf,method="rf",verbose=F)
model_svm <- svm(RETPLASMA~., betadf)
model_lm <- glm(RETPLASMA~., betadf,family="gaussian")

Sau đó, Nhi cũng tạo ra 3 object metadata cho 3 mô hình này bằng hàm explain.

ex_rf <- explain(model_rf,data=betadf[,-14],y=betadf$RETPLASMA,label="RandomForest")
ex_svm <- explain(model_svm,data=betadf[,-14],y=betadf$RETPLASMA,label="SVM")
ex_lm <- explain(model_lm,data=betadf[,-14],y=betadf$RETPLASMA,label="Linear")

Hàm variable_importance sẽ đánh giá vai trò của từng biến trong mô hình bằng phương pháp permutation của Fisher và công sự (Journal of Computational and Graphical Statistics. http://arxiv.org/abs/1801.01489 ; 2018) :

Cú pháp hàm này:

variable_importance(explainer, loss_function = function(observed, predicted) sum((observed - predicted)^2), …, type = “raw”, n_sample = 1000)

Nhiều phiên bản mô hình sẽ được dựng, trong đó từng biến sẽ được loại bỏ, kết quả của hàm loss_function sẽ lần lượt được tính cho : mô hình bão hòa, mô hình rỗng (cơ bản), và các phiên bản còn lại.

Đối với mô hình tuyến tính, variable importance được suy diễn trực tiếp từ hệ số hồi quy trong mô hình.

Cũng như hàm predict, nội dung hàm loss_function cũng có thể tùy chỉnh theo ý thích. Mặc định là root_mean_square

vi_rf <- variable_importance(ex_rf,type="difference")
vi_svm <- variable_importance(ex_svm,type="difference")
vi_lm <- variable_importance(ex_lm,type="difference")

plot(vi_rf,vi_svm,vi_lm)+theme_bw(8)

7 Diễn giải hiệu ứng bộ phận của biến rời rạc/phân nhóm

Trong thí dụ minh họa ở trên, Nhi đã khảo sát quan hệ phụ thuộc riêng phần/bộ phận giữa kết quả định lượng và predictor định lượng.

Trong trường hợp predictor là biến rời rạc hoặc phân nhóm, phương pháp “Factor Mering Path” của Sitko và Biece (2017) sẽ được áp dụng, với tùy chỉnh: type = “factor”. Kết quả là một biểu đồ so sánh hiệu hứng riêng phần cho mỗi bậc trong factor, hơn nữa 1 dendogram sẽ gom các levels có hiệu ứng tương đương nhau thành 1 cụm (cluster). Ý tưởng này gần giống như phân tích tương phản /post-hoc test trong các thiết kế ANOVA

plot(variable_response(ex_rf, "VITUSE",type="factor")
)+theme_bw(6)+ggtitle("Vitamin A use")

plot(variable_response(ex_svm, "VITUSE",type="factor") )+theme_bw(6)+ggtitle("Vitamin A use")

plot(variable_response(ex_lm, "VITUSE",type="factor")
)+theme_bw(6)+ggtitle("Vitamin A use")

Ta thử khảo sát quan hệ bộ phận giữa Retinol và 2 biến định lượng khác: Cholesterol và Alcohol

plot(single_variable(ex_rf, "ALCOHOL"),
     single_variable(ex_svm, "ALCOHOL"),
     single_variable(ex_lm, "ALCOHOL")
)+theme_bw()+ggtitle("ALCOHOL")

plot(variable_response(ex_rf, "CHOLESTEROL"),
     variable_response(ex_svm, "CHOLESTEROL"),
     variable_response(ex_lm, "CHOLESTEROL")
)+theme_bw()+ggtitle("CHOLESTEROL")

8 Diễn giải mô hình phân loại h2o

Để kết thúc bài này, Nhi sẽ thử áp dụng DALEX cho một mô hình GBM dựng bằng h2o. Thí dụ minh họa này dùng Dữ liệu Thoracic Surgery từ nghiên cứu của Tomczak và cộng sự (Balan) nhằm xây dựng mô hình tiên lượng nguy cơ tử vong trong vòng 1 năm ở bệnh nhân phẫu thuật lồng ngực. Dữ liệu gồm 16 predictors bao gồm triệu chứng chức năng tiền phẫu, chỉ số hô hấp ký, thông tin về chẩn đoán, bệnh lý đi kèm và tình trạng hút thuốc. Nghiên cứu gốc sử dụng giải thuật Support vector machine.

Reference: Tomczak et al. (2013). Boosted SVM for extracting rules from imbalanced data in application to prediction of the post-operative life expectancy in the lung cancer patients. Applied Soft Computing.

Nhi dựng một mô hình GBM bằng package h2o từ 90% dữ liệu gốc

library(h2o)

df=read.csv("https://www.openml.org/data/get_csv/1761654/phpGuu4iR.csv")%>%as_tibble()

colnames(df)<-c("DIAG","FVC","FEV1","ZUBROD","PAIN",
                "HEMOP","DYSPN","COUGH","WEAK","TGRAD",
                "DIAB","MI","PAD","SMOKE","ASTH","AGE","SURV")

df$SURV%<>%as.factor()%>%
  recode_factor(.,`TRUE` = "Dead",`FALSE` = "Survived")

library(rsample)

set.seed(206)

data_split <- initial_split(df, strata = "SURV", prop=0.9)
trainset <- training(data_split)
testset <- testing(data_split)

h2o.init(nthreads = -1,max_mem_size ="4g")

train_hf <- as.h2o(trainset)
test_hf <- as.h2o(testset)

features=setdiff(colnames(train_hf),"SURV")
label="SURV"

gbmod=h2o.gbm(x = features,
               y = label,
               training_frame = train_hf,nfolds=10,
               categorical_encoding="Enum",
               balance_classes = TRUE,
               ntrees =500, max_depth = 10,min_rows = 30,sample_rate=0.8,
               keep_cross_validation_fold_assignment = F, 
               keep_cross_validation_predictions=F,
               score_each_iteration = F,
               seed=206)

Để có thể áp dụng các hàm của package DALEX, ta phải viết 1 hàm predict cho phép xuất kết quả là xác suất của 1 nhãn, thí dụ “Dead”,

Bạn có thể test thử hàm này trên 10 cases đầu tiên của testset nếu bạn thích

prob_predict <- function(model,newdata)  {
  newdata_h2o <- as.h2o(newdata)
  res <- as.data.frame(h2o.predict(model, newdata_h2o))
  return(as.numeric(res$Dead))
}

prob_predict(gbmod,testset[c(1:10),])

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

##  [1] 2.318625e-03 7.419177e-01 8.655132e-07 8.346448e-05 6.471125e-06
##  [6] 1.415870e-03 1.778604e-03 5.450894e-05 6.755577e-06 2.619766e-04

Giả sử ta muốn khảo sát quan hệ riêng phần giữa dung tích sống (FVC) tiền phẫu với xác suất tử vong, ta sẽ có kết quả như sau:

exp_gbm_train<- explain(model = gbmod, 
                             data = trainset[,-17],  
                             y = trainset$SURV,
                             predict_function =prob_predict,
                             label = "h2o_GBM")

pdp_FVC=variable_response(exp_gbm_train, variable = "FVC")

plot(pdp_FVC)+theme_bw()

Tương tự, ta có thể khảo sát quan hệ bộ phận giữa Diagnosis là một biến rời rạc nhiều giá trịvà nguy cơ tử vong

mppDiag=variable_response(exp_gbm_train, variable = "DIAG", type = "factor")

plot(mppDiag)+theme_bw()

9 Diễn giải cho từng cá thể

Sau cùng, DALEX cho phép diễn giải mô hình cho mỗi cá thể bằng phương pháp phân rã mô hình như trong packae breakdown mà Nhi đã giới thiệu ở bài trước, với hàm single_prediction( ).

Thí dụ ta diễn giải kết quả tiên lượng của mô hình cho 1 trường hợp tử vong

new_case <- testset[2,]
brk_gbm1 <- prediction_breakdown(exp_gbm_train,observation = new_case)

plot(brk_gbm1)+theme_bw()

Theo kết quả này, các biến : phân loại T trong TNM, FEV1 thấp, mã chẩn đoán DGN2, điểm Zubrod loại 2, có triệu chứng ho, có hút thuốc lần lượt có ảnh hưởng nhất định đến tiên lượng tử vong, trong khi sự vắng mặt của triệu chứng khó thở, và FVC=3.98 làm giảm nguy cơ tử vong ở bệnh nhân này.

Ta có thể làm tương tự cho 3 mô hình hồi quy ước lượng giá trị Retinol plasma ở trên: Tuy nhiên kết quả diễn giải có vẻ không tương đồng giữa 2 mô hình SVM và Random Forest

bkd_rf=prediction_breakdown(ex_rf,observation=betadf[50,])

bkd_svm=prediction_breakdown(ex_svm,observation=betadf[50,])

plot(bkd_rf,bkd_svm)+theme_bw(8)

10 Nhận xét

Sự góp mặt của package DALEX vào danh sách những công cụ diễn giải nội dung mô hình là một tín hiệu lạc quan cho thấy chúng ta đang tiến rất gần đến một thời kì mới, trong đó mô hình hồi quy tuyến tính không còn là sự lựa chọn duy nhất cho nghiên cứu Y học. Nhiều mô hình khác vốn thuộc về trường phái Statistical learning hay Machine learning cũng có thể được dùng như công cụ để suy diễn thống kê và khai thác thông tin từ dữ liệu thực nghiệm. Thậm chí những mô hình vốn được xem là hộp đen cũng có thể được phân tích về cấu trúc bên trong, cũng như giải thích được cơ chế hoạt động một cách trực quan.