当我们说机器学习的的时候，我们在说什么？

上图解释了完整的机器学习流程，包括构建任务、准备训练数据集及测试数据集、选择学习方法（leaner）、构建模型、预测、评估、重复执行并可视化学习情况。

下面所提及到的有监督机器学习各个方法只是Learner的一种方法而已。

有监督机器学习领域中包含许多可用于分类的方法，如逻辑回归、决策树、随机森林、支持向量机、神经网络等。

有监督学习基于一组包含预测变量值和输出变量值的样本单元。将全部数据分为一个训练集和一个验证集，其中训练集用于建立预测模型，验证集用于测试模型的准确性。对训练集和验证集的划分尤其重要，因为任何分类技术都会最大化给定数据的预测效果。用训练集建立模型并测试模型会使得模型的有效性被过分夸大，而用单独的验证集来测试基于训练集得到的模型则可使得估计更准确、更切合实际。得到一个有效的预测模型后，就可以预测那些只知道预测变量值的样本单元对应的输出值了。

通过 rpart 、 rpart.plot 和 party 包来实现决策树模型及其可视化，通过randomForest包拟合随机森林，通过 e1071包构造支持向量机，通过R中的基本函数 glm()实现逻辑回归。

机器学习的数据准备

下面给出R中数据准备流程。数据从UCI数据库中抽取，并随机分出训练集和验证集，其中训练集中包含499个样本单元（占70%），其中良性样本单元329个，恶性160个；验证集中包含210个样本单元（占30%），其中良性129个，恶性81个。

数据集中的ID变量不作考虑，类别变量为输出变量用作判断。

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)

predict函数对应一系列的函数用于预测模型，predict.glm, predict.lm, predict.loess, predict.nls, predict.poly, predict.princomp, predict.smooth.spline.
SafePrediction for prediction from (univariable) polynomial and spline fits.
For time-series prediction, predict.ar, predict.Arima, predict.arima0, predict.HoltWinters, predict.StructTS.

logistic回归模型分类

fit.logit <- glm(class~., data = df.train, family = binomial())
summary(fit.logit)
prob <- predict.glm(fit.logit, newdata = df.validate, type = "response") 
# 注意type= 参数
# type  
# the type of prediction required. The default is on the scale of the linear
# predictors; the alternative "response" is on the scale of the response variable.
# Thus for a default binomial model the default predictions are of log-odds
# (probabilities on logit scale) and type = "response" gives the predicted
# probabilities. The "terms" option returns a matrix giving the fitted values of
# each term in the model formula on the linear predictor scale.logit.prob <- factor(prob > 0.5, levels = c(FALSE, TRUE, labels = c("Benign", "Malignant"))) # 关注一下这里factor()进行差别并标注标签的妙用
table(df.validate$class, logit.prob, dnn = c("Actual", "Predicted")) # table()的交叉验证对象长度需要一致

解释：
1. 首先，以类别为响应变量，其余变量为预测变量。 2. 基于 df.train 数据框中的数据构造逻辑回归模型, 接着给出了模型中的系数。 3. 接着，采用基于 df.train 建立的模型来对 df.validate 数据集中的样本单元分类。predict() 函数默认输出肿瘤为恶性的对数概率，指定参数 type="response" 即可得到预测肿瘤为恶性的概率。样本单元中，概率大于0.5的被分为恶性肿瘤类，概率小于等于0.5的被分为良性肿瘤类。 4. 最后给出预测与实际情况对比的交叉表（即混淆矩阵，confusion matrix）。

df.validate 数据集中有10个样本单元因包含缺失数据而无法判别。在验证集上，正确分类的模型（即准确率，accuracy）为(76+118)/200=97%。下一步是验证模型的有效性。
同时要注意的是，模型中有三个预测变量（ sizeUniformity 、 shapeUniformity 和singleEpithelialCellSize ）的系数未通过显著性检验（即p值大于0.1）。从预测的角度来说，我们一般不会将这些变量纳入最终模型。当这类不包含相关信息的变量特别多时，可以直接将其认定为模型中的噪声。在这种情况下，可用逐步逻辑回归生成一个包含更少解释变量的模型，其目的是通过增加或移除变量来得到一个更小的AIC值。精简后的模型在验证集上的误差相对全变量模型更小。此处需要联系到之间学习的logistic回归模型如何减少纳入考虑的无效变量，参广义线性模型_高级回归分析.md。

