1 基于树的模型

https://bonsai.tidymodels.org/articles/bonsai.html

Code
library(tidymodels)
#> ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
#> ✔ broom        1.0.6     ✔ rsample      1.2.1
#> ✔ dials        1.2.1     ✔ tune         1.2.1
#> ✔ infer        1.0.7     ✔ workflows    1.1.4
#> ✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
#> ✔ parsnip      1.2.1     ✔ yardstick    1.3.1
#> ✔ recipes      1.1.0
library(vip)

1.1 决策树 (分类)

Code
df <- read_csv("data/breast-cancer-wisconsin.data",col_names = F,na = c("","NA","?"))
#> Rows: 699 Columns: 11
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> dbl (11): X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

colnames(df) <- c("id","肿块厚度","细胞大小均匀性","细胞形状均匀性","边际附着力","单个上皮细胞大小",
               "裸核","bland_chromatin","正常核","有丝分裂","class")

df <- df |>
    select(-1) |>
    mutate(class = factor(class, levels = c(2, 4), labels = c("良性", "恶性"))) |>
    drop_na()
glimpse(df)
#> Rows: 683
#> Columns: 10
#> $ 肿块厚度         <dbl> 5, 5, 3, 6, 4, 8, 1, 2, 2, 4, 1, 2, 5, 1, 8, 7, 4, 4,…
#> $ 细胞大小均匀性   <dbl> 1, 4, 1, 8, 1, 10, 1, 1, 1, 2, 1, 1, 3, 1, 7, 4, 1, 1…
#> $ 细胞形状均匀性   <dbl> 1, 4, 1, 8, 1, 10, 1, 2, 1, 1, 1, 1, 3, 1, 5, 6, 1, 1…
#> $ 边际附着力       <dbl> 1, 5, 1, 1, 3, 8, 1, 1, 1, 1, 1, 1, 3, 1, 10, 4, 1, 1…
#> $ 单个上皮细胞大小 <dbl> 2, 7, 2, 3, 2, 7, 2, 2, 2, 2, 1, 2, 2, 2, 7, 6, 2, 2,…
#> $ 裸核             <dbl> 1, 10, 2, 4, 1, 10, 10, 1, 1, 1, 1, 1, 3, 3, 9, 1, 1,…
#> $ bland_chromatin  <dbl> 3, 3, 3, 3, 3, 9, 3, 3, 1, 2, 3, 2, 4, 3, 5, 4, 2, 3,…
#> $ 正常核           <dbl> 1, 2, 1, 7, 1, 7, 1, 1, 1, 1, 1, 1, 4, 1, 5, 3, 1, 1,…
#> $ 有丝分裂         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 4, 1, 1, 1,…
#> $ class            <fct> 良性, 良性, 良性, 良性, 良性, 恶性, 良性, 良性, 良性,…
Code
# 2=良性 4=恶性
table(df$class)
#> 
#> 良性 恶性 
#>  444  239
# 拆分训练集和测试集         ####
set.seed(100)
split <- initial_split(df, prop = 0.70, strata = class)

split
#> <Training/Testing/Total>
#> <477/206/683>
train <- training(split)
test  <-  testing(split)

table(train$class)
#> 
#> 良性 恶性 
#>  310  167
table(test$class)
#> 
#> 良性 恶性 
#>  134   72
Code
class_tree_spec <- decision_tree() %>%
    set_engine("rpart") %>%
    set_mode("classification") 

cdtree <- class_tree_spec |> fit(class ~ . ,data = train)
cdtree
#> parsnip model object
#> 
#> n= 477 
#> 
#> node), split, n, loss, yval, (yprob)
#>       * denotes terminal node
#> 
#>  1) root 477 167 良性 (0.64989518 0.35010482)  
#>    2) 细胞大小均匀性< 2.5 293  11 良性 (0.96245734 0.03754266)  
#>      4) 裸核< 5.5 285   4 良性 (0.98596491 0.01403509) *
#>      5) 裸核>=5.5 8   1 恶性 (0.12500000 0.87500000) *
#>    3) 细胞大小均匀性>=2.5 184  28 恶性 (0.15217391 0.84782609)  
#>      6) 细胞形状均匀性< 2.5 16   3 良性 (0.81250000 0.18750000) *
#>      7) 细胞形状均匀性>=2.5 168  15 恶性 (0.08928571 0.91071429)  
#>       14) 细胞大小均匀性< 4.5 46  12 恶性 (0.26086957 0.73913043)  
#>         28) 裸核< 2.5 10   3 良性 (0.70000000 0.30000000) *
#>         29) 裸核>=2.5 36   5 恶性 (0.13888889 0.86111111) *
#>       15) 细胞大小均匀性>=4.5 122   3 恶性 (0.02459016 0.97540984) *

