mlr3實戰(zhàn) | 基于臨床參數(shù)的肝病患者分類(7種常用的機器學習方法)

圖又掛了，原文鏈接：https://mp.weixin.qq.com/s/bbzCEV7vSubxTIxBOhSrTw

序言

下面的例子是慕尼黑大學機器學習入門講座的一部分內(nèi)容。該項目的目標是為手頭的問題創(chuàng)建并比較一個或幾個機器學習管道，同時進行探索性分析并對結(jié)果進行闡述邓尤。

準備

mlr3的詳細指南見：

mlr3 book (https://mlr3book.mlr-org.com/index.html)

## 安裝與加載所需包
install.packages('mlr3verse')
install.packages('DataExplorer')
install.packages('gridExtra')
library(mlr3verse)
library(dplyr)
library(tidyr)
library(DataExplorer)
library(ggplot2)
library(gridExtra)

用一個固定的種子來初始化隨機數(shù)發(fā)生器，以保證可重復性唧龄，并減少記錄器的冗長性局冰，以保持輸出的清晰表現(xiàn)。

set.seed(7832)
lgr::get_logger("mlr3")$set_threshold("warn")
lgr::get_logger("bbotk")$set_threshold("warn")

在這個示例中训貌，作者研究了機器學習算法和學習器在肝病檢測方面的具體應用制肮。因此，該任務是一個二元分類任務递沪，根據(jù)一些常見的診斷測量結(jié)果預測病人是否患有肝病豺鼻。

印度肝病數(shù)據(jù)

# Importing data
data("ilpd", package = "mlr3data")

它包含了在印度的安得拉邦東北部收集的583名患者的數(shù)據(jù)。根據(jù)病人是否有肝病款慨，觀察結(jié)果被分為兩類儒飒。除了我們的目標變量外，還提供了十個主要是數(shù)字的特征檩奠。為了更詳細地描述這些特征桩了，下表列出了數(shù)據(jù)集中的變量。

顯然埠戳，一些測量值是其它變量的一部分井誉。例如，血清總膽紅素是直接膽紅素和間接膽紅素水平的總和整胃；而白蛋白的數(shù)量則用于計算血清總蛋白以及白蛋白-球蛋白比率的數(shù)值颗圣。因此，一些特征是彼此高度相關(guān)的屁使，下面會進行處理在岂。

數(shù)據(jù)預處理

單變量分布

接下來，研究每個變量的單變量分布蛮寂。從目標變量和唯一的離散特征--性別開始蔽午，它們都是二元變量的。

##  所有離散變量的頻率分布
plot_bar(ilpd,ggtheme = theme_bw())

[圖片上傳失敗...(image-1454a9-1649651130723)]

可以看到共郭，目標變量（即肝病與非肝病患者）的分布是相當不平衡的祠丝，如柱狀圖所示：有肝病和無肝病的患者數(shù)量分別為416和167疾呻。一個類別的代表性不足，可能會使ML模型的性能惡化写半。為了研究這個問題岸蜗，作者還在一個數(shù)據(jù)集上擬合了模型，在這個數(shù)據(jù)集上叠蝇，隨機地對少數(shù)人類別進行了過度抽樣璃岳，結(jié)果是一個完全平衡的數(shù)據(jù)集。此外悔捶，我們還應用了分層抽樣铃慷，以確保在交叉驗證過程中保持各類的比例。唯一的離散特征gender也是相當不平衡的蜕该。

## 查看所有連續(xù)變量的頻率分布直方圖
plot_histogram(ilpd,ggtheme = theme_mlr3())

[圖片上傳失敗...(image-6a00fa-1649651130723)]

可以看到犁柜，一些指標特征是極度右偏的，包含幾個極端值堂淡。為了減少離群值的影響馋缅，并且由于一些模型假設(shè)了特征的正態(tài)性，我們對這些變量進行了log轉(zhuǎn)換绢淀。

特征分組

為了描繪目標和特征之間的關(guān)系萤悴，我們按類別分析了特征的分布情況。首先皆的，我們研究了離散特征性別覆履。

plot_bar(ilpd,by = 'diseased',ggtheme = theme_mlr3())

[圖片上傳失敗...(image-ddfd77-1649651130723)]

在 "疾病 "類中，男性的比例略高费薄，但總體而言硝全，差異不大。除此之外楞抡，正如我們之前提到的柳沙，在兩個類別中都可以觀察到性別不平衡的現(xiàn)象。

為了看到連續(xù)特征的差異拌倍，我們比較了以下的boxplots，其中右偏的特征還沒有進行對數(shù)轉(zhuǎn)換噪径。

## View bivariate continuous distribution based on `diseased`
plot_boxplot(ilpd,by = 'diseased')

[圖片上傳失敗...(image-d03894-1649651130723)]

可以看到除了total_protein柱恤，對于每一個特征，我們都得到了兩個類的中位值之間的差異找爱。值得注意的是梗顺，在強右偏的特征中，"疾病 "類包含的極端值遠遠多于 "無疾病 "類车摄，這可能是因為其規(guī)模較大寺谤。

從下面的圖中可以看出仑鸥，這種影響在對數(shù)轉(zhuǎn)換后會被削弱。此外变屁，這些特征在 "疾病 "類中的分散性更大眼俊，正如箱線圖的長度所示∷诠兀總的來說疮胖，這些特征似乎與目標相關(guān)，所以將它們用于這項任務并建立它們與目標的關(guān)系模型是有意義的闷板。

對部分特征進行l(wèi)og轉(zhuǎn)換

ilpd_log = ilpd %>%
  mutate(
    # Log for features with skewed distributions
    alanine_transaminase = log(alanine_transaminase),
    total_bilirubin =log(total_bilirubin),
    alkaline_phosphatase = log(alkaline_phosphatase),
    aspartate_transaminase = log(aspartate_transaminase),
    direct_bilirubin = log(direct_bilirubin)
  )
plot_histogram(ilpd_log,ggtheme = theme_mlr3(),ncol = 3)
plot_boxplot(ilpd_log,by = 'diseased')

[圖片上傳失敗...(image-e00fb9-1649651130723)]

[圖片上傳失敗...(image-dd3dd1-1649651130723)]

可以看到log轉(zhuǎn)換后的數(shù)據(jù)分布改善了許多澎灸。

相關(guān)分析

正如我們在數(shù)據(jù)描述中提到的，有些特征是由另一個特征間接測量的遮晚。這表明它們是高度相關(guān)的性昭。我們要比較的一些模型假設(shè)是獨立的特征，或者有多重共線性的問題县遣。因此糜颠，我們檢查了特征之間的相關(guān)性。

plot_correlation(ilpd)

[圖片上傳失敗...(image-bdd2df-1649651130723)]

可以看到艺玲，其中四對有非常高的相關(guān)系數(shù)括蝠。看一下這些特征饭聚，很明顯它們是相互影響的忌警。由于模型的復雜性應該最小化，并且由于多重共線性的考慮秒梳，我們決定每對特征中只取一個法绵。在決定保留哪些特征時，我們選擇了那些關(guān)于肝病的更具體和相關(guān)的特征酪碘。因此朋譬，我們選擇了白蛋白，而不是白蛋白和球蛋白的比例兴垦，也不是蛋白質(zhì)的總量徙赢。同樣的觀點也適用于使用直接膽紅素的量而不是總膽紅素。關(guān)于天門冬氨酸轉(zhuǎn)氨酶和丙氨酸轉(zhuǎn)氨酶探越，我們沒有注意到這兩個特征的數(shù)據(jù)有任何根本性的差異狡赐，所以我們?nèi)我膺x擇了天冬氨酸轉(zhuǎn)氨酶。

