> For the complete documentation index, see [llms.txt](https://book.ncrnalab.org/teaching/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://book.ncrnalab.org/teaching/part-iv.-machine-learning/3.machine-learning-with-python.md). # 3.Machine Learning with Python [scikit-learn](https://scikit-learn.org/)等python package对很多机器学习模型提供了易于使用的接口。但是，如果想用机器学习方法解决实际问题，我们需要考虑的步骤不只是拟合模型，还包括数据预处理，特征选择，调参，模型评估等等。我们只有对数据有比较深入理解，才能选择合适的方法。本节中我们提供一个用基于python的机器学习工具解决实际的分类问题的例子。我们希望读者读者结合本章内容，能通过实践掌握相关工具的使用。 ## 1) 数据说明 我们这里用到的还是 [BreastCancer.csv](https://cloud.tsinghua.edu.cn/d/e3215eb0090f4b3dac1e/?p=%2FFiles%2FPART_IV\&mode=list) 数据集。 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 | ## 2) Load python packages ```python 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) 数据预处理 接下来我们进行数据预处理: * 读取数据 ```python 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](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html) 方法来随机的将20%的数据预留为一个用于最终评估模型泛化能力的数据集，剩下80%数据为可以见到label的Discovery set。 ```python 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。 ```python 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只会保留一部分特征: ```python 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个特征的组合 ```python 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类，把能使用到的特征的标号也作为一个超参数: ```python 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评估每个特征组合的性能: ```python 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 ```python y_pred_proba = clf.predict_proba(X_validation)[:,1] fpr, tpr,_ = roc_curve(y_validation,y_pred_proba) AUROC = auc(fpr, tpr) ``` * 绘制ROC曲线 ```python 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](/files/EagauizNmAA1IYCSC22o) ## 7) Homework ### 7.1) 我们提供了一个qPCR数据集[qPCR\_data.csv](https://cloud.tsinghua.edu.cn/d/ad22768345664924b202/?p=%2FFiles%2FPART_IV\&mode=list)，第1列为sample id，第2-12列为特征(11个基因的表达量)，第13列为样本标签(负例为健康人:NC,正例为肝癌病人:HCC)。请大家完成以下任务: * 数据预处理 * 用均值或中位数补全空缺值 * 对数据进行scaling * 数据可视化 * 用PCA对数据进行可视化 * 数据集划分：预留20%数据用于评估模型泛化能力，剩下的用于模型拟合 * 模型选择和模型拟合 * 任选一种分类器即可 * 调参: 根据所选的分类器对相应的超参数进行搜索 * 特征选择: 用RFE/穷举或其他方式均可，特征数量不限 * 模型评估 * 计算预留数据集上的AUROC * 绘制ROC曲线 * 请提交代码，必要的文字解释，数据可视化结果及ROC曲线 ### 7.2) 随机森林是生物信息经常使用的一个分类器。请大家查阅资料，回答以下两个问题: * 随机森林中树的数量是不是一个 超参数?为什么 * 请问什么是随机森林的out-of-bag (OOB) error?它和bootstrapping有什么关系?