精品教学案例 | 基于分类算法的肝病诊断

2020-05-19 17:41:56 浏览数 (1)

查看本案例完整的数据、代码和报告请登录数据酷客(http://cookdata.cn)案例板块。

本案例适合作为大数据专业数据科学导引或机器学习实践课程的配套教学案例。通过本案例,能够达到以下教学效果:

  • 培养学生对真实数据的探索分析能力。案例中使用Pandas对数据进行清洗,对特征进行了年龄性别做交叉分析,并对特征做相关分析,了解各特征之间的关系。
  • 帮助学生掌握数据绘图工具。例如在案例中使用Matplotlib做柱状图,清晰的了解各特征的分布,并且用Seaborn将相关矩阵做可视化呈现。
  • 提高学生机器学习建模的能力。案例中使用了两种思想建模,一种是单一的分类器,另一种是集成学习,通过使用两种方法建模增加学生对机器学习模型的理解,并提高建模的实践能力。

肝脏是人体内最大的消化腺,帮助分解体内毒素,是人体物质能量代谢的中心站和人体最大的解毒器官,更是维持生命活动必不可少的器官。当下,由于人们生活水平的提升,饮食、环境等方面的改变,长期过量饮酒,有害气体和受污染的食物的摄入,熬夜过度等不良生活习惯增加,致使肝病患者不断增加,肝病的诊断与治疗成为医疗行业的一大重要问题。

常见的肝功能诊断中,主要包括三大类的指标:血清酶、胆红素和血清蛋白。其中,血清酶中的医学指标主要包括丙氨酸氨基转移酶、天冬氨酸氨基转移酶和碱性磷酸酶等,当肝脏细胞被破坏时,酶会被大量释放到血液中,引起指标上升。胆红素指标包括总胆红素、直接胆红素和间接胆红素等,它们反映了胆红素的代谢情况,当肝细胞变性坏死,胆红素代谢出现障碍时,胆红素指标会升高。血清蛋白指标反映了肝脏的合成功能,其包含白蛋白、球蛋白、总蛋白等,可用于检测慢性肝损伤、机体免疫等情况。将患者的医疗检测数据与预测算法相结合,有助于帮助医生更精确地做出诊断,也可减轻医生的负担。

1.数据说明与预处理

我们知道患者的患病情况不仅与其医疗检测指标有关,还与其自身的生理特征有关。数据显示,肝病的分布呈现“重男轻女”的现象,男女患病比例约为5:1。另外,肝病患者的年龄大多集中在45岁以上,即在中老年人群中发病率较高。

本案例使用了UCI公开的肝病患者诊断数据,基于患者的生理指标(性别、年龄)和医疗检测指标来对患病情况进行分析,数据共计583条,其中患肝病人数为416,未患肝病人数为167。对于医疗检测指标,命名较长,后面将考虑用其缩写来代替,而生理指标的缩写与原列名保持相同。各数据指标含义如下表所示:

首先,我们读入数据,并将数据的列名替换成各个特征的缩写形式。

代码语言:javascript复制
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
patient = pd.read_csv("./input/indian_liver_patient.csv")
name = ["Age","Gender","TBIL","DBIL","ALP","ALT","AST","TP","ALB","A/G","Dataset"]
patient.columns = name
patient.head()
代码语言:javascript复制
patient.info()

从特征的情况来看,只有性别(Gender)为定性的特征,下面将考虑将其转化为数值型特征,即男性设为0,女性设为1。另外,数据中白蛋白和球蛋白比率(A/G)存在数据缺失,需要进行处理。

代码语言:javascript复制
#Gender特征转换
patient["Gender"].replace("Female",1,inplace = True)
patient["Gender"].replace("Male",0,inplace = True)
patient["Gender"] = patient["Gender"].astype("int8")
代码语言:javascript复制
#缺失值处理
patient.isnull().sum(axis = 0)
代码语言:javascript复制
print("The proportion of missing values is %.5f  "%(patient.isnull().sum(axis = 0).max()/patient.shape[0]))  #缺失值占比

观察到缺失数据的占比仅为0.686%,占比较小,故考虑直接做删除处理,这里使用dropna()方法。

代码语言:javascript复制
patient.dropna(inplace = True)
print("The shape of data is ",patient.shape)

2.探索性分析

由于数据特征较多,我们需要探究一下各特征之间的关系,这将有助于后面建立合理有效的模型。下面将从几个方面进行入手:

  • 首先,采用直方图观察各个医疗检测指标的分布情况;
  • 其次,使用分组柱状图分析患病与性别之间的关系;
  • 然后,使用堆积直方图分析患病与年龄之间的关系;
  • 最后,使用特征的相关系数图分析各个特征之间的相关性。

2.1 医疗指标的分布情况

代码语言:javascript复制
fig,ax = plt.subplots(2,4,figsize = (15,7))
for i in range(1,9):
    sub = int("24" str(i))
    plt.subplot(sub)
    plt.hist(patient[name[i 1]],bins = 15)
    plt.xlabel(name[i 1],fontsize = 13)
    plt.ylabel("Count",fontsize = 13)  
plt.tight_layout()
plt.show()   

观察到总胆红素(TBIL)直接胆红素(DBIL)碱性磷酸酶(ALP)谷丙转氨酶(ALT)天冬氨酸氨基转移酶(AST)这5个特征呈现明显的右偏分布,这是由于个别严重肝病患者的这些指标会明显高于大众水平。其中,总胆红素(TBIL)直接胆红素(DBIL)这两个特征的分布虽有一定的右偏,但数据的极差较小,暂不做处理;碱性磷酸酶(ALP)谷丙转氨酶(ALT)天冬氨酸氨基转移酶(AST)这三个特征的最大值与最小值的差距在千倍以上,这里考虑将它们进行对数变换,压缩数据的尺度,使数据更为平稳,这将有助于后续的建模,使用的方法为NumPy中的log()

代码语言:javascript复制
#对数变换
patient[["ALP","ALT","AST"]] = np.log(patient[["ALP","ALT","AST"]])

下面,我们查看一下对数变换后的碱性磷酸酶(ALP)谷丙转氨酶(ALT)天冬氨酸氨基转移酶(AST)这三个特征的数据分布情况。

代码语言:javascript复制
fig,ax = plt.subplots(1,3,figsize = (12,4))
for i in range(1,4):
    sub = int("13" str(i))
    plt.subplot(sub)
    plt.hist(patient[name[i 3]],bins = 15)
    plt.xlabel(name[i 3],fontsize = 13)
    plt.ylabel("Count",fontsize = 13)  
plt.tight_layout()
plt.show()   

可以发现,对数变换后的数据分布范围缩小,且更接近于正态,这将有助于后续的分析。

2.2 患病与性别的分布情况

代码语言:javascript复制
fig, ax = plt.subplots(figsize = (7,5))
name_list = ['Male','Fmale'] 
gen_count = patient.groupby(['Gender','Dataset']).Age.count()
ill_count = [gen_count[0][1],gen_count[1][1]]  
notill_count = [gen_count[0][2],gen_count[1][2]]  
x = list(range(len(name_list)))  
total_width, n = 0.8, 2  
width = total_width / n  

plt.bar(x,ill_count,width = width,label = 'ill',color = 'lightsalmon')  
for i in range(len(x)):  
    x[i] = x[i]   width  
plt.bar(x,notill_count, width = width, label = 'notill',tick_label = name_list,color = 'lightskyblue')  
plt.legend(loc = "upper right")  
plt.xlabel("Gender",fontsize = 13)
plt.ylabel("Count",fontsize = 13)
plt.title('Patient by gender')
plt.show()  

从图中可以看出,患病人群中的男性人数约为女性的三倍,这与现实中的肝病人群分布略有不同。其中,男性中患肝病人数与未患肝病人数的比例约为3:1,女性中患肝病人数与未患肝病人数的比例约为2:1。性别对是否患肝病可能会有一定的影响。

2.3 患病与年龄的分布情况

代码语言:javascript复制
fig, ax = plt.subplots(figsize = (7,5))
x = [patient[patient["Dataset"] == 1]["Age"],patient[patient["Dataset"] == 2]["Age"]]
bins = range(0, 91, 10)

plt.hist(x, bins = bins, color = ['lightsalmon','lightskyblue'], histtype = "bar", rwidth = 10,stacked = True, label = ["ill", "notill"])
plt.xlabel("Age",fontsize = 13)
plt.ylabel("Count",fontsize = 13)
plt.title("Patient by Age")
plt.legend(loc = "upper right")
plt.show() 

从堆积直方图中可以看出,收集的数据涵盖了各个年龄段的人群,其中以30岁-80岁为主。观察到,40-60岁和70-80岁的人群中,患肝病比例较高,这与现实情况是相符的。现代人由于长期不良的生活习惯的影响,致使疾病朝着“年轻化”的方向发展,中年人群已成为疾病发展的高危人群,而老年人群由于自身免疫力的下降也极易受疾病的困扰。因此,是否患肝病与年龄之间可能存在着密切的关系。

2.4 特征间的相关性

在Python中,可以直接对数据框使用.corr()来计算相关性矩阵。

代码语言:javascript复制
#特征相关系数图
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False
fig, ax = plt.subplots(figsize = (8,8))
corr=patient.corr()
sns.heatmap(corr,annot = True, vmax = 1,vmin=-1,square = True, cmap = "Blues")
plt.ylabel('Feature',fontsize = 15)
plt.xlabel('Feature',fontsize = 15)
plt.title("The correlation of features",fontsize = 15)
plt.show()

从特征之间的相关系数图可以看出,颜色越深,特征之间的正相关性越强;颜色越浅,特征之间的负相关性越强。其中,总胆红素(TBIL)直接胆红素(DBIL)的相关系数为0.87,它们均是衡量血液中的胆红素情况的指标;谷丙转氨酶(ALT)天冬氨酸氨基转移酶(AST)的相关系数为0.84,它们衡量的是血液中血清酶的含量;白蛋白(ALB)总蛋白(TP)白蛋白和球蛋白比率(A/G)的相关系数均在0.6以上,它们衡量的是血清蛋白的情况。整体来看,特征之间有着较强的相关性,在后续进行建模时,需要重点考虑模型的特征选择问题。

3.分类建模

判断患者是否患病本质上是一个二分类问题,在后文将考虑分别使用逻辑回归、决策树这两种单一的分类器和随机森林、AdaBoost这两种集成方法进行建模。首先,将数据集按照3:1划分成训练集和测试集,使用的方法为train_test_split;之后定义一个模型评价函数evaluate,用于输出模型的预测报告、混淆矩阵和分类正确率。

代码语言:javascript复制
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score   

##训练集测试集的划分
x = patient.drop("Dataset",axis = 1)
y = patient["Dataset"]
y.replace(2,0,inplace = True)
x_train, x_test, y_train, y_test = train_test_split(x, y,test_size = 0.25, random_state = 2)
代码语言:javascript复制
#定义一个模型评价函数
def evaluate(method,y_test,y_pred):
    print(classification_report(y_test,y_pred))
    print(confusion_matrix(y_test, y_pred))
    print("The acc of %s is %.5f"%(method,accuracy_score(y_test,y_pred)))

3.1 逐步逻辑回归

由特征间的相关系数图可知,部分特征之间有较强的相关性,故而不能直接将所有的特征都放入逻辑回归模型中进行训练,需要进行筛选。在回归模型中,可以考虑使用逐步回归的方法,筛选出对因变量影响最显著的那些特征。

在逐步逻辑回归中,每次只引入或剔除一个特征,具体步骤为: (1) 初始时,模型中无特征,我们选取偏回归平方和最大的那个特征进行F检验,如果检验的p_value小于0.1,则说明该特征对因变量有显著影响,将被选入模型; (2) 模型引入新特征后,我们选取偏回归平方和最小的那个特征进行F检验,如果检验的p_value大于0.15,则说明该特征对因变量的影响不显著,故将该特征从模型中剔除; (3) 重复上述过程,直到模型中没有特征被引入或剔除为止。

由于Python中没有直接的API可以进行逐步逻辑回归,所以需要借助StatsModels中Logit类进行函数编写,再根据回归结果中的各特征显著性进行选择和剔除特征。

代码语言:javascript复制
#自定义逐步回归的函数
import statsmodels.api as sm
def stepwise(x,y,alpha_in = 0.1,alpha_out = 0.15):
    '''x为所有可能的特征构成的DataFrame,
      y为响应变量,
      alpha_in为引入特征的显著性水平,
      alpha_out为剔除特征的显著性水平'''
    included = []
    while True:
        changed = False
        excluded = list(set(x.columns)-set(included))
        p_val = pd.Series(index = excluded)
        for new_column in excluded:
            model = sm.Logit(y,sm.add_constant(x[included [new_column]])).fit()
            p_val[new_column] = model.pvalues[new_column]
        min_pval = p_val.min()
        #forward step
        if min_pval < alpha_in:
            changed = True
            add_feature = p_val.idxmin()
            included.append(add_feature)
            print("Add {} with p_value   {:.6}".format(add_feature,min_pval))
        #backward step
        model = sm.Logit(y,sm.add_constant(x[included])).fit()
        p_val = model.pvalues.iloc[1:]
        max_pval = p_val.max()
        if max_pval > alpha_out:
            changed = True
            drop_feature = p_val.idxmax()
            included.remove(drop_feature)
            print("Drop {} with p_value   {:.6}".format(drop_feature,max_pval))
        if not changed:
            break
    return included
代码语言:javascript复制
result = stepwise(x_train,y_train) #这里的结果较长,暂不做展示

在逐步逻辑回归的过程中,依次被纳入模型的特征为天冬氨酸氨基转移酶(AST)年龄(Age)谷丙转氨酶(ALT)直接胆红素(DBIL),由于直接胆红素(DBIL)的引入,导致天冬氨酸氨基转移酶(AST)对因变量的影响不再显著,以p_value等于0.345被剔除模型。之后被纳入模型的特征为碱性磷酸酶(ALP),故最终的逻辑回归模型中将包含四个特征。下面我们将查看逻辑回归模型的结果。

代码语言:javascript复制
new_columns = result
model = sm.Logit(y_train,sm.add_constant(x_train[new_columns])).fit()
model.summary2()

从逐步逻辑回归的结果来看,年龄(Age)谷丙转氨酶(ALT)直接胆红素(DBIL)碱性磷酸酶(ALP)这四个特征被引入模型,且所有特征的p_value均小于0.1,特征显著,模型的整体p_value小于0.01,模型有效。另外,除截距项之外,各个特征的系数均为正值,说明在控制其它因素不变的情况下,上述四个特征之一增加,均会导致患者患肝病的概率增加。下面在测试集上进行预测。

代码语言:javascript复制
y_pred = model.predict(sm.add_constant(x_test[new_columns]))
y_pred1 = [1 if y>0.6 else 0 for y in y_pred]
evaluate("StepLogit",y_test,y_pred1)

逐步逻辑回归的 F1值为0.69,预测正确率为67.6%左右。

3.2 决策树

决策树的建模过程主要包括:特征选择、树的生成和剪枝。在特征选择的过程中,我们可以选择使用信息增益、信息增益比和基尼系数三种方法来进行特征选择,每次选择出对训练数据最有分类能力的特征。在Python建决策树的过程中,最常使用的方法为基尼系数,其定义为:

表示数据中共有K个类,pk 为样本属于第k个类的概率。

在Python中,使用sklearn.tree的DecisionTreeClassifier类进行建模,建模过程较为简单,但重在参数调试方面,决策树DecisionTreeClassifier常用参数说明如下:

(1) critrion : 特征选择的准则,"gini" or "entropy"(default = "gini"),即选择基尼系数还是信息增益;

(2) splitter : 特征划分标准,"best" or "random" (default = "best"),即选择全局最优划分还是局部最优划分;

(3) max_depth : 树的最大深度,当特征和数据较少时,可以使用默认(default = None),一般设为5-20,深度过大易出现过拟合;

(4) min_samples_split : 节点的再划分需要的最小样本数(default = 2);

(5) min_impurity_decrease :节点再划分的最小不纯度(default = 0);

(6) min_samples_leaf :叶子节点最小样本数(default = 1);

(7) random_state:随机状态,可以设定一个值,保证每次运行的结果是一样的;

(8) class_weight:指定各样本的类别权重,"balanced" or None(default = None),或者自行设定,当样本不均衡时,常用"balanced"平衡。

这里,由于特征个数并不是很多,所以我们设置树的最大深度为20,由于患病和不患病两种样本的比例约为3:1,样本本身不均衡,所以需设class_weight = "balanced",建模结果如下:

代码语言:javascript复制
## 决策树
from sklearn.tree import DecisionTreeClassifier
dt = DecisionTreeClassifier(criterion = 'gini',max_depth = 20,random_state = 50,class_weight = "balanced").fit(x_train,y_train)
y_pred = dt.predict(x_test)
evaluate("DecisionTree",y_test,y_pred)

从决策树的分类结果来看,F1值为0.68,预测正确率为67.6%,与逻辑回归的预测性能几乎相同。

在集成方法中,主要有两种:Bagging和Boosting.

(1) Bagging方法的思想为 :对训练样本进行M次Bootstrap抽样,得到M个样本容量为n的训练集,对每个训练集分别训练一个基分类器,各分类器之间相互独立,最后组合分类器的分类结果,采取“多数表决”的原则得到最终的样本分类结果。Bagging的典型代表方法为随机森林,其主要思想为:通过每次随机抽取特征的方式来建立M棵相互独立的树,再将各树的结果进行平均得到最终结果。

(2) Boosting方法的思想为 : 训练M个基分类器,各分类器之间相互关联,以一定的方式进行组合。Boosting的典型代表方法为AdaBoost,其主要思想为:每次提升前一个分类器中分类错误样本的权值,降低分类正确的样本权值;在进行各分类器的组合时,采用“加权多数表决”的原则,加大分类错误率小的分类器的权值,减小分类错误率大的分类器的权值。

3.3 随机森林

在Python中使用sklearn.ensemble中的RandomForestClassifier类进行随机森林建模,其主要参数包括:

(1) n_estimators : 森林中树的数量 (default = 10);

(2) max_features :树的最大特征数(default = "auto"),因为随机森林的一大特点是每棵树随机抽取特征,"auto"表示sqrt(n_features), "sqrt"(同"auto"),"log2"表示log2(n_features),None表示max_features = n_features;

(3) bootstrap : 是否选用bootstrap来进行抽取数据集进行建树(default = True); (4) 其它类似于决策树中的参数。

代码语言:javascript复制
## 随机森林
from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators = 100, max_depth = 20,max_features = "auto", random_state = 20).fit(x_train,y_train)
y_pred = rf.predict(x_test)
evaluate("RandomForest",y_test,y_pred)