最終數(shù)據(jù)集

## Reducing, transforming and scaling dataset
ilpd = ilpd %>%
  select(-total_bilirubin, -alanine_transaminase, -total_protein,
         -albumin_globulin_ratio) %>%
  mutate(
    # Recode gender
    gender = as.numeric(ifelse(gender == "Female", 1, 0)),
     # Remove labels for class
    diseased = factor(ifelse(diseased == "yes", 1, 0)),
     # Log for features with skewed distributions
    alkaline_phosphatase = log(alkaline_phosphatase),
    aspartate_transaminase = log(aspartate_transaminase),
    direct_bilirubin = log(direct_bilirubin)
  )
## 標準化
po_scale = po("scale")
po_scale$param_set$values$affect_columns =
  selector_name(c("age", "direct_bilirubin", "alkaline_phosphatase",
                  "aspartate_transaminase", "albumin"))
task_liver = as_task_classif(ilpd_m, target = "diseased", positive = "1")
ilpd_f = po_scale$train(list(task_liver))[[1]]$data()

最后钦幔，我們對所有的連續(xù)變量特征進行了標準化枕屉，這對k-NN模型尤其重要。下表顯示了最終的數(shù)據(jù)集和我們應用的轉(zhuǎn)換鲤氢。注意：與對數(shù)或其他轉(zhuǎn)換不同搀擂，縮放取決于數(shù)據(jù)本身西潘。在數(shù)據(jù)被分割之前對數(shù)據(jù)進行縮放會導致數(shù)據(jù)泄露(詳見：Nature Reviews Genetics | 在基因組學中應用機器學習的常見陷阱），因為訓練集和測試集的信息是共享的哨颂。由于數(shù)據(jù)泄露會導致更高的性能喷市，縮放應該總是單獨應用于ML工作流程所引起的每個數(shù)據(jù)分割。因此咆蒿，我們強烈建議在這種情況下使用PipeOpScale东抹。

學習器和調(diào)參

首先，我們需要定義一個task沃测，其中包含最終的數(shù)據(jù)集和一些元信息缭黔。此外，我們還需要指定正類蒂破，因為軟件包默認將第一個正類作為正類馏谨。正類的指定對后面的評估有影響。

## Task definition
task_liver = as_task_classif(ilpd_f, target = "diseased", positive = "1")

下面我們將對logistic regression, linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), naive Bayes, k-nearest neighbour (k-NN), classification trees (CART) and random forest的二元分類目標進行評估附迷。