1.1.1 模型可视化

Code
rpart::plotcp(cdtree$fit)

Code
cdtree %>%
    extract_fit_engine() %>%
    rpart.plot::rpart.plot(roundint = F)

1.1.2 模型性能评估

Code
augment(cdtree, new_data = test) %>%
    accuracy(truth = class, estimate = .pred_class)
#> # A tibble: 1 × 3
#>   .metric  .estimator .estimate
#>   <chr>    <chr>          <dbl>
#> 1 accuracy binary         0.971

augment(cdtree, new_data = test) %>%
    conf_mat(truth = class, estimate = .pred_class)
#>           Truth
#> Prediction 良性 恶性
#>       良性  132    4
#>       恶性    2   68

1.2 随机森林 (分类)

随机森林是装袋法的一种扩展,它不仅对数据进行随机抽样,还对特征进行随机抽样,以此增加模型的多样性和泛化能力。随机森林由大量决策树组成,每棵树都是在一个随机抽取的样本和特征子集上训练的。

Code
cf_spec_class <-
    rand_forest(#mtry = .cols(), 
        trees = 500 ,min_n = 1) %>%
    set_engine('randomForest', importance = TRUE) %>%
    set_mode('classification')

class_rf_fit <- cf_spec_class |> 
    fit(class ~ . , data = train)
class_rf_fit
#> parsnip model object
#> 
#> 
#> Call:
#>  randomForest(x = maybe_data_frame(x), y = y, ntree = ~500, nodesize = min_rows(~1,      x), importance = ~TRUE) 
#>                Type of random forest: classification
#>                      Number of trees: 500
#> No. of variables tried at each split: 3
#> 
#>         OOB estimate of  error rate: 3.77%
#> Confusion matrix:
#>      良性 恶性 class.error
#> 良性  300   10  0.03225806
#> 恶性    8  159  0.04790419

1.2.1 特征重要性:基于Gini系数的减少

OOB ,out of bag 袋外预测误差

这种方法主要用于分类任务。

Mean Decrease in Gini (MDG) 这种方法通过衡量某个特征对分类纯度的贡献来计算其重要性。具体步骤如下:

  1. 训练模型:使用所有特征训练随机森林模型。

  2. 计算Gini系数:在决策树中,每次节点分裂都会计算Gini系数减少量。Gini系数用于衡量数据集的纯度,越低表示越纯。

  3. 累加Gini减少量:在每棵树中,计算每个特征在分裂过程中带来的Gini减少量,并将这些减少量累加起来。

  4. 计算平均值:对所有树的累加值取平均值,作为该特征的重要性得分。

Code
class_rf_fit %>% vip::vi()
#> # A tibble: 9 × 2
#>   Variable         Importance
#>   <chr>                 <dbl>
#> 1 裸核                  25.6 
#> 2 肿块厚度              21.0 
#> 3 bland_chromatin       19.2 
#> 4 细胞大小均匀性        18.4 
#> 5 细胞形状均匀性        17.9 
#> 6 边际附着力            14.0 
#> 7 正常核                13.9 
#> 8 单个上皮细胞大小      10.6 
#> 9 有丝分裂               5.88
class_rf_fit %>% vip::vip()

Code
class_rf_fit$fit$importance
#>                         良性        恶性 MeanDecreaseAccuracy MeanDecreaseGini
#> 肿块厚度         0.045985557 0.041538273          0.044359259        12.465397
#> 细胞大小均匀性   0.049172066 0.076704399          0.058509661        47.927558
#> 细胞形状均匀性   0.009876762 0.096802558          0.040209397        41.331213
#> 边际附着力       0.016551632 0.035342949          0.023040888         6.945234
#> 单个上皮细胞大小 0.012280597 0.010654858          0.011723405        16.171318
#> 裸核             0.061713383 0.066225145          0.063075655        40.588065
#> bland_chromatin  0.017570370 0.058096398          0.031839179        31.183405
#> 正常核           0.027576658 0.021034748          0.025191390        17.830550
#> 有丝分裂         0.003095108 0.001130822          0.002400936         1.946241
Code

augment(class_rf_fit, new_data = test) %>%
    accuracy(truth = class, estimate = .pred_class)
#> # A tibble: 1 × 3
#>   .metric  .estimator .estimate
#>   <chr>    <chr>          <dbl>
#> 1 accuracy binary         0.990

