基于逻辑回归及随机森林算法的冠心病预测与分析

本文是一个课程报告，由我和另外一位同学合作完成。自我感觉做的还行决定放上来。

 数据集来源：Cardiovascular Study Dataset | Kaggle

目录

1.项目背景... 3

1.1项目说明... 3

1.2需求分析... 3

2.数据挖掘准备... 3

2.1数据字段含义介绍... 3

2.2基础统计分析... 4

3.数据挖掘过程... 5

3.1数据预处理... 5

3.1.1文字型变量数值化... 5

3.1.2缺失值处理... 6

3.1.3异常值处理... 8

3.1.4数据规范化... 10

3.2数据挖掘与可视化分析... 10

3.2.1人口统计信息分析... 11

3.2.2疾病史与亚健康状态分析... 13

3.2.3重要影响指标分析... 18

3.3特征工程... 19

3.3.1创建新特征... 19

3.3.2计算特征相关性... 20

3.3.3降维... 22

3.4建模    23

4.结果展示与评价... 25

4.1可视化结果展示... 25

4.2 模型结果展示... 26

5.总结... 28

5.1问题及解决方案... 28

5.2设计方案的优点及不足... 28

5.2.1优点... 28

5.2.2不足... 28

6.参考资料... 29

附录1：完整代码... 30

1.项目背景

1.1项目说明

冠心病，即冠状动脉性心脏病（cardio vascular disease），是一种常见的心脏病，也是夺走众多生命的重大疾病之一，世界卫生组织统计，一半的死亡是由心血管疾病造成的。截止至2021年，推算的中国冠心病现患人数为1139万[3]，且中国CVD患病率处于持续上升阶段。患有高血压、糖尿病等疾病，以及过度肥胖、不良生活习惯等是诱发该病的主要因素，以老年人为主要发病群体。

提前监测身体健康指标并预测心血管疾病的发作，可以最大程度上减少CVD的患病与伤亡数；同时，对已患病人群尽早干预，使其回归正确生活方式，可以帮助高危患者减少并发症。

综上所述，本研究旨在通过数据可视化与数据挖掘的方法，找出与冠心病病最相关的因素，探索数据之间潜在的联系，使用逻辑回归预测患病风险，来帮助相关人员预防并控制冠心病。

1.2需求分析

通过数据可视化分析与数据挖掘预测，深层探索数据之间的潜在联系，为冠心病的检测与预防发展提供依据。一方面分析个人的基本信息、生活习惯、病史等和冠心病风险的关系，判断哪些人群是患病高危人群，需要重点关注预防，并且如果和生活习惯有关，可以借此提出习惯改善建议；另一方面分析临床指标和冠心病风险的关系，用于辅助冠心病的临床诊断。再结合两方面的分析，对某个特定的人的冠心病风险进行预测，实现“早发现、早预防”。

2.数据挖掘准备

2.1数据字段含义介绍

数据集在Kaggle网站上公开发布，它来自马萨诸塞州，弗雷明汉镇的研究人员正在进行的心血管研究。数据集共有17列，其名称与含义如表格1：

表格 1 数据列语义信息

2.2基础统计分析

该数据集一共有3390条数据，16个特征，一个预测目标。其中，id字段用于作为每条数据的唯一标识，并无现实意义；sex和is_smoking两个字段为文字型变量，其余字段均为数值型变量。

对数据集进行缺失情况统计，统计结果如下表所示

表格 2 缺失值统计结果

glucose字段缺失值较多，缺失条目占到了总条目的近10%，需要进行缺失值处理。各列的数据统计情况分析如表格2与表格3。

表格 3 连续型数据统计情况

表格 4 离散型数据统计情况

3.数据挖掘过程

3.1数据预处理

3.1.1文字型变量数值化

将sex，is_smoking两个字段数值化：sex中F设为1，M设为-1，使得两者距离原点对称，对于原点具有相同的重要性；is_smoking中YES设为1，NO设为0。代码及结果如下

def numerical(data):
    """
    将文本变量进行编码
    """
    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缺失值处理

考虑到葡萄糖水平和是否患糖尿病可能有较强的相关性，计算葡萄糖水平和糖尿病的皮尔逊相关系数，并和其它指标进行比较，验证猜想

 corr = data['glucose'].corr(data['diabetes'])
    print('糖尿病和葡萄糖水平相关系数：', corr)
    corr = data['glucose'].corr(data['BMI'])
    print('肥胖程度和葡萄糖水平相关系数：', corr)
    corr = data['glucose'].corr(data['prevalentHyp'])
    print('高血压和葡萄糖水平相关系数：', corr)