# detect overfitting
install.packages('e1071')
install.packages('kknn')
learners = list(
  learner_logreg = lrn("classif.log_reg", predict_type = "prob",
                       predict_sets = c("train", "test")),
  learner_lda = lrn("classif.lda", predict_type = "prob",
                    predict_sets = c("train", "test")),
  learner_qda = lrn("classif.qda", predict_type = "prob",
                    predict_sets = c("train", "test")),
  learner_nb = lrn("classif.naive_bayes", predict_type = "prob",
                   predict_sets = c("train", "test")),
  learner_knn = lrn("classif.kknn", scale = FALSE,
                    predict_type = "prob"),
  learner_rpart = lrn("classif.rpart",
                      predict_type = "prob"),
  learner_rf = lrn("classif.ranger", num.trees = 1000,
                   predict_type = "prob")
)

調(diào)參

為了找到最佳的超參數(shù)惧互，我們使用隨機搜索來更好地覆蓋超參數(shù)空間。我們定義了要調(diào)整的超參數(shù)喇伯。我們只調(diào)整了k-NN喊儡、CART和隨機森林的超參數(shù)，因為其他方法有很強的假設(shè)稻据，并作為基線艾猜。

對于k-NN，我們選擇3作為k（鄰居數(shù)量）的下限捻悯，50作為上限匆赃。太小的k會導致過度擬合。我們還嘗試了不同的距離測量方法（Manhattan distance為1, Euclidean distance為2）和內(nèi)核今缚。對于CART算柳，我們調(diào)整了超參數(shù)cp（復雜度參數(shù)）和minsplit（為了嘗試分割，一個節(jié)點中的最小觀察數(shù)）姓言。cp控制了tree的大兴蚕睢：小的值會導致過擬合，而大的值會導致欠擬合何荚。我們還調(diào)整了隨機森林的終端節(jié)點的最小尺寸和每次分裂時隨機抽樣作為候選變量的數(shù)量（從1到特征數(shù)）的參數(shù)滥壕。

