kaggle專案:基於隨機森林模型的心臟病患者預測分類!

語言: CN / TW / HK

公眾號:尤而小屋
作者:Peter
編輯:Peter

大家好,我是Peter~

新年的第一個專案實踐 ~給大家分享一個新的kaggle案例: 基於隨機森林模型(RandomForest)的心臟病人預測分類 。本文涉及到的知識點主要包含:

  • 資料預處理和型別轉化

  • 隨機森林模型建立與解釋

  • 決策樹如何視覺化

  • 基於混淆矩陣的分類評價指標

  • 部分依賴圖PDP的繪製和解釋

  • AutoML機器學習SHAP庫的使用和解釋(待提升)

導讀

Of all the applications of machine-learning, diagnosing any serious disease using a black box is always going to be a hard sell. If the output from a model is the particular course of treatment (potentially with side-effects), or surgery, or the absence of treatment, people are going to want to know why .

在機器學習的所有應用中,使用黑匣子診斷任何嚴重疾病總是很難的。如果模型的輸出是特定的治療過程(可能有副作用)、手術或是否有療效,人們會想知道為什麼。

This dataset gives a number of variables along with a target condition of having or not having heart disease. Below, the data is first used in a simple random forest model, and then the model is investigated using ML explainability tools and techniques.

該資料集提供了許多變數以及患有或不患有心臟病的目標條件。下面,資料首先用於一個簡單的隨機森林模型,然後使用 ML 可解釋性工具和技術對該模型進行研究。

感興趣的可以參考原文學習,notebook地址為:https://www.kaggle.com/tentotheminus9/what-causes-heart-disease-explaining-the-model

匯入庫

本案例中涉及到多個不同方向的庫:

  • 資料預處理

  • 多種視覺化繪圖;尤其是 shap的視覺化 ,模型可解釋性的使用(後面會專門寫這個庫)

  • 隨機森林模型

  • 模型評價等

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns 
from sklearn.ensemble import RandomForestClassifier 
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_graphviz 
from sklearn.metrics import roc_curve, auc 
from sklearn.metrics import classification_report 
from sklearn.metrics import confusion_matrix 
from sklearn.model_selection import train_test_split 
import eli5 
from eli5.sklearn import PermutationImportance
import shap 
from pdpbox import pdp, info_plots 
np.random.seed(123) 

pd.options.mode.chained_assignment = None  

資料探索EDA

1、匯入資料

2、缺失值情況

資料比較完美,沒有任何缺失值!

欄位含義