从随机森林的分类结果来看,F1值为0.72,预测正确率为72.4%,相较于单棵的决策树和逻辑回归方法,预测正确率提升了近5个百分点,足以体现集成方法的优势。

3.4 AdaBoost

在Python中使用sklearn.ensemble的AdaBoostClassifier类进行分类建模,其主要参数包括:

(1) base_estimator : 基分类器,默认是决策树,最大深度为1;

(2) n_estimators : 训练分类器的数量(default = 50);

(3) learning_rate : 学习率(default = 1);

(4) algorithm : Boosting算法,也就是模型提升准则,有两种方式SAMME, 和SAMME.R(default = SAMME.R),二者的区别主要是弱学习器权重的度量方式不同,前者是对样本集预测错误的概率进行划分的,后者是对样本集的错分率进行划分的;

(5) random_state : 随机状态,可以设定一个值,保证每次运行的结果是一样的。

代码语言:javascript复制
from sklearn.ensemble import AdaBoostClassifier
ab = AdaBoostClassifier(n_estimators = 200,learning_rate = 1,algorithm = 'SAMME.R', random_state = 20).fit(x_train,y_train)
y_pred = ab.predict(x_test)
evaluate("AdaBoost",y_test,y_pred)

从AdaBoost分类结果来看,F1值为0.73,预测正确率为73.8%。正确率相较于决策树提升近6个百分点,相较于随机森林提升近1个百分点,预测效果最好。