tune_ps_knn = ps(
  k = p_int(lower = 3, upper = 50), # Number of neighbors considered
  distance = p_dbl(lower = 1, upper = 3),
  kernel = p_fct(levels = c("rectangular", "gaussian", "rank", "optimal"))
)
tune_ps_rpart = ps(
  # Minimum number of observations that must exist in a node in order for a
  # split to be attempted
  minsplit = p_int(lower = 10, upper = 40),
  cp = p_dbl(lower = 0.001, upper = 0.1) # Complexity parameter
)
tune_ps_rf = ps(
  # Minimum size of terminal nodes
  min.node.size = p_int(lower = 10, upper = 50),
  # Number of variables randomly sampled as candidates at each split
  mtry = p_int(lower = 1, upper = 6)
)

下一步是將mlr3tuning中的AutoTuner實例化。我們對嵌套重采樣的內(nèi)循環(huán)采用了5-fold交叉驗證法兽泣。評價次數(shù)被設(shè)定為100次作為停止標準。我們使用AUC作為評價指標胁孙，唠倦。

如前所述称鳞，由于目標類別不平衡，我們選擇了完美平衡類稠鼻。通過使用mlr3pipelines冈止，我們可以在以后應用基準函數(shù)。

# Oversampling minority class to get perfectly balanced classes
po_over = po("classbalancing", id = "oversample", adjust = "minor",
             reference = "minor", shuffle = FALSE, ratio = 416/167)
table(po_over$train(list(task_liver))$output$truth()) # Check class balance

# Learners with balanced/oversampled data
learners_bal = lapply(learners, function(x) {
  GraphLearner$new(po_scale %>>% po_over %>>% x)
})
lapply(learners_bal, function(x) x$predict_sets = c("train", "test"))

模型擬合和基準設(shè)定

在定義了學習器候齿、選擇了嵌套重采樣的內(nèi)部方法和設(shè)置了調(diào)整器之后熙暴，我們開始選擇外部重采樣方法。我們選擇了分層的5倍交叉驗證法慌盯，以保持目標變量的分布周霉，不受過度采樣的影響。然而亚皂，事實證明俱箱，沒有分層的正常交叉驗證法也會產(chǎn)生非常相似的結(jié)果。

# 5-fold cross-validation
resampling_outer = rsmp(id = "cv", .key = "cv", folds = 5L)

# Stratification
task_liver$col_roles$stratum = task_liver$target_names

為了對不同的學習器進行排名灭必，并最終決定哪一個最適合手頭的任務狞谱，我們使用了基準測試(benchmarking)。下面的代碼塊執(zhí)行了我們對所有學習者的基準測試禁漓。

design = benchmark_grid(
  tasks = task_liver,
  learners = c(learners, learners_bal),
  resamplings = resampling_outer
)

bmr = benchmark(design, store_models = FALSE) ## 耗時較長

如上所述跟衅，我們選擇了分層的5折交叉驗證法。這意味著性能被確定為五個模型評估的平均值播歼，train-test-split為80%和20%伶跷。此外，性能指標的選擇對于不同學習器的排名至關(guān)重要荚恶。雖然每一個都有其特定的使用情況撩穿，但我們選擇了AUC，一個同時考慮了敏感性和特異性的性能指標谒撼，我們也使用它來進行超參數(shù)調(diào)整食寡。

我們首先通過AUC對所有學習者進行了比較，包括有無超采樣廓潜，以及訓練和測試數(shù)據(jù)抵皱。

measures = list(
  msr("classif.auc", predict_sets = "train", id = "auc_train"),
  msr("classif.auc", id = "auc_test")
)

tab = bmr2$aggregate(measures)
tab_1 = tab[,c('learner_id','auc_train','auc_test')]
print(tab_1)

