3.Machine Learning with Python
scikit-learn等python package对很多机器学习模型提供了易于使用的接口。但是，如果想用机器学习方法解决实际问题，我们需要考虑的步骤不只是拟合模型，还包括数据预处理，特征选择，调参，模型评估等等。我们只有对数据有比较深入理解，才能选择合适的方法。本节中我们提供一个用基于python的机器学习工具解决实际的分类问题的例子。我们希望读者读者结合本章内容，能通过实践掌握相关工具的使用。
- 1.文件包括11个列，第1列为样本编号，第2-10列为特征，11列为标签（benign为良性,malignant为恶性)
- 2.数据集有458个良性（benign）样本 和241个恶性（malignant）样本
- 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 |
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对于交叉验证调参过程的封装，可以简化编程实现
接下来我们进行数据预处理:
- 读取数据
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
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
- 简单起见我们这里只考虑logistic regression。其他分类器的实现也非常类似。
- 我们这里用两种方式进行特征选择
- recursive feature elimination (RFE): 用交叉验证评估特征在discovery set上的分类效果，可以用sklearn提供的RFECV类实现
- exhaustive enumeration: 枚举所有三个特征的组合，挑选出discovery set上交叉验证效果最好的组合。我们的数据集只有9个特征，所有三个特征的组合总共有84种，这还是可以接受的。当特征数目非常多事，使用这种穷举法几乎是不可能的。
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]
我们首先枚举出所有的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']
- 计算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
我们提供了一个qPCR数据集qPCR_data.csv，第1列为sample id，第2-12列为特征(11个基因的表达量)，第13列为样本标签(负例为健康人:NC,正例为肝癌病人:HCC)。请大家完成以下任务:
- 数据预处理
- 用均值或中位数补全空缺值
- 对数据进行scaling
- 数据可视化
- 用PCA对数据进行可视化
- 数据集划分：预留20%数据用于评估模型泛化能力，剩下的用于模型拟合
- 模型选择和模型拟合
- 任选一种分类器即可
- 调参: 根据所选的分类器对相应的超参数进行搜索
- 特征选择: 用RFE/穷举或其他方式均可，特征数量不限
- 模型评估
- 计算预留数据集上的AUROC
- 绘制ROC曲线
- 请提交代码，必要的文字解释，数据可视化结果及ROC曲线
随机森林是生物信息经常使用的一个分类器。请大家查阅资料，回答以下两个问题:
- 随机森林中树的数量是不是一个需要通过交叉验证调整的超参数?为什么?
- 请问什么是随机森林的out-of-bag (OOB) error?它和bootstrapping有什么关系?