决策树分类

基本思想是对预测变量进行二元分离，从而构造一棵可用于预测新样本单元所属类别的树。此处将介绍两类决策树：经典树和条件推断树。

经典决策树

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

不过，上述算法通常会得到一棵过大的树，从而出现过拟合现象。结果就是，对于训练集外单元的分类性能较差。为解决这一问题，可采用10折交叉验证法选择预测误差最小的树。这一剪枝后的树即可用于预测。
R中的 rpart 包支持 rpart() 函数构造决策树， prune()函数对决策树进行剪枝。

# rpart包主要有两个函数组成，分别介绍如下：
rpart(formula, data, weight s, subset, na. action = na. rpart, method, model= FALSE, x= FALSE,y= TRUE, parms, control, cost, . . . )

• fomula 回归方程形式: 例如 y~ x 1+ x2+ x3。
• data 数据: 包含前面方程中变量的数据框( data frame) 。
• na.action 缺失数据的处理办法: 默认办法是删除因变量缺失的观测而保留自变量缺失的观测。
• method 根据树末端的数据类型选择相应变量分割方法,本参数有四种取值: 连续型>anova; 离散型>class; 计数型( 泊松过程)>poisson; 生存分析型>exp。程序会根据因变量的类型自动选择方法, 但一般情况下最好还是指明本参数, 以便让程序清楚做哪一种树模型。
• 连续性method=“anova”,离散型method=“class”,计数型method=“poisson”,生存分布型method=“exp”
• parms 对应于method的选择进行设置，是关于如何分割（splitting）数据的optional选项。在classification splitting中，可供选择进行设置的三个参数包括: 向量的先验概率(component prior)、损失矩阵(component loss)或者分割指数（splitting index）(component split)。
  • component priors必须为正值且和为1，component loss必须对角为0且非对角为正数，component split可以是gini（基尼系数）或者information（信息增益）。
  • 默认值如下，prior为数据计数的占比情况，loss为1，split为gini
  • 如果method = "anova", 则没有参数；
  • 如果method = "poisson", 分割(split)有一个参数，先验分布(prior)变异系数的比率，默认为1；
  • 如果method = "exp", 其参数设置与method = "poisson"一致；
  • 如果为离散型数据，则需要分别设置prior、loss、split；

• control 控制每个节点上的最小样本量、交叉验证的次数、复杂性参量: 即cp: complexity pamemeter, 这个参数意味着对每一步拆分, 模型的拟合优度必须提高的程度, 等等。在rpart的生成对象中即包含cptable

# 交叉验证指，比如xval是10折交叉验证：将数据集分为10组，进行10次拟合，第i次拟合用除了第i组以外的数据训练，用第i组进行预测；目的是减少misclaassification rate。prune(tree, cp, . . . )

• tree 一个回归树对象, 常是rpart()的结果对象。
• cp 复杂性参量, 指定剪枝采用的阈值。

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.81764706      0 1.00000000 1.0000000 0.06194645
# 2 0.04117647      1 0.18235294 0.1823529 0.03169642
# 3 0.01764706      3 0.10000000 0.1588235 0.02970979
# 4 0.01000000      4 0.08235294 0.1235294 0.02637116
plotcp(dtree) # 借助 plotcp() 函数可画出交叉验证误差与复杂度参数的关系图

