我们使用包randomForest并利用鸢尾花数据建立一个预测模型。包里面的randomForest()函数有两点不足:第一,它不能处理缺失值,使得用户必须在使用该函数之前填补这些缺失值;第二,每个分类属性的最大数量不能超过32个,如果属性超过32个,那么在使用randomForest()之前那些属性必须被转化。 也可以通过另外一个包’cforest’建立随机森林,并且这个包里面的函数并不受属性的最大数量约束,尽管如此,高维的分类属性会使得它在建立随机森林的时候消耗大量的内存和时间。

R中与决策树有关的Package:

单棵决策树:rpart/tree/C50
随机森林:  randomforest/ranger
梯度提升树:gbm/xgboost
树的可视化:rpart.plot

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,]
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)   

randomForest小结

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()查看结果。图中数据点的边距为正确分类的比例减去被归到其他类别的最大比例。一般来说,边距为正数说明该数据点划分正确。