评估变量贡献度及模型预测——利用随机森林模型或广义线性模型
参考资料：
随机森林(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等等。需要注意，如果是非连续变量，需要转化为因子进行分析。
