R语言实现随机森林(RF)以及绘制ROC曲线和混淆矩阵

R语言实现随机森林(RF)以及绘制ROC曲线和混淆矩阵

随机森林是一种集成学习方法，通过构建多个决策树并综合它们的预测结果来进行分类或回归。它在处理非线性关系、鲁棒性和抗噪性方面具有显著优势，特别适合处理高维数据和包含缺失值的数据集。本文将详细介绍如何在R语言中实现随机森林算法，并绘制ROC曲线和混淆矩阵来评估模型性能。

随机森林算法简介

随机森林（Random Forest）是一种集成学习方法，由多个决策树构成的模型。它通过对训练数据进行自助采样（bootstrap sampling）和特征随机选择（random feature selection）来构建多个决策树，并最终通过投票或平均预测结果来进行分类或回归。

随机森林的优势

• 处理非线性关系：随机森林可以捕捉更复杂的非线性关系。它由多个决策树组成，每个决策树可以学习和表示不同的特征组合，从而更好地处理非线性关系。而逻辑回归是线性模型，只能捕捉线性关系。
• 鲁棒性和抗噪性：随机森林对于噪声和异常值具有一定的鲁棒性。由于采用自助采样和多个决策树的投票或平均，随机森林可以减少个别样本的影响，从而更具鲁棒性。逻辑回归相对较为脆弱，受异常值和噪声的影响较大。
• 处理缺失值：随机森林可以有效地处理含有缺失值的数据。随机森林利用随机特征选择的方式，可以在不使用缺失值的情况下对样本进行预测。逻辑回归在处理缺失值时需要对缺失值进行特殊处理。

但是逻辑回归也具有其自身的优势，如模型解释性好、计算效率高等。在实际应用中，选择使用哪种方法应根据数据特点、问题需求以及性能评估等因素进行综合考虑。

在R中实现随机森林

测试数据和代码

百度云链接：https://pan.baidu.com/s/14xTrOn3R7ClU7i4wZ7ZYew 提取码：6lli

# 导入必要的包,没有安装的可以先安装一下  
library(dplyr) #数据处理使用  
library(data.table) #数据读取使用  
library(randomForest) #RF模型使用  
library(caret) # 调参和计算模型评价参数使用  
library(pROC) #绘图使用  
library(ggplot2) #绘图使用  
library(ggpubr) #绘图使用  
library(ggprism) #绘图使用  

# 读取数据  
data <- fread("./RF_data.txt",data.table = F)  # 替换为你的数据文件名或路径  

数据长这个样子，一共35727行，214列。每一行代表一个样本，第一列是样本标签（malignant或normal），后面213列是213个特征。我们想根据213个特征，使用RF训练出一个能够对样本进行精准分类的模型。

构建RF模型

# 分割数据为训练集和测试集  
set.seed(123)  # 设置随机种子，保证结果可复现  
split <- sample.split(data$type, SplitRatio = 0.8)  # 将数据按照指定比例分割  
train_data <- subset(data, split == TRUE)  # 训练集  
test_data <- subset(data, split == FALSE)  # 测试集  

# 定义训练集特征和目标变量  
X_train <- train_data[, -1]  
y_train <- as.factor(train_data[, 1])  

# 创建随机森林分类模型  
model <- randomForest(x = X_train, y = y_train, ntree = 100)  

# 输出默认参数下的模型性能  
print(model)  

这是一个没有经过调参的模型结果，尽管看起来模型已经很不错了，但我们还是继续进行调参，看一下模型效果能够上升多少。

使用caret包进行调参

有至少两个参数需要进行测序，分别是mtry和ntree。

• mtry参数：它表示每棵决策树在进行节点分裂时考虑的特征数量。默认情况下，mtry的取值是平方根（对于分类问题）或总特征数的三分之一（对于回归问题）。mtry的选择影响了随机森林的多样性和模型的复杂度。较小的mtry值可能导致过拟合；较大的mtry值可能导致欠拟合。通常，可以通过交叉验证或网格搜索等技术来选择最佳的mtry值。
• ntree参数：它表示随机森林中决策树的数量。ntree的值决定了随机森林中决策树的个数，也可以说是集成模型中基学习器的数量。增加ntree的值可以增加随机森林的稳定性和预测准确性，但也会增加计算开销。通常，可以通过交叉验证或学习曲线等方法来选择合适的ntree值。

由于caret包只提供了mtry参数的调节，关于ntree参数的调节我们这里手动进行。

mtry参数调节

# 进行参数调优  
# 创建训练控制对象  
ctrl <- trainControl(method = "cv", number = 5) #使用五折交叉验证，也可以选择10折交叉验证。  

# 定义参数网格  
grid <- expand.grid(mtry = c(2, 4, 6))  # 每棵树中用于分裂的特征数量，这里只是随便给的测试，主要为了介绍如何调参，并非最优选择。  

# 使用caret包进行调参  
rf_model <- train(x = X_train, y = y_train,  
                  method = "rf",  
                  trControl = ctrl,  
                  tuneGrid = grid)  

# 输出最佳模型和参数  
print(rf_model)  

使用准确性来选择最佳模型。该模型最终mtry值为mtry = 6。

ntree参数调节

