# 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有什么关系? --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://book.ncrnalab.org/teaching/part-iv.-machine-learning/3.machine-learning-with-python.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.