R-机器学习秘籍-三-

R 机器学习秘籍（三）

原文：annas-archive.org/md5/167951e4c9873bc214756c25a63518c8

译者：飞龙

协议：CC BY-NC-SA 4.0

第七章：模型评估

在本章中，我们将涵盖以下主题：

• 使用 k 折交叉验证估计模型性能

• 使用 e1071 包进行交叉验证

• 使用caret包进行交叉验证

• 使用caret包对变量重要性进行排名

• 使用rminer包对变量重要性进行排名

• 使用caret包寻找高度相关的特征

• 使用caret包选择特征

• 测量回归模型的性能

• 使用混淆矩阵测量预测性能

• 使用 ROCR 测量预测性能

• 使用caret包比较 ROC 曲线

• 使用caret包测量模型之间的性能差异

简介

模型评估是为了确保拟合的模型能够准确预测未来或未知主体的响应。如果没有模型评估，我们可能会训练出在训练数据上过度拟合的模型。为了防止过度拟合，我们可以使用如caret、rminer和rocr等包来评估拟合模型的性能。此外，模型评估有助于选择最佳模型，该模型更稳健，并能准确预测未来主体的响应。

在下一章中，我们将讨论如何实现一个简单的 R 脚本或使用一个包（例如caret或rminer）来评估拟合模型的性能。

使用 k 折交叉验证估计模型性能

k 折交叉验证技术是一种常用的技术，用于估计分类器的性能，因为它克服了过度拟合的问题。对于 k 折交叉验证，该方法不使用整个数据集来构建模型，而是将数据分割成训练数据集和测试数据集。因此，使用训练数据集构建的模型可以用来评估模型在测试数据集上的性能。通过执行 n 次 k 折验证，我们可以使用 n 个准确率的平均值来真正评估构建的模型性能。在本菜谱中，我们将说明如何执行 k 折交叉验证。

准备工作

在这个菜谱中，我们将继续使用 telecom churn数据集作为输入数据源来训练支持向量机。对于那些尚未准备数据集的人，请参阅第五章，分类（I）-树、懒惰和概率性，以获取详细信息。

如何操作...

执行以下步骤以交叉验证 telecom churn数据集：

1. 使用cut函数将索引分割成 10 折：

   > ind = cut(1:nrow(churnTrain), breaks=10, labels=F)
   
   

2. 接下来，使用for循环执行 10 折交叉验证，重复 10 次：

   > accuracies = c()
   > for (i in 1:10) {
   +   fit = svm(churn ~., churnTrain[ind != i,])
   +   predictions = predict(fit, churnTrain[ind == i, ! names(churnTrain) %in% c("churn")])
   +   correct_count = sum(predictions == churnTrain[ind == i,c("churn")])
   +   accuracies = append(correct_count / nrow(churnTrain[ind == i,]), accuracies)
   + }
   
   

3. 你可以打印出准确率：

   > accuracies
    [1] 0.9341317 0.8948949 0.8978979 0.9459459 0.9219219 0.9281437 0.9219219 0.9249249 0.9189189 0.9251497
   
   

4. 最后，你可以使用mean函数生成平均准确率：

   > mean(accuracies)
   [1] 0.9213852
   
   

它是如何工作的...

在本例中，我们实现了一个简单的脚本，执行 10 折交叉验证。我们首先使用cut函数生成 10 个折的索引。然后，我们实现一个for循环，执行 10 次 10 折交叉验证。在循环中，我们首先将svm应用于9个数据折作为训练集。然后，我们使用拟合的模型预测剩余数据（测试数据集）的标签。接下来，我们使用正确预测的标签总和来生成准确率。因此，循环存储了 10 个生成的准确率。最后，我们使用mean函数检索准确率的平均值。

还有更多...

如果您希望使用其他模型执行 k 折交叉验证，只需替换生成变量 fit 的行，以您偏好的分类器为准。例如，如果您想使用 10 折交叉验证来评估朴素贝叶斯模型，只需将调用函数从svm替换为naiveBayes：

> for (i in 1:10) {
+   fit = naiveBayes(churn ~., churnTrain[ind != i,])
+   predictions = predict(fit, churnTrain[ind == i, ! names(churnTrain) %in% c("churn")])
+   correct_count = sum(predictions == churnTrain[ind == i,c("churn")])
+   accuracies = append(correct_count / nrow(churnTrain[ind == i,]), accuracies)
+ }

使用 e1071 包进行交叉验证

除了实现一个loop函数来执行 k 折交叉验证外，您还可以在e1071包中使用tuning函数（例如，tune.nnet、tune.randomForest、tune.rpart、tune.svm和tune.knn）来获取最小误差值。在本例中，我们将说明如何使用tune.svm执行 10 折交叉验证并获得最佳分类模型。

准备工作

在本例中，我们继续使用电信churn数据集作为输入数据源执行 10 折交叉验证。

如何操作...

执行以下步骤以使用交叉验证检索最小估计误差：

1. 在训练数据集trainset上应用tune.svm，使用 10 折交叉验证作为调整控制。（如果您发现错误消息，例如could not find function predict.func，请清除工作区，重新启动 R 会话并重新加载e1071库）：

   > tuned = tune.svm(churn~., data = trainset, gamma = 10^-2, cost = 10², tunecontrol=tune.control(cross=10))
   
   

2. 接下来，您可以获取模型的摘要信息，调整：

   > summary(tuned)
   
   Error estimation of 'svm' using 10-fold cross validation: 0.08164651
   
   

3. 然后，您可以访问调整后模型的性能细节：

   > tuned$performances
    gamma cost      error dispersion
   1  0.01  100 0.08164651 0.02437228
   
   

4. 最后，您可以使用最佳模型生成分类表：

   > svmfit = tuned$best.model
   > table(trainset[,c("churn")], predict(svmfit))
   
    yes   no
    yes  234  108
    no    13 1960
   
   

工作原理...

e1071包提供了构建和评估模型的各种函数，因此，您无需重新发明轮子来评估拟合模型。在本例中，我们使用tune.svm函数使用给定的公式、数据集、gamma、成本和控制函数调整 svm 模型。在tune.control选项中，我们将选项配置为cross=10，在调整过程中执行 10 折交叉验证。调整过程最终将返回最小估计误差、性能细节以及调整过程中的最佳模型。因此，我们可以获得调整的性能指标，并进一步使用最佳模型生成分类表。

参考信息

• 在e1071包中，tune函数使用网格搜索来调整参数。对于那些对其他调整函数感兴趣的人，请使用帮助函数查看tune文档：

  > ?e1071::tune
  
  

使用 caret 包进行交叉验证

Caret（分类和回归训练）包包含许多关于回归和分类问题训练过程的函数。类似于e1071包，它也包含一个执行 k 折交叉验证的函数。在本菜谱中，我们将演示如何使用caret包执行 k 折交叉验证。

准备工作

在本菜谱中，我们将继续使用电信churn数据集作为输入数据源来执行 k 折交叉验证。

如何操作...

执行以下步骤以使用caret包执行 k 折交叉验证：

1. 首先，设置控制参数以进行 10 折交叉验证，重复 3 次：

   > control = trainControl(method="repeatedcv", number=10, repeats=3)
   
   

2. 然后，你可以使用rpart在电信客户流失数据上训练分类模型：

   > model = train(churn~., data=trainset, method="rpart", preProcess="scale", trControl=control)
   
   

3. 最后，你可以检查生成的模型的输出：

   > model
   CART 
   
   2315 samples
    16 predictor
    2 classes: 'yes', 'no' 
   
   Pre-processing: scaled 
   Resampling: Cross-Validated (10 fold, repeated 3 times) 
   
   Summary of sample sizes: 2084, 2083, 2082, 2084, 2083, 2084, ... 
   
   Resampling results across tuning parameters:
   
    cp      Accuracy  Kappa  Accuracy SD  Kappa SD
    0.0556  0.904     0.531  0.0236       0.155 
    0.0746  0.867     0.269  0.0153       0.153 
    0.0760  0.860     0.212  0.0107       0.141 
   
   Accuracy was used to select the optimal model using the largest value.
   The final value used for the model was cp = 0.05555556.
   
   

工作原理...

在本菜谱中，我们展示了使用caret包进行 k 折交叉验证是多么方便。在第一步中，我们设置了训练控制并选择了在三次重复中执行 10 折交叉验证的选项。重复 k 折验证的过程称为重复 k 折验证，用于测试模型的稳定性。如果模型稳定，应该得到相似的训练结果。然后，我们使用rpart在训练数据集上应用，并选择缩放数据以及使用之前步骤中配置的选项来训练模型。

训练过程完成后，模型输出三个重采样结果。在这些结果中，cp=0.05555556的模型具有最大的准确值（0.904），因此被选为分类的最佳模型。

参考以下内容

• 你可以在trainControl中配置resampling函数，其中你可以指定boot、boot632、cv、repeatedcv、LOOCV、LGOCV、none、oob、adaptive_cv、adaptive_boot或adaptive_LGOCV。要查看如何选择重采样方法的更详细信息，请查看trainControl文档：

  > ?trainControl
  
  

使用 caret 包对变量重要性进行排名

在构建监督学习模型后，我们可以估计特征的重要性。这种估计采用敏感性分析来衡量当输入变化时对给定模型输出的影响。在本菜谱中，我们将向您展示如何使用caret包对变量重要性进行排名。

准备工作

你需要完成之前的菜谱，并将拟合的rpart对象存储在model变量中。

如何操作...

执行以下步骤以使用caret包对变量重要性进行排名：

1. 首先，你可以使用varImp函数估计变量重要性：

   > importance = varImp(model, scale=FALSE)
   > importance
   rpart variable importance
   
    Overall
   number_customer_service_calls 116.015
   total_day_minutes             106.988
   total_day_charge              100.648
   international_planyes          86.789
   voice_mail_planyes             25.974
   total_eve_charge               23.097
   total_eve_minutes              23.097
   number_vmail_messages          19.885
   total_intl_minutes              6.347
   total_eve_calls                 0.000
   total_day_calls                 0.000
   total_night_charge              0.000
   total_intl_calls                0.000
   total_intl_charge               0.000
   total_night_minutes             0.000
   total_night_calls               0.000
   
   

2. 然后，您可以使用plot函数生成可变重要性图：

   > plot(importance)
   
   

   

   图 1：使用 caret 包的可变重要性可视化

工作原理...

在本食谱中，我们首先使用varImp函数检索可变重要性并获取摘要。整体结果显示了每个属性的敏感性度量。接下来，我们按排名绘制可变重要性，这表明number_customer_service_calls属性在敏感性度量中是最重要的变量。

更多内容...

在某些分类包中，例如rpart，从训练模型生成的对象包含变量重要性。我们可以通过访问输出对象来检查变量重要性：

> library(rpart)
> model.rp = rpart(churn~., data=trainset)
> model.rp$variable.importance
 total_day_minutes              total_day_charge 
 111.645286                    110.881583 
number_customer_service_calls            total_intl_minutes 
 58.486651                     48.283228 
 total_intl_charge              total_eve_charge 
 47.698379                     47.166646 
 total_eve_minutes            international_plan 
 47.166646                     42.194508 
 total_intl_calls         number_vmail_messages 
 36.730344                     19.884863 
 voice_mail_plan             total_night_calls 
 19.884863                      7.195828 
 total_eve_calls            total_night_charge 
 3.553423                      1.754547 
 total_night_minutes               total_day_calls 
 1.754547                      1.494986 

使用 rminer 包对可变重要性进行排名

除了使用caret包生成可变重要性外，您还可以使用rminer包生成分类模型的可变重要性。在下面的食谱中，我们将说明如何使用rminer获取拟合模型的可变重要性。

准备工作

在本食谱中，我们将继续使用电信churn数据集作为输入数据源来对可变重要性进行排序。

如何操作...

执行以下步骤以使用rminer对可变重要性进行排名：

1. 安装并加载包，rminer：

   > install.packages("rminer")
   > library(rminer)
   
   

2. 使用训练集拟合 svm 模型：

   > model=fit(churn~.,trainset,model="svm")
   
   

3. 使用Importance函数获取可变重要性：

   > VariableImportance=Importance(model,trainset,method="sensv")
   
   

4. 按方差绘制可变重要性图

   > L=list(runs=1,sen=t(VariableImportance$imp),sresponses=VariableImportance$sresponses)
   > mgraph(L,graph="IMP",leg=names(trainset),col="gray",Grid=10)
   
   

   

   图 2：使用rminer包的可变重要性可视化

工作原理...

与caret包类似，rminer包也可以生成分类模型的可变重要性。在本食谱中，我们首先使用fit函数在训练数据集trainset上训练 svm 模型。然后，我们使用Importance函数使用敏感性度量对可变重要性进行排名。最后，我们使用mgraph绘制可变重要性的排名。与使用caret包获得的结果相似，number_customer_service_calls是在敏感性度量中最重要的变量。

参见

• rminer包提供许多分类模型供用户选择。如果您对使用 svm 以外的模型感兴趣，可以使用以下命令查看这些选项：

  > ?rminer::fit
  
  

使用 caret 包寻找高度相关的特征

在执行回归或分类时，如果移除高度相关的属性，某些模型的表现会更好。caret包提供了findCorrelation函数，可以用来查找彼此高度相关的属性。在本食谱中，我们将演示如何使用caret包查找高度相关的特征。

准备工作

在本食谱中，我们将继续使用电信churn数据集作为输入数据源来寻找高度相关的特征。

如何操作...

执行以下步骤以找到高度相关的属性：

1. 移除未用数值字符编码的特征：

   > new_train = trainset[,! names(churnTrain) %in% c("churn", "international_plan", "voice_mail_plan")]
   
   

2. 然后，你可以获得每个属性的关联性：

   >cor_mat = cor(new_train)
   
   

3. 接下来，我们使用findCorrelation搜索截止值为 0.75 的高度相关属性：

   > highlyCorrelated = findCorrelation(cor_mat, cutoff=0.75)
   
   

4. 然后，我们获得高度相关属性的名称：

   > names(new_train)[highlyCorrelated]
   [1] "total_intl_minutes"  "total_day_charge"    "total_eve_minutes"   "total_night_minutes"
   
   

它是如何工作的...

在本菜谱中，我们使用caret包搜索高度相关的属性。为了检索每个属性的关联性，首先应移除非数值属性。然后，我们执行关联性分析以获得关联矩阵。接下来，我们使用findCorrelation找到高度相关的属性，截止值为 0.75。我们最终获得高度相关（相关系数超过 0.75）的属性名称，分别是total_intl_minutes、total_day_charge、total_eve_minutes和total_night_minutes。你可以考虑移除一些高度相关的属性，保留一个或两个属性以获得更好的准确性。

相关内容

• 除了caret包外，你还可以使用subselect包中的leaps、genetic和anneal函数达到相同的目的