# 调整Caret没有提供的参数  
# 如果我们想调整的参数Caret没有提供，可以用下面的方式自己手动调参。  
# 用刚刚调参的最佳mtry值固定mtry  
grid <- expand.grid(mtry = c(6))  # 每棵树中用于分裂的特征数量  

# 定义模型列表，存储每一个模型评估结果  
modellist <- list()  

# 调整的参数是决策树的数量  
for (ntree in c(100,200, 300)) {  
  set.seed(123)  
  fit <- train(x = X_train, y = y_train, method="rf",   
               metric="Accuracy", tuneGrid=grid,   
               trControl=ctrl, ntree=ntree)  
  key <- toString(ntree)  
  modellist[[key]] <- fit  
}  

# compare results  
results <- resamples(modellist)  

# 输出最佳模型和参数  
summary(results)  

从准确性可以看出，ntree = 200是最佳的。这样我们就完成了调参，最佳的参数组合是mtry = 6，ntree = 200。

使用最佳参数训练模型

# 使用最佳参数训练最终模型  
final_model <- randomForest(x = X_train, y = y_train,mtry = 6,ntree = 200)  

# 输出最终模型  
print(final_model)  

从结果可以看出，经过调参的模型比初始模型好了一点点。

应用于测试集

这里使用caret包中的函数来输出模型的评价指标，想手动计算可以参考逻辑回归(LR)的推文。

# 在测试集上进行预测  
X_test <- test_data[, -1]  
y_test <- as.factor(test_data[, 1])  
test_predictions <- predict(final_model, newdata = test_data)  

# 计算模型指标  
confusion_matrix <- confusionMatrix(test_predictions, y_test)  
accuracy <- confusion_matrix$overall["Accuracy"]  
precision <- confusion_matrix$byClass["Pos Pred Value"]  
recall <- confusion_matrix$byClass["Sensitivity"]  
f1_score <- confusion_matrix$byClass["F1"]  

# 输出模型指标  
print(confusion_matrix)  
print(paste("Accuracy:", accuracy))  
print(paste("Precision:", precision))  
print(paste("Recall:", recall))  
print(paste("F1 Score:", f1_score))  

从测试集中看，模型表现的也不错。下面我们来绘制一下混淆矩阵，ROC曲线。

绘制混淆矩阵

混淆矩阵提供了对分类模型性能的全面评估。它展示了实际类别和预测类别之间的对应关系，可以清晰地看到模型的预测结果中真正例、真反例、假正例和假反例的数量或比例。绘制混淆矩阵可以将复杂的分类结果以直观的方式展示出来，使得结果更易于理解和解释。

# 绘制混淆矩阵热图  
# 将混淆矩阵转换为数据框  
confusion_matrix_df <- as.data.frame.matrix(confusion_matrix$table)  
colnames(confusion_matrix_df) <- c("cluster1","cluster2")  
rownames(confusion_matrix_df) <- c("cluster1","cluster2")  
draw_data <- round(confusion_matrix_df / rowSums(confusion_matrix_df),2)  
draw_data$real <- rownames(draw_data)  
draw_data <- melt(draw_data)  

ggplot(draw_data, aes(real,variable, fill = value)) +  
  geom_tile() +  
  geom_text(aes(label = scales::percent(value))) +  
  scale_fill_gradient(low = "#F0F0F0", high = "#3575b5") +  
  labs(x = "True", y = "Guess", title = "Confusion matrix") +  
  theme_prism(border = T)+  
  theme(panel.border = element_blank(),  
        axis.ticks.y = element_blank(),  
        axis.ticks.x = element_blank(),  
        legend.position="none")  

绘制ROC曲线

ROC（Receiver Operating Characteristic）曲线和AUC（Area Under the Curve）是评估二分类模型性能常用的指标。

• ROC曲线是以模型的真阳性率（True Positive Rate，也称为灵敏度或召回率）为纵坐标，以模型的假阳性率（False Positive Rate）为横坐标所绘制的曲线。ROC曲线能够综合考虑分类模型在不同阈值下的预测能力。
• AUC是ROC曲线下的面积，代表模型的整体性能。AUC值范围在0.5到1之间，其中0.5表示模型的预测性能等同于随机预测，而1表示模型的预测性能完美。通常来说，AUC值越接近1，模型的性能越好。

# 绘制ROC曲线需要将预测结果以概率的形式输出  
test_predictions <- predict(final_model, newdata = test_data,type = "prob")  

# 计算ROC曲线的参数  
roc_obj <- roc(response = y_test, predictor = test_predictions[, 2])  
roc_auc <- auc(roc_obj)  

# 将ROC对象转换为数据框  
roc_data <- data.frame(1 - roc_obj$specificities, roc_obj$sensitivities)  

# 绘制ROC曲线  
ggplot(roc_data, aes(x = 1 - roc_obj$specificities, y = roc_obj$sensitivities)) +  
  geom_line(color = "#0073C2FF", size = 1.5) +  
  geom_segment(aes(x = 0, y = 0, xend = 1, yend = 1), linetype = "dashed", color = "gray") +  
  geom_text(aes(x = 0.8, y = 0.2, label = paste("AUC =", round(roc_auc, 2))), size = 4, color = "black") +  
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +  
  theme_pubr() +  
  labs(x = "1 - Specificity", y = "Sensitivity") +  
  ggtitle("ROC Curve") +  
  theme(plot.title = element_text(size = 14, face = "bold"))+  
  theme_prism(border = T)