rpart() 返回的 cptable 值中包括不同大小的树对应的预测误差，因此可用于辅助设定最终的树的大小。其中，复杂度参数（ cp ）用于惩罚过大的树；树的大小即分支数（ nsplit ），有n个分支的树将有n+1个终端节点； rel error 栏即训练集中各种树对应的误差；交叉验证误差（ xerror ）即基于训练样本所得的10折交叉验证误差； xstd 栏为交叉验证误差的标准差。
对于所有交叉验证误差在最小交叉验证误差一个标准差范围内的树，最小的树即最优的树
虚线是基于一个标准差准则得到的上限（0.1235295+1×0.02637116=0.149）。从图像来看，应选择虚线下最左侧 cp 值对应的树

# 根据本例中，最小的交叉验证误差为0.123，标准误差为0.026
# 故最优的树为交叉验证误差在0.123±0.026（0.097和0.149）之间的树
# cp选择0.01，即dtree$cptable中4次分割对应的复杂度符合
dtree.pruned <- prune(dtree, cp = 0.01) library(rpart.plot)
rpart.plot::prp(dtree.pruned, type = 2, extra = 104, fallen.leaves = TRUE,main = "Decision Tree")
# type=2 可画出每个节点下分割的标签
# extra=104 可画出每一类的概率以及每个节点处的样本占比
# fallen.leaves=TRUE 可在图的底端显示终端节点

# 验证模型
dtree.pred <- predict(dtree.pruned, newdata = df.validate, type = "class") 
# 这里与Logistic的type参数不同table(df.validate$class, dtree.pred, dnn = c("Actual", "Predicted"))
###                 Predicted
### Actual      Benign Malignant
### Benign       129        10
### Malignant      2        69

predict() 函数用来对验证集中的观测点分类。
整体来看，验证集中的准确率达到了96%。与逻辑回归不同的是，验证集中的210个样本单元都可由最终树来分类。值得注意的是，对于水平数很多或缺失值很多的预测变量，决策树可能会有偏。

条件推断树

经典决策树的一种重要变体，即条件推断树（conditional inference tree）。条件推断树与经典决策树类似，但变量和分割的选取是基于显著性检验的，而不是纯净度或同质性一类的度量。显著性检验是置换检验。

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

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

library(party)
fit.party <- ctree(class~., data = df.train)
plot(fit.party, main = "Conditinal Inference Tree")df.pred <- predict(fit.party, newdata = df.validate, type = "response")
table(df.validate$class, df.pred, dnn = c("Actual", "Predicted"))
#                  Predicted
# Actual      Benign Malignant
#  Benign       131         8
#  Malignant      4        67

值得注意的是，对于条件推断树来说，剪枝不是必需的, 其生成过程相对更自动化一些。另外， party 包也提供了许多图像参数。每个节点中的阴影区域代表这个节点对应的恶性肿瘤比例。

如果想通过 rpart() 函数得到一棵经典决策树，但想要以上图的形式展示这棵决策树，则可借助 partykit 包。安装并载入这个包后，可通过 plot(as.party(an.rpart.tree))绘制想要的图。

随机森林分类

• 随机森林（random forest）是一种组成式的有监督学习方法。在随机森林中，我们同时生成多个预测模型，并将模型的结果汇总以提升分类准确率。
• 随机森林的算法涉及对样本单元和变量进行抽样，从而生成大量决策树。对每个样本单元来说，所有决策树依次对其进行分类。所有决策树预测类别中的众数类别即为随机森林所预测的这一样本单元的类别。

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

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

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

library(randomForest)
set.seed(1234)
fit.rf <- randomForest(class~., data = df.train, importance = TRUE,na.action = na.roughfix)importance(fit.rf, type = 2)df.pred <- predict(fit.rf, newdata = df.validate)
table(df.validate$class, df.pred)

• 首先， randomForest() 函数从训练集中有放回地随机抽取489个观测点，在每棵树的每个节点随机抽取3个变量，从而生成了500棵传统决策树
• na.action=na.roughfix 参数可将数值变量中的缺失值替换成对应列的中位数，类别变量中的缺失值替换成对应列的众数类（若有多个众数则随机选一个）。
• 随机森林可度量变量重要性，通过设置 information=TRUE 参数得到，并通过 importance()函数输出。由 type=2 参数得到的变量相对重要性就是分割该变量时节点不纯度（异质性）的下降总量对所有树取平均。节点不纯度由Gini系数定义。
• 最后，再通过随机森林算法对验证集中的样本单元进行分类，并计算预测准确率。分类时剔除验证集中有缺失值的单元。