使用 caret 包进行特征选择

特征选择方法搜索具有最小预测误差的特征子集。我们可以应用特征选择来识别构建准确模型所需的属性。caret包提供了一个递归特征消除函数rfe，可以帮助自动选择所需的特征。在下面的菜谱中，我们将演示如何使用caret包进行特征选择。

准备工作

在本菜谱中，我们将继续使用电信churn数据集作为特征选择的输入数据源。

如何操作...

执行以下步骤以选择特征：

1. 将训练数据集trainset中名为international_plan的特征转换为intl_yes和intl_no：

   > intl_plan = model.matrix(~ trainset.international_plan - 1, data=data.frame(trainset$international_plan))
   > colnames(intl_plan) = c("trainset.international_planno"="intl_no", "trainset.international_planyes"= "intl_yes")
   
   

2. 将训练数据集trainset中名为voice_mail_plan的特征转换为voice_yes和voice_no：

   > voice_plan = model.matrix(~ trainset.voice_mail_plan - 1, data=data.frame(trainset$voice_mail_plan))
   > colnames(voice_plan) = c("trainset.voice_mail_planno" ="voice_no", "trainset.voice_mail_planyes"="voidce_yes")
   
   

3. 移除international_plan和voice_mail_plan属性，并将训练数据集trainset与数据框intl_plan和voice_plan合并：

   > trainset$international_plan = NULL
   > trainset$voice_mail_plan = NULL
   > trainset = cbind(intl_plan,voice_plan, trainset)
   
   

4. 将测试数据集testset中名为international_plan的特征转换为intl_yes和intl_no：

   > intl_plan = model.matrix(~ testset.international_plan - 1, data=data.frame(testset$international_plan))
   > colnames(intl_plan) = c("testset.international_planno"="intl_no", "testset.international_planyes"= "intl_yes")
   
   

5. 将训练数据集trainset中名为voice_mail_plan的特征转换为voice_yes和voice_no：

   > voice_plan = model.matrix(~ testset.voice_mail_plan - 1, data=data.frame(testset$voice_mail_plan))
   > colnames(voice_plan) = c("testset.voice_mail_planno" ="voice_no", "testset.voice_mail_planyes"="voidce_yes")
   
   

6. 移除international_plan和voice_mail_plan属性，并将测试数据集testset与数据框intl_plan和voice_plan合并：

   > testset$international_plan = NULL
   > testset$voice_mail_plan = NULL
   > testset = cbind(intl_plan,voice_plan, testset)
   
   

7. 然后，我们使用线性判别分析创建一个特征选择算法：

   > ldaControl = rfeControl(functions = ldaFuncs, method = "cv")
   
   

8. 接下来，我们在训练数据集trainset上使用从 1 到 18 的子集进行向后特征选择：

   > ldaProfile = rfe(trainset[, !names(trainset) %in% c("churn")], trainset[,c("churn")],sizes = c(1:18), rfeControl = ldaControl)
   > ldaProfile
   
   Recursive feature selection
   
   Outer resampling method: Cross-Validated (10 fold) 
   
   Resampling performance over subset size:
   
    Variables Accuracy  Kappa AccuracySD KappaSD Selected
    1   0.8523 0.0000   0.001325 0.00000 
    2   0.8523 0.0000   0.001325 0.00000 
    3   0.8423 0.1877   0.015468 0.09787 
    4   0.8462 0.2285   0.016593 0.09610 
    5   0.8466 0.2384   0.020710 0.09970 
    6   0.8466 0.2364   0.019612 0.09387 
    7   0.8458 0.2315   0.017551 0.08670 
    8   0.8458 0.2284   0.016608 0.09536 
    9   0.8475 0.2430   0.016882 0.10147 
    10   0.8514 0.2577   0.014281 0.08076 
    11   0.8518 0.2587   0.014124 0.08075 
    12   0.8544 0.2702   0.015078 0.09208        *
    13   0.8544 0.2721   0.015352 0.09421 
    14   0.8531 0.2663   0.018428 0.11022 
    15   0.8527 0.2652   0.017958 0.10850 
    16   0.8531 0.2684   0.017897 0.10884 
    17   0.8531 0.2684   0.017897 0.10884 
    18   0.8531 0.2684   0.017897 0.10884 
   
   The top 5 variables (out of 12):
    total_day_charge, total_day_minutes, intl_no, number_customer_service_calls, total_eve_charge
   
   

9. 接下来，我们可以绘制选择结果：

   > plot(ldaProfile, type = c("o", "g"))
   
   

   

   图 3：特征选择结果

10. 然后，我们可以检查变量的最佳子集：

    > ldaProfile$optVariables
     [1] "total_day_charge" 
     [2] "total_day_minutes" 
     [3] "intl_no" 
     [4] "number_customer_service_calls"
     [5] "total_eve_charge" 
     [6] "total_eve_minutes" 
     [7] "voidce_yes" 
     [8] "total_intl_calls" 
     [9] "number_vmail_messages" 
    [10] "total_intl_charge" 
    [11] "total_intl_minutes" 
    [12] "total_night_minutes" 
    
    

11. 现在，我们可以检查拟合的模型：

    > ldaProfile$fit
    Call:
    lda(x, y)
    
    Prior probabilities of groups:
     yes        no 
    0.1477322 0.8522678 
    
    Group means:
     total_day_charge total_day_minutes   intl_no
    yes         35.00143          205.8877 0.7046784
    no          29.62402          174.2555 0.9351242
     number_customer_service_calls total_eve_charge
    yes                      2.204678         18.16702
    no                       1.441460         16.96789
     total_eve_minutes voidce_yes total_intl_calls
    yes          213.7269  0.1666667         4.134503
    no           199.6197  0.2954891         4.514445
     number_vmail_messages total_intl_charge
    yes              5.099415          2.899386
    no               8.674607          2.741343
     total_intl_minutes total_night_minutes
    yes           10.73684            205.4640
    no            10.15119            201.4184
    
    Coefficients of linear discriminants:
     LD1
    total_day_charge               0.715025524
    total_day_minutes             -0.130486470
    intl_no                        2.259889324
    number_customer_service_calls -0.421997335
    total_eve_charge              -2.390372793
    total_eve_minutes              0.198406977
    voidce_yes                     0.660927935
    total_intl_calls               0.066240268
    number_vmail_messages         -0.003529233
    total_intl_charge              2.315069869
    total_intl_minutes            -0.693504606
    total_night_minutes           -0.002127471
    
    

12. 最后，我们可以计算跨样本的性能：

    > postResample(predict(ldaProfile, testset[, !names(testset) %in% c("churn")]), testset[,c("churn")])
    Accuracy     Kappa
    0.8605108 0.2672027
    
    

工作原理...

在这个菜谱中，我们使用caret包进行特征选择。由于数据集中存在因子编码的属性，我们首先使用一个名为model.matrix的函数将因子编码的属性转换为多个二元属性。因此，我们将international_plan属性转换为intl_yes和intl_no。此外，我们将voice_mail_plan属性转换为voice_yes和voice_no。

接下来，我们使用交叉验证方法cv和线性判别函数ldaFuncs设置训练的控制参数。然后，我们使用control函数ldaFuncs进行特征选择，使用递归特征消除rfe。rfe函数生成特征选择摘要，其中包含对子集大小和顶级变量的重采样性能。

然后，我们可以使用获得模型信息来绘制变量数量与准确率的关系图。从图 3 中可以看出，使用 12 个特征可以获得最佳准确率。除此之外，我们还可以检索出拟合模型中最佳变量子集（总共有 12 个变量）。最后，我们可以计算跨样本的性能，得到准确率为 0.86 和 kappa 值为 0.27。

参见

• 为了指定用于控制特征选择的算法，可以在rfeControl中更改指定的控制函数。以下是一些你可以使用的选项：

  caretFuncs      SVM (caret)
  lmFuncs     lm (base)
  rfFuncs         RF(randomForest)
  treebagFuncs     DT (ipred)
  ldaFuncs       lda(base)
  nbFuncs       NB(klaR)
  gamFuncs      gam(gam)
  
  

衡量回归模型的表现

为了衡量回归模型的表现，我们可以计算预测输出与实际输出之间的距离，作为模型性能的量化指标。在这里，我们通常使用均方根误差（RMSE）、相对平方误差（RSE）和 R-Square 作为常见的测量指标。在下面的菜谱中，我们将展示如何从一个构建的回归模型中计算这些测量值。

准备工作

在这个菜谱中，我们将使用包含四个回归数据集的Quartet数据集作为我们的输入数据源。

如何做...

执行以下步骤来衡量回归模型的表现：

1. 从car包中加载Quartet数据集：

   > library(car)
   > data(Quartet)
   
   

2. 使用lm函数将属性y3与 x 绘制出来：

   > plot(Quartet$x, Quartet$y3)
   > lmfit = lm(Quartet$y3~Quartet$x)
   > abline(lmfit, col="red")
   
   

   

   图 4：线性回归图

3. 你可以使用predict函数检索预测值：

   > predicted= predict(lmfit, newdata=Quartet[c("x")])
   
   

4. 现在，你可以计算均方根误差：

   > actual = Quartet$y3
   > rmse = (mean((predicted - actual)²))⁰.5
   > rmse
   [1] 1.118286
   
   

5. 你可以计算相对平方误差：

   > mu = mean(actual)
   > rse = mean((predicted - actual)²) / mean((mu - actual)²) 
   > rse
   [1] 0.333676
   
   

6. 此外，你还可以使用 R-Square 作为测量指标：

   > rsquare = 1 - rse
   > rsquare
   [1] 0.666324
   
   

7. 然后，你可以使用 MASS 包中的rlm函数将属性 y3 与 x 绘制出来：

   > library(MASS)
   > plot(Quartet$x, Quartet$y3)
   > rlmfit = rlm(Quartet$y3~Quartet$x)
   > abline(rlmfit, col="red")
   
   

   

   图 5：Quartet 数据集上的稳健线性回归图

8. 然后，你可以使用predict函数检索预测值：

   > predicted = predict(rlmfit, newdata=Quartet[c("x")])
   
   

9. 接下来，你可以使用预测值和实际值之间的距离来计算根均方误差：

   > actual = Quartet$y3
   > rmse = (mean((predicted - actual)²))⁰.5
   > rmse
   [1] 1.279045
   
   

10. 计算预测标签和实际标签之间的相对平方误差：

    > mu = mean(actual)
    > rse =mean((predicted - actual)²) / mean((mu - actual)²) 
    > rse
    [1] 0.4365067
    
    

11. 现在，你可以计算 R-Square 值：

    > rsquare = 1 - rse
    > rsquare
    [1] 0.5634933
    
    

工作原理...

回归模型性能的测量使用预测值和实际值之间的距离。我们通常使用这三个度量，即均方根误差、相对平方误差和 R-Square，作为回归模型性能的量化指标。在这个菜谱中，我们首先从 car 包中加载 Quartet 数据。然后，我们使用 lm 函数拟合线性模型，并在 x 变量对 y3 变量的散点图上添加回归线。接下来，我们使用预测函数计算预测值，并开始计算构建模型的 均方根误差（RMSE）、相对平方误差（RSE）和 R-Square。

由于这个数据集在 x=13 处有一个异常值，我们希望量化异常值对性能测量的影响。为了实现这一点，我们首先使用 MASS 包中的 rlm 函数训练一个回归模型。类似于前面的步骤，然后我们生成根平方均方误差、相对误差和 R-Square 的性能测量。从输出测量来看，显然 lm 模型的均方误差和相对平方误差小于由 rlm 构建的模型，R-Square 的得分显示使用 lm 构建的模型具有更大的预测能力。然而，对于实际场景，我们应该移除 x=13 处的异常值。这个比较表明异常值可能偏向于性能测量，并可能导致我们选择错误的模型。

更多内容...

如果你想要对一个线性回归模型进行交叉验证，你可以在 e1071 包中使用 tune 函数：

> tune(lm, y3~x, data = Quartet)
Error estimation of 'lm' using 10-fold cross validation: 2.33754

除了 e1071 包之外，你还可以使用 caret 包中的 train 函数来进行交叉验证。此外，你还可以使用 DAAG 包中的 cv.lm 来达到相同的目的。

使用混淆矩阵衡量预测性能

为了衡量分类模型的性能，我们首先可以根据我们的预测标签和实际标签生成一个分类表。然后，我们可以使用混淆矩阵来获得性能度量，如精确度、召回率、特异性和准确率。在这个菜谱中，我们将演示如何使用 caret 包检索混淆矩阵。

准备工作

在这个菜谱中，我们将继续使用电信 churn 数据集作为我们的示例数据集。

如何做...

执行以下步骤以生成分类测量：

1. 使用训练数据集训练 svm 模型：

   > svm.model= train(churn ~ .,
   +                   data = trainset,
   +                   method = "svmRadial")
   
   

2. 然后，你可以使用拟合的模型 svm.model 来预测标签：

   > svm.pred = predict(svm.model, testset[,! names(testset) %in% c("churn")])
   
   

3. 接下来，你可以生成一个分类表：

   > table(svm.pred, testset[,c("churn")])
   
   svm.pred yes  no
    yes  73  16
    no   68 861
   
   

4. 最后，你可以使用预测结果和测试数据集中的实际标签生成一个混淆矩阵：

   > confusionMatrix(svm.pred, testset[,c("churn")])
   Confusion Matrix and Statistics
   
    Reference
   Prediction yes  no
    yes  73  16
    no   68 861
   
    Accuracy : 0.9175 
    95% CI : (0.8989, 0.9337)
    No Information Rate : 0.8615 
    P-Value [Acc > NIR] : 2.273e-08 
   
    Kappa : 0.5909 
    Mcnemar's Test P-Value : 2.628e-08 
   
    Sensitivity : 0.51773 
    Specificity : 0.98176 
    Pos Pred Value : 0.82022 
    Neg Pred Value : 0.92680 
    Prevalence : 0.13851 
    Detection Rate : 0.07171 
    Detection Prevalence : 0.08743 
    Balanced Accuracy : 0.74974 
   
    'Positive' Class : yes 
   
   

工作原理...

在这个菜谱中，我们展示了如何获取混淆矩阵来衡量分类模型的性能。首先，我们使用caret包中的train函数来训练 svm 模型。接下来，我们使用predict函数使用测试数据集提取 svm 模型的预测标签。然后，我们执行table函数以根据预测和实际标签获取分类表。最后，我们使用caret包中的confusionMatrix函数生成一个混淆矩阵来衡量分类模型的性能。

相关内容

• 如果你感兴趣，可以在train函数中使用的方法，可以参考这个网站：topepo.github.io/caret/modelList.html

使用 ROCR 测量预测性能

