本文是一個課程報告，由我和另外一位同學(xué)合作完成。自我感覺做的還行決定放上來。

?數(shù)據(jù)集來源：Cardiovascular Study Dataset | Kaggle

目錄

1.項目背景... 3

1.1項目說明... 3

1.2需求分析... 3

2.數(shù)據(jù)挖掘準(zhǔn)備... 3

2.1數(shù)據(jù)字段含義介紹... 3

2.2基礎(chǔ)統(tǒng)計分析... 4

3.數(shù)據(jù)挖掘過程... 5

3.1數(shù)據(jù)預(yù)處理... 5

3.1.1文字型變量數(shù)值化... 5

3.1.2缺失值處理... 6

3.1.3異常值處理... 8

3.1.4數(shù)據(jù)規(guī)范化... 10

3.2數(shù)據(jù)挖掘與可視化分析... 10

3.2.1人口統(tǒng)計信息分析... 11

3.2.2疾病史與亞健康狀態(tài)分析... 13

3.2.3重要影響指標(biāo)分析... 18

3.3特征工程... 19

3.3.1創(chuàng)建新特征... 19

3.3.2計算特征相關(guān)性... 20

3.3.3降維... 22

3.4建模??? 23

4.結(jié)果展示與評價... 25

4.1可視化結(jié)果展示... 25

4.2 模型結(jié)果展示... 26

5.總結(jié)... 28

5.1問題及解決方案... 28

5.2設(shè)計方案的優(yōu)點及不足... 28

5.2.1優(yōu)點... 28

5.2.2不足... 28

6.參考資料... 29

附錄1：完整代碼... 30

1.項目背景

1.1項目說明

冠心病，即冠狀動脈性心臟病（cardio vascular disease），是一種常見的心臟病，也是奪走眾多生命的重大疾病之一，世界衛(wèi)生組織統(tǒng)計，一半的死亡是由心血管疾病造成的。截止至2021年，推算的中國冠心病現(xiàn)患人數(shù)為1139萬[3]，且中國CVD患病率處于持續(xù)上升階段。患有高血壓、糖尿病等疾病，以及過度肥胖、不良生活習(xí)慣等是誘發(fā)該病的主要因素，以老年人為主要發(fā)病群體。

提前監(jiān)測身體健康指標(biāo)并預(yù)測心血管疾病的發(fā)作，可以最大程度上減少CVD的患病與傷亡數(shù)；同時，對已患病人群盡早干預(yù)，使其回歸正確生活方式，可以幫助高?；颊邷p少并發(fā)癥。

綜上所述，本研究旨在通過數(shù)據(jù)可視化與數(shù)據(jù)挖掘的方法，找出與冠心病病最相關(guān)的因素，探索數(shù)據(jù)之間潛在的聯(lián)系，使用邏輯回歸預(yù)測患病風(fēng)險，來幫助相關(guān)人員預(yù)防并控制冠心病。

1.2需求分析

通過數(shù)據(jù)可視化分析與數(shù)據(jù)挖掘預(yù)測，深層探索數(shù)據(jù)之間的潛在聯(lián)系，為冠心病的檢測與預(yù)防發(fā)展提供依據(jù)。一方面分析個人的基本信息、生活習(xí)慣、病史等和冠心病風(fēng)險的關(guān)系，判斷哪些人群是患病高危人群，需要重點關(guān)注預(yù)防，并且如果和生活習(xí)慣有關(guān)，可以借此提出習(xí)慣改善建議；另一方面分析臨床指標(biāo)和冠心病風(fēng)險的關(guān)系，用于輔助冠心病的臨床診斷。再結(jié)合兩方面的分析，對某個特定的人的冠心病風(fēng)險進(jìn)行預(yù)測，實現(xiàn)“早發(fā)現(xiàn)、早預(yù)防”。

2.數(shù)據(jù)挖掘準(zhǔn)備

2.1數(shù)據(jù)字段含義介紹

數(shù)據(jù)集在Kaggle網(wǎng)站上公開發(fā)布，它來自馬薩諸塞州，弗雷明漢鎮(zhèn)的研究人員正在進(jìn)行的心血管研究。數(shù)據(jù)集共有17列，其名稱與含義如表格1：

表格 1 數(shù)據(jù)列語義信息

2.2基礎(chǔ)統(tǒng)計分析

該數(shù)據(jù)集一共有3390條數(shù)據(jù)，16個特征，一個預(yù)測目標(biāo)。其中，id字段用于作為每條數(shù)據(jù)的唯一標(biāo)識，并無現(xiàn)實意義；sex和is_smoking兩個字段為文字型變量，其余字段均為數(shù)值型變量。

對數(shù)據(jù)集進(jìn)行缺失情況統(tǒng)計，統(tǒng)計結(jié)果如下表所示

