【Python机器学习项目】项目一：心脏病二分类问题

使用机器学习预测心脏病

根据一些病理学属性预测心脏病

特别说明：

1. 开新坑啦！本系列共2个项目，难度不大，特别适合新手入坑

2. 由于本项目只是系列课程的第一个项目，所以很多细节不深挖，仅做示范，在第二个项目中再完善。

以下为整体思路概述

1. 问题定义

给定一个病人的临床诊断，能否预测他们是否患有心脏病？

2. 数据来源

https://archive.ics.uci.edu/ml/datasets/Heart+Disease

3. 评估

期望准确率达到95%

4. 特征和标签

数据字典

1. age: age in years
2. sex: sex (1 = male; 0 = female)
3. cp: chest pain type

• – Value 0: typical angina
• – Value 1: atypical angina
• – Value 2: non-anginal pain
• – Value 3: asymptomatic

4. trestbps: resting blood pressure (in mm Hg on admission to the hospital)
5. chol: serum cholestoral in mg/dl
6. fbs: (fasting blood sugar > 120 mg/dl) (1 = true; 0 = false)
7. restecg: resting electrocardiographic results

• – Value 0: normal
• – Value 1: having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV)
• – Value 2: showing probable or definite left ventricular hypertrophy by Estes’ criteria

8. thalach: maximum heart rate achieved
9. exang: exercise induced angina (1 = yes; 0 = no)
10. oldpeak = ST depression induced by exercise relative to rest
11. slope: the slope of the peak exercise ST segment

• – Value 0: upsloping
• – Value 1: flat
• – Value 2: downsloping

12. ca: number of major vessels (0-3) colored by flourosopy
13. thal: 0 = normal; 1 = fixed defect; 2 = reversable defect
14. target: 0 = no disease, 1 = disease
    
    

0. 导包

# EDA
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
sns.set()
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
%config InlineBackend.figure_config = 'svg'

# sklearn模型
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

# 模型评估
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.model_selection import RandomizedSearchCV, GridSearchCV
from sklearn.metrics import confusion_matrix, classification_report
from sklearn.metrics import precision_score, recall_score, f1_score
from sklearn.metrics import plot_roc_curve
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

载入数据

hd_df = pd.read_csv('heart-disease.csv')
hd_df.shape
1
2

(303, 14)
1

1. EDA

了解更多有关这个数据集的信息，成为该数据集的懂王

1. 要解决什么问题？
2. 都有些什么数据，要怎么处理？
3. 有无缺失值，如何处理？
4. 有无异常值，如何处理？
5. 如何通过创建衍生特征、处理和筛选现有特征得到更多信息？

hd_df.head()
1

hd_df.tail()
1

# 查看样本分布
targets = hd_df['target'].value_counts()
targets
1
2
3

1    165
0    138
Name: target, dtype: int64
1
2
3

targets.plot(
    kind='bar', 
    color=['salmon', 'lightblue'],
    figsize=(10,6)
)
plt.xticks(rotation=0)
plt.show()
1
2
3
4
5
6
7

​

hd_df.info()
1

