Chapter 9 第九章 机器学习建模入门

学习目标

  • 📌 了解机器学习建模中数据准备过程

  • 📌 掌握监督和非监督学习建模

  • 📌 掌握模型优化基本方法

  • 📌 掌握模型评估基本方法

9.1 数据准备

9.1.1 数据划分

9.1.2 重采样

9.1.3 数据平衡

9.1.4 特征工程

9.1.4.1 因变量特征工程

9.1.4.2 缺失值

9.1.4.3 特征筛选

9.1.4.4 自变量特征工程

9.2 算法类型

9.2.1 监督学习

9.2.2 非监督学习

9.2.2.1 降维

9.2.2.2 聚类

9.3 模型优化

9.4 模型评估

9.5 建模工具

9.6 课堂任务

任务9.1 降维技术：主成份分析与因子分析

主成分分析

因子分析

# 尝试不同因子个数
fac.3 = factanal(metrics_SCA[, -c(1:10,25)], 
                  factors = 3,
                  scores = "regression")
fac.3 
# fac.4 = factanal(metrics_SCA[, -c(1:10,25)], 
#                   factors = 4,
#                   scores = "regression")
# 
# fac.5 = factanal(metrics_SCA[, -c(1:10,25)], 
#                   factors = 5,
#                   scores = "regression")

# 方法1：用 1 - uniqueness
communality <- 1 - fac.3$uniquenesses

# 方法2：用载荷平方和（验证一致性）
loadings_matrix <- fac.3$loadings[]
communality_from_loadings <- rowSums(loadings_matrix^2)

# 比较两种方法
comparison <- data.frame(
  变量 = names(fac.3$uniquenesses),
  独特方差 = round(fac.3$uniquenesses, 4),
  方法1_公因子方差 = round(communality, 4),
  方法2_公因子方差 = round(communality_from_loadings, 4),
  是否相等 = round(communality, 4) == round(communality_from_loadings, 4)
)
print(comparison)


loadings = loadings(fac.3) 
plot(loadings, type = "n", xlim = c(-0.4, 1)) 
text(loadings, rownames(loadings), cex = 0.8)



# 尝试不同因子旋转方式
# 默认是正交 varimax
# 斜交
fac.3.rotate = factanal(metrics_SCA[, -c(1:10,25)], 
                  factors = 3,
           rotation = "promax",
                  scores = "regression")
fac.3.rotate

loadings.rot = loadings(fac.3.rotate) 
plot(loadings.rot, type = "n", xlim = c(-0.4, 1)) 
text(loadings.rot, rownames(loadings), cex = 0.8)


# 提取因子分数
# sc.fac.5=factanal(metrics_SCA[, -c(1:10,25)], 
#                   factors = 5,
#                   scores = "regression")
# 
# sc.fac.5$scores
# 
# loadings(sc.fac.5) 
# loadings(sc.fac.4) 
# loadings(sc.fac.3) 


# factor rotation
# sc.fac.rot = factanal(metrics_SCA[, -c(1:10,25)], 
#                    factors = 5, rotation = "promax") 
# 
# loadings.rot = loadings(sc.fac.rot) 
# plot(loadings.rot, type = "n", xlim = c(-0.4, 1)) 
# text(loadings.rot, rownames(loadings)) 
# abline(h = -0.1, col = "darkgrey")

任务9.2 降维技术：对应分析与多维尺度分析

对应分析

任务9.3 聚类算法

层次聚类

分割优化

计算文体学

任务9.4 监督学习-回归任务

数据集划分示例

pictureNaming <- read_csv("data/ch9/pictureNaming.new.csv")%>%
  filter(AoA_r != "#NULL!")%>%
  drop_na()%>%
  mutate(AoA_r = as.numeric(AoA_r),
         AoA_o = as.numeric(AoA_o))%>%
  select("LEN","VAR","IMG","FAM","VC","FREQ", 
         "NAgree","H_value","CA","AoA_r","AoA_o", "RT_harm")%>%
  mutate(LEN = as.factor(LEN))

# Using base R 
set.seed(123) # for reproducibility  
index_1 <- sample(1:nrow(pictureNaming), round(nrow(pictureNaming) * 0.7)) 
train_1 <- pictureNaming[index_1, ] 
test_1 <- pictureNaming[-index_1, ]  