表格 2 缺失值統(tǒng)計結(jié)果

glucose字段缺失值較多，缺失條目占到了總條目的近10%，需要進(jìn)行缺失值處理。各列的數(shù)據(jù)統(tǒng)計情況分析如表格2與表格3。

表格 3 連續(xù)型數(shù)據(jù)統(tǒng)計情況

表格 4 離散型數(shù)據(jù)統(tǒng)計情況

3.數(shù)據(jù)挖掘過程

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

3.1.1文字型變量數(shù)值化

將sex，is_smoking兩個字段數(shù)值化：sex中F設(shè)為1，M設(shè)為-1，使得兩者距離原點對稱，對于原點具有相同的重要性；is_smoking中YES設(shè)為1，NO設(shè)為0。代碼及結(jié)果如下

def numerical(data):
    """
    將文本變量進(jìn)行編碼
    """
    data.loc[data.sex == 'F', 'sex'] = 1
    data.loc[data.sex == 'M', 'sex'] = -1
    data.loc[data.is_smoking == 'YES', 'is_smoking'] = 1
    data.loc[data.is_smoking == 'NO', 'is_smoking'] = 0
    return data

3.1.2缺失值處理

考慮到葡萄糖水平和是否患糖尿病可能有較強的相關(guān)性，計算葡萄糖水平和糖尿病的皮爾遜相關(guān)系數(shù)，并和其它指標(biāo)進(jìn)行比較，驗證猜想

 corr = data['glucose'].corr(data['diabetes'])
    print('糖尿病和葡萄糖水平相關(guān)系數(shù)：', corr)
    corr = data['glucose'].corr(data['BMI'])
    print('肥胖程度和葡萄糖水平相關(guān)系數(shù)：', corr)
    corr = data['glucose'].corr(data['prevalentHyp'])
    print('高血壓和葡萄糖水平相關(guān)系數(shù)：', corr)

可以看到，是否患有糖尿病與葡萄糖水平具有較高的關(guān)聯(lián)，可以將人群分為患病和不患病兩個類別，進(jìn)一步分析，作出兩者與葡萄糖關(guān)聯(lián)的箱線圖

圖 1 糖尿病與葡萄糖水平

可以看到，患有糖尿病人群的葡萄糖水平分布更廣，方差更大，均值更高，而未患病人群中葡萄糖水平的異常值更多，可能對均值產(chǎn)生較大影響，因此兩者均采用眾數(shù)進(jìn)行填補。

另外對缺失值較少的字段離散值用眾數(shù)填補，連續(xù)值用均值填補，代碼和填補結(jié)果如下

# 填充葡萄糖
    a = float(data[data.diabetes == 1]['glucose'].mode())
    b = float(data[data.diabetes == 0]['glucose'].mode())
    data1 = data.loc[data.diabetes == 1]
    data2 = data.loc[data.diabetes == 0]
    data1['glucose'].fillna(a, inplace=True)
    data2['glucose'].fillna(b, inplace=True)
    data = pd.concat([data1, data2])

    data['education'].fillna(int(data['education'].mode()), inplace=True)
    data['cigsPerDay'].fillna(int(data['cigsPerDay'].mode()), inplace=True)
    data['BPMeds'].fillna(int(data['BPMeds'].mode()), inplace=True)
    data['totChol'].fillna(data['totChol'].median(), inplace=True)
    data['BMI'].fillna(data['BMI'].median(), inplace=True)
    data['heartRate'].fillna(data['heartRate'].median(), inplace=True)

至此，缺失值全部處理完畢

?

3.1.3異常值處理

作出連續(xù)屬性的箱線圖

?

圖 2 連續(xù)屬性箱線圖

可以看到，除cigsPerDay外，每個屬性均有較多的值處于異常范圍，若將其直接刪掉，會丟失大量信息，影響數(shù)據(jù)分析的準(zhǔn)確度。但同時每個指標(biāo)的取值是否處于異常范圍或許與是否患病有關(guān)，因此在此處計算獲得每個連續(xù)指標(biāo)的非異常范圍，為后續(xù)特征分析做準(zhǔn)備，代碼和結(jié)果如下

def find_outlier(data):
    """
    畫箱線圖，找出異常范圍
    """
    constant = ['cigsPerDay', 'totChol', 'sysBP', 'diaBP', 'BMI', 'heartRate', 'glucose']

    plt.figure(figsize=(20, 15))
    for i in range(7):
        plt.subplot(2, 4, i + 1)
        bp = plt.boxplot(x=data[constant[i]])
        lower = [item.get_ydata()[1] for item in bp['whiskers']][0]
        upper = [item.get_ydata()[1] for item in bp['whiskers']][1]
        print('非異常范圍：', [lower, upper])

    # plt.show()