RangeIndex: 303 entries, 0 to 302
Data columns (total 14 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   age       303 non-null    int64  
 1   sex       303 non-null    int64  
 2   cp        303 non-null    int64  
 3   trestbps  303 non-null    int64  
 4   chol      303 non-null    int64  
 5   fbs       303 non-null    int64  
 6   restecg   303 non-null    int64  
 7   thalach   303 non-null    int64  
 8   exang     303 non-null    int64  
 9   oldpeak   303 non-null    float64
 10  slope     303 non-null    int64  
 11  ca        303 non-null    int64  
 12  thal      303 non-null    int64  
 13  target    303 non-null    int64  
dtypes: float64(1), int64(13)
memory usage: 33.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

# 查看缺失值
hd_df.isna().sum()
1
2

age         0
sex         0
cp          0
trestbps    0
chol        0
fbs         0
restecg     0
thalach     0
exang       0
oldpeak     0
slope       0
ca          0
thal        0
target      0
dtype: int64
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 查看描述性统计信息
hd_df.describe([0.01, 0.25, 0.5, 0.75, 0.99]).T
1
2

查看性别和标签之间的关系

hd_df['sex'].value_counts()
1

1    207
0     96
Name: sex, dtype: int64
1
2
3

# cross_tab改进版函数
def to_cross_tab(origin_df, index_name, col_name):
    df = pd.crosstab(origin_df[index_name], origin_df[col_name])
    df['rate'] = df.iloc[:,1] / (df.iloc[:,0] + df.iloc[:,1])
    return df
1
2
3
4
5

sex_target_df = to_cross_tab(hd_df, 'target', 'sex')
sex_target_df
1
2

# 方便绘图的函数
def to_plot(df, title, xlabel, ylabel, legend):
    df.plot(
    kind='bar', 
    color=['lightblue', 'salmon'],
    figsize=(10,6)
)
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.xticks(rotation=0)
    plt.legend(legend)
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13

to_plot(sex_target_df[[0,1]], '按性别统计的心脏病概率', '0 = 女生， 1 = 男生', '总人数', ['未得病', '得病'])
1

​

​

明显女性发病率高得多

查看得病/未得病两类人中年龄和最大心率的关系

plt.figure(figsize=(10,6))

# 查看得病人群
plt.scatter(hd_df['age'][hd_df['target']==1],
            hd_df['thalach'][hd_df['target']==1],
            c='salmon'
)

# 查看未得病人群
plt.scatter(hd_df['age'][hd_df['target']==0],
            hd_df['thalach'][hd_df['target']==0],
            c='lightblue'
)

# 说明
plt.title('根据是否得心脏病分成两类人群来查看年龄和最大心率')
plt.xlabel('年龄')
plt.ylabel('最大心率')
plt.legend(['得病', '未得病'])

plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

​

​

# 查看年龄分布
hd_df['age'].hist()
1
2

1

​

# 做正态性检验
stats.normaltest(hd_df['age'])
1
2

NormaltestResult(statistic=8.74798581312778, pvalue=0.012600826063683705)
1

年龄符合正态分布

查看心绞痛类型和标签之间的关系

3. cp: chest pain type

• – Value 0: typical angina
• – Value 1: atypical angina
• – Value 2: non-anginal pain
• – Value 3: asymptomatic

cp_target_df = to_cross_tab(hd_df, 'cp', 'target')
cp_target_df
1
2

to_plot(cp_target_df[[0,1]], '按心绞痛类型统计的心脏病人数', '心绞痛类型', '总人数', ['未得病', '得病'])
1

​

​

# 相关系数矩阵
corr_matrix = hd_df.corr()
corr_matrix
1
2
3

plt.figure(figsize=(14, 10))
sns.heatmap(
    corr_matrix, 
    vmin=-1, 
    annot=True, 
    linewidth=5, 
    fmt='.2f', 
    cmap='YlGnBu'
)
plt.show()
1
2
3
4
5
6
7
8
9
10

​

​

这个相关性看起来还是比较好的，大部分特征和标签之间都有一定的相关性，且特征之间也没有相关性>0.8的需要排除。当然，真的想看相关性还得分类别变量和连续值变量，连续值变量又得做正态检验。

3. 建模

hd_df.head()
1

X = hd_df.drop(columns=['target'])
y = hd_df['target']
1
2

# 设置随机种子，便于其他人重复实验
np.random.seed(13)

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
1
2
3
4
5

依次使用逻辑斯蒂回归、KNN、随机森林

# 创建字典
models = {
    'lr': LogisticRegression(),
    'knn': KNeighborsClassifier(),
    'rf': RandomForestClassifier()
}

# 一个简单的试探性fit和score的函数
def fit_and_score(models, X_train, X_test, y_train, y_test):
    np.random.seed(13)
    model_score = {}
    for name, model in models.items():
        model.fit(X_train, y_train)
        model_score[name] = model.score(X_test, y_test)
    return model_score
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

model_scores = fit_and_score(models, X_train, X_test, y_train, y_test)
model_scores
1
2

 {'lr': 0.8360655737704918, 'knn': 0.639344262295082, 'rf': 0.819672131147541}
1

模型比较

model_compare = pd.DataFrame(model_scores, index=['正确率'])
model_compare.T.plot(kind='bar')
plt.xticks(rotation=0)
plt.show()
1
2
3
4

​

接下来做什么？

• 超参数优化
• 特征重要性
• 混淆矩阵
• 交叉验证
• 精确率
• 召回率
• F1 score
• 分类报告
• ROC
• AUC

# knn调参（假装不会GSCV和RSCV）
train_scores = []
test_scores = []

neighbors = range(1, 21)

knn = KNeighborsClassifier()
for i in n_neighbors:
    knn.set_params(n_neighbors=i)
    knn.fit(X_train, y_train)
    train_scores.append(knn.score(X_train, y_train))
    test_scores.append(knn.score(X_test, y_test))    
1
2
3
4
5
6
7
8
9
10
11
12

train_scores
1

[1.0,
 0.8016528925619835,
 0.8057851239669421,
 0.7603305785123967,
 0.768595041322314,
 0.7355371900826446,
 0.7396694214876033,
 0.71900826446281,
 0.7024793388429752,
 0.6900826446280992,
 0.7107438016528925,
 0.6859504132231405,
 0.7024793388429752,
 0.6776859504132231,
 0.6942148760330579,
 0.6859504132231405,
 0.6694214876033058,
 0.6859504132231405,
 0.7024793388429752,
 0.7066115702479339]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

test_scores
1

[0.6065573770491803,
 0.4426229508196721,
 0.5737704918032787,
 0.5409836065573771,
 0.639344262295082,
 0.6557377049180327,
 0.6065573770491803,
 0.6721311475409836,
 0.6557377049180327,
 0.6557377049180327,
 0.6885245901639344,
 0.6885245901639344,
 0.6885245901639344,
 0.7377049180327869,
 0.7213114754098361,
 0.7213114754098361,
 0.7213114754098361,
 0.7049180327868853,
 0.7377049180327869,
 0.7377049180327869]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

plt.plot(neighbors, train_scores, label='Train score')
plt.plot(neighbors, test_scores, label='Test score')
plt.xticks(range(1,21,1))
plt.xlabel('n_neighbors参数值')
plt.ylabel('正确率')
plt.legend()
plt.show()
1
2
3
4
5
6
7

​

​

knn最高分也没达到80%正确率，放弃

使用RandomizedSearchCV调参

# 逻辑斯蒂回归
# 由于主要是想找最优C值，其他参数就不设置了，并且这里使用np.logspace故意把C值分布得开一些，因为完全不知道在哪里取得最优值
log_reg_grid = {
    'C':np.logspace(-4, 4, 20),
    'solver': ['liblinear']
}

# 随机森林
rf_grid = {
    'n_estimators': np.arange(10, 1000, 50),
    'max_depth': [None, 3, 5, 10],
    'min_samples_split': np.arange(2, 20, 2),
    'min_samples_leaf': np.arange(1, 20, 2)
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14

np.random.seed(13)

# 实例化RSCV对象
rs_log_reg = RandomizedSearchCV(
    LogisticRegression(),
    param_distributions=log_reg_grid,
    cv=5,
    n_iter=20,
    verbose=True
)
# fit
rs_log_reg.fit(X_train, y_train)
1
2
3
4
5
6
7
8
9
10
11
12

Fitting 5 folds for each of 20 candidates, totalling 100 fits
1

RandomizedSearchCV(cv=5, estimator=LogisticRegression(), n_iter=20,
               param_distributions={'C': array([1.00000000e-04, 2.63665090e-04, 6.95192796e-04, 1.83298071e-03,
   4.83293024e-03, 1.27427499e-02, 3.35981829e-02, 8.85866790e-02,
   2.33572147e-01, 6.15848211e-01, 1.62377674e+00, 4.28133240e+00,
   1.12883789e+01, 2.97635144e+01, 7.84759970e+01, 2.06913808e+02,
   5.45559478e+02, 1.43844989e+03, 3.79269019e+03, 1.00000000e+04]),
                                    'solver': ['liblinear']},
               verbose=True)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. 
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

RandomizedSearchCV(cv=5, estimator=LogisticRegression(), n_iter=20,
               param_distributions={'C': array([1.00000000e-04, 2.63665090e-04, 6.95192796e-04, 1.83298071e-03,
   4.83293024e-03, 1.27427499e-02, 3.35981829e-02, 8.85866790e-02,
   2.33572147e-01, 6.15848211e-01, 1.62377674e+00, 4.28133240e+00,
   1.12883789e+01, 2.97635144e+01, 7.84759970e+01, 2.06913808e+02,
   5.45559478e+02, 1.43844989e+03, 3.79269019e+03, 1.00000000e+04]),
                                    'solver': ['liblinear']},
               verbose=True)

LogisticRegression()

LogisticRegression()

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14