今天在读一篇文章的时候，其中对细胞聚类的判断，为了exclude low confidence assignments by clustering we used a classifier for final assignment of cells to clusters，作者采用了一个叫Gradient boosting的方法（R包xgboost，这是一个机器学习的算法，大家熟悉的还有随机森林），applied to build a cell to cluster classifier for each of the top-level seven cluster types and subsequently to each of the 88 low-level clusters（构建聚类的分类器），对于这个，我们需要深入看看了。

关于XGboost

1. 什么是XGBoost

XGBoost是陈天奇等人开发的一个开源机器学习项目，高效地实现了GBDT算法并进行了算法和工程上的许多改进，被广泛应用在Kaggle竞赛及其他许多机器学习竞赛中并取得了不错的成绩。

说到XGBoost，不得不提GBDT(Gradient Boosting Decision Tree)。因为XGBoost本质上还是一个GBDT，但是力争把速度和效率发挥到极致，所以叫X (Extreme) GBoosted。包括前面说过，两者都是boosting方法。

关于GBDT，这里不再提，可以查看我前一篇的介绍，点此跳转。

1.1 XGBoost树的定义

先来举个例子，我们要预测一家人对电子游戏的喜好程度，考虑到年轻和年老相比，年轻更可能喜欢电子游戏，以及男性和女性相比，男性更喜欢电子游戏，故先根据年龄大小区分小孩和大人，然后再通过性别区分开是男是女，逐一给各人在电子游戏喜好程度上打分，如下图所示。

image

就这样，训练出了2棵树tree1和tree2，类似之前gbdt的原理，两棵树的结论累加起来便是最终的结论，所以小孩的预测分数就是两棵树中小孩所落到的结点的分数相加：2 + 0.9 = 2.9。爷爷的预测分数同理：-1 + （-0.9）= -1.9。具体如下图所示：

image

事实上，如果不考虑工程实现、解决问题上的一些差异，XGBoost与GBDT比较大的不同就是目标函数的定义。XGBoost的目标函数如下图所示：

image

其中：

• 红色箭头所指向的L 即为损失函数（比如平方损失函数：(,i)=(−))
• 红色方框所框起来的是正则项（包括L1正则、L2正则）
• 红色圆圈所圈起来的为常数项
• 对于f(x)，XGBoost利用泰勒展开三项，做一个近似。f(x)表示的是其中一颗回归树。

看到这里可能有些读者会头晕了，这么多公式，我在这里只做一个简要式的讲解，具体的算法细节和公式求解请查看这篇博文，讲得很仔细：通俗理解kaggle比赛大杀器xgboost

XGBoost的核心算法思想不难，基本就是：

1. 不断地添加树，不断地进行特征分裂来生长一棵树，每次添加一个树，其实是学习一个新函数f(x)，去拟合上次预测的残差。
2. 当我们训练完成得到k棵树，我们要预测一个样本的分数，其实就是根据这个样本的特征，在每棵树中会落到对应的一个叶子节点，每个叶子节点就对应一个分数
3. 最后只需要将每棵树对应的分数加起来就是该样本的预测值。

显然，我们的目标是要使得树群的预测值^′尽量接近真实值。，而且有尽量大的泛化能力。类似之前GBDT的套路，XGBoost也是需要将多棵树的得分累加得到最终的预测得分（每一次迭代，都在现有树的基础上，增加一棵树去拟合前面树的预测结果与真实值之间的残差）。

image

那接下来，我们如何选择每一轮加入什么 f 呢？答案是非常直接的，选取一个 f 来使得我们的目标函数尽量最大地降低。这里 f 可以使用泰勒展开公式近似。

image

实质是把样本分配到叶子结点会对应一个obj，优化过程就是obj优化。也就是分裂节点到叶子不同的组合，不同的组合对应不同obj，所有的优化围绕这个思想展开。到目前为止我们讨论了目标函数中的第一个部分：训练误差。接下来我们讨论目标函数的第二个部分：正则项，即如何定义树的复杂度。