augment(class_rf_fit, new_data = test) %>%
    conf_mat(truth = class, estimate = .pred_class)
#>           Truth
#> Prediction 良性 恶性
#>       良性  132    0
#>       恶性    2   72

1.2.2 基于表达数据的应用

Code
df <- dendextend::khan

df$train.classes
#>  [1] EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS   
#> [11] EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS   
#> [21] EWS    EWS    EWS    BL-NHL BL-NHL BL-NHL BL-NHL BL-NHL BL-NHL BL-NHL
#> [31] BL-NHL NB     NB     NB     NB     NB     NB     NB     NB     NB    
#> [41] NB     NB     NB     RMS    RMS    RMS    RMS    RMS    RMS    RMS   
#> [51] RMS    RMS    RMS    RMS    RMS    RMS    RMS    RMS    RMS    RMS   
#> [61] RMS    RMS    RMS    RMS   
#> Levels: EWS BL-NHL NB RMS
train <- t(df$train) |> bind_cols(tibble(class=df$train.classes)) |> 
    relocate(class, .before = 1) |> 
    mutate(
        class=factor(class,levels = c("EWS", "BL-NHL", "NB","RMS"))
    )
str(train$class)
#>  Factor w/ 4 levels "EWS","BL-NHL",..: 1 1 1 1 1 1 1 1 1 1 ...
table(train$class)
#> 
#>    EWS BL-NHL     NB    RMS 
#>     23      8     12     21
Code
df$test.classes
#>  [1] Normal Normal Normal NB     RMS    Normal Normal NB     EWS    RMS   
#> [11] BL-NHL EWS    RMS    EWS    EWS    EWS    RMS    BL-NHL RMS    NB    
#> [21] NB     NB     NB     BL-NHL EWS   
#> Levels: EWS BL-NHL NB RMS Normal

test <- t(df$test) |> bind_cols(tibble(class=df$test.classes)) |> 
    relocate(class, .before = 1) |> 
    mutate(
        class=factor(class,levels = c("EWS", "BL-NHL", "NB","RMS","Normal"))
    )
str(test$class)
#>  Factor w/ 5 levels "EWS","BL-NHL",..: 5 5 5 3 4 5 5 3 1 4 ...
table(test$class)
#> 
#>    EWS BL-NHL     NB    RMS Normal 
#>      6      3      6      5      5
Code
#
dt <- class_tree_spec |> fit(class ~ . ,data = train)
rpart::plotcp(dt$fit)

Code

dt%>%
    extract_fit_engine() %>%
    rpart.plot::rpart.plot(roundint = F)

Code
#
rf_fit_eg <- cf_spec_class |> 
    fit(class ~ . , data = train)
rf_fit_eg
#> parsnip model object
#> 
#> 
#> Call:
#>  randomForest(x = maybe_data_frame(x), y = y, ntree = ~500, nodesize = min_rows(~1,      x), importance = ~TRUE) 
#>                Type of random forest: classification
#>                      Number of trees: 500
#> No. of variables tried at each split: 17
#> 
#>         OOB estimate of  error rate: 0%
#> Confusion matrix:
#>        EWS BL-NHL NB RMS class.error
#> EWS     23      0  0   0           0
#> BL-NHL   0      8  0   0           0
#> NB       0      0 12   0           0
#> RMS      0      0  0  21           0

1.3 随机森林 (回归)

Code
library(ggplot2)
library(dplyr)
library(tidymodels)
data <- haven::read_sav("data/抑郁随机森林模型变量重要性.sav") 


# glimpse(data)
# 
# attributes(data$Gender)
# attributes(data$ZD)
data <- data %>% 
    mutate(
        Gender = factor(Gender),
        XK = factor(XK),
        DQ = factor(DQ),
        SYDLX = factor(SYDLX),
        
    )

set.seed(123)
split <- initial_split(data, prop = 0.7)
train_data <- training(split)
test_data <- testing(split)
Code
# 构建随机森林模型
rf_sepc_reg <- rand_forest(mtry = 5, trees = 500) %>%
    set_engine("randomForest" , importance = TRUE
               ) %>%
    set_mode("regression")

# 创建配方
rf_recipe <- recipe(ZD ~ ., data = train_data) 

# 工作流程
rf_workflow <- workflow() %>%
    add_recipe(rf_recipe) %>%
    add_model(rf_sepc_reg)

# 训练模型
reg_rf_fit <- rf_workflow %>%
    fit(data = train_data)
