「R」逻辑回归、决策树、随机森林

有监督学习基于一组包含预测变量和输出变量的样本单元。将全部数据分为一个训练数据集和一个验证数据集，其中训练集用于建立预测模型，验证集用于测试模型的准确性。

这部分通过 rpart 、 rpart.plot 和 party 包来实现决策树模型及其可视化，通过 randomForest 包拟合随机森林，通过 e1071 包构造支持向量机，通过R中的基本函数 glm() 实现逻辑回归。在探索之前，先安装好相应的包。

pkgs <- c("rpart", "rpart.plot", "party",
          "randomForest", "e1071")
install.packages(pkgs, dependencies = TRUE)

这里的例子来源于UCI机器学习数据库中的威斯康星州乳腺癌数据。数据分析的目的是根据细胞组织细针抽吸活检所反应的特征，来判断被捡者是否患有乳腺癌。

数据准备

该数据集是逗号分隔的txt文件，包含699个样本蛋白，其中458个良性，241个为恶性。数据集中有11个变量，表中未标明变量名。其中16个样本单元中有缺失数据并用问号（?）表示。

变量为

• ID
• 肿块厚度
• 细胞大小的均匀性
• 细胞形状的均匀性
• 边际附着力
• 单个上皮细胞的大小
• 裸核
• 乏味染色体
• 正常核
• 有丝分裂
• 类别

ID不纳入数据分析，最后一个变量是输出变量（良性=2，恶性=4）。

对于每一个样本来说，另外九个变量是与判别恶性肿瘤相关的细胞特征，并加以记录。这些细胞特征得分为1（最接近良性）到10（最接近病变）之间的整数。任一变量都不能单独作为判别良性或恶性的标准，建模的目的是找到九个细胞特征的某种组合，从而实现对恶性肿瘤的准确预测。

## Data preparation
loc <- "http://archive.ics.uci.edu/ml/machine-learning-databases/"
ds  <- "breast-cancer-wisconsin/breast-cancer-wisconsin.data" 
url <- paste(loc, ds, sep="")
breast <- read.table(url, sep = ",", header = FALSE, na.strings = "?")
names(breast) <- c("ID", "clumpThickness", "sizeUniformity", "shapeUniformity",
                   "maginalAdhesion", "singleEpithelialCellSize", "bareNuclei",
                   "blandChromatin", "normalNucleoli", "mitosis", "class")
df <- breast[-1]
df$class <- factor(df$class, levels = c(2,4),
                   labels = c("benign", "malignant"))
set.seed(1234)
train <- sample(nrow(df), 0.7*nrow(df))
df.train <- df[train,]
df.validate <- df[-train,]
table(df.train$class)
table(df.validate$class)

逻辑回归

逻辑回归 是广义线性模型的一种，它根据一组数值变量预测二元输出（之前在广义模型中有介绍）。基本函数 glm() 可以用于拟合逻辑回归模型。

## logistic regression
fit.logit <- glm(class ~ ., data=df.train, family=binomial()) # 拟合逻辑回归
summary(fit.logit)                            #　检查模型
Call:
glm(formula = class ~ ., family = binomial(), data = df.train)
Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.75813  -0.10602  -0.05679   0.01237   2.64317  
Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -10.42758    1.47602  -7.065 1.61e-12 ***
clumpThickness             0.52434    0.15950   3.287  0.00101 ** 
sizeUniformity            -0.04805    0.25706  -0.187  0.85171    
shapeUniformity            0.42309    0.26775   1.580  0.11407    
maginalAdhesion            0.29245    0.14690   1.991  0.04650 *  
singleEpithelialCellSize   0.11053    0.17980   0.615  0.53871    
bareNuclei                 0.33570    0.10715   3.133  0.00173 ** 
blandChromatin             0.42353    0.20673   2.049  0.04049 *  
normalNucleoli             0.28888    0.13995   2.064  0.03900 *  
mitosis                    0.69057    0.39829   1.734  0.08295 .  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
prob <- predict(fit.logit, df.validate, type="response")　　#　建立预测模型,对训练集集外样本单元分类
logit.pred <- factor(prob > .5, levels=c(FALSE, TRUE),
                     labels = c("benign", "malignant"))　　
logit.perf <- table(df.validate$class, logit.pred,
                    dnn=c("Actual", "Predicted"))　　　　　　# 评估预测准确性
logit.perf
           Predicted
Actual      benign malignant
  benign       118         2
  malignant      4        76