接收者操作特征（ROC）曲线是说明二元分类器系统性能的图表，它将不同截断点的真正例率和假正例率进行绘图。我们最常用这个图表来计算曲线下面积（AUC），以衡量分类模型的性能。在这个菜谱中，我们将展示如何绘制 ROC 曲线并计算 AUC 来衡量分类模型的性能。

准备工作

在这个菜谱中，我们将继续使用电信churn数据集作为我们的示例数据集。

如何操作...

执行以下步骤以生成具有不同成本的两种不同的分类示例：

1. 首先，你应该安装并加载ROCR包：

   > install.packages("ROCR")
   > library(ROCR)
   
   

2. 使用训练数据集并设置概率等于TRUE来训练 svm 模型：

   > svmfit=svm(churn~ ., data=trainset, prob=TRUE)
   
   

3. 在测试数据集上基于训练模型进行预测，并将概率设置为TRUE：

   >pred=predict(svmfit,testset[, !names(testset) %in% c("churn")], probability=TRUE)
   
   

4. 使用yes标签获取概率：

   > pred.prob = attr(pred, "probabilities") 
   > pred.to.roc = pred.prob[, 2] 
   
   

5. 使用prediction函数生成预测结果：

   > pred.rocr = prediction(pred.to.roc, testset$churn)
   
   

6. 使用performance函数获取性能度量：

   > perf.rocr = performance(pred.rocr, measure = "auc", x.measure = "cutoff") 
   > perf.tpr.rocr = performance(pred.rocr, "tpr","fpr") 
   
   

7. 使用plot函数可视化 ROC 曲线：

   > plot(perf.tpr.rocr, colorize=T,main=paste("AUC:",(perf.rocr@y.values)))
   
   

   

   图 6：svm 分类器性能的 ROC 曲线

工作原理...

在本例中，我们展示了如何生成 ROC 曲线来展示二元分类器的性能。首先，我们应该安装并加载库，ROCR。然后，我们使用来自e1071包的 svm 来训练一个分类模型，并使用该模型对测试数据集进行预测。接下来，我们使用预测函数（来自ROCR包）生成预测结果。然后，我们将性能函数调整为获得真正例率对假正例率的性能测量。最后，我们使用plot函数可视化 ROC 图，并在标题上添加 AUC 值。在这个例子中，AUC 值是 0.92，这表明 svm 分类器在电信用户流失数据集的分类中表现良好。

参考内容

• 对于对 ROC 的概念和术语感兴趣的人，您可以参考en.wikipedia.org/wiki/Receiver_operating_characteristic

使用 caret 包比较 ROC 曲线

在前面的章节中，我们介绍了许多分类方法；每种方法都有其自身的优缺点。然而，当涉及到如何选择最佳拟合模型的问题时，您需要比较来自不同预测模型生成的所有性能指标。为了使比较变得容易，caret包允许我们生成并比较模型的性能。在本例中，我们将使用caret包提供的函数来比较在同一数据集上训练的不同算法模型。

准备工作

在这里，我们将继续使用 telecom 数据集作为我们的输入数据源。

如何做...

执行以下步骤以生成每个拟合模型的 ROC 曲线：

1. 安装并加载库，pROC：

   > install.packages("pROC")
   > library("pROC")
   
   

2. 设置训练控制，使用 3 次重复的 10 折交叉验证：

   > control = trainControl(method = "repeatedcv",
   +                            number = 10,
   +                            repeats = 3,
   +                            classProbs = TRUE,
   +                            summaryFunction = twoClassSummary)
   
   

3. 然后，您可以使用glm在训练数据集上训练一个分类器：

   > glm.model= train(churn ~ .,
   +                     data = trainset,
   +                     method = "glm",
   +                     metric = "ROC",
   +                     trControl = control)
   
   

4. 此外，您还可以使用svm在训练数据集上训练一个分类器：

   > svm.model= train(churn ~ .,
   +                   data = trainset,
   +                   method = "svmRadial",
   +                   metric = "ROC",
   +                   trControl = control)
   
   

5. 为了查看rpart在训练数据上的表现，我们使用rpart函数：

   > rpart.model= train(churn ~ .,
   +                   data = trainset,
   +                   method = "rpart",
   +                   metric = "ROC",
   +                   trControl = control)
   
   

6. 您可以根据不同的训练模型分别进行预测：

   > glm.probs = predict(glm.model, testset[,! names(testset) %in% c("churn")], type = "prob")
   > svm.probs = predict(svm.model, testset[,! names(testset) %in% c("churn")], type = "prob")
   > rpart.probs = predict(rpart.model, testset[,! names(testset) %in% c("churn")], type = "prob")
   
   

7. 您可以为每个模型生成 ROC 曲线，并在同一图上绘制曲线：

   > glm.ROC = roc(response = testset[,c("churn")],
   +                predictor =glm.probs$yes,
   +                levels = levels(testset[,c("churn")]))
   > plot(glm.ROC, type="S", col="red") 
   
   Call:
   roc.default(response = testset[, c("churn")], predictor = glm.probs$yes,     levels = levels(testset[, c("churn")]))
   
   Data: glm.probs$yes in 141 controls (testset[, c("churn")] yes) > 877 cases (testset[, c("churn")] no).
   Area under the curve: 0.82
   
   > svm.ROC = roc(response = testset[,c("churn")],
   +                predictor =svm.probs$yes,
   +                levels = levels(testset[,c("churn")]))
   > plot(svm.ROC, add=TRUE, col="green") 
   
   Call:
   roc.default(response = testset[, c("churn")], predictor = svm.probs$yes,     levels = levels(testset[, c("churn")]))
   
   Data: svm.probs$yes in 141 controls (testset[, c("churn")] yes) > 877 cases (testset[, c("churn")] no).
   Area under the curve: 0.9233
   
   > rpart.ROC = roc(response = testset[,c("churn")],
   +                predictor =rpart.probs$yes,
   +                levels = levels(testset[,c("churn")]))
   > plot(rpart.ROC, add=TRUE, col="blue")
   
   Call:
   roc.default(response = testset[, c("churn")], predictor = rpart.probs$yes,     levels = levels(testset[, c("churn")]))
   
   Data: rpart.probs$yes in 141 controls (testset[, c("churn")] yes) > 877 cases (testset[, c("churn")] no).
   Area under the curve: 0.7581
   
   

   

   图 7：三个分类器性能的 ROC 曲线

它是如何工作的...

在这里，我们通过展示如何在一张图上绘制它们的 ROC 曲线来说明我们如何比较拟合模型。首先，我们通过在twoClassSummary中进行性能评估，使用 10 折交叉验证和 3 次重复来设置训练过程的控制。在设置好训练过程的控制后，我们接着在训练数据集上应用glm、svm和rpart算法来拟合分类模型。接下来，我们可以基于每个生成的模型进行预测并分别绘制 ROC 曲线。在生成的图中，我们发现由 svm 训练的模型具有最大的曲线下面积，为 0.9233（用绿色绘制），glm模型的 AUC 为 0.82，rpart模型的 AUC 为 0.7581。从图 7中可以看出，svm在此训练数据集上（无需调整）的所有拟合模型中表现最佳。

参见

• 我们使用另一个 ROC 可视化包pROC，它可以用来显示和分析 ROC 曲线。如果您想了解更多关于这个包的信息，请使用help函数：

  > help(package="pROC")
  
  

使用 caret 包测量模型之间的性能差异

在之前的步骤中，我们介绍了如何为每个生成的模型生成 ROC 曲线，并将曲线绘制在同一张图上。除了使用 ROC 曲线外，还可以使用重采样方法来生成每个拟合模型在 ROC、敏感性和特异性指标中的统计数据。因此，我们可以使用这些统计数据来比较每个模型之间的性能差异。在接下来的步骤中，我们将介绍如何使用caret包来测量拟合模型之间的性能差异。

准备工作

需要完成之前的步骤，将glm拟合模型、svm拟合模型和rpart拟合模型分别存储到glm.model、svm.model和rpart.model中。

如何做...

执行以下步骤来测量每个拟合模型之间的性能差异：

1. 重采样三个生成的模型：

   > cv.values = resamples(list(glm = glm.model, svm=svm.model, rpart = rpart.model))
   
   

2. 然后，您可以获得重采样结果的摘要：

   > summary(cv.values)
   
   Call:
   summary.resamples(object = cv.values)
   
   Models: glm, svm, rpart 
   Number of resamples: 30 
   
   ROC 
    Min. 1st Qu. Median   Mean 3rd Qu.   Max. NA's
   glm   0.7206  0.7847 0.8126 0.8116  0.8371 0.8877    0
   svm   0.8337  0.8673 0.8946 0.8929  0.9194 0.9458    0
   rpart 0.2802  0.7159 0.7413 0.6769  0.8105 0.8821    0
   
   Sens 
    Min. 1st Qu. Median   Mean 3rd Qu.   Max. NA's
   glm   0.08824  0.2000 0.2286 0.2194  0.2517 0.3529    0
   svm   0.44120  0.5368 0.5714 0.5866  0.6424 0.7143    0
   rpart 0.20590  0.3742 0.4706 0.4745  0.5929 0.6471    0
   
   Spec 
    Min. 1st Qu. Median   Mean 3rd Qu.   Max. NA's
   glm   0.9442  0.9608 0.9746 0.9701  0.9797 0.9949    0
   svm   0.9442  0.9646 0.9746 0.9740  0.9835 0.9949    0
   rpart 0.9492  0.9709 0.9797 0.9780  0.9848 0.9949    0
   
   

3. 使用dotplot在 ROC 指标中绘制重采样结果：

   > dotplot(cv.values, metric = "ROC")
   
   

   

   图 8：ROC 指标中重采样结果的点图

4. 此外，您还可以使用箱线图来绘制重采样结果：

   > bwplot(cv.values, layout = c(3, 1))
   
   

   

   图 9：重采样结果的箱线图

它是如何工作的...

在这个菜谱中，我们展示了如何使用重采样方法来衡量三个拟合模型之间的性能差异。首先，我们使用 resample 函数生成每个拟合模型（svm.model、glm.model 和 rpart.model）的统计数据。然后，我们可以使用 summary 函数获取这三个模型在 ROC、敏感性和特异性指标中的统计数据。接下来，我们可以在重采样结果上应用 dotplot 来查看 ROC 在不同模型之间的变化情况。最后，我们使用箱线图在重采样结果上展示不同模型在 ROC、敏感性和特异性指标上的箱线图。

参见

• 除了使用 dotplot 和 bwplot 来衡量性能差异外，还可以使用 densityplot、splom 和 xyplot 来可视化每个拟合模型在 ROC、敏感性和特异性指标中的性能差异。

第八章：集成学习

在本章中，我们将涵盖以下主题：

• 使用 Bagging 方法分类数据

• 使用 Bagging 方法进行交叉验证

• 使用 Boosting 方法分类数据

• 使用 Boosting 方法进行交叉验证

• 使用梯度提升分类数据

• 计算分类器的边缘

• 计算集成方法的误差演变

• 使用随机森林分类数据

• 估计不同分类器的预测误差

简介

集成学习是一种将不同学习器产生的结果组合成一种格式的方法，目的是产生更好的分类结果和回归结果。在前面的章节中，我们讨论了几种分类方法。这些方法采取不同的方法，但它们都有相同的目标，即找到最优的分类模型。然而，单个分类器可能是不完美的，它可能在某些类别中错误地分类数据。由于并非所有分类器都是不完美的，因此更好的方法是通过对结果进行平均投票。换句话说，如果我们对每个分类器相同的输入的预测结果进行平均，我们可能创建一个比使用单个方法更好的模型。

在集成学习中，Bagging、Boosting 和随机森林是三种最常用的方法：

• Bagging 是一种投票方法，它首先使用 Bootstrap 生成不同的训练集，然后使用该训练集来制作不同的基础学习器。Bagging 方法采用基础学习器的组合来做出更好的预测。

• Boosting 与 Bagging 方法类似。然而，使 Boosting 与众不同的地方在于它首先按顺序构建基础学习器，其中每个后续学习器都是为了构建前一个学习器的预测残差。通过创建互补学习者的手段，它使用先前学习者的错误来训练下一个基础学习器。

• 随机森林使用来自许多分类树的投票结果。这个想法很简单；单个分类树将使用单个输入向量获得单个分类结果。然而，随机森林生长了许多分类树，从单个输入中获得多个结果。因此，随机森林将使用所有决策树中的多数投票来分类数据或使用平均输出进行回归。

在下面的菜谱中，我们将讨论如何使用 Bagging 和 Boosting 来分类数据。然后我们可以进行交叉验证来估计每个分类器的错误率。除此之外，我们还将介绍使用边缘来衡量模型确定性的方法。接下来，我们将介绍类似于 Bagging 和 Boosting 方法的随机森林，并介绍如何训练模型来分类数据以及如何使用边缘来估计模型确定性。最后，我们将演示如何估计每个分类器的错误率，并使用错误率来比较不同分类器的性能。

使用 bagging 方法对数据进行分类

adabag包实现了提升和 bagging 方法。对于 bagging 方法，该包实现了 Breiman 的 Bagging 算法，该算法首先生成多个分类器的多个版本，然后获得一个聚合分类器。在这个菜谱中，我们将说明如何使用adabag中的 bagging 方法来使用电信churn数据集生成分类模型。

准备工作

在这个菜谱中，我们继续使用电信churn数据集作为 bagging 方法的输入数据源。对于那些尚未准备数据集的人，请参阅第五章, 分类（I） – 树、懒惰和概率分类，获取详细信息。

如何操作...

执行以下步骤以生成电信churn数据集的分类模型：

1. 首先，你需要安装并加载adabag包（安装adabag可能需要一段时间）：

   > install.packages("adabag")
   > library(adabag)
   
   

2. 接下来，你可以使用bagging函数来训练一个训练数据集（训练过程中结果可能会有所不同）：

   > set.seed(2)
   > churn.bagging = bagging(churn ~ ., data=trainset, mfinal=10)
   
   

3. 从 bagging 结果中获取变量重要性：

   > churn.bagging$importance
    international_plan number_customer_service_calls 
    10.4948380                    16.4260510 
    number_vmail_messages               total_day_calls 
    0.5319143                     0.3774190 
    total_day_charge             total_day_minutes 
    0.0000000                    28.7545042 
    total_eve_calls              total_eve_charge 
    0.1463585                     0.0000000 
    total_eve_minutes              total_intl_calls 
    14.2366754                     8.7733895 
    total_intl_charge            total_intl_minutes 
    0.0000000                     9.7838256 
    total_night_calls            total_night_charge 
    0.4349952                     0.0000000 
    total_night_minutes               voice_mail_plan 
    2.3379622                     7.7020671 
   
   

4. 在生成分类模型后，你可以使用测试数据集的预测结果：

   > churn.predbagging= predict.bagging(churn.bagging, newdata=testset)
   
   