可以看到，是否患有糖尿病与葡萄糖水平具有较高的关联，可以将人群分为患病和不患病两个类别，进一步分析，作出两者与葡萄糖关联的箱线图

图 1 糖尿病与葡萄糖水平

可以看到，患有糖尿病人群的葡萄糖水平分布更广，方差更大，均值更高，而未患病人群中葡萄糖水平的异常值更多，可能对均值产生较大影响，因此两者均采用众数进行填补。

另外对缺失值较少的字段离散值用众数填补，连续值用均值填补，代码和填补结果如下

# 填充葡萄糖
    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异常值处理

作出连续属性的箱线图

图 2 连续属性箱线图

可以看到，除cigsPerDay外，每个属性均有较多的值处于异常范围，若将其直接删掉，会丢失大量信息，影响数据分析的准确度。但同时每个指标的取值是否处于异常范围或许与是否患病有关，因此在此处计算获得每个连续指标的非异常范围，为后续特征分析做准备，代码和结果如下

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连续值正常范围

并且经过查询可知，该方法得出的正常值范围与一般医学上认定的正常值范围较为接近，认为该方法有效且该数据集正常。

3.1.4数据规范化

对临床指标进行Z值规范化处理，以消除单位不统一带来的量纲的影响，至此，数据预处理完成

3.2数据挖掘与可视化分析

联系数据列含义与问题背景，将从4个维度对数据进行分析。

3.2.1人口统计信息分析

绘制数据集的基本信息的面板图，如图1所示。其中，每一行代表着不同的教育水平，从上到下依次升序；每一列代表着不同的年龄，从最左列开始，以步长为5依次递增；面板中每个条形图表示着性别的分布情况，第一列为女性，第二列为男性；用深红色的表示该部分人吸烟，浅黄色则为不吸烟；最后，用面板的颜色深浅表示该部分人的平均患病风险情况。

图 3 人口统计数据面板图

可以看到，（1）数据的来源以中低教育水平、女性居多；（2）抽烟人群，以低教育水平，中年人居多，性别之间基本持平；（3）抽烟比例随年龄的上升，呈下降趋势，以女性尤为明显；（4）患病风险随年龄上升而显著增加。需要注意的是，在年龄分布的边缘区域偶见反常值，为样本数据过少而不可避免的偶然性，不应作为参考。

进一步，绘制抽烟比例与年龄之间的堆积百分比条形图，如图2所示，其中，左栏为女性，右栏为男性：

图 4 抽烟比例-年龄百分比堆积图

可以看到，随着年龄的增加，有抽烟习惯的人显著减少，女性部分的变化大于男性，与面板图情况一致。

3.2.2疾病史与亚健康状态分析

冠心病的发作与过往疾病史和身体的亚健康状态息息相关，通过分别绘制有无患病风险的人群的数据分布情况，观察其影响因素与大小。

绘制疾病记录与年龄桑基图，图5图6分别是无风险与有风险人群。其中，左侧的堆积条形图表示疾病史的分类情况，右侧的堆积条形图表示年龄的堆积情况，中间用曲线表示两部分的联系，曲线的粗细表明记录数的多少。

图 5 无风险人群患病与年龄桑基图

 图 6有风险人群患病与年龄桑基图

可以看到，（1）两类人中有过高血压经历的人均明显占多，数据呈有偏分布；（2）两类人中，年龄越大，患病的经历就越多。不同的是，（3）有风险人群有中风和糖尿病的比例更多些；（4）有风险人群中老年人群所占的比例更大些。这一点与常识相符。

桑基图在边缘地区产生畸变，同样是因为样本偶然性原因，不应考虑。

值得注意的是，在两类图中，年龄以中间多两边少的趋势呈近似正态分布，绘制年龄与患病记录直方图如图7。

 图 7年龄-患病记录直方图

    可以得到结果：（1）列1中，整体呈左偏分布，说明年龄越大，没有患病记录的人就越少了，这两个因素之间为正相关；（2）列2及列3，整体呈右偏分布，进一步证实相关性；（3）且随着患病记录的增加，有风险人群的占比越来越大，说明患病会增加冠心病发作的概率。

各种指标是身体状态的直接表示，考虑有风险与无风险之间的指标差异对比，绘制箱型图如图8。

 图 8 有无风险人群各身体指标箱型图

