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

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

随机森林(RF)优势

这部分可以直接跳过，直接看后面R语言实现随机森林(RF)代码。

上次介绍了逻辑回归(LR)，这里简单说一下随机森林的优势,随机森林在生物信息学中应用场景和逻辑回归类似。关于逻辑回归的内容可以点击下面链接查看。

R语言实现逻辑回归(LR)以及绘制ROC曲线和混淆矩阵

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

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

下面讲一下如何在R中实现随机森林

R语言实现RF

测试数据和代码

百度云链接：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)

本文章转载微信公众号@Bio小菜鸟