1.2 正则项：树的复杂度

XGBoost对树的复杂度包含了两个部分：

• 一个是树里面叶子节点的个数T
• 一个是树上叶子节点的得分w的L2模平方（对w进行L2正则化，相当于针对每个叶结点的得分增加L2平滑，目的是为了避免过拟合）

image

我们再来看一下XGBoost的目标函数（损失函数揭示训练误差 + 正则化定义复杂度）：

图片.png

正则化公式也就是目标函数的后半部分，对于上式而言，^'是整个累加模型的输出，正则化项∑kΩ(ft)是则表示树的复杂度的函数，值越小复杂度越低，泛化能力越强。

1.3 树该怎么长

很有意思的一个事是，我们从头到尾了解了xgboost如何优化、如何计算，但树到底长啥样，我们却一直没看到。很显然，一棵树的生成是由一个节点一分为二，然后不断分裂最终形成为整棵树。那么树怎么分裂的就成为了接下来我们要探讨的关键。对于一个叶子节点如何进行分裂，XGBoost作者在其原始论文中给出了一种分裂节点的方法：枚举所有不同树结构的贪心法

不断地枚举不同树的结构，然后利用打分函数来寻找出一个最优结构的树，接着加入到模型中，不断重复这样的操作。这个寻找的过程使用的就是贪心算法。选择一个feature分裂，计算loss function最小值，然后再选一个feature分裂，又得到一个loss function最小值，你枚举完，找一个效果最好的，把树给分裂，就得到了小树苗。

总而言之，XGBoost使用了和CART回归树一样的想法，利用贪婪算法，遍历所有特征的所有特征划分点，不同的是使用的目标函数不一样。具体做法就是分裂后的目标函数值比单子叶子节点的目标函数的增益，同时为了限制树生长过深，还加了个阈值，只有当增益大于该阈值才进行分裂。从而继续分裂，形成一棵树，再形成一棵树，每次在上一次的预测基础上取最优进一步分裂/建树。

1.4 如何停止树的循环生成

凡是这种循环迭代的方式必定有停止条件，什么时候停止呢？简言之，设置树的最大深度、当样本权重和小于设定阈值时停止生长以防止过拟合。具体而言，则

1. 当引入的分裂带来的增益小于设定阀值的时候，我们可以忽略掉这个分裂，所以并不是每一次分裂loss function整体都会增加的，有点预剪枝的意思，阈值参数为（即正则项里叶子节点数T的系数）；
2. 当树达到最大深度时则停止建立决策树，设置一个超参数max_depth，避免树太深导致学习局部样本，从而过拟合；
3. 样本权重和小于设定阈值时则停止建树。什么意思呢，即涉及到一个超参数-最小的样本权重和min_child_weight，和GBM的 min_child_leaf 参数类似，但不完全一样。大意就是一个叶子节点样本太少了，也终止同样是防止过拟合；

2. XGBoost与GBDT有什么不同

除了算法上与传统的GBDT有一些不同外，XGBoost还在工程实现上做了大量的优化。总的来说，两者之间的区别和联系可以总结成以下几个方面。

1. GBDT是机器学习算法，XGBoost是该算法的工程实现。
2. 在使用CART作为基分类器时，XGBoost显式地加入了正则项来控制模 型的复杂度，有利于防止过拟合，从而提高模型的泛化能力。
3. GBDT在模型训练时只使用了代价函数的一阶导数信息，XGBoost对代 价函数进行二阶泰勒展开，可以同时使用一阶和二阶导数。
4. 传统的GBDT采用CART作为基分类器，XGBoost支持多种类型的基分类 器，比如线性分类器。
5. 传统的GBDT在每轮迭代时使用全部的数据，XGBoost则采用了与随机 森林相似的策略，支持对数据进行采样。
6. 传统的GBDT没有设计对缺失值进行处理，XGBoost能够自动学习出缺 失值的处理策略。

3. 为什么XGBoost要用泰勒展开，优势在哪里？