浅绿色的列为无风险，橙色的列为有风险。可以看到，差距仅仅是有风险列的方差过大。说明冠心病是由多种因素并发导致的，风险是否与单个指标的关系并不明显。

基于此，收集六大指标在正常范围之外的信息，并将其表注为亚健康记录。

表格 6 亚健康记录划分标准

身体亚健康与不良习惯有关，即每日抽烟的根数。探索每日香烟根数的分布情况，患病记录与亚健康记录之间的关系，绘制香烟根数-亚健康记录-患病记录直方图，如图9。

 图 9 香烟根数-亚健康记录-患病记录直方图

绿色条状图为吸烟根数的分布，小部分吸烟人群每日香烟数集中在18-24根；红色折线图表示人群平均的患病记录，黄色折线表示人群平均的亚健康记录。两条折线之间仅有很小的相似性，折线与香烟根数之间的相关性也有限，判断该分析以香烟为维度结果不显著。以香烟为不良习惯的代表有很大局限性。

但患病记录与亚健康记录之间的关系确实存在，绘制瀑布图，如图10。

图 10 亚健康记录-患病记录瀑布图

每有一个指标为亚健康，患病记录的增加成指数型增长。

综上所述，单个身体健康指标对患冠心病风险的作用不大，但多个健康指标会综合影响患病记录，患病记录再增加冠心病风险的可能。

3.2.3重要影响指标分析

冠心病为经典的老年病症，且与肥胖程度高度相关，综合考虑BMI与年龄对冠心病风险的影响，画体重-年龄矩阵密度图，如图11，其中，浅绿色的点为无风险，橙色的点为有风险，颜色深浅表示落在这个坐标的样本数的多少。

图 11 体重-年龄矩阵密度图

    与无风险人群相比，有风险人群大多落在第一、第二、第四象限，有着高龄高体重的特征，且有风险人群的BMI平均值大于无风险人群。

3.3特征工程

3.3.1创建新特征

根据以上分析，在Python中创建“亚健康计数”特征

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计算特征相关性

计算每个特征的皮尔逊相关系数，画出归一化和未归一化的热力图，进行特征比较

 图 12 热力图-未归一化

图 13 热力图-归一化

不管是否归一化，是否有十年冠心病风险与单个特征的相关性都很小，且和年龄相关性最大，和tableau分析相符。且归一化后削弱了各个特征之间的相关性，减弱了线性回归中的多重共线性问题。但可以看到，是否患有糖尿病和收缩压、亚健康记录数和是否服用过降压药等特征之间仍有很强的相关性，因此进行人工特征选择降维或者进行PCA降维。

3.3.3降维

根据以上分析，初步人工选择特征为”age””prevalentHyp”“subHealth”三个特征。进行PCA降维。先降为两维，画出聚类散点图

def pca(data):
    """
    pca降维，作图，找出各成分和主成分关系
    """
    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("降维后的数据为：", a)
    print("各主成分贡献率为：", pca.explained_variance_ratio_)

    from sklearn.cluster import KMeans
    model = KMeans(n_clusters=2)
    model.fit(a)
    yhat = model.predict(a)  # 得出每个点所属类别的预测
    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 聚类散点图

当只保留两个主成分时，主成分累计贡献率已经达到了96%以上，可以只保留两个主成分进行后续建模。两个主成分呈现出条带形的关系，展现出一种既有序又难以用函数描述的关联。对散点图出现原因进行推测，可能是原始维度中既有离散型变量又有连续型变量，连续变量使得降维后的散点图呈条带状，离散型变量使得各个条带分立开来。同时可以注意到，聚类结果中一个条带上的点均属于同一类别，说明离散特征是影响是否有患病风险的重要特征，这也与之前的分析相符。

3.4建模

将上文中不同的数据处理方法和不同的机器学习方法相结合，建立12种组合模型，如下表所示

表格 7 组合模型参数

对随机森林模型进行调参

图 15 随机森林最大深度调整

可以看到最优深度为2，模型再深时，训练集准确率一直上升，测试集准确率下降，模型过拟合。接着建立模型。部分代码如下，完整代码见附录

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.结果展示与评价

4.1可视化结果展示

   通过上述可视化图表，制作仪表盘与故事，呈现数据分析过程与结论。

图 16 故事预览图

4.2 模型结果展示

