# 随机森林多分类模型 # packages 安装 —— install.packages("") setwd("D:/Rpackages") # packages 载入 —— library() library(randomForest) library(tidyverse) library(skimr) library(DataExplorer) library(caret) library(pROC) library(ggplot2) library(splines) library(reshape2) library(scales) library(ggprism) library(ggpubr) # 加载数据 boston <- read.csv(file.choose()) # 数据框转换 boston <- as.data.frame(boston) # boston <- na.pass(boston) # boston <- na.omit(boston) # na.action = na.omit 返回删除所有包含NA值的行的数据框或矩阵 # na.action = na.pass 将NA值传递给后续的计算，可能会导致结果中也出现NA # 查看数据概貌 skim(boston) # 查看数据缺失情况 plot_missing(boston) boston[is.na(boston)] <- 0 #将matrix（此处为boston）中的NA值替换为0，按需运行 # 数据类型修正，将boston中的i列转为factor，i 将依次取c()中的值 for (i in c(1:2)) { boston[,i] <- factor(boston[,i]) } # 因变量分布情况 table(boston$Class) # 拆分数据 set.seed(1) trains <- createDataPartition( y = boston$Class, p = 0.7, list = F ) traindata <- boston[trains,] testdata <- boston[-trains,] # # test：自选训练集测试集（无特殊可不理 # traindata <- read.csv(file.choose()) # traindata <- as.data.frame(traindata) # for (i in c(1:5)) { # traindata[,i] <- factor(traindata[,i]) # } # # testdata <- read.csv(file.choose()) # testdata <- as.data.frame(testdata) # for (i in c(1:5)) { # testdata[,i] <- factor(testdata[,i]) # } # 拆分后因变量分布 table(traindata$Class) table(testdata$Class) dim(boston) # # 构建公式，自变量列名~a到b列（colnames[a:b]） # # colnames(boston) 用于获取数据框的列名， # form_clsm <- as.formula( # paste0( # "class~", # paste(colnames(traindata)[3:112],collapse = "+") # ) # ) # form_clsm # # # 构建模型 # set.seed(100001) # fit_rf_clsm <- randomForest( # form_clsm, # data = traindata, # ntree = 500, # 默认500，error & trees不稳定需增加树的数量 # mtry = 6, # 每个节点可提供选择的变量数目（指定节点中用于二叉树的变量个数 # # 默认情况下数据集变量个数的二次方根(分类模型)或三分之一(预测模型)） # importance = T # importance = TRUE: 启用特征重要性计算 # # 随机森林会根据每个特征对模型性能的影响，计算出每个特征的重要性 # ) # print(fit_rf_clsm) # 定义训练集特征和目标变量 X_train <- traindata[, -(1:2)] #traindata中除第1、2列外的列为自变量 x y_train <- as.factor(traindata[, 2]) #traindata中第2列为因变量 y # # 创建随机森林分类模型（基准模型 # model <- randomForest(x = X_train, y = y_train, ntree = 500) # # # 创建训练控制对象 # ctrl <- trainControl(method = "repeatedcv", number = 5, repeats = 10) # # k 折交叉验证作为模型评估方法,小数据做重复10次5折交叉验证 # # # 进行参数调优 # # mtry参数调节 # grid1 <- expand.grid(mtry = c(48:52)) # 定义mtry范围 # # grid1 <- expand.grid(mtry = c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) # # 进阶 # # num_features <- ncol(X_train) # # mtry_range <- floor(sqrt(num_features)) # 取平方根并向下取整 # # grid_mtry <- expand.grid(mtry = seq(max(2, mtry_range - 5), # # min(mtry_range + 5, num_features))) # # # 使用caret包进行调参 # rf_model <- train(x = X_train, y = y_train, # 自变量 x， 因变量 y # method = "rf", # 使用的模型方法为随机森林（rf） # trControl = ctrl, # tuneGrid = grid1,) # # 输出结果 # print(rf_model) # # # # 定义最佳mtry参数 # grid2 <- expand.grid(mtry = c(30)) # # # 定义模型列表，存储每一个模型评估结果 # modellist <- list() # # 调整的参数是决策树的数量 # for (ntree in seq(100, 2000, by=100)) { # seq函数构建一个集合，间距by = 100 # set.seed(101) # fit <- train(x = X_train, y = y_train, method="rf", # metric="Accuracy", tuneGrid=grid2, # trControl=ctrl, ntree=ntree) # key <- toString(ntree) # modellist[[key]] <- fit # } # # # compare results # results <- resamples(modellist) # # # 输出最佳模型和参数 # summary(results) # # # 可视化模型性能，accuracy越近1越好，kappa越近1越好 # bwplot(results) # 箱线图 # dotplot(results) # 点图 # densityplot(results) # 密度图 # 使用最佳参数训练模型 set.seed(1001) fit_rf_clsm <- randomForest(x = X_train, y = y_train, mtry = 50, ntree = 500, importance = T) print(fit_rf_clsm) # ntree参数与error之间的关系图示 plot(fit_rf_clsm,main = "ERROR & TREES") legend("top", legend = colnames(fit_rf_clsm$err.rate), lty = 1:6, col = 1:6, horiz = T, cex = 0.9) plot(randomForest::margin(fit_rf_clsm), main = '观测值被判断正确的概率图') # 变量重要性 varImpPlot(fit_rf_clsm,main ="varImpPlot") varImpPlot(fit_rf_clsm,main = "varImpPlot",type = 1) varImpPlot(fit_rf_clsm,main = "varImpPlot",type = 2) importance_genus <- data.frame(importance(fit_rf_clsm)) importance_genus <- importance_genus[order(importance_genus$MeanDecreaseGini, decreasing=TRUE),] importance_genus <- importance_genus[order(importance_genus$MeanDecreaseAccuracy, decreasing=TRUE),] head(importance_genus) write.table(importance_genus,"importance_genus_DaChuang_SiteClass_unexposed20250209.txt", sep = '\t',col.names = NA,quote = FALSE) # 保存重要特征到文件 # # 预测 # # 训练集预测概率 # trainpredprob <- predict(fit_rf_clsm,newdata = traindata,type = "prob") # # 训练集ROC # multiclass.roc(response = traindata$class,predictor = trainpredprob) # # 训练集预测分类 # trainpredlab <- predict(fit_rf_clsm,newdata = traindata,type = "Class") # # 训练集混淆矩阵 # confusionMatrix_train <- confusionMatrix(data = trainpredlab, # reference = traindata$Class, # mode = "everything") # # 训练集综合结果 # multiClassSummary( # data.frame(obs = traindata$Class,pred=trainpredlab), # lev = levels(traindata$Class) # ) # # # # 作图展示 top N 重要的 OTUs # varImpPlot(fit_rf_clsm, n.var = min(20, nrow(fit_rf_clsm$importance)), # main = 'Top 20 - variable importance') # # # 测试集预测概率 # # X_test <- testdata[, -(1:2)] #testdata中除第1、2列外的列为自变量 x # # y_test <- as.factor(testdata[, 2]) #testdata中第2列为因变量 y # testpredprob <- predict(fit_rf_clsm,newdata = testdata,type = "prob") # # 测试集ROC # multiclass.roc(response = testdata$Class,predictor = testpredprob) # # 测试集预测分类 # testpredlab <- predict(fit_rf_clsm,newdata = testdata,type = "Class") # # 测试集混淆矩阵 # confusionMatrix_test <- confusionMatrix(data = testpredlab, # reference = testdata$Class, # mode = "everything") # confusionMatrix_test # # 测试集综合结果 # multiClassSummary( # data.frame(obs = testdata$Class,pred=testpredlab), # lev = levels(testdata$Class) # ) # 交叉验证帮助选择特定数量的特征 # 5次重复十折交叉验证 set.seed(10001) otu_train.cv <- replicate(5, rfcv(traindata[,-(1:2)], # 除去第a到b列(a:b) traindata$Class, cv.fold = 10, step = 1.5), simplify = FALSE) otu_train.cv <- data.frame(sapply(otu_train.cv, '[[', 'error.cv')) otu_train.cv$otus <- rownames(otu_train.cv) otu_train.cv <- reshape2::melt(otu_train.cv, id = 'otus') otu_train.cv$otus <- as.numeric(as.character(otu_train.cv$otus)) otu_train.cv.mean <- aggregate(otu_train.cv$value, by = list(otu_train.cv$otus), FUN = mean) head(otu_train.cv.mean, 18) # 绘图观察拐点 p <- ggplot(otu_train.cv,aes(otus,value)) + geom_smooth(se = FALSE, method = 'glm',formula = y~ns(x,6)) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black',fill = 'transparent')) + labs(title = '',x='Number of genus',y='Cross-validation error') p # 在横坐标(xintecept)绘制竖线 p + geom_vline(xintercept = 160) # # 备用 # p2 <- ggplot(otu_train.cv,aes(otus,value)) + # geom_line() + # theme(panel.grid = element_blank(), # panel.background = element_rect(color = 'black', fill = 'transparent')) + # labs(title = '',x = 'Number of OTUs', y = 'Cross-validation error') # p2 # # p2 + geom_vline(xintercept = 30) # 大约提取前 N个重要的特征 importance_genus[1:160,] # importance_Ngenus <- importance_genus[1:160,] # # 输出表格 # write.table(importance_genus[1:70, ], # 'importance_genus_top70_of_zhangxiaofeng_human_drowning.txt', # sep = '\t', col.names = NA, quote = FALSE) # # 变量重要性 # varImpPlot(fit_rf_clsm,main ="varImpPlot") # varImpPlot(fit_rf_clsm,main = "varImpPlot",type = 1) # varImpPlot(fit_rf_clsm,main = "varImpPlot",type = 2) # # varImpPlot(fit_rf_clsm, n.var = min(160, nrow(fit_rf_clsm$importance)), # min(num,)处的num为图中的菌属数量 # main = 'Top 19 - variable importance',type = 1) # 简约分类器(只取部分高影响的自变量) # 选择 top N 重要的 OTUs，例如上述已经根据“Mean Decrease Accuracy”排名获得 genus_select <- rownames(importance_genus)[1:160] # 数据子集的训练集和测试集 genus_train_top <- traindata[ ,c(genus_select, 'Class')] genus_test_top<- testdata[ ,c(genus_select, 'Class')] # set.seed(10001) # form_clsm1 <- as.formula( # paste0( # "class~", # paste(colnames(genus_train_top)[1:10],collapse = "+") # ) # ) # 构建模型 # fit_rf_clsm1 <- randomForest( # form_clsm1, # data = genus_train_top, # ntree = 500, # mtry = 6, # importance = T # ) x_train1 <- genus_train_top[, -160 - 1] # 自变量 x y_train1 <- as.factor(genus_train_top[, 160 + 1]) # 因变量 y # fit_rf_clsm_test1 <- randomForest(x = x_train1, # y = y_train1, # ntree = 500, # 增加树的数量以提高稳定性 # importance = TRUE # 启用特征重要性计算 # ) # # fit_rf_clsm_test1 # # # 5 折交叉验证，重复 10 次 # ctrl1 <- trainControl(method = "repeatedcv", number = 5, repeats = 10) # # # 定义 mtry 和 ntree 的参数范围 # grid3 <- expand.grid(mtry = c(2:15)) # 定义mtry范围 # # # 进阶 # # num_features <- ncol(X_train) # # mtry_range <- floor(sqrt(num_features)) # 取平方根并向下取整 # # grid_mtry <- expand.grid(mtry = seq(max(2, mtry_range - 5), # # min(mtry_range + 5, num_features))) # # # 使用caret包进行调参 # fit_rf_clsm_test2 <- train(x = x_train1, y = y_train1, # 自变量 x， 因变量 y # method = "rf", # 使用的模型方法为随机森林（rf） # trControl = ctrl1, # tuneGrid = grid3,) # # 输出结果 # print(fit_rf_clsm_test2) # # # 定义最佳mtry参数 # grid4 <- expand.grid(mtry = c(16)) # # # 定义模型列表，存储每一个模型评估结果 # modellist1 <- list() # # 调整的参数是决策树的数量 # for (ntree in seq(100, 2000, by=100)) { # seq函数构建一个集合，间距by = 100 # set.seed(100003) # fit1 <- train(x = x_train1, y = y_train1, method="rf", # metric="Accuracy", tuneGrid=grid4, # trControl=ctrl1, ntree=ntree) # key1 <- toString(ntree) # modellist1[[key1]] <- fit1 # } # # # compare results # results1 <- resamples(modellist1) # # # 输出最佳模型和参数 # summary(results1) # # # 可视化模型性能，accuracy越近进1越好，kappa越近1越好 # bwplot(results1) # 箱线图 # dotplot(results1) # 点图 # densityplot(results1) # 密度图 # 使用最佳参数训练模型 set.seed(1) fit_rf_clsm1 <- randomForest(x = x_train1, y = y_train1, mtry = 12, ntree = 500, importance = T) print(fit_rf_clsm1) # # "ERROR & TREES" # plot(fit_rf_clsm1,main = "ERROR & TREES") # # # plot(randomForest::margin(fit_rf_clsm1), main = '观测值被判断正确的概率图') # # 预测 # # 训练集预测概率 # trainpredprob <- predict(fit_rf_clsm1,newdata = genus_train_top,type = "prob") # # 训练集ROC # multiclass.roc(response = genus_train_top$Class,predictor = trainpredprob) # # 训练集预测分类 #trainpredlab <- predict(fit_rf_clsm1,newdata = genus_train_top,type = "Class") # # 训练集混淆矩阵 # confusionMatrix(data = trainpredlab, # reference = genus_train_top$Class, # mode = "everything") # 预测 # 测试集预测概率 # x_test1 <- genus_test_top[, -11] # 自变量 x # y_test1 <- as.factor(genus_test_top[, 11]) # 因变量 y testpredprob <- predict(fit_rf_clsm1,newdata = genus_test_top,type = "prob") write.table(testpredprob, file = "D:/Rpackages/testpredprob.txt", sep = "\t", row.names = FALSE, col.names = TRUE) # 测试集ROC multiclass.roc(response = genus_test_top$Class,predictor = testpredprob) # 测试集预测分类 testpredlab <- predict(fit_rf_clsm1,newdata = genus_test_top,type = "Class") # 测试集混淆矩阵 confusion_matrix <- confusionMatrix(data = testpredlab, reference = genus_test_top$Class, mode = "everything") # 测试集综合结果 multiClassSummary( data.frame(obs = genus_test_top$Class,pred=testpredlab), lev = levels(genus_test_top$Class) ) # 查看样本预测结果 results <- data.frame(Actual = genus_test_top$Class, Predicted = testpredlab) # 测试集预测分类 # testpredlab <- predict(fit_rf_clsm1,newdata = testdata,type = "class") # t <- table(testpredlab,testdata$class) # acc = sum(diag(t))/nrow(testdata)*100 # print(paste("模型准确率为：",round(acc,4),sep='')) # 绘制混淆矩阵热图(内容复杂，好汉谨慎处之) # confusion_matrix是混淆矩阵对象 # 转换混淆矩阵为数据框 roc1 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,1:6]) roc1 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,1]) plot(roc1$rocs[[1]],col="#1f77b4",print.auc = TRUE,print.auc.x=0.8,print.auc.y=0.8) roc2 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,2]) plot(roc2$rocs[[1]],add=TRUE,col="#ff7f0e",print.auc = TRUE,print.auc.x=0.6,print.auc.y=0.6) roc3 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,3]) plot(roc3$rocs[[1]],add=TRUE,col="#2ca02c",print.auc=TRUE,print.auc.x=0.5,print.auc.y=0.5) roc4 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,4]) plot(roc4$rocs[[1]],add=TRUE,col="#d62728",print.auc=TRUE,print.auc.x=0.4,print.auc.y=0.4) roc5 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,5]) plot(roc5$rocs[[1]],add=TRUE,col="#9467bd",print.auc=TRUE,print.auc.x=0.3,print.auc.y=0.3) roc6 <- multiclass.roc(response = genus_test_top$Class,predictor =testpredprob[,6]) plot(roc1$rocs[[6]], add = TRUE, col = "#8c564b", print.auc = TRUE, print.auc.x = 0.2, print.auc.y = 0.2) # confusion_matrix_df <- as.data.frame.matrix(confusion_matrix$table) # colnames(confusion_matrix_df) <- c("F","H") # rownames(confusion_matrix_df) <- c("F","H") # # #c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15") # colnames(confusion_matrix_df) <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14", "15") # rownames(confusion_matrix_df) <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14", "15") # # 计算归一化值 # draw_data <- round(confusion_matrix_df / colSums(confusion_matrix_df),2) # draw_data <- as.data.frame(t(draw_data)) #使用t()函数转换draw_data的行列，原本列为真实类别，现转为行 # draw_data$real <- rownames(draw_data) #提取行名作为新生成列real的类别 # draw_data <- melt(draw_data) #将宽格式的数据框转换为长格式 # # # 绘制矩阵热图 # confusion1 <- 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 = "Prediction", y = "Reference", title = "Confusion matrix") + # theme_prism(border = T) + # theme(panel.border = element_blank(), # axis.ticks.y = element_blank(), # axis.ticks.x = element_blank(), # legend.position="right", # plot.title = element_text(hjust = 0.5)) # # confusion1 # # 二分类ROC曲线绘制 # # 计算ROC曲线的参数 # roc_obj <- roc(response = genus_test_top$class, predictor = testpredprob[, 2]) # roc_auc <- auc(roc_obj) # # 将ROC对象转换为数据框 # roc_data <- data.frame(1 - roc_obj$specificities, roc_obj$sensitivities) # # 绘制ROC曲线 # ROC1 <- ggplot(roc_data, aes(x = 1 - roc_obj$specificities, y = roc_obj$sensitivities)) + # geom_line(color = "#0073C2FF", size = 1) + # annotate("segment", x = 0, y = 0, xend = 1, yend = 1, linetype = "dashed", color = "gray") + # annotate("text", x = 0.8, y = 0.2, label = paste("AUC =", round(roc_auc, 3)), 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) # ROC1 # # # 计算 ROC 和 AUC # roc_obj <- roc(response = genus_test_top$class, predictor = testpredprob[, 2]) # roc_auc <- auc(roc_obj) # # # 将 ROC 对象转换为数据框 # roc_data <- data.frame( # FPR = 1 - roc_obj$specificities, # TPR = roc_obj$sensitivities # ) # # # 平滑处理 # smooth_roc <- data.frame( # FPR = spline(roc_data$FPR, n = 500)$x, # TPR = spline(roc_data$TPR, n = 500)$y # ) # # # 绘制平滑后的 ROC 曲线 # ROC1 <- ggplot(smooth_roc, aes(x = FPR, y = TPR)) + # geom_line(color = "#0073C2FF", size = 1) + # annotate("segment", x = 0, y = 0, xend = 1, yend = 1, linetype = "dashed", color = "gray") + # annotate("text", 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("Smoothed ROC Curve") + # theme(plot.title = element_text(size = 14, face = "bold")) + # theme_prism(border = T) # # ROC1 # # geom_smooth(se = FALSE, method = 'glm',formula = y~ns(x,6)) # # # 保存模型 # saveRDS(fit_rf_clsm, # file = "D:/Documents/R.data/fit_rf_clsm1_UnSimplifiedSiteClass_ConcernGender_UnexposedTop19_240209.rds") # # # 读取模型 # fit_rf_clsm1 <- readRDS("D:/Documents/R.data/fit_rf_clsm1_hand_Simplified_240102.rds") # # # 读取待分类数据 # testdata1 <- read.csv(file.choose()) # testdata1 <- as.data.frame(testdata) # for (i in c(1:2)) { # testdata1[,i] <- factor(testdata1[,i]) # } # # # 应用 # # 待分类数据集预测概率 # testpredprob <- predict(fit_rf_clsm1, newdata = testdata1, type = "prob") # # 测试集ROC # multiclass.roc(response = testdata1$class, predictor = testpredprob) # # 待分类数据集预测分类 # testpredlab <- predict(fit_rf_clsm1, newdata = testdata1,type = "class") # # 待分类数据集混淆矩阵 # confusion_matrix <- confusionMatrix(data = testpredlab, # reference = testdata1$class, # mode = "everything") # # 待分类数据集综合结果 # multiClassSummary( # data.frame(obs = testdata1$class,pred=testpredlab), # lev = levels(testdata1$class) # ) # # # 查看样本分类结果 # results <- data.frame(Actual = testdata1$class, Predicted = testpredlab) （这是整个的代码，testpredprob <- predict(fit_rf_clsm1,newdata = genus_test_top[1:5, ],type = "prob")，跑不下来）