美文网首页临床数据的的机器模型构建
MachineLearning 9. 癌症诊断机器学习之梯度提升

MachineLearning 9. 癌症诊断机器学习之梯度提升

作者: 桓峰基因 | 来源:发表于2022-05-12 12:40 被阅读0次

    前     言

    梯度提升机是一个强大的机器学习技术家族,在广泛的实际应用中显示了相当大的成功。它们可以根据应用程序的特定需求进行高度定制,就像学习不同的损失函数一样。这篇文章提供了一个教程,介绍梯度提升方法的方法论,重点关注建模的机器学习方面。理论信息是补充描述性的例子和插图,涵盖梯度推进模型设计的所有阶段。讨论了处理模型复杂性的注意事项。给出了三个梯度助推应用实例,并进行了综合分析。

    基本原理

    梯度提升法的主要思想是,先建立一个某种形式的初始模型(线性、样条、树或其他),称为基学习器;然后检查残差,在残差的基础上围绕损失函数拟合模型。损失函数测量模型和现实之间的差别,例如,在回归问题中可以用误差的平方。一直继续这个过程,直到满足某个特定的结束条件。这与下面的情形有点相似:一个学生进行模拟考试,100道题中错了30道,然后只研究那30道错题;在下次模考中,30道题中又错了10道,然后只研究那10道题,以此类推。

    实例解析

    1. 软件安装

    if (!require(xgboost)) install.packages("xgboost")

    if (!require(caret)) install.packages("caret")

    library(xgboost)
    library(caret)

    2. 数据读取

    数据来源《机器学习与R语言》书中,具体来自UCI机器学习仓库。地址:http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/ 下载wbdc.data和wbdc.names这两个数据集,数据经过整理,成为面板数据。查看数据结构,其中第一列为id列,无特征意义,需要删除。第二列diagnosis为响应变量(B,M),字符型,一般在R语言中分类任务都要求响应变量为因子类型,因此需要做数据类型转换。剩余的为预测变量,数值类型。查看数据维度,568个样本,32个特征(包括响应特征)。

    BreastCancer <- read.csv("wisc_bc_data.csv", stringsAsFactors = FALSE)
    dim(BreastCancer)
    ## [1] 568 32
    table(BreastCancer$diagnosis)
    ##
    ## B M
    ## 357 211

    3. 数据分布

    比较恶性和良性之间的差距

    library(reshape2)
    bc <- BreastCancer[, -1]
    bc.melt <- melt(bc, id.var = "diagnosis")
    head(bc.melt)
    ## diagnosis variable value
    ## 1 M radius_mean 20.57
    ## 2 M radius_mean 19.69
    ## 3 M radius_mean 11.42
    ## 4 M radius_mean 20.29
    ## 5 M radius_mean 12.45
    ## 6 M radius_mean 18.25
    ggplot(data = bc.melt, aes(x = diagnosis, y = log(value + 1), fill = diagnosis)) +
    geom_boxplot() + theme_bw() + facet_wrap(~variable, ncol = 8)

    比较变量之间的相关性,如下:

    library(tidyverse)
    data <- select(BreastCancer, -1) %>%
    mutate_at("diagnosis", as.factor)
    data$diagnosis = ifelse(data$diagnosis == "B", "B", "M")
    sum(is.na(data))
    ## [1] 0
    corrplot::corrplot(cor(data[, -1]))

    4. 数据分割

    当我们只有一套数据的时候,可以将数据分为训练集和测试集,具体怎么分割可以看公众号的专题:Topic 5. 样本量确定及分割我们将整个数据进行分割,分为训练集和测试集,并保证其正负样本的比例,如下:

    # 数据分割 install.packages('sampling')
    library(sampling)
    set.seed(123)
    # 每层抽取70%的数据
    train_id <- strata(data, "diagnosis", size = rev(round(table(data$diagnosis) * 0.7)))$ID_unit
    # 训练数据
    train_data <- data[train_id, ]
    # 测试数据
    test_data <- data[-train_id, ]

    # 查看训练、测试数据中正负样本比例
    prop.table(table(train_data$diagnosis))
    ##
    ## B M
    ## 0.6281407 0.3718593

    prop.table(table(test_data$diagnosis))
    ##
    ## B M
    ## 0.6294118 0.3705882

    5. 模型构建

    参数配置

    在这一节要使用已经加载的xgboost包。因为这种方法的效果如雷贯耳,所以我们直接用它解决最有挑战性的乳腺癌诊断性问题。正如在提升算法简介中所说,我们要调整一些参数。(1) nrounds:最大迭代次数(最终模型中树的数量)。(2) colsample_bytree:建立树时随机抽取的特征数量,用一个比率表示,默认值为1(使用100%的特征)。(3) min_child_weight:对树进行提升时使用的最小权重,默认为1o(4) eta:学习率,每棵树在最终解中的贡献,默认为0.3。(5) gamma:在树中新增一个叶子分区时所需的最小减损。(6)subsample:子样本数据占整个观测的比例,默认值为1(100%)。(7) max_depth:单个树的最大深度。

    使用expand.grid()函数可以建立实验网格,以运行caret包的训练过程。对于前面列出的参数,如果没有设定具体值,那么即使有默认值,运行函数时也会收到出错信息。下面的参数取值是基于我以前的一些训练迭代而设定的。建议各位亲自实验参数调整过程。使用以下命令建立网格:

    ###参数选择
    grid = expand.grid(
    nrounds = c(75, 100),
    colsample_bytree = 1,
    min_child_weight = 1,
    eta = c(0.01, 0.1, 0.3), #0.3 is default,
    gamma = c(0.5, 0.25),
    subsample = 0.5,
    max_depth = c(2, 3)
    )
    grid
    ## nrounds colsample_bytree min_child_weight eta gamma subsample max_depth
    ## 1 75 1 1 0.01 0.50 0.5 2
    ## 2 100 1 1 0.01 0.50 0.5 2
    ## 3 75 1 1 0.10 0.50 0.5 2
    ## 4 100 1 1 0.10 0.50 0.5 2
    ## 5 75 1 1 0.30 0.50 0.5 2
    ## 6 100 1 1 0.30 0.50 0.5 2
    ## 7 75 1 1 0.01 0.25 0.5 2
    ## 8 100 1 1 0.01 0.25 0.5 2
    ## 9 75 1 1 0.10 0.25 0.5 2
    ## 10 100 1 1 0.10 0.25 0.5 2
    ## 11 75 1 1 0.30 0.25 0.5 2
    ## 12 100 1 1 0.30 0.25 0.5 2
    ## 13 75 1 1 0.01 0.50 0.5 3
    ## 14 100 1 1 0.01 0.50 0.5 3
    ## 15 75 1 1 0.10 0.50 0.5 3
    ## 16 100 1 1 0.10 0.50 0.5 3
    ## 17 75 1 1 0.30 0.50 0.5 3
    ## 18 100 1 1 0.30 0.50 0.5 3
    ## 19 75 1 1 0.01 0.25 0.5 3
    ## 20 100 1 1 0.01 0.25 0.5 3
    ## 21 75 1 1 0.10 0.25 0.5 3
    ## 22 100 1 1 0.10 0.25 0.5 3
    ## 23 75 1 1 0.30 0.25 0.5 3
    ## 24 100 1 1 0.30 0.25 0.5 3

    交叉验证

    使用car包的train()函数之前,我要创建一个名为cntrl的对象,来设定trainControl的参数。这个对象会保存我们要使用的方法,以训练调优参数。我们使用10折交叉验证,代码如下所示:

    cntrl = trainControl(method = "cv", number = 10, verboseIter = FALSE, returnData = FALSE,
    returnResamp = "final")

    根据配置好的参数,进行建模

    set.seed(123)
    train.xgb = train(x = train_data[, -1], y = train_data[, 1], trControl = cntrl, tuneGrid = grid,
    method = "xgbTree")

    train.xgb
    ## eXtreme Gradient Boosting
    ##
    ## No pre-processing
    ## Resampling: Cross-Validated (10 fold)
    ## Summary of sample sizes: 359, 358, 358, 358, 358, 359, ...
    ## Resampling results across tuning parameters:
    ##
    ## eta max_depth gamma nrounds Accuracy Kappa
    ## 0.01 2 0.25 75 0.9448077 0.8815793
    ## 0.01 2 0.25 100 0.9448718 0.8825558
    ## 0.01 2 0.50 75 0.9448077 0.8811347
    ## 0.01 2 0.50 100 0.9448077 0.8814371
    ## 0.01 3 0.25 75 0.9498718 0.8926807
    ## 0.01 3 0.25 100 0.9473077 0.8868425
    ## 0.01 3 0.50 75 0.9448077 0.8815793
    ## 0.01 3 0.50 100 0.9473077 0.8868425
    ## 0.10 2 0.25 75 0.9623718 0.9187054
    ## 0.10 2 0.25 100 0.9598718 0.9133019
    ## 0.10 2 0.50 75 0.9548077 0.9022024
    ## 0.10 2 0.50 100 0.9623718 0.9192980
    ## 0.10 3 0.25 75 0.9648718 0.9247034
    ## 0.10 3 0.25 100 0.9673718 0.9298262
    ## 0.10 3 0.50 75 0.9548077 0.9022182
    ## 0.10 3 0.50 100 0.9598718 0.9128120
    ## 0.30 2 0.25 75 0.9623718 0.9190002
    ## 0.30 2 0.25 100 0.9648718 0.9241266
    ## 0.30 2 0.50 75 0.9623077 0.9184364
    ## 0.30 2 0.50 100 0.9598077 0.9133136
    ## 0.30 3 0.25 75 0.9649359 0.9247785
    ## 0.30 3 0.25 100 0.9699359 0.9357355
    ## 0.30 3 0.50 75 0.9748718 0.9457522
    ## 0.30 3 0.50 100 0.9748718 0.9457522
    ##
    ## Tuning parameter 'colsample_bytree' was held constant at a value of 1
    ##
    ## Tuning parameter 'min_child_weight' was held constant at a value of 1
    ##
    ## Tuning parameter 'subsample' was held constant at a value of 0.5
    ## Accuracy was used to select the optimal model using the largest value.
    ## The final values used for the model were nrounds = 75, max_depth = 3, eta
    ## = 0.3, gamma = 0.5, colsample_bytree = 1, min_child_weight = 1 and subsample
    ## = 0.5.

    最优参数设置

    由此可以得到最优的参数组合来建立模型。模型在训练数据上的正确率是0.975,Kappa值是0.946。下面要做的事情有点复杂,但我认为这才是最佳实践。首先创建一个参数列表,供xgboost包的训练函数xgb.train()使用。然后将数据框转换为一个输入特征矩阵,以及一个带标号的数值型结果列表(其中的值是0和1)。接着,将特征矩阵和标号列表组合成符合要求的输入,即一个xgb.Dmatrix对象。代码如下:

    param <- list(objective = "binary:logistic", booster = "gbtree", eval_metric = "error",
    eta = 0.1, max_depth = 2, subsample = 0.5, colsample_bytree = 1, gamma = 0.5)

    x <- as.matrix(train_data[, -1])
    y <- ifelse(train_data$diagnosis == "B", 0, 1)
    train.mat <- xgb.DMatrix(data = x, label = y)

    优化模型

    set.seed(123)
    xgb.fit <- xgb.train(params = param, data = train.mat, nrounds = 100)
    xgb.fit
    ## ##### xgb.Booster
    ## raw: 87 Kb
    ## call:
    ## xgb.train(params = param, data = train.mat, nrounds = 100)
    ## params (as set within xgb.train):
    ## objective = "binary:logistic", booster = "gbtree", eval_metric = "error", eta = "0.1", max_depth = "2", subsample = "0.5", colsample_bytree = "1", gamma = "0.5", validate_parameters = "TRUE"
    ## xgb.attributes:
    ## niter
    ## callbacks:
    ## cb.print.evaluation(period = print_every_n)
    ## # of features: 30
    ## niter: 100
    ## nfeatures : 30

    可视化变量

    从图上我们可以出前几个重要的变量占比几乎达到0.3,这个阈值是在之前grid网格中设置:eta = c(0.01, 0.1, 0.3), #0.3 is default。

    impMatrix <- xgb.importance(feature_names = dimnames(x)[[2]], model = xgb.fit)
    impMatrix
    ## Feature Gain Cover Frequency
    ## 1: concave_points_worst 0.2967527144 0.147574780 0.074766355
    ## 2: concave_points_mean 0.2341084307 0.144084708 0.079439252
    ## 3: radius_worst 0.0776151887 0.075881804 0.056074766
    ## 4: perimeter_worst 0.0742989918 0.080783308 0.065420561
    ## 5: area_se 0.0707442192 0.075966262 0.065420561
    ## 6: texture_worst 0.0431688078 0.084473059 0.107476636
    ## 7: area_worst 0.0349753102 0.054806620 0.056074766
    ## 8: concavity_worst 0.0274044659 0.045700578 0.060747664
    ## 9: texture_mean 0.0242682068 0.042618830 0.056074766
    ## 10: smoothness_worst 0.0180421347 0.030139745 0.046728972
    ## 11: area_mean 0.0137108659 0.022756702 0.018691589
    ## 12: compactne_mean 0.0105402518 0.012610045 0.028037383
    ## 13: symmetry_se 0.0104209161 0.020799297 0.042056075
    ## 14: smoothness_mean 0.0090815731 0.020561037 0.032710280
    ## 15: radius_se 0.0084606195 0.029712377 0.018691589
    ## 16: compactne_se 0.0080685708 0.014326525 0.032710280
    ## 17: symmetry_worst 0.0076659203 0.027927952 0.037383178
    ## 18: fractal_dimension_mean 0.0057942369 0.008967406 0.014018692
    ## 19: smoothness_se 0.0050382309 0.018076694 0.023364486
    ## 20: fractal_dimension_se 0.0043507637 0.007250846 0.018691589
    ## 21: fractal_dimension_worst 0.0029493268 0.006075975 0.009345794
    ## 22: symmetry_mean 0.0026459958 0.006109480 0.014018692
    ## 23: concavity_mean 0.0023260374 0.003932079 0.009345794
    ## 24: perimeter_mean 0.0022726213 0.003359381 0.004672897
    ## 25: compactne_worst 0.0019021532 0.003878213 0.009345794
    ## 26: concave_points_se 0.0011308883 0.002550426 0.004672897
    ## 27: radius_mean 0.0010039937 0.002045465 0.004672897
    ## 28: perimeter_se 0.0008485220 0.002095077 0.004672897
    ## 29: texture_se 0.0004100425 0.004935332 0.004672897
    ## Feature Gain Cover Frequency
    xgb.plot.importance(impMatrix, main = "Gain by Feature")

    模型性能评估

    library(InformationValue)
    pred <- predict(xgb.fit, x)
    optimalCutoff(y, pred)
    ## [1] 0.2692855

    测试集验证

    testMat <- as.matrix(test_data[, -1])
    xgb.test <- predict(xgb.fit, testMat)
    y.test <- ifelse(test_data$diagnosis == "B", 0, 1)
    optimalCutoff(y.test, xgb.test)
    ## [1] 0.1793153

    混淆矩阵计算假阳和假阴个数,如下:

    confusionMatrix(y.test, xgb.test, threshold = 0.29)
    ## 0 1
    ## 0 103 0
    ## 1 4 63
    1 - misClassError(y.test, xgb.test, threshold = 0.29)
    ## [1] 0.9765

    绘制ROC曲线

    最后就是模型的准确性评估,这里我们使用的是InformationValue软件包中的plotROC,绘制ROC曲线,如下:

    ###
    library(InformationValue)
    plotROC(y.test, xgb.test)

    结果解读

    我们再回忆一下对乳腺癌数据我们都用过哪些机器学习方法。

    基于乳腺癌的数据我们已经做过四种类型的机器学习算法,如下:

    K-邻近算法(KNN)准确率为0.9471

    支持向量机(SVM)准确率为 0.9765

    分类随机森林RT准确率为0.969

    这期使用梯度提升准确率为0.997,准确率嗖的一下就上去了,这模型应该是让老板十万分满意。

    其实这说明做模型时,不应单一只是用一种算法,需要多种算法比较,找到最优的选择!

    还有就是注意绘制ROC的方法,在做分类随机森林时我们使用的是ROSE软件包,而在做回归随机森林时我们使用的是ROCR,需要注意使用时方法的选择。在这里我们使用InformationValue软件包里面的plotROC

    References:

    1. Frontiers in Neurorobotics, Gradient boosting machines,a tutorial,Natekin A., Knoll A.(2013)

    本文使用 文章同步助手 同步

    相关文章

      网友评论

        本文标题:MachineLearning 9. 癌症诊断机器学习之梯度提升

        本文链接:https://www.haomeiwen.com/subject/gfbsurtx.html