randomForest 包根据传统决策树生成随机森林，而 party 包中的 cforest() 函数则可基于条件推断树生成随机森林。当预测变量间高度相关时，基于条件推断树的随机森林可能效果更好。
相较于其他分类方法，随机森林的分类准确率通常更高。
随机森林算法可处理大规模问题（即多样本单元、多变量），可处理训练集中有大量缺失值的数据，也可应对变量远多于样本单元的数据。可计算袋外预测误差（OOB error）、度量变量重要性也是随机森林的两个明显优势。
随机森林的一个明显缺点是分类方法（此例中相当于500棵决策树）较难理解和表达。另外，我们需要存储整个随机森林以对新样本单元分类。

支持向量机分类

支持向量机（SVM）是一类可用于分类和回归的有监督机器学习模型。其流行归功于两个方面：一方面，他们可输出较准确的预测结果；另一方面，模型基于较优雅的数学理论。

SVM旨在在多维空间中找到一个能将全部样本单元分成两类的最优平面，这一平面应使两类中距离最近的点的间距（margin）尽可能大，在间距边界上的点被称为支持向量（support vector，它们决定间距），分割的超平面位于间距的中间。

SVM可以通过R中 kernlab 包的 ksvm() 函数和 e1071 包中的 svm() 函数实现。 ksvm() 功能更强大，但 svm() 相对更简单。

library(e1071)
fit.svm <- svm(class~., df.train)
fit.svmsvm.pred <- predict(fit.svm, newdata = na.omit(df.validate))
# svm不支持缺失值
table(df.validate$calss, svm.pred, dnn = c("Actual", "Predicted"))
###                 Predicted
### Actual      Benign Malignant
### Benign       136         8
### Malignant      1        65

由于方差较大的预测变量通常对SVM的生成影响更大， svm() 函数默认在生成模型前对每个变量标准化，使其均值为0、标准差为1。

选择调和参数

svm() 函数默认通过径向基函数（Radial Basis Function，RBF）将样本单元投射到高维空间。

在用带RBF核的SVM拟合样本时，两个参数可能影响最终结果：gamma和成本（cost）。

gamma是核函数的参数，控制分割超平面的形状。gamma越大，通常导致支持向量越多。我们也可将gamma看作控制训练样本“到达范围”的参数，即gamma越大意味着训练样本到达范围越广，而越小则意味着到达范围越窄。gamma必须大于0。

成本参数代表犯错的成本。一个较大的成本意味着模型对误差的惩罚更大，从而将生成一个更复杂的分类边界，对应的训练集中的误差也会更小，但也意味着可能存在过拟合问题，即对新样本单元的预测误差可能很大。较小的成本意味着分类边界更平滑，但可能会导致欠拟合。与gamma一样，成本参数也恒为正。

svm() 函数默认设置gamma为预测变量个数的倒数，成本参数为1。不过gamma与成本参数的不同组合可能生成更有效的模型。在建模时，我们可以尝试变动参数值建立不同的模型，但利用格点搜索法可能更有效。可以通过tune.svm()对每个参数设置一个候选范围， tune.svm()函数对每一个参数组合生成一个SVM模型，并输出在每一个参数组合上的表现。

set.seed(1234)
tuned <- tune.svm(class~., data = df.train, gamma = 10^(-6:1), cost = 10^(-10:10))set.seed(1234)
tuned <- tune.svm(class~., data = df.train, gamma = 10^(-6:1), cost = 10^(-10:10))
tuned
Parameter tuning of ‘svm’:
# sampling method: 10-fold cross validation 
# best parameters:
# gamma cost
# 0.001   10
# best performance: 0.0327381 fit.svm2 <- svm(class~., data = df.train, gamma = 0.001, cost = 10)
svm.pred2 <- predict(fit.svm2, newdata = na.omit(df.validate))
table(df.validate$class, svm.pred2)
#           Predicted
# Actual      Benign Malignant
#  Benign       138         6
#  Malignant      1        65

