3.Machine Learning with Python

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

1) 数据说明

我们这里用到的还是 BreastCancer.csv 数据集。

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

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

• 数据可视化

  • 用PCA对数据进行可视化

• 数据集划分：预留20%数据用于评估模型泛化能力，剩下的用于模型拟合

• 模型选择和模型拟合

  • 任选一种分类器即可

  • 调参: 根据所选的分类器对相应的超参数进行搜索

  • 特征选择: 用RFE/穷举或其他方式均可，特征数量不限

• 模型评估

  • 计算预留数据集上的AUROC

  • 绘制ROC曲线

• 请提交代码，必要的文字解释，数据可视化结果及ROC曲线

7.2)

随机森林是生物信息经常使用的一个分类器。请大家查阅资料，回答以下两个问题:

• 随机森林中树的数量是不是一个 超参数?为什么

• 请问什么是随机森林的out-of-bag (OOB) error?它和bootstrapping有什么关系?