reg_rf_fit %>% extract_fit_parsnip()
#> parsnip model object
#> 
#> 
#> Call:
#>  randomForest(x = maybe_data_frame(x), y = y, ntree = ~500, mtry = min_cols(~5,      x), importance = ~TRUE) 
#>                Type of random forest: regression
#>                      Number of trees: 500
#> No. of variables tried at each split: 5
#> 
#>           Mean of squared residuals: 48.45618
#>                     % Var explained: 87.48

1.3.1 特征重要性:基于均方误差(MSE)的减少

这种方法主要用于回归任务。

Mean Decrease in Accuracy (MDA) 这种方法通过衡量某个特征对整体模型预测准确性的贡献来计算其重要性。具体步骤如下:

  1. 训练模型:使用所有特征训练随机森林模型。

  2. 计算基线误差:使用训练好的模型在验证集上计算基线误差(例如,均方误差)。

  3. 扰动特征值:对于每个特征,随机打乱验证集中该特征的值,从而破坏其与目标变量的关系。 重新计算误差:使用扰动后的数据再次计算模型误差。

  4. 计算重要性:特征重要性得分等于扰动后的误差与基线误差之差。误差增加越多,说明该特征对模型预测的贡献越大。

Code
# 查看特征重要性
importance <- reg_rf_fit %>%
    extract_fit_parsnip() %>%
    vip::vi()
vip::vi(reg_rf_fit$fit$fit) 
#> # A tibble: 26 × 2
#>    Variable Importance
#>    <chr>         <dbl>
#>  1 HL             34.9
#>  2 NE_A           29.4
#>  3 Sport          21.9
#>  4 SA             14.6
#>  5 I              14.1
#>  6 SMML           13.6
#>  7 FC             13.5
#>  8 BECKNC         13.2
#>  9 YS             12.8
#> 10 MVS            12.6
#> # ℹ 16 more rows


# 可视化特征重要性
ggplot(importance, aes(x = reorder(Variable, Importance), y = Importance)) + 
    geom_col(fill = "skyblue") +
    coord_flip() +
    labs(title = "特征重要性", x = "特征", y = "重要性得分")

1.4 装袋法

装袋法(Bagging)或称自助聚合(Bootstrap Aggregation)是基于树的模型的一种集成学习技术。

装袋法是一种集成学习技术,它通过构建多个模型(通常是决策树)并将其预测结果进行平均(对于回归任务)或投票(对于分类任务)来提高模型的准确性和稳健性。其主要步骤如下:

  1. 数据抽样:从原始训练数据集中通过自助法(Bootstrap)随机有放回地抽取多个子集。每个子集的大小与原始数据集相同,但由于是有放回地抽样,因此每个子集中可能包含重复的样本。

  2. 模型训练:对每个子集训练一个模型(通常是决策树模型)。由于每个子集的样本可能不同,训练得到的每个模型也可能不同。

  3. 模型集成:在进行预测时,将所有模型的预测结果进行整合。对于分类任务,采用多数投票法,即选择出现次数最多的类别作为最终预测结果;对于回归任务,则采用平均法,即取各模型预测值的平均值作为最终预测结果。

1.5 梯度提升树

梯度提升树(Gradient Boosting Trees, GBT)是一种强大的集成学习方法,它通过逐步构建多个决策树,并将它们的预测结果进行加权组合来提高模型的预测性能。与装袋法(Bagging)不同,梯度提升树是一个迭代的过程,在每一步中都试图纠正前一步模型的错误。

梯度提升树的基本思想是通过逐步构建一系列的弱学习器(通常是决策树)来逼近目标函数。每一个新的树都在先前树的基础上进行改进,使整体模型的预测误差逐步减小。其主要步骤如下:

  1. 初始化模型:首先用一个简单的模型(如常数值模型)初始化预测值。

  2. 计算残差:计算初始模型的预测值与实际值之间的差值(残差),这些残差代表了当前模型的误差。

  3. 训练新树:基于残差训练一个新的决策树,目的是学习如何纠正当前模型的误差。

  4. 更新模型:将新树的预测结果加权加入到当前模型中,从而更新整体模型。更新公式通常为:

\[ F_m(x)=F_{m-1}(x)+ηh_m(x)F_{m}(x) \] 其中,\(F_{m}(x)\) 是第 m 次迭代的模型,\(\eta\) 是学习率(通常在0和1之间),\(h_m(x)\) 是第 m 棵树的预测值。

  1. 重复迭代:重复步骤2-4,直到达到预定的树的数量或误差收敛。

1.6 贝叶斯相加回归树

Back to top