Diễn giải mô hình logistic ở cấp độ cá thể

1 Đặt vấn đề

Hồi quy logistic là một phương pháp rất phổ biến trong nghiên cứu Y học. Đây là một dạng cá biệt của mô hình tuyến tính tổng quát (GLM) với outcome là một biến số nhị phân giới hạn trong khoảng 0-1, giả định có phân bố binomial, được liên kết với công thức mô hình bằng một hàm logit.

Tuy mô hình logistic được xem như một trong những giải thuật cho bài toán xếp loại nhị phân (binary classification) với bản chất hồi quy tuyến tính, hiếm khi hồi quy logistic được vận dụng cho mục tiêu xây dựng quy luật chẩn đoán trong y học, nhưng phần lớn nghiên cứu chỉ sử dụng hồi quy logistic như một công cụ suy diễn thống kê, với mục đích diễn giải vai trò của một hay nhiều biến số (đại lượng, yếu tố ) gây ảnh hưởng lên nguy cơ của một bệnh lý.

Hiếm gặp hơn, hồi quy logistic được sử dụng nhằm xây dựng một thang điểm lâm sàng, nhưng gần như không bao giờ mô hình logistic được trực tiếp ứng dụng một cách thủ công hay tự động trong thực hành lâm sàng để đưa ra quyết định.

Tuy công dụng chính của hồi quy logistic là để diễn giải về vai trò, ảnh hưởng của yếu tố nguy cơ/triệu chứng lâm sàng lên nguy cơ hiện diện của bệnh lý, cho đến nay việc diễn giải này thực chất chỉ có ý nghĩa tổng quát (ở cấp độ quần thể) nhưng không thể áp dụng được cho từng cá thể riêng biệt. Thí dụ, dựa vào mô hình logistic, ta thường diễn giải vai trò của một yếu tố theo kiểu : Sự hiện diện của yếu tố nguy cơ X làm tăng nguy cơ mắc bệnh Y gấp 2-3 lần so với những người không có X, hoặc : mỗi đơn vị gia tăng của biomarker X sẽ làm tăng nguy cơ hiện diện bệnh lý Y từ 2 đến 3 lần.

Việc diễn giải mang tính chất chung chung này tuy cũng hữu ích (nó cho phép xác định được yếu tố nguy cơ, định lượng được nguy cơ) từ đó đưa ra quy tắc về phòng ngừa, trị liệu, chẩn đoán …), tuy nhiên cách làm này cũng có những giới hạn, thí dụ :

• Việc diễn giải dựa trên giả định « đây là một mô hình chính xác, phù hợp với hiện thực ». Trên thực tế, mô hình có thể sai khi áp dụng lên một quần thể đích khác.

• Mỗi predictor được diễn giải độc lập với những predictor còn lại, với giả định là giá trị của chúng giữ nguyên. Giả định này phi lý, vì nếu các predictor thực sự là những biến độc lập, và ngẫu nhiên, mỗi cá thể có thể là một tổ hợp ngẫu nhiên từ bất cứ giá trị nào cho từng predictor, cũng như rất khó ước lượng hiệu quả can thiệp trị liệu lên từng yếu tố riêng lẻ.

2 Mục tiêu bài thực hành

Trong bài thực hành hôm nay, Nhi sẽ giới thiệu với các bạn 2 phương pháp hiện đại cho phép diễn giải nội dung mô hình logistic ở cấp độ cá thể (một có từ năm 2017, một mới xuất hiện vào tháng 4 năm 2018 ), và cả 2 đều có thể thực hiện dễ dàng trong R.

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

Thí dụ minh họa trong bài là một nghiên cứu giả định với nội dung: “Tiên lượng nguy cơ tử vong ở bệnh nhân bị nhiễm khuẩn huyết dựa vào giá trị của 10 chỉ số xét nghiệm sinh hóa”. Dữ liệu gồm có 200 trường hợp bao gồm 107 bệnh nhân sống sót và 93 bệnh nhân tử vong trong thời gian điều trị tại bệnh viện.

Có 10 biến định lượng liên tục, bao gồm:

1. Gamma globulin (U/l)

2. ASAT (U/l)

3. ALAT (U/l)

4. Bilirubin (μmol/l)

5. Urea (mmol/l)

6. creatinine (μmol/l)

7. creatinine clearance (ml/min)

8. erythrocyte sedimentation rate) (mm)

9. c-reactive protein (mg/l)

10. lượng tế bào bạch cầu (×10^9 /l)

library(tidyverse)

data<-read.csv("https://raw.githubusercontent.com/kinokoberuji/R-Tutorials/master/Sepsis.csv",sep=";",dec = ",")%>%
  as_tibble()%>%na.omit()

colnames(data)<-c("Tuvong","GAM","ASAT","ALAT","BIL",
                  "UREA","CREA","CLEA","ESR","CRP","LEUCOS")

data$Tuvong%<>%as.factor()%>%
  recode_factor(.,`0` = "Song",`1` = "Tuvong")

table(data$Tuvong)

## 
##   Song Tuvong 
##    107     93