表格 5連續(xù)值正常范圍

并且經(jīng)過查詢可知，該方法得出的正常值范圍與一般醫(yī)學(xué)上認(rèn)定的正常值范圍較為接近，認(rèn)為該方法有效且該數(shù)據(jù)集正常。

3.1.4數(shù)據(jù)規(guī)范化

對臨床指標(biāo)進(jìn)行Z值規(guī)范化處理，以消除單位不統(tǒng)一帶來的量綱的影響，至此，數(shù)據(jù)預(yù)處理完成

3.2數(shù)據(jù)挖掘與可視化分析

聯(lián)系數(shù)據(jù)列含義與問題背景，將從4個維度對數(shù)據(jù)進(jìn)行分析。

?

3.2.1人口統(tǒng)計信息分析

繪制數(shù)據(jù)集的基本信息的面板圖，如圖1所示。其中，每一行代表著不同的教育水平，從上到下依次升序；每一列代表著不同的年齡，從最左列開始，以步長為5依次遞增；面板中每個條形圖表示著性別的分布情況，第一列為女性，第二列為男性；用深紅色的表示該部分人吸煙，淺黃色則為不吸煙；最后，用面板的顏色深淺表示該部分人的平均患病風(fēng)險情況。

?

圖 3 人口統(tǒng)計數(shù)據(jù)面板圖

可以看到，（1）數(shù)據(jù)的來源以中低教育水平、女性居多；（2）抽煙人群，以低教育水平，中年人居多，性別之間基本持平；（3）抽煙比例隨年齡的上升，呈下降趨勢，以女性尤為明顯；（4）患病風(fēng)險隨年齡上升而顯著增加。需要注意的是，在年齡分布的邊緣區(qū)域偶見反常值，為樣本數(shù)據(jù)過少而不可避免的偶然性，不應(yīng)作為參考。

進(jìn)一步，繪制抽煙比例與年齡之間的堆積百分比條形圖，如圖2所示，其中，左欄為女性，右欄為男性：

?

圖 4 抽煙比例-年齡百分比堆積圖

可以看到，隨著年齡的增加，有抽煙習(xí)慣的人顯著減少，女性部分的變化大于男性，與面板圖情況一致。

3.2.2疾病史與亞健康狀態(tài)分析

冠心病的發(fā)作與過往疾病史和身體的亞健康狀態(tài)息息相關(guān)，通過分別繪制有無患病風(fēng)險的人群的數(shù)據(jù)分布情況，觀察其影響因素與大小。

繪制疾病記錄與年齡桑基圖，圖5圖6分別是無風(fēng)險與有風(fēng)險人群。其中，左側(cè)的堆積條形圖表示疾病史的分類情況，右側(cè)的堆積條形圖表示年齡的堆積情況，中間用曲線表示兩部分的聯(lián)系，曲線的粗細(xì)表明記錄數(shù)的多少。

圖 5 無風(fēng)險人群患病與年齡?；鶊D

?圖 6有風(fēng)險人群患病與年齡?；鶊D

可以看到，（1）兩類人中有過高血壓經(jīng)歷的人均明顯占多，數(shù)據(jù)呈有偏分布；（2）兩類人中，年齡越大，患病的經(jīng)歷就越多。不同的是，（3）有風(fēng)險人群有中風(fēng)和糖尿病的比例更多些；（4）有風(fēng)險人群中老年人群所占的比例更大些。這一點與常識相符。

桑基圖在邊緣地區(qū)產(chǎn)生畸變，同樣是因為樣本偶然性原因，不應(yīng)考慮。

值得注意的是，在兩類圖中，年齡以中間多兩邊少的趨勢呈近似正態(tài)分布，繪制年齡與患病記錄直方圖如圖7。

?圖 7年齡-患病記錄直方圖

??? 可以得到結(jié)果：（1）列1中，整體呈左偏分布，說明年齡越大，沒有患病記錄的人就越少了，這兩個因素之間為正相關(guān)；（2）列2及列3，整體呈右偏分布，進(jìn)一步證實相關(guān)性；（3）且隨著患病記錄的增加，有風(fēng)險人群的占比越來越大，說明患病會增加冠心病發(fā)作的概率。

各種指標(biāo)是身體狀態(tài)的直接表示，考慮有風(fēng)險與無風(fēng)險之間的指標(biāo)差異對比，繪制箱型圖如圖8。

?

?圖 8 有無風(fēng)險人群各身體指標(biāo)箱型圖

淺綠色的列為無風(fēng)險，橙色的列為有風(fēng)險?？梢钥吹?，差距僅僅是有風(fēng)險列的方差過大。說明冠心病是由多種因素并發(fā)導(dǎo)致的，風(fēng)險是否與單個指標(biāo)的關(guān)系并不明顯。