> print(tab_1)
                               learner_id auc_train  auc_test
 1:                       classif.log_reg 0.7548382 0.7485372
 2:                           classif.lda 0.7546522 0.7487159
 3:                           classif.qda 0.7683438 0.7441634
 4:                   classif.naive_bayes 0.7539374 0.7498427
 5:                    classif.kknn.tuned 0.8652143 0.7150679
 6:                   classif.rpart.tuned 0.7988561 0.6847818
 7:                  classif.ranger.tuned 0.9871615 0.7426650
 8:      scale.oversample.classif.log_reg 0.7540066 0.7497002
 9:          scale.oversample.classif.lda 0.7537952 0.7489675
10:          scale.oversample.classif.qda 0.7679012 0.7481963
11:  scale.oversample.classif.naive_bayes 0.7536208 0.7503436
12:   scale.oversample.classif.kknn.tuned 0.9982251 0.6870297
13:  scale.oversample.classif.rpart.tuned 0.8903927 0.6231100
14: scale.oversample.classif.ranger.tuned 1.0000000 0.7409655

從上面的結(jié)果可以看出，無論是否應用了超采樣辩蛋，邏輯回歸呻畸、LDA、QDA和NB在訓練和測試數(shù)據(jù)上的表現(xiàn)非常相似悼院。另一方面伤为，k-NN、CART和隨機森林在訓練數(shù)據(jù)上的預測效果要好得多，這表明過度擬合绞愚。

此外叙甸，過度取樣使所有學習器的AUC性能幾乎沒有變化。

下面的箱線圖展示了所有學習器的5折交叉驗證的AUC性能位衩。

# boxplot of AUC values across the 5 folds
autoplot(bmr2, measure = msr("classif.auc"))

[圖片上傳失敗...(image-175248-1649651130723)]

autoplot(bmr2,type = "roc")+
  scale_color_discrete() +
  theme_bw()

[圖片上傳失敗...(image-f107be-1649651130723)]

隨后裆蒸，輸出每個學習器的敏感性、特異性糖驴、假陰性率（FNR）和假陽性率（FPR）僚祷。

tab2 = bmr2$aggregate(msrs(c('classif.auc', 'classif.sensitivity','classif.specificity',
                            'classif.fnr',  'classif.fpr')))
tab2 = tab2[,c('learner_id','classif.auc','classif.sensitivity','classif.specificity',
               'classif.fnr',   'classif.fpr')]
print(tab2)

> print(tab2)
                               learner_id classif.auc classif.sensitivity
 1:                       classif.log_reg   0.7485372           0.8917097
 2:                           classif.lda   0.7487159           0.9037005
 3:                           classif.qda   0.7441634           0.6779116
 4:                   classif.naive_bayes   0.7498427           0.6250430
 5:                    classif.kknn.tuned   0.7180074           0.8509180
 6:                   classif.rpart.tuned   0.6987046           0.8679289
 7:                  classif.ranger.tuned   0.7506405           0.9447504
 8:      scale.oversample.classif.log_reg   0.7475678           0.6008893
 9:          scale.oversample.classif.lda   0.7489090           0.5841652
10:          scale.oversample.classif.qda   0.7431096           0.5529547
11:  scale.oversample.classif.naive_bayes   0.7494055           0.5505164
12:   scale.oversample.classif.kknn.tuned   0.6924480           0.6948078
13:  scale.oversample.classif.rpart.tuned   0.6753005           0.7090075
14: scale.oversample.classif.ranger.tuned   0.7393948           0.7427424
    classif.specificity classif.fnr classif.fpr
 1:           0.2516934  0.10829030   0.7483066
 2:           0.1855615  0.09629948   0.8144385
 3:           0.6946524  0.32208835   0.3053476
 4:           0.7488414  0.37495697   0.2511586
 5:           0.2581105  0.14908204   0.7418895
 6:           0.3108734  0.13207114   0.6891266
 7:           0.1554367  0.05524957   0.8445633
 8:           0.7663102  0.39911073   0.2336898
 9:           0.8023173  0.41583477   0.1976827