5. 从预测结果中，你可以获得一个分类表：

   > churn.predbagging$confusion
    Observed Class
   Predicted Class yes  no
    no   35 866
    yes 106  11
   
   

6. 最后，你可以检索 bagging 结果的平均误差：

   > churn.predbagging$error
   [1] 0.0451866
   
   

如何工作...

Bagging 来源于 Bootstrap aggregating（自助聚合）的名字，这是一种稳定、准确且易于实现的数据分类和回归模型。Bagging 的定义如下：给定一个大小为n的训练数据集，bagging 执行 Bootstrap 抽样并生成m个新的训练集，Di，每个大小为n。最后，我们可以将m个 Bootstrap 样本拟合到m个模型中，并通过平均输出（对于回归）或投票（对于分类）来组合结果：

bagging 方法的示意图

使用 bagging 的优点是它是一种强大的学习方法，易于理解和实现。然而，这种技术的缺点是结果难以分析。

在这个菜谱中，我们使用adabag中的提升方法对电信客户流失数据进行分类。类似于前面章节中讨论的其他分类方法，你可以使用公式和训练数据集来训练一个提升分类器。此外，你可以在mfinal参数中将迭代次数设置为 10。一旦构建了分类模型，你可以检查每个属性的重要性。按重要性对属性进行排序显示，客户服务电话的数量在分类模型中起着至关重要的作用。

接下来，使用拟合好的模型，您可以应用predict.bagging函数来预测测试数据集的标签。因此，您可以使用测试数据集的标签和预测结果来生成分类表并获取平均误差，在这个例子中为 0.045。

还有更多...

除了adabag，ipred包还提供了一个用于分类树的袋装方法。在这里，我们演示如何使用ipred包的袋装方法来训练一个分类模型：

1. 首先，您需要安装并加载ipred包：

   > install.packages("ipred")
   > library(ipred)
   
   

2. 然后，您可以使用bagging方法来拟合分类方法：

   > churn.bagging = bagging(churn ~ ., data = trainset, coob = T)
   > churn.bagging
   
   Bagging classification trees with 25 bootstrap replications 
   
   Call: bagging.data.frame(formula = churn ~ ., data = trainset, coob = T)
   
   Out-of-bag estimate of misclassification error:  0.0605 
   
   

3. 获取错误分类的出袋估计：

   > mean(predict(churn.bagging) != trainset$churn)
   [1] 0.06047516
   
   

4. 然后，您可以使用predict函数来获取测试数据集的预测标签：

   > churn.prediction = predict(churn.bagging, newdata=testset, type="class")
   
   

5. 从测试数据集的标签和预测结果中获取分类表：

   > prediction.table = table(churn.prediction, testset$churn)
   
   churn.prediction yes  no
    no   31 869
    yes 110   8
   
   

使用袋装方法进行交叉验证

为了评估分类器的预测能力，您可以运行交叉验证方法来测试分类模型的鲁棒性。在这个菜谱中，我们将介绍如何使用bagging.cv来执行带有袋装方法的交叉验证。

准备工作

在这个菜谱中，我们继续使用 telecom churn数据集作为输入数据源，使用袋装方法进行 k 折交叉验证。

如何做...

通过使用袋装方法进行交叉验证来执行以下步骤以检索最小估计误差：

1. 首先，我们使用bagging.cv在训练数据集上进行 10 次迭代的 10 折分类：

   > churn.baggingcv = bagging.cv(churn ~ ., v=10, data=trainset, mfinal=10)
   
   

2. 然后，您可以从交叉验证结果中获取混淆矩阵：

   > churn.baggingcv$confusion
    Observed Class
   Predicted Class  yes   no
    no   100 1938
    yes  242   35
   
   

3. 最后，您可以从交叉验证结果中检索最小估计误差：

   > churn.baggingcv$error
   [1] 0.05831533
   
   

工作原理...

adabag包提供了一个函数来执行带有袋装或提升方法的 k 折验证。在这个例子中，我们使用bagging.cv来使用袋装方法进行 k 折交叉验证。我们首先通过指定v=10和mfinal=10进行 10 次迭代的 10 折交叉验证。请注意，由于迭代次数较多，这相当耗时。交叉验证过程完成后，我们可以从交叉验证结果中获得混淆矩阵和平均误差（在这种情况下为 0.058）。

参见

• 对于那些想要调整bagging.cv参数的人来说，请使用help函数查看bagging.cv文档：

  > help(bagging.cv)
  
  

使用提升方法进行数据分类

与袋装法类似，提升法从简单或弱分类器开始，通过重新加权误分类样本逐渐改进它。因此，新的分类器可以从前一个分类器中学习。adabag包提供了AdaBoost.M1和SAMME算法的实现。因此，可以使用adabag中的提升方法进行集成学习。在本菜谱中，我们将使用adabag中的提升方法对电信churn数据集进行分类。

准备工作

在本菜谱中，我们将继续使用电信churn数据集作为输入数据源，使用提升法进行分类。此外，在开始菜谱之前，您需要将adabag包加载到 R 中。

如何操作...

执行以下步骤，使用提升法对电信churn数据集进行分类：

1. 您可以使用adabag包中的提升函数来训练分类模型：

   > set.seed(2)
   > churn.boost = boosting(churn ~.,data=trainset,mfinal=10, coeflearn="Freund", boos=FALSE , control=rpart.control(maxdepth=3))
   
   

2. 然后，您可以根据提升模型和测试数据集进行预测：

   > churn.boost.pred = predict.boosting(churn.boost,newdata=testset)
   
   

3. 接下来，您可以从预测结果中检索分类表：

   > churn.boost.pred$confusion
    Observed Class
   Predicted Class yes  no
    no   41 858
    yes 100  19
   
   

4. 最后，您可以从预测结果中获得平均误差：

   > churn.boost.pred$error
   [1] 0.0589391
   
   

它是如何工作的...

提升法的思想是将弱学习器（例如单个决策树）提升为强学习器。假设我们的训练数据集中有n个点，我们可以为每个点分配一个权重，Wi（0 <= i < n）。然后，在迭代学习过程中（我们假设迭代次数为m），我们可以根据每次迭代的分类结果重新加权每个点。如果点被正确分类，我们应该减少权重。否则，我们增加该点的权重。当迭代过程完成后，我们可以获得m个拟合模型，fi（0 <= i < n）。最后，我们可以通过每个树的预测的加权平均获得最终预测，其中权重 b 基于每个树的质量：

提升法的示意图

无论是袋装法（bagging）还是提升法（boosting），它们都是集成方法，将每个单个学习器的预测能力结合成一个强大的学习器。袋装法和提升法之间的区别在于，袋装法结合了独立的模型，而提升法通过迭代过程，通过后续模型预测先前模型的错误来减少错误。

在这个菜谱中，我们展示了如何在提升方法内拟合分类模型。类似于袋装法，必须指定用于训练分类模型的公式和训练数据集。此外，还可以指定参数，例如迭代次数（mfinal）、权重更新系数（coeflearn）、每个观测值的权重（boos）以及rpart（单个决策树）的控制。在这个菜谱中，我们将迭代次数设置为 10，使用Freund（Adaboost.M1 算法实现的方法）作为coeflearn，将boos设置为 false，并将rpart配置的最大深度设置为3。

我们使用提升方法拟合分类模型，并将其保存在churn.boost中。然后，我们可以使用prediction函数获取预测标签。此外，我们可以使用table函数根据预测标签和测试数据集标签检索分类表。最后，我们可以获取预测结果的平均误差。

更多...

除了在adabag包中使用提升函数外，还可以使用caret包通过提升方法进行分类：

1. 首先，加载mboost和pROC包：

   > library(mboost)
   > install.packages("pROC")
   > library(pROC)
   
   

2. 然后，我们可以使用trainControl函数设置训练控制，并使用train函数使用 adaboost 训练分类模型：

   > set.seed(2)
   > ctrl = trainControl(method = "repeatedcv", repeats = 1, classProbs = TRUE, summaryFunction = twoClassSummary)
   > ada.train = train(churn ~ ., data = trainset, method = "ada", metric = "ROC", trControl = ctrl)
   
   

3. 使用summary函数获取分类模型的详细信息：

   > ada.train$result
    nu maxdepth iter       ROC      Sens        Spec      ROCSD     SensSD      SpecSD
   1 0.1        1   50 0.8571988 0.9152941 0.012662155 0.03448418 0.04430519 0.007251045
   4 0.1        2   50 0.8905514 0.7138655 0.006083679 0.03538445 0.10089887 0.006236741
   7 0.1        3   50 0.9056456 0.4036134 0.007093780 0.03934631 0.09406015 0.006407402
   2 0.1        1  100 0.8550789 0.8918487 0.015705276 0.03434382 0.06190546 0.006503191
   5 0.1        2  100 0.8907720 0.6609244 0.009626724 0.03788941 0.11403364 0.006940001
   8 0.1        3  100 0.9077750 0.3832773 0.005576065 0.03601187 0.09630026 0.003738978
   3 0.1        1  150 0.8571743 0.8714286 0.016720505 0.03481526 0.06198773 0.006767313
   6 0.1        2  150 0.8929524 0.6171429 0.011654617 0.03638272 0.11383803 0.006777465
   9 0.1        3  150 0.9093921 0.3743697 0.007093780 0.03258220 0.09504202 0.005446136
   
   

4. 使用plot函数在不同迭代中绘制 ROC 曲线：

   > plot(ada.train)
   
   

   

   重复交叉验证图

5. 最后，我们可以使用predict函数进行预测，并查看分类表：

   > ada.predict = predict(ada.train, testset, "prob")
   > ada.predict.result = ifelse(ada.predict[1] > 0.5, "yes", "no")
   
   > table(testset$churn, ada.predict.result)
    ada.predict.result
    no yes
    yes  40 101
    no  872   5
   
   

使用提升方法进行交叉验证

类似于bagging函数，adabag为提升方法提供了一个交叉验证函数，名为boosting.cv。在这个菜谱中，我们将演示如何使用adabag包中的boosting.cv进行交叉验证。

准备工作

在这个菜谱中，我们继续使用电信churn数据集作为输入数据源，使用boosting方法进行 k 折交叉验证。

如何做...

通过使用boosting方法进行交叉验证，执行以下步骤以检索最小估计误差：

1. 首先，你可以使用boosting.cv对训练数据集进行交叉验证：

   > churn.boostcv = boosting.cv(churn ~ ., v=10, data=trainset, mfinal=5,control=rpart.control(cp=0.01))
   
   

2. 然后，你可以从提升结果中获得混淆矩阵：

   > churn.boostcv$confusion
    Observed Class
   Predicted Class  yes   no
    no   119 1940
    yes  223   33
   
   

3. 最后，你可以检索提升方法的平均误差：

   > churn.boostcv$error
   [1] 0.06565875
   
   

它是如何工作的...

类似于bagging.cv，我们可以使用boosting.cv在提升方法中进行交叉验证。如果将v设置为10，将mfinal设置为5，则提升方法将执行 10 次交叉验证，每次迭代 5 次。此外，还可以在参数内设置rpart拟合的控制。在这个例子中，我们可以将复杂性参数设置为 0.01。一旦训练完成，将获得提升结果的混淆矩阵和平均误差。

参见

• 对于需要更多关于调整boosting.cv参数信息的人，请使用help函数查看boosting.cv文档：

  > help(boosting.cv)
  
  

使用梯度提升法分类数据

梯度提升集成弱学习器并创建一个与损失函数的负梯度最大相关的新基学习器。可以将此方法应用于回归或分类问题，并且在不同数据集上表现良好。在本菜谱中，我们将介绍如何使用gbm对电信churn数据集进行分类。

准备工作

在这个菜谱中，我们继续使用电信churn数据集作为bagging方法的输入数据源。对于那些尚未准备数据集的人，请参阅第五章，分类（I）-树、懒惰和概率性，以获取详细信息。

如何做

执行以下步骤以使用梯度提升法计算和分类数据：

1. 首先，安装并加载gbm包：

   > install.packages("gbm")
   > library(gbm)
   
   

2. gbm函数只使用从0到1的响应；因此，你应该将是/否响应转换为数值响应（0/1）：

   > trainset$churn = ifelse(trainset$churn == "yes", 1, 0)
   
   

3. 接下来，你可以使用gbm函数来训练一个训练数据集：

   > set.seed(2)
   > churn.gbm = gbm(formula = churn ~ .,distribution = "bernoulli",data = trainset,n.trees = 1000,interaction.depth = 7,shrinkage = 0.01, cv.folds=3)
   
   

4. 然后，你可以从拟合模型中获取摘要信息：

   > summary(churn.gbm)
    var    rel.inf
   total_day_minutes          total_day_minutes 28.1217147
   total_eve_minutes                total_eve_minutes 16.8097151
   number_customer_service_calls number_customer_service_calls 12.7894464
   total_intl_minutes             total_intl_minutes  9.4515822
   total_intl_calls                   total_intl_calls  8.1379826
   international_plan               international_plan  8.0703900
   total_night_minutes             total_night_minutes  4.0805153
   number_vmail_messages         number_vmail_messages  3.9173515
   voice_mail_plan                  voice_mail_plan  2.5501480
   total_night_calls              total_night_calls  2.1357970
   total_day_calls                     total_day_calls  1.7367888
   total_eve_calls                     total_eve_calls  1.4398047
   total_eve_charge                 total_eve_charge  0.5457486
   total_night_charge              total_night_charge  0.2130152
   total_day_charge                total_day_charge  0.0000000
   total_intl_charge                 total_intl_charge  0.0000000
   
   

   

   拟合模型的相对影响图

5. 你可以使用交叉验证获得最佳迭代次数：

   > churn.iter = gbm.perf(churn.gbm,method="cv")
   
   

   

   性能测量图

6. 然后，你可以检索伯努利损失函数返回的对数中的奇数值：

   > churn.predict = predict(churn.gbm, testset, n.trees = churn.iter)
   > str(churn.predict)
    num [1:1018] -3.31 -2.91 -3.16 -3.47 -3.48 ...
   
   

7. 接下来，你可以绘制 ROC 曲线并获取具有最大准确性的最佳截止值：

   > churn.roc = roc(testset$churn, churn.predict)
   > plot(churn.roc)
   Call:
   roc.default(response = testset$churn, predictor = churn.predict)
   Data: churn.predict in 141 controls (testset$churn yes) > 877 cases (testset$churn no).
   Area under the curve: 0.9393
   
   

   

   拟合模型的 ROC 曲线

8. 你可以使用coords函数检索最佳截止值并使用此截止值来获取预测标签：

   > coords(churn.roc, "best")
    threshold specificity sensitivity 
    -0.9495258   0.8723404   0.9703535 
   > churn.predict.class = ifelse(churn.predict > coords(churn.roc, "best")["threshold"], "yes", "no")
   
   

