参考资料:
随机森林(2)R实战][https://www.yuque.com/biotrainee/biostat/chapter3-24]
优雅R 微信公众号:https://mp.weixin.qq.com/s/fHDxqOOTW2CbKbSDVOD0bg
关于决策树及随机森林模型的理论基础可以参考优雅R 微信公众号文章:https://mp.weixin.qq.com/s/fHDxqOOTW2CbKbSDVOD0bg
1、数据
> data1
success_rate body.mass..g..x Migration.distance..km.. Mass...average...g Length...average...mm
1 0.0000000 30.0 6052.5 36.00 195.0
2 0.0000000 8.6 19.8 8.00 140.0
3 0.9453125 1141.0 1678.2 1082.00 NA
4 0.0000000 54.5 1238.7 57.75 181.0
5 0.9981481 354.5 0.0 358.70 NA
...
success_rate:繁殖的成功率
body.mass:体重
Migration.distance:年迁徙距离
Length:体长
其中包含了NA缺失值
建立模型的意图为:1、计算各变量对模型的贡献度。2、用训练好的模型预测测试集。
2、缺失值处理
(2.1)删除缺失值(2.1或2.2都可,取决于具体样本情况)
> data1<-na.omit(data1)
(2.2)利用随机森林模型预测补全NA值
install.packages("randomForest")
library(randomForest)
data.imputed <- rfImpute(success_rate ~ ., data = data1, iter=10,ntree=500)
使用500棵决策树组成的随机森林进行补全,重复10次。用success_rate来拟合其他所有缺失值。“success_rate ~ .”中的“.”意为所有参数。
3、建立模型
(3.1)广义线性模型(glm()函数)(3.1或3.2择一)
glm_model<-glm(success_rate~ ., data=data.imputed,family=binomial(link = "logit"))
glm()函数是R自带的函数。“success_rate~ .”意为用其他参数作为变量因素来拟合success_rate。family参数为模型的选择,默认为正态分布(高加索分布),其他的选择还有:
binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")
结果查看:
> summary(glm_model)
Call:
glm(formula = success_rate ~ ., family = binomial(link = "logit"),
data = data.imputed)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5494 -0.9805 -0.4454 0.9659 1.6085
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7536342 0.6148484 -1.226 0.2203
body.mass..g..x 0.0001538 0.0004608 0.334 0.7385
Migration.distance..km.. -0.0003455 0.0001704 -2.028 0.0426 *
Mass...average...g -0.0002686 0.0004168 -0.644 0.5193
Length...average...mm 0.0026700 0.0018083 1.477 0.1398
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 97.757 on 89 degrees of freedom
Residual deviance: 88.492 on 85 degrees of freedom
AIC: 123.9
Number of Fisher Scoring iterations: 4
Estimate表示模型中各参数的系数,正负值即为相关性的正负,绝对值大小代表对模型的贡献程度。Pr(>|z|) 显示了因素对模型的影响是否显著,一般来说应该保留显著的,去掉不显著的因素,重新跑模型(这里只Migration.distance有一个显著的,就不继续做了)。
(3.2)随机森林模型
library(randomForest)
model_rf <- randomForest(success_rate ~ ., data=data.imputed,ntree=600)
model_rf 为训练好的模型
plot(model_rf)
绘制标准差变化图,可以看到随着决策树计算的个数增多,标准差下降到稳定值。
查看因素重要性并绘图
> model_rf$importance
IncNodePurity
body.mass..g..x 3.868491
Migration.distance..km.. 3.640498
Mass...average...g 4.019628
Length...average...mm 3.728895
#IncNodePurity越高,表明Gini值越低,因素贡献率越高
> varImpPlot(model_rf)
用模型预测测试集:
假设有一个预测集testset,数据结构与data1相同,通过model_rf预测相关参数:
model.prediction <- predict(model.rf, testset)
此外,随机森林还可以用于利用SNPs预测杂交categories等等。需要注意,如果是非连续变量,需要转化为因子进行分析。