我们使用包randomForest并利用鸢尾花数据建立一个预测模型。包里面的randomForest()函数有两点不足:第一,它不能处理缺失值,使得用户必须在使用该函数之前填补这些缺失值;第二,每个分类属性的最大数量不能超过32个,如果属性超过32个,那么在使用randomForest()之前那些属性必须被转化。 也可以通过另外一个包’cforest’建立随机森林,并且这个包里面的函数并不受属性的最大数量约束,尽管如此,高维的分类属性会使得它在建立随机森林的时候消耗大量的内存和时间。
R中与决策树有关的Package:
单棵决策树:rpart/tree/C50
随机森林: randomforest/ranger
梯度提升树:gbm/xgboost
树的可视化:rpart.plot
data(iris)
ind <- sample(1:2, nrow(iris), replace=TRUE, prob=c(0.7, 0.3))
trainData <- iris[ind==1,]
testData <- iris[ind==2,]
str(trainData)
## 'data.frame': 105 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.7 5 5.4 4.9 5.4 4.3 5.8 5.7 5.4 ...
## $ Sepal.Width : num 3.5 3.2 3.6 3.9 3.1 3.7 3 4 4.4 3.9 ...
## $ Petal.Length: num 1.4 1.3 1.4 1.7 1.5 1.5 1.1 1.2 1.5 1.3 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.4 0.1 0.2 0.1 0.2 0.4 0.4 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
str(testData)
## 'data.frame': 45 obs. of 5 variables:
## $ Sepal.Length: num 4.9 4.6 4.6 5 4.4 4.8 4.8 5.1 5.4 4.6 ...
## $ Sepal.Width : num 3 3.1 3.4 3.4 2.9 3.4 3 3.8 3.4 3.6 ...
## $ Petal.Length: num 1.4 1.5 1.4 1.5 1.4 1.6 1.4 1.5 1.7 1 ...
## $ Petal.Width : num 0.2 0.2 0.3 0.2 0.2 0.2 0.1 0.3 0.2 0.2 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
选取randomforest – mtry节点值,对应误差最小为2,一般可默认。通常也是2记得。mtry指定节点中用于二叉树的变量个数,默认情况下数据集变量个数的二次方根(分类模型)或三分之一(预测模型)。
library(randomForest)
n <- length(names(trainData))
set.seed(100)
for(i in 1:(n-1)){
mtry_fit <- randomForest(Species~., data=trainData, mtry=i)
err <- mean(mtry_fit$err.rate)
print(err)
}
## [1] 0.05697697
## [1] 0.05508195
## [1] 0.05608843
## [1] 0.05356275
之后选择ntree值,ntree指定随机森林所包含的决策树数目,默认为500;.在400左右时,模型内误差基本稳定,故取ntree=400。
library(randomForest)
set.seed(100)
ntree_fit <- randomForest(Species~.,data=trainData, mtry=2, ntree=1000)
plot(ntree_fit)
由上图的结果可知,OOB误差为2.8%,同时在随机森林中,第二类和第三类仍然有误差,会被误判,也可以通过输入plot(rf)绘制每一棵树的误判率的图。
重要性查看
library(randomForest)
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
importance(ntree_fit)
## MeanDecreaseGini
## Sepal.Length 6.882673
## Sepal.Width 2.093474
## Petal.Length 31.676958
## Petal.Width 28.588881
importance(ntree_fit,type=1) # 重要性评分
##
## Sepal.Length
## Sepal.Width
## Petal.Length
## Petal.Width
importance(ntree_fit,type=2) # Gini指数
## MeanDecreaseGini
## Sepal.Length 6.882673
## Sepal.Width 2.093474
## Petal.Length 31.676958
## Petal.Width 28.588881
varImpPlot(ntree_fit) # 可视化
利用iris数据,可以看到这四个变量的重要性排序基本是一样的。
模型的预测功能:predict中有多种参数,比如Nodes,Proximity,predict.all。
predict(object, newdata, type="response",
norm.votes=TRUE, predict.all=FALSE, proximity=FALSE, nodes=FALSE,
cutoff, ...)
Nodes判断是否是终点。Proximity判断是否需要进行近邻测量。predict.all判断是否保留所有的预测器。
举例,以前面的随机森林模型进行建模。
predict.all会输出一个150*150的字符矩阵,代表每一颗树的150个预测值(前面预设了ntree=100); Nodes输出100颗树的节点情况。
prediction <- predict(Randommodel, data[,1:5],type="class") # 还有response回归类型
table(observed =data$Species,predicted=prediction)
library(randomForest)
data(iris)
ind <- sample(1:2, nrow(iris), replace=TRUE, prob=c(0.7, 0.3))
trainData <- iris[ind==1,]
testData <- iris[ind==2,]
train_data = trainData
test_data = trainData[,-5]
numtree = 300
# train_data = trainData[,1:4]
# test_data = trainData[,-c(4:5)]
# numtree = 300
set.seed(100)
formula <- as.formula(paste(names(train_data)[ncol(train_data)], " ~ ", paste(names(train_data)[-ncol(train_data)], collapse= "+")))
train_rf <- randomForest(formula, data=train_data, ntree=numtree, proximity=TRUE, importance=TRUE)
result <- predict(train_rf, newdata=test_data)
# result
print(train_rf)
##
## Call:
## randomForest(formula = formula, data = train_data, ntree = numtree, proximity = TRUE, importance = TRUE)
## Type of random forest: classification
## Number of trees: 300
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 6.67%
## Confusion matrix:
## setosa versicolor virginica class.error
## setosa 33 0 0 0.00000000
## versicolor 0 33 3 0.08333333
## virginica 0 4 32 0.11111111
attributes(train_rf)
## $names
## [1] "call" "type" "predicted"
## [4] "err.rate" "confusion" "votes"
## [7] "oob.times" "classes" "importance"
## [10] "importanceSD" "localImportance" "proximity"
## [13] "ntree" "mtry" "forest"
## [16] "y" "test" "inbag"
## [19] "terms"
##
## $class
## [1] "randomForest.formula" "randomForest"
# rf[1:length(train_rf)]
plot(train_rf)
importance(train_rf)
## setosa versicolor virginica MeanDecreaseAccuracy
## Sepal.Length 5.363700 2.802764 4.7755814 6.794226
## Sepal.Width 3.755962 2.384586 -0.5129371 3.436602
## Petal.Length 17.231612 20.527279 16.2581403 22.687567
## Petal.Width 18.162923 25.496244 21.7707477 27.602349
## MeanDecreaseGini
## Sepal.Length 6.205470
## Sepal.Width 3.105701
## Petal.Length 28.410449
## Petal.Width 31.523618
# importance(train_rf,type=1)
# importance(train_rf,type=2)
varImpPlot(train_rf)
irisPred <- predict(train_rf, testData)
table(irisPred, testData$Species)
##
## irisPred setosa versicolor virginica
## setosa 17 0 0
## versicolor 0 13 0
## virginica 0 1 14
# table(irisPred, testData[,ncol(testData)])
# table(type=testData[,ncol(testData)], predict=irisPred)
# 绘制每一个观测值被判断正确的概率图
plot(margin(train_rf, testData$Species),main="观测值被判断正确的概率图")
fun_randomForest <- function(train_data, test_data, numtree){
set.seed(100)
formula <- as.formula(paste(names(train_data)[ncol(train_data)], " ~ ", paste(names(train_data)[-ncol(train_data)], collapse= "+")))
train_rf <- randomForest(formula, data=train_data, ntree=numtree, proximity=TRUE, importance=TRUE)
result <- predict(train_rf, newdata=test_data)
result
}
fun_randomForest(train_data, test_data, numtree)
## 1 2 3 5 7 8
## setosa setosa setosa setosa setosa setosa
## 10 11 13 15 19 21
## setosa setosa setosa setosa setosa setosa
## 22 23 24 25 27 28
## setosa setosa setosa setosa setosa setosa
## 29 30 31 32 34 36
## setosa setosa setosa setosa setosa setosa
## 37 38 40 43 44 45
## setosa setosa setosa setosa setosa setosa
## 46 48 49 54 55 57
## setosa setosa setosa versicolor versicolor versicolor
## 58 59 61 62 63 65
## versicolor versicolor versicolor versicolor versicolor versicolor
## 66 67 68 70 73 75
## versicolor versicolor versicolor versicolor versicolor versicolor
## 76 77 78 79 80 81
## versicolor versicolor versicolor versicolor versicolor versicolor
## 83 84 85 87 88 89
## versicolor versicolor versicolor versicolor versicolor versicolor
## 90 92 93 94 95 96
## versicolor versicolor versicolor versicolor versicolor versicolor
## 97 98 99 101 103 104
## versicolor versicolor versicolor virginica virginica virginica
## 106 107 108 109 110 111
## virginica virginica virginica virginica virginica virginica
## 112 114 115 116 117 119
## virginica virginica virginica virginica virginica virginica
## 120 122 123 124 125 127
## virginica virginica virginica virginica virginica virginica
## 129 131 132 134 135 138
## virginica virginica virginica virginica virginica virginica
## 140 141 142 144 145 147
## virginica virginica virginica virginica virginica virginica
## 148 149 150
## virginica virginica virginica
## Levels: setosa versicolor virginica
代码解读:
randomForset,执行建模,x参数设定自变量数据集,y参数设定因变量数据列,importance设定是否输出因变量在模型中的重要性,如果移除某个变量,模型方差增加的比例是它判断变量重要性的标准之一,proximity参数用于设定是否计算模型的临近矩阵,ntree用于设定随机森林的树数,最后一句输出模型在训练集上的效果。
print:输出模型在训练集上的效果,可以看出错误率为3.33%,维持在比较低的水平。
mtry:选择的分裂属性的个数
proximity=TRUE:表示生成临近矩阵
importance=TRUE:输出分裂属性的重要性
最后使用测试集进行测试,用table()和margin()查看结果。图中数据点的边距为正确分类的比例减去被归到其他类别的最大比例。一般来说,边距为正数说明该数据点划分正确。