Đầu tiên, Nhi dùng package rsample để cắt ngẫu nhiên dữ liệu gốc thành 2 phần, phần lớn (90%) sẽ được dụng để dựng mô hình logistic, một phần nhỏ (n=19) sẽ được dùng cho phần diễn giải cá thể.

library(rsample)

## Loading required package: broom

## 
## Attaching package: 'rsample'

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

set.seed(2705)

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

table(trainset$Tuvong)

## 
##   Song Tuvong 
##     97     84

table(testset$Tuvong)

## 
##   Song Tuvong 
##     10      9

Để đảm bảo bí mật của testset đến phút cuối, ta sẽ không nhìn vào tập testset, 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 Sống sót/Tử vong

trainset%>%gather(GAM:LEUCOS,key="Feature",value="Score")%>%
  ggplot(aes(x=Tuvong,y=Score,col=Tuvong))+
  geom_jitter(alpha=0.5,size=2)+
  theme_bw(8)+coord_flip()+
  facet_wrap(~Feature,scales="free",ncol=2)+
  scale_color_manual(values=c("blue","red"))

trainset%>%gather(GAM:LEUCOS,key="Feature",value="Score")%>%
  ggplot(aes(x=Score,fill=Tuvong))+
  geom_density(alpha=0.5,col="black")+
  theme_bw(8)+
  facet_wrap(~Feature,scales="free",ncol=2)+
  scale_fill_manual(values=c("blue","red"))

Toàn bộ các biến số đều cho thấy sự tương phản rõ giữa 2 phân nhóm, các trường hợp tử vong có giá trị cao hơn hẳn so với các trường hợp sống sót

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

Bước tiếp theo Nhi sẽ sơ chế dữ liệu một chút, do chúng ta sử dụng hồi quy tuyến tính, Nhi muốn đưa trung bình của toàn bộ predictor về 0; cũng như loại bỏ những biến số tương quan lẫn nhau. Để 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. Tiếp theo, đưa trung bình về Zero (centering)

3. 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.

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

library(recipes)