10:           0.8139037  0.44704532   0.1860963
11:           0.8381462  0.44948365   0.1618538
12:           0.5811052  0.30519220   0.4188948
13:           0.5449198  0.29099254   0.4550802
14:           0.5509804  0.25725760   0.4490196

事實證明，在沒有超采樣的情況下贮缕，邏輯回歸辙谜、LDA、k-NN跷睦、CART和隨機森林在敏感性方面得分很高筷弦，而在特異性方面得分相當?shù)停涣硪环矫嬉种睿琎DA和天真貝葉斯在特異性方面得分相對較高烂琴，但在敏感性方面卻沒有那么高。根據(jù)定義蜕乡，高靈敏度（特異性）源于低的假陰性（陽性）率奸绷，這在數(shù)據(jù)中也有體現(xiàn)。

提取單個模型

## 提取隨機森林模型
bmr_rf = bmr2$clone(deep = TRUE)$filter(learner_ids = 'classif.ranger.tuned')
## ROC
autoplot(bmr_rf,type = "roc")+
  scale_color_discrete() +
  theme_bw()
## PRC
autoplot(bmr_rf, type = "prc")+
  scale_color_discrete() +
  theme_bw()

[圖片上傳失敗...(image-4c97eb-1649651130723)]

[圖片上傳失敗...(image-a2ecaa-1649651130723)]

關(guān)于哪種學習器效果最好层玲，也包括是否應該使用超量取樣号醉，在很大程度上取決于敏感性和特異性的現(xiàn)實意義。就實際的重要性而言辛块，兩者中的一個可能會超過另一個很多倍畔派。想想典型的HIV快速診斷測試的例子，以低特異性為代價的高靈敏度可能會引起（不必要的）震驚润绵，但除此之外并不危險线椰，而低靈敏度則是非常危險的。正如通常的情況一樣尘盼，這里不存在黑白分明的 "最佳模型"憨愉。回顧一下卿捎，即使有超額取樣配紫，我們的模型沒有一個在靈敏度和特異性方面表現(xiàn)良好。在我們的案例中午阵，我們需要思考：以低靈敏度為代價的高特異性的后果是什么躺孝，這意味著告訴許多肝病患者他們是健康的；而以低特異性為代價的高靈敏度的后果是什么，這意味著告訴許多健康患者他們有肝病括细。在沒有進一步的特定主題信息的情況下伪很，我們只能說明在所選擇的特定性能指標上表現(xiàn)最好的學習器。如上所述奋单，基于AUC的隨機森林表現(xiàn)最好。此外猫十，隨機森林是靈敏度得分最高（FNR最低）的學習器览濒，而樸素貝葉斯是特異性最好（FPR最低）的學習器。

然而拖云，我們進行的分析決不是詳盡的贷笛。在特征層面上，雖然我們在分析過程中幾乎只關(guān)注了機器學習和統(tǒng)計分析方面宙项，但也可以更深入地挖掘?qū)嶋H的主題（肝卜唷），并嘗試更徹底地理解變量以及潛在的相關(guān)性和互動性尤筐。這可能也意味著要再次考慮已經(jīng)刪除的變量汇荐。此外，可以對數(shù)據(jù)集進行特征工程和數(shù)據(jù)預處理盆繁，例如使用主成分分析掀淘。關(guān)于超參數(shù)的調(diào)整，可以考慮使用更大的超參數(shù)空間和評估數(shù)量的不同超參數(shù)油昂。此外革娄，調(diào)整也可以應用于那些被我們標記為基線學習者的一些學習器。最后冕碟，還有更多的分類器存在拦惋，特別是梯度提升和支持向量機可以另外應用于這項任務，并有可能產(chǎn)生更好的結(jié)果安寺。

參考

• （mlr3gallery: Liver Patient Classification Based on Diagnostic Measures ）(https://mlr3gallery.mlr-org.com/posts/2020-09-11-liver-patient-classification/)