# Using caret package 
library(caret)
set.seed(123) # for reproducibility  
index_2 <- createDataPartition(pictureNaming$RT_harm, p = 0.7, list = FALSE) 
train_2 <- pictureNaming[index_2, ] 
test_2 <- pictureNaming[-index_2, ]  

# Using rsample package 
library(rsample) # for data splitting  
set.seed(123) # for reproducibility  
split_1 <- initial_split(pictureNaming, prop = 0.7) 
train_3 <- training(split_1) 
test_3 <- testing(split_1)  

# Using h2o package  
# library(h2o) # for resampling and model training  
# h2o set-up 
# h2o.no_progress() # turn off h2o progress bars 
# h2o.init() # launch h2o  
# pictureNaming.h2o <- as.h2o(pictureNaming)
# split_2 <- h2o.splitFrame(pictureNaming, ratios = 0.7, seed = 123) 
# train_4 <- split_2[[1]] 
# test_4 <- split_2[[2]]

# 分层抽样
# orginal response distribution  
table(pictureNaming$LEN) %>% prop.table() 

# stratified sampling with the rsample package 
set.seed(123)  
split_strat <- initial_split(pictureNaming, prop = 0.7, strata = "LEN") 
train_strat <- training(split_strat) 
test_strat <- testing(split_strat)  
# consistent response ratio between train & test 
table(train_strat$LEN) %>% prop.table() 
table(test_strat$LEN) %>% prop.table()

回归任务建模

# 数据分为训练集和测试集
set.seed(123)
pictureNaming_split <- initial_split(pictureNaming, prop = 0.7) 
pictureNaming_train <- training(pictureNaming_split) 
pictureNaming_test <- testing(pictureNaming_split)



pictureNaming_m1 <- rpart(formula = RT_harm ~ ., 
                          data = pictureNaming_train, method = "anova" )


rpart.plot(pictureNaming_m1)

plotcp(pictureNaming_m1)

# rpart cross validation results 
pictureNaming_m1$cptable

# 不同特征的重要性
vip(pictureNaming_m1, num_features = 7, bar = FALSE)

# 1. 使用模型进行预测
test_predictions <- predict(pictureNaming_m1, newdata = pictureNaming_test)

# 2. 获取测试集真实值
test_actual <- pictureNaming_test$RT_harm

# 3. 创建结果数据框
results <- data.frame(
  Actual = test_actual,
  Predicted = test_predictions,
  Residual = test_actual - test_predictions
)


# 基础回归评估指标
model_performance <- data.frame(
  # 1. 误差指标
  MAE = mae(test_actual, test_predictions),          # 平均绝对误差
  RMSE = rmse(test_actual, test_predictions),        # 均方根误差
  MAPE = mape(test_actual, test_predictions) * 100,  # 平均绝对百分比误差
  
  # 2. 拟合优度指标
  R_squared = cor(test_actual, test_predictions)^2,  # R平方
  Adjusted_R2 = 1 - (1 - cor(test_actual, test_predictions)^2) * 
    (nrow(pictureNaming_test) - 1) / 
    (nrow(pictureNaming_test) - ncol(pictureNaming_test) - 1),
  
  # 3. 相对指标
  NRMSE = rmse(test_actual, test_predictions) / mean(test_actual),  # 标准化RMSE
  
  # 4. 样本统计
  Sample_Size = nrow(pictureNaming_test),
  Predictor_Count = ncol(pictureNaming_test) - 1  # 减去因变量
)

任务9.5 监督学习分类任务

9.6.0.1 决策树模型

## model 1
datasetAE.tree1 <- rpart(formula = focus_type ~ ., 
                         data = datasetAE_train)

rpart.plot(datasetAE.tree1)

plotcp(datasetAE.tree1)

# 模型评估
# 对测试集进行预测
predictions <- predict(datasetAE.tree1, 
                       newdata = datasetAE_test, 
                       type = "class")
probabilities <- predict(datasetAE.tree1, 
                         newdata = datasetAE_test, 
                         type = "prob")



# 创建混淆矩阵
confusion_matrix <- confusionMatrix(predictions, 
                                   datasetAE_test$focus_type)
print(confusion_matrix)

