随机森林回归模型

1. 分为测试集和验证集 

library(randomForest)
data("mtcars")
data=mtcars
set.seed(123)
train <- sample(nrow(data), nrow(data)*0.7)
data_train <- data[train, ]
data_test <- data[-train, ]

2. 初步建立随机森林模型

set.seed(123)
data_train.forest <- randomForest(mpg~., data = data_train, importance = TRUE)
data_train.forest

Call:
 randomForest(formula = mpg ~ ., data = data_train, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 3

          Mean of squared residuals: 4.936502
                    % Var explained: 84.44

结果中，% Var explained体现了预测变量对响应变量有关方差的整体解释率。

查看该模型的预测性能，可以看到具有较高的精度。

#使用训练集，查看预测精度
data_predict <- predict(data_train.forest,data_train)

plot(data_train$mpg, data_predict, main = '训练集',
    xlab = 'x', ylab = 'Predict')
abline(1, 1)
 
#使用测试集，评估预测性能
data_predict <- predict(data_train.forest, data_test)
 
plot(data_test$mpg, data_predict, main = '测试集',
    xlab = 'x', ylab = 'Predict')
abline(1, 1)

 

但是，并非这所有的变量都对回归的精度具有可观的贡献。有些变量可能在回归中产生较大噪声，对模型精度带来较高的误差。因此，最好将低贡献的变量去除。

3. 评估预测变量的重要性

基于已经构建好的随机森林回归模型，可以从中评估变量的重要性。将变量按重要程度高低进行排名后，选择排名靠前的一部分变量，这些重要的变量是明显的与植物生长时间密切关联的一些细菌类群。

 

#查看表示每个预测变量重要性的得分
#summary(data_train.forest)
importance_data <- data_train.forest$importance
head(importance_data)

#或者使用函数 importance()
importance_data <- data.frame(importance(data_train.forest), check.names = FALSE)
head(importance_data)
        %IncMSE IncNodePurity
cyl   6.3586844      96.76923
disp  6.8292862     134.34787
hp    7.1451720     128.48986
drat  0.9911803      38.03542
wt   11.3510772     162.70969
qsec  1.0467295      19.05801

#作图展示 top30 
varImpPlot(data_train.forest, n.var = min(30, nrow(data_train.forest$importance)),
    main = 'Top 30 - variable importance')

“%IncMSE”即increase in mean squared error，通过对每一个预测变量随机赋值，如果该预测变量更为重要，那么其值被随机替换后模型预测的误差会增大。因此，该值越大表示该变量的重要性越大；

“IncNodePurity”即increase in node purity，通过残差平方和来度量，代表了每个变量对分类树每个节点上观测值的异质性的影响，从而比较变量的重要性。该值越大表示该变量的重要性越大。

对于“%IncMSE”或“IncNodePurity”，二选一作为判断预测变量重要性的指标。需注意的是，二者的排名存在一定的差异。

#可以根据某种重要性的高低排个序，例如根据“IncNodePurity”指标
importance_data <- importance_data[order(importance_data$IncNodePurity, decreasing = TRUE), ]
head(importance_data)

4. 交叉验证

最终选择的变量：可通过执行十折交叉验证，根据交叉验证曲线对变量进行取舍。交叉验证法的作用就是尝试利用不同的训练集/验证集划分来对模型做多组不同的训练/验证，来应对单独测试结果过于片面以及训练数据不足的问题。此处使用训练集本身进行交叉验证。

##交叉验证辅助评估选择特定数量的 variable
#5 次重复十折交叉验证
set.seed(123)
data_train.cv <- replicate(5, rfcv(data_train[-ncol(data_train)], data_train$mpg, cv.fold = 10, step = 1.5), simplify = FALSE)
data_train.cv

#提取验证结果绘图
data_train.cv <- data.frame(sapply(data_train.cv, '[[', 'error.cv'))
data_train.cv$id <- rownames(data_train.cv)
data_train.cv <- reshape2::melt(data_train.cv, id = 'id')
data_train.cv$id <- as.numeric(as.character(data_train.cv$id))
 
data_train.cv.mean <- aggregate(data_train.cv$value, by = list(data_train.cv$id), FUN = mean)
head(data_train.cv.mean, 10)
 
#拟合线图
library(ggplot2)
 
ggplot(data_train.cv.mean, aes(Group.1, x)) +
geom_line() +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +
labs(title = '',x = 'Number of variabls', y = 'Cross-validation error')

选择低点，示例选4个

#提示保留 4个重要的 OTU，可以使随机森林回归的精度最大化
#首先根据某种重要性的高低排个序，例如根据“IncNodePurity”指标
importance_data<- importance_data[order(importance_data$IncNodePurity, decreasing = TRUE), ]

#然后取出排名靠前的变量，例如 top4 最重要的 OTU
importance_data.select <- importance_data[1:10, ]
importance_data.select

#不妨简单查看下这些重要的 变量与observation的关系
data_id.select <- rownames(importance_data.select)
data.select <- data[ ,c(data_id.select, 'mpg')]
data.select <- reshape2::melt(data.select, id = 'mpg')
 
ggplot(data.select, aes(x = mpg, y = value)) +
geom_point() +
geom_smooth() +
facet_wrap(~variable, ncol = 3, scale = 'free_y') +
labs(title = '',x = 'mpg', y = 'Relative abundance')

 5. 最终模型的确定

##只包含 4个重要预测变量的简约回归
data.select <- data[ ,c(data_id.select, 'mpg')]
 
#为了方便后续评估随机森林模型的性能，将总数据集分为训练集（占 70%）和测试集（占 30%）
set.seed(123)
train <- sample(nrow(data.select), nrow(data.select)*0.7)
data_train.select <- data.select[train, ]
data_test.select <- data.select[-train, ]
 
#随机森林计算（默认生成 500 棵决策树），详情 ?randomForest
set.seed(123)
data_train.select.forest <- randomForest(mpg~., data = data_train.select, importance = TRUE)
data_train.select.forest
 
#使用训练集，查看预测精度
data_predict <- predict(data_train.select.forest, data_train.select)
 
plot(data_train.select$mpg, data_predict, main = '训练集',
    xlab = 'mpg', ylab = 'Predict')
abline(1, 1)
 
#使用测试集，评估预测性能
data_predict <- predict(data_train.select.forest, data_test.select)
 
plot(data_test.select$mpg,data_predict, main = '测试集',
    xlab = 'mpg', ylab = 'Predict')
abline(1, 1)

如果不涉及预测，只是评估变量的重要性

上文过程中，为了评估随机森林回归模型的预测性能，将数据拆分为训练集和测试集分开讨论。但如果目的不涉及预测，只是为了寻找重要的预测变量，完全可以直接使用所有数据代入回归，无需再区分训练集与测试集，这也是很多文献中经常提到的方法。