9. 最后，你可以从预测结果中获得分类表：

   > table( testset$churn,churn.predict.class)
    churn.predict.class
    no yes
    yes  18 123
    no  851  26
   
   

它是如何工作的...

梯度提升算法涉及的过程，首先计算每个分区的残差偏差，然后确定每个阶段的最佳数据分区。接下来，连续的模型将拟合前一阶段的残差并构建一个新的模型以减少残差方差（一个误差）。残差方差的减少遵循功能梯度下降技术，通过沿着其导数下降来最小化残差方差，如下所示：

梯度下降法

在这个菜谱中，我们使用gbm中的梯度提升方法对电信客户流失数据集进行分类。为了开始分类，我们首先安装并加载gbm包。然后，我们使用gbm函数来训练分类模型。在这里，由于我们的预测目标是churn属性，这是一个二元结果，因此我们在distribution参数中将分布设置为bernoulli。此外，我们在n.tree参数中将 1000 棵树设置为拟合，在interaction.depth中将变量交互的最大深度设置为7，在shrinkage中将步长减少的学习率设置为 0.01，在cv.folds中将交叉验证的次数设置为3。模型拟合后，我们可以使用summary函数从表中获取每个变量的相对影响信息，并从图表中获取。相对影响显示了每个变量在平方误差和中的减少量。在这里，我们可以找到total_day_minutes是减少损失函数的最有影响力的一个变量。

接下来，我们使用gbm.perf函数找到最佳迭代次数。在这里，我们通过指定method参数为cv来进行交叉验证，以估计最佳迭代次数。该函数进一步生成两个图表，其中黑色线条表示训练误差，绿色线条表示验证误差。这里的误差测量是一个bernoulli分布，我们已经在训练阶段定义了它。图表上的蓝色虚线显示了最佳迭代的位置。

然后，我们使用predict函数从伯努利损失函数返回的每个测试案例中获取对数概率值。为了获得最佳的预测结果，可以将n.trees参数设置为最佳迭代次数。然而，由于返回值是对数概率值的奇数，我们仍然需要确定最佳的截止点来确定标签。因此，我们使用roc函数生成 ROC 曲线，并获取具有最大准确率的截止点。

最后，我们可以使用coords函数检索最佳截止阈值，并使用ifelse函数从对数概率值中确定类别标签。现在，我们可以使用table函数生成分类表，以查看分类模型的准确性。

还有更多...

除了在gbm包中使用提升函数外，还可以使用mboost包通过梯度提升方法进行分类：

1. 首先，安装并加载mboost包：

   > install.packages("mboost")
   > library(mboost)
   
   

2. mboost函数仅使用数值响应；因此，应将是/否响应转换为数值响应（0/1）：

   > trainset$churn = ifelse(trainset$churn == "yes", 1, 0)
   
   

3. 此外，还应删除非数值属性，例如voice_mail_plan和international_plan：

   > trainset$voice_mail_plan = NULL
   > trainset$international_plan = NULL
   
   

4. 然后，我们可以使用mboost来训练分类模型：

   > churn.mboost = mboost(churn ~ ., data=trainset,  control = boost_control(mstop = 10))
   
   

5. 使用summary函数获取分类模型的详细信息：

   > summary(churn.mboost)
   
    Model-based Boosting
   
   Call:
   mboost(formula = churn ~ ., data = trainset, control = boost_control(mstop = 10))
   
    Squared Error (Regression) 
   
   Loss function: (y - f)² 
   
   Number of boosting iterations: mstop = 10 
   Step size:  0.1 
   Offset:  1.147732 
   Number of baselearners:  14 
   
   Selection frequencies:
    bbs(total_day_minutes) bbs(number_customer_service_calls) 
    0.6                                0.4 
   
   

6. 最后，使用plot函数绘制每个属性的局部贡献图：

   > par(mfrow=c(1,2))
   > plot(churn.mboost)
   
   

   

   重要属性的局部贡献图

计算分类器的边缘

边缘是分类确定性的度量。此方法计算正确类别的支持与错误类别的最大支持之间的差异。在本配方中，我们将演示如何计算生成的分类器的边缘。

准备工作

您需要完成前面的配方，将拟合的 bagging 模型存储在变量 churn.bagging 和 churn.predbagging 中。同时，将拟合的提升分类器放入 churn.boost 和 churn.boost.pred 中。

如何操作...

执行以下步骤以计算每个集成学习者的边缘：

1. 首先，使用 margins 函数计算提升分类器的边缘：

   > boost.margins = margins(churn.boost, trainset)
   > boost.pred.margins = margins(churn.boost.pred, testset)
   
   

2. 然后，您可以使用 plot 函数绘制提升分类器的边缘累积分布图：

   > plot(sort(boost.margins[[1]]), (1:length(boost.margins[[1]]))/length(boost.margins[[1]]), type="l",xlim=c(-1,1),main="Boosting: Margin cumulative distribution graph", xlab="margin", ylab="% observations", col = "blue")
   > lines(sort(boost.pred.margins[[1]]), (1:length(boost.pred.margins[[1]]))/length(boost.pred.margins[[1]]), type="l", col = "green")
   > abline(v=0, col="red",lty=2)
   
   

   

   使用提升方法的边缘累积分布图

3. 然后，您可以计算负边缘匹配训练错误的百分比和负边缘匹配测试错误的百分比：

   > boosting.training.margin = table(boost.margins[[1]] > 0)
   > boosting.negative.training = as.numeric(boosting.training.margin[1]/boosting.training.margin[2])
   > boosting.negative.training
    [1] 0.06387868
   
   > boosting.testing.margin = table(boost.pred.margins[[1]] > 0)
   > boosting.negative.testing = as.numeric(boosting.testing.margin[1]/boosting.testing.margin[2])
   > boosting.negative.testing
   [1] 0.06263048
   
   

4. 此外，您还可以计算 bagging 分类器的边缘。您可能会看到显示 "no non-missing argument to min" 的警告信息。该信息仅表示 min/max 函数应用于长度为 0 的参数的数值：

   > bagging.margins = margins(churn.bagging, trainset)
   > bagging.pred.margins = margins(churn.predbagging, testset)
   
   

5. 然后，您可以使用 plot 函数绘制 bagging 分类器的边缘累积分布图：

   > plot(sort(bagging.margins[[1]]), (1:length(bagging.margins[[1]]))/length(bagging.margins[[1]]), type="l",xlim=c(-1,1),main="Bagging: Margin cumulative distribution graph", xlab="margin", ylab="% observations", col = "blue")
   
   > lines(sort(bagging.pred.margins[[1]]), (1:length(bagging.pred.margins[[1]]))/length(bagging.pred.margins[[1]]), type="l", col = "green")
   > abline(v=0, col="red",lty=2)
   
   

   

   bagging 方法的边缘累积分布图

6. 最后，然后您可以计算负边缘匹配训练错误的百分比和负边缘匹配测试错误的百分比：

   > bagging.training.margin = table(bagging.margins[[1]] > 0)
   > bagging.negative.training = as.numeric(bagging.training.margin[1]/bagging.training.margin[2])
   > bagging.negative.training
   [1] 0.1733401
   
   > bagging.testing.margin = table(bagging.pred.margins[[1]] > 0)
   > bagging.negative.testing = as.numeric(bagging.testing.margin[1]/bagging.testing.margin[2])
   > bagging.negative.testing
   [1] 0.04303279
   
   

工作原理...

边缘是分类确定性的度量；它是通过正确类别的支持与错误类别的最大支持来计算的。边缘的公式可以表示为：

在这里，xi 样本的边缘等于正确分类样本的支持（c 表示正确类别）减去被分类到类别 j 的样本的最大支持（其中 j≠c 且 j=1…k）。因此，正确分类的示例将具有正边缘，错误分类的示例将具有负边缘。如果边缘值接近 1，则表示正确分类的示例具有高度的置信度。另一方面，不确定分类的示例将具有较小的边缘。

margins函数计算 AdaBoost.M1、AdaBoost-SAMME 或袋装分类器的边缘，并返回一个边缘向量。为了可视化边缘分布，可以使用边缘累积分布图。在这些图中，x 轴表示边缘，y 轴表示边缘小于或等于 x 轴边缘值的观测值的百分比。如果每个观测值都被正确分类，则图表将在边缘等于 1 的位置显示一条垂直线（其中边缘=1）。

对于提升分类器的边缘累积分布图，我们可以看到图上有两条线，其中绿色线表示测试数据集的边缘，蓝色线表示训练集的边缘。图显示大约 6.39%的负边缘与训练错误匹配，6.26%的负边缘与测试错误匹配。另一方面，我们可以在袋装分类器的边缘累积分布图中找到，17.33%的负边缘与训练错误匹配，4.3%的负边缘与测试错误匹配。通常，匹配训练错误的负边缘百分比应该接近匹配测试错误的负边缘百分比。因此，我们应该检查为什么匹配训练错误的负边缘百分比远高于匹配测试错误的负边缘百分比。

参考内容

• 如果您对边缘分布图有更多细节感兴趣，请参阅以下来源：Kuncheva LI (2004), Combining Pattern Classifiers: Methods and Algorithms, John Wiley & Sons。

计算集成方法的错误演化

adabag包提供了一个errorevol函数，用户可以根据迭代次数估计集成方法错误。在这个配方中，我们将演示如何使用errorevol来展示每个集成分类器错误的演化。

准备工作

您需要完成上一个配方，将拟合的袋装模型存储在变量churn.bagging中。同时，将拟合的提升分类器放入churn.boost。

如何做...

执行以下步骤来计算每个集成学习者的错误演化：

1. 首先，使用errorevol函数来计算提升分类器的错误演化：

   > boosting.evol.train = errorevol(churn.boost, trainset)
   > boosting.evol.test = errorevol(churn.boost, testset)
   > plot(boosting.evol.test$error, type = "l", ylim = c(0, 1),
   +       main = "Boosting error versus number of trees", xlab = "Iterations",
   +       ylab = "Error", col = "red", lwd = 2)
   > lines(boosting.evol.train$error, cex = .5, col = "blue", lty = 2, lwd = 2)
   > legend("topright", c("test", "train"), col = c("red", "blue"), lty = 1:2, lwd = 2)
   
   

   

   提升错误与树的数量对比

2. 接下来，使用errorevol函数来计算袋装分类器的错误演化：

   > bagging.evol.train = errorevol(churn.bagging, trainset)
   > bagging.evol.test = errorevol(churn.bagging, testset)
   > plot(bagging.evol.test$error, type = "l", ylim = c(0, 1),
   +       main = "Bagging error versus number of trees", xlab = "Iterations",
   +       ylab = "Error", col = "red", lwd = 2)
   > lines(bagging.evol.train$error, cex = .5, col = "blue", lty = 2, lwd = 2)
   > legend("topright", c("test", "train"), col = c("red", "blue"), lty = 1:2, lwd = 2)
   
   

   

   袋装错误与树的数量对比

它是如何工作的...

errorest函数计算 AdaBoost.M1、AdaBoost-SAMME 或袋装分类器的错误演化，并返回一个错误演化的向量。在这个配方中，我们使用提升和袋装模型生成错误演化向量，并绘制错误与树的数量对比图。

生成的图表揭示了每次迭代的错误率。错误率的变化趋势可以帮助衡量错误减少的速度，同时迭代次数增加。此外，图表可能显示模型是否过拟合。

相关内容

• 如果集成模型过拟合，你可以使用predict.bagging和predict.boosting函数来修剪集成模型。更多信息，请使用帮助函数参考predict.bagging和predict.boosting：

  > help(predict.bagging)
  > help(predict.boosting)
  
  

使用随机森林进行数据分类

随机森林是另一种有用的集成学习方法，在训练过程中会生成多个决策树。每个决策树将输出其对应的输入预测结果。森林将使用投票机制来选择得票最多的类别作为预测结果。在本教程中，我们将说明如何使用randomForest包进行数据分类。

准备工作

在本教程中，我们将继续使用电信churn数据集作为输入数据源，使用随机森林方法进行分类。

如何操作...

执行以下步骤以使用随机森林进行数据分类：

1. 首先，你必须安装并加载randomForest包；

   > install.packages("randomForest")
   > library(randomForest)
   
   

2. 然后，你可以使用训练集拟合随机森林分类器：

   > churn.rf = randomForest(churn ~ ., data = trainset, importance = T)
   > churn.rf
   
   Call:
    randomForest(formula = churn ~ ., data = trainset, importance = T) 
    Type of random forest: classification
    Number of trees: 500
   No. of variables tried at each split: 4
   
    OOB estimate of  error rate: 4.88%
   Confusion matrix:
    yes   no class.error
   yes 247   95 0.277777778
   no   18 1955 0.009123163
   
   

3. 接下来，根据拟合模型和测试数据集进行预测：

   > churn.prediction = predict(churn.rf, testset)
   
   

4. 与其他分类方法类似，你可以获得分类表：

   > table(churn.prediction, testset$churn)
   
   churn.prediction yes  no
    yes 110   7
    no   31 870
   
   

5. 你可以使用plot函数绘制森林对象的均方误差：

   > plot(churn.rf)
   
   

   

   随机森林的均方误差

6. 然后，你可以检查拟合分类器中每个属性的重要性：

   > importance(churn.rf)
    yes         no
   international_plan            66.55206691 56.5100647
   voice_mail_plan               19.98337191 15.2354970
   number_vmail_messages         21.02976166 14.0707195
   total_day_minutes             28.05190188 27.7570444
   
   

7. 接下来，你可以使用varImpPlot函数获取变量重要性的绘图：

   > varImpPlot(churn.rf)
   
   

   

   变量重要性的可视化

8. 你还可以使用margin函数来计算边缘并绘制边缘累积分布图：

   > margins.rf=margin(churn.rf,trainset)
   > plot(margins.rf)
   
   

   

   随机森林方法的边缘累积分布图

9. 此外，你可以使用直方图来可视化随机森林的边缘分布：

   > hist(margins.rf,main="Margins of Random Forest for churn dataset")
   
   

   

   边缘分布的直方图

10. 你还可以使用boxplot通过类别可视化随机森林的边缘：

    > boxplot(margins.rf~trainset$churn, main="Margins of Random Forest for churn dataset by class")
    
    

    

    随机森林的类别边缘

工作原理...

随机森林的目的是将弱学习器（例如单个决策树）集成到强学习器中。开发随机森林的过程与袋装法非常相似，假设我们有一个包含N个样本和M个特征的训练集。这个过程首先进行自助采样，随机抽取N个案例，作为每个单个决策树的训练数据集。接下来，在每个节点中，过程首先随机选择m个变量（其中m << M），然后找到在 m 个变量中提供最佳分割的预测变量。接下来，过程在不剪枝的情况下生长完整的树。最后，我们可以从每棵树中获得示例的预测结果。因此，我们可以通过对输出取平均或加权平均（对于回归）或取多数投票（对于分类）来获得预测结果：