# 查看详细的评估指标
confusion_matrix$overall  # 总体准确率
confusion_matrix$byClass  # 各类别的详细指标（精确率、召回率、F1等）


# model2
datasetAE.tree2 <- rpart( formula = focus_type ~ ., 
                         data = datasetAE_train,
                         control = list(cp = 0, xval = 4) )
rpart.plot(datasetAE.tree2)
plotcp(datasetAE.tree2)

# model 3
# caret cross validation results 
datasetAE.tree3 <- train( focus_type ~ ., 
                          data = datasetAE_train,
                          method = "rpart", 
                   trControl = trainControl(method = "cv", number = 10), 
                   tuneLength = 20 )  

ggplot(datasetAE.tree3)

# feature importance
vip(datasetAE.tree3, num_features = 21, bar = FALSE)

# Construct partial dependence plots  
# p1 <- partial(datasetAE.tree3, pred.var = "Fallexc") %>% autoplot() 
# p2 <- partial(datasetAE.tree3, pred.var = "Falldur") %>% autoplot() 
# p3 <- partial(datasetAE.tree3, pred.var = c("Fallexc", "Falldur")) %>% 
# 
# plotPartial(levelplot = FALSE, zlab ="yhat", 
#               drape = TRUE, colorkey = TRUE, 
#               screen = list(z = -20, x = -60))  
# # Display plots side by side  
# gridExtra::grid.arrange(p1, p2, p3, ncol = 3)

random forest

## package 1 = randomForest


set.seed(123)

RF_AE = randomForest(focus_type ~ .,
                    ntree = 500, 
                    data = datasetAE_train, 
                    importance = TRUE)
importance(RF_AE)   

varImpPlot(RF_AE, n.var= 5, main ="AE")

# 在构建随机森林的过程中，每棵决策树都是通过对原始训练数据进行有放回抽样（Bootstrap Sampling）得到的子样本训练出来的。剩下样本则没有被该树用于训练，这些样本就是袋外样本。利用那些没有进行训练的数据来进行预测。

as.data.frame(RF_AE$confusion)%>%
mutate(across(where(is.numeric), ~ round(.x, digits = 2)))%>%

  flextable()%>%
  set_caption(caption = "Table x .")


# 模型评估

RF_predictions <- predict(RF_AE, newdata = datasetAE_test)
RF_probabilities <- predict(RF_AE, newdata = datasetAE_test, type = "prob")

RF_cm <- confusionMatrix(RF_predictions, datasetAE_test$focus_type)
# 4. 可视化混淆矩阵
conf_matrix_df <- as.data.frame(RF_cm$table)
ggplot(conf_matrix_df, aes(x = Reference, y = Prediction, fill = Freq)) +
  geom_tile(color = "white") +
  geom_text(aes(label = Freq), color = "black", size = 4) +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(title = "混淆矩阵热力图",
       x = "真实类别",
       y = "预测类别",
       fill = "计数") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


# 提取每个类别的指标
class_metrics <- RF_cm$byClass
print(round(class_metrics, 4))



# 创建汇总表
metrics_summary <- data.frame(
  指标 = c("总体准确率", "Kappa系数", "微平均F1", "宏平均F1", "加权平均F1"),
  值 = c(
    RF_cm$overall["Accuracy"],
    RF_cm$overall["Kappa"],
    # 微平均F1（等于准确率）
    RF_cm$overall["Accuracy"],
    # 宏平均F1
    mean(class_metrics[, "F1"], na.rm = TRUE),
    # 加权平均F1
    weighted.mean(class_metrics[, "F1"], table(datasetAE_test$focus_type), na.rm = TRUE)
  )
)
print(metrics_summary)
## ranger package

# number of features  
n_features <- length(setdiff(names(datasetAE_train), "focus_type"))  

# train a default random forest model 
datasetAE_rf1 <- ranger( focus_type ~ ., data = datasetAE_train, 
                    mtry = floor(n_features / 3), 
                    respect.unordered.factors = "order",
                    seed = 123 )  
# get OOB RMSE  
default_rmse <- sqrt(datasetAE_rf1$prediction.error)
# create hyperparameter grid 
hyper_grid <- expand.grid( mtry = floor(n_features * c(.05, .15, .25, .333, .4)), 
                           min.node.size = c(1, 3, 5, 10), 
                           replace = c(TRUE, FALSE), 
                           sample.fraction = c(.5, .63, .8), rmse = NA )  