在這裡重點介紹下各個欄位的含義。Peter近期匯出的資料集中的額欄位和原notebook中的欄位名字寫法稍有差異(時間原因導致),還好Peter已經為大家做了一一對應的關係,下面是具體的中文含義:

  1. age:年齡

  2. sex 性別 1=male  0=female

  3. cp  胸痛型別;4種取值情況

  • 1:典型心絞痛

  • 2:非典型心絞痛

  • 3:非心絞痛

  • 4:無症狀

  • trestbps 靜息血壓

  • chol 血清膽固醇

  • fbs 空腹血糖 >120mg/dl  :1=true;0=false

  • restecg 靜息心電圖(值0,1,2)

  • thalach 達到的最大心率

  • exang 運動誘發的心絞痛(1=yes;0=no)

  • oldpeak 相對於休息的運動引起的ST值(ST值與心電圖上的位置有關)

  • slope 運動高峰ST段的坡度

    • 1:upsloping向上傾斜

    • 2:flat持平

    • 3:downsloping向下傾斜

  • ca  The number of major vessels(血管) (0-3)

  • thal A blood disorder called thalassemia ,一種叫做地中海貧血的血液疾病(3 = normal;6 = fixed defect;;7 = reversable defect)

  • target 生病沒有(0=no;1=yes)

  • 原notebook中的英文含義;

    下面是Peter整理的對應關係。本文中以當前的版本為標準:

    欄位轉化

    轉化編碼

    對部分欄位進行一一的轉化。以sex欄位為例:將資料中的0變成female,1變成male

    # 1、sex
    
    df["sex"][df["sex"] == 0] = "female"
    df["sex"][df["sex"] == 1] = "male"
    

    欄位型別轉化

    # 指定資料型別
    df["sex"] = df["sex"].astype("object")
    df["cp"] = df["cp"].astype("object")
    df["fbs"] = df["fbs"].astype("object")
    df["restecg"] = df["restecg"].astype("object")
    df["exang"] = df["exang"].astype("object")
    df["slope"] = df["slope"].astype("object")
    df["thal"] = df["thal"].astype("object")
    

    生成啞變數

    # 生成啞變數
    df = pd.get_dummies(df,drop_first=True)
    df
    

    隨機森林RandomForest

    切分資料

    # 生成特徵變數資料集和因變數資料集
    X = df.drop("target",1)
    y = df["target"]
    
    # 切分比例為8:2
    X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2,random_state=10)
    
    X_train
    

    建模

    rf = RandomForestClassifier(max_depth=5)
    rf.fit(X_train, y_train)
    

    3個重要屬性

    隨機森林中3個重要的屬性:

    • 檢視森林中樹的狀況 :estimators_

    • 袋外估計準確率得分:oob_score_,必須是 oob_score 引數選擇True的時候才可用
    • 變數的重要性 :feature_importances_

    決策樹視覺化

    在這裡我們選擇的第二棵樹的視覺化過程:

    # 檢視第二棵樹的狀況
    estimator = rf.estimators_[1]
    
    # 全部屬性
    feature_names = [i for i in X_train.columns]
    #print(feature_names)
    
    # 指定資料型別
    y_train_str = y_train.astype('str')
    # 0-no 1-disease
    y_train_str[y_train_str == '0'] = 'no disease'
    y_train_str[y_train_str == '1'] = 'disease'
    # 訓練資料的取值
    y_train_str = y_train_str.values
    y_train_str[:5]
    

    繪圖的具體程式碼為:

    # 繪圖顯示
    
    export_graphviz(
        estimator,   # 傳入第二顆樹
        out_file='tree.dot',   # 匯出檔名
        feature_names = feature_names,  # 屬性名
        class_names = y_train_str,  # 最終的分類資料
        rounded = True, 
        proportion = True, 
        label='root',
        precision = 2, 
        filled = True
    )
    
    from subprocess import call
    call(['dot', '-Tpng', 'tree.dot', '-o', 'tree.png', '-Gdpi=600'])
    
    from IPython.display import Image
    Image(filename = 'tree.png')
    

    決策樹的視覺化過程能夠讓我們看到具體的分類過程,但是並不能解決哪些特徵或者屬性比較重要。後面會對部分屬性的特徵重要性進行探索

    模型得分驗證

    關於混淆矩陣和使用特異性(specificity)以及靈敏度(sensitivity)這兩個指標來描述分類器的效能:

    # 模型預測
    y_predict = rf.predict(X_test)
    y_pred_quant = rf.predict_proba(X_test)[:,1]
    y_pred_bin = rf.predict(X_test)
    
    # 混淆矩陣
    confusion_matrix = confusion_matrix(y_test,y_pred_bin)
    confusion_matrix
    
    # 計算sensitivity and specificity 
    total=sum(sum(confusion_matrix))
    sensitivity = confusion_matrix[0,0]/(confusion_matrix[0,0]+confusion_matrix[1,0])
    specificity = confusion_matrix[1,1]/(confusion_matrix[1,1]+confusion_matrix[0,1])
    

    繪製ROC曲線

    fpr, tpr, thresholds = roc_curve(y_test, y_pred_quant)
    
    fig, ax = plt.subplots()
    ax.plot(fpr, tpr)
    
    ax.plot([0,1],[0,1],
            transform = ax.transAxes,
            ls = "--",
            c = ".3"
           )
    
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.0])
    
    plt.rcParams['font.size'] = 12
    
    # 標題
    plt.title('ROC Curve')
    # 兩個軸的名稱
    plt.xlabel('False Positive Rate (1 - Specificity)')
    plt.ylabel('True Positive Rate (Sensitivity)')
    # 網格線
    plt.grid(True)
    

    本案例中的ROC曲線值:

    auc(fpr, tpr)
    # 結果
    0.9076923076923078
    

    根據一般ROC曲線的評價標準,案例的表現結果還是不錯的:

    • 0.90 - 1.00 = excellent

    • 0.80 - 0.90 = good

    • 0.70 - 0.80 = fair

    • 0.60 - 0.70 = poor

    • 0.50 - 0.60 = fail

    補充知識點:分類器的評價指標

    考慮一個二分類的情況,類別為1和0,我們將1和0分別作為正類(positive)和負類(negative),根據實際的結果和預測的結果,則最終的結果有4種,表格如下:

    常見的評價指標:

    1、ACC:classification accuracy,描述分類器的分類準確率

    計算公式為:ACC=(TP+TN)/(TP+FP+FN+TN)

    2、BER:balanced error rate 計算公式為:BER=1/2*(FPR+FN/(FN+TP))

    3、TPR:true positive rate,描述識別出的所有正例佔所有正例的比例 計算公式為:TPR=TP/ (TP+ FN)

    4、FPR:false positive rate,描述將負例識別為正例的情況佔所有負例的比例 計算公式為:FPR= FP / (FP + TN)

    5、TNR:true negative rate,描述識別出的負例佔所有負例的比例 計算公式為:TNR= TN / (FP + TN)

    6、PPV:Positive predictive value計算公式為:PPV=TP / (TP + FP)

    7、NPV:Negative predictive value計算公式:NPV=TN / (FN + TN)

    其中 TPR 即為敏感度(sensitivity), TNR 即為特異度(specificity)。

    來自維基百科的經典圖形:

    可解釋性

    排列重要性-Permutation Importance

    下面的內容是關於機器學習模型的結果可解釋性。首先考察的是每個變數對模型的重要性。重點考量的排列重要性Permutation Importance:

    部分依賴圖( Partial dependence plots ,PDP)

    一維PDP

    Partial Dependence就是用來解釋某個特徵和目標值y的關係的,一般是通過畫出Partial Dependence Plot(PDP)來體現。也就是說PDP在X1的值,就是把訓練集中第一個變數換成X1之後,原模型預測出來的平均值。

    重點:檢視單個特徵和目標值的關係

    欄位ca

    base_features = df.columns.values.tolist()
    base_features.remove("target")
    
    feat_name = 'ca'  # ca-num_major_vessels 原文
    pdp_dist = pdp.pdp_isolate(
        model=rf,  # 模型
        dataset=X_test,  # 測試集
        model_features=base_features,  # 特徵變數;除去目標值 
        feature=feat_name  # 指定單個欄位
    )
    
    pdp.pdp_plot(pdp_dist, feat_name)  # 傳入兩個引數
    plt.show()
    

    通過下面的圖形我們觀察到:當ca欄位增加的時候,患病的機率在下降。ca欄位的含義是血管數量(num_major_vessels),也就是說當血管數量增加的時候,患病率隨之降低

    欄位age

    feat_name = 'age'
    
    pdp_dist = pdp.pdp_isolate(
        model=rf, 
        dataset=X_test, 
        model_features=base_features, 
        feature=feat_name)
    
    pdp.pdp_plot(pdp_dist, feat_name)
    plt.show()
    

    關於年齡欄位age,原文的描述:

    That's a bit odd. The higher the age, the lower the chance of heart disease? Althought the blue confidence regions show that this might not be true (the red baseline is within the blue zone).
    

    翻譯:這有點奇怪。年齡越大,患心臟病的機率越低?儘管藍色置信區間表明這可能不是真的(紅色基線在藍色區域內)

    欄位oldpeak

    feat_name = 'oldpeak'
    
    pdp_dist = pdp.pdp_isolate(
        model=rf, 
        dataset=X_test, 
        model_features=base_features, 
        feature=feat_name)
    
    pdp.pdp_plot(pdp_dist, feat_name)
    plt.show()
    

    oldpeak欄位同樣表明:取值越大,患病機率越低。

    這個變數稱之為“相對休息運動引起的ST壓低值”。正常的狀態下, 該值越高,患病機率越高 。但是上面的影象卻顯示了相反的結果。

    作者推斷:造成這個結果的原因除了the depression amount,可能還和slope type有關係。原文摘錄如下,於是作者繪製了2D-PDP圖形

    Perhaps it's not just the depression amount that's important, but the interaction with the slope type? Let's check with a 2D PDP
    

    2D-PDP圖

    檢視的是 slope_upsloping 、slope_flat和 oldpeak的關係:

    inter1  =  pdp.pdp_interact(
        model=rf,  # 模型
        dataset=X_test,  # 特徵資料集
        model_features=base_features,  # 特徵
        features=['slope_upsloping', 'oldpeak'])
    
    pdp.pdp_interact_plot(
        pdp_interact_out=inter1, 
        feature_names=['slope_upsloping', 'oldpeak'], 
        plot_type='contour')
    plt.show()
    
    ## ------------
    
    inter1  =  pdp.pdp_interact(
        model=rf, 
        dataset=X_test,
        model_features=base_features, 
        features=['slope_flat', 'oldpeak']
    )
    
    pdp.pdp_interact_plot(
        pdp_interact_out=inter1, 
        feature_names=['slope_flat', 'oldpeak'], 
        plot_type='contour')
    plt.show()
    

    從兩張圖形中我們可以觀察到:在oldpeak取值較低的時候,患病機率都比較高(黃色),這是一個奇怪的現象。於是作者進行了如下的SHAP視覺化探索:針對單個變數進行分析。

    SHAP視覺化

    關於SHAP的介紹可以參考文章:https://zhuanlan.zhihu.com/p/83412330    和   https://blog.csdn.net/sinat_26917383/article/details/115400327

    SHAP是Python開發的一個"模型解釋"包,可以解釋任何機器學習模型的輸出。下面SHAP使用的部分功能:

    Explainer

    在SHAP中進行模型解釋之前需要先建立一個 explainer ,SHAP支援很多型別的explainer,例如deep, gradient, kernel, linear, tree, sampling)。在這個案例我們以tree為例:

    # 傳入隨機森林模型rf
    explainer = shap.TreeExplainer(rf)  
    # 在explainer中傳入特徵值的資料,計算shap值
    shap_values = explainer.shap_values(X_test)  
    shap_values
    
    image-20220131173928357

    Feature Importance

    取每個特徵 SHAP值的絕對值的平均值 作為該特徵的重要性,得到一個標準的條形圖(multi-class則生成堆疊的條形圖:

    結論:能夠直觀地觀察到ca欄位的SHAP值最高

    summary_plot

    summary plot 為每個樣本繪製其每個特徵的SHAP值,這可以更好地理解整體模式,並允許發現預測異常值。

    • 每一行代表一個特徵,橫座標為SHAP值

    • 一個點代表一個樣本,顏色表示特徵值的高低(紅色高,藍色低)

    個體差異

    檢視單個病人的不同特徵屬性對其結果的影響,原文描述為:

    Next, let's pick out individual patients and see how the different variables are affecting their outcomes
    
    def heart_disease_risk_factors(model, patient):
    
        explainer = shap.TreeExplainer(model)  # 建立explainer
        shap_values = explainer.shap_values(patient)  # 計算shape值
        shap.initjs()
        return shap.force_plot(
          explainer.expected_value[1],
          shap_values[1], 
          patient)
    

    從兩個病人的結果中顯示:

    • P1:預測準確率高達29%(baseline是57%),更多的因素集中在ca、thal_fixed_defect、oldpeak等藍色部分。

    • P3:預測準確率高達82%,更多的影響因素在sel_male=0,thalach=143等

    通過對比不同的患者,我們是可以觀察到不同病人之間的預測率和主要影響因素。

    dependence_plot

    為了理解單個feature如何影響模型的輸出,我們可以將該feature的SHAP值與資料集中所有樣本的feature值進行比較:

    多樣本視覺化探索

    下面的圖形是針對多患者的預測和影響因素的探索。在jupyter notebook的互動式作用下,能夠觀察不同的特徵屬性對前50個患者的影響:

    shap_values = explainer.shap_values(X_train.iloc[:50])
    shap.force_plot(explainer.expected_value[1], 
                    shap_values[1], 
                    X_test.iloc[:50])
    

    資料集獲取方式 :公眾號 後臺 ,回覆  HeartDisease   即可領取~

    乾貨學習, 三連