随机森林使用两个参数：ntree（树的数量）和mtry（用于找到最佳特征的特性数量），而袋装法只使用 ntree 作为参数。因此，如果我们设置 mtry 等于训练数据集中的特性数量，那么随机森林就等同于袋装法。

随机森林的主要优点是它易于计算，可以高效地处理数据，并且对缺失或不平衡数据具有容错性。随机森林的主要缺点是它不能预测超出训练数据集范围的值。此外，它容易对噪声数据进行过拟合。

在这个菜谱中，我们使用来自randomForest包的随机森林方法来拟合分类模型。首先，我们在 R 会话中安装并加载randomForest。然后，我们使用随机森林方法训练分类模型。我们设置importance = T，这将确保评估预测变量的重要性。

与袋装法和提升法类似，一旦模型拟合完成，就可以在测试数据集上使用拟合的模型进行预测，并且还可以获得分类表。

为了评估每个属性的重要性，randomForest包提供了重要性和varImpPlot函数，既可以列出拟合模型中每个属性的重要性，也可以使用平均降低准确度或平均降低gini来可视化重要性。

与包含计算袋装法和提升法边缘方法的adabag类似，randomForest提供了margin函数来计算森林对象的边缘。使用plot、hist和boxplot函数，你可以从不同方面可视化正确分类观察值的比例。

还有更多...

除了randomForest包之外，party包还提供了一个随机森林的实现。在以下步骤中，我们将展示如何在party包中使用cforest函数来进行分类：

1. 首先，安装并加载party包：

   > install.packages("party")
   > library(party)
   
   

2. 然后，你可以使用cforest函数来拟合分类模型：

   > churn.cforest = cforest(churn ~ ., data = trainset, controls=cforest_unbiased(ntree=1000, mtry=5))
   > churn.cforest
   
    Random Forest using Conditional Inference Trees
   
   Number of trees:  1000 
   
   Response:  churn 
   Inputs:  international_plan, voice_mail_plan, number_vmail_messages, total_day_minutes, total_day_calls, total_day_charge, total_eve_minutes, total_eve_calls, total_eve_charge, total_night_minutes, total_night_calls, total_night_charge, total_intl_minutes, total_intl_calls, total_intl_charge, number_customer_service_calls 
   Number of observations:  2315 
   
   

3. 你可以根据构建的模型和测试数据集进行预测：

   > churn.cforest.prediction = predict(churn.cforest, testset, OOB=TRUE, type = "response")
   
   

4. 最后，从预测标签和测试数据集的标签中获取分类表：

   > table(churn.cforest.prediction, testset$churn)
   
   churn.cforest.prediction yes  no
    yes  91   3
    no   50 874
   
   

估计不同分类器的预测误差

在本章开头，我们讨论了为什么使用集成学习以及它如何与仅使用单个分类器相比提高预测性能。我们现在通过比较每种方法的性能来验证集成模型是否比单个决策树表现更好。为了比较不同的分类器，我们可以对每种分类方法进行 10 折交叉验证，使用ipred包中的erroreset来估计测试误差。

准备工作

在这个菜谱中，我们将继续使用电信churn数据集作为输入数据源来估计不同分类器的预测误差。

如何操作...

执行以下步骤来估计每种分类方法的预测误差：

1. 你可以估计袋装模型的错误率：

   > churn.bagging= errorest(churn ~ ., data = trainset, model = bagging)
   > churn.bagging
   
   Call:
   errorest.data.frame(formula = churn ~ ., data = trainset, model = bagging)
   
    10-fold cross-validation estimator of misclassification error 
   
   Misclassification error:  0.0583 
   
   

2. 然后，你可以估计提升方法的错误率：

   > install.packages("ada")
   > library(ada)
   > churn.boosting= errorest(churn ~ ., data = trainset, model = ada)
   > churn.boosting
   
   Call:
   errorest.data.frame(formula = churn ~ ., data = trainset, model = ada)
   
    10-fold cross-validation estimator of misclassification error 
   
   Misclassification error:  0.0475 
   
   

3. 接下来，估计随机森林模型的错误率：

   > churn.rf= errorest(churn ~ ., data = trainset, model = randomForest)
   > churn.rf
   
   Call:
   errorest.data.frame(formula = churn ~ ., data = trainset, model = randomForest)
   
    10-fold cross-validation estimator of misclassification error 
   
   Misclassification error:  0.051 
   
   

4. 最后，使用churn.predict创建一个预测函数，然后使用该函数估计单个决策树的错误率：

   > churn.predict = function(object, newdata) {predict(object, newdata = newdata, type = "class")}
   > churn.tree= errorest(churn ~ ., data = trainset, model = rpart,predict = churn.predict)
   > churn.tree
   
   Call:
   errorest.data.frame(formula = churn ~ ., data = trainset, model = rpart, 
    predict = churn.predict)
   
    10-fold cross-validation estimator of misclassification error 
   
   Misclassification error:  0.0674 
   
   

工作原理...

在这个菜谱中，我们使用ipred包中的errorest函数来估计四种不同分类器的错误率。我们比较了提升、袋装和随机森林方法，以及单个决策树分类器。errorest函数对每个分类器进行 10 折交叉验证，并计算误分类误差。四个选择模型的估计结果显示，提升方法具有最低的错误率（0.0475）。随机森林方法具有第二低的错误率（0.051），而袋装方法的错误率为 0.0583。单个决策树分类器rpart在四种方法中表现最差，错误率等于 0.0674。这些结果表明，提升、袋装和随机森林这三种集成学习方法都优于单个决策树分类器。

相关内容

• 在这个菜谱中，我们提到了ada包，它包含一个执行随机提升的方法。对那些对这个包感兴趣的人，请参阅：Friedman 等（2000）的《加性逻辑回归：提升的统计视角》。

第九章。聚类

在本章中，我们将涵盖以下主题：

• 使用层次聚类聚类数据

• 将树木切割成簇

• 使用 k-means 方法聚类数据

• 绘制双变量簇图

• 比较聚类方法

• 从聚类中提取轮廓信息

• 获取 k-means 的最佳簇

• 使用基于密度的方法聚类数据

• 使用基于模型的方法聚类数据

• 可视化差异矩阵

• 外部验证簇

简介

聚类是一种用于将相似对象（在距离上接近）分组到同一组（簇）中的技术。与前面章节中涵盖的监督学习方法（例如，分类和回归）不同，聚类分析不使用任何标签信息，而是简单地使用数据特征之间的相似性将它们分组到簇中。

聚类可以广泛应用于商业分析。例如，营销部门可以使用聚类根据个人属性对客户进行细分。因此，可以设计针对不同类型客户的差异化营销活动。

最常见的四种聚类方法为层次聚类、k-means 聚类、基于模型聚类和基于密度的聚类：

• 层次聚类：它创建簇的层次结构，并以树状图的形式展示。这种方法在开始时不需要指定簇的数量。

• k-means 聚类：也称为平面聚类。与层次聚类不同，它不创建簇的层次结构，并且需要一个输入的簇数量。然而，它的性能比层次聚类快。

• 基于模型聚类：层次聚类和 k-means 聚类都使用启发式方法构建簇，并且不依赖于正式模型。基于模型聚类假设数据模型，并应用 EM 算法来找到最可能的模型组件和簇的数量。

• 基于密度的聚类：它根据密度测量构建簇。这种方法中的簇比数据集的其余部分具有更高的密度。

在下面的菜谱中，我们将讨论如何使用这四种聚类技术来聚类数据。我们讨论了如何使用簇内平方和、平均轮廓宽度和外部真实情况来内部验证簇。

使用层次聚类聚类数据

层次聚类采用聚合或划分方法来构建簇的层次结构。无论采用哪种方法，两者首先都使用距离相似度度量来合并或分割簇。递归过程持续进行，直到只剩下一个簇或无法再分割簇。最终，我们可以使用树状图来表示簇的层次结构。在这个菜谱中，我们将演示如何使用层次聚类对客户进行聚类。

准备工作

在这个菜谱中，我们将对客户数据进行层次聚类，这涉及到将客户分割成不同的组。您可以从这个 GitHub 页面下载数据：github.com/ywchiu/ml_R_cookbook/tree/master/CH9。

如何做...

执行以下步骤将客户数据聚类到簇的层次结构中：

1. 首先，您需要从customer.csv加载数据并将其保存到customer中：

   > customer= read.csv('customer.csv', header=TRUE)
   > head(customer)
    ID Visit.Time Average.Expense Sex Age
   1  1          3             5.7   0  10
   2  2          5            14.5   0  27
   3  3         16            33.5   0  32
   4  4          5            15.9   0  30
   5  5         16            24.9   0  23
   6  6          3            12.0   0  15
   
   

2. 然后，您可以检查数据集的结构：

   > str(customer)
   'data.frame':  60 obs. of  5 variables:
    $ ID             : int  1 2 3 4 5 6 7 8 9 10 ...
    $ Visit.Time     : int  3 5 16 5 16 3 12 14 6 3 ...
    $ Average.Expense: num  5.7 14.5 33.5 15.9 24.9 12 28.5 18.8 23.8 5.3 ...
    $ Sex            : int  0 0 0 0 0 0 0 0 0 0 ...
    $ Age            : int  10 27 32 30 23 15 33 27 16 11 ...
   
   

3. 接下来，您应该将客户数据归一化到相同的尺度：

   > customer = scale(customer[,-1])
   
   

4. 您可以使用聚合层次聚类对客户数据进行聚类：

   > hc = hclust(dist(customer, method="euclidean"), method="ward.D2")
   > hc
   
   Call:
   hclust(d = dist(customer, method = "euclidean"), method = "ward.D2")
   
   Cluster method   : ward.D2 
   Distance         : euclidean 
   Number of objects: 60
   
   

5. 最后，您可以使用plot函数绘制树状图：

   > plot(hc, hang = -0.01, cex = 0.7)
   
   

   

   使用"ward.D2"的层次聚类树状图

6. 此外，您可以使用单方法进行层次聚类，并查看生成的树状图与之前的区别：

   > hc2 = hclust(dist(customer), method="single")
   > plot(hc2, hang = -0.01, cex = 0.7)
   
   

   

   使用"single"的层次聚类树状图

如何工作...

层次聚类是一种聚类技术，它试图通过迭代构建簇的层次结构。通常，有两种方法来构建层次簇：

• 聚合层次聚类：这是一种自下而上的方法。每个观测值最初都在自己的簇中。然后我们可以计算每个簇之间的相似度（或距离），并在每次迭代中将两个最相似的簇合并，直到只剩下一个簇。

• 划分层次聚类：这是一种自上而下的方法。所有观测值最初都在一个簇中，然后我们递归地将簇分割成两个最不相似的簇，直到每个观测值对应一个簇：

  层次聚类的示意图

在执行层次聚类之前，我们需要确定两个簇之间的相似度。在这里，我们列出了一些常用的相似度测量距离函数：

• 单连接：这指的是每个簇中两点之间的最短距离：

• 完全连接：这指的是每个簇中两点之间的最长距离：

• 平均链接：这指的是每个簇中两点之间的平均距离（其中是簇的大小，而是簇的大小）：

• 沃德方法：这指的是每个点到合并簇均值的平方距离之和（其中是的均值向量）：

在本菜谱中，我们对客户数据进行层次聚类。首先，我们从customer.csv加载数据，然后将其加载到客户数据框中。在数据中，我们找到五个客户账户信息的变量，它们是 ID、访问次数、平均消费、性别和年龄。由于每个变量的尺度不同，我们使用尺度函数来归一化尺度。

在对所有属性的尺度进行归一化后，我们使用hclust函数执行层次聚类。我们使用欧几里得距离作为距离度量，并使用沃德的最小方差方法进行聚合聚类。

最后，我们使用plot函数绘制层次簇的树状图。我们指定hang以在树状图的底部显示标签，并使用cex将标签缩小到正常大小的 70%。为了比较使用ward.D2和single方法生成的簇层次结构，我们在前面的图中使用single绘制另一个树状图（步骤 6）。

还有更多...

在执行层次聚类时，您可以选择不同的距离度量和方法。有关更多详细信息，您可以参考dist和hclust的文档：

> ? dist
> ? hclust

在本菜谱中，我们使用hclust执行聚合层次聚类；如果您想执行分裂层次聚类，可以使用diana函数：

1. 首先，您可以使用diana执行分裂层次聚类：

   > dv = diana(customer, metric = "euclidean")
   
   

2. 然后，您可以使用summary获取摘要信息：

   > summary(dv)
   
   

3. 最后，您可以使用plot函数绘制树状图和横幅：

   > plot(dv)
   
   

如果您想绘制水平树状图，可以使用dendextend包。使用以下步骤生成水平树状图：

1. 首先，安装并加载dendextend和magrittr包（如果您的 R 版本是 3.1 及以上，您不需要安装和加载magrittr包）：

   > install.packages("dendextend")
   > library(dendextend)
   > install.packages("margrittr")
   > library(magrittr)
   
   

2. 设置树状图：

   > dend = customer %>% dist %>% hclust %>% as.dendrogram
   
   

3. 最后，绘制水平树状图：

   dend %>% plot(horiz=TRUE, main = "Horizontal Dendrogram")
   
   

   

   水平树状图

将树切割成簇

在树状图中，我们可以看到簇的层次结构，但我们还没有将数据分组到不同的簇中。然而，我们可以确定树状图中有多少簇，并在某个树高切割树状图以将数据分离到不同的组中。在这个菜谱中，我们展示了如何使用cutree函数将数据分离成给定的簇数。

准备工作

为了执行cutree函数，你需要完成之前的步骤，通过生成 hclust 对象hc。

如何做...

执行以下步骤将簇的层次结构切割成给定的簇数：

1. 首先，将数据分类为四个组：

   > fit = cutree(hc, k = 4)
   
   

2. 然后，你可以检查数据的簇标签：

   > fit
    [1] 1 1 2 1 2 1 2 2 1 1 1 2 2 1 1 1 2 1 2 3 4 3 4 3 3 4 4 3 4
   [30] 4 4 3 3 3 4 4 3 4 4 4 4 4 4 4 3 3 4 4 4 3 4 3 3 4 4 4 3 4
   [59] 4 3
   
   

3. 计算每个簇内的数据数量：

   > table(fit)
   fit
    1  2  3  4 
   11  8 16 25 
   
   

4. 最后，你可以使用红色矩形边框可视化数据的聚类情况：

   > plot(hc)
   > rect.hclust(hc, k = 4 , border="red")
   
   

   

   使用红色矩形边框来区分树状图中不同的簇

它是如何工作的...