3.5 树模型特征选择

由于树模型可以自行做特征选择,在树的生成过程中,每次将会选择对数据最有分类能力的特征,下面将分别查看3种树模型中各个特征的重要性。

代码语言:javascript复制
model = [dt,rf,ab]
name = ["Decision Tree","Random Forest","AdaBoost"]
fig,ax = plt.subplots(1,3,figsize = (15,6))
for i in range(3):
    sub = int("13" str(i 1))
    x = list(patient.columns[:-1])
    y = list(model[i].feature_importances_)
    plt.subplot(sub,sharey = ax[0])
    plt.bar(x,y) 
    plt.xlabel("Feature",fontsize = 13)
    plt.ylabel("Feature_importance",fontsize = 13)
    plt.title("The features' importance of %s "%name[i])  
plt.tight_layout()
plt.show()   

上图分别展示了3种树模型的特征重要性,不同模型中,各个特征的重要程度不同。其中,决策树模型中总胆红素(TBIL)的重要性最高,年龄(Age)第二;随机森林模型中三种血清酶的重要性最高;AdaBoost模型中,碱性磷酸酶(ALP)谷丙转氨酶(ALT)的重要性较高。事实上,树模型在选择最优特征建树时,并没有考虑特征之间的相关性,故其选出的重要特征中可能具有高相关性,这也是树模型的一大缺点。另外,观察到,性别(Gender)在所有模型中的重要性均最低,这说明它对患者肝病诊断的帮助最小。下面考虑将该特征从模型中剔除,重新建立树模型。