基于此，收集六大指標(biāo)在正常范圍之外的信息，并將其表注為亞健康記錄。

表格 6 亞健康記錄劃分標(biāo)準(zhǔn)

身體亞健康與不良習(xí)慣有關(guān)，即每日抽煙的根數(shù)。探索每日香煙根數(shù)的分布情況，患病記錄與亞健康記錄之間的關(guān)系，繪制香煙根數(shù)-亞健康記錄-患病記錄直方圖，如圖9。

?圖 9 香煙根數(shù)-亞健康記錄-患病記錄直方圖

綠色條狀圖為吸煙根數(shù)的分布，小部分吸煙人群每日香煙數(shù)集中在18-24根；紅色折線圖表示人群平均的患病記錄，黃色折線表示人群平均的亞健康記錄。兩條折線之間僅有很小的相似性，折線與香煙根數(shù)之間的相關(guān)性也有限，判斷該分析以香煙為維度結(jié)果不顯著。以香煙為不良習(xí)慣的代表有很大局限性。

但患病記錄與亞健康記錄之間的關(guān)系確實存在，繪制瀑布圖，如圖10。

圖 10 亞健康記錄-患病記錄瀑布圖

每有一個指標(biāo)為亞健康，患病記錄的增加成指數(shù)型增長。

綜上所述，單個身體健康指標(biāo)對患冠心病風(fēng)險的作用不大，但多個健康指標(biāo)會綜合影響患病記錄，患病記錄再增加冠心病風(fēng)險的可能。

3.2.3重要影響指標(biāo)分析

冠心病為經(jīng)典的老年病癥，且與肥胖程度高度相關(guān)，綜合考慮BMI與年齡對冠心病風(fēng)險的影響，畫體重-年齡矩陣密度圖，如圖11，其中，淺綠色的點為無風(fēng)險，橙色的點為有風(fēng)險，顏色深淺表示落在這個坐標(biāo)的樣本數(shù)的多少。

圖 11 體重-年齡矩陣密度圖

??? 與無風(fēng)險人群相比，有風(fēng)險人群大多落在第一、第二、第四象限，有著高齡高體重的特征，且有風(fēng)險人群的BMI平均值大于無風(fēng)險人群。

3.3特征工程

3.3.1創(chuàng)建新特征

根據(jù)以上分析，在Python中創(chuàng)建“亞健康計數(shù)”特征

def create_sub_health(data):
    data['subHealth'] = [0]*data.shape[0]
    for i in range(data.shape[0]):
        count = 0
        if data['totChol'][i] < 240 or data['totChol'][i] > 560:
            count+=1
        if data['sysBP'][i] < 90 or data['sysBP'][i] > 140:
            count+=1
        if data['diaBP'][i] < 60 or data['diaBP'][i] > 90:
            count+=1
        if data['BMI'][i] < 18 or data['BMI'][i] > 24:
            count += 1
        if data['heartRate'][i] < 60 or data['heartRate'][i] > 100:
            count += 1
        if data['glucose'][i] < 39 or data['glucose'][i] > 61:
            count += 1
        data['subHealth'][i] = count

    return data

3.3.2計算特征相關(guān)性

計算每個特征的皮爾遜相關(guān)系數(shù)，畫出歸一化和未歸一化的熱力圖，進(jìn)行特征比較

?圖 12 熱力圖-未歸一化

?

圖 13 熱力圖-歸一化

不管是否歸一化，是否有十年冠心病風(fēng)險與單個特征的相關(guān)性都很小，且和年齡相關(guān)性最大，和tableau分析相符。且歸一化后削弱了各個特征之間的相關(guān)性，減弱了線性回歸中的多重共線性問題。但可以看到，是否患有糖尿病和收縮壓、亞健康記錄數(shù)和是否服用過降壓藥等特征之間仍有很強的相關(guān)性，因此進(jìn)行人工特征選擇降維或者進(jìn)行PCA降維。

3.3.3降維

根據(jù)以上分析，初步人工選擇特征為”age””prevalentHyp”“subHealth”三個特征。進(jìn)行PCA降維。先降為兩維，畫出聚類散點圖