predict() 函数默认输出肿瘤为恶性的对数概率，指定参数 type="response" 即可以得到预测肿瘤为恶性的概率。在样本单元中，概率大于.5的被分为恶性肿瘤类，概率小于等于.5的被分为良性肿瘤类。

最后给出预测与实际情况对比的交叉表（混淆矩阵，confusion matrix）。数据集中有１０个单元包含缺失数据而无法判别。

在验证集上，正确分类的模型（准确率，accuracy）为(76+118)/200=97%。看起来还是非常准确的哈～

值得注意的是，模型中有三个预测变量的系数未能通过显著性检验，一般而言可以将它们去除从而精简模型。

当然，可以用逐步逻辑回归生成一个包含更少解释变量的模型，其目的是通过增加或移除变量来得到一个更小的AIC值。例如本例可以用

logit.fit.reduced <- step(fit.logit) 来得到一个更为精简的模型。

决策树

决策树是数据挖掘领域中常用模型。其 基本思想 是对预测变量进行二元分离，从而构造一颗可以预测新样本单元所属类别的树。这里介绍两类决策树：经典树和条件推断树。

这里讲的基本思想有点精悍，似懂非懂的样子。我们具体来看看它们究竟是什么吧。

经典决策树

经典决策树以一个二元输出变量和一组预测变量为基础。其具体算法如下：

1. 选定一个最佳预测变量将全部样本单元分为两类，实现两类中的纯度最大化（即一类中良性样本单元尽可能多，另一类恶性样本单元尽可能多）。如果预测变量连续，则选定一个分割点进行分类，使得两类纯度最大化；如果预测变量为分类变量，则对各类别进行合并再分类。
2. 对每个子类别继续执行步骤1。
3. 重复步骤1~2，直到子类别中所含的样本单元树过少，或者没有分类能将不纯度下降到一个给定阈值以下。最终集中的子类别即 终端节点 。根据每一个终端节点中样本单元的类别数众数来判别这一终端节点的所属类别。
4. 对任一样本单元执行决策树，得到其终端节点，即可以根据步骤3得到模型预测的所属类别。

这一过程就类似一棵树生长不断形成分支，这些分支的生成是依赖具体的算法要求（这里就是让它们纯度最大化）。

上述算法构建的树过大，容易出现过度拟合现象。可采用10折交叉验证法预测误差最小的树，然后用它进行预测。

R中的 rpart 包支持 rpart() 函数构造决策树， prune() 函数对决策树进行剪枝。下面给出针对数据集的算法实现。

library(rpart)
set.seed(1234)
dtree <- rpart(class ~ ., data=df.train, method="class",
              parms=list(split="information"))     # 生成树
dtree$cptable
        CP nsplit rel error xerror       xstd
1 0.800000      0   1.00000 1.0000 0.06484605
2 0.046875      1   0.20000 0.2750 0.03954867
3 0.012500      3   0.10625 0.1500 0.02985779
4 0.010000      4   0.09375 0.1625 0.03101007
plotcp(dtree)
dtree.pruned <- prune(dtree, cp=.0125)          # 剪枝
library(rpart.plot)
prp(dtree.pruned, type=2, extra=104,
    fallen.leaves = TRUE, main="Decision Tree")
dtree.pred <- predict(dtree.pruned, df.validate, type="class")   # 对训练集外样本单元分类
dtree.perf <- table(df.validate$class, dtree.pred,
                    dnn=c("Actual", "Predicted"))
dtree.perf
           Predicted
Actual      benign malignant
  benign       122         7
  malignant      2        79

dtree.png

rpart() 返回的 cptable 值中包括不同大小的树对应的预测误差，因此可以用于辅助设定最终的树的大小。其中，复杂度参数 cp 用于惩罚过大的树；树的大小即分支数 nsplit ，有n个分支的树将有n+1个终端节点； rel error 栏即各种训练集中各种树对应的误差；交叉验证误差 xerror 即基于训练样本所得的10折交叉验证误差； xstd 栏为交叉验证误差的标准差。

借助 plotcp() 函数可画出交叉验证误差与复杂度参数的关系图（上图）。对于所有交叉验证误差在最小交叉验证误差一个标准差范围内的树，最小的树即最优的树。

由代码中的 cptable 表可以知道，四个终端节点（三次分割）的树满足要求。

deci_tree.png

在完整树的基础上， prune() 函数根据复杂度参数减掉最不重要的枝，从而将树的大小控制在理想范围内。从代码中的 cptable 内容中可以看到，三次分割对应的复杂度参数是0.0125，从而 prune(dtree, cp=0.0125) 可得到一个理想大小的树。