代码语言:javascript复制
del x_train["Gender"]
del x_test["Gender"]
代码语言:javascript复制
dt1 = DecisionTreeClassifier(criterion = 'gini',max_depth = 20,random_state = 50,class_weight = "balanced").fit(x_train,y_train)
y_pred = dt1.predict(x_test)
print("-------------The result of DecisionTree-------------n")
evaluate("DecisionTree",y_test,y_pred)

rf1 = RandomForestClassifier(n_estimators = 100, max_depth = 20,max_features = "auto", random_state = 20).fit(x_train,y_train)
y_pred = rf1.predict(x_test)
print("n-------------The result of RandomForset-------------n")
evaluate("RandomForest",y_test,y_pred)

ab1 = AdaBoostClassifier(n_estimators = 200,learning_rate = 1,algorithm = 'SAMME.R', random_state = 20)
ab1.fit(x_train,y_train)
y_pred = ab1.predict(x_test)
print("n---------------The result of AdaBoost---------------n")
evaluate("AdaBoost",y_test,y_pred)

观察到,在本文使用的三种树模型中,删除性别(Gender)后,决策树模型的F1值上升0.02,分类正确率上升2%,整体性能提升;随机森林模型的F1值减小0.05,分类正确率下降5%,整体性能下降,这可能是由于随机森林在建树时随机选择特征带来的影响;AdaBoost模型的性能保持不变,这是由于性别(Gender)在AdaBoost中的重要性几乎为0,当基分类器的最大深度为1时,建树基本上不会选到性别(Gender),故新模型与原模型结果一致。另外,由该结果可知,并不是所有的集成方法都会优于单一的分类器,集成的结果与特征选择、基分类器的组成、参数设定等多个方面有关。所以在建立预测模型时,需要结合各个方面来选择合适的分类器。

4.总结

本案例基于肝病诊断数据进行了分析和预测,通过数据预处理、探索性分析和分类建模三个方面对数据进行了深入挖掘。在分类建模过程中,分别使用了单一分类器和集成方法进行预测,通过对比F1值和预测正确率对模型进行评估,AdaBoost的预测效果最好。另外,从模型的预测结果来看,年龄、总胆红素、三种血清酶等特征对肝病的诊断起着重要作用,因而在实际肝病诊断时,可重点参考这些指标。

0 人点赞