# execute full cartesian grid search 
for(i in seq_len(nrow(hyper_grid))) {  
  # fit model for ith hyperparameter combination 
fit <- ranger(  formula = focus_type ~ ., 
                  data = datasetAE_train, 
                  num.trees = n_features * 10, 
                  mtry = hyper_grid$mtry[i], 
                  min.node.size = hyper_grid$min.node.size[i], 
                  replace = hyper_grid$replace[i], 
                  sample.fraction = hyper_grid$sample.fraction[i], 
                  verbose = FALSE, seed = 123, respect.unordered.factors = "order", )  
  # export OOB error  
hyper_grid$rmse[i] <- sqrt(fit$prediction.error) 
  }
 
# assess top 10 models 
hyper_grid %>%  
  arrange(rmse) %>%  
  mutate(perc_gain = (default_rmse - rmse) / default_rmse * 100) %>% 
  head(10)%>% mutate(across(where(is.numeric), ~ round(.x, digits = 2)))%>%
  flextable()%>%
  set_caption(caption = "Table x .")

# re-run model with impurity-based variable importance
rf_impurity <- ranger(
  formula = focus_type ~ ., 
  data = datasetAE_train, 
  num.trees = 2000,
  mtry = 20,
  min.node.size = 1,
  sample.fraction = .80,
  replace = FALSE,
  importance = "impurity",
  respect.unordered.factors = "order",
  verbose = FALSE,
  seed  = 123
)

# re-run model with permutation-based variable importance
rf_permutation <- ranger(
  formula = focus_type ~ ., 
  data = datasetAE_train, 
  num.trees = 2000,
  mtry = 20,
  min.node.size = 1,
  sample.fraction = .80,
  replace = FALSE,
  importance = "permutation",
  respect.unordered.factors = "order",
  verbose = FALSE,
  seed  = 123
)

p1 <- vip::vip(rf_impurity, num_features = 25, bar = FALSE)
p2 <- vip::vip(rf_permutation, num_features = 25, bar = FALSE)

gridExtra::grid.arrange(p1, p2, nrow = 1)
# pacakge 3 = h2o

h2o.no_progress()
h2o.init(max_mem_size = "5g")

# convert training data to h2o object
train_h2o <- as.h2o(datasetAE_train)

# set the response column to Sale_Price
response <- "focus_type"

# set the predictor names
predictors <- setdiff(colnames(datasetAE_train), response)

h2o_rf1 <- h2o.randomForest(
    x = predictors, 
    y = response,
    training_frame = train_h2o, 
    ntrees = n_features * 10,
    seed = 123
)

h2o_rf1

# hyperparameter grid
# hyper_grid <- list(
#   mtries = floor(n_features * c(.05, .15, .25, .333, .4)),
#   min_rows = c(1, 3, 5, 10),
#   max_depth = c(10, 20, 30),
#   sample_rate = c(.55, .632, .70, .80)
# )

# random grid search strategy
# search_criteria <- list(
#   strategy = "RandomDiscrete",
#   stopping_metric = "mse",
#   stopping_tolerance = 0.001,   # stop if improvement is < 0.1%
#   stopping_rounds = 10,         # over the last 10 models
#   max_runtime_secs = 60*5      # or stop search after 5 min.
# )


# # perform grid search 
# random_grid <- h2o.grid(
#   algorithm = "randomForest",
#   grid_id = "rf_random_grid",
#   x = predictors, 
#   y = response, 
#   training_frame = train_h2o,
#   hyper_params = hyper_grid,
#   ntrees = n_features * 10,
#   seed = 123,
#   stopping_metric = "RMSE",   
#   stopping_rounds = 10,           # stop if last 10 trees added 
#   stopping_tolerance = 0.005,     # don't improve RMSE by 0.5%
#   search_criteria = search_criteria
# )


# collect the results and sort by our model performance metric 
# of choice
# random_grid_perf <- h2o.getGrid(
#   grid_id = "rf_random_grid", 
#   sort_by = "mse", 
#   decreasing = FALSE
# )
# 
# random_grid_perf

