R : 随机森林（测试版1）

# 清空当前环境中的所有对象
rm (list = ls ())

# 设置工作目录
setwd("C:\\Users\\Administrator\\Desktop\\随机森林4")

library(randomForest) #随机森林
library(tidyverse) #数据分析和可视化
library(skimr) #生成数据摘要统计分析
library(DataExplorer) #探索性数据分析
library(caret) #分类和回归模型的函数
library(pROC) #生成ROC曲线和计算AUC
library(caTools) #数据分割和抽样

#加载数据，指定第一行包含列名（变量名）
otu <- read.table("otutable.txt", header = TRUE, sep = "\t")

#因变量分布情况
table(otu$gene)

# 创建训练集和测试集的索引
split_index <- sample.split(otu$gene, SplitRatio = 0.8)

# 根据索引划分数据
train_data <- subset(otu, split_index == TRUE)
test_data <- subset(otu, split_index == FALSE)

#拆分后因变量分布
table(train_data$gene)
table(test_data$gene)

#因变量自变量构建公式
colnames(otu) #获取列名
form_cls <- as.formula(
  paste0(
    "gene ~",   #构建分类模型的公式
    paste0(colnames(train_data)[2:554],collapse = "+")
  )
)
form_cls

#构建模型
set.seed(1)
# 将 gene 列转换为因子（因为这是一个分类问题）
train_data$gene <- factor(train_data$gene)
fit_rf_cls <- randomForest(
  form_cls,
  data = train_data,
  ntree = 500,
  mtry = 23,
  importance = T #是否计算变量的重要性
)
#模型概况
fit_rf_cls

#绘制 nreww 参数与误差之间的关系图
plot(fit_rf_cls, main = "ERROR & TREES")
# 添加图例
legend("top",
       legend = colnames(fit_rf_cls$err.rate),
       lty = 1:3,
       col = 1:3,
       horiz = T)

#获取变量重要性
importance(fit_rf_cls)
# 绘制变量重要性图（类型为默认值）
varImpPlot(fit_rf_cls,main = "varImpPlot")
# 绘制变量重要性图（类型为1）
varImpPlot(fit_rf_cls,main = "varImpPlot",type = 1)
# 绘制变量重要性图（类型为2）
varImpPlot(fit_rf_cls,main = "varImpPlot",type = 2)

#偏相关图
partialPlot(x=fit_rf_cls, #指定使用的模型
            pred.data = train_data,  #预测数据集
            x.var = Ramlibacter, #要绘制偏相关图的自变量
            which.class = "B73", #指定类别
            ylab = "B73") #指定 y 轴标签
# 计算变量 Dyadobacter 在类别 B73 中的分布比例
prop.table(table(train_data$Ramlibacter,train_data$gene),margin = 1)

#预测
#训练集预训练概率,返回概率而不是类别
trainpredprob <- predict(fit_rf_cls, newdata = train_data , type = "prob")
# 计算 ROC 曲线
trainroc <- roc(response = train_data$gene,
                predictor = trainpredprob[, 2])

#绘制训练集ROC曲线
plot(trainroc,
     print.auc = TRUE, #打印 AUC（Area Under the Curve）值
     auc.polygon = TRUE, #在曲线下方填充多边形，用于突出 AUC 区域
     grid = T, #显示网格线
     max.auc.polygon = T, #在最大 AUC 区域填充多边形
     auc.polygon.col = "skyblue", #指定填充多边形的颜色
     print.thres = T, #在图中标注阈值
     legacy.axes = T, #使用旧版坐标轴
     bty = "l" #指定绘图的边界类型
     )
#约登法则,计算最佳 Youden Index 时的分类阈值
bastp <- trainroc$thresholds[
  which.max(trainroc$sensitivities + trainroc$specificities - 1)
]#找到具有最大 Youden Index 的分类阈值的索引
bastp

#训练集预测分类
trainpredlab <- as.factor(
  ifelse(trainpredprob[, 2] > bastp, "Mo17", "B73")
)

#训练集混淆矩阵
confusionMatrix(data = trainpredlab, #预测类别
                reference = train_data$gene, #实际类别
                positive = "Mo17", # 正类别标签
                mode = "everything" # 显示所有评估指标
                )

#测试集预测概率
testpredprob <- predict(fit_rf_cls, newdata = test_data, type = "prob")

#测试集预测分类
testpredlab <- as.factor(
  ifelse(testpredprob[, 2] > bastp, "Mo17","B73")
)

# 将 test_data$gene 转换为因子
test_data$gene <- factor(test_data$gene)
#测试集混淆矩阵
confusionMatrix(data = testpredlab, #预测类别
                reference = test_data$gene, #实际类别
                positive = "Mo17", # 正类别标签
                mode = "everything" # 显示所有评估指标
                )
#测试集ROC
testroc <- roc(response = test_data$gene, #实际类别
                predictor = testpredprob[, 2]) #预测概率
#训练集、测试集ROC曲线叠加
plot(trainroc,
     print.auc = TRUE,
     grid = c(0.1,0.2),
     auc.polygon = F,
     max.auc.polygon = T,
     main = "随机森林--ROC",
     grid.col=c("green","red"))
plot(testroc,
     print.auc = TRUE,
     print.auc.y = 0.4,
     add = T,
     col = "red")
legend("bottomright",
       legend = c("train_data","test_data"),
       col = c(par("fg"),"red"),
       lwd = 2,
       cex = 0.9)

##############################


# 生成混淆矩阵
cm <- confusionMatrix(data = testpredlab, reference = test_data$gene, positive = "Mo17", mode = "everything")

# 可视化混淆矩阵
fourfoldplot(cm$table, color = c("#CC6666", "#99CC99"), conf.level = 0, margin = 1, main = "Confusion Matrix")