3.Machine Learning with Python
scikit-learn等python package对很多机器学习模型提供了易于使用的接口。但是,如果想用机器学习方法解决实际问题,我们需要考虑的步骤不只是拟合模型,还包括数据预处理,特征选择,调参,模型评估等等。我们只有对数据有比较深入理解,才能选择合适的方法。本节中我们提供一个用基于python的机器学习工具解决实际的分类问题的例子。我们希望读者读者结合本章内容,能通过实践掌握相关工具的使用。

1) 数据说明

我们这里用到的还是BreastCancer.csv数据集。
  1. 1.
    文件包括11个列,第1列为样本编号,第2-10列为特征,11列为标签(benign为良性,malignant为恶性)
  2. 2.
    数据集有458个良性(benign)样本和241个恶性(malignant)样本
  3. 3.
    数据集有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
1000025
5
1
1
1
2
1
3
1
1
benign
1002945
5
4
4
5
7
10
3
2
1
benign
1015425
3
1
1
1
2
2
3
1
1
benign
1016277
6
8
8
1
3
4
3
7
1
benign
1017023
4
1
1
3
2
1
3
1
1
benign
1017122
8
10
10
8
7
10
9
7
1
malignant

2) Load python packages

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) 数据预处理

接下来我们进行数据预处理:
  • 读取数据
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。
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类如何获取特征重要性的信息
    • step参数控制每轮筛选移除多少个特征
    • RFECV对象的support_属性记录了有哪些特征被保留了下来,ranking_记录了特征在倒数第几轮筛选中被移除。最终保留的特征ranking_为1。
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只会保留一部分特征:
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个特征的组合
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类,把能使用到的特征的标号也作为一个超参数:
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评估每个特征组合的性能:
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) 评估模型效果

  • 计算AUROC
y_pred_proba = clf.predict_proba(X_validation)[:,1]
fpr, tpr,_ = roc_curve(y_validation,y_pred_proba)
AUROC = auc(fpr, tpr)
  • 绘制ROC曲线
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()
png

7) Homework

7.1

我们提供了一个qPCR数据集qPCR_data.csv,第1列为sample id,第2-12列为特征(11个基因的表达量),第13列为样本标签(负例为健康人:NC,正例为肝癌病人:HCC)。请大家完成以下任务:
  • 数据预处理
    • 用均值或中位数补全空缺值
    • 必要时对数据进行scaling
  • 数据集划分:预留20%数据用于评估模型泛化能力,剩下的用于模型拟合
  • 模型选择和模型拟合
    • 任选一种分类器即可
    • 调参: 根据所选的分类器对相应的超参数进行搜索
    • 特征选择: 用RFE/穷举或其他方式均可,特征数量不限
  • 模型评估
    • 计算预留数据集上的AUROC
    • 绘制ROC曲线
  • 请提交代码,必要的文字解释和ROC曲线

7.2

随机森林是生物信息经常使用的一个分类器。请大家查阅资料,回答以下两个问题:
  • 随机森林中树的数量是不是一个需要通过交叉验证调整的超参数?为什么?
  • 请问什么是随机森林的out-of-bag (OOB) error?它和bootstrapping有什么关系?