XGBoost使用了一阶和二阶偏导, 二阶导数有利于梯度下降的更快更准. 使用泰勒展开取得函数做自变量的二阶导数形式, 可以在不选定损失函数具体形式的情况下, 仅仅依靠输入数据的值就可以进行叶子分裂优化计算, 本质上也就把损失函数的选取和模型算法优化/参数选择分开了. 这种去耦合增加了XGBoost的适用性, 使得它按需选取损失函数, 可以用于分类, 也可以用于回归。

那么，xgboost在单细胞数据的运用是什么？？

For each of the seven top-level cell types, a separate classifier was trained to predict each cell type separately (one-versusall), in a 5-fold cross-validation scheme（对细胞类型进行训练）. Next, using the probability scores of the held-out test-set we identified an optimal cutoff for each class based on an ROC analysis comparing the true positive rate (TPR = true positives divided by all positive predictions) to the false positive rate (FPR = true negative divided by all negatives) and selecting the point at which the ROC curve intersects with the diagonal（常见的检验方法）. Cells that were ambiguously assigned in this way to more than one cluster were removed as potential doublets.（多细胞是这么识别的，有意思）。

Next, a similar classification training scheme was applied separately to cells from each top-level cell type (epithelial, stromal, mast, B, plasma, myeloid, and T cells). We used 5-fold cross-validation and ROC analysis to select thresholds. In cases where a cell was assigned to more than one subtype, we used the assignment with the higher predictive score（阈值的调整这里还是很重要的）. Cells that could not be assigned confidently by any classifier were excluded from further analysis.（无法进行分类的细胞进一步分析）。

这个方法还用来Classifying malignant cells by gradient boosting（识别恶性细胞）

Adjacent normal tissue, which was sampled distantly from the tumor (e.g., ~10cm apart), is expected to be tumor-cell free. We used gradient boosting to train a classifier predicting malignant from non-malignant epithelial cells based on the source channel type (tumor versus adjacent normal), in a 5-fold cross validation scheme（根据样本的来源进行训练）. We separately trained two classifiers, one predicting isTumor and another predicting isNormal, and used the geometric mean of the resulting probabilities as the final statistic. In subsequent analyses, we considered epithelial cells from tumor channels with a predicted score greater than 0.75 to be malignant, and cells from normal channels with a predicted score < 0.25 to be normal epithelial cells（细胞的划分）. Overall, by this measure ~95% of tumor channel epithelial cells were predicted to be malignant, and 98% of normal channel epithelial cells were predicted to be non-malignant cells. The classifier predictions were highly concordant with those made by inferred copy number alterations with only ~11% of likely malignant cells showing no substantial copy number differences from normal (8% for MMRp, and 15% for MMRd), and 2% of likely normal cells showing copy number differences (data not shown). Copy number alterations were only determined for epithelial cells.看来分类器的作用了，就是为了识别恶性细胞。

R语言实现xgboost

一、XGBoost包

我们可以直接调用R语言的xgboost包中的xgboost()函数实现该算法，xgboost包已经上传到CRAN上，可直接通过install.packages("xgboost")进行安装。下面首先对几个主要参数进行介绍

一般参数

• silent : 默认值是0。您需要指定0连续打印消息,静默模式1。

• booster : 默认值是gbtree。你需要指定要使用的上升模型:gbtree(树)或gblinear(线性函数)。

• num_pbuffer : 这是由xgboost自动设置,不需要由用户设定。阅读xgboost文档的更多细节。

• num_feature : 这是由xgboost自动设置,不需要由用户设定。

辅助参数

具体参数树状图：

• eta：默认值设置为0.3。您需要指定用于更新步长收缩来防止过度拟合。每个提升步骤后,我们可以直接获得新特性的权重。实际上 eta 收缩特征权重的提高过程更为保守。范围是0到1。低η值意味着模型过度拟合更健壮。

• gamma:默认值设置为0。您需要指定最小损失减少应进一步划分树的叶节点。
  更大,更保守的算法。范围是0到∞。γ越大算法越保守。