Gradient Boosting

h2o.init(max_mem_size = "10g")

train_h2o <- as.h2o(datasetAE_train)
response <- "focus_type"
predictors <- setdiff(colnames(datasetAE_train), response)

# run a basic GBM model
set.seed(123)  # for reproducibility

datasetAE_gbm1 <- gbm(
  formula = focus_type ~ .,
  data = datasetAE_train,
  distribution = "gaussian",  # SSE loss function
  n.trees = 5000,
  shrinkage = 0.1,
  interaction.depth = 3,
  n.minobsinnode = 10,
  cv.folds = 10
)

# find index for number trees with minimum CV error
best <- which.min(datasetAE_gbm1 $cv.error)

# get MSE and compute RMSE
sqrt(datasetAE_gbm1$cv.error[best])

# plot error curve
gbm.perf(datasetAE_gbm1, method = "cv")


# create grid search
hyper_grid <- expand.grid(
  learning_rate = c(0.3, 0.1, 0.05, 0.01, 0.005),
  RMSE = NA,
  trees = NA,
  time = NA
)

# execute grid search
for(i in seq_len(nrow(hyper_grid))) {

  # fit gbm
  set.seed(123)  # for reproducibility
  train_time <- system.time({
    m <- gbm(
      formula = focus_type ~ .,
      data = datasetAE_train,
      distribution = "gaussian",
      n.trees = 5000, 
      shrinkage = hyper_grid$learning_rate[i], 
      interaction.depth = 3, 
      n.minobsinnode = 10,
      cv.folds = 10 
   )
  })
  
  # add SSE, trees, and training time to results
  hyper_grid$RMSE[i]  <- sqrt(min(m$cv.error))
  hyper_grid$trees[i] <- which.min(m$cv.error)
  hyper_grid$Time[i]  <- train_time[["elapsed"]]

}

# results
arrange(hyper_grid, RMSE)

# second try
# search grid
hyper_grid <- expand.grid(
  n.trees = 6000,
  shrinkage = 0.01,
  interaction.depth = c(3, 5, 7),
  n.minobsinnode = c(5, 10, 15)
)

# create model fit function
model_fit <- function(n.trees, shrinkage, interaction.depth, n.minobsinnode) {
  set.seed(123)
  m <- gbm(
    formula = focus_type ~ .,
    data = datasetAE_train,
    distribution = "gaussian",
    n.trees = n.trees,
    shrinkage = shrinkage,
    interaction.depth = interaction.depth,
    n.minobsinnode = n.minobsinnode,
    cv.folds = 10
  )
  # compute RMSE
  sqrt(min(m$cv.error))
}

# perform search grid with functional programming
hyper_grid$rmse <- purrr::pmap_dbl(
  hyper_grid,
  ~ model_fit(
    n.trees = ..1,
    shrinkage = ..2,
    interaction.depth = ..3,
    n.minobsinnode = ..4
    )
)

# results
arrange(hyper_grid, rmse)

# implement a stochastic GBM

# refined hyperparameter grid
hyper_grid <- list(
  sample_rate = c(0.5, 0.75, 1),              # row subsampling
  col_sample_rate = c(0.5, 0.75, 1),          # col subsampling for each split
  col_sample_rate_per_tree = c(0.5, 0.75, 1)  # col subsampling for each tree
)

# random grid search strategy
search_criteria <- list(
  strategy = "RandomDiscrete",
  stopping_metric = "mse",
  stopping_tolerance = 0.001,   
  stopping_rounds = 10,         
  max_runtime_secs = 60*60      
)

# perform grid search 
grid <- h2o.grid(
  algorithm = "gbm",
  grid_id = "gbm_grid",
  x = predictors, 
  y = response,
  training_frame = train_h2o,
  hyper_params = hyper_grid,
  ntrees = 6000,
  learn_rate = 0.01,
  max_depth = 7,
  min_rows = 5,
  nfolds = 10,
  stopping_rounds = 10,
  stopping_tolerance = 0,
  search_criteria = search_criteria,
  seed = 123
)

# collect the results and sort by our model performance metric of choice
grid_perf <- h2o.getGrid(
  grid_id = "gbm_grid", 
  sort_by = "mse", 
  decreasing = FALSE
)