我们可以从前一个图中的树状图中确定簇的数量。在这个菜谱中，我们确定树中应该有四个簇。因此，我们在cutree函数中将簇数指定为4。除了使用簇数来切割树，你还可以指定height作为切割树的参数。

接下来，我们可以输出数据的簇标签，并使用table函数计算每个簇内的数据数量。从计数表中，我们发现大部分数据都在簇 4 中。最后，我们可以使用rect.hclust函数在簇周围绘制红色矩形，以显示数据如何被分类到四个簇中。

还有更多...

除了在所有层次簇周围绘制矩形，你还可以在某个簇周围放置红色矩形：

> rect.hclust(hc, k = 4 , which =2, border="red")

在某个簇周围绘制红色矩形。

此外，你可以使用dendextend包通过在簇周围添加红色矩形以不同颜色着色簇。你必须完成前一个菜谱中还有更多部分中概述的说明，并执行以下步骤：

1. 根据所属的簇给分支上色：

   > dend %>% color_branches(k=4) %>% plot(horiz=TRUE, main = "Horizontal Dendrogram")
   
   

2. 然后，你可以在簇周围添加红色矩形：

   > dend %>% rect.dendrogram(k=4,horiz=TRUE)
   
   

   

   在水平树状图中绘制簇周围的红色矩形

3. 最后，你可以添加一条线来显示树木切割的位置：

   > abline(v = heights_per_k.dendrogram(dend)["4"] + .1, lwd = 2, lty = 2, col = "blue")
   
   

   

   在水平树状图中绘制切割线

使用 k-means 方法聚类数据

k-means 聚类是一种平面聚类技术，它只产生一个包含 k 个聚类的分区。与不需要用户在开始时确定聚类数量的层次聚类不同，k-means 方法需要首先确定这一点。然而，由于构建层次树非常耗时，k-means 聚类比层次聚类快得多。在本食谱中，我们将演示如何在客户数据集上执行 k-means 聚类。

准备工作

在本食谱中，我们将继续使用客户数据集作为输入数据源来执行 k-means 聚类。

如何操作...

按以下步骤使用 k-means 方法对 customer 数据集进行聚类：

1. 首先，您可以使用 kmeans 对客户数据进行聚类：

   > set.seed(22)
   > fit = kmeans(customer, 4)
   > fit
   K-means clustering with 4 clusters of sizes 8, 11, 16, 25
   
   Cluster means:
    Visit.Time Average.Expense        Sex        Age
   1  1.3302016       1.0155226 -1.4566845  0.5591307
   2 -0.7771737      -0.5178412 -1.4566845 -0.4774599
   3  0.8571173       0.9887331  0.6750489  1.0505015
   4 -0.6322632      -0.7299063  0.6750489 -0.6411604
   
   Clustering vector:
    [1] 2 2 1 2 1 2 1 1 2 2 2 1 1 2 2 2 1 2 1 3 4 3 4 3 3 4 4 3
   [29] 4 4 4 3 3 3 4 4 3 4 4 4 4 4 4 4 3 3 4 4 4 3 4 3 3 4 4 4
   [57] 3 4 4 3
   
   Within cluster sum of squares by cluster:
   [1]  5.90040 11.97454 22.58236 20.89159
    (between_SS / total_SS =  74.0 %)
   
   Available components:
   
   [1] "cluster"      "centers"      "totss" 
   [4] "withinss"     "tot.withinss" "betweenss" 
   [7] "size"         "iter"         "ifault
   
   

2. 您可以使用 barplot 检查每个聚类的中心：

   > barplot(t(fit$centers), beside = TRUE,xlab="cluster", ylab="value")
   
   

   

   四个聚类中不同属性的中心的条形图

3. 最后，您可以在数据上绘制散点图，并根据聚类对点进行着色：

   > plot(customer, col = fit$cluster)
   
   

   

   显示数据根据其聚类标签着色的散点图

工作原理...

k-means 聚类是一种分区聚类方法。算法的目标是将 n 个对象划分为 k 个聚类，其中每个对象属于最近的均值所在的聚类。算法的目标是最小化 聚类内平方和（WCSS）。假设 x 是给定的观察值集合，S = 表示 k 个分区，而 是 的均值，那么我们可以将 WCSS 函数表述如下：

k-means 聚类的过程可以通过以下五个步骤来表示：

1. 指定 k 个聚类的数量。

2. 随机创建 k 个分区。

3. 计算分区中心。

4. 将对象与聚类中心最接近的关联起来。

5. 重复步骤 2、3 和 4，直到 WCSS 变化很小（或是最小化）。

在本食谱中，我们展示了如何使用 k-means 聚类对客户数据进行聚类。与层次聚类不同，k-means 聚类需要用户输入 K 的数量。在本例中，我们使用 K=4。然后，拟合模型的输出显示了每个聚类的规模、四个生成聚类的聚类均值、每个数据点的聚类向量、聚类内的平方和以及其他可用组件。

此外，您可以在条形图中绘制每个聚类的中心，这将提供更多关于每个属性如何影响聚类的细节。最后，我们在散点图中绘制数据点，并使用拟合的聚类标签根据聚类标签分配颜色。

相关内容

• 在 k-means 聚类中，您可以指定用于执行聚类分析的算法。您可以指定 Hartigan-Wong、Lloyd、Forgy 或 MacQueen 作为聚类算法。有关更多详细信息，请使用help函数参考kmeans函数的文档：

  >help(kmeans)
  
  

绘制双变量聚类图

在上一个配方中，我们使用了 k-means 方法将数据拟合到聚类中。然而，如果有超过两个变量，就无法在二维中显示数据的聚类方式。因此，您可以使用双变量聚类图首先将变量降低为两个成分，然后使用轴和圆等成分作为聚类来显示数据的聚类方式。在本配方中，我们将说明如何创建双变量聚类图。

准备工作

在本配方中，我们将继续使用customer数据集作为输入数据源来绘制双变量聚类图。

如何做...

执行以下步骤以绘制双变量聚类图：

1. 安装并加载聚类包：

   > install.packages("cluster")
   > library(cluster)
   
   

2. 然后，您可以绘制双变量聚类图：

   > clusplot(customer, fit$cluster, color=TRUE, shade=TRUE)
   
   

   

   客户数据集的双变量聚类图

3. 您还可以放大双变量聚类图：

   > par(mfrow= c(1,2))
   > clusplot(customer, fit$cluster, color=TRUE, shade=TRUE)
   > rect(-0.7,-1.7, 2.2,-1.2, border = "orange", lwd=2)
   > clusplot(customer, fit$cluster, color = TRUE, xlim = c(-0.7,2.2), ylim = c(-1.7,-1.2))
   
   

   

   双变量聚类图的放大

它是如何工作的...

在本配方中，我们绘制双变量聚类图以显示数据的聚类方式。要绘制双变量聚类图，我们首先需要安装cluster包并将其加载到 R 中。然后，我们使用clusplot函数从客户数据集中绘制双变量聚类图。在clustplot函数中，我们可以将shade设置为TRUE，将color设置为TRUE以显示带有颜色和阴影的聚类。根据前面的图（步骤 2），我们发现双变量使用两个成分作为 x 轴和 y 轴，这两个成分解释了 85.01%的点变异性。然后，数据点根据成分 1 和成分 2 在图上散布。同一聚类内的数据以相同颜色和阴影的圆圈表示。

除了在单个图中绘制四个聚类外，您还可以使用rect在给定 x 轴和 y 轴范围内的特定区域内添加一个矩形。然后，您可以使用clusplot函数中的xlim和ylim来放大图表，以检查每个聚类内的数据。

更多内容

clusplot函数使用princomp和cmdscale将原始特征维度降低到主成分。因此，可以看到数据如何在这两个成分作为 x 轴和 y 轴的单个图中聚类。要了解更多关于princomp和cmdscale的信息，可以使用help函数查看相关文档：

> help(cmdscale)
> help(princomp)

对于那些对如何使用cmdscale进行降维感兴趣的人，请执行以下步骤：

> mds = cmdscale(dist(customer), k = 2)
> plot(mds, col = fit$cluster)

关于缩放维度的数据散点图

比较聚类方法

使用不同的聚类方法将数据拟合到簇中后，您可能希望测量聚类的准确性。在大多数情况下，您可以使用簇内或簇间指标作为测量标准。我们现在介绍如何使用fpc包中的cluster.stat比较不同的聚类方法。

准备工作

为了执行聚类方法比较，需要完成前面的配方，生成customer数据集。

如何操作...

执行以下步骤以比较聚类方法：

1. 首先，安装并加载fpc包：

   > install.packages("fpc")
   > library(fpc)
   
   

2. 然后，您需要使用single方法进行层次聚类以聚类客户数据并生成对象hc_single：

   > single_c =  hclust(dist(customer), method="single")
   > hc_single = cutree(single_c, k = 4)
   
   

3. 使用complete方法进行层次聚类以聚类客户数据并生成对象hc_complete：

   > complete_c =  hclust(dist(customer), method="complete")
   > hc_complete =  cutree(complete_c, k = 4)
   
   

4. 然后，您可以使用 k-means 聚类将客户数据进行聚类并生成对象km：

   > set.seed(22)
   > km = kmeans(customer, 4)
   
   

5. 接下来，检索任何聚类方法的聚类验证统计信息：

   > cs = cluster.stats(dist(customer), km$cluster)
   
   

6. 通常，我们关注使用within.cluster.ss和avg.silwidth来验证聚类方法：

   > cs[c("within.cluster.ss","avg.silwidth")]
   $within.cluster.ss
   [1] 61.3489
   
   $avg.silwidth
   [1] 0.4640587
   
   

7. 最后，我们可以生成每种聚类方法的聚类统计信息并将它们列在表中：

   > sapply(list(kmeans = km$cluster, hc_single = hc_single, hc_complete = hc_complete), function(c) cluster.stats(dist(customer), c)[c("within.cluster.ss","avg.silwidth")])
    kmeans    hc_single hc_complete
   within.cluster.ss 61.3489   136.0092  65.94076
   avg.silwidth      0.4640587 0.2481926 0.4255961
   
   

如何工作...

在本配方中，我们展示了如何验证聚类。为了验证聚类方法，我们通常采用两种技术：簇间距离和簇内距离。在这些技术中，簇间距离越高，越好；簇内距离越低，越好。为了计算相关统计信息，我们可以将fpc包中的cluster.stat应用于拟合的聚类对象。

从输出中，within.cluster.ss测量表示簇内平方和，而avg.silwidth表示平均轮廓宽度。within.cluster.ss测量表示簇内对象的相关性；值越小，簇内相关对象越紧密。另一方面，轮廓是一个考虑簇内对象的相关性和簇之间分离程度的测量。数学上，我们可以定义每个点x的轮廓宽度如下：

在前面的公式中，a(x)是x与簇内所有其他点的平均距离，而b(x)是x与其他簇中点的平均距离的最小值。轮廓值通常在0到1之间；接近1的值表明数据聚类得更好。

最后一步生成的摘要表显示，完全层次聚类方法在within.cluster.ss和avg.silwidth方面优于单一层次聚类方法和 k-means 聚类。

相关内容

• kmeans函数也输出统计信息（例如，withinss和betweenss），供用户验证聚类方法：

  > set.seed(22)
  > km = kmeans(customer, 4)
  > km$withinss
  [1]  5.90040 11.97454 22.58236 20.89159
  > km$betweenss
  [1] 174.6511
  
  

从聚类中提取轮廓信息

轮廓信息是用于验证数据聚类的测量指标。在之前的步骤中，我们提到聚类的测量涉及计算数据在每个聚类内聚类的紧密程度，以及衡量不同聚类之间的距离。轮廓系数结合了簇内和簇间距离的测量。输出值通常在 0 到 1 之间；越接近 1，聚类越好。在本步骤中，我们将介绍如何计算轮廓信息。

准备工作

为了从聚类中提取轮廓信息，你需要完成之前的步骤，通过生成 customer 数据集。

如何操作...

执行以下步骤来计算轮廓信息：

1. 使用 kmeans 生成一个 k-means 对象，km：

   > set.seed(22)
   > km = kmeans(customer, 4)
   
   

2. 然后，你可以计算轮廓信息：

   > kms = silhouette(km$cluster,dist(customer))
   > summary(kms)
   Silhouette of 60 units in 4 clusters from silhouette.default(x = km$cluster, dist = dist(customer)) :
    Cluster sizes and average silhouette widths:
    8        11        16        25 
   0.5464597 0.4080823 0.3794910 0.5164434 
   Individual silhouette widths:
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.1931  0.4030  0.4890  0.4641  0.5422  0.6333 
   
   

3. 接下来，你可以绘制轮廓信息图：

   > plot(kms)
   
   

   

   k-means 聚类结果的轮廓图

它是如何工作的...

在本步骤中，我们展示了如何使用轮廓图来验证聚类。你可以首先检索轮廓信息，它显示了聚类大小、平均轮廓宽度和个体轮廓宽度。轮廓系数是一个介于 0 到 1 之间的值；越接近 1，聚类质量越好。

最后，我们使用 plot 函数绘制轮廓图。图的左侧显示了水平线的数量，代表聚类的数量。图的右侧列显示了其自身聚类图与下一个相似聚类图平均相似度的差值。平均轮廓宽度在图的底部显示。

相关内容

• 对于那些对如何计算轮廓感兴趣的人，请参考维基百科上的轮廓值条目：[zh.wikipedia.org/wiki/轮廓 _%28 聚类%29](zh.wikipedia.org/wiki/轮廓 _%28 聚类%29)

获取 k-means 的最佳聚类数量

虽然 k-means 聚类速度快且易于使用，但它需要在开始时输入 k 值。因此，我们可以使用平方和来确定哪个 k 值最适合找到 k-means 的最佳聚类数量。在下面的步骤中，我们将讨论如何找到 k-means 聚类方法的最佳聚类数量。

准备工作

为了找到最佳聚类数量，你需要完成之前的步骤，通过生成 customer 数据集。

如何操作...

执行以下步骤以找到 k-means 聚类最佳聚类数量：

1. 首先，计算不同数量聚类的内部平方和 (withinss)：

   > nk = 2:10
   > set.seed(22)
   > WSS = sapply(nk, function(k) {
   +     kmeans(customer, centers=k)$tot.withinss
   + })
   > WSS
   [1] 123.49224  88.07028  61.34890  48.76431  47.20813
   [6]  45.48114  29.58014  28.87519  23.21331
   
   

2. 然后，你可以使用线图绘制不同 k 值的内部平方和：

   > plot(nk, WSS, type="l", xlab= "number of k", ylab="within sum of squares")
   
   

   

   关于不同 k 值的内部平方和的线图

3. 接下来，你可以计算不同数量聚类的平均轮廓宽度（avg.silwidth）：

   > SW = sapply(nk, function(k) {
   +   cluster.stats(dist(customer), kmeans(customer, centers=k)$cluster)$avg.silwidth
   + })
   > SW
   [1] 0.4203896 0.4278904 0.4640587 0.4308448 0.3481157
   [6] 0.3320245 0.4396910 0.3417403 0.4070539
   
   

