Python基于逻辑回归的糖尿病视网膜病变检测（数据集messidor_features.arff）

一. 引言
本项目基于逻辑回归理论，运用Python语言对数据集messidor_features.arff进行分析，实现对糖尿病视网膜病变的检测。糖尿病视网膜病变（DR）是糖尿病最常见的微血管并发症之一，是慢性进行性糖尿病导致的视网膜微血管渗漏和阻塞从而引起一系列的眼底病变，如微血管瘤、硬性渗出、黄班水肿甚至视网膜脱离。DR检测对于糖尿病人群筛查、糖尿病患者早期治疗具有重要意义。

二. 数据集描述
1. 下载地址messidor_features.arff
2. 数据集开头有一些描述信息，训练是用不到的，我选择删掉，方便用pd.read_csv()函数读取。

当然不删也是可以的，有专门的函数读取.arff文件

from scipy.io import arffimport pandas as pddf = arff.loadarff("messidor_features.arff") #读取出来是一个元组dataframe = pd.DataFrame(df[0])

3.数据集messidor_features.arff包含从 Messidor 图像集中提取的特征，用于预测图像是否包含糖尿病视网膜病变的迹象。所有特征都代表检测到的病变、解剖部位的描述特征或图像级描述符。该数据集有20条属性，类标签是最后一条，如图。

4. 对数据集有大概了解后，来简单探索性分析一下它。首先查看data.info(),看看有没有缺失值及数据类型。