grid_perf


# Grab the model_id for the top model, chosen by cross validation error
best_model_id <- grid_perf@model_ids[[1]]
best_model <- h2o.getModel(best_model_id)

# Now let’s get performance metrics on the best model
h2o.performance(model = best_model, xval = TRUE)

#### XGBoost

xgb_prep <- recipe(focus_type ~ ., data = datasetAE_train) %>%
  step_integer(all_nominal()) %>%
  prep(training = datasetAE_train, retain = TRUE) %>%
  juice()

X <- as.matrix(xgb_prep[setdiff(names(xgb_prep), "focus_type")])
Y <- xgb_prep$focus_type

set.seed(123)
datasetAE_xgb <- xgb.cv(
  data = X,
  label = Y,
  nrounds = 6000,
  objective = "reg:linear",
  early_stopping_rounds = 50, 
  nfold = 10,
  params = list(
    eta = 0.1,
    max_depth = 3,
    min_child_weight = 3,
    subsample = 0.8,
    colsample_bytree = 1.0),
  verbose = 0
)  

# minimum test CV RMSE
min(datasetAE_xgb$evaluation_log$test_rmse_mean)

# hyperparameter grid
hyper_grid <- expand.grid(
  eta = 0.01,
  max_depth = 3, 
  min_child_weight = 3,
  subsample = 0.5, 
  colsample_bytree = 0.5,
  gamma = c(0, 1, 10, 100, 1000),
  lambda = c(0, 1e-2, 0.1, 1, 100, 1000, 10000),
  alpha = c(0, 1e-2, 0.1, 1, 100, 1000, 10000),
  rmse = 0,          # a place to dump RMSE results
  trees = 0          # a place to dump required number of trees
)

# optimal parameter list
params <- list(
  eta = 0.01,
  max_depth = 3,
  min_child_weight = 3,
  subsample = 0.5,
  colsample_bytree = 0.5
)

# train final model
xgb.fit.final <- xgboost(
  params = params,
  data = X,
  label = Y,
  nrounds = 3944,
  objective = "reg:linear",
  verbose = 0
)

SVM

# Tune an SVM with radial basis kernel 
set.seed(1854) # for reproducibility 
datasetAE_svm <- train( focus_type ~ ., 
                    data = datasetAE_train, 
                    method = "svmRadial", 
                    preProcess = c("center", "scale"), 
                    trControl = trainControl(method = "cv", number = 10), 
                    tuneLength = 10 )

# Plot results  
ggplot(datasetAE_svm ) + 
  theme_light()

as.data.frame(datasetAE_svm$results)%>%
mutate(across(where(is.numeric), ~ round(.x, digits = 2)))%>%
  flextable()%>%
  set_caption(caption = "Table x .")

# Control params for SVM 
ctrl <- trainControl( method = "cv", number = 10, 
                      classProbs = TRUE, 
                      summaryFunction = multiClassSummary # also needed for AUC/ROC 
                      )  
# Tune an SVM 
set.seed(5628) # for reproducibility 
datasetAE_svm_auc <- train( focus_type ~ .,
                        data = datasetAE_train, 
                        method = "svmRadial", 
                        preProcess = c("center", "scale"),  
                        metric = "ROC", # area under ROC curve (AUC)
                        trControl = ctrl, tuneLength = 10 )

confusionMatrix(datasetAE_svm_auc)

# Variable importance plot 
set.seed(2827) # for reproducibility  
# prob_bro <- function(object, newdata) { predict(object, 
#                                                 newdata = newdata, 
#                                                 type = "prob")[, "bro"] }
# vip(datasetAE_svm_auc, 
#     method = "permute", 
#     nsim = 5, 
#     train = datasetAE_train, 
#     target = "focus_type", 
#     metric = "roc_auc", 
#     reference_class = "bro", 
#     pred_wrapper = prob_bro)

vip(
  datasetAE_svm_auc, 
  method = "permute", 
  nsim = 5, 
  train = datasetAE_train, 
  target = "focus_type", 
  metric = yardstick::accuracy_vec, 
  pred_wrapper = function(object, newdata) predict(object, newdata = newdata),
  smaller_is_better = FALSE
)

# vip::list_metrics()
# 
# levels(datasetAE_train$focus_type)