scikit-learn 等python package对很多机器学习模型提供了易于使用的接口。但是，如果想用机器学习方法解决实际问题，我们需要考虑的步骤不只是拟合模型，还包括数据预处理，特征选择，调参，模型评估等等。我们只有对数据有比较深入理解，才能选择合适的方法。本节中我们提供一个用基于python的机器学习工具解决实际的分类问题的例子。我们希望读者读者结合本章内容，能通过实践掌握相关工具的使用。
1) 数据说明
我们这里用到的还是BreastCancer.csv 数据集。
文件包括11个列，第1列为样本编号，第2-10列为特征，11列为标签（benign为良性,malignant为恶性)
数据集有458个良性（benign）样本和241个恶性（malignant）样本
数据集有Cl.thickness，Cell.size，Cell.shape，Marg.adhesion，Epith.c.size，Bare.nuclei，Bl.cromatin，Normal.nucleoli，Mitoses9个特征，每个特征取值在为0-10之间的整数。
样例数据如下：
Id
Cl.thickness
Cell.size
Cell.shape
Marg.adhesion
Epith.c.size
Bare.nuclei
Bl.cromatin
Normal.nucleoli
Mitoses
Class
2) Load python packages
Copy import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_curve,auc LogisticRegression是sklearn为逻辑斯蒂回归实现的一个类
GridSearchCV类是sklearn对于交叉验证调参过程的封装，可以简化编程实现
3) 数据预处理
接下来我们进行数据预处理:
Copy data = pd.read_csv('BreastCancer.csv',sep=',',index_col=0)
label_lut = pd.Series({"benign":0,"malignant":1})
y = label_lut.loc[data["Class"]].values # 将样本标签由字符串转化为整数表示
data = data.iloc[:,:-1] #选取feature对应的列
data = data.fillna(data.mean(axis=0)) #用均值填充缺失值
X = data.values
X = StandardScaler().fit_transform(X) # Z score scaling 4) 数据集划分
这里我们首先使用train_test_split 方法来随机的将20%的数据预留为一个用于最终评估模型泛化能力的数据集，剩下80%数据为可以见到label的Discovery set。
Copy X_discovery, X_validation, y_discovery, y_validation = train_test_split(X, y, test_size=0.2, random_state=666)
print(f"{X_discovery.shape[0]} samples in discovery set")
print(f"{X_validation.shape[0]} samples in validation set")
# 559 samples in discovery set
# 140 samples in validation set 5) 特征选择
简单起见我们这里只考虑logistic regression。其他分类器的实现也非常类似。
我们这里用两种方式进行特征选择
recursive feature elimination (RFE): 用交叉验证评估特征在discovery set上的分类效果，可以用sklearn提供的RFECV类实现
exhaustive enumeration: 枚举所有三个特征的组合，挑选出discovery set上交叉验证效果最好的组合。我们的数据集只有9个特征，所有三个特征的组合总共有84种，这还是可以接受的。当特征数目非常多事，使用这种穷举法几乎是不可能的。
5.1) RFE
GridSearchCV是sklearn封装的一个用于调参的类
为了初始化这个类，我们需要提供分类器以及包含调参中参数取值范围的字典param_grid
我们可以像一个普通的分类器一样使用它，但是它的内部会对我们提供的每一个超参数组合进行交叉验证，自动选出最好的超参数组合，并用所有的输入数据和最好的超参数重新拟合模型(refit参数默认为True)
最好的模型被保存在best_estimator_属性中。
RFECV是sklearn封装的一个用RFE做特征选择的类
我们同样可以像一个普通的分类器一样使用它，但是它的内部会根据拟合的分类器提供的coef_或feature_importances_，每次去除一些最不重要的特征，直到交叉验证的性能下降
为了初始化这个类，我们需要提供一个分类器，这个分类器需要有coef_或feature_importances_属性，如果没有则需要定义一个函数传给importance_getter参数，告诉RFECV类如何获取特征重要性的信息
RFECV对象的support_属性记录了有哪些特征被保留了下来，ranking_记录了特征在倒数第几轮筛选中被移除。最终保留的特征ranking_为1。
Copy clf = GridSearchCV(LogisticRegression(penalty="l2"),
param_grid={"C":[0.01,0.1,1,10]}, cv=5,
scoring="roc_auc",refit=True,verbose=4)
selector = RFECV(clf, step=1,
min_features_to_select=3,cv=5,
importance_getter=lambda clf:clf.best_estimator_.coef_,
scoring="roc_auc",verbose=4)
selector = selector.fit(X_discovery, y_discovery)
print(selector.support_)
# array([ True, True, True, True, True, True, True, True, True])
print(selector.ranking_)
# [1 1 1 1 1 1 1 1 1] 这里我们用到的LogisticRegression回归默认使用l2正则化，相当于自带特征选择功能，所以所有的特征都被保留了下来
如果我们人为使用比较弱的正则化(把超参数C的搜索空间限制到比较大的值)，我们会发现RFECV只会保留一部分特征:
Copy clf = GridSearchCV(LogisticRegression(penalty="l2"),
param_grid={"C":[100,1000]}, cv=5,
scoring="roc_auc",refit=True,verbose=4)
selector = RFECV(clf, step=1,
min_features_to_select=3,cv=5,
importance_getter=lambda clf:clf.best_estimator_.coef_,
scoring="roc_auc",verbose=4)
selector = selector.fit(X_discovery, y_discovery)
print(selector.support_)
# [ True False True True False True True False True]
print(selector.ranking_)
# [1 2 1 1 4 1 1 3 1] 5.2) exhaustive enumeration
我们首先枚举出所有的3个特征的组合
Copy from itertools import combinations
feature_combinations=[]
for i in combinations(list(range(9)), 3):
feature_combinations.append(list(i))
print(len(feature_combinations))
# 84 然后我们定义一个MaskedLogisticRegression类，把能使用到的特征的标号也作为一个超参数:
Copy from sklearn.base import ClassifierMixin, BaseEstimator
class MaskedLogisticRegression(BaseEstimator, ClassifierMixin):
def __init__(self,feature_indices=None,**params):
self.feature_indices = feature_indices
self.estimator = LogisticRegression(**params)
def mask(self,X):
if self.feature_indices is None:
return X
else:
return X[:,self.feature_indices]
def fit(self, X, y=None):
self.classes_ = np.unique(y)
return self.estimator.fit(self.mask(X),y)
def predict(self, X):
return self.estimator.predict(self.mask(X))
def predict_proba(self, X):
return self.estimator.predict_proba(self.mask(X)) 用GridSearchCV评估每个特征组合的性能:
Copy clf = GridSearchCV(MaskedLogisticRegression(),
param_grid={"feature_indices":feature_combinations}, cv=5,
scoring="roc_auc",refit=True,verbose=4)
clf = clf.fit(X_discovery, y_discovery)
print(list(data.columns[clf.best_params_['feature_indices']]))
# ['Cl.thickness', 'Cell.shape', 'Bare.nuclei'] 6) 评估模型效果
Copy y_pred_proba = clf.predict_proba(X_validation)[:,1]
fpr, tpr,_ = roc_curve(y_validation,y_pred_proba)
AUROC = auc(fpr, tpr)
Copy plt.figure(figsize=(4,4))
plt.plot(fpr, tpr, '-', color='b', label='Validation AUC of {:.4f}'.format(AUROC), lw=2)
plt.plot([0, 1], [0, 1], '--', color=(0.6, 0.6, 0.6), label='Random Chance')
plt.xlim([-0.01, 1.01])
plt.ylim([-0.01, 1.01])
plt.title('ROC curve of test data')
plt.xlabel('FPR')
plt.ylabel('TPR')
plt.legend(loc='best',fontsize='small')
plt.tight_layout()
plt.show()
plt.close() 7) Homework
7.1)
我们提供了一个qPCR数据集qPCR_data.csv ，第1列为sample id，第2-12列为特征(11个基因的表达量)，第13列为样本标签(负例为健康人:NC,正例为肝癌病人:HCC)。请大家完成以下任务:
数据集划分：预留20%数据用于评估模型泛化能力，剩下的用于模型拟合
模型选择和模型拟合
特征选择: 用RFE/穷举或其他方式均可，特征数量不限
请提交代码，必要的文字解释，数据可视化结果及ROC曲线
7.2)
随机森林是生物信息经常使用的一个分类器。请大家查阅资料，回答以下两个问题:
请问什么是随机森林的out-of-bag (OOB) error?它和bootstrapping有什么关系?
Last updated 8 months ago