import pandas as pdpath="E:/Python_file/zuoye/messidor_features.arff"Cnames = ["x0", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9",      "x10", "x11", "x12", "x13", "x14", "x15", "x16", "x17", "x18", "y"]#删掉与数据集无关的内容data=pd.read_csv(path,header=None,names=Cnames)print("数据集基础信息：")print(data.info())

可以看到数据里没有缺失值。再统计一下，类标签为‘1’的数据有611条，占比为53.1%；标签为‘0’的数据有540条，占比为46.9%。说明正例数据与负例数据的数量分布是均衡的。
5. 这里推荐使用pandas-profiling库，可以一键生成对数据集的分析报告，非常好用。没有安装的话用pip install pandas_profiling 命令安装一下。

import pandas_profilingreport= pandas_profiling.ProfileReport(data)report.to_file("output_file.html")

运行后生成一个可交互的.html文件，通常包含对数据的类型检测；计算唯一值、缺失值；分位数统计如最小最大值、四分位数、中位数等；描述统计如平均数、众数、峰度偏度等；变量间相关系数的热力图，等等。
我们来看看该数据集下各属性间用Spearman秩相关系数（ρ）描述的单调相关的度量热力图，ρ的值介于-1和+1之间，-1表示完全负单调相关，0表示没有单调相关，1表示完全正单调相关。

三. 方法介绍
逻辑回归的原理有很多博主写的很好，我就不班门弄斧了，主要说说代码。利用sklearn库提供的LogisticRegression()可以很方便的完成训练和预测。

from sklearn.linear_model import LogisticRegressionfrom sklearn.model_selection import train_test_splitX=data[data.columns[0:19]] #提取特征，不要标签y=data["y"]                #train_size=0.8，80%的训练集占比x_train,x_test,y_train,y_test=train_test_split(X,y,train_size=0.8,random_state=90)lr=LogisticRegression(max_iter=3000)clm=lr.fit(x_train,y_train)print("对测试集的预测结果：")#输出预测结果、预测结果的结构类型及尺寸result=clm.predict(x_test)

LogisticRegression()的参数很多，但需要设置的不多。我们的数据集是分布均衡的，参数类别权重class_weight不需要设置；关于参数优化算法solver用默认的‘liblinear’就好，因为这是二分类问题（只看有没有病变），而且我们是小数据集，也用不到面向大数据集的‘sag’和‘saga’；而max_iter是设置迭代次数，如果小了，可能模型没收敛就运行结束了，这里我设置为3000次。

四. 结果和模型评价及可视化
1.测试集的检测结果如图，1代表有病变，0代表没有。要注意的是划分训练集和测试集时，random_state等于不同的值，会得到不同的测试集，我这里是random_state=90，改成其他数，预测结果就和我不同，但对模型评价没有影响。

2. 模型评价的指标有很多，比如召回率、精度、准确率、F统计量、决定系数R²、ROC曲线的包络面积AUC等。可以用classification_report()一键生成评估报告。

from sklearn.metrics import classification_reportprint("性能报告;")print(classification_report(y_test,result))confusion = metrics.confusion_matrix(y_test, result)

3. 以FPR为横轴，TPR为纵轴，绘制ROC曲线，并由曲线计算得AUC=0.77。

from sklearn.metrics import roc_curve,aucfrom matplotlib import pyplot as pltfpr, tpr, thr = roc_curve(y_test, result, drop_intermediate=False)fpr, tpr = [0] + list(fpr), [0] + list(tpr)plt.plot(fpr, tpr)plt.title("ROC curve for diabetes classifier")plt.xlabel("False Positive Rate (1 - Specificity)")plt.ylabel("True Positive Rate (Sensitivity)")plt.grid(True)plt.show()print("AUC:"+ str(auc(fpr,tpr)))

该模型的检测能力还是可以的。
4. 从逻辑回归模型中导出各个变量的回归系数，由此作出重要程度的条形图。

print("逻辑回归各变量系数：")print(clm.coef_)coef_lr = pd.DataFrame({"var" : x_test.columns,                        "coef" : clm.coef_.flatten()                        })index_sort =  np.abs(coef_lr["coef"]).sort_values().indexcoef_lr_sort = coef_lr.loc[index_sort,:]# 水平柱形图绘图fig,ax=plt.subplots()x, y = coef_lr_sort["var"], coef_lr_sort["coef"]rects = plt.barh(x, y, color="dodgerblue")plt.grid(linestyle="-.", axis="y", alpha=0.4)plt.tight_layout()#添加数据标签for rect in rects:    w = rect.get_width()    ax.text(w, rect.get_y()+rect.get_height()/2,"%.2f" %w,ha="left",va="center")

逻辑回归就是把线性回归的结果输入到了sigmoid函数，所以各变量的系数还是有的。

可以看到对视网膜病变检测结果影响最大的五个变量分别是x14,x1,x2,x0,x15，它们的含义参照上文。可以确定这五个变量是导致糖尿病视网膜病变的主要因素。

五. 完整代码

import pandas as pdfrom sklearn.linear_model import LogisticRegressionfrom sklearn.model_selection import train_test_splitfrom sklearn.metrics import accuracy_score, precision_score, recall_score, f1_scorefrom sklearn.metrics import classification_reportfrom sklearn import metricsfrom sklearn.metrics import roc_curve,aucfrom matplotlib import pyplot as pltimport numpy as nppath="/content/drive/MyDrive/messidor_features.arff"Cnames = ["x0", "x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9",      "x10", "x11", "x12", "x13", "x14", "x15", "x16", "x17", "x18", "y"]#首先删掉与数据集无关的内容data=pd.read_csv(path,header=None,names=Cnames)X=data[data.columns[0:19]] y=data["y"]x_train,x_test,y_train,y_test=train_test_split(X,y,train_size=0.8,random_state=90)lr=LogisticRegression(max_iter=3000)clm=lr.fit(x_train,y_train)print("对测试集的预测结果：")#输出预测结果、预测结果的结构类型及尺寸result=clm.predict(x_test)print(result,type(result),result.shape) print("模型评分："+ str(clm.score(x_test,y_test))) #用决定系数来打分print("性能报告;")print(classification_report(y_test,result))confusion = metrics.confusion_matrix(y_test, result)print("混淆矩阵：")print(confusion)TP = confusion[1, 1]TN = confusion[0, 0]FP = confusion[0, 1]FN = confusion[1, 0]print("TPR(正确识别的正例数据在实际正例数据中的占比):"+str(TP/(TP + FN)))print("TNR(正确识别的负例数据在实际负例数据中的占比):"+str(TN/(TN + FP)))print("Accuracy score: ", format(accuracy_score(y_test, result)))print("Precision score: ", format(precision_score(y_test, result)))print("Recall score: ", format(recall_score(y_test, result)))print("F1 score: ", format(f1_score(y_test, result)))fpr, tpr, thr = roc_curve(y_test, result, drop_intermediate=False)fpr, tpr = [0] + list(fpr), [0] + list(tpr)plt.plot(fpr, tpr)plt.title("ROC curve for diabetes classifier")plt.xlabel("False Positive Rate (1 - Specificity)")plt.ylabel("True Positive Rate (Sensitivity)")plt.grid(True)plt.show()print("AUC:"+ str(auc(fpr,tpr)))print("逻辑回归各变量系数：")print(clm.coef_)coef_lr = pd.DataFrame({"var" : x_test.columns,                        "coef" : clm.coef_.flatten()                        })index_sort =  np.abs(coef_lr["coef"]).sort_values().indexcoef_lr_sort = coef_lr.loc[index_sort,:]# 水平柱形图绘图fig,ax=plt.subplots()x, y = coef_lr_sort["var"], coef_lr_sort["coef"]rects = plt.barh(x, y, color="dodgerblue")plt.grid(linestyle="-.", axis="y", alpha=0.4)plt.tight_layout()#添加数据标签for rect in rects:    w = rect.get_width()    ax.text(w, rect.get_y()+rect.get_height()/2,"%.2f" %w,ha="left",va="center")

如果有错，还望指正。