4. 然后，你可以使用线图来绘制不同k值的平均轮廓宽度：

   > plot(nk, SW, type="l", xlab="number of clusers", ylab="average silhouette width")
   
   

   

   与不同数量的 k 相关的平均轮廓宽度的线图

5. 获取最大聚类数量：

   > nk[which.max(SW)]
   [1] 4
   
   

它是如何工作的...

在这个配方中，我们通过迭代地获取平方和与平均轮廓值之和来演示如何找到最佳聚类数量。对于平方和，较低的值表示质量更好的聚类。通过绘制不同数量k的平方和，我们发现图表的肘部在k=4。

另一方面，我们使用cluster.stats根据不同数量的聚类计算平均轮廓宽度。我们还可以使用线图来绘制不同数量聚类对应的平均轮廓宽度。前面的图（步骤 4）显示最大平均轮廓宽度出现在k=4。最后，我们使用which.max来获取 k 的值，以确定最大平均轮廓宽度的位置。

参考信息

• 对于那些对平方和如何计算感兴趣的人，请参阅维基百科上的K-means 聚类条目：en.wikipedia.org/wiki/K-means_clustering

使用基于密度的方法进行聚类

作为距离测量的替代，你可以使用基于密度的测量来聚类数据。这种方法找到一个比剩余区域密度更高的区域。最著名的方法之一是 DBSCAN。在下面的配方中，我们将演示如何使用 DBSCAN 进行基于密度的聚类。

准备工作

在这个配方中，我们将使用由mlbench包生成的模拟数据。

如何操作...

执行以下步骤以执行基于密度的聚类：

1. 首先，安装并加载fpc和mlbench包：

   > install.packages("mlbench")
   > library(mlbench)
   > install.packages("fpc")
   > library(fpc)
   
   

2. 然后，你可以使用mlbench库绘制 Cassini 问题图：

   > set.seed(2)
   > p = mlbench.cassini(500)
   > plot(p$x)
   
   

   

   Cassini 问题图

3. 接下来，你可以根据密度的测量来对数据进行聚类：

   > ds = dbscan(dist(p$x),0.2, 2, countmode=NULL, method="dist")
   > ds
   dbscan Pts=500 MinPts=2 eps=0.2
    1   2   3
   seed  200 200 100
   total 200 200 100
   
   

4. 以不同的聚类标签作为颜色绘制散点图中的数据：

   > plot(ds, p$x)
   
   

   

   根据聚类标签着色的数据散点图

5. 你还可以使用dbscan来预测数据点属于哪个聚类。在这个例子中，首先在矩阵p中输入三个输入：

   > y = matrix(0,nrow=3,ncol=2)
   > y[1,] = c(0,0)
   > y[2,] = c(0,-1.5)
   > y[3,] = c(1,1)
   > y
    [,1] [,2]
   [1,]    0  0.0
   [2,]    0 -1.5
   [3,]    1  1.0
   
   

6. 然后，你可以预测数据属于哪个聚类：

   > predict(ds, p$x, y)
   [1] 3 1 2
   
   

它是如何工作的...

基于密度的聚类利用密度可达性和密度连通性的概念，这使得它在发现非线性形状的聚类中非常有用。在讨论基于密度的聚类过程之前，必须解释一些重要的背景概念。基于密度的聚类考虑两个参数：eps和MinPts。eps代表邻域的最大半径；MinPts表示eps邻域内点的最小数量。有了这两个参数，我们可以定义核心点为在eps邻域内有超过MinPts个点的点。此外，我们还可以定义边界点为点数少于MinPts，但位于核心点邻域内的点。然后，我们可以定义核心对象为如果点p的eps邻域内的点数超过MinPts。

此外，我们必须定义两点之间的可达性。我们可以这样说，如果点q位于点p的eps邻域内，并且p是一个核心对象，那么点p是直接从点q密度可达的。然后，我们可以定义，如果存在一个点的链，p[1]，p[2]，...，p[n]，其中 p[1] = q，p[n] = p，并且对于 1 <= i <= n，p[i]+1 相对于 Eps 和MinPts直接从 pi 密度可达，那么点p是通用且从点q密度可达的：

点 p 和 q 是密度可达的

在对基于密度的聚类有一个初步概念之后，我们可以通过以下步骤来阐述最流行的基于密度的聚类算法 DBSCAN 的过程，如图所示：

1. 随机选择一个点，p。

2. 根据Eps和MinPts检索所有从p密度可达的点。

3. 如果p是一个核心点，那么就形成了一个簇。否则，如果它是边界点，并且没有点从p密度可达，则过程将标记该点为噪声，并继续访问下一个点。

4. 重复此过程，直到所有点都被访问。

在这个菜谱中，我们展示了如何使用基于密度的 DBSCAN 方法来聚类客户数据。首先，我们必须安装和加载mlbench和fpc库。mlbench包提供了许多方法来生成不同形状和大小的模拟数据。在这个例子中，我们生成一个 Cassini 问题图。

接下来，我们在 Cassini 数据集上执行dbscan以聚类数据。我们指定可达距离为 0.2，点的最小可达数为 2，进度报告为 null，并使用距离作为测量标准。聚类方法成功地将数据聚为三个大小分别为 200、200 和 100 的簇。通过在图上绘制点和簇标签，我们看到 Cassini 图的三部分被不同颜色分开。

fpc包还提供了一个predict函数，您可以使用此函数预测输入矩阵的聚类标签。点 c(0,0)被分类到聚类 3，点 c(0, -1.5)被分类到聚类 1，点 c(1,1)被分类到聚类 2。

参见

• fpc包包含灵活的聚类过程，并具有有用的聚类分析函数。例如，您可以使用plotcluster函数生成判别投影图。有关更多信息，请参阅以下文档：

  > help(plotcluster)
  
  

使用基于模型的方法进行聚类数据

与使用启发式方法且不依赖于正式模型的层次聚类和 k-means 聚类相比，基于模型聚类技术假设各种数据模型，并应用 EM 算法以获得最可能的模型，并进一步使用该模型推断最可能的聚类数量。在本配方中，我们将演示如何使用基于模型的方法来确定最可能的聚类数量。

准备工作

为了执行基于模型的方法对客户数据进行聚类，您需要完成之前的配方，生成客户数据集。

如何操作...

执行以下步骤以执行基于模型聚类：

1. 首先，请安装并加载库mclust：

   > install.packages("mclust")
   > library(mclust)
   
   

2. 然后，您可以在customer数据集上执行基于模型聚类：

   > mb = Mclust(customer)
   > plot(mb)
   
   

3. 然后，您可以按 1 键获取 BIC 与组件数量的关系：

   BIC 与组件数量的关系图

4. 然后，您可以按 2 键显示关于不同特征组合的分类：

   根据不同特征组合的分类图

5. 按 3 键显示关于不同特征组合的分类不确定性：

   根据不同特征组合的分类不确定性图

6. 接下来，按 4 键绘制密度估计：

   密度估计图

7. 然后，您可以按 0 键绘制密度并退出绘图菜单。

8. 最后，使用summary函数获取最可能的模型和聚类数量：

   > summary(mb)
   ----------------------------------------------------
   Gaussian finite mixture model fitted by EM algorithm 
   ----------------------------------------------------
   
   Mclust VII (spherical, varying volume) model with 5 components:
   
    log.likelihood  n df       BIC       ICL
    -218.6891 60 29 -556.1142 -557.2812
   
   Clustering table:
    1  2  3  4  5
   11  8 17 14 10
   
   

如何工作...

与基于启发式方法构建聚类不同，基于模型聚类使用基于概率的方法。基于模型聚类假设数据是由潜在的概率分布生成的，并试图从数据中恢复分布。一种常见的基于模型的方法是使用有限混合模型，它为概率分布的分析提供了一个灵活的建模框架。有限混合模型是组件概率分布的线性加权总和。假设数据*y=(y[1],y[2]…y[n])*包含 n 个独立的多变量观测值；G 是组件数量；有限混合模型的似然可以表示为：

和 是混合模型中第 k 个成分的密度和参数，而 (![如何工作...] 和 ) 是观察属于第 k 个成分的概率。

基于模型的聚类过程有几个步骤：首先，过程选择组件概率分布的数量和类型。然后，它拟合一个有限混合模型并计算组件成员的后验概率。最后，它将每个观察的成员分配给概率最大的组件。

在这个菜谱中，我们展示了如何使用基于模型的聚类方法对数据进行聚类。我们首先在 R 中安装并加载Mclust库。然后，我们使用Mclust函数将客户数据拟合到基于模型的方法中。

数据拟合到模型后，我们根据聚类结果绘制模型。有四种不同的图：BIC 图、分类图、不确定性和密度图。BIC 图显示了 BIC 值，人们可以使用这个值来选择簇的数量。分类图显示了数据如何根据不同的维度组合进行聚类。不确定图显示了根据不同维度组合的分类不确定性。密度图显示了轮廓中的密度估计。

你也可以使用summary函数来获取最可能的模型和最可能的簇数量。对于这个例子，最可能的簇数量是五个，BIC 值为-556.1142。

参见

• 对于那些对Mclust如何工作感兴趣的人，请参阅以下来源：C. Fraley, A. E. Raftery, T. B. Murphy 和 L. Scrucca (2012). mclust Version 4 for R: Normal Mixture Modeling for Model-Based Clustering, Classification, and Density Estimation. 技术报告第 597 号，华盛顿大学统计学系。

可视化相似性矩阵

相似性矩阵可以用作衡量聚类质量的一个指标。为了可视化矩阵，我们可以在距离矩阵上使用热图。在图中，低相似性（或高相似性）的条目以较深的颜色绘制，这有助于识别数据中的隐藏结构。在这个菜谱中，我们将讨论一些可视化相似性矩阵的有用技术。

准备工作

为了可视化相似性矩阵，你需要完成之前的菜谱，生成客户数据集。此外，还需要生成并存储在变量km中的 k-means 对象。

如何操作...

执行以下步骤以可视化相似性矩阵：

1. 首先，安装并加载seriation包：

   > install.packages("seriation")
   > library(seriation)
   
   

2. 然后，你可以使用dissplot在热图上可视化相似性矩阵：

   > dissplot(dist(customer), labels=km$cluster, options=list(main="Kmeans Clustering With k=4"))
   
   

   

   k-means 聚类的相似性图

3. 接下来，在热图中应用dissplot于层次聚类：

   > complete_c =  hclust(dist(customer), method="complete")
   > hc_complete =  cutree(complete_c, k = 4)
   > dissplot(dist(customer), labels=hc_complete, options=list(main="Hierarchical Clustering"))
   
   

   

   层次聚类的相似度图

它是如何工作的...

在这个菜谱中，我们使用不相似度图来可视化不相似度矩阵。我们首先安装并加载seriation包，然后对 k-means 聚类输出应用dissplot函数，生成前面的图（步骤 2）。

这显示，相似度高的聚类用较深的颜色表示，而不同组合的聚类用较浅的颜色表示。因此，我们可以看到对应聚类（如聚类 4 到聚类 4）是斜对角且颜色较深的。另一方面，彼此不相似的聚类用较浅的颜色表示，并且远离对角线。

同样，我们可以在层次聚类的输出上应用dissplot函数。图中的生成图（步骤 3）显示了每个聚类的相似性。

还有更多...

除了使用dissplot来可视化不相似度矩阵外，还可以使用dist和image函数来可视化距离矩阵。在生成的图中，密切相关项用红色表示。不太相关项则更接近白色：

> image(as.matrix(dist(customer)))

客户数据集的距离矩阵图

为了绘制树状图和热图以显示数据的聚类方式，您可以使用heatmap函数：

> cd=dist(customer)
> hc=hclust(cd)
> cdt=dist(t(customer))
> hcc=hclust(cdt)
> heatmap(customer, Rowv=as.dendrogram(hc), Colv=as.dendrogram(hcc))

列和行侧带有树状图的热图

外部验证聚类

除了生成统计信息来验证生成的聚类的质量外，您还可以使用已知的数据聚类作为基准来比较不同的聚类方法。在这个菜谱中，我们将演示聚类方法在已知聚类数据上的差异。

准备工作

在这个菜谱中，我们将继续使用手写数字作为聚类输入；您可以在作者的 GitHub 页面上找到该图：github.com/ywchiu/ml_R_cookbook/tree/master/CH9。

如何做...

执行以下步骤以使用不同的聚类技术对数字进行聚类：

1. 首先，您需要安装并加载png包：

   > install.packages("png")
   > library(png)
   
   

2. 然后，请从handwriting.png读取图像，并将读取的数据转换为散点图：

   > img2 = readPNG("handwriting.png", TRUE)
   > img3 = img2[,nrow(img2):1]
   > b = cbind(as.integer(which(img3 < -1) %% 28), which(img3 < -1) / 28)
   > plot(b, xlim=c(1,28), ylim=c(1,28))
   
   

   

   手写数字的散点图

3. 对手写数字执行 k-means 聚类方法：

   > set.seed(18)
   > fit = kmeans(b, 2)
   > plot(b, col=fit$cluster)
   > plot(b, col=fit$cluster,  xlim=c(1,28), ylim=c(1,28))
   
   

   

   手写数字的 k-means 聚类结果

4. 接下来，对手写数字执行dbscan聚类方法：

   > ds = dbscan(b, 2)
   > ds
   dbscan Pts=212 MinPts=5 eps=2
    1   2
   seed  75 137
   total 75 137
   > plot(ds, b,  xlim=c(1,28), ylim=c(1,28))
   
   

手写数字的 DBSCAN 聚类结果

它是如何工作的...

在这个菜谱中，我们展示了不同的聚类方法在处理手写数据集时的效果。聚类的目的是将 1 和 7 分离到不同的簇中。我们执行不同的技术来观察数据在 k-means 和 DBSCAN 方法下的聚类情况。

为了生成数据，我们使用 Windows 应用程序 paint.exe 创建一个 28 x 28 像素的 PNG 文件。然后我们使用 readPNG 函数读取 PNG 数据，并将读取的 PNG 数据点转换为散点图，该图显示了 17 中的手写数字。

在数据读取之后，我们对手写数字进行聚类技术处理。首先，我们在数据集上执行 k-means 聚类，其中 k=2。由于 k-means 聚类使用距离度量，构建的簇覆盖了 1 和 7 数字区域。然后我们对数据集执行 DBSCAN。由于 DBSCAN 是一种基于密度的聚类技术，它成功地将数字 1 和数字 7 分离到不同的簇中。

参见

• 如果你对如何在 R 中读取各种图形格式感兴趣，你可以参考以下文档：

  > help(package="png")