如何使用python sklearn 中的LeaveOneGroupOu？特例：leave one pair out的实现多重比较校正？

引言

一般我们会在研究中验证某个特征的有效性，会使用该特征对病例个样本进行二分类，通常回选择支持向量机（SVM）,这里的目的并非一定要得到很高的正确率，更多地只是想说明该特征可以区分case组和control组。

留一法 or 留一对法？

因为样本量的原因，一般是使用留一法进行多重比较校正，但是有时候也会被要求使用留一对的方法，简单的就是每一次从case中拿出一例，然后从control中拿出一例，使用剩余的（m+n）- 2的样本进行训练，对取出的样本（一个来自case一个来自control）进行测试。

但是你可能会想case组有m个样本，control组有n个样本，两组人的数量不一样（m≠n）要怎么处理呢？

首先k=m<n?m:n（选取二者中较小的一个数）
一共进行k次随机取样本，训练模型，分类，得到分类正确率
然后重复步骤（1,2）N次，就会得到多个准确率和多条ROC曲线。
你可能会问为什么要重复这多（N）次？因为这m和n不相等，进行k次抽样，肯定会导致有一些样本（样本量大的那一组有样本无法被抽中成为测试集），因此通过重复N次来较少这个因素。

具体如何实现留一对？leave one pair out?

在我以前的博客中已经介绍了如何使用留一法的SVM,具体见留一法SVM，可以借助sklearn.model_selection中的LeaveOneGroupOut实现，该函数可以根据组别（可以理解为每一个样本有一个所属组的编号）将数据分成与组别相同的组数，比如我有20个样本，这20个样本根据性别可以分为2组。在后续的多重比较校正中每一次取出其中的一组作为测试集，剩余组作为训练集即可。leave one pair out 是上述过程的一个特例（因为在对两组人进行分类时，可能会有需求每次将类别1和类别2同时拿出来一个作为测试集），因此接祖leaveOneGroupOut,我们只需要对两组人进行分别编号，然后编号对应时作为一组，具体实现过程如下：

# 该方法是每次排除一个组，怎样才能每次排除一对？那就是正常组和病人组对应配对？
pair_data1 = [i for i in range(np.size(data1_label,0))]
pair_data2 = [i for i in range(np.size(data2_label,0))]
# 随机打乱分组
random.shuffle(pair_data1)
random.shuffle(pair_data2)
pair = pair_data1+pair_data2
loo = LeaveOneGroupOut()
loo.get_n_splits(dataset,groups=pair)
...
for train_index, test_index in loo.split(dataset,groups=pair):X_train, X_test = dataset[train_index], dataset[test_index]Y_train, Y_test = np.array(datalabels)[train_index], np.array(datalabels)[test_index]...

上述代码说明

首先将data1进行编号0：M-1,然后对N进行编号0：N-1，数字为组别
之后对组别进行打乱，打乱的组别（后续说明为什么打乱）
然后将分组编号合并（前面为data1的组别，后面为data2的组别）
然后将样本特征及其对应的分组编号作为参数输入到LeaveOneGroupOut的get_n_splits中，其实上述过程会将数据分成max(M，N)份，如果M与N的大小不一样的话，整个过程会出现leave one out 的情况（后续会对此进行处理）

特殊处理说明

因为编号是由小到大，因此当编号小于min(M，N)时程序都能正常的leave one pair out，当编号大于等于min(M，N)时只能得到一个一个样本，此时在代码中加入判断条件跳出循环，就可以了。代码如下：

print(test_index)if np.size(X_test,0)==1:print('条件不满足，跳出循环！')break

为什么需要随机打乱？因为两组人的数量可能不同，分组的方式也不唯一，因此需要随机打乱来排除这种认为分组方式对结果的影响，因此每次代码的运行分组的方式都是不一样的，可以通过多运行几次来减小分组带来的误差，同时可以借此判断结果是否稳定。

完整实现代码