rec_obj<-recipe(Tuvong~., data = trainset)%>%
  step_center(all_numeric(),-all_outcomes())%>%
  step_corr(all_numeric(),-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':    181 obs. of  11 variables:
##  $ Tuvong: Factor w/ 2 levels "Song","Tuvong": 1 1 1 1 1 1 1 1 1 1 ...
##  $ GAM   : num  -180 -186 -170 -165 -177 ...
##  $ ASAT  : num  -171 -173 -159 -160 -161 ...
##  $ ALAT  : num  -154 -155 -156 -148 -166 ...
##  $ BIL   : num  -59.3 -58.3 -57.3 -57.3 -57.3 ...
##  $ UREA  : num  -11.67 -13.07 -9.47 -9.07 -8.97 ...
##  $ CREA  : num  -103.4 -125.4 -134.4 -116.4 -97.4 ...
##  $ CLEA  : num  -30.9 -31.9 -35.9 -29.9 -39.9 ...
##  $ ESR   : num  -37 -32 -31 -33 -30 ...
##  $ CRP   : num  -21.1 -20.1 -19.1 -18.1 -17.1 ...
##  $ LEUCOS: num  -7.05 -6.05 -8.05 -5.05 -6.05 ...

Như vậy, sau khi sơ chế không có biến nào bị loại bỏ, toàn bộ 11 biến đi vào quy trình dựng mô hình logistic.

Centering không làm thay đổi đặc tính phân phối, nó chỉ thay đổi thang đo. Công dụng của Centering giúp cho việc diễn giải hệ số hồi quy chính xác hơn.

train_df%>%gather(GAM:LEUCOS,key="Feature",value="Score")%>%
  ggplot(aes(x=Tuvong,y=Score,col=Tuvong))+
  geom_jitter(alpha=0.5,size=2)+
  theme_bw(8)+coord_flip()+
  facet_wrap(~Feature,scales="free",ncol=2)+
  scale_color_manual(values=c("blue","red"))

5 Bước 3: Dựng một mô hình logistic

Tiếp theo, Nhi dùng caret để dựng một mô hình logistic bằng quy trình step-wise kiểm chứng chéo 5 blocks:

# A logistic model by Step-wise

library(caret)

tct<- trainControl(method = "cv", number = 5, returnResamp = "all",
                       classProbs = TRUE, 
                       summaryFunction = twoClassSummary)

log_mod1 <- caret::train(x=train_df[,-1], y=train_df$Tuvong, 
                             trControl = tct,
                             method = "glmStepAIC", 
                             verbose=FALSE,
                             metric = "ROC", 
                             preProc = c("center", "scale"))

## Start:  AIC=39.73
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CLEA    1   17.740 37.740
## - UREA    1   17.742 37.742
## - CRP     1   17.745 37.745
## - ALAT    1   18.083 38.083
## - CREA    1   18.197 38.197
## - ASAT    1   18.469 38.469
## - GAM     1   18.512 38.512
## <none>        17.735 39.735
## - ESR     1   20.043 40.043
## - BIL     1   26.337 46.337
## - LEUCOS  1   32.417 52.417
## 
## Step:  AIC=37.74
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + ESR + CRP + 
##     LEUCOS
## 
##          Df Deviance    AIC
## - UREA    1   17.749 35.749
## - CRP     1   17.750 35.750
## - ALAT    1   18.089 36.089
## - CREA    1   18.437 36.437
## - ASAT    1   18.473 36.473
## - GAM     1   18.513 36.513
## <none>        17.740 37.740
## - ESR     1   20.043 38.043
## - BIL     1   26.588 44.588
## - LEUCOS  1   33.276 51.276
## 
## Step:  AIC=35.75
## .outcome ~ GAM + ASAT + ALAT + BIL + CREA + ESR + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CRP     1    17.90  33.90
## - ALAT    1    18.20  34.20
## - GAM     1    18.51  34.51
## - ASAT    1    18.56  34.56
## <none>         17.75  35.75
## - ESR     1    20.69  36.69
## - CREA    1    22.66  38.66
## - BIL     1    27.76  43.76
## - LEUCOS  1   504.61 520.61
## 
## Step:  AIC=33.9
## .outcome ~ GAM + ASAT + ALAT + BIL + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## - ALAT    1   18.354 32.354
## - GAM     1   18.574 32.574
## - ASAT    1   18.592 32.592
## <none>        17.902 33.902
## - ESR     1   20.783 34.783
## - CREA    1   29.125 43.125
## - BIL     1   31.971 45.971
## - LEUCOS  1   36.616 50.616
## 
## Step:  AIC=32.35
## .outcome ~ GAM + ASAT + BIL + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## - ASAT    1   18.787 30.787
## - GAM     1   19.068 31.068
## <none>        18.354 32.354
## - ESR     1   20.947 32.947
## - CREA    1   29.393 41.393
## - BIL     1   33.111 45.111
## - LEUCOS  1   37.369 49.369
## 
## Step:  AIC=30.79
## .outcome ~ GAM + BIL + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## - GAM     1   19.122 29.122
## <none>        18.787 30.787
## - ESR     1   21.523 31.523
## - CREA    1   29.594 39.594
## - BIL     1   33.291 43.291
## - LEUCOS  1   37.827 47.827
## 
## Step:  AIC=29.12
## .outcome ~ BIL + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## <none>        19.122 29.122
## - ESR     1   21.553 29.553
## - CREA    1   29.709 37.709
## - LEUCOS  1   38.311 46.311
## - BIL     1   63.974 71.974
## Start:  AIC=22
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ASAT    1    0.000 20.000
## - BIL     1    0.000 20.000
## - ESR     1    0.000 20.000
## - CREA    1    0.000 20.000
## - ALAT    1    0.000 20.000
## - GAM     1    0.000 20.000
## - UREA    1    0.000 20.000
## - LEUCOS  1    0.000 20.000
## - CLEA    1    0.000 20.000
## <none>         0.000 22.000
## - CRP     1   12.954 32.954
## 
## Step:  AIC=20
## .outcome ~ GAM + ALAT + BIL + UREA + CREA + CLEA + ESR + CRP + 
##     LEUCOS
## 
##          Df Deviance     AIC
## - BIL     1    0.000  18.000
## - CREA    1    0.000  18.000
## - ESR     1    0.000  18.000
## - ALAT    1    0.000  18.000
## - UREA    1    0.000  18.000
## - GAM     1    0.000  18.000
## - LEUCOS  1    0.000  18.000
## <none>         0.000  20.000
## - CRP     1   13.373  31.373
## - CLEA    1  288.349 306.349
## 
## Step:  AIC=18
## .outcome ~ GAM + ALAT + UREA + CREA + CLEA + ESR + CRP + LEUCOS
## 
##          Df Deviance   AIC
## - ESR     1     0.00 16.00
## - CREA    1     0.00 16.00
## - UREA    1     0.00 16.00
## - LEUCOS  1     0.00 16.00
## - GAM     1     0.00 16.00
## - ALAT    1     0.00 16.00
## - CLEA    1     0.00 16.00
## <none>          0.00 18.00
## - CRP     1    16.14 32.14
## 
## Step:  AIC=16
## .outcome ~ GAM + ALAT + UREA + CREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CREA    1     0.00  14.00
## - GAM     1     0.00  14.00
## - ALAT    1     0.00  14.00
## <none>          0.00  16.00
## - CRP     1    16.28  30.28
## - UREA    1   432.52 446.52
## - CLEA    1   504.61 518.61
## - LEUCOS  1   648.79 662.79
## 
## Step:  AIC=14
## .outcome ~ GAM + ALAT + UREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - GAM     1    0.000 12.000
## - ALAT    1    0.000 12.000
## - LEUCOS  1    0.000 12.000
## - UREA    1    0.000 12.000
## - CLEA    1    0.000 12.000
## <none>         0.000 14.000
## - CRP     1   16.506 28.506
## 
## Step:  AIC=12
## .outcome ~ ALAT + UREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - LEUCOS  1    0.000 10.000
## - CLEA    1    0.000 10.000
## <none>         0.000 12.000
## - UREA    1    6.145 16.145
## - CRP     1   20.762 30.762
## - ALAT    1   37.201 47.201
## 
## Step:  AIC=10
## .outcome ~ ALAT + UREA + CLEA + CRP
## 
##        Df Deviance    AIC
## <none>       0.000 10.000
## - CLEA  1    5.309 13.309
## - UREA  1   19.346 27.346
## - ALAT  1   37.829 45.829
## - CRP   1   69.427 77.427
## Start:  AIC=52.13
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ASAT    1   30.138 50.138
## - ALAT    1   30.141 50.141
## - CREA    1   30.173 50.173
## - CLEA    1   30.494 50.494
## - UREA    1   30.827 50.827
## - ESR     1   31.277 51.277
## - GAM     1   31.580 51.580
## <none>        30.132 52.132
## - CRP     1   32.552 52.552
## - BIL     1   33.725 53.725
## - LEUCOS  1   40.479 60.479
## 
## Step:  AIC=50.14
## .outcome ~ GAM + ALAT + BIL + UREA + CREA + CLEA + ESR + CRP + 
##     LEUCOS
## 
##          Df Deviance    AIC
## - ALAT    1   30.157 48.157
## - CREA    1   30.178 48.178
## - CLEA    1   30.509 48.509
## - UREA    1   30.830 48.830
## - ESR     1   31.318 49.318
## <none>        30.138 50.138
## - GAM     1   32.185 50.185
## - CRP     1   32.557 50.557
## - BIL     1   33.796 51.796
## - LEUCOS  1   40.543 58.543
## 
## Step:  AIC=48.16
## .outcome ~ GAM + BIL + UREA + CREA + CLEA + ESR + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CREA    1   30.192 46.192
## - CLEA    1   30.518 46.518
## - UREA    1   30.856 46.856
## - ESR     1   31.446 47.446
## <none>        30.157 48.157
## - CRP     1   32.693 48.693
## - BIL     1   34.403 50.403
## - GAM     1   35.038 51.038
## - LEUCOS  1   41.301 57.301
## 
## Step:  AIC=46.19
## .outcome ~ GAM + BIL + UREA + CLEA + ESR + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CLEA    1   30.555 44.555
## - UREA    1   30.937 44.937
## - ESR     1   31.561 45.561
## <none>        30.192 46.192
## - CRP     1   32.694 46.694
## - BIL     1   34.620 48.620
## - GAM     1   35.728 49.728
## - LEUCOS  1   41.304 55.304
## 
## Step:  AIC=44.55
## .outcome ~ GAM + BIL + UREA + ESR + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ESR     1   31.626 43.626
## - UREA    1   32.137 44.137
## <none>        30.555 44.555
## - CRP     1   32.946 44.946
## - BIL     1   34.956 46.956
## - GAM     1   36.004 48.004
## - LEUCOS  1   43.575 55.575
## 
## Step:  AIC=43.63
## .outcome ~ GAM + BIL + UREA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## <none>        31.626 43.626
## - CRP     1   34.244 44.244
## - BIL     1   34.972 44.972
## - UREA    1   35.645 45.645
## - GAM     1   38.240 48.240
## - LEUCOS  1   44.463 54.463
## Start:  AIC=41.7
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ALAT    1   19.701 39.701
## - ESR     1   19.740 39.740
## - ASAT    1   19.809 39.809
## - CREA    1   19.842 39.842
## - CLEA    1   20.590 40.590
## - GAM     1   21.146 41.146
## <none>        19.701 41.701
## - UREA    1   21.814 41.814
## - BIL     1   22.382 42.382
## - CRP     1   23.219 43.219
## - LEUCOS  1   33.667 53.667
## 
## Step:  AIC=39.7
## .outcome ~ GAM + ASAT + BIL + UREA + CREA + CLEA + ESR + CRP + 
##     LEUCOS
## 
##          Df Deviance    AIC
## - ESR     1   19.741 37.741
## - ASAT    1   19.840 37.840
## - CREA    1   19.847 37.847
## - CLEA    1   20.657 38.657
## <none>        19.701 39.701
## - GAM     1   21.790 39.790
## - UREA    1   21.822 39.822
## - BIL     1   22.406 40.406
## - CRP     1   23.283 41.283
## - LEUCOS  1   33.831 51.831
## 
## Step:  AIC=37.74
## .outcome ~ GAM + ASAT + BIL + UREA + CREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ASAT    1   19.892 35.892
## - CREA    1   19.894 35.894
## - CLEA    1   20.666 36.666
## <none>        19.741 37.741
## - GAM     1   21.837 37.837
## - UREA    1   22.484 38.484
## - BIL     1   22.557 38.557
## - CRP     1   23.750 39.750
## - LEUCOS  1   36.088 52.088
## 
## Step:  AIC=35.89
## .outcome ~ GAM + BIL + UREA + CREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CREA    1   20.024 34.024
## - CLEA    1   20.793 34.793
## <none>        19.892 35.892
## - UREA    1   22.527 36.527
## - BIL     1   22.564 36.564
## - CRP     1   23.844 37.844
## - GAM     1   26.611 40.611
## - LEUCOS  1   36.191 50.191
## 
## Step:  AIC=34.02
## .outcome ~ GAM + BIL + UREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CLEA    1   21.107 33.107
## <none>        20.024 34.024
## - UREA    1   22.758 34.758
## - BIL     1   22.991 34.991
## - CRP     1   23.861 35.861
## - GAM     1   27.242 39.242
## - LEUCOS  1   36.194 48.194
## 
## Step:  AIC=33.11
## .outcome ~ GAM + BIL + UREA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## <none>        21.107 33.107
## - BIL     1   24.489 34.489
## - CRP     1   24.774 34.774
## - UREA    1   27.349 37.349
## - GAM     1   27.881 37.881
## - LEUCOS  1   36.566 46.566
## Start:  AIC=46.43
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CRP     1   24.483 44.483
## - ALAT    1   24.495 44.495
## - CLEA    1   24.501 44.501
## - ASAT    1   24.890 44.890
## - GAM     1   24.970 44.970
## - CREA    1   25.847 45.847
## <none>        24.433 46.433
## - BIL     1   26.629 46.629
## - UREA    1   29.104 49.104
## - ESR     1   29.473 49.473
## - LEUCOS  1   49.677 69.677
## 
## Step:  AIC=44.48
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     LEUCOS
## 
##          Df Deviance    AIC
## - CLEA    1   24.522 42.522
## - ALAT    1   24.554 42.554
## - ASAT    1   24.962 42.962
## - GAM     1   24.989 42.989
## - CREA    1   25.848 43.848
## <none>        24.483 44.483
## - BIL     1   26.715 44.715
## - ESR     1   29.477 47.477
## - UREA    1   30.598 48.598
## - LEUCOS  1   51.849 69.849
## 
## Step:  AIC=42.52
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## - ALAT    1   24.592 40.592
## - ASAT    1   24.985 40.985
## - GAM     1   25.038 41.038
## - CREA    1   26.019 42.019
## <none>        24.522 42.522
## - BIL     1   26.771 42.771
## - ESR     1   29.479 45.479
## - UREA    1   30.977 46.977
## - LEUCOS  1   52.340 68.340
## 
## Step:  AIC=40.59
## .outcome ~ GAM + ASAT + BIL + UREA + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## - ASAT    1   24.994 38.994
## - GAM     1   25.038 39.038
## - CREA    1   26.061 40.061
## <none>        24.592 40.592
## - BIL     1   27.390 41.390
## - ESR     1   29.527 43.527
## - UREA    1   30.978 44.978
## - LEUCOS  1   52.700 66.700
## 
## Step:  AIC=38.99
## .outcome ~ GAM + BIL + UREA + CREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## - CREA    1   26.319 38.319
## <none>        24.994 38.994
## - BIL     1   27.971 39.971
## - ESR     1   30.208 42.208
## - UREA    1   31.071 43.071
## - GAM     1   34.413 46.413
## - LEUCOS  1   53.070 65.070
## 
## Step:  AIC=38.32
## .outcome ~ GAM + BIL + UREA + ESR + LEUCOS
## 
##          Df Deviance    AIC
## <none>        26.319 38.319
## - BIL     1   28.594 38.594
## - ESR     1   30.610 40.610
## - UREA    1   34.272 44.272
## - GAM     1   34.413 44.413
## - LEUCOS  1   53.073 63.073
## Start:  AIC=55.16
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CREA + CLEA + ESR + 
##     CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CREA    1   33.167 53.167
## - ALAT    1   33.168 53.168
## - ASAT    1   33.191 53.191
## - ESR     1   33.458 53.458
## - CLEA    1   33.826 53.826
## - GAM     1   34.423 54.423
## - UREA    1   34.694 54.694
## <none>        33.164 55.164
## - CRP     1   35.660 55.660
## - BIL     1   37.146 57.146
## - LEUCOS  1   53.298 73.298
## 
## Step:  AIC=53.17
## .outcome ~ GAM + ASAT + ALAT + BIL + UREA + CLEA + ESR + CRP + 
##     LEUCOS
## 
##          Df Deviance    AIC
## - ALAT    1   33.171 51.171
## - ASAT    1   33.194 51.194
## - ESR     1   33.495 51.495
## - CLEA    1   34.163 52.163
## - GAM     1   34.463 52.463
## - UREA    1   34.993 52.993
## <none>        33.167 53.167
## - CRP     1   35.693 53.693
## - BIL     1   37.509 55.509
## - LEUCOS  1   53.310 71.310
## 
## Step:  AIC=51.17
## .outcome ~ GAM + ASAT + BIL + UREA + CLEA + ESR + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ASAT    1   33.215 49.215
## - ESR     1   33.504 49.504
## - CLEA    1   34.164 50.164
## - GAM     1   34.822 50.822
## - UREA    1   35.065 51.065
## <none>        33.171 51.171
## - CRP     1   35.766 51.766
## - BIL     1   38.166 54.166
## - LEUCOS  1   53.488 69.488
## 
## Step:  AIC=49.22
## .outcome ~ GAM + BIL + UREA + CLEA + ESR + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - ESR     1   33.581 47.581
## - CLEA    1   34.234 48.234
## - UREA    1   35.068 49.068
## <none>        33.215 49.215
## - CRP     1   35.768 49.768
## - BIL     1   38.191 52.191
## - GAM     1   39.040 53.040
## - LEUCOS  1   53.791 67.791
## 
## Step:  AIC=47.58
## .outcome ~ GAM + BIL + UREA + CLEA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## - CLEA    1   34.338 46.338
## <none>        33.581 47.581
## - CRP     1   37.362 49.362
## - UREA    1   38.188 50.188
## - BIL     1   38.237 50.237
## - GAM     1   39.617 51.617
## - LEUCOS  1   56.433 68.433
## 
## Step:  AIC=46.34
## .outcome ~ GAM + BIL + UREA + CRP + LEUCOS
## 
##          Df Deviance    AIC
## <none>        34.338 46.338
## - CRP     1   38.169 48.169
## - BIL     1   39.426 49.426
## - GAM     1   39.884 49.884
## - UREA    1   42.461 52.461
## - LEUCOS  1   57.380 67.380

Kết quả của mô hình như sau:

log_mod1

## Generalized Linear Model with Stepwise Feature Selection 
## 
## 181 samples
##  10 predictor
##   2 classes: 'Song', 'Tuvong' 
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 145, 144, 145, 145, 145 
## Resampling results:
## 
##   ROC        Sens       Spec     
##   0.9683185  0.9278947  0.9044118

log_mod1$finalModel%>%tidy()%>%knitr::kable()

term	estimate	std.error	statistic	p.value
(Intercept)	1.345704	1.281448	1.050143	0.2936522
GAM	6.227955	3.254551	1.913614	0.0556694
BIL	3.245595	1.398497	2.320773	0.0202991
UREA	-3.764874	1.453637	-2.589968	0.0095985
CRP	-2.222855	1.010863	-2.198967	0.0278803
LEUCOS	8.088358	2.595259	3.116590	0.0018296

logitsum=log_mod1$finalModel%>%tidy()%>%mutate(OR=exp(estimate),
                                    UL=exp(estimate+1.96*std.error),
                                    LL=exp(estimate-1.96*std.error))

logitsum%>%dplyr::select(term,OR,LL,UL,p.value)%>%knitr::kable()

term	OR	LL	UL	p.value
(Intercept)	3.8408891	0.3116315	4.733933e+01	0.2936522
GAM	506.7183637	0.8598784	2.986044e+05	0.0556694
BIL	25.6769748	1.6562241	3.980784e+02	0.0202991
UREA	0.0231705	0.0013415	4.002181e-01	0.0095985
CRP	0.1082995	0.0149335	7.853992e-01	0.0278803
LEUCOS	3256.3366089	20.1187037	5.270582e+05	0.0018296

Nếu ta xem mô hình trên đây là 1 quy luật chẩn đoán, nó rất chính xác: ROC AUC = 0.96, SENS=0.93 và SPEC=0.91

Mô hình step-wise chỉ giữ lại 5 biến là GAM, BIL, UREA, CRP và LEUCOS.

Từ mô hình gốc, ta có thể tính Odds-ratio và 95%CI của chúng. Việc diễn giải Odds-ratio thì các bạn đã quen nên Nhi không nhắc lại nữa. Ta sẽ đi tiếp sang phần chính của bài đó là diễn giải cục bộ hay cá thể mô hình này.

6 Diễn giải cá thể bằng phương pháp LIME

Năm 2016, Tulio Ribeiro và cs. Giới thiệu phương pháp “Locally Interpretable Model-agnostic Explanations (LIME)”, tạm dịch: “Diễn giải cục bộ cho mô hình bất khả tri”, cho phép diễn giải cơ chế hoạt động của bất cứ mô hình hồi quy và phân loại nào ở cấp độ cá thể, cho mọi dữ liệu đầu vào là văn bản, hình ảnh hay ma trận.

Phương pháp này ban đầu chỉ có trong Python, nhưng sau đó được Thomas Lin Pedersen trang bị cho R qua package lime. Nhi đã trình bày về LIME trong 2 bài tutorial gần đây để diễn giải cho Random Forest và Deep neural network:

http://rpubs.com/lengockhanhi/391615

http://rpubs.com/lengockhanhi/314349

Trong bài thực hành này, chúng ta sẽ áp dụng LIME cho một trường hợp cá biệt là mô hình logistic.

Ta phân chia 19 bệnh nhân trong test_dataframe thành 2 nhóm: 9 trường hợp tử vong và 10 trường hợp sống sót

s_df<-test_df%>%filter(Tuvong=="Song")
d_df<-test_df%>%filter(Tuvong!="Song")

library(lime)

explainer <- lime(train_df,log_mod1)
explanation <- explain(s_df[c(1:6),],explainer,n_features=5,n_labels = 1)

plot_features(explanation)

explainer <- lime(train_df,log_mod1)
explanation <- explain(d_df[c(1:6),],explainer,n_features=5,n_labels = 1)

plot_features(explanation)

explanation2<-explain(
  test_df,
  explainer = explainer, 
  n_labels = 2, 
  n_features =5)

plot_explanations(explanation2)+theme_bw(8)

Do phương pháp LIME dựa vào việc chuyển giá trị định lượng tại một vùng dữ liệu cục bộ thành vector nhị phân rồi dựng một mô hình Whitebox (có thể diễn giải được) từ dữ liệu mô phỏng quanh vị trí này, các biểu đồ diễn giải trình bày mỗi biến số theo quy tắc nhị phân, thí dụ: Bạch cầu có cao hơn 6x 10^9 không ? Creatinin có nằm trong khoảng -114 đến -99 không ?…

Màu xanh tương ứng với “Chứng cứ ủng hộ cho kết quả phân loại”, màu đỏ ngược lại, đại diện cho chứng cứ Chống lại kết quả phân loại. Kích thước của biểu đồ thanh tỉ lệ với hiệu ứng mà mỗi biến số đã góp phần vào việc ủng hộ hay chống lại kết quả mà mô hình đưa ra. Khi giá trị thực của biến số trái với quy tắc nhị phân dành cho nó, hiệu ứng nó gây ra là nghịch với chẩn đoán (làm giảm xác suất cho nhãn giá trị đích), và hiển thị bằng một màu đỏ. Khi giá trị thực của biến số Xi trên cá thể này làm tăng xác suất dự báo cho nhãn giá trị Sống hay tử vong, hiệu ứng của biến số được xem là Ủng hộ, thuận chiều, và biểu thị bằng màu xanh.

7 Diễn giải cá thể bằng phương pháp Breakdown

Vào tháng 4 năm 2018, Mateusz Staniak và Przemysław Biecek thuộc nhóm data scientist MI2 của Ba Lan giới thiệu một phương pháp khác để diễn giải mô hình ở cấp độ cá thể, có tên là « Model agnostic break down » , tạm dịch : « Phân rã mô hình bất khả tri », trong R package breakDown. Package này hoạt động theo một cơ chế khác so với LIME, vì nó nhắm đến mục tiêu phân tích một giá trị prediction xác định (của 1 cá thể) thành từng phần tử nhỏ tương ứng với những predictors (features) đã tham gia trong mô hình. Vai trò đóng góp của mỗi predictors có thể được diễn giải cho cá thể này qua một mô hình tuyến tính :

Dữ liệu của một cá thể (1 hàng trong ma trận test_df) có dạng 1 vector x_new gồm p phần tử từ x1 đến xp, tương ứng giá trị của p predictors trong mô hình. Giá trị prediction là kết quả của hàm f(x_new) với f là model cần diễn giải. Giá trị f(x_new) này được xem như tương đương với một mô hình hồi quy tuyến tính có dạng :

\[f(x^{new}) = (1,x^{new})(\mu , \beta)^{T} = \mu + x_{1}^{new} \beta _{1} + ... x_{p}^{new} \beta _{p}\]

Như vậy các tham số beta_i cho phép ước lượng sự đóng góp / vai trò / ảnh hưởng của từng predictor i trong mô hình logistic (bao gồm Intercept, các dummy variables cho biến định tính và hàm đa thức, nếu có).

Sau khi tải package breakDown vào R, ta viết 1 hàm predict.fun có chức năng trích xuất giá trị xác suất cho 1 cá thể từ mô hình logistic cho trước

Khi test hàm này trên 1 cá thể số 9 từ dữ liệu 10 người sống sót, và cá thể thứ 5 trong số những case tử vong, ta có kết quả như sau:

# Predict function

library("breakDown")

predict.fun <- function(model,instance) predict(model,instance, type = "prob")[,2]

predict.fun(log_mod1, s_df[9,])

## [1] 0.0005495391

predict.fun(log_mod1, d_df[5,])

## [1] 0.9999418

Tiếp theo, ta áp dụng hàm broken trên 1 cá thể tử vong đầu tiên, kết quả như sau:

ex1 <- broken(log_mod1, d_df[1,], 
                   data = train_df[,-1],
                   predict.function = predict.fun)
ex1

##                             contribution
## (Intercept)                        0.464
## + GAM = 700.475138121547           0.535
## + BIL = 225.67955801105            0.001
## + LEUCOS = 10.9502762430939        0.000
## + UREA = 29.9292817679558          0.000
## + CRP = 79.8618784530387           0.000
## + ASAT = 565.309392265193          0.000
## + ALAT = 667.883977900553          0.000
## + CREA = 339.558011049724          0.000
## + CLEA = 72.0718232044199          0.000
## + ESR = 69.9834254143646           0.000
## final_prognosis                    1.000
## baseline:  0

exlist=list(ex1)

Kết quả thô có 2 phần: Cột bên trái là giá trị thực tế của 10 biến của cá thể này (tạm gọi là D1): cột thứ 2 là những hệ số cho biết mức độ đóng góp, hiệu ứng của từng biến lên xác suất dự báo cuối cùng (prob=1, tử vong !). Kết quả cho thấy Gamma là biến quan trọng nhất quyết định kết quả tiên lượng tử vong cho cá thể này, tiếp theo là Bilirubin, ngoài ra chỉ có Intercept của mô hình. Vì mô hình không chứa các biến như ASAT, ALAT, CLEA, ESR…, hiệu ứng của chúng dĩ nhiên = zero.

Ta làm tương tự cho 1 case Sống sót đầu tiên:

ex2 <- broken(log_mod1, s_df[1,], 
              data = train_df[,-1],
              predict.function = predict.fun)
ex2

##                              contribution
## (Intercept)                         0.464
## + LEUCOS = -6.04972375690608       -0.358
## + GAM = -165.524861878453          -0.091
## + BIL = -58.3204419889503          -0.014
## + ASAT = -168.690607734807          0.000
## + ALAT = -175.116022099448          0.000
## + CREA = -125.441988950276          0.000
## + CLEA = -44.9281767955801          0.000
## + ESR = -24.0165745856354           0.000
## + CRP = -16.1381215469613           0.000
## + UREA = -11.0707182320442          0.001
## final_prognosis                     0.001
## baseline:  0

exlist2=list(ex2)

Kết quả diễn giải cho thấy xác suất dự báo sau cùng cho nhãn " tử vong " là 0.001, chứng tỏ mô hình tiên lượng chính xác, bệnh nhân sống. Những biến số như Leucos, Gam, Bil đã góp phần kéo xác suất tử vong về giá trị thấp, trong khi Urea làm tăng xác suất này, nhưng hiệu ứng quá yếu: 0.001; lưu ý là hiệu ứng của Intercept không đổi (xem như baseline).

Tiếp theo, ta viết 2 vòng lặp for_loop để áp dụng quy trình trên hàng loạt cho 8 case tử vong và 9 cases sống sót còn lại:

for(i in 2:nrow(d_df)){
  temp=broken(log_mod1, d_df[i,], 
           data = train_df[,-1],
           predict.function = predict.fun)
  temp=list(temp)
  exlist=c(exlist,temp)
}

for(i in 2:nrow(s_df)){
  temp=broken(log_mod1, s_df[i,], 
              data = train_df[,-1],
              predict.function = predict.fun)
  temp=list(temp)
  exlist2=c(exlist2,temp)
}

Bây giờ ta sẽ vẽ biểu đồ diễn giải cho toàn bộ 19 bệnh nhân :

exlist%>%map(.,plot)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

exlist2%>%map(.,plot)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

19 biểu đồ trên đây là hình ảnh của kết quả diễn giải mô hình logistic của chúng ta cho từng bệnh nhân. Bạn có thể đọc các biểu đồ này để nhận ra ý nghĩa, vai trò của mỗi biến số ở từng cá thể.

Một cách tổng quát, những thông tin sau đây có thể rút ra từ mỗi biểu đồ:

1. Mô hình có hoạt động chính xác hay không cho bệnh nhân này ? (Dựa vào giá trị xác suất tử vong, so với ngưỡng 0.5)

2. Xác suất tử vong của mỗi cá thể là bao nhiêu ?

3. Giá trị của 10 predictors

4. Vai trò đóng góp của mỗi predictor làm tăng/giảm xác suất tử vong ?

5. Những predictor nào quan trọng nhất , ảnh hưởng mạnh nhất quyết địnhcho kết quả tiên lượng ?

6. Những predictor nào làm đảo chiều tiên lượng ?

8 Bàn luận

Nếu chúng ta có thể diễn giải một mô hình logistic cho từng cá thể, đặc biệt là trên cá thể hoàn toàn mới (chưa từng có trong dữ liệu gốc dùng để dựng mô hình), nhiều câu hỏi thú vị hơn về thực hành lâm sàng có thể được trả lời, thí dụ :

1. Quy luật chẩn đoán trong y văn liệu có chính xác/phù hợp cho bệnh nhân Nguyễn Văn A này không ?

2. Nếu có, tôi muốn biết vai trò/ảnh hưởng của từng yếu tố nguy cơ/bệnh cảnh của cá thể này đã làm ảnh hưởng chẩn đoán/loại trừ bệnh lý

3. Bệnh nhân muốn biết dựa vào chứng cứ nào tôi đã chẩn đoán anh ta có nguy cơ tử vong cao ?

4. Tôi phải đưa ra bằng chứng định lượng để biện luận về quyết định chẩn đoán/trị liệu với trưởng khoa và đồng nghiệp

5. Tôi phải ưu tiên can thiệp vấn đề nào cho bệnh nhân này để đảo ngược nguy cơ bệnh lý ? , dựa theo mô hình đã được sử dụng ?

6. Nếu mô hình phạm sai lầm cho bệnh nhân này, tôi muốn biết lý do vì sao ? Bộ phận (biến số) nào trong mô hình đã gây ra sai lầm này ? Hiệu ứng của mỗi predictor có phù hợp với quy luật sinh lý/sinh lý bệnh hay không ?

7. Tôi gặp một trường hợp đặc biệt khó chẩn đoán, và tất cả mô hình trong y văn đều cho ra kết quả ở ranh giới nghi ngờ, tôi muốn báo cáo điều này trong 1 case report

8. Tương tự, tôi gặp một trường hợp kì lạ mà mô hình trong y văn cho ra kết quả chẩn đoán Dương tính trong khi bệnh cảnh lâm sàng không điển hình, tôi muốn kiểm tra vì sao và báo cáo điều này trong 1 case report.

(Cần lưu ý rằng cho đến nay, mặc định là case study và case report không có phân tích thống kê nào được làm cả)

Chúc các bạn thực hành vui !