可以得到以上十二个模型的ROC曲线和测试集上的预测准确率如表格8和图17。可以看到，最简单的模型1，即不进行任何数据处理，直接采用逻辑回归进行预测，在测试集上的准确率是最高的。而在特征选择上，采用人工特征选择和PCA降维并无较大差异，两者相互印证，说明年龄、亚健康记录数、是否服用降压药确实与冠心病风险密切相关。另外从ROC曲线来看，随机森林的表现效果比逻辑回归要更好，该模型可能在其它类似数据集上预测准确率更稳定。

表格 8 测试准确率

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.总结

5.1问题及解决方案

1. 因人口信息维度的数量多，不知道从哪个视角开始分析，才不会遗漏重要的相关性。后综合考虑，决定统一做一张面板图：通过新建列与行的离散字段，获得条形图面板与面积图面板，在仪表盘叠加，最终得到了有详细信息的人口信息统计图：人口信息的任意两个维度在此图上都能够有机会得到呈现相关性的机会。
2. 在分析患病风险与疾病史与亚健康状态时，最开始因六大指标与三大疾病史在有无风险人群中的分布区别不够明显，找不到能呈现显著结果分析维度，无法作数据分析，甚至怀疑冠心病是否没有主要的致病因素。后利用疾病史与亚健康记录的累加值新建特征，得到了较好的结果，发现单一指标或维度不能导致患病风险的提升，患病是所有因素综合作用的结果。
3. 对缺失值进行填充时，葡萄糖水平缺失过多，如果只是用均值填充感觉会引入大量噪声。猜测葡萄糖水平应该和是否患糖尿病有关联，于是用是否患糖尿病将人群分类再填充，希望能填充的更准确一些。
4. 在tableau可视化分析的基础上做Python建模时，开始有些迷惑不知道该怎样把两个部分整合起来形成一个整体。后来用可视化分析出的“亚健康记录”这个指标作为突破口，在此基础上进行特征工程，将两个部分结合起来

5.2设计方案的优点及不足

5.2.1优点

1. 首先通过探索性数据分析，观察分析各类图的可视化结果，得到了可能相关的初步因素与相关性结论。然后在通过数据挖掘的方法，用机器学习算法进行降维、预测等操作，得到最终结论。
2. 采用多种分析工具相结合，能够综合各个工具的优势进行分析，且得出的结论可以相互印证
3. 发现了“亚健康记录数”这一重要指标，将各个连续的特征联系了起来。
4. 使用了多种模型进行模型消融实验，并采用测试集准确率和ROC两种方法进行模型评价，能够选出效果最好的模型

5.2.2不足

1. 特征之间的关联挖掘不够深入。还可以查阅相关文献资料，找到各个特征间可能的相互关系，再进行可视化分析进行验证分析。
2. 模型具有局限性。可以看到在测试集的预测准确率上，各个模型相差并不大，可以采用更多的模型进行尝试，找到预测效果更好的模型。
3. BP Meds数据始终没能用上，尝试过多种数据维度分析，效果都不够显著。再有机会进行深度挖掘时，该变量可能会派上用场

6.参考资料 

1. Cardiovascular Risk Prediction | Kaggle
2. 冠心病 - 医学百科 (yixue.com)
3. 《中国心血管健康与疾病报告2021》要点解读[J].中国心血管杂志,2022,27(04):305-318.

附录1：完整代码文章来源地址https://www.toymoban.com/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):
    """
    统计各字段的缺失值
    """
    print('缺失值统计：', data.isnull().sum())
    print(data.info())
    print('数值分布统计', data.describe())


def numerical(data):
    """
    将文本变量进行编码
    """
    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):
    """
    填充缺失值
    """
    # 计算其它指标和葡萄糖水平的相关系数
    corr = data['glucose'].corr(data['diabetes'])
    print('糖尿病和葡萄糖水平相关系数：', corr)
    corr = data['glucose'].corr(data['BMI'])
    print('肥胖程度和葡萄糖水平相关系数：', corr)
    corr = data['glucose'].corr(data['prevalentHyp'])
    print('高血压和葡萄糖水平相关系数：', 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  # 实现z-score标准化
    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降维，作图，找出各成分和主成分关系
    """
    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("降维后的数据为：", a)
    print("各主成分贡献率为：", pca.explained_variance_ratio_)

    from sklearn.cluster import KMeans
    model = KMeans(n_clusters=2)
    model.fit(a)
    yhat = model.predict(a)  # 得出每个点所属类别的预测
    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()

到了这里，关于基于逻辑回归及随机森林算法的冠心病预测与分析的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！