"""
Created on Wed Sep  9 10:06:14 2020@author: Administer
"""import numpy as np
from sklearn.model_selection import LeaveOneOut
from sklearn.model_selection import LeaveOneGroupOut
import scipy.io as scio
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
import random###################功能性函数###################################################
# 将一个任意嵌套的列表整理为一个列表
def flatten(nested):try:for sublist in nested:for element in flatten(sublist):yield elementexcept TypeError:yield nested# 读取一个txt文件中的数字并将其转换为一个列表
def read_txt2number(path):significant_index = []with open(path) as f:for line in f.readlines():line = line.strip('\n')      #去掉列表中每一个元素的换行符significant_index.append(line)number_list = [int(i) for i in significant_index]   # 将字符类型转换为数字return number_list  # 获取一个预测结果的TP TN FP FN
def get_TpTN_FpFn(list1,list2):# list1 为真实的label list2 为预测的labelreallabel_list=list(flatten(list1))predictlabel_list=list(flatten(list2))TP_count = 0TN_count = 0FP_count = 0FN_count = 0for i in range(len(reallabel_list)):if reallabel_list[i] == 1 and predictlabel_list[i] ==1:TP_count += 1if reallabel_list[i] == -1 and predictlabel_list[i] ==-1:TN_count += 1if reallabel_list[i] == -1 and predictlabel_list[i] == 1:FP_count +=1if reallabel_list[i] == 1 and predictlabel_list[i] == -1:FN_count += 1return TP_count, TN_count, FP_count, FN_count
##################数据的读取与整理了#############################################path2 = r'C:\Users\Administer\Desktop\GTCS_HC\20201013data_and_result\6_gene expression class_SVM\SVM class\dataset'MSN_HC = scio.loadmat(path2 + '\\average_MSr_dataset1_HC.mat')
train_MSN_HC = MSN_HC['result']
MSN_GTCS = scio.loadmat(path2 + '\\average_MSr_dataset1_GTCS.mat')
train_MSN_GTCS = MSN_GTCS['result']train_data1 = train_MSN_HC;
train_data2 = train_MSN_GTCS;
dataset =(np.hstack((train_data1,train_data2))).Tdata1_label = list(-1 for i in range(np.size(train_data1,1)));
data2_label = list(1 for i in range(np.size(train_data2,1)))
datalabels = data1_label + data2_label#####################################LOOVC#####################################
# loo = LeaveOneOut()
# loo.get_n_splits(dataset)# 该方法是每次排除一个组，怎样才能每次排除一对？那就是正常组和病人组对应配对？
pair_data1 = [i for i in range(np.size(data1_label,0))]
pair_data2 = [i for i in range(np.size(data2_label,0))]
# 随机打乱分组
random.shuffle(pair_data1)
random.shuffle(pair_data2)
pair = pair_data1+pair_data2
loo = LeaveOneGroupOut()
loo.get_n_splits(dataset,groups=pair)predictlabel_list = []
reallabel_list = []
coef_weight_list = []
Y_score_list = []
sum_of_coef = np.zeros((1,np.size(dataset,1)))
#clf = make_pipeline(StandardScaler(), SVC(C=0.9, kernel='linear', gamma='auto'))
#clf = SVC(C=0.9, kernel='linear', gamma='auto')
#clf = svm.LinearSVC()
scaler = StandardScaler()
dataset = scaler.fit_transform(dataset)
clf = SVC(C=1, kernel='linear', gamma='auto')
count_right_label = 0
count = 0                     # 循环次数
# 用留一法进行验证for train_index, test_index in loo.split(dataset,groups=pair):X_train, X_test = dataset[train_index], dataset[test_index]Y_train, Y_test = np.array(datalabels)[train_index], np.array(datalabels)[test_index]print(test_index)if np.size(X_test,0)==1:print('条件不满足，跳出循环！')breakclf.fit(X_train,Y_train)Y_score_temp = clf.decision_function([X_test[0,:]])            # 得到的结果为该类到超平面的距离temp1 = clf.coef_Y_score_temp1 = clf.decision_function([X_test[1,:]])            # 得到的结果为该类到超平面的距离temp2 = clf.coef_sum_of_coef += np.abs(temp1)sum_of_coef += np.abs(temp2)coef_weight_list.append(temp1)coef_weight_list.append(temp2)predictlabel_list.append(clf.predict([X_test[0,:]]))predictlabel_list.append(clf.predict([X_test[1,:]]))reallabel_list.append(Y_test[0])reallabel_list.append(Y_test[1])Y_score_list.append(Y_score_temp)Y_score_list.append(Y_score_temp1)if Y_test[0] == clf.predict([X_test[0,:]]):count_right_label +=1 if Y_test[1] == clf.predict([X_test[1,:]]):count_right_label +=1 count += 1print('第{}次循环'.format(count))
accurancy = count_right_label/len(reallabel_list)
print('******循环结束！************')
print('准确率为：%.2f%%' %(accurancy*100))
print('******运行结束！************')average_weight = list(flatten(sum_of_coef/2*count)) TP_count, TN_count, FP_count, FN_count = get_TpTN_FpFn(reallabel_list,predictlabel_list)
F1_score = (2 * TP_count)/(2*TP_count + FP_count +FN_count)
F2_score = (2 * TN_count)/(2*TN_count + FP_count +FN_count)
ACC = (TP_count + TN_count)/(TP_count + FN_count + TN_count + FP_count)
SEN = TP_count/(TP_count + FN_count)
SPE = TN_count / (TN_count + FP_count)print('F1_SCORE为：%.2f%%' %(F1_score*100))
print('F2_SCORE为：%.2f%%' %(F2_score*100))
print('ACC(准确率)为：%.2f%%' %(ACC*100))
print('SEN（敏感度）为：%.2f%%' %(SEN*100))
print('SPE(特异性)为：%.2f%%' %(SPE*100))  #print('ROI权重为：{}'.format(average_weight))import matplotlib.pyplot as plt# ROC AUCfrom sklearn.metrics import roc_curve, auc
# 为每个类别计算ROC曲线和AUC
real = np.array(list(flatten(reallabel_list)))
predict = np.array(list(flatten(predictlabel_list)))
roc_auc = dict()
fpr, tpr, threshold = roc_curve(real,Y_score_list)
roc_auc = auc(fpr, tpr)plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange',lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
order=10
# plt.savefig('./loocv_ROC_roc'+str(order)+'.png',dpi=400) # 以400大批保存图片
plt.show()
# scio.savemat('fpr_fdr'+str(order)+'.mat',{'fpr':fpr,'tpr':tpr} )