• max_depth:默认值设置为6。您需要指定一个树的最大深度。参数范围是1到∞。

• min_child_weight:默认值设置为1。您需要在子树中指定最小的（海塞）实例权重的和，然后这个构建过程将放弃进一步的分割。在线性回归模式中，在每个节点最少所需实例数量将简单的同时部署。更大,更保守的算法。参数范围是0到∞。

• max_delta_step：默认值设置为0。max_delta_step 允许我们估计每棵树的权重。如果该值设置为0,这意味着没有约束。
  如果它被设置为一个正值,它可以帮助更新步骤更为保守。通常不需要此参数,但是在逻辑回归中当分类是极为不均衡时需要用到。将其设置为1 - 10的价值可能有助于控制更新。参数范围是0到∞。

• subsample： 默认值设置为1。您需要指定训练实例的子样品比。
  设置为0.5意味着XGBoost随机收集一半的数据实例来生成树来防止过度拟合。参数范围是0到1。

• colsample_bytree : 默认值设置为1。在构建每棵树时,您需要指定列的子样品比。范围是0到1。

线性上升具体参数

• lambda and alpha : 这些都是正则化项权重。λ默认值假设是1和α= 0。

• lambda_bias : L2正则化项在偏差上的默认值为0。

任务参数

• base_score : 默认值设置为0.5。您需要指定初始预测分数作为全局偏差。

• objective : 默认值设置为reg:linear。您需要指定你想要的类型的学习者,包括线性回归、逻辑回归、泊松回归等。

• eval_metric : 您需要指定验证数据的评估指标,一个默认的指标分配根据客观(rmse回归,错误分类,意味着平均精度等级

• seed : 随机数种子，确保重现数据相同的输出。

数据预处理

xgboost包中的xgboost函数对于数据格式有所要求，在使用xgboost函数之前，我们需要先对数据进行预处理。

library(Matrix)
####训练集的数据预处理
# 将自变量转化为矩阵
traindata1 <- data.matrix(train_data[,c(1:11)]) 
# 利用Matrix函数，将sparse参数设置为TRUE，转化为稀疏矩阵
traindata2 <- Matrix(traindata1,sparse=T) 
traindata3 <- train_data[,13]
# 将自变量和因变量拼接为list
traindata4 <- list(data=traindata2,label=traindata3) 
# 构造模型需要的xgb.DMatrix对象，处理对象为稀疏矩阵
dtrain <- xgb.DMatrix(data = traindata4$data, label = traindata4$label) 

####测试集的数据预处理
# 将自变量转化为矩阵
testset1 <- data.matrix(test_data[,c(1:11)]) 
# 利用Matrix函数，将sparse参数设置为TRUE，转化为稀疏矩阵
testset2 <- Matrix(testset1,sparse=T) 
# 将因变量转化为numeric
testset3 <- test_data[,13]
# 将自变量和因变量拼接为list
testset4 <- list(data=testset2,label=testset3) 
# 构造模型需要的xgb.DMatrix对象，处理对象为稀疏矩阵
dtest <- xgb.DMatrix(data = testset4$data, label = testset4$label) 
xgb <- xgboost(data = dtrain,max_depth=6, eta=0.5,  objective='binary:logistic', nround=25)

ROC曲线和AUC值

#在测试集上预测
pre_xgb = round(predict(xgb,newdata = dtest))
#输出混淆矩阵
table(test_data$等级,pre_xgb,dnn=c("真实值","预测值"))
xgboost_roc <- roc(test_data$等级,as.numeric(pre_xgb))
#绘制ROC曲线和AUC值
plot(xgboost_roc, print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2),grid.col=c("green", "red"), max.auc.polygon=TRUE,auc.polygon.col="skyblue", print.thres=TRUE,main='xgboost模型ROC曲线')

图片.png

关于分类器，是机器学习很深的一部分，大家如果对数学感兴趣，好好学习。

生活很好，有你更好