• 首先，对不同的gamma和成本拟合一个带RBF核的SVM模型。我们一共将尝试八个不同的gamma （从0.000 001到10）以及21个成本参数（从0.01到1010）。总体来说，我们共拟合了168 （8×21）个模型，并比较了其结果。训练集中10折交叉验证误差最小的模型所对应的参数为gamm=0.1，成本参数为1。
• 基于这一参数值组合，我们对全部训练样本拟合出新的SVM模型，然后用这一模型对验证集中的样本单元进行预测，并给出错分个数。在本例中，调和后的模型轻微减少了错分个数（从7减少到6）。一般来说，为SVM模型选取调和参数通常可以得到更好的结果。

选择预测效果最好的解

涉及一个分类器的敏感度（sensitivity）、特异性（sensitivity）、正例命中率（positive predictive power）和负例命中率（negative predictive power）。

通过比较上述值在不同模型中的表现，可以得出最佳模型。

predict() 函数可以估计一个样本单元为恶性组织的概率。如果这一概率值大于0.5，则分类器会把这一样本单元判为恶性。这个0.5即阈值（threshold）或门槛值（cutoff value）。通过变动这一阈值，我们可以通过牺牲分类器的特异性来增加其敏感度。这同样适用于决策树、随机森林和支持向量机（尽管语句写法上会有差别）。

变动阈值可能带来的影响可以通过ROC（Receiver Operating Characteristic）曲线来进一步观察。ROC曲线可对一个区间内的门槛值画出特异性和敏感度之间的关系，然后我们就能针对特定问题选择特异性和敏感度的最佳组合。许多R包都可以画ROC曲线，如 ROCR 、 pROC 等。

ROC & AUC curve with pROC

The best threshold is of as high as possible true positive rate while as low as possible false positive rate. This decision made through the graph is the point the farest away from the green line.

The better method decided by the AUC value, the higher the better. Through the graph, the red one is better than the green one.

pROC的roc()函数的两个主要参数response和predictor由原始数据或者待测数据的分类结果表示，默认为0为control而1为case。如果为多分类变量，则要用levels指明control和case。predictor则为probability分布值。在logistic regression中为fit$fitted.value，在randomForest中为fit$votes[,]的具体列。另外，需要指出的是这两个参数必须长度相等。而回归分析往往对NA值敏感，会自动剔除，导致长度不等，因此需要注意。

library(pROC)
library(randomForest)
fit.logit <- glm(class~., data = df.train, family = binomial())
# the variable values and fitted.value stored in fit.logit into the lines() functionfit.rf <- randomForest(class~., data = df.train, importance = TRUE,na.action = na.roughfix)
roc(df.train$class, fit.logit$fitted.value,plot=TRUE, legacy.axes = TRUE, percent = TRUE,print.auc = TRUE, col = "red")
# 这里理论上会报错，因为进行glm（）回归时，含有NA值的那一行会自动忽略
# 因此最后response和predictor的length是不一样的。
# 这种情况下应该用df.train[complete.cases(df.train),]$class
# legacy.axes = TRUE表示1-specificity
# percent = TRUE表示用百分数表示
# print.auc = TRUE显示曲线下面积# 比较两种方法的ROC曲线，绘制于同一图像上
plot.roc(df.train$class, fit.rf$votes[,1], percent = TRUE, print.auc = TRUE, add = TRUE, levels = c(control = "Benign", case = "Malignant"), col = "blue")
# randomForest对NA不敏感
# print.auc默认打印在中间，如果要调整，选用print.auc.y或者print.auc.x = 40(百分比位置)
# levels = 设置class的类型便于取舍
# 多种方法时要区分col# legend
legend("bottomright", legend = c("Logistic regression", "Random Forest"),col = c("red", "blue"), lwd =4) # 添加Legend时保持颜色一致

对于更改两边margin，可以用par(pty = "s")。默认“m”

更新于

2020-10-27