def pca(data):
    """
    pca降維，作圖，找出各成分和主成分關(guān)系
    """
    from sklearn.decomposition import PCA
    data = data.iloc[:, 1:-1]
    x = np.array(data)
    pca = PCA(n_components=2)
    a = pca.fit_transform(x)
    print("降維后的數(shù)據(jù)為：", a)
    print("各主成分貢獻(xiàn)率為：", pca.explained_variance_ratio_)

    from sklearn.cluster import KMeans
    model = KMeans(n_clusters=2)
    model.fit(a)
    yhat = model.predict(a)  # 得出每個點所屬類別的預(yù)測
    label = np.unique(yhat)
    plt.rcParams['font.sans-serif'] = ["SimHei"]  # 使其能顯示中文
    plt.rcParams['axes.unicode_minus'] = False
    # print(label)
    for i in label:
        row_indx = np.where(yhat == i)
        plt.scatter(a[row_indx, 0], a[row_indx, 1], label="第" + str(i + 1) + "類", alpha=0.5, s=30)

    plt.ylabel('dim2')
    plt.xlabel('dim1')
    plt.legend()
    # plt.savefig()
    plt.show()

?

圖 14 聚類散點圖

當(dāng)只保留兩個主成分時，主成分累計貢獻(xiàn)率已經(jīng)達(dá)到了96%以上，可以只保留兩個主成分進(jìn)行后續(xù)建模。兩個主成分呈現(xiàn)出條帶形的關(guān)系，展現(xiàn)出一種既有序又難以用函數(shù)描述的關(guān)聯(lián)。對散點圖出現(xiàn)原因進(jìn)行推測，可能是原始維度中既有離散型變量又有連續(xù)型變量，連續(xù)變量使得降維后的散點圖呈條帶狀，離散型變量使得各個條帶分立開來。同時可以注意到，聚類結(jié)果中一個條帶上的點均屬于同一類別，說明離散特征是影響是否有患病風(fēng)險的重要特征，這也與之前的分析相符。

3.4建模

將上文中不同的數(shù)據(jù)處理方法和不同的機器學(xué)習(xí)方法相結(jié)合，建立12種組合模型，如下表所示

表格 7 組合模型參數(shù)

對隨機森林模型進(jìn)行調(diào)參

?

圖 15 隨機森林最大深度調(diào)整

可以看到最優(yōu)深度為2，模型再深時，訓(xùn)練集準(zhǔn)確率一直上升，測試集準(zhǔn)確率下降，模型過擬合。接著建立模型。部分代碼如下，完整代碼見附錄