rpart.plo 包中的 prp() 函数可用于画出最终的决策树，它有很多的可供选择参数，如 type=2 可画出每个节点下分割的标签， extra=104 可画出每一类的概率以及每个节点处的样本占比， fallen.leaves=TRUE 可在图的底端显示终端节点。

对观测点分类时，从树的顶端开始，若满足条件则从左枝往下，否则右枝往下，重复这个过程知道碰到一个终端节点为止。该终端节点即为这一观测点的所属类别。

最后 predict() 函数用来对验证集中的观测点分类。代码内容中给出了实际类别与预测类别的交叉表。整体来看，准确率还是非常高的。

条件推断树

条件推断树与传统决策树类似，但变量和分割的选取是基于显著性检验的，而不是纯净度或同质性一类的度量。显著性检验是置换检验。

条件推断的算法如下：

1. 对输出变量与每个预测变量间的关系计算p值。
2. 选取p值最小的变量。
3. 在因变量与被选中的变量间尝试所有可能的二元分割（通过排列检验），并选取最显著的分割。
4. 将数据集分成两群，并对每个子群重复上述步骤。
5. 重复直至所有分割都不显著或已经达到最小节点为止。

条件推断树可由 party 包中的 ctree() 函数获得。

library(party)
fit.ctree <- ctree(class ~ ., data=df.train)
plot(fit.ctree, main="Conditional Inference Tree")
ctree.pred <- predict(fit.ctree, df.validate, type="response")
ctree.perf <- table(df.validate$class, ctree.pred, 
                    dnn=c("Actual", "Predicted"))
ctree.perf

con_infer_tree.png

值得注意的是，对于条件推断树来说，剪枝不是必需的，其生成过程相对更自动化一些。另外， party 包也提供了许多图像参数。

随机森林

随机森林 是一种组成式的有监督学习方法。在随机森林中，我们同时生成多个预测模型，并将模型的结果汇总以提升分类准确率。http://mng.bz/7Nul上有关于随机森林的详尽介绍。

随机森林的算法涉及对样本单元和变量的抽样，从而生成大量决策树。对每个样本单元来说，所有的决策树依次对其进行分类。所有决策树预测类别中的众数类别即为随机森林所预测的这一样本的类别。

假设训练集中共有N个样本单元，M个变量，则随机森林算法如下：

1. 从训练集中随机有放回地抽取N个样本单元，生成大量决策树。
2. 在每一个节点随机地抽取m<M个变量，将其作为分割节点的候选变量。每一个节点处的变量数应一致。
3. 完整生成所有决策树，无需剪枝。
4. 终端节点的所属类别由节点对应的众数类别决定。
5. 对于新的观测点，用所有的树对其进行分类，其类别由多数决定原则生成。

生成树时没有用到的样本点所对应的类别可以由生成的树估计，与其真实类别比较即可得到袋外预测（out-of-bag, OOB）误差。无法获得验证集时，这是随机森林的一大优势。随机森林算法可以计算变量的相对重要程度。

randomForest 包中的 randomForest() 函数可以用于生成随机森林。函数默认生成500棵树，并且默认在每个节点处抽取 sqrt(M) 个变量，最小节点为1。

library(randomForest)
set.seed(1234)
fit.forest <- randomForest(class ~ ., data=df.train,
                           na.action = na.roughfix,
                           importance = TRUE)   # 生成森林
fit.forest
importance(fit.forest, type=2)    # 给出变量重要性
forest.pred <- predict(fit.forest, df.validate)    # 对训练集外样本进行分类
forest.perf <- table(df.validate$class, forest.pred,
                     dnn = c("Actual", "Predicted"))
forest.perf

几个输出结果

> fit.forest
Call:
 randomForest(formula = class ~ ., data = df.train, importance = TRUE,      na.action = na.roughfix) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 3
        OOB estimate of  error rate: 3.68%
Confusion matrix:
          benign malignant class.error
benign       319        10  0.03039514
malignant      8       152  0.05000000
> importance(fit.forest, type=2)
                         MeanDecreaseGini
clumpThickness                  12.504484
sizeUniformity                  54.770143
shapeUniformity                 48.662325
maginalAdhesion                  5.969580
singleEpithelialCellSize        14.297239
bareNuclei                      34.017599
blandChromatin                  16.243253
normalNucleoli                  26.337646
mitosis                          1.814502