def model1(data):
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model = model.fit(X_train, y_train)
    test_accuracy = model.score(X_test, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='Logistic Regression',color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy

4.結(jié)果展示與評價

4.1可視化結(jié)果展示

?? 通過上述可視化圖表，制作儀表盤與故事，呈現(xiàn)數(shù)據(jù)分析過程與結(jié)論。

?

圖 16 故事預(yù)覽圖

4.2 模型結(jié)果展示

可以得到以上十二個模型的ROC曲線和測試集上的預(yù)測準(zhǔn)確率如表格8和圖17。可以看到，最簡單的模型1，即不進(jìn)行任何數(shù)據(jù)處理，直接采用邏輯回歸進(jìn)行預(yù)測，在測試集上的準(zhǔn)確率是最高的。而在特征選擇上，采用人工特征選擇和PCA降維并無較大差異，兩者相互印證，說明年齡、亞健康記錄數(shù)、是否服用降壓藥確實與冠心病風(fēng)險密切相關(guān)。另外從ROC曲線來看，隨機森林的表現(xiàn)效果比邏輯回歸要更好，該模型可能在其它類似數(shù)據(jù)集上預(yù)測準(zhǔn)確率更穩(wěn)定。

表格 8 測試準(zhǔn)確率

def select_model(data):
    # 畫ROC曲線
    acclist = []
    plt.figure(figsize=(5, 15))
    plt.subplot(6, 2, 1)
    acclist.append(model1(train))
    plt.subplot(6, 2, 2)
    acclist.append(model2(train))
    plt.subplot(6, 2, 3)
    acclist.append(model3(train))
    plt.subplot(6, 2, 4)
    acclist.append(model4(train))
    plt.subplot(6, 2, 5)
    acclist.append(model5(train))
    plt.subplot(6, 2, 6)
    acclist.append(model6(train))
    plt.subplot(6, 2, 7)
    acclist.append(model7(train))
    plt.subplot(6, 2, 8)
    acclist.append(model8(train))
    plt.subplot(6, 2, 9)
    acclist.append(model9(train))
    plt.subplot(6, 2, 10)
    acclist.append(model10(train))
    plt.subplot(6, 2, 11)
    acclist.append(model11(train))
    plt.subplot(6, 2, 12)
    acclist.append(model12(train))
    plt.show()
    print(acclist)

?

圖 17 ROC曲線

5.總結(jié)

5.1問題及解決方案

1. 因人口信息維度的數(shù)量多，不知道從哪個視角開始分析，才不會遺漏重要的相關(guān)性。后綜合考慮，決定統(tǒng)一做一張面板圖：通過新建列與行的離散字段，獲得條形圖面板與面積圖面板，在儀表盤疊加，最終得到了有詳細(xì)信息的人口信息統(tǒng)計圖：人口信息的任意兩個維度在此圖上都能夠有機會得到呈現(xiàn)相關(guān)性的機會。
2. 在分析患病風(fēng)險與疾病史與亞健康狀態(tài)時，最開始因六大指標(biāo)與三大疾病史在有無風(fēng)險人群中的分布區(qū)別不夠明顯，找不到能呈現(xiàn)顯著結(jié)果分析維度，無法作數(shù)據(jù)分析，甚至懷疑冠心病是否沒有主要的致病因素。后利用疾病史與亞健康記錄的累加值新建特征，得到了較好的結(jié)果，發(fā)現(xiàn)單一指標(biāo)或維度不能導(dǎo)致患病風(fēng)險的提升，患病是所有因素綜合作用的結(jié)果。
3. 對缺失值進(jìn)行填充時，葡萄糖水平缺失過多，如果只是用均值填充感覺會引入大量噪聲。猜測葡萄糖水平應(yīng)該和是否患糖尿病有關(guān)聯(lián)，于是用是否患糖尿病將人群分類再填充，希望能填充的更準(zhǔn)確一些。
4. 在tableau可視化分析的基礎(chǔ)上做Python建模時，開始有些迷惑不知道該怎樣把兩個部分整合起來形成一個整體。后來用可視化分析出的“亞健康記錄”這個指標(biāo)作為突破口，在此基礎(chǔ)上進(jìn)行特征工程，將兩個部分結(jié)合起來

5.2設(shè)計方案的優(yōu)點及不足

5.2.1優(yōu)點

1. 首先通過探索性數(shù)據(jù)分析，觀察分析各類圖的可視化結(jié)果，得到了可能相關(guān)的初步因素與相關(guān)性結(jié)論。然后在通過數(shù)據(jù)挖掘的方法，用機器學(xué)習(xí)算法進(jìn)行降維、預(yù)測等操作，得到最終結(jié)論。
2. 采用多種分析工具相結(jié)合，能夠綜合各個工具的優(yōu)勢進(jìn)行分析，且得出的結(jié)論可以相互印證
3. 發(fā)現(xiàn)了“亞健康記錄數(shù)”這一重要指標(biāo)，將各個連續(xù)的特征聯(lián)系了起來。
4. 使用了多種模型進(jìn)行模型消融實驗，并采用測試集準(zhǔn)確率和ROC兩種方法進(jìn)行模型評價，能夠選出效果最好的模型

5.2.2不足

1. 特征之間的關(guān)聯(lián)挖掘不夠深入。還可以查閱相關(guān)文獻(xiàn)資料，找到各個特征間可能的相互關(guān)系，再進(jìn)行可視化分析進(jìn)行驗證分析。
2. 模型具有局限性?？梢钥吹皆跍y試集的預(yù)測準(zhǔn)確率上，各個模型相差并不大，可以采用更多的模型進(jìn)行嘗試，找到預(yù)測效果更好的模型。
3. BP Meds數(shù)據(jù)始終沒能用上，嘗試過多種數(shù)據(jù)維度分析，效果都不夠顯著。再有機會進(jìn)行深度挖掘時，該變量可能會派上用場

6.參考資料?

1. Cardiovascular Risk Prediction | Kaggle
2. 冠心病 - 醫(yī)學(xué)百科 (yixue.com)
3. 《中國心血管健康與疾病報告2021》要點解讀[J].中國心血管雜志,2022,27(04):305-318.

附錄1：完整代碼文章來源地址http://www.zghlxwxcb.cn/news/detail-801481.html

# main.py
import pandas as pd
import utils
import warnings
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, roc_curve
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
import numpy as np

warnings.filterwarnings("ignore")

train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')

pd.set_option('display.max_columns', None)


def model1(data):
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model = model.fit(X_train, y_train)
    test_accuracy = model.score(X_test, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='Logistic Regression',color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model2(data):
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    # utils.adjust_rfc(x, y)
    model = RandomForestClassifier(max_depth=2)
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='RandomForestClassifier', color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model3(data):
    y = data['TenYearCHD']
    x = data.loc[:, ['age', 'prevalentHyp', 'subHealth']]
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr,color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model4(data):
    y = data['TenYearCHD']
    x = data.loc[:, ['age', 'prevalentHyp', 'subHealth']]
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = RandomForestClassifier(max_depth=2)
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='RandomForestClassifier', color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model5(data):
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    x = np.array(x)
    pca = PCA(n_components=2)
    x = pca.fit_transform(x)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr,color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model6(data):
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    x = np.array(x)
    pca = PCA(n_components=2)
    x = pca.fit_transform(x)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = RandomForestClassifier(max_depth=2)
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr,color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model7(data):
    data = utils.standardize(data)
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='Logistic Regression',color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model8(data):
    data = utils.standardize(data)
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    # utils.adjust_rfc(x, y)
    model = RandomForestClassifier(max_depth=2)
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='RandomForestClassifier', color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model9(data):
    data = utils.standardize(data)
    y = data['TenYearCHD']
    x = data.loc[:, ['age', 'prevalentHyp', 'subHealth']]
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr,color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model10(data):
    data = utils.standardize(data)
    y = data['TenYearCHD']
    x = data.loc[:, ['age', 'prevalentHyp', 'subHealth']]
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = RandomForestClassifier(max_depth=2)
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr, label='RandomForestClassifier', color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model11(data):
    data = utils.standardize(data)
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    x = np.array(x)
    pca = PCA(n_components=2)
    x = pca.fit_transform(x)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = LogisticRegression()
    model.fit(X_train, y_train)
    test_class_preds = model.predict(X_test)
    test_accuracy = accuracy_score(test_class_preds, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr,color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def model12(data):
    data = utils.standardize(data)
    y = data['TenYearCHD']
    x = data.drop('TenYearCHD', axis=1)
    x = np.array(x)
    pca = PCA(n_components=2)
    x = pca.fit_transform(x)
    X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0)
    model = RandomForestClassifier(max_depth=2)
    model.fit(X_train, y_train)
    test_accuracy = model.score(X_test, y_test)
    y_lr_predict_pro = model.predict_proba(X_test)[:, 1]
    fpr, tpr, thresholds = roc_curve(y_test, y_lr_predict_pro)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.plot(fpr, tpr,color='r')
    plt.xlabel('TPR')
    plt.ylabel('FPR')
    return test_accuracy


def select_model(data):
    # 畫ROC曲線
    acclist = []
    plt.figure(figsize=(5, 15))
    plt.subplot(6, 2, 1)
    acclist.append(model1(train))
    plt.subplot(6, 2, 2)
    acclist.append(model2(train))
    plt.subplot(6, 2, 3)
    acclist.append(model3(train))
    plt.subplot(6, 2, 4)
    acclist.append(model4(train))
    plt.subplot(6, 2, 5)
    acclist.append(model5(train))
    plt.subplot(6, 2, 6)
    acclist.append(model6(train))
    plt.subplot(6, 2, 7)
    acclist.append(model7(train))
    plt.subplot(6, 2, 8)
    acclist.append(model8(train))
    plt.subplot(6, 2, 9)
    acclist.append(model9(train))
    plt.subplot(6, 2, 10)
    acclist.append(model10(train))
    plt.subplot(6, 2, 11)
    acclist.append(model11(train))
    plt.subplot(6, 2, 12)
    acclist.append(model12(train))
    plt.show()
    print(acclist)



if __name__ == '__main__':
    # utils.data_statics(train)
    train = utils.numerical(train)
    train = utils.fill_nan(train)
    # print(train.head())
    print(train.isnull().sum())
    # utils.find_outlier(train)
    train = utils.create_sub_health(train)
    # utils.draw_therm(train)
    # select_model(train)

# utils.py
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

def data_statics(data):
    """
    統(tǒng)計各字段的缺失值
    """
    print('缺失值統(tǒng)計：', data.isnull().sum())
    print(data.info())
    print('數(shù)值分布統(tǒng)計', data.describe())


def numerical(data):
    """
    將文本變量進(jìn)行編碼
    """
    data.loc[data.sex == 'F', 'sex'] = 1
    data.loc[data.sex == 'M', 'sex'] = -1
    data.loc[data.is_smoking == 'YES', 'is_smoking'] = 1
    data.loc[data.is_smoking == 'NO', 'is_smoking'] = 0
    return data


def fill_nan(data):
    """
    填充缺失值
    """
    # 計算其它指標(biāo)和葡萄糖水平的相關(guān)系數(shù)
    corr = data['glucose'].corr(data['diabetes'])
    print('糖尿病和葡萄糖水平相關(guān)系數(shù)：', corr)
    corr = data['glucose'].corr(data['BMI'])
    print('肥胖程度和葡萄糖水平相關(guān)系數(shù)：', corr)
    corr = data['glucose'].corr(data['prevalentHyp'])
    print('高血壓和葡萄糖水平相關(guān)系數(shù)：', corr)

    # 畫出葡萄糖
    # palette = sns.color_palette("BrBG", 2)
    # sns.boxplot(data=data, x='diabetes', y=data['glucose'], palette=palette)
    # # plt.show()

    # 填充葡萄糖
    a = float(data[data.diabetes == 1]['glucose'].mode())
    b = float(data[data.diabetes == 0]['glucose'].mode())
    data1 = data.loc[data.diabetes == 1]
    data2 = data.loc[data.diabetes == 0]
    data1['glucose'].fillna(a, inplace=True)
    data2['glucose'].fillna(b, inplace=True)
    data = pd.concat([data1, data2])

    data['education'].fillna(int(data['education'].mode()), inplace=True)
    data['cigsPerDay'].fillna(int(data['cigsPerDay'].mode()), inplace=True)
    data['BPMeds'].fillna(int(data['BPMeds'].mode()), inplace=True)
    data['totChol'].fillna(data['totChol'].median(), inplace=True)
    data['BMI'].fillna(data['BMI'].median(), inplace=True)
    data['heartRate'].fillna(data['heartRate'].median(), inplace=True)

    return data


def find_outlier(data):
    """
    畫箱線圖，找出異常范圍
    """
    constant = ['cigsPerDay', 'totChol', 'sysBP', 'diaBP', 'BMI', 'heartRate', 'glucose']

    plt.figure(figsize=(20, 15))
    for i in range(7):
        plt.subplot(2, 4, i + 1)
        bp = plt.boxplot(x=data[constant[i]])
        lower = [item.get_ydata()[1] for item in bp['whiskers']][0]
        upper = [item.get_ydata()[1] for item in bp['whiskers']][1]
        print('非異常范圍：', [lower, upper])

    # plt.show()


def standardize(data):
    from sklearn.preprocessing import StandardScaler  # 實現(xiàn)z-score標(biāo)準(zhǔn)化
    feature = ['totChol', 'sysBP', 'diaBP', 'BMI', 'heartRate', 'glucose']
    for i in range(6):
        x = np.array(data[feature[i]]).reshape(-1,1)
        ss = StandardScaler()
        x = ss.fit_transform(x)
        for j in range(data.shape[0]):
            data[feature[i]][j] = x[j][0]
    # print(data.head())
    return data


def create_sub_health(data):
    data['subHealth'] = [0]*data.shape[0]
    for i in range(data.shape[0]):
        count = 0
        if data['totChol'][i] < 240 or data['totChol'][i] > 560:
            count+=1
        if data['sysBP'][i] < 90 or data['sysBP'][i] > 140:
            count+=1
        if data['diaBP'][i] < 60 or data['diaBP'][i] > 90:
            count+=1
        if data['BMI'][i] < 18 or data['BMI'][i] > 24:
            count += 1
        if data['heartRate'][i] < 60 or data['heartRate'][i] > 100:
            count += 1
        if data['glucose'][i] < 39 or data['glucose'][i] > 61:
            count += 1
        data['subHealth'][i] = count

    return data


def draw_therm(data):
    """
    畫熱力圖
    """
    plt.figure()
    sns.heatmap(data.corr(), annot=True, center=0)
    plt.show()


def pca(data):
    """
    pca降維，作圖，找出各成分和主成分關(guān)系
    """
    from sklearn.decomposition import PCA
    data = data.iloc[:, 1:-1]
    x = np.array(data)
    pca = PCA(n_components=2)
    a = pca.fit_transform(x)
    print("降維后的數(shù)據(jù)為：", a)
    print("各主成分貢獻(xiàn)率為：", pca.explained_variance_ratio_)

    from sklearn.cluster import KMeans
    model = KMeans(n_clusters=2)
    model.fit(a)
    yhat = model.predict(a)  # 得出每個點所屬類別的預(yù)測
    label = np.unique(yhat)
    plt.rcParams['font.sans-serif'] = ["SimHei"]  # 使其能顯示中文
    plt.rcParams['axes.unicode_minus'] = False
    # print(label)
    for i in label:
        row_indx = np.where(yhat == i)
        plt.scatter(a[row_indx, 0], a[row_indx, 1], label="第" + str(i + 1) + "類", alpha=0.5, s=30)

    plt.ylabel('dim2')
    plt.xlabel('dim1')
    plt.legend()
    # plt.savefig()
    plt.show()


def adjust_rfc(X, y):
    from sklearn.model_selection import train_test_split
    from sklearn.ensemble import RandomForestClassifier
    from sklearn.model_selection import cross_val_score
    Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.3, random_state=0)
    tr = []
    te = []
    for i in range(10):
        clf = RandomForestClassifier(max_depth=i + 1)
        clf = clf.fit(Xtrain, ytrain)
        score_tr = clf.score(Xtrain, ytrain)
        score_te = cross_val_score(clf, X, y, cv=10).mean()
        tr.append(score_tr)
        te.append(score_te)
    print(max(te))
    plt.plot(range(1, 11), tr, label="train")
    plt.plot(range(1, 11), te, label="test")
    plt.xticks(range(1, 11))
    plt.legend()
    plt.show()

到了這里，關(guān)于基于邏輯回歸及隨機森林算法的冠心病預(yù)測與分析的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！