R-机器学习第二版-三-

R 机器学习第二版（三）

原文：annas-archive.org/md5/e79782597ad902714a518953c5d1e70b

译者：飞龙

协议：CC BY-NC-SA 4.0

第八章：发现模式—使用关联规则进行市场篮子分析

回想一下你上一次冲动购买的经历。也许你在超市结账时买了一包口香糖或一块巧克力棒。也许你在深夜去买尿布和奶粉时顺便拿了一瓶含咖啡因的饮料或一六瓶啤酒。你可能甚至是在书店推荐下冲动地买了这本书。这些冲动购买可不是偶然的，零售商通过使用复杂的数据分析技术来识别那些能够驱动零售行为的模式。

在过去的几年里，这种推荐系统主要依赖于营销专业人员和库存经理或采购员的主观直觉。近年来，随着条形码扫描仪、计算机化库存系统和在线购物趋势积累了大量的交易数据，机器学习被越来越多地应用于学习购买模式。由于这种技术经常应用于超市数据，因此这一实践通常被称为市场篮子分析。

尽管该技术起源于购物数据，但它在其他背景下也同样有用。当你读完这一章时，你将能够将市场篮子分析技术应用到你自己的任务中，无论这些任务是什么。一般而言，这项工作包括：

• 使用简单的性能指标在大型数据库中寻找关联

• 理解交易数据的特殊性

• 知道如何识别有用且可操作的模式

市场篮子分析的结果是可操作的模式。因此，当我们应用这项技术时，即使你与零售链没有任何关联，也很可能会识别出对你的工作有用的应用。

理解关联规则

市场篮子分析的基本构件是可能出现在任何给定交易中的商品。一个或多个商品组成的组被括号括起来，表示它们形成了一个集合，或者更具体地说，是一个在数据中有规律出现的项集。交易是通过项集来指定的，例如在典型的杂货店中可能找到的以下交易：

市场篮子分析的结果是一组关联规则，这些规则指定了在商品项集之间的关系中发现的模式。关联规则总是由项集的子集组成，并通过将规则的左侧（LHS）中的一个项集与右侧（RHS）中的另一个项集联系起来来表示。LHS 是需要满足的条件，以触发规则，而 RHS 是满足该条件后预期的结果。一个从示例交易中识别出的规则可以表示为：

用简单的话来说，这个关联规则表示，如果花生酱和果冻一起购买，那么面包也很可能会被购买。换句话说，“花生酱和果冻意味着面包。”

在零售交易数据库的背景下，关联规则不是用来进行预测，而是用于在大型数据库中进行无监督的知识发现。这与前几章介绍的分类和数值预测算法不同。尽管如此，你会发现关联规则学习器与第五章中的分类规则学习器紧密相关，并共享许多相同的特征，分而治之——使用决策树和规则进行分类。

由于关联规则学习器是无监督的，因此算法无需进行训练；数据无需提前标注。程序只是被释放到数据集上，希望能够发现有趣的关联。当然，缺点是，除了通过定性有用性评估规则学习器外，没有简单的方法来客观衡量其性能——通常需要某种形式的人工检查。

虽然关联规则最常用于市场篮分析，但它们在发现多种不同类型数据的模式中也很有帮助。其他潜在应用包括：

• 在癌症数据中搜索有趣且频繁出现的 DNA 和蛋白质序列模式

• 寻找与信用卡或保险欺诈使用相关联的购买模式或医疗索赔模式

• 识别导致客户取消手机服务或升级有线电视套餐的行为组合

关联规则分析用于在大量元素中搜索有趣的关联。人类能够凭直觉进行这样的洞察，但通常需要专家级的知识或大量经验，才能完成规则学习算法在几分钟甚至几秒钟内就能完成的任务。此外，一些数据集对于人类来说过于庞大和复杂，难以从中找到“针尖”。

关联规则学习的 Apriori 算法

正如人类感到具有挑战性一样，交易数据使得关联规则挖掘对于机器来说也是一项挑战。交易数据集通常非常庞大，无论是交易数量还是监控的商品或特征数量都很庞大。问题在于，潜在项集的数量随着特征数量的增加呈指数级增长。假设有k个商品可能会出现在集合中或不出现在集合中，那么有2^k个潜在的项集可能成为规则。一个只卖 100 种商品的零售商，可能需要评估大约2¹⁰⁰ = 1.27e+30个项集——这看似是不可能完成的任务。

与其逐一评估这些项集，不如利用更智能的规则学习算法，利用许多潜在的商品组合在实践中很少出现这一事实。例如，即使商店同时出售汽车配件和女性化妆品，*{机油, 口红}*这种组合也可能极其不常见。通过忽略这些稀有（而且可能不太重要）组合，可以将规则搜索的范围限制到一个更易管理的大小。

已经做了大量工作来识别启发式算法，以减少需要搜索的项集数量。也许最广泛使用的高效搜索大数据库中规则的方法是Apriori算法。该算法由 Rakesh Agrawal 和 Ramakrishnan Srikant 于 1994 年提出，此后已成为关联规则学习的代名词。其名称来源于该算法利用一种简单的先验信念（即a priori）关于频繁项集的属性。

在我们深入讨论之前，值得注意的是，这种算法与所有学习算法一样，并非没有优缺点。以下是其中的一些优缺点：

优势	劣势

|

• 能够处理大量的交易数据

• 产生易于理解的规则

• 有助于“数据挖掘”和发现数据库中的意外知识

|

• 对于小型数据集不太有用

• 需要努力将真正的洞察力与常识区分开

• 容易从随机模式中得出虚假的结论

|

如前所述，Apriori 算法通过简单的a priori信念来减少关联规则搜索空间：频繁项集的所有子集也必须是频繁的。这一启发式方法被称为Apriori 属性。利用这一精明的观察，我们可以大幅限制需要搜索的规则数量。例如，*{机油, 口红}集只有在{机油}和{口红}*都频繁出现时，才能是频繁的。因此，如果机油或口红出现频率较低，那么包含这些物品的任何集合都可以从搜索中排除。

注意

如需了解 Apriori 算法的更多细节，请参阅：Agrawal R, Srikant R. 快速挖掘关联规则的算法。第 20 届国际大型数据库会议论文集，1994:487-499。

为了了解这一原则如何应用于更现实的场景，让我们考虑一个简单的交易数据库。下表显示了一个假想医院礼品店的五笔已完成交易：

交易编号	购买的物品
1	{花, 康复卡片, 苏打水}
2	{毛绒玩具熊, 花, 气球, 巧克力棒}
3	{康复卡片, 巧克力棒, 花}
4	{毛绒玩具熊, 气球, 苏打水}
5	{花, 康复卡片, 苏打水}

通过查看购买集，可以推断出一些典型的购买模式。拜访生病的朋友或家人时，人们往往会购买一张祝早日康复的卡片和鲜花，而拜访新妈妈的人则倾向于购买毛绒玩具熊和气球。这些模式之所以引人注目，是因为它们频繁出现，足以引起我们的兴趣；我们只需运用一些逻辑和专业经验来解释这些规则。

类似地，Apriori 算法通过使用项集“有趣性”的统计度量来定位更大事务数据库中的关联规则。在接下来的章节中，我们将探讨 Apriori 如何计算这些有趣性的度量，以及如何将它们与 Apriori 属性结合起来，以减少需要学习的规则数量。

衡量规则兴趣 – 支持度与置信度

是否将某个关联规则视为有趣，取决于两个统计度量：支持度和置信度。通过为每个度量提供最小阈值，并应用 Apriori 原则，可以轻松地大幅度限制报告的规则数量，甚至可能只识别出显而易见或常识性的规则。因此，重要的是要仔细理解在这些标准下被排除的规则类型。

项集或规则的支持度衡量它在数据中出现的频率。例如，项集 {祝早日康复的卡片，鲜花} 在医院礼品店数据中的支持度为 3 / 5 = 0.6。类似地，{祝早日康复的卡片} → {鲜花} 的支持度也是 0.6。支持度可以计算任何项集，甚至是单一项；例如，{糖果棒} 的支持度为 2 / 5 = 0.4，因为糖果棒出现在 40% 的购买中。定义项集 X 的支持度的函数可以如下定义：

这里，N 是数据库中事务的数量，而 count(X) 是包含项集 X 的事务数量。

一条规则的置信度是其预测能力或准确度的度量。它被定义为包含 X 和 Y 的项集的支持度除以仅包含 X 的项集的支持度：

本质上，置信度告诉我们，在多少比例的交易中，物品或项集 X 的存在导致了物品或项集 Y 的存在。记住，X 导致 Y 的置信度和 Y 导致 X 的置信度是不同的。例如，{鲜花} → {祝贺卡} 的置信度是 0.6 / 0.8 = 0.75。相比之下，{祝贺卡} → {鲜花} 的置信度是 0.6 / 0.6 = 1.0。这意味着购买鲜花的交易有 75% 的可能性会同时购买祝贺卡，而购买祝贺卡的交易则有 100% 的可能性会购买鲜花。这些信息对于礼品店管理层非常有用。

提示

你可能已经注意到支持度、置信度和贝叶斯概率规则之间的相似性，这些内容在第四章，概率学习——使用朴素贝叶斯分类中有提到。事实上，support(A, B) 就是 P(A∩B)，而 confidence(A → B) 就是 P(B|A)。只是上下文不同。

类似 {祝贺卡} → {鲜花} 的规则被称为强规则，因为它们具有较高的支持度和置信度。找到更多强规则的一种方法是检查礼品店中每一对物品的所有可能组合，衡量支持度和置信度的值，然后仅报告那些符合特定兴趣水平的规则。然而，如前所述，这种策略通常只适用于最小的数据集。

在下一部分中，你将看到 Apriori 算法如何利用最小支持度和置信度水平，并结合 Apriori 原则，通过减少规则数量到一个更易管理的水平，快速找到强规则。

使用 Apriori 原则构建规则集

回忆一下，Apriori 原则指出，频繁项集的所有子集也必须是频繁的。换句话说，如果 {A, B} 是频繁的，那么 {A} 和 {B} 都必须是频繁的。还要记住，根据定义，支持度表示项集在数据中出现的频率。因此，如果我们知道 {A} 未达到所需的支持度阈值，那么就没有理由考虑 {A, B} 或任何包含 {A} 的项集；它不可能是频繁的。

Apriori 算法利用这种逻辑，在实际评估潜在的关联规则之前，先排除不可能的规则。创建规则的实际过程分为两个阶段：

1. 确定所有满足最小支持度阈值的项集。

2. 使用满足最小置信度阈值的项集来创建规则。

第一个阶段在多个迭代中进行。每次迭代都涉及评估一组越来越大的项集的支持度。例如，第一次迭代评估 1 项项集（1-itemsets），第二次迭代评估 2 项项集，以此类推。每次迭代i的结果是满足最低支持度阈值的所有i-项集的集合。

第i次迭代中的所有项集将被合并，以生成第i + 1次迭代的候选项集进行评估。但是，Apriori 原理甚至可以在下一轮开始之前就排除其中的一些项集。如果在第一次迭代中*{A}、{B}和{C}是频繁项集，而{D}不是，那么第二次迭代将只考虑{A, B}、{A, C}和{B, C}。因此，算法只需要评估三个项集，而不是评估包含D的六个项集，如果没有通过先验*排除这些项集的话。

继续这个思路，假设在第二次迭代中，发现*{A, B}和{B, C}是频繁项集，但{A, C}不是。尽管第三次迭代通常会通过评估{A, B, C}的支持度来开始，但并不强制要求这一步一定发生。为什么不呢？Apriori 原理指出，由于子集{A, C}不是频繁的，{A, B, C}*不可能是频繁项集。因此，在第三次迭代中没有生成新的项集，算法可能会停止。

此时，Apriori 算法的第二阶段可以开始。给定频繁项集的集合，从所有可能的子集生成关联规则。例如，*{A, B}将生成{A} → {B}和{B} → {A}*的候选规则。这些规则将根据最低置信度阈值进行评估，任何不满足所需置信度水平的规则都会被淘汰。

示例 – 使用关联规则识别常购商品

如本章介绍所述，市场篮子分析在许多实体店和在线零售商的推荐系统背后默默发挥着作用。通过学习得到的关联规则能够指示出经常一起购买的商品组合。了解这些模式为杂货连锁店提供了新的洞察力，帮助优化库存、进行促销广告或组织商店的物理布局。例如，如果顾客经常购买咖啡或橙汁与早餐点心一起食用，那么通过将点心移近咖啡和橙汁的位置，可能会提高利润。

在本教程中，我们将对来自一家杂货店的交易数据进行市场篮子分析。然而，这些技术可以应用于许多不同类型的问题，从电影推荐到约会网站，再到寻找药物之间的危险相互作用。在这个过程中，我们将看到 Apriori 算法如何高效地评估一个潜在的庞大的关联规则集。

第一步 – 收集数据

我们的市场篮分析将利用来自一家现实世界杂货店运营一个月的购买数据。数据包含 9,835 笔交易，约合每天 327 笔交易（在一个 12 小时的营业日中大约 30 笔交易），表明该零售商规模既不特别大，也不特别小。

注意

这里使用的数据集改编自arules R 包中的Groceries数据集。有关更多信息，请参见：Hahsler M, Hornik K, Reutterer T. 概率数据建模对挖掘关联规则的影响。载：Gaul W, Vichi M, Weihs C, ed. 分类学、数据分析与知识组织研究：从数据与信息分析到知识工程。纽约：Springer; 2006:598–605。

一般的杂货店提供种类繁多的商品。可能会有五种不同品牌的牛奶、十几种不同类型的洗衣粉和三种品牌的咖啡。鉴于零售商的规模适中，我们假设他们并不太关心寻找只适用于某一特定品牌牛奶或洗衣粉的规则。考虑到这一点，可以从购买数据中去除所有品牌名称。这样，商品种类减少为 169 种，更易于管理，涵盖了鸡肉、冷冻餐、玛琪琳和苏打水等广泛类别。

提示

如果你希望识别非常具体的关联规则——例如，客户是否喜欢将葡萄酒或草莓果酱与花生酱搭配——你将需要大量的事务数据。大型连锁零售商使用数百万笔交易的数据库来发现特定品牌、颜色或口味的物品之间的关联。

你有什么猜测关于哪些类型的物品可能会一起购买吗？葡萄酒和奶酪会是常见的搭配吗？面包和黄油？茶和蜂蜜？让我们深入分析这些数据，看看是否能确认我们的猜测。

第 2 步——探索与准备数据

事务数据以与我们之前使用的格式略有不同的方式存储。我们以前的大多数分析使用的是矩阵形式的数据，其中行表示示例实例，列表示特征。由于矩阵格式的结构，所有示例都必须具有完全相同的特征集。

相比之下，事务数据的形式更加自由。像往常一样，数据中的每一行表示一个单独的示例——在这种情况下，是一次交易。然而，与其有固定数量的特征，每条记录由一个用逗号分隔的物品列表组成，数量从一个到多个不等。实质上，特征可能会因示例而异。

提示

为了跟随本次分析，请从 Packt Publishing 网站下载groceries.csv文件，并将其保存在你的 R 工作目录中。

grocery.csv 文件的前五行如下：

citrus fruit,semi-finished bread,margarine,ready soups
tropical fruit,yogurt,coffee
whole milk
pip fruit,yogurt,cream cheese,meat spreads
other vegetables,whole milk,condensed milk,long life bakery product

这些行表示五个独立的杂货店交易。第一次交易包含了四个项目：柑橘类水果、半成品面包、玛格丽琳和即食汤。相比之下，第三次交易仅包含一个项目：全脂牛奶。

假设我们像之前的分析那样，尝试使用 read.csv() 函数加载数据。R 会欣然执行，并将数据读取成如下的矩阵形式：

你会注意到，R 创建了四列来存储事务数据中的项目：V1、V2、V3 和 V4。虽然这看起来合情合理，但如果我们使用这种格式的数据，后续会遇到问题。R 选择创建四个变量是因为第一行恰好有四个用逗号分隔的值。然而，我们知道，杂货购买可能包含超过四个项目；在四列设计中，这样的交易会被拆分到矩阵中的多行。我们可以尝试通过将包含最多项目的交易放在文件的顶部来解决这个问题，但这忽略了另一个更为棘手的问题。

通过这种方式构建数据，R 创建了一组特征，不仅记录了交易中的商品，还记录了它们出现的顺序。如果我们把学习算法想象成是在尝试找出 V1、V2、V3 和 V4 之间的关系，那么 V1 中的全脂牛奶可能会与 V2 中的全脂牛奶有所不同。相反，我们需要的是一个数据集，它不会把一笔交易视为一组需要填写（或不填写）特定商品的位置，而是将其视为一个市场购物篮，里面包含或不包含每一个特定商品。

数据准备 – 为事务数据创建稀疏矩阵

解决这个问题的方法是使用一种名为 稀疏矩阵 的数据结构。你可能还记得，我们在第四章中，概率学习 – 使用朴素贝叶斯分类，使用了稀疏矩阵来处理文本数据。和前面的数据集一样，稀疏矩阵中的每一行表示一次交易。然而，稀疏矩阵为每个可能出现在某人购物袋中的商品提供了一个列（即特征）。由于我们的杂货店数据中有 169 种不同的商品，所以我们的稀疏矩阵将包含 169 列。

为什么不直接像我们在大多数分析中那样将其存储为数据框（data frame）呢？原因是随着更多的交易和商品被添加，传统的数据结构很快会变得太大，无法适应可用的内存。即使是这里使用的相对较小的交易数据集，矩阵也包含近 170 万个单元格，其中大多数是零（因此称为“稀疏”矩阵——非零值非常少）。由于存储所有这些零值没有任何意义，稀疏矩阵实际上并不会将整个矩阵存储在内存中；它只存储那些由项目占用的单元格。这使得该结构比同等大小的矩阵或数据框更节省内存。

为了从交易数据中创建稀疏矩阵数据结构，我们可以使用arules包提供的功能。通过执行install.packages("arules")和library(arules)命令来安装和加载该包。

注意

欲了解更多有关 arules 包的信息，请参考：Hahsler M, Gruen B, Hornik K. arules – a computational environment for mining association rules and frequent item sets. Journal of Statistical Software. 2005; 14。

由于我们正在加载交易数据，不能简单地使用之前的read.csv()函数。相反，arules提供了一个read.transactions()函数，类似于read.csv()，但其结果是适用于交易数据的稀疏矩阵。sep = ","参数指定输入文件中的项目由逗号分隔。要将groceries.csv数据读入名为groceries的稀疏矩阵，输入以下代码：

> groceries <- read.transactions("groceries.csv", sep = ",")

要查看我们刚刚创建的groceries矩阵的基本信息，可以对该对象使用summary()函数：

> summary(groceries)
transactions as itemMatrix in sparse format with
 9835 rows (elements/itemsets/transactions) and
 169 columns (items) and a density of 0.02609146

输出中的第一部分信息（如前所示）提供了我们创建的稀疏矩阵的摘要。输出中的9835 rows表示交易的数量，169 columns表示可能出现在购物篮中的 169 个不同商品。矩阵中的每个单元格，如果相应的交易中购买了该商品，则为1，否则为0。

0.02609146（2.6%）的密度值表示矩阵中非零单元格的比例。由于矩阵中共有9,835 * 169 = 1,662,115个位置，我们可以计算出在商店的 30 天运营期间，共购买了1,662,115 * 0.02609146 = 43,367个商品（忽略了同一商品可能重复购买的情况）。通过额外的步骤，我们可以确定平均每笔交易包含了43,367 / 8,835 = 4.409个不同的杂货商品。当然，如果我们进一步查看输出结果，会看到每笔交易的商品平均数已经在输出中给出。

summary()输出的下一个区块列出了在交易数据中最常见的项目。由于2,513 / 9,835 = 0.2555，我们可以确定全脂牛奶出现在 25.6％的交易中。其他常见的项目包括其他蔬菜、卷/小圆面包、汽水和酸奶，具体如下：

most frequent items:
 whole milk other vegetables       rolls/buns
 2513             1903             1809
 soda           yogurt          (Other)
 1715             1372            34055

最后，我们展示了一组关于交易大小的统计信息。共有 2,159 个交易仅包含单一项目，而有一个交易包含 32 个项目。第一四分位数和中位数购买大小分别为 2 和 3 个项目，这意味着 25％的交易包含两个或更少的项目，并且交易在包含少于三个项目的和包含更多项目的交易之间平均分配。每个交易 4.409 个项目的平均值与我们手动计算的结果一致。

element (itemset/transaction) length distribution:
sizes
 1    2    3    4    5    6    7    8    9   10   11   12
2159 1643 1299 1005  855  645  545  438  350  246  182  117
 13   14   15   16   17   18   19   20   21   22   23   24
 78   77   55   46   29   14   14    9   11    4    6    1
 26   27   28   29   32
 1    1    1    3    1

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 1.000   2.000   3.000   4.409   6.000  32.000

arules包包含一些有用的功能，用于检查交易数据。要查看稀疏矩阵的内容，请将inspect()函数与向量操作符结合使用。可以按如下方式查看前五个交易：

> inspect(groceries[1:5])
 items 
1 {citrus fruit, 
 margarine, 
 ready soups, 
 semi-finished bread} 
2 {coffee, 
 tropical fruit, 
 yogurt} 
3 {whole milk} 
4 {cream cheese, 
 meat spreads, 
 pip fruit, 
 yogurt} 
5 {condensed milk, 
 long life bakery product,
 other vegetables, 
 whole milk}

这些交易与我们查看原始 CSV 文件时的数据一致。要检查特定的项目（即一列数据），可以使用[row, column]矩阵概念。将其与itemFrequency()函数结合使用，可以查看包含该项目的交易比例。例如，这使我们能够查看groceries数据中前 3 个项目的支持度：

> itemFrequency(groceries[, 1:3])
abrasive cleaner artif. sweetener   baby cosmetics
 0.0035587189     0.0032536858     0.0006100661

请注意，稀疏矩阵中的项目按字母顺序按列排序。磨砂清洁剂和人工甜味剂出现在大约 0.3％的交易中，而婴儿化妆品出现在大约 0.06％的交易中。

可视化项目支持度 - 项目频率图

为了直观展示这些统计信息，可以使用itemFrequencyPlot()函数。这允许你生成一张条形图，显示包含特定项目的交易比例。由于交易数据包含大量项目，你通常需要限制显示在图表中的项目，以便生成清晰的图表。

如果你希望这些项目出现在交易的最低比例中，可以使用带有support参数的itemFrequencyPlot()：

> itemFrequencyPlot(groceries, support = 0.1)

如下图所示，这导致了一个显示groceries数据中至少有 10％支持度的八个项目的直方图：

如果你希望限制图表中的项目数量，可以使用带有topN参数的itemFrequencyPlot()：

> itemFrequencyPlot(groceries, topN = 20)

然后，直方图按支持度递减排序，如下图所示，展示了groceries数据中的前 20 个项目：

可视化交易数据 - 绘制稀疏矩阵

除了查看商品外，还可以可视化整个稀疏矩阵。为此，请使用image()函数。显示前五笔交易的稀疏矩阵的命令如下：

> image(groceries[1:5])

结果图表描绘了一个包含 5 行和 169 列的矩阵，表示我们请求的 5 笔交易和 169 个可能的商品。矩阵中的单元格用黑色填充，表示在某笔交易（行）中购买了某个商品（列）。

尽管前面的图表较小且可能略显难以阅读，但你可以看到第一、第四和第五行的交易每行都包含四个项目，因为它们的行中有四个单元格被填充。你还可以看到第三、第五、第二和第四行在右侧图表中有一个共同的项目。

这种可视化可以成为探索数据的有用工具。首先，它可能有助于识别潜在的数据问题。完全填充的列可能表明在每一笔交易中都有购买某个商品——这可能是一个问题，例如，如果零售商的名称或识别号不小心包含在交易数据集中。

此外，图表中的模式可能有助于揭示交易和商品中的可操作见解，尤其是当数据以有趣的方式排序时。例如，如果交易按日期排序，黑色点的模式可能揭示购买的商品数量或类型的季节性变化。或许在圣诞节或光明节时，玩具更为常见；而在万圣节时，糖果可能变得受欢迎。如果商品也被分类，这种可视化可能尤其强大。然而，在大多数情况下，图表看起来将相当随机，就像电视屏幕上的雪花。

请记住，这种可视化对于极大的交易数据库来说并不那么有用，因为单元格会小到难以辨认。尽管如此，通过将其与sample()函数结合使用，你可以查看一个随机抽样的交易集的稀疏矩阵。创建 100 笔交易随机选择的命令如下：

> image(sample(groceries, 100))

这将生成一个有 100 行和 169 列的矩阵图：

一些列看起来填充得相当密集，表明商店里有一些非常受欢迎的商品。但总体来说，点的分布似乎相当随机。既然没有其他值得注意的内容，我们继续进行分析。

步骤 3 – 在数据上训练模型

数据准备完成后，我们现在可以开始寻找购物车物品之间的关联了。我们将使用arules包中的 Apriori 算法实现，这个包我们之前已经用来探索和准备杂货数据。如果你还没有安装和加载这个包，记得先操作。以下表格显示了使用apriori()函数创建规则集的语法：

尽管运行apriori()函数非常直接，但有时可能需要一定的试错过程来找到产生合理数量关联规则的support和confidence参数。如果你设定这些值过高，你可能找不到规则，或者找到的规则过于宽泛，无法提供实际价值。另一方面，阈值设定得过低可能会导致规则数量过多，甚至更糟，可能会导致学习阶段运行时间过长或内存不足。

在这种情况下，如果我们尝试使用默认设置support = 0.1和confidence = 0.8，最终将得到一组零规则：

> apriori(groceries)
set of 0 rules

显然，我们需要稍微扩大一下搜索范围。

提示

如果你仔细想一想，这个结果其实并不令人惊讶。因为support = 0.1是默认设置，要生成规则，一个项必须至少出现在0.1 * 9,385 = 938.5次交易中。由于在我们的数据中只有八个项出现得如此频繁，因此我们没能找到任何规则也就不足为奇了。

设定最小支持度阈值的一种方法是思考在什么样的最小交易数下，你会认为一个模式是有趣的。例如，你可以认为，如果某个项每天被购买两次（在一个月的数据中大约是 60 次），它可能是一个有趣的模式。从这里开始，你可以计算出需要的支持度水平，以便仅找到匹配至少那么多次交易的规则。由于 60/9,835 等于 0.006，我们首先尝试将支持度设为这个值。

设置最小置信度是一个微妙的平衡。一方面，如果置信度太低，我们可能会被大量不可靠的规则所淹没——例如，几十个规则显示电池通常会与其他物品一起购买。那么，我们怎么知道该将广告预算投向哪里呢？另一方面，如果我们将置信度设定得过高，那么规则将仅限于那些显而易见或不可避免的情况——类似于烟雾探测器总是与电池一起购买的事实。在这种情况下，将烟雾探测器放得更靠近电池可能不会带来额外的收益，因为这两种物品已经几乎总是一起购买了。

提示

适当的最小置信度水平在很大程度上取决于你分析的目标。如果你从保守的值开始，若未找到可操作的情报，你总是可以减少置信度，扩大搜索范围。

我们将从 0.25 的置信度阈值开始，这意味着为了被包括在结果中，规则必须至少正确 25% 的时间。这将剔除最不可靠的规则，同时为我们通过有针对性的促销调整行为提供一定的空间。

我们现在可以生成一些规则。除了最小的 support 和 confidence 参数外，设置 minlen = 2 也很有帮助，以排除那些包含少于两个项目的规则。这可以防止一些无趣的规则被创建，比如仅因为某个项目购买频繁而生成的规则，如 {} → whole milk。这条规则符合最小支持度和置信度，因为全脂牛奶在超过 25% 的交易中被购买，但它并不是一个非常有行动价值的洞察。

使用 Apriori 算法查找关联规则的完整命令如下：

> groceryrules <- apriori(groceries, parameter = list(support =
 0.006, confidence = 0.25, minlen = 2))

这将我们的规则保存在 rules 对象中，可以通过键入其名称来查看：

> groceryrules
set of 463 rules

我们的 groceryrules 对象包含了一组 463 条关联规则。为了确定它们是否有用，我们需要深入挖掘。

步骤 4 – 评估模型性能

为了获得关联规则的高层概览，我们可以按如下方式使用 summary()。规则长度分布告诉我们每条规则包含的项目数量。在我们的规则集中，150 条规则只有两个项目，297 条规则有三个项目，16 条规则有四个项目。与此分布相关的 summary 统计数据也会提供：

> summary(groceryrules)
set of 463 rules

rule length distribution (lhs + rhs):sizes
 2   3   4
150 297  16

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 2.000   2.000   3.000   2.711   3.000   4.000

提示

如前述输出所示，规则的大小是通过计算规则的左侧（lhs）和右侧（rhs）的总和来得出的。这意味着像 {bread} → {butter} 这样的规则是两个项，而 {peanut butter, jelly} → {bread} 是三个。

接下来，我们看到规则质量度量的汇总统计数据：支持度、置信度和提升度。支持度和置信度的度量不应该令人惊讶，因为我们使用这些作为规则的选择标准。如果大多数或所有规则的支持度和置信度接近最低阈值，我们可能会感到警觉，因为这意味着我们可能设置的标准过高。但这里并非如此，因为有许多规则的每个值都远高于这些最低值。

summary of quality measures:
 support           confidence          lift 
 Min.   :0.006101   Min.   :0.2500   Min.   :0.9932 
 1st Qu.:0.007117   1st Qu.:0.2971   1st Qu.:1.6229 
 Median :0.008744   Median :0.3554   Median :1.9332 
 Mean   :0.011539   Mean   :0.3786   Mean   :2.0351 
 3rd Qu.:0.012303   3rd Qu.:0.4495   3rd Qu.:2.3565 
 Max.   :0.074835   Max.   :0.6600   Max.   :3.9565 

第三列是我们尚未考虑的度量。规则的 lift 衡量的是在已知另一个项目或项目集被购买的情况下，某个项目或项目集的购买可能性相对于其通常购买率的提高程度。这个值通过以下方程定义：

提示

与置信度不同，提升度 lift(X → Y) 与 lift(Y → X) 是相同的。

例如，假设在一家超市里，大多数人都购买牛奶和面包。仅凭偶然，我们就可以预期会有许多同时购买牛奶和面包的交易。然而，如果 lift(牛奶 → 面包) 大于 1，这意味着这两种商品比偶然情况下更常一起出现。因此，较高的提升值是规则重要性的强烈指示，反映了商品之间的真实联系。

在 summary() 输出的最后部分，我们会看到挖掘信息，告诉我们规则是如何选择的。在这里，我们看到包含 9,835 个交易的 groceries 数据被用来构建规则，最小支持度为 0.0006，最小置信度为 0.25：

mining info:
 data ntransactions support confidence
 groceries          9835   0.006       0.25

我们可以使用 inspect() 函数查看特定规则。例如，可以通过以下方式查看 groceryrules 对象中的前三条规则：

> inspect(groceryrules[1:3])

第一个规则可以用简单的语言表述为：“如果顾客购买盆栽植物，他们也会购买全脂牛奶。”通过支持度为 0.007 和置信度为 0.400，我们可以确定这个规则覆盖了 0.7% 的交易，并且在涉及盆栽植物的购买中正确率为 40%。提升值告诉我们，相比于平均顾客，购买盆栽植物的顾客购买全脂牛奶的可能性有多大。由于我们知道大约 25.6% 的顾客购买了全脂牛奶（support），而 40% 购买盆栽植物的顾客购买了全脂牛奶（confidence），因此我们可以计算提升值为 0.40 / 0.256 = 1.56，这与显示的值一致。

注意

请注意，标有 support 的列表示的是规则的支持值，而不是 lhs 或 rhs 的支持值单独的支持值。

尽管置信度和提升值很高，{盆栽植物} → {全脂牛奶} 看起来仍然像是一个非常有用的规则吗？可能不是，因为似乎没有逻辑理由说明某人购买盆栽植物时更有可能购买牛奶。然而，我们的数据却暗示了不同的结果。我们如何解释这一事实呢？

一种常见的方法是将关联规则分为以下三类：

• 可操作的

• 微不足道的

• 无法解释的

显然，市场篮分析的目标是找到可操作的规则，这些规则能提供明确且有用的洞察。有些规则很明确，有些规则很有用；同时具备这两者的组合则较为少见。

所谓的微不足道规则是指那些过于明显、没有价值的规则——它们很清晰，但没有用处。假设你是一名市场咨询师，受雇于客户，获取丰厚的报酬去识别新的交叉销售机会。如果你报告发现 {纸尿裤} → {奶粉}，你可能不会被邀请回来进行下一次咨询工作。

提示

微不足道的规则也可能伪装成更有趣的结果。例如，假设你发现某个特定品牌的儿童谷物和某部 DVD 电影之间存在关联。如果这部电影的主角出现在谷物盒的封面上，那么这个发现就没有太大意义。

如果规则之间的关联不明确，导致无法或几乎无法理解如何使用这些信息，则该规则是无法解释的。这个规则可能仅仅是数据中的一个随机模式，例如，一个规则表明*{腌黄瓜}* → {巧克力冰淇淋}，可能是因为某个顾客的孕妇妻子经常对奇怪的食物组合有特殊的渴望。

最佳规则是隐藏的宝石——那些似乎一旦被发现就显而易见的模式洞察。给定足够的时间，可以评估每一条规则来发现这些宝石。然而，我们（执行市场购物篮分析的人）可能并不是判断规则是否具有可操作性、微不足道或无法解释的最佳人选。在接下来的部分中，我们将通过使用排序和共享已学规则的方法来提高我们的工作效用，从而使最有趣的结果能够浮到最上面。

第 5 步 – 提高模型表现

专家可能能够非常迅速地识别有用的规则，但让他们评估成百上千的规则将是浪费他们时间。因此，能够根据不同标准对规则进行排序，并将它们从 R 中导出，以便与营销团队共享并深入研究是很有用的。通过这种方式，我们可以通过使结果更具可操作性来提高规则的表现。

排序关联规则集合

根据市场购物篮分析的目标，最有用的规则可能是那些具有最高support、confidence或lift值的规则。arules包包含一个sort()函数，可以用来重新排序规则列表，使得具有最高或最低质量度量值的规则排在最前面。

要重新排序groceryrules对象，我们可以在by参数中指定一个"support"、"confidence"或"lift"值并应用sort()函数。通过将sort函数与向量操作符结合使用，我们可以获得一个特定数量的有趣规则。例如，按照 lift 统计值排序，可以使用以下命令检查最佳的五条规则：

> inspect(sort(groceryrules, by = "lift")[1:5])

输出如下所示：

这些规则似乎比我们之前查看的规则更有趣。第一条规则的lift约为 3.96，意味着购买香草的人比普通顾客更可能购买根菜类蔬菜，其可能性是普通顾客的近四倍——也许是为了做某种炖菜？第二条规则也很有趣。打发奶油在购物车中与浆果一起出现的概率是其他购物车的三倍多，这或许暗示了一种甜点搭配？

提示

默认情况下，排序顺序是降序的，意味着最大的值排在前面。要反转此顺序，只需添加一行 parameterdecreasing = FALSE。

获取关联规则的子集

假设根据前面的规则，营销团队对创建一个广告来推广当前正当季的浆果感到兴奋。然而，在最终确定活动之前，他们要求你调查一下浆果是否经常与其他商品一起购买。为了回答这个问题，我们需要找出所有包含浆果的规则。

subset() 函数提供了一种方法来查找交易、项目或规则的子集。要使用它查找任何包含 berries 的规则，可以使用以下命令。它会将规则存储在一个名为 berryrules 的新对象中：

> berryrules <- subset(groceryrules, items %in% "berries")

然后，我们可以像处理更大数据集时那样检查这些规则：

> inspect(berryrules)

结果是以下一组规则：

涉及浆果的规则有四条，其中两条似乎足够有趣，可以称为可执行的。除了鲜奶油，浆果也经常与酸奶一起购买——这种搭配既适合早餐或午餐，也可以作为甜点。

subset() 函数非常强大。选择子集的标准可以通过多个关键字和操作符来定义：

• 前面解释过的关键字 items 匹配规则中任何位置出现的项目。为了将子集限制为只匹配出现在左侧或右侧的项目，可以改用 lhs 和 rhs。

• 操作符 %in% 表示在你定义的列表中，至少有一个项目必须存在。如果你想要匹配浆果或酸奶的任何规则，可以写 items %in%c("berries", "yogurt")。

• 还提供了额外的操作符，用于部分匹配（%pin%）和完全匹配（%ain%）。部分匹配允许你通过一个查询同时找到柑橘类水果和热带水果：items %pin% "fruit"。完全匹配要求所有列出的项目都必须出现。例如，items %ain% c("berries", "yogurt") 只会找到包含 berries 和 yogurt 的规则。

• 子集还可以通过 support、confidence 或 lift 来限制。例如，confidence > 0.50 将限制为那些置信度大于 50% 的规则。

• 匹配标准可以与标准的 R 逻辑运算符结合使用，如与（&）、或（|）和非（!）。

使用这些选项，你可以将规则的选择限制得尽可能具体或宽泛。

将关联规则保存到文件或数据框

为了分享你的市场购物篮分析结果，你可以使用 write() 函数将规则保存为 CSV 文件。这将生成一个可以在大多数电子表格程序（包括 Microsoft Excel）中使用的 CSV 文件：

> write(groceryrules, file = "groceryrules.csv",
 sep = ",", quote = TRUE, row.names = FALSE)

有时，将规则转换为 R 数据框也是很方便的。可以通过 as() 函数轻松实现，如下所示：

> groceryrules_df <- as(groceryrules, "data.frame")

这将创建一个数据框，其中包含按因子格式表示的规则，以及支持度、置信度和提升度的数值向量：

> str(groceryrules_df)
'data.frame':463 obs. of 4 variables:
 $ rules     : Factor w/ 463 levels "{baking powder} => {other vegetables}",..: 340 302 207 206 208 341 402 21 139 140 ...
 $ support   : num  0.00691 0.0061 0.00702 0.00773 0.00773 ...
 $ confidence: num  0.4 0.405 0.431 0.475 0.475 ...
 $ lift      : num  1.57 1.59 3.96 2.45 1.86 ...

如果你想对规则执行进一步处理或需要将它们导出到另一个数据库中，你可以选择这样做。

总结

关联规则通常用于在大型零售商的海量交易数据库中提供有用的洞察。作为一种无监督学习过程，关联规则学习器能够从大型数据库中提取知识，而无需提前了解要寻找什么样的模式。问题在于，它需要一定的努力才能将大量信息缩减为一个更小、更易管理的结果集。我们在本章中研究的 Apriori 算法通过设置最小兴趣阈值，并仅报告满足这些标准的关联，从而实现了这一目标。

我们在进行一个中型超市一个月交易数据的市场购物篮分析时，运用了 Apriori 算法。即使在这个小示例中，也发现了大量的关联。在这些关联中，我们注意到了一些可能对未来营销活动有用的模式。我们应用的相同方法也被用于更大规模的零售商，分析的数据集大小是这些的多倍。

在下一章中，我们将研究另一种无监督学习算法。与关联规则类似，它也旨在在数据中找到模式。但与寻找特征间模式的关联规则不同，下一章中的方法更关注在示例之间发现联系。

第九章：数据分组—使用 k-means 进行聚类

您是否曾花时间观察过一大群人？如果有，您可能已经见过一些反复出现的个性。也许某种类型的人，穿着刚熨好的西装，手拿公文包，成为了典型的“肥猫”商界高管。一个穿着紧身牛仔裤、法兰绒衬衫和太阳镜的二十多岁年轻人可能被称为“嬉皮士”，而一位从小面包车里抱出孩子的女性则可能被标记为“足球妈妈”。

当然，这些刻板印象在应用到个体时是危险的，因为没有两个完全相同的人。然而，作为描述一个集体的方式，这些标签捕捉到了群体内个体之间某些潜在的相似性。

正如您很快会学到的，聚类或在数据中识别模式的过程，与在人群中识别模式并没有太大不同。在这一章中，您将学到：

• 聚类任务与我们之前研究的分类任务有何不同

• 聚类如何定义一个组，以及如何通过 k-means 这一经典且易于理解的聚类算法识别这些组

• 将聚类应用于现实世界任务的步骤，比如在青少年社交媒体用户中识别营销细分

在开始实际操作之前，我们将深入探讨聚类的具体内容。

理解聚类

聚类是一种无监督的机器学习任务，它会自动将数据分成聚类，即一组组相似的项目。它在不知道如何预先划分组的情况下进行。这就像我们可能甚至不知道自己在寻找什么，聚类用于知识发现而非预测。它为我们提供了数据中自然分组的洞察。

在没有预先了解一个聚类由什么组成的情况下，计算机如何可能知道一个组何时结束，另一个组何时开始呢？答案很简单。聚类是由一个原则指导的：聚类内部的项应该彼此非常相似，而与聚类外部的项有很大不同。相似性的定义可能会因应用场景的不同而有所变化，但基本思想始终相同——将数据分组，使相关元素聚集在一起。

生成的聚类可以用于实际操作。例如，您可能会发现聚类方法应用于以下场景：

• 将客户按相似的人口统计信息或购买模式进行分组，以进行定向营销活动

• 通过识别使用模式，发现超出已知聚类范围的异常行为，例如未经授权的网络入侵

• 通过将具有相似值的特征分组为更少的同质类别，从而简化极其庞大的数据集

总的来说，聚类在数据多样且可以通过更少的组来概括时非常有用。它能产生有意义且可操作的数据结构，减少复杂性并为关系模式提供洞察。

聚类作为机器学习任务

聚类与我们之前探讨的分类、数值预测和模式检测任务有所不同。在这些任务中，结果是一个模型，将特征与结果或特征之间的关系进行关联；从概念上讲，模型描述了数据中的现有模式。相比之下，聚类会创造新的数据。未标记的示例会被赋予一个聚类标签，该标签完全是根据数据内部的关系推断出来的。因此，有时你会看到聚类任务被称为无监督分类，因为从某种意义上来说，它是对未标记的示例进行分类。

关键在于，来自无监督分类器的类别标签没有内在意义。聚类将告诉你哪些示例组之间有紧密的关系——例如，它可能会返回 A、B 和 C 组——但你需要为这些组应用一个可操作且有意义的标签。为了了解这对聚类任务的影响，我们来看一个假设的例子。

假设你正在组织一个关于数据科学的会议。为了促进专业的社交和合作，你计划根据三种研究专长之一将与会者分组：计算机和/或数据库科学、数学和统计学、以及机器学习。不幸的是，在发送会议邀请后，你意识到忘记了包含一项调查，询问与会者希望与哪个学科的人员坐在一起。

一时灵感迸发，你意识到可以通过检查每个学者的出版历史来推断他们的研究专长。为此，你开始收集每个与会者在计算机科学相关期刊上发表的文章数量，以及在数学或统计学相关期刊上发表的文章数量。通过收集几位学者的数据，你绘制了一个散点图：

正如预期的那样，似乎确实存在一个模式。我们可能会猜测，位于左上角的群体，代表那些在计算机科学方面有许多出版物，但在数学方面文章较少的人，可能是计算机科学家的聚类。按照这一逻辑，右下角可能是数学家的群体。同样，右上角那些既有数学又有计算机科学经验的人，可能是机器学习专家。

我们的分组是通过视觉方式形成的；我们只是简单地将数据点作为紧密聚集的群体进行识别。然而，尽管这些分组看起来显而易见，但不幸的是，我们没有办法知道这些分组是否真正同质化，因为我们无法逐个询问每位学者的学术专长。我们所应用的标签要求我们对可能属于该组的人进行定性、假设性的判断。因此，你可以把群体标签理解为不确定的术语，如下所示：

与其主观地定义群体边界，不如使用机器学习来客观地定义它们。考虑到前面图中的轴向分割，我们的问题似乎是第五章中所描述的决策树的明显应用，分治法 – 使用决策树和规则进行分类。这可能会为我们提供一个规则，形式为“如果学者的数学出版物较少，那么他/她就是计算机科学专家。”不幸的是，这个计划存在问题。由于我们没有每个点的真实类别数据，监督学习算法无法学习到这样的模式，因为它无法知道哪些分割能形成同质化的群体。

另一方面，聚类算法使用的过程与我们通过视觉检查散点图所做的非常相似。通过衡量样本之间的关系程度，可以识别出同质的群体。在接下来的部分中，我们将开始探讨聚类算法是如何实现的。

提示

这个例子突出了聚类的一个有趣应用。如果你从未标记的数据开始，你可以使用聚类来创建类别标签。从那里，你可以应用监督学习算法，如决策树，来找出这些类别最重要的预测因素。这被称为半监督学习。

k-means 聚类算法

k-means 算法可能是最常用的聚类方法。经过数十年的研究，它为许多更复杂的聚类技术奠定了基础。如果你理解了它使用的简单原理，你就具备了理解今天几乎所有聚类算法所需的知识。许多此类方法列在以下网站上，即CRAN 聚类任务视图：cran.r-project.org/web/views/Cluster.html。

注意

随着 k-means 的不断发展，该算法有许多不同的实现。一个常见的方法描述在：Hartigan JA, Wong MA. A k-means clustering algorithm. Applied Statistics. 1979; 28:100-108。

尽管聚类方法自 k-means 诞生以来已有了进展，但这并不意味着 k-means 已经过时。事实上，这种方法现在可能比以往任何时候都更受欢迎。以下表格列出了一些 k-means 仍然广泛使用的原因：

优势	弱点

|

• 使用简单的原则，可以用非统计学术语进行解释

• 高度灵活，可以通过简单的调整适应，解决几乎所有的缺点

• 在许多实际应用场景中表现良好

|

• 比不上更现代的聚类算法复杂

• 由于它使用了随机因素，因此无法保证找到最优的聚类集

• 需要合理猜测数据中自然存在多少个聚类

• 对于非球形聚类或密度差异较大的聚类不理想

|

如果“k-means”这个名字让你觉得熟悉，你可能是在回忆第三章中讨论的 k-NN 算法，懒学习 – 使用最近邻分类。正如你将很快看到的，k-means 与 k 最近邻的相似之处远不止字母 k。

k-means 算法将每个n个样本分配到k个聚类中的一个，k是一个事先确定的数字。目标是最小化每个聚类内的差异，并最大化聚类之间的差异。

除非k和n非常小，否则计算所有可能组合的最优聚类是不可行的。相反，算法使用启发式过程来寻找局部最优解。简单来说，这意味着它从初始聚类分配开始，然后稍微修改分配，看看这些变化是否能提高聚类内部的一致性。

我们将很快深入讨论这个过程，但该算法基本上分为两个阶段。首先，它将样本分配到初始的k个聚类中。然后，它通过根据当前落入聚类中的样本调整聚类边界来更新分配。更新和分配的过程会重复几次，直到变化不再改善聚类拟合为止。此时，过程停止，聚类最终确定。

提示

由于 k-means 的启发式特性，仅通过轻微改变初始条件，你可能会得到略有不同的最终结果。如果结果差异很大，这可能意味着存在问题。例如，数据可能没有自然分组，或者k的值选择不当。考虑到这一点，最好多次尝试聚类分析，以测试结果的稳健性。

为了了解分配和更新的过程如何在实践中工作，让我们重新审视一下假设的数据科学会议的案例。虽然这是一个简单的例子，但它将展示 k-means 在背后如何运作的基础。

使用距离分配和更新聚类

与 k-NN 一样，k-means 将特征值视为多维特征空间中的坐标。对于这次会议数据，只有两个特征，因此我们可以将特征空间表示为如前所示的二维散点图。

k-means 算法首先通过在特征空间中选择 k 个点作为聚类中心。这些中心是推动剩余样本归类的催化剂。通常，这些点通过从训练数据集中随机选择 k 个样本来确定。假设我们希望识别三个聚类，按照这种方法，k = 3 个点将被随机选择。这些点在下图中分别用星号、三角形和菱形表示：

值得注意的是，尽管前图中的三个聚类中心恰好相距较远，但这并不一定总是如此。由于它们是随机选择的，这三个中心也有可能是相邻的三个点。由于 k-means 算法对聚类中心的初始位置高度敏感，这意味着随机因素可能对最终的聚类结果产生重大影响。

为了解决这个问题，k-means 可以通过不同的方法来选择初始中心。例如，一种变体会在特征空间中的任意位置选择随机值（而不仅仅是从数据中选择已观察到的值）。另一种选择是完全跳过这一步；通过将每个样本随机分配到一个聚类中，算法可以立即跳到更新阶段。每种方法都会对最终的聚类集添加特定的偏差，你可能可以利用这些偏差来改进你的结果。

注意

2007 年，提出了一种名为 k-means++ 的算法，它提供了一种选择初始聚类中心的替代方法。该方法声称是一种高效的方式，能够在减少随机因素影响的同时，接近最优的聚类结果。欲了解更多信息，请参阅 Arthur D, Vassilvitskii S。k-means++: 精心初始化的优势。第十八届 ACM-SIAM 离散算法年会论文集，2007：1027–1035。

选择初始聚类中心后，其他样本会根据距离函数被分配到最近的聚类中心。你应该还记得我们在学习 k-最近邻时研究了距离函数。传统上，k-means 使用欧几里得距离，但也有时使用曼哈顿距离或闵可夫斯基距离。

回想一下，如果 n 表示特征的数量，则示例 x 和示例 y 之间的欧几里得距离公式为：

例如，如果我们要比较一位有五篇计算机科学论文和一篇数学论文的访客与一位没有计算机科学论文但有两篇数学论文的访客，我们可以在 R 中按如下方式计算：

> sqrt((5 - 0)² + (1 - 2)²)
[1] 5.09902

使用该距离函数，我们计算每个示例与每个聚类中心之间的距离。然后，示例被分配到最近的聚类中心。

提示

请记住，由于我们正在使用距离计算，所有特征需要是数值型的，且值应提前标准化到一个标准范围。第三章中讨论的方法，惰性学习——使用最近邻分类，对这个任务会很有帮助。

如下图所示，三个聚类中心将示例分为三个部分，分别标记为聚类 A、聚类 B和聚类 C。虚线表示由聚类中心创建的Voronoi 图的边界。Voronoi 图显示了离某个聚类中心比其他任何中心都近的区域；三个边界相交的顶点是离所有三个聚类中心最远的点。利用这些边界，我们可以轻松看到每个初始 k-means 种子所占据的区域：

现在初始分配阶段已经完成，k-means 算法进入更新阶段。更新聚类的第一步是将初始中心移动到新的位置，称为质心，其位置是当前分配给该聚类的点的平均位置。下图展示了当聚类中心移动到新的质心时，Voronoi 图中的边界也会发生变化，而原本在聚类 B中的一个点（由箭头表示）被加入到了聚类 A中：

由于此次重新分配，k-means 算法将继续进行另一个更新阶段。在移动聚类中心、更新聚类边界并重新分配点到新聚类之后（如箭头所示），图形如下所示：

由于又有两个点被重新分配，因此必须进行另一次更新，这将移动中心并更新聚类边界。然而，由于这些变化没有导致任何重新分配，k-means 算法停止。聚类分配现在是最终的：

最终的聚类结果可以通过两种方式之一来报告。首先，你可以简单地报告每个样本的聚类分配，如 A、B 或 C。或者，你可以报告最终更新后聚类中心的坐标。无论哪种报告方式，你都可以通过计算中心点或将每个样本分配给其最近的聚类来定义聚类边界。

选择适当数量的聚类

在 k-means 的介绍中，我们了解到该算法对随机选择的聚类中心非常敏感。实际上，如果我们在前一个例子中选择了不同的三个起始点，可能会得到与我们预期不同的数据分组。同样，k-means 对聚类数量也非常敏感；这个选择需要一个微妙的平衡。将 k 设置得非常大将提高聚类的同质性，但同时也有过拟合数据的风险。

理想情况下，你应该具有 先验 知识（即先前的信念）关于真实的分组情况，并可以利用这些信息来选择聚类数量。例如，如果你在对电影进行聚类，你可能会将 k 设置为奥斯卡奖提名的类型数量。在我们之前讨论的数据科学会议座位问题中，k 可能反映了受邀的学术领域数量。

有时候，聚类的数量由业务需求或分析的动机决定。例如，会议厅中的桌子数量可能决定了应该从数据科学参与者名单中创建多少个小组。将这个思路扩展到另一个商业案例，如果市场部门只有限制资源来创建三种不同的广告活动，那么将 k = 3 设置为将所有潜在客户分配给三个吸引点中的一个可能是合理的。

在没有任何先验知识的情况下，有一个经验法则建议将 k 设置为 (n / 2) 的平方根，其中 n 是数据集中的样本数量。然而，对于大型数据集来说，这个经验法则可能会导致聚类数量过多而难以处理。幸运的是，还有其他统计方法可以帮助找到合适的 k-means 聚类集。

一种名为 肘部法则 的技术试图衡量在不同的 k 值下，聚类内的同质性或异质性是如何变化的。正如下图所示，随着聚类数量的增加，聚类内部的同质性预期会增加；同样，异质性也会随着聚类数量的增加而继续减少。尽管你可以继续观察到每个样本被分配到其自己的聚类，目标不是最大化同质性或最小化异质性，而是找到一个 k，在这个点之后，收益递减。这种 k 值被称为 肘部点，因为它看起来像一个肘部。

有许多统计方法可以衡量聚类内部的同质性和异质性，这些方法可以与肘部法（以下信息框提供了详细的引用）一起使用。然而，在实际应用中，并不总是可行的反复测试大量k值。部分原因是，聚类大数据集本身就可能非常耗时；而反复进行数据聚类则更加浪费时间。无论如何，需要精确的最优聚类集的应用相对较少。在大多数聚类应用中，选择一个方便的k值而非严格的性能要求通常就足够了。

注意

对于关于聚类性能度量的大量综述，请参考：Halkidi M, Batistakis Y, Vazirgiannis M。关于聚类验证技术。智能信息系统杂志。2001；17:107-145。

设置k的过程本身有时会带来有趣的洞见。通过观察随着k变化，聚类的特征如何变化，可能会推测出数据自然的边界所在。聚集得更紧密的组变化较小，而异质性较大的组则会随着时间的推移不断形成和解散。

一般来说，花费很少的时间去担心k是否完全准确是明智的。下一个例子将展示，即使是从一部好莱坞电影中借来的一点点专业知识，也能用于设定k，以便发现可操作且有趣的聚类。由于聚类是无监督的，任务的本质实际上是你如何理解它；真正的价值在于从算法发现中获得的洞察。

示例——使用 k-means 聚类寻找青少年市场细分

与朋友在社交网络服务（SNS）上互动，例如 Facebook、Tumblr 和 Instagram，已经成为全球青少年的一项通行仪式。这些青少年通常拥有相对较多的可支配收入，因此成为了企业争相吸引的群体，企业希望通过他们销售零食、饮料、电子产品和卫生用品。

使用这些网站的数百万青少年消费者已经吸引了市场营销人员的关注，他们在日益竞争激烈的市场中努力寻找竞争优势。一种获得这种优势的方法是识别出具有相似品味的青少年群体，以便客户避免将广告投放给那些对所售产品没有兴趣的青少年。例如，运动服饰可能很难成功地销售给对体育不感兴趣的青少年。

给定青少年 SNS 页面的文本，我们可以识别出一些具有共同兴趣的群体，例如体育、宗教或音乐。聚类可以自动化发现这一人群中的自然分段过程。然而，是否认为这些聚类有趣，以及如何利用它们进行广告投放，仍然需要我们自己决定。让我们从头到尾尝试这一过程。

第 1 步 – 收集数据

对于本次分析，我们将使用一个代表 2006 年在一个著名 SNS 上有个人资料的 30,000 名美国高中生的随机样本数据集。为了保护用户的匿名性，SNS 的名称将保持不公开。然而，在数据收集时，这个 SNS 是美国青少年常用的网络平台。因此，可以合理假设这些个人资料代表了 2006 年美国青少年群体的一个较为广泛的横截面。

提示

该数据集是 Brett Lantz 在大学进行青少年身份的社会学研究时编制的。如果你用于研究目的，请引用这本书的章节。完整数据集可在 Packt Publishing 网站上下载，文件名为snsdata.csv。为了进行互动操作，本章假设你已将该文件保存到你的 R 工作目录中。

数据在四个高中毕业年份（2006 年至 2009 年）之间均匀抽样，代表了数据收集时的高年级、低年级、二年级和一年级学生。使用自动化的网络爬虫，下载了 SNS 个人资料的完整文本，并记录了每个青少年的性别、年龄和 SNS 好友数量。

使用了一个文本挖掘工具，将其余的 SNS 页面内容分割成单词。从所有页面中出现的前 500 个单词中，选择了 36 个单词来代表五类兴趣：即课外活动、时尚、宗教、浪漫和反社会行为。选中的 36 个单词包括足球、性感、亲吻、圣经、购物、死亡和毒品等。最终的数据集显示了每个人在其 SNS 个人资料中每个单词出现的次数。

第 2 步 – 探索和准备数据

我们可以使用read.csv()的默认设置将数据加载到数据框中：

> teens <- read.csv("snsdata.csv")

让我们也快速查看一下数据的具体情况。str()输出的前几行如下：

> str(teens)
'data.frame':  30000 obs. of  40 variables:
 $ gradyear    : int  2006 2006 2006 2006 2006 2006 2006 2006 ...
 $ gender      : Factor w/ 2 levels "F","M": 2 1 2 1 NA 1 1 2 ...
 $ age         : num  19 18.8 18.3 18.9 19 ...
 $ friends     : int  7 0 69 0 10 142 72 17 52 39 ...
 $ basketball  : int  0 0 0 0 0 0 0 0 0 0 ...

正如我们预期的那样，数据包括 30,000 名青少年，四个变量表示个人特征，36 个单词表示兴趣。

你是否注意到gender行中的异常情况？如果你仔细查看，可能已经注意到NA值，它与1和2的值不太一致。NA是 R 用来告诉我们记录缺少值的方式——我们不知道这个人的性别。到目前为止，我们还没有处理缺失数据，但对于许多类型的分析来说，缺失数据可能是一个重要问题。

让我们看看这个问题的严重程度。一个选择是使用table()命令，如下所示：

> table(teens$gender)
 F     M
22054  5222

尽管此命令告诉我们有多少F和M值存在，但table()函数排除了NA值，而不是将其视为一个单独的类别。为了包括NA值（如果有的话），我们只需添加一个额外的参数：

> table(teens$gender, useNA = "ifany")
 F     M  <NA>
22054  5222  2724

在这里，我们看到有 2,724 条记录（占 9%）缺少性别数据。有趣的是，SNS 数据中女性的数量是男性的四倍多，这表明男性使用 SNS 网站的倾向不如女性。

如果你检查数据框中的其他变量，你会发现除了gender，只有age存在缺失值。对于数值型数据，summary()命令会告诉我们缺失的NA值的数量：

> summary(teens$age)
 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
 3.086  16.310  17.290  17.990  18.260 106.900    5086

总共有 5,086 条记录（占 17%）缺少年龄数据。令人担忧的是，最小和最大值似乎不合理；一个 3 岁或 106 岁的学生不太可能上高中。为了确保这些极端值不会对分析造成问题，我们需要在继续之前清理它们。

高中生的合理年龄范围应包括至少 13 岁但不满 20 岁的学生。任何超出此范围的年龄值应视为缺失数据——我们无法信任提供的年龄。为了重新编码年龄变量，我们可以使用ifelse()函数，当年龄至少为 13 岁且小于 20 岁时，将teen$age的值设为teen$age；否则，赋值为NA：

> teens$age <- ifelse(teens$age >= 13 & teens$age < 20,
 teens$age, NA)

通过重新检查summary()输出，我们发现现在的年龄范围呈现出一个更像实际高中学生的分布：

> summary(teens$age)
 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
 13.03   16.30   17.26   17.25   18.22   20.00    5523

不幸的是，现在我们制造了一个更大的缺失数据问题。在继续分析之前，我们需要找到一种处理这些缺失值的方法。

数据准备——虚拟编码缺失值

处理缺失值的一个简单方法是排除任何缺失值的记录。然而，如果你考虑这种做法的后果，可能会在做之前三思而后行——仅仅因为这种方法简单，并不意味着它是个好主意！这种方法的问题在于，即使缺失值的数量不多，你也可能会轻易排除大量数据。

例如，假设在我们的数据中，缺失性别值的人群与缺失年龄数据的人群完全不同。这意味着，如果你排除了缺失性别或年龄的记录，你将排除*9% + 17% = 26%*的数据，或者超过 7,500 条记录。而这仅仅是针对两个变量的缺失数据！缺失值数量越多，任何给定记录被排除的可能性就越大。很快，你将只剩下一个非常小的子集数据，或者更糟，剩下的记录将系统性地不同或不具代表性，无法代表整体人群。

对于像性别这样的分类变量，另一种选择是将缺失值视为单独的类别。例如，而不是仅限于女性和男性，我们可以为未知性别添加一个额外的类别。这允许我们使用虚拟编码，这在第三章中有讲解，懒惰学习 – 使用最近邻分类。

如果你还记得，虚拟编码涉及为名义特征的每个级别创建一个单独的二进制（1 或 0）值的虚拟变量，除了一个级别，该级别被保留作为参考组。可以排除一个类别的原因是因为可以从其他类别推断出它的状态。例如，如果某人既不是女性也不是未知性别，他们必须是男性。因此，在这种情况下，我们只需为女性和未知性别创建虚拟变量：

> teens$female <- ifelse(teens$gender == "F" &
 !is.na(teens$gender), 1, 0)
> teens$no_gender <- ifelse(is.na(teens$gender), 1, 0)

正如你可能预期的那样，is.na() 函数测试性别是否等于 NA。因此，第一条语句如果性别等于 F 且性别不等于 NA，则为 teens$female 赋值 1；否则，赋值 0。在第二条语句中，如果 is.na() 返回 TRUE，表示性别缺失，则将 teens$no_gender 变量赋值为 1；否则，赋值为 0。为了确认我们的工作是否正确，让我们将构建的虚拟变量与原始的 gender 变量进行比较：

> table(teens$gender, useNA = "ifany")
 F     M  <NA>
22054  5222  2724
> table(teens$female, useNA = "ifany")
 0     1
 7946 22054
> table(teens$no_gender, useNA = "ifany")
 0     1
27276  2724

对于 teens$female 和 teens$no_gender 中的 1 值数量与 F 和 NA 值的数量匹配，因此我们应该能够信任我们的工作。

数据准备 – 填补缺失值

接下来，让我们排除 5,523 个缺失的年龄值。由于年龄是数值型的，为未知值创建一个额外的类别没有意义——相对于其他年龄，你会将"未知"排名在哪里呢？相反，我们将使用一种称为插补的不同策略，它涉及用真实值的猜测填补缺失数据。

你能想到我们如何能够利用 SNS 数据来推断青少年的年龄吗？如果你考虑使用毕业年份，那么你有正确的想法。一个毕业队列中的大多数人在一个日历年内出生。如果我们能够确定每个队列的典型年龄，我们将有一个相当合理的估计来描述该毕业年份的学生的年龄。

找到典型值的一种方法是计算平均值或均值。如果我们尝试应用 mean() 函数，就像我们之前分析过的那样，会有一个问题：

> mean(teens$age)
[1] NA

问题在于对包含缺失数据的向量计算均值是未定义的。由于我们的年龄数据包含缺失值，mean(teens$age) 返回一个缺失值。我们可以通过添加额外的参数在计算均值之前删除缺失值来纠正这一点：

> mean(teens$age, na.rm = TRUE)
[1] 17.25243

这表明我们数据中的平均学生年龄大约为 17 岁。这只帮助我们达成了一部分目标；我们实际上需要每个毕业年份的平均年龄。你可能会想计算四次均值，但 R 的一个优点是通常有方法避免重复。在这种情况下，aggregate() 函数就是合适的工具。它计算数据子组的统计信息。在这里，它计算每个毕业年份的平均年龄，并在去除 NA 值后进行：

> aggregate(data = teens, age ~ gradyear, mean, na.rm = TRUE)
 gradyear      age
1     2006 18.65586
2     2007 17.70617
3     2008 16.76770
4     2009 15.81957

平均年龄每变化一次毕业年份大约差一年。这一点并不令人惊讶，但对于确认我们的数据合理性是一个有用的发现。

aggregate() 输出的是一个数据框。这对某些目的很有帮助，但需要额外的工作将其合并回原始数据。作为替代方案，我们可以使用 ave() 函数，它返回一个包含组均值的向量，这样结果的长度与原始向量相同：

> ave_age <- ave(teens$age, teens$gradyear, FUN =
 function(x) mean(x, na.rm = TRUE))

为了将这些均值填充到缺失值上，我们需要再调用一次 ifelse()，只有当原始年龄值为 NA 时才使用 ave_age 值：

> teens$age <- ifelse(is.na(teens$age), ave_age, teens$age)

summary() 结果显示，现在缺失值已被消除：

> summary(teens$age)
 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 13.03   16.28   17.24   17.24   18.21   20.00

数据准备好进行分析后，我们可以深入到这个项目的有趣部分。让我们看看我们的努力是否得到了回报。

第 3 步 – 在数据上训练模型

为了将青少年分为不同的营销群体，我们将使用 stats 包中的 k-means 实现，该包应该默认包含在你的 R 安装中。如果你碰巧没有这个包，可以像安装其他包一样安装它，并使用 library(stats) 命令加载它。尽管各种 R 包中有很多 k-means 函数，但 stats 包中的 kmeans() 函数被广泛使用，并提供了该算法的标准实现。

kmeans() 函数需要一个仅包含数值数据的数据框以及一个指定所需簇数的参数。如果你准备好了这两样东西，实际构建模型的过程就很简单。问题在于，选择合适的数据和簇的组合有时是一门艺术；这通常需要经过大量的反复试验。

我们将通过仅考虑代表青少年 SNS 个人资料中各种兴趣出现次数的 36 个特征来开始我们的聚类分析。为方便起见，我们创建一个仅包含这些特征的数据框：

> interests <- teens[5:40]

如果你记得第三章，懒学习 - 使用最近邻分类，在进行任何距离计算分析之前的常见做法是对特征进行归一化或 z-score 标准化，以确保每个特征使用相同的范围。通过这样做，你可以避免某些特征仅因为它们具有比其他特征更大的数值范围而主导结果的问题。

z-score 标准化过程会重新缩放特征，使其均值为零，标准差为一。这种转换改变了数据的解释方式，这在此处可能会有用。具体来说，如果某人在其个人资料中提到足球三次，在没有更多信息的情况下，我们无法判断这是否意味着他们比其他人更喜欢足球。另一方面，如果 z-score 为三，我们就知道他们提到足球的次数远远超过了平均水平的青少年。

要对interests数据框应用 z-score 标准化，我们可以使用lapply()配合scale()函数，如下所示：

> interests_z <- as.data.frame(lapply(interests, scale))

由于lapply()返回的是矩阵，因此必须使用as.data.frame()函数将其转换回数据框形式。

我们最后的决定是决定使用多少个簇来对数据进行分段。如果我们使用太多的簇，可能会发现它们过于具体，无法发挥作用；相反，选择太少的簇可能导致分组不均。你应该敢于尝试不同的k值。如果你不喜欢结果，可以轻松尝试另一个值并重新开始。

提示

如果你对分析群体有所了解，选择簇的数量会更加容易。对自然分组的真实数量有直觉可以帮助你节省一些试错的时间。

为了帮助我们预测数据中的簇的数量，我想引用我最喜欢的电影之一，《早餐俱乐部》，这是一部 1985 年上映的成长喜剧片，由约翰·休斯执导。电影中的青少年角色按照五种刻板印象进行分类：学霸、运动员、怪胎、公主和罪犯。考虑到这些身份在流行的青少年小说中常常出现，五个似乎是k的一个合理起点。

要使用 k-means 算法将青少年的兴趣数据分成五个簇，我们在interests数据框上使用kmeans()函数。由于 k-means 算法使用随机起始点，因此使用set.seed()函数来确保结果与以下示例中的输出一致。如果你记得前几章的内容，这个命令初始化了 R 的随机数生成器，设置为特定的序列。如果没有这个命令，每次运行 k-means 算法时结果都会不同：

> set.seed(2345)
> teen_clusters <- kmeans(interests_z, 5)

k-means 聚类过程的结果是一个名为teen_clusters的列表，存储了五个聚类的各项属性。让我们深入了解一下，看看算法是如何将青少年的兴趣数据划分的。

提示

如果你发现你的结果与这里展示的不同，请确保在运行kmeans()函数之前，立即执行set.seed(2345)命令。

步骤 4 – 评估模型性能

评估聚类结果可能是有一定主观性的。最终，模型的成功或失败取决于聚类是否能为其预期的目的提供帮助。由于本次分析的目标是识别具有相似兴趣的青少年群体，以便用于市场营销，我们将主要从定性角度来衡量成功。对于其他聚类应用，可能需要更多定量的成功衡量标准。

评估一组聚类有效性的最基本方法之一是检查每个组中示例的数量。如果这些组太大或太小，它们可能不会非常有用。要获取kmeans()聚类的大小，请使用teen_clusters$size组件，如下所示：

> teen_clusters$size
[1]   871   600  5981  1034 21514

在这里，我们看到了我们请求的五个聚类。最小的聚类有 600 个青少年（占 2%），而最大的聚类有 21,514 个（占 72%）。虽然最大和最小聚类之间人数差距较大，这有点令人担忧，但在没有更仔细检查这些组的情况下，我们无法知道这是否表示存在问题。也许，聚类大小的差异反映了某些实际情况，例如一大群有相似兴趣的青少年，或者它可能是由初始的 k-means 聚类中心引起的随机巧合。随着我们开始查看每个聚类的同质性，我们会了解更多。

提示

有时，k-means 可能会找到极小的聚类——有时甚至只有一个点。如果初始聚类中心恰好落在一个远离其他数据的离群点上，就可能发生这种情况。是否将这样的极小聚类视为一个真实的发现，代表一个极端案例的聚类，或者视为由随机机会引起的问题，并不总是很明确。如果遇到这个问题，可以考虑使用不同的随机种子重新运行 k-means 算法，看看这个小聚类是否对不同的起始点具有稳健性。

要更深入地了解这些聚类，我们可以通过teen_clusters$centers组件检查聚类中心的坐标，以下是前四个兴趣的情况：

> teen_clusters$centers
 basketball   football      soccer    softball
1  0.16001227  0.2364174  0.10385512  0.07232021
2 -0.09195886  0.0652625 -0.09932124 -0.01739428
3  0.52755083  0.4873480  0.29778605  0.37178877
4  0.34081039  0.3593965  0.12722250  0.16384661
5 -0.16695523 -0.1641499 -0.09033520 -0.11367669

输出的行（标记为1到5）表示五个聚类，而每行中的数字表示该聚类在每列顶部列出的兴趣项的平均值。由于这些值是经过 z 分数标准化的，正值表示该兴趣项高于所有青少年总体均值，负值表示低于总体均值。例如，第三行在篮球这一列中的值最高，这意味着聚类3在所有聚类中对篮球的平均兴趣最高。

通过检查聚类在每个兴趣类别上是否高于或低于平均水平，我们可以开始注意到区分聚类的模式。实际上，这涉及打印聚类中心并搜索其中的模式或极端值，类似于一个数字版的单词搜索谜题。以下截图显示了五个聚类在 36 个青少年兴趣中的 19 个兴趣的突出模式：

根据这一子集的兴趣数据，我们已经能够推断出一些聚类的特征。聚类 3在所有体育项目上的兴趣水平都显著高于平均水平。这表明这可能是早餐俱乐部刻板印象中的运动员群体。聚类 1包括最多提到“啦啦队”的内容、词汇“热”和高于平均水平的足球兴趣。这些是所谓的公主吗？

通过继续以这种方式检查聚类，我们可以构建一张表，列出每个群体的主要兴趣。在以下表格中，显示了每个聚类与其他聚类最具区别性的特征，以及最能准确描述该群体特点的早餐俱乐部身份。

有趣的是，聚类 5的特点在于其成员在每项活动中的兴趣水平都低于平均值，且它是成员数量最多的单一群体。一种可能的解释是，这些用户创建了一个网站个人资料，但从未发布任何兴趣。

提示

在分享分段分析结果时，通常应用具有信息量的标签可以帮助简化并捕捉群体的本质，如此处应用的早餐俱乐部类型学。加入此类标签的风险在于，它们可能通过刻板印象掩盖群体的细微差别。由于此类标签可能会偏向我们的思维，如果将标签视为全部真理，可能会错过一些重要的模式。

根据表格，营销主管将能清晰地描绘出五种类型的青少年社交网站访客。基于这些档案，主管可以向相关产品的企业出售有针对性的广告展示。在接下来的部分中，我们将看到如何将群体标签应用回原始人群以用于此类用途。

第五步 – 改进模型性能

由于聚类创造了新信息，聚类算法的表现至少在某种程度上取决于群体本身的质量以及如何利用这些信息。在前面的部分中，我们已经展示了这五个群体为青少年的兴趣提供了有用且新颖的见解。以此衡量，算法的表现相当不错。因此，我们现在可以将精力集中在将这些见解转化为行动上。

我们将首先将群体应用回完整数据集。由kmeans()函数创建的teen_clusters对象包含一个名为cluster的组件，该组件包含样本中 30,000 个个体的群体分配信息。我们可以通过以下命令将其作为一列添加到teens数据框中：

> teens$cluster <- teen_clusters$cluster

鉴于这一新数据，我们可以开始检查群体分配与个人特征之间的关系。例如，以下是 SNS 数据中前五个青少年的个人信息：

> teens[1:5, c("cluster", "gender", "age", "friends")]
 cluster gender    age friends
1       5      M 18.982       7
2       3      F 18.801       0
3       5      M 18.335      69
4       5      F 18.875       0
5       4   <NA> 18.995      10

使用aggregate()函数，我们还可以查看各群体的不同人口特征。群体间的平均年龄变化不大，这并不令人惊讶，因为这些青少年的身份通常在上高中之前就已确定。情况如下所示：

> aggregate(data = teens, age ~ cluster, mean)
 cluster      age
1       1 16.86497
2       2 17.39037
3       3 17.07656
4       4 17.11957
5       5 17.29849

另一方面，不同群体中女性所占比例存在一些显著差异。这是一个非常有趣的发现，因为我们并未使用性别数据来创建群体，但群体仍然能预测性别：

> aggregate(data = teens, female ~ cluster, mean)
 cluster    female
1       1 0.8381171
2       2 0.7250000
3       3 0.8378198
4       4 0.8027079
5       5 0.6994515

回想一下，SNS 用户中大约 74%是女性。群体 1，即所谓的公主群体，女性比例接近 84%，而群体 2和群体 5的女性比例仅约为 70%。这些差异表明青少年男孩和女孩在社交网络页面上讨论的兴趣有所不同。

鉴于我们在预测性别上的成功，您可能也会怀疑群体能否预测用户拥有的朋友数量。这个假设似乎得到了数据的支持，具体如下：

> aggregate(data = teens, friends ~ cluster, mean)
 cluster  friends
1       1 41.43054
2       2 32.57333
3       3 37.16185
4       4 30.50290
5       5 27.70052

平均而言，公主拥有最多的朋友（41.4），其次是运动员（37.2）和学霸（32.6）。最低的是罪犯（30.5）和精神病患者（27.7）。与性别一样，青少年朋友数量与其预测聚类之间的联系非常显著，尽管我们并未将朋友数据作为聚类算法的输入。另外有趣的是，朋友数量似乎与每个群体在高中受欢迎程度的刻板印象相关；那些在刻板印象中受欢迎的群体往往拥有更多的朋友。

群体成员、性别和朋友数量之间的关联表明，聚类可以是行为的有用预测指标。通过这种方式验证它们的预测能力，可以使聚类在向营销团队推销时更具吸引力，从而最终提升算法的表现。

总结

我们的发现支持了那句流行的格言：“物以类聚，人以群分。”通过使用机器学习方法将青少年与有相似兴趣的人进行聚类，我们能够开发出一种青少年身份的分类法，该分类法能够预测个人特征，如性别和朋友数量。这些相同的方法也可以应用于其他情境，并取得类似的结果。

本章仅介绍了聚类的基本概念。作为一种非常成熟的机器学习方法，k-means 算法有许多变体，还有许多其他聚类算法，它们为任务带来了独特的偏差和启发式方法。基于本章的基础，你将能够理解并应用其他聚类方法来解决新的问题。

在下一章，我们将开始探讨衡量学习算法成功的方法，这些方法适用于许多机器学习任务。虽然我们的过程一直在评估学习的成功，但为了获得最高的性能水平，能够在最严格的术语下定义和衡量它是至关重要的。

第十章：评估模型性能

当只有富人能负担得起教育时，测试和考试并没有评估学生的潜力。相反，教师是根据家长的要求来评判的，家长们希望知道他们的孩子是否学到了足够的知识，以证明教员的薪水。显然，随着时间的推移，这种情况发生了变化。现在，这些评估被用来区分高成就和低成就的学生，并将他们筛选到职业和其他机会中。

鉴于这个过程的重要性，投入了大量精力来开发准确的学生评估。公平的评估有大量的问题，覆盖广泛的主题，奖励真实的知识而不是运气猜测。它们还要求学生思考他们之前从未遇到过的问题。因此，正确的回答表明学生能够更广泛地概括他们的知识。

评估机器学习算法的过程与评估学生的过程非常相似。由于算法有不同的优缺点，测试应该能够区分不同的学习者。预测学习者在未来数据上的表现也同样重要。

本章提供了评估机器学习者所需的信息，例如：

• 为什么预测准确率不足以衡量性能，以及你可以使用的其他性能度量

• 确保性能度量合理反映模型预测或预测未见过案例的能力的方法

• 如何使用 R 语言将这些更有用的度量和方法应用到前几章中讨论的预测模型上

就像学习某个主题的最佳方式是尝试将其教授给别人一样，教学和评估机器学习者的过程将为你提供更多关于迄今为止所学方法的洞察。

衡量分类的性能

在前几章中，我们通过将正确预测的比例除以预测的总数来衡量分类器的准确性。这表示学习者在多少情况下是正确的或错误的。例如，假设在 100,000 名新生儿中，有 99,990 名婴儿的基因缺陷是否携带被分类器正确预测。这样的话，准确率将是 99.99%，错误率仅为 0.01%。

乍一看，这似乎是一个极其准确的分类器。然而，在将你孩子的生命交给该测试之前，最好先收集更多的信息。如果这种基因缺陷仅在每 10 万个婴儿中有 10 个发现，怎么办？无论在什么情况下，始终预测没有缺陷的测试对 99.99%的所有案例都是正确的，但对 100%最重要的案例却是错误的。换句话说，尽管预测非常准确，但这个分类器对于防止可治疗的出生缺陷并没有多大帮助。

提示

这是类别不平衡问题的一个后果，指的是数据中大多数记录属于同一类别所带来的问题。

尽管有许多方法可以衡量分类器的性能，但最好的衡量标准总是能够捕捉分类器在其预期目标上是否成功的标准。定义性能度量时，关键是要以效用为导向，而非单纯的准确率。为此，我们将开始探索从混淆矩阵中衍生出的各种替代性性能度量方法。然而，在我们开始之前，我们需要考虑如何准备分类器进行评估。

在 R 中处理分类预测数据

评估分类模型的目标是更好地理解其性能如何推断到未来的案例。由于在实际环境中测试一个尚未验证的模型通常是不可行的，我们通常通过要求模型对一个包含类似未来任务的案例的数据集进行分类，从而模拟未来的情况。通过观察学习者对这一检验的回应，我们可以了解其优点和缺点。

尽管我们在前面的章节中评估了分类器，但值得反思一下我们所拥有的数据类型：

• 实际类别值

• 预测类别值

• 预测的估计概率

实际值和预测的类别值可能是显而易见的，但它们是评估的关键。就像老师用答案解析来评估学生的答案一样，我们需要知道机器学习者预测的正确答案。目标是保持两个数据向量：一个存储正确或实际的类别值，另一个存储预测的类别值。这两个向量必须包含相同数量的值，并按相同的顺序排列。预测值和实际值可以存储为独立的 R 向量，或者在一个 R 数据框中作为列存储。

获取这些数据很容易。实际的类别值直接来自测试数据集中的目标特征。预测类别值是通过基于训练数据构建的分类器来获取的，并应用于测试数据。对于大多数机器学习包来说，这通常涉及对模型对象和测试数据框应用predict()函数，例如：predicted_outcome <- predict(model, test_data)。

到目前为止，我们仅仅使用这两个数据向量来检查分类预测。然而，大多数模型可以提供另一个有用的信息。即使分类器对每个样本做出单一的预测，它对于某些决策的信心可能会高于其他决策。例如，分类器可能有 99% 的把握认为包含“免费”和“铃声”字样的短信是垃圾短信，但对含有“今晚”字样的短信只有 51% 的把握是垃圾短信。在这两种情况下，分类器都会将消息归类为垃圾短信，但对其中一个决策的信心远高于另一个。

研究这些内部预测概率可以提供有用的数据来评估模型的表现。如果两个模型犯了相同数量的错误，但其中一个能够更准确地评估其不确定性，那么这个模型更为智能。理想情况下，应该找到一个在做出正确预测时非常自信，而在面对不确定时则保持谨慎的学习者。信心与谨慎之间的平衡是模型评估的关键部分。

不幸的是，获取内部预测概率可能有些棘手，因为不同的分类器获取预测概率的方法不同。通常，对于大多数分类器，predict()函数用于指定所需的预测类型。要获取单一预测类别（如垃圾邮件或正常邮件），通常需要将type = "class"参数设置为该值。要获取预测概率，type参数应根据所使用的分类器设置为"prob"、"posterior"、"raw"或"probability"之一。

提示

本书中介绍的几乎所有分类器都会提供预测概率。在每个模型的语法框中都会包含type参数。

例如，要输出在第五章中构建的 C5.0 分类器的预测概率，可以使用predict()函数，并设置type = "prob"，如下所示：

> predicted_prob <- predict(credit_model, credit_test, type = "prob")

为了进一步说明评估学习算法的过程，让我们更详细地看看在第四章中开发的 SMS 垃圾邮件分类模型的表现，概率学习 – 使用朴素贝叶斯分类。要输出朴素贝叶斯的预测概率，可以使用predict()函数，并设置type = "raw"，如下所示：

> sms_test_prob <- predict(sms_classifier, sms_test, type = "raw")

在大多数情况下，predict()函数会为每个结果类别返回一个概率。例如，在像 SMS 分类器这样的二分类模型中，预测的概率可能是一个矩阵或数据框，如下所示：

> head(sms_test_prob)
 ham         spam
[1,] 9.999995e-01 4.565938e-07
[2,] 9.999995e-01 4.540489e-07
[3,] 9.998418e-01 1.582360e-04
[4,] 9.999578e-01 4.223125e-05
[5,] 4.816137e-10 1.000000e+00
[6,] 9.997970e-01 2.030033e-04

输出中的每一行显示了分类器对垃圾邮件和正常邮件的预测概率，这两个概率的总和始终为 1，因为这是唯一的两个可能结果。在构建评估数据集时，确保你使用的是与所关注类别级别相符的正确概率非常重要。为了避免混淆，在二分类情况下，甚至可以考虑去掉其中一个类别的预测向量。

为了方便评估过程，可以构建一个数据框，包含预测的类别值、实际类别值，以及感兴趣的估计概率。

提示

构建评估数据集的步骤为了简洁起见已被省略，但它们包含在本章的代码中，可以在 Packt Publishing 网站上找到。要跟随本示例操作，请下载sms_results.csv文件，并使用sms_results <- read.csv("sms_results.csv")命令将其加载到数据框中。

sms_results数据框非常简单，它包含四个向量，包含 1,390 个值。一个向量包含表示实际短信类型（spam或ham）的值，一个向量表示朴素贝叶斯模型的预测类型，第三个和第四个向量分别表示消息是spam或ham的概率：

> head(sms_results)
 actual_type predict_type prob_spam prob_ham
1         ham          ham   0.00000  1.00000
2         ham          ham   0.00000  1.00000
3         ham          ham   0.00016  0.99984
4         ham          ham   0.00004  0.99996
5        spam         spam   1.00000  0.00000
6         ham          ham   0.00020  0.99980

对于这六个测试案例，预测值与实际的短信类型一致；模型正确地预测了它们的状态。此外，预测的概率表明模型对这些预测极其自信，因为它们的值都接近零或一。

当预测值和实际值远离零和一时，会发生什么？使用subset()函数，我们可以找出一些这样的记录。以下输出显示了模型预测的spam概率介于 40 到 60 百分比之间的测试案例：

> head(subset(sms_results, prob_spam > 0.40 & prob_spam < 0.60))
 actual_type predict_type prob_spam prob_ham
377         spam          ham   0.47536  0.52464
717          ham         spam   0.56188  0.43812
1311         ham         spam   0.57917  0.42083

根据模型自己的说明，这些是正确预测几乎等同于掷硬币的情况。然而，所有三个预测都是错误的——这是一个不幸的结果。让我们再看看几个模型错误的案例：

> head(subset(sms_results, actual_type != predict_type))
 actual_type predict_type prob_spam prob_ham
53         spam          ham   0.00071  0.99929
59         spam          ham   0.00156  0.99844
73         spam          ham   0.01708  0.98292
76         spam          ham   0.00851  0.99149
184        spam          ham   0.01243  0.98757
332        spam          ham   0.00003  0.99997

这些案例说明了一个重要的事实：模型可以非常有信心，但也可能极其错误。所有这六个测试案例都是spam，而分类器认为它们被判定为ham的概率不低于 98%。

尽管存在这些错误，这个模型是否仍然有用呢？我们可以通过对评估数据应用各种错误度量来回答这个问题。事实上，许多这样的度量是基于我们在前几章中已广泛使用的工具。

更深入地看混淆矩阵

混淆 矩阵是一个表格，用于根据预测值是否与实际值匹配来对预测进行分类。表格的一个维度表示预测值的可能类别，另一个维度表示实际值的类别。尽管我们至今只见过 2 x 2 的混淆矩阵，但也可以为预测任何类别值的模型创建矩阵。下图展示了熟悉的二分类模型的混淆矩阵，以及三类模型的 3 x 3 混淆矩阵。

当预测值与实际值相同，说明是正确分类。正确的预测会出现在混淆矩阵的对角线上（用O表示）。对角线外的矩阵单元格（用X表示）表示预测值与实际值不同的情况，这些是错误预测。分类模型的性能度量是基于这些表格中对角线上的预测数和对角线外的预测数：

最常见的性能衡量指标考虑的是模型区分一个类别与所有其他类别的能力。关注类别被称为正类，而所有其他类别被称为负类。

提示

使用“正类”和“负类”这些术语并不意味着任何价值判断（即好与坏），也不一定表示结果是存在或不存在（例如，出生缺陷与否）。正类的选择甚至可以是任意的，比如在模型预测“晴天与雨天”或“狗与猫”等类别的情况下。

正类和负类预测之间的关系可以通过一个 2 x 2 的混淆矩阵来表示，矩阵记录了预测是否属于以下四个类别之一：

• 真正（TP）：正确地分类为关注类别

• 真负（TN）：正确地分类为非关注类别

• 假正（FP）：错误地分类为关注类别

• 假负（FN）：错误地分类为非关注类别

对于垃圾邮件分类器，正类是 spam，因为这是我们希望检测的结果。我们可以将混淆矩阵想象为以下示意图所示：

以这种方式呈现的混淆矩阵是许多重要模型性能指标的基础。在接下来的部分，我们将使用该矩阵更好地理解准确率的含义。

使用混淆矩阵衡量性能

使用 2 x 2 混淆矩阵，我们可以形式化定义预测的准确率（有时称为成功率）为：

在此公式中，TP、TN、FP 和 FN 指的是模型预测落入这些类别的次数。因此，准确率是一个比例，表示真正例和真负例的数量，除以总预测数量。

错误率或错误分类示例的比例定义为：

请注意，错误率可以通过 1 减去准确率来计算。直观上，这是有道理的；一个 95% 正确的模型，其错误率为 5%。

将分类器的预测结果汇总到混淆矩阵中，一种简便的方法是使用 R 的 table() 函数。创建 SMS 数据混淆矩阵的命令如下所示。该表格中的计数值可以用于计算准确率和其他统计数据：

> table(sms_results$actual_type, sms_results$predict_type)

 ham spam
 ham  1203    4
 spam   31  152

如果你想创建一个更具信息性的混淆矩阵输出，gmodels包中的CrossTable()函数提供了一个可定制的解决方案。如果你还记得，我们在第二章，管理与理解数据中首次使用了这个函数。如果你当时没有安装这个包，你需要使用install.packages("gmodels")命令进行安装。

默认情况下，CrossTable()的输出包括每个单元格中的比例，表示该单元格的计数占表格的行、列或总体总计的百分比。输出结果还包括行和列的总计。如下面的代码所示，语法与table()函数类似：

> library(gmodels)
> CrossTable(sms_results$actual_type, sms_results$predict_type)

结果是一个包含大量附加细节的混淆矩阵：

我们在前几章中已经使用了CrossTable()，所以现在你应该对输出结果比较熟悉。如果你忘记了如何解读输出结果，只需参考关键部分（标记为Cell Contents），它提供了表格单元格中每个数字的定义。

我们可以使用混淆矩阵来获得准确率和误差率。由于准确率是*(TP + TN) / (TP + TN + FP + FN)*，我们可以使用以下命令来计算它：

> (152 + 1203) / (152 + 1203 + 4 + 31)
[1] 0.9748201

我们还可以计算误差率*(FP + FN) / (TP + TN + FP + FN)*，方法如下：

> (4 + 31) / (152 + 1203 + 4 + 31)
[1] 0.02517986

这与准确度的补集相同：

> 1 - 0.9748201
[1] 0.0251799

尽管这些计算看起来很简单，但重要的是要练习思考混淆矩阵的各个组成部分是如何相互关联的。在接下来的章节中，你将看到如何将这些组件以不同方式组合，从而创建各种附加的性能度量。

除了准确度——其他的性能度量

无数的性能度量已经为特定目的在诸如医学、信息检索、市场营销和信号检测理论等领域开发并使用。涵盖所有这些度量将填满数百页，因此在这里进行全面描述是不可行的。相反，我们将只考虑机器学习文献中最常用和最常引用的一些度量。

Max Kuhn 的分类与回归训练包caret包括计算许多此类性能度量的函数。该包提供了大量的工具，用于准备、训练、评估和可视化机器学习模型和数据。除了在这里的使用外，我们还将在第十一章，提升模型性能中广泛使用caret。在继续之前，你需要使用install.packages("caret")命令安装该包。

提示

关于caret的更多信息，请参考：Kuhn M. 使用 caret 包在 R 中构建预测模型。统计学期刊软件。2008 年；28。

caret 包提供了另一个函数来创建混淆矩阵。如下命令所示，其语法与 table() 类似，但有一个小的差异。因为 caret 提供了考虑到分类正类能力的模型性能度量，所以应指定 positive 参数。在本例中，由于 SMS 分类器旨在检测 spam，我们将设置 positive = "spam"，如下所示：

> library(caret)
> confusionMatrix(sms_results$predict_type,
 sms_results$actual_type, positive = "spam")

这将产生如下输出：

输出顶部是一个混淆矩阵，类似于 table() 函数生成的矩阵，但它被转置了。输出还包括一组性能度量。其中一些，如准确性，是我们熟悉的，而许多其他度量则是新的。让我们看看几个最重要的指标。

Kappa 统计量

Kappa 统计量（在之前的输出中标记为Kappa）通过考虑仅凭随机猜测就能做出正确预测的可能性来调整准确性。这对于具有严重类别不平衡的数据集尤为重要，因为分类器只需始终猜测最频繁的类别就能获得高准确率。Kappa 统计量只有在分类器的正确率超过这种简单策略时，才会给予奖励。

Kappa 值的范围从 0 到最大值 1，表示模型预测与真实值之间的完美协议。值小于 1 表示协议不完全。根据模型的使用方式，Kappa 统计量的解释可能有所不同。以下是常见的解释：

• 差的协议 = 小于 0.20

• 公平协议 = 0.20 到 0.40

• 中等协议 = 0.40 到 0.60

• 良好的协议 = 0.60 到 0.80

• 非常好的协议 = 0.80 到 1.00

需要注意的是，这些类别是主观的。虽然“良好的协议”可能足以预测某人最喜欢的冰淇淋口味，但如果目标是识别出生缺陷，单凭“非常好的协议”可能不足够。

提示

有关前述量表的更多信息，请参阅：Landis JR, Koch GG. The measurement of observer agreement for categorical data. Biometrics. 1997; 33:159-174.

以下是计算 Kappa 统计量的公式。在这个公式中，Pr(a) 指的是实际协议的比例，而 Pr(e) 指的是在假设随机选择的情况下，分类器与真实值之间的期望协议：

提示

定义 Kappa 统计量的方法不止一种。这里描述的最常见方法使用 Cohen 的 Kappa 系数，该方法在论文中有所阐述：Cohen J. A coefficient of agreement for nominal scales. Education and Psychological Measurement. 1960; 20:37-46.

这些比例可以通过混淆矩阵轻松获得，一旦你知道该从哪里查找。让我们考虑使用CrossTable()函数创建的 SMS 分类模型的混淆矩阵，为了方便起见，这里重复显示：

请记住，每个单元格底部的值表示所有实例中落入该单元格的比例。因此，计算观察到的一致性*Pr(a)时，我们只需将预测类型与实际短信类型一致的所有实例的比例相加。这样，我们可以计算Pr(a)*如下：

> pr_a <- 0.865 + 0.109
> pr_a
[1] 0.974

对于这个分类器，观察值和实际值有 97.4%的时间是一致的——你会注意到这与准确度是相同的。kappa 统计量根据预期一致性*Pr(e)*调整了准确度，*Pr(e)*是指在假设两者都是根据观察到的比例随机选择的前提下，仅凭运气，预测值和实际值匹配的概率。

为了找到这些观察到的比例，我们可以使用我们在第四章中学到的概率规则，概率学习 – 使用朴素贝叶斯分类。假设两个事件是独立的（意味着一个事件不会影响另一个事件），概率规则指出，两者同时发生的概率等于每个事件发生概率的乘积。例如，我们知道两者都选择正常邮件的概率是：

Pr(实际类型是正常邮件) * Pr(预测类型是正常邮件)

两者都选择垃圾邮件的概率是：

Pr(实际类型是垃圾邮件) * Pr(预测类型是垃圾邮件)

预测类型或实际类型是垃圾邮件（spam）或正常邮件（ham）的概率可以从行或列的总计中获得。例如，Pr(实际类型是正常邮件) = 0.868 和 Pr(预测类型是正常邮件) = 0.888。

Pr(e)是通过计算预测值和实际值因运气而一致的概率之和来计算的，无论消息是垃圾邮件还是正常邮件。回想一下，对于互斥事件（不能同时发生的事件），发生任意一个的概率等于它们各自概率的总和。因此，为了获得最终的Pr(e)，我们只需将两个乘积相加，如以下命令所示：

> pr_e <- 0.868 * 0.888 + 0.132 * 0.112
> pr_e
[1] 0.785568

由于*Pr(e)*是 0.786，单纯依靠运气，我们预期观察值与实际值大约有 78.6%的时间是一致的。

这意味着我们现在拥有了完成 kappa 公式所需的所有信息。将*Pr(a)和Pr(e)*值代入 kappa 公式中，我们得到：

> k <- (pr_a - pr_e) / (1 - pr_e)
> k
[1] 0.8787494

kappa 大约是 0.88，这与caret的confusionMatrix()输出一致（小差异是由于四舍五入）。根据建议的解释，我们注意到分类器的预测与实际值之间有很好的一致性。

有几个 R 函数可以自动计算 kappa。Visualizing Categorical Data (vcd) 包中的 Kappa() 函数（请注意大写的 'K'）使用预测值和实际值的混淆矩阵。在安装该包后（使用命令 install.packages("vcd")），可以使用以下命令获取 kappa 值：

> library(vcd)
> Kappa(table(sms_results$actual_type, sms_results$predict_type))
 value        ASE
Unweighted 0.8825203 0.01949315
Weighted   0.8825203 0.01949315

我们关心的是不带权的 kappa。值 0.88 与我们预期的相符。

提示

加权 kappa 用于存在不同程度一致性的情况。例如，使用“冷、凉、温暖、热”这样的尺度时，温暖与热的值更为接近，而与冷的值差异较大。在二分类事件中，例如垃圾邮件和正常邮件，带权和不带权的 kappa 统计量将是相同的。

Inter-Rater Reliability (irr) 包中的 kappa2() 函数可以用来计算存储在数据框中的预测值和实际值向量的 kappa。安装该包（使用 install.packages("irr") 命令）后，可以使用以下命令获取 kappa 值：

> kappa2(sms_results[1:2])
 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 1390 
 Raters = 2 
 Kappa = 0.883 

 z = 33 
 p-value = 0

Kappa() 和 kappa2() 函数报告相同的 kappa 统计量，因此你可以选择更熟悉的函数。

提示

小心不要使用内置的 kappa() 函数。它与之前报告的 kappa 统计量完全无关！

灵敏度和特异性

寻找一个有用的分类器通常涉及在过于保守和过于激进的预测之间做出平衡。例如，一个电子邮件过滤器可以通过激进地过滤几乎所有的正常邮件来确保删除每一封垃圾邮件。另一方面，为了确保不误过滤正常邮件，可能需要允许不可接受的垃圾邮件通过过滤器。一对性能度量捕捉了这种权衡：灵敏度和特异性。

模型的灵敏度（也叫做真正率）衡量的是正例中被正确分类的比例。因此，如下公式所示，它是通过将真正例数除以所有正例的总数（包括正确分类的真正例和错误分类的假负例）来计算的：

模型的特异性（也叫做真负率）衡量的是负例中被正确分类的比例。与灵敏度类似，这个值是通过将真负例数除以所有负例的总数（包括真负例和假正例）来计算的：

给定短信分类器的混淆矩阵，我们可以轻松手动计算这些度量。假设垃圾邮件为正类，我们可以确认 confusionMatrix() 输出中的数字是正确的。例如，灵敏度的计算公式如下：

> sens <- 152 / (152 + 31)
> sens
[1] 0.8306011

同样，对于特异性，我们可以计算：

> spec <- 1203 / (1203 + 4)
> spec
[1] 0.996686

caret 包提供了从预测值和实际值的向量直接计算灵敏度和特异性（sensitivity and specificity）的函数。请小心地指定 positive 或 negative 参数，如下所示：

> library(caret)
> sensitivity(sms_results$predict_type, sms_results$actual_type,
 positive = "spam")
[1] 0.8306011

> specificity(sms_results$predict_type, sms_results$actual_type,
 negative = "ham")
[1] 0.996686

灵敏度和特异性范围从 0 到 1，接近 1 的值更为理想。当然，找到两者之间的适当平衡是很重要的——这一任务通常是特定于上下文的。

例如，在这种情况下，0.831 的灵敏度意味着 83.1%的垃圾短信被正确分类。类似地，0.997 的特异性意味着 99.7%的非垃圾短信被正确分类；或者，0.3%的有效短信被误判为垃圾短信。拒绝 0.3%的有效短信可能是不可接受的，或者考虑到垃圾短信减少，这可能是一个合理的权衡。

灵敏度和特异性为思考这种权衡提供了工具。通常，通过对模型进行修改并测试不同的模型，直到找到一个满足期望的灵敏度和特异性阈值的模型。可视化工具（如本章稍后讨论的内容）也有助于理解灵敏度和特异性之间的权衡。

精准度与召回率

与灵敏度和特异性密切相关的，还有两个与分类妥协相关的性能度量：精准度和召回率。它们主要用于信息检索领域，旨在提供模型结果的相关性和有趣性指示，或者判断预测是否被无意义的噪声稀释。

精准度（也称为正预测值）定义为真正的正例所占比例；换句话说，当模型预测为正类时，它有多大的准确性？一个精准的模型只会在非常可能为正的情况下预测为正类，因此它非常值得信赖。

想象一下，如果模型非常不精准，会发生什么。随着时间的推移，结果的可信度会降低。在信息检索的背景下，这类似于一个像 Google 这样的搜索引擎返回无关的结果。最终，用户会转向像 Bing 这样的竞争对手。在短信垃圾过滤器的案例中，高精准度意味着模型能够精确地识别垃圾短信，同时忽略非垃圾短信。

另一方面，召回率是衡量结果完整性的一个指标。如以下公式所示，召回率定义为真阳性数占所有阳性数的比例。你可能已经认识到这与灵敏度相同。然而，在这种情况下，解释略有不同。一个具有高召回率的模型捕捉到了大量的正例，这意味着它具有广泛的覆盖范围。例如，一个具有高召回率的搜索引擎会返回大量与搜索查询相关的文档。同样，短信垃圾邮件过滤器具有高召回率时，意味着大多数垃圾短信都能被正确识别。

我们可以从混淆矩阵中计算精确度和召回率。再次假设spam是正类，精确度为：

> prec <- 152 / (152 + 4)
> prec
[1] 0.974359

召回率为：

> rec <- 152 / (152 + 31)
> rec
[1] 0.8306011

caret包可以用来从预测和实际类别的向量中计算这些度量之一。精确度使用posPredValue()函数：

> library(caret)
> posPredValue(sms_results$predict_type, sms_results$actual_type,
 positive = "spam")
[1] 0.974359

召回率使用我们之前使用的sensitivity()函数：

> sensitivity(sms_results$predict_type, sms_results$actual_type, 
 positive = "spam")
[1] 0.8306011

类似于灵敏度和特异性之间的固有权衡，对于大多数实际问题而言，很难构建一个同时具有高精确度和高召回率的模型。如果你只针对简单的、容易分类的样本，保持精确度就变得容易。同样，如果一个模型通过使用一个非常宽泛的筛选标准来捕捉尽可能多的正例，那么它的召回率会很高。在这种情况下，模型会过于激进地识别正样本。相比之下，同时具备高精度和高召回率非常具有挑战性。因此，测试多种模型以找到精度和召回率的最佳组合，以满足项目需求是至关重要的。

F-measure

一种结合了精确度和召回率的模型性能度量方法称为F-measure（有时也叫F[1]分数或F-score）。F-measure 通过调和平均数将精确度和召回率结合起来，调和平均数是一种用于变化率的平均值类型。由于精确度和召回率都表示为介于零和一之间的比例，可以解释为比率，因此使用调和平均数而非常见的算术平均数。以下是 F-measure 的公式：

要计算 F-measure，请使用之前计算的精确度和召回率值：

> f <- (2 * prec * rec) / (prec + rec)
> f
[1] 0.8967552

这个计算结果与使用混淆矩阵中的计数值完全相同：

> f <- (2 * 152) / (2 * 152 + 4 + 31)
> f
[1] 0.8967552

由于 F 度量在一个数字中描述了模型的性能，它提供了一种方便的方式来并排比较多个模型。然而，这假设了精确度和召回率应该赋予相同的权重，这一假设并不总是有效。可以使用不同的权重来计算 F 分数，但选择权重可能在最好的情况下比较棘手，最坏的情况下则显得任意。更好的做法是将像 F 分数这样的度量与更全面考虑模型优缺点的方法结合使用，如下一节中描述的方法。

可视化性能权衡

可视化有助于更详细地理解机器学习算法的性能。当统计数据如敏感性和特异性，或精确度和召回率试图将模型性能简化为一个数字时，可视化则描绘了学习者在各种条件下的表现。

由于学习算法有不同的偏差，两个模型可能在准确率相似的情况下，在实现准确率的方式上存在巨大的差异。有些模型可能在一些预测上遇到困难，而其他模型则轻松完成这些预测，同时对于其他模型无法正确预测的情况表现得游刃有余。可视化提供了一种理解这些权衡的方法，通过将多个学习者并排比较在一个图表中。

ROCR包提供了一套易于使用的函数，用于可视化分类模型的性能。它包括用于计算常见性能度量和可视化的大量函数。ROCR官网 rocr.bioinf.mpi-sb.mpg.de/ 列出了完整的功能集以及多个可视化功能示例。继续之前，请使用install.packages("ROCR")命令安装该包。

提示

有关 ROCR 开发的更多信息，请参见：Sing T, Sander O, Beerenwinkel N, Lengauer T. ROCR：在 R 中可视化分类器性能。生物信息学。2005；21：3940-3941。

要使用ROCR创建可视化，需要两个数据向量。第一个必须包含预测的类别值，第二个必须包含正类的估计概率。这些数据用于创建预测对象，然后可以通过ROCR的绘图功能进行检查。

SMS 分类器的预测对象需要分类器的垃圾邮件概率估计值和实际类别标签。这些数据通过prediction()函数结合在以下几行中：

> library(ROCR)
> pred <- prediction(predictions = sms_results$prob_spam,
 labels = sms_results$actual_type)

接下来，performance()函数将允许我们从刚刚创建的prediction对象中计算性能度量，然后可以使用 R 的plot()函数进行可视化。通过这三步，可以创建多种有用的可视化图。

ROC 曲线

接收者操作特征（ROC）曲线通常用于检验在避免假阳性的同时，检测真实阳性的权衡。正如你可能从名称中猜到的，ROC 曲线最初由通信领域的工程师开发。在二战时期，雷达和无线电操作员使用 ROC 曲线来衡量接收器区分真实信号和假警报的能力。今天，这一技术在可视化机器学习模型的有效性时依然非常有用。

典型 ROC 图的特征如下面的图所示。曲线定义在一个图上，纵轴表示真实阳性比例，横轴表示假阳性比例。由于这些值分别等同于灵敏度和（1 – 特异性），因此该图也被称为灵敏度/特异性图：

组成 ROC 曲线的点表示在不同假阳性阈值下的真实阳性率。为了创建这些曲线，分类器的预测结果按模型对正类的估计概率排序，最大的值排在前面。从原点开始，每个预测对真实阳性率和假阳性率的影响将导致曲线向上（对于正确预测）或向右（对于错误预测）延伸。

为了说明这一概念，前面的图中对比了三种假设的分类器。首先，从图的左下角到右上角的对角线代表一个没有预测价值的分类器。这种分类器以相同的速度检测到真实阳性和假阳性，意味着分类器无法区分二者。这是其他分类器评判的基准线。接近此线的 ROC 曲线表示模型没有太大用处。完美分类器的曲线通过一个点，表示 100%的真实阳性率和 0%的假阳性率。它能够在错误分类任何负结果之前正确识别所有的正结果。大多数现实世界的分类器与测试分类器类似，它们的表现位于完美分类器和无用分类器之间。

曲线越接近完美分类器，越能更好地识别正值。这可以通过一个统计量来衡量，称为ROC 曲线下面积（简称AUC）。AUC 将 ROC 图作为一个二维方形，并测量 ROC 曲线下的总面积。AUC 的值范围从 0.5（表示分类器没有预测价值）到 1.0（表示完美分类器）。解释 AUC 分数的惯例使用一个类似于学术字母评分的系统：

• A：优秀 = 0.9 到 1.0

• B：优秀/良好 = 0.8 到 0.9

• C：可接受/一般 = 0.7 到 0.8

• D：差 = 0.6 到 0.7

• E：无区分能力 = 0.5 到 0.6

与大多数类似的量表一样，某些任务可能比其他任务更适合使用这些等级；这种分类是有一定主观性的。

提示

还值得注意的是，两个 ROC 曲线的形状可能截然不同，但 AUC 却相同。正因为如此，单独使用 AUC 可能会产生误导。最佳做法是结合 AUC 和 ROC 曲线的定性检查一起使用。

使用ROCR包创建 ROC 曲线涉及从我们之前计算的prediction对象中构建一个performance对象。由于 ROC 曲线绘制的是真实正例率与假正例率之间的关系，我们只需调用performance()函数并指定tpr和fpr这两个度量值，如下代码所示：

> perf <- performance(pred, measure = "tpr", x.measure = "fpr")

使用perf对象，我们可以通过 R 的plot()函数可视化 ROC 曲线。如以下代码所示，可以使用许多标准参数来调整可视化效果，例如main（添加标题）、col（改变线条颜色）和lwd（调整线条宽度）：

> plot(perf, main = "ROC curve for SMS spam filter",
 col = "blue", lwd = 3)

尽管plot()命令足以创建有效的 ROC 曲线，但添加参考线来指示一个没有预测价值的分类器的表现会更有帮助。

为了绘制这样的曲线，我们将使用abline()函数。这个函数可以用来指定一个斜截式方程，其中a是截距，b是斜率。由于我们需要一条通过原点的单位线，我们将截距设置为a=0，斜率设置为b=1，如下图所示。lwd参数调整线条的粗细，而lty参数调整线条的类型。例如，lty = 2表示虚线：

> abline(a = 0, b = 1, lwd = 2, lty = 2)

最终结果是带有虚线参考线的 ROC 图：

从定性上看，我们可以看到这条 ROC 曲线似乎占据了图表的左上角区域，这表明它比表示无用分类器的虚线更接近一个完美的分类器。为了定量验证这一点，我们可以使用 ROCR 包来计算 AUC。为此，我们首先需要创建另一个performance对象，这次指定measure = "auc"，如下代码所示：

> perf.auc <- performance(pred, measure = "auc")

由于perf.auc是一个 R 对象（具体来说是 S4 对象），我们需要使用特殊的符号来访问其中存储的值。S4 对象在被称为槽位的位置存储信息。可以使用str()函数查看一个对象的所有槽位：

> str(perf.auc)
Formal class 'performance' [package "ROCR"] with 6 slots
 ..@ x.name      : chr "None"
 ..@ y.name      : chr "Area under the ROC curve"
 ..@ alpha.name  : chr "none"
 ..@ x.values    : list()
 ..@ y.values    :List of 1
 .. ..$ : num 0.984
 ..@ alpha.values: list()

请注意，槽位前面有@符号。在访问存储在y.values槽位中的 AUC 值时，我们可以使用@符号和unlist()函数，后者将列表简化为数值向量：

> unlist(perf.auc@y.values)
[1] 0.9835862

SMS 分类器的 AUC 为 0.98，这非常高。但是我们怎么知道该模型是否同样能够在另一个数据集上表现良好呢？为了回答这些问题，我们需要更好地理解我们能够将模型的预测结果从测试数据外推的范围。

估计未来表现

一些 R 语言机器学习包在构建模型过程中会显示混淆矩阵和性能度量。这些统计数据的目的是提供关于模型重新替换误差的见解，这种误差发生在即使模型直接从训练数据构建，训练数据仍被错误预测时。这些信息可以用作粗略的诊断工具，以识别明显表现不佳的模型。

重新替换误差并不是一个非常有用的未来性能指标。例如，一个通过死记硬背完美分类每个训练实例并且零重新替换误差的模型，将无法将其预测泛化到从未见过的数据上。因此，训练数据上的错误率可能对模型未来的表现过于乐观。

与其依赖重新替换误差，更好的做法是评估模型在它尚未见过的数据上的表现。我们在前几章中使用了这种方法，将可用数据分为训练集和测试集。然而，在某些情况下，创建训练集和测试集并不总是理想的。例如，在只有一小部分数据的情况下，你可能不希望再进一步减少样本量。

幸运的是，还有其他方法可以估计模型在未见数据上的表现。我们用来计算性能度量的caret包也提供了多个函数来估计未来的表现。如果你正在跟随 R 语言代码示例并且尚未安装caret包，请安装它。你还需要将该包加载到 R 会话中，使用library(caret)命令。

保留法

我们在前几章中使用的将数据划分为训练集和测试集的过程被称为保留法。如下面的图示所示，训练集用于生成模型，然后将其应用于测试集以生成预测结果并进行评估。通常，大约三分之一的数据用于测试，三分之二用于训练，但这个比例可以根据可用数据的多少而有所不同。为了确保训练数据和测试数据没有系统性的差异，它们的样本会被随机地分配到两个组中。

为了使保留法真正准确地估计未来的性能，在任何时候都不应让测试数据集上的表现影响模型。很容易不知不觉地违反这一规则，通过反复测试选择最佳模型。例如，假设我们在训练数据上构建了多个模型，并选择了在测试数据上表现最好的那个。因为我们已经挑选了最佳结果，所以测试性能并不能公正地衡量在未见数据上的表现。

为了避免这个问题，最好将原始数据分割，以便除了训练数据集和测试数据集之外，还可以有一个验证数据集。验证数据集将用于迭代和优化所选择的模型，留出测试数据集仅在最后一步使用，用于报告未来预测的估计误差率。训练、测试和验证的典型分割比例是 50%、25%和 25%。

提示

一位敏锐的读者会注意到，前几章中使用了保留测试数据集来评估模型并提高模型性能。这样做是为了说明问题，但确实违反了前面提到的规则。因此，所展示的模型性能统计数据并不是对未来在未见数据上的性能的有效估计，整个过程应该更准确地称为验证。

创建保留样本的一种简单方法是使用随机数生成器将记录分配到不同的分区中。这个技术最早在第五章，分而治之 – 使用决策树和规则进行分类中，用于创建训练数据集和测试数据集。

提示

如果你希望跟随以下的示例进行操作，可以从 Packt 出版网站下载credit.csv数据集，并使用credit <- read.csv("credit.csv")命令将其加载到数据框中。

假设我们有一个名为 credit 的数据框，包含 1000 行数据。我们可以按如下方式将其分为三个分区。首先，使用runif()函数创建一个随机排序的行 ID 向量，范围从 1 到 1000。runif()函数默认生成指定数量的 0 到 1 之间的随机值。runif()函数的名称来源于随机均匀分布，这一点在第二章，管理和理解数据中已有讨论。

> random_ids <- order(runif(1000))

这里使用的order()返回一个向量，表示 1,000 个随机数的排名顺序。例如，order(c(0.5, 0.25, 0.75, 0.1))返回序列4 2 1 3，因为最小的数（0.1）排在第四，第二小的（0.25）排在第二，依此类推。

我们可以使用生成的随机 ID 将credit数据框分成包含 500、250 和 250 条记录的训练、验证和测试数据集：

> credit_train <- credit[random_ids[1:500], ]
> credit_validate <- credit[random_ids[501:750], ]
> credit_test <- credit[random_ids[751:1000], ]

留出抽样的一个问题是每个分区可能具有某些类的较大或较小比例。在某些情况下，特别是某些类在数据集中占很小比例的情况下，这可能导致某个类在训练数据集中被省略。这是一个重大问题，因为模型将无法学习该类。

为了减少这种情况发生的机会，可以使用一种称为分层随机抽样的技术。尽管从长远来看，随机样本应该包含与整个数据集相同比例的每个类值，但分层随机抽样确保随机分区几乎与整个数据集中每个类的比例相同，即使某些类很小。

caret包提供了一个createDataPartition()函数，该函数基于分层留出抽样创建分区。以下是为credit数据集创建分层样本的代码示例。要使用该函数，必须指定一个类值向量（这里的default指的是贷款是否违约），以及一个参数p，指定要包含在分区中的实例比例。list = FALSE参数防止结果以列表格式存储：

> in_train <- createDataPartition(credit$default, p = 0.75,
 list = FALSE)
> credit_train <- credit[in_train, ]
> credit_test <- credit[-in_train, ]

in_train向量指示包含在训练样本中的行号。我们可以使用这些行号从credit_train数据框中选择示例。类似地，通过使用负号，我们可以使用未在in_train向量中找到的行号来获取credit_test数据集。

尽管分层抽样可以均匀分布类别，但不能保证其他类型的代表性。一些样本可能包含太多或太少的难例、易于预测的案例或异常值。对于较小的数据集尤其如此，这些案例可能没有足够大的部分可以在训练和测试集之间划分。

除了可能存在偏倚的样本外，留出法的另一个问题是必须保留大量数据用于测试和验证模型。由于这些数据在评估模型性能之前不能用于训练模型，因此性能估计可能过于保守。

提示

由于通常在较大数据集上训练的模型表现更好，一种常见的做法是在选择和评估最终模型后，重新在整套数据（即训练加测试和验证）上训练模型。

一种叫做重复留出法的技术有时被用来减轻随机组成的训练数据集所带来的问题。重复留出法是留出法的一种特殊情况，它通过使用几个随机留出样本的平均结果来评估模型的表现。由于使用了多个留出样本，因此模型训练或测试时使用非代表性数据的可能性较小。我们将在下一节中详细讨论这个概念。

交叉验证

重复留出法是一个叫做k 折交叉验证（或k-fold CV）技术的基础，这已经成为估计模型表现的行业标准。与反复进行随机抽样（可能会多次使用相同的记录）不同，k 折交叉验证将数据随机划分为k个完全独立的随机分区，称为折叠。

尽管k可以设置为任何数字，但迄今为止，最常见的约定是使用10 折交叉验证（10-fold CV）。为什么是 10 个折叠？原因在于经验数据表明，使用更多折叠的好处并不显著。对于每个 10 折中的折叠（每个折叠包含总数据的 10%），机器学习模型会在剩余的 90%数据上进行训练。然后，使用与折叠匹配的 10%样本进行模型评估。经过 10 次训练和评估（10 种不同的训练/测试组合）后，报告所有折叠的平均表现。

提示

k 折交叉验证的一个极端情况是留一法，它使用数据的每个示例作为一个折叠来进行 k 折交叉验证。这确保了最大限度地利用数据来训练模型。尽管这看起来很有用，但由于计算成本极高，它在实际中很少使用。

数据集的交叉验证可以通过caret包中的createFolds()函数来创建。与分层随机留出采样类似，这个函数会尝试在每个折叠中保持与原始数据集相同的类别平衡。以下是创建 10 个折叠的命令：

> folds <- createFolds(credit$default, k = 10)

createFolds()函数的结果是一个向量列表，存储着每个请求的k = 10个折叠的行号。我们可以通过str()查看其内容：

> str(folds)
List of 10
 $ Fold01: int [1:100] 1 5 12 13 19 21 25 32 36 38 ...
 $ Fold02: int [1:100] 16 49 78 81 84 93 105 108 128 134 ...
 $ Fold03: int [1:100] 15 48 60 67 76 91 102 109 117 123 ...
 $ Fold04: int [1:100] 24 28 59 64 75 85 95 97 99 104 ...
 $ Fold05: int [1:100] 9 10 23 27 29 34 37 39 53 61 ...
 $ Fold06: int [1:100] 4 8 41 55 58 103 118 121 144 146 ...
 $ Fold07: int [1:100] 2 3 7 11 14 33 40 45 51 57 ...
 $ Fold08: int [1:100] 17 30 35 52 70 107 113 129 133 137 ...
 $ Fold09: int [1:100] 6 20 26 31 42 44 46 63 79 101 ...
 $ Fold10: int [1:100] 18 22 43 50 68 77 80 88 106 111 ...

在这里，我们看到第一个折叠被命名为Fold01，并存储了100个整数，表示第一个折叠中的信用数据框中的 100 行。为了创建训练集和测试集以构建和评估模型，还需要额外的步骤。以下命令展示了如何为第一个折叠创建数据。我们将选定的 10%分配给测试数据集，并使用负号将剩余的 90%分配给训练数据集：

> credit01_test <- credit[folds$Fold01, ]
> credit01_train <- credit[-folds$Fold01, ]

要执行完整的 10 次交叉验证（10-fold CV），此步骤需要重复 10 次；每次都要构建模型并计算模型的表现。最后，所有表现指标会被平均，以获得整体性能。幸运的是，我们可以通过应用之前学到的几种技术来自动化这一过程。

为了演示该过程，我们将使用 10 次交叉验证估计信用数据的 C5.0 决策树模型的卡帕统计量。首先，我们需要加载一些 R 包：caret（用于创建折叠）、C50（用于决策树）和 irr（用于计算卡帕）。后两个包是为了示例目的选择的；如果你愿意，也可以使用其他模型或其他性能指标，沿用相同的步骤。

> library(caret)
> library(C50)
> library(irr)

接下来，我们将像之前一样创建一个包含 10 个折叠的列表。这里使用 set.seed() 函数来确保如果再次运行相同的代码，结果是一致的：

> set.seed(123)
> folds <- createFolds(credit$default, k = 10)

最后，我们将使用 lapply() 函数对折叠列表应用一系列相同的步骤。如以下代码所示，由于没有现成的函数可以完美满足我们的需求，我们必须定义自己的函数并传递给 lapply()。我们自定义的函数将信用数据框分为训练数据和测试数据，使用 C5.0() 函数在训练数据上构建决策树，从测试数据生成一组预测，并使用 kappa2() 函数比较预测值和实际值：

> cv_results <- lapply(folds, function(x) {
 credit_train <- credit[-x, ]
 credit_test <- credit[x, ]
 credit_model <- C5.0(default ~ ., data = credit_train)
 credit_pred <- predict(credit_model, credit_test)
 credit_actual <- credit_test$default
 kappa <- kappa2(data.frame(credit_actual, credit_pred))$value
 return(kappa)
 })

结果的卡帕统计量被编译成一个列表，存储在 cv_results 对象中，我们可以使用 str() 函数检查它：

> str(cv_results)
List of 10
 $ Fold01: num 0.343
 $ Fold02: num 0.255
 $ Fold03: num 0.109
 $ Fold04: num 0.107
 $ Fold05: num 0.338
 $ Fold06: num 0.474
 $ Fold07: num 0.245
 $ Fold08: num 0.0365
 $ Fold09: num 0.425
 $ Fold10: num 0.505

10 次交叉验证过程只剩下最后一步：我们必须计算这 10 个值的平均值。虽然你可能会倾向于输入 mean(cv_results)，但是由于 cv_results 不是一个数值向量，这样会导致错误。相反，应该使用 unlist() 函数，它可以消除列表结构，将 cv_results 转化为数值向量。然后，我们就可以像预期的那样计算平均卡帕值：

> mean(unlist(cv_results))
[1] 0.283796

这个卡帕统计量相对较低，对应于解释尺度中的“公平”，这表明信用评分模型的表现仅略优于随机猜测。在下一章中，我们将研究基于 10 次交叉验证的自动化方法，帮助我们提升该模型的表现。

提示

也许目前最可靠的模型性能估计方法是 重复的 k 折交叉验证（repeated k-fold CV）。正如你从名字中猜到的那样，这涉及到反复应用 k 折交叉验证，并平均结果。一种常见的策略是进行 10 次 10 折交叉验证。尽管这种方法计算量大，但它提供了一个非常稳健的估计。

自助抽样

一种使用频率稍低的替代方法是自助采样（bootstrap sampling），简称bootstrap或自助法。一般而言，这些指的是使用数据的随机样本来估计较大数据集属性的统计方法。当这个原理应用于机器学习模型的性能时，它意味着创建几个随机选择的训练集和测试集，然后用它们来估计性能统计数据。来自不同随机数据集的结果会被平均，最终得到对未来性能的估计。

那么，这个过程与 k 折交叉验证有何不同呢？交叉验证将数据划分为不同的分区，其中每个示例只能出现一次，而自助法通过有放回采样允许示例被多次选择。这意味着，从原始的n个示例的数据集中，自助法程序将创建一个或多个新的训练数据集，这些数据集也包含n个示例，其中一些是重复的。相应的测试数据集则从未被选入相应训练数据集的示例集中构建。

使用前面描述的带有替代的采样方法，任何给定实例被包含在训练数据集中的概率为 63.2%。因此，任何实例出现在测试数据集中的概率为 36.8%。换句话说，训练数据仅代表了 63.2%的可用示例，其中一些是重复的。与使用 90%示例用于训练的 10 折交叉验证（10-fold CV）相比，自助法样本对完整数据集的代表性较差。

因为仅在 63.2%的训练数据上训练的模型可能会比在更大训练集上训练的模型表现差，所以自助法的性能估计可能比后来在完整数据集上训练模型时的结果低得多。自助法的一个特例，称为0.632 自助法，通过将最终的性能度量计算为训练数据（过于乐观）和测试数据（过于悲观）性能的函数来考虑这一点。最终的误差率可以通过以下方式估算：

自助法相较于交叉验证的一个优点是它在非常小的数据集上表现得更好。此外，自助采样不仅用于性能测量，还有其他应用。特别是在下一章中，我们将学习如何利用自助采样的原理来提高模型性能。

摘要

本章介绍了评估机器学习分类模型性能的几种常见度量和技术。尽管准确度提供了一种简单的方法来检查模型的正确性，但在稀有事件的情况下，这可能会误导，因为此类事件的实际成本可能与它们出现的频率成反比。

基于混淆矩阵的多个度量方法更好地捕捉了各种错误类型成本之间的平衡。仔细审视灵敏度与特异性，或者精确度与召回率之间的权衡，可能是思考现实世界中错误影响的一个有用工具。像 ROC 曲线这样的可视化工具也有助于这一点。

还值得一提的是，有时评估一个模型性能的最佳标准是考虑它在满足或未能满足其他目标方面的表现。例如，你可能需要用简单的语言解释模型的逻辑，这将排除一些模型的考虑范围。此外，即使一个模型表现非常好，如果它太慢或难以在生产环境中扩展，那么它也是完全没有用的。

测量性能的一个明显扩展是找到自动化方法，为特定任务寻找最佳模型。在下一章中，我们将基于目前的工作，研究通过系统地迭代、优化和结合学习算法来构建更智能模型的方法。

第十一章：改进模型性能

当一支运动队未能达到其目标时——无论目标是获得奥运金牌、联赛冠军还是世界纪录——它必须寻找可能的改进方向。假设你是该队的教练，你会如何安排训练？或许你会指示运动员更加努力训练或改变训练方式，以最大化他们的潜力。或者，你可能会强调更好的团队合作，更聪明地利用运动员的长处和短处。

现在，假设你正在训练一款世界级的机器学习算法。也许你希望参加数据挖掘竞赛，比如 Kaggle 上发布的竞赛（www.kaggle.com/competitions）。或者你仅仅是希望改善商业成果。你该从哪里开始？尽管背景不同，但提升运动队表现的策略同样可以用来提升统计学习器的表现。

作为教练，你的任务是找到训练技巧和团队协作技能的组合，以帮助你实现性能目标。本章在本书中所覆盖的内容的基础上，介绍了一组提高机器学习器预测性能的技巧。你将学习：

• 如何通过系统地寻找最佳训练条件的组合来自动化模型性能调优

• 将模型组合成利用团队合作解决困难学习任务的方法

• 如何应用一种变体的决策树，这种决策树因其出色的表现而迅速流行

这些方法并非对每个问题都有效。然而，从机器学习竞赛的获奖作品来看，你可能会发现其中至少使用了某种方法。为了具有竞争力，你也需要将这些技能纳入你的技能库。

调优标准模型以提高性能

一些学习问题非常适合前几章中介绍的标准模型。在这种情况下，可能不需要花费太多时间反复调整和优化模型；它可能已经足够好。然而，另一方面，有些问题本质上更为复杂。需要学习的基本概念可能极为复杂，涉及许多微妙的关系，或者它可能受到随机变化的影响，使得在噪音中定义信号变得困难。

开发在困难问题上表现极佳的模型既是一门艺术，也是一门科学。在尝试找出性能提升的方向时，直觉有时是有帮助的。在其他情况下，找到提升的方法可能需要一种蛮力式的反复试验方法。当然，借助自动化程序，搜索可能的改进方法的过程可以得到帮助。

在第五章，分治法——使用决策树和规则进行分类中，我们尝试了解决一个难题：识别可能进入违约的贷款。尽管我们通过性能调优方法获得了大约 82%的可接受分类准确率，但在第十章，模型性能评估中仔细检查后，我们发现高准确率有些误导。尽管准确率合理，但 Kappa 统计量只有大约 0.28，表明模型的实际表现并不理想。在这一节中，我们将重新审视信用评分模型，看看是否可以改善结果。

提示

要跟随示例，下载 Packt Publishing 网站上的 credit.csv 文件，并将其保存到你的 R 工作目录中。使用命令 credit <- read.csv("credit.csv") 将文件加载到 R 中。

你可能还记得，我们首先使用了一个标准的 C5.0 决策树来构建信用数据的分类器。然后，我们尝试通过调整 trials 参数来增加提升迭代次数，从而提高模型性能。通过将迭代次数从默认值 1 增加到 10，我们成功提高了模型的准确度。这个调整模型选项以找出最佳拟合的过程称为参数 调优。

参数调优不限于决策树。例如，当我们搜索最佳 k 值时，我们对 k-NN 模型进行了调优。当我们调整神经网络和支持向量机的节点数或隐藏层数，或者选择不同的核函数时，我们也进行了调优。大多数机器学习算法允许调整至少一个参数，而最复杂的模型提供了大量调节模型拟合的方法。尽管这使得模型能够更好地适应学习任务，但所有可能选项的复杂性可能会让人感到压倒。此时，更系统的方法是必要的。

使用 caret 进行自动化参数调优

与其为每个模型的参数选择任意的值——这不仅是繁琐的，而且有些不科学——不如通过搜索多个可能的参数值来找到最佳组合。

caret 包，我们在第十章，模型性能评估中广泛使用，提供了帮助自动化参数调优的工具。核心功能由 train() 函数提供，该函数作为一个标准化接口，支持超过 175 种不同的机器学习模型，用于分类和回归任务。通过使用这个函数，可以通过选择不同的评估方法和指标，自动化地搜索最优模型。

提示

不要被大量模型吓到——我们在前面的章节中已经介绍了很多模型。其他的模型只是基础概念的简单变体或扩展。考虑到你目前所学的内容，你应该有信心能理解所有可用的方法。

自动调优参数需要你考虑三个问题：

• 应该在数据上训练什么类型的机器学习模型（以及具体的实现）？

• 哪些模型参数可以调整，应该如何调节这些参数以找到最佳设置？

• 应该使用什么标准来评估模型，以找到最佳候选模型？

回答第一个问题需要在机器学习任务和 175 个模型之间找到一个合适的匹配。显然，这需要对机器学习模型的广度和深度有一定了解。进行排除法也会有所帮助。根据任务是分类还是数值预测，几乎一半的模型可以被排除；其他的可以根据数据的格式或是否需要避免使用黑箱模型来排除。无论如何，也没有理由不能尝试多种方法并比较每种方法的最佳结果。

解决第二个问题在很大程度上取决于模型的选择，因为每个算法使用一组独特的参数。本书中涵盖的预测模型的可用调优参数列在下表中。请记住，虽然一些模型可能有未显示的额外选项，但caret仅支持表中列出的选项进行自动调优。

模型	学习任务	方法名称	参数
k 近邻算法	分类	knn	k
朴素贝叶斯	分类	nb	fL, usekernel
决策树	分类	C5.0	model, trials, winnow
OneR 规则学习器	分类	OneR	无
RIPPER 规则学习器	分类	JRip	NumOpt
线性回归	回归	lm	无
回归树	回归	rpart	cp
模型树	回归	M5	pruned, smoothed, rules
神经网络	双重用途	nnet	size, decay
支持向量机（线性核）	双重用途	svmLinear	C
支持向量机（径向基核）	双重用途	svmRadial	C, sigma
随机森林	双重用途	rf	mtry

提示

要查看caret所涵盖的模型及其调优参数的完整列表，请参考包作者 Max Kuhn 提供的表格：topepo.github.io/caret/modelList.html。

如果你忘记了某个模型的调优参数，可以使用modelLookup()函数来查找它们。只需提供方法名称，以下是 C5.0 模型的示例：

> modelLookup("C5.0")
 model parameter                 label forReg forClass probModel
1  C5.0    trials # Boosting Iterations  FALSE     TRUE      TRUE
2  C5.0     model            Model Type  FALSE     TRUE      TRUE
3  C5.0    winnow                Winnow  FALSE     TRUE      TRUE

自动调整的目标是搜索一个候选模型集合，这些模型由一组参数组合的矩阵或网格构成。由于不现实去搜索所有可能的组合，因此只使用部分可能性来构建网格。默认情况下，caret 每个 p 参数最多搜索三个值。这意味着最多会测试 3^p 个候选模型。例如，默认情况下，k-最近邻的自动调整会比较 3¹ = 3 个候选模型，分别对应 k=5、k=7 和 k=9。类似地，调整决策树将比较最多 27 个不同的候选模型，这些模型由 3³ = 27 个 model、trials 和 winnow 设置的组合构成。然而，实际上，只会测试 12 个模型。这是因为 model 和 winnow 参数只能取两个值（分别是 tree 与 rules 和 TRUE 与 FALSE），所以网格大小为 3 * 2 * 2 = 12。

提示

由于默认的搜索网格可能并不适合你的学习问题，caret 允许你通过一个简单的命令提供自定义的搜索网格，我们将在后面介绍。

自动调整模型的第三步也是最后一步，涉及从候选模型中识别最佳模型。这个过程使用了在第十章中讨论的方法，即通过选择重采样策略来创建训练集和测试集，并使用模型性能统计量来衡量预测准确性。

所有我们学到的重采样策略和许多性能统计量都得到了 caret 的支持。这些包括精度、Kappa（对于分类器）以及 R-squared 或 RMSE（对于数值模型）等统计量。如果需要，还可以使用如灵敏度、特异性和 ROC 曲线下的面积（AUC）等成本敏感度度量。

默认情况下，caret 会选择具有最大性能度量值的候选模型。由于这种做法有时会选择那些通过大幅增加模型复杂度来实现边际性能提升的模型，因此提供了替代的模型选择函数。

鉴于有各种各样的选项，许多默认设置是合理的，这一点非常有帮助。例如，caret 将使用在自助法样本上的预测准确度来选择最佳的分类模型表现者。通过这些默认值开始，我们可以调整 train() 函数来设计各种各样的实验。

创建一个简单的调整模型

为了说明调整模型的过程，我们首先观察一下当我们尝试使用 caret 包的默认设置来调整信用评分模型时会发生什么。接下来，我们将根据需要调整选项。

调整学习器的最简单方法只需通过method参数指定模型类型即可。由于我们之前在credit模型中使用了 C5.0 决策树，因此我们将继续通过优化此学习器来进行后续工作。使用默认设置调整 C5.0 决策树的基本train()命令如下：

> library(caret)
> set.seed(300)
> m <- train(default ~ ., data = credit, method = "C5.0")

首先，使用set.seed()函数初始化 R 的随机数生成器，以一个设定的起始位置。你可能还记得我们在之前的章节中使用过这个函数。通过设置seed参数（在这里设为任意的数字 300），随机数将按照预定的序列生成。这使得使用随机抽样的模拟可以重复进行并得到相同的结果——这对于共享代码或重复先前的结果非常有帮助。

接下来，我们使用 R 公式接口定义一个树模型，表示为default ~ .。这个模型使用credit数据框中的所有其他特征来建模贷款违约状态（yes或no）。参数method = "C5.0"告诉caret使用 C5.0 决策树算法。

在输入前述命令后，可能会有一个较长的延迟（具体取决于你计算机的性能），因为调整过程正在进行。尽管这是一个相对较小的数据集，但仍需要进行大量的计算。R 必须反复生成随机数据样本，构建决策树，计算性能统计量并评估结果。

实验的结果保存在一个名为m的对象中。如果你想检查该对象的内容，str(m)命令会列出所有相关数据，但这可能会显得有些繁杂。相反，只需输入对象的名称，即可获得一个简洁的结果概览。例如，输入m会产生以下输出（注意，为了清晰起见，已添加标签）：

标签突出显示了输出中的四个主要组成部分：

1. 输入数据集的简要描述：如果你对自己的数据比较熟悉，并且正确应用了train()函数，这些信息应该不会让你感到惊讶。

2. 应用的预处理和重采样方法报告：在这里，我们看到使用了 25 个 bootstrap 样本，每个样本包括 1,000 个示例，用于训练模型。

3. 评估的候选模型列表：在这一部分，我们可以确认测试了 12 个不同的模型，这些模型是基于三个 C5.0 调整参数的组合——model、trials和winnow。每个候选模型的准确性和卡帕统计量的平均值和标准差也在这里显示。

4. 最佳模型的选择：正如脚注所述，选择了准确性最高的模型。这个模型使用了一个有 20 次试验的决策树，并设置winnow = FALSE。

在确定最佳模型后，train() 函数使用其调优参数在完整的输入数据集上构建模型，并将其存储在 m 列表对象中的 m$finalModel。在大多数情况下，你不需要直接使用 finalModel 子对象。相反，只需按如下方式使用 predict() 函数和 m 对象：

> p <- predict(m, credit)

生成的预测向量按预期工作，使我们能够创建一个比较预测值和实际值的混淆矩阵：

> table(p, credit$default)

p      no yes
 no  700   2
 yes   0 298

在用于训练最终模型的 1,000 个示例中，只有两个被错误分类。然而，非常重要的是要注意，由于模型是基于训练数据和测试数据构建的，因此该准确性是乐观的，不能作为未见数据的性能指标。总结输出中的 73% 自助法估计值是对未来表现的更现实估计。

使用 train() 和 predict() 函数除了自动参数调优外，还提供了几个额外的好处。

首先，train() 函数应用的任何数据准备步骤将同样应用于用于生成预测的数据。这包括像中心化、缩放以及缺失值填充等变换。让 caret 处理数据准备可以确保那些对最佳模型性能有贡献的步骤在模型部署时得以保留。

其次，predict() 函数提供了一个标准化接口，用于获取预测的类别值和类别概率，即使是那些通常需要额外步骤才能获取这些信息的模型类型。预测类别默认会被提供：

> head(predict(m, credit))
[1] no  yes no  no  yes no
Levels: no yes

要获得每个类别的估计概率，请使用 type = "prob" 参数：

> head(predict(m, credit, type = "prob"))
 no        yes
1 0.9606970 0.03930299
2 0.1388444 0.86115561
3 1.0000000 0.00000000
4 0.7720279 0.22797208
5 0.2948062 0.70519385
6 0.8583715 0.14162851

即使在底层模型使用不同字符串表示预测概率的情况下（例如，naiveBayes 模型使用 "raw"），predict() 函数也会在幕后将 type = "prob" 转换为适当的字符串。

自定义调优过程

我们之前创建的决策树展示了 caret 包在最小干预下生成优化模型的能力。默认设置使得优化模型的创建变得简单。然而，也可以更改默认设置，使其更符合学习任务的具体需求，这有助于解锁更高水平的性能。

模型选择过程中的每个步骤都可以定制。为了说明这一灵活性，让我们修改在信用决策树中的工作，使其与我们在第十章中使用的过程相似，评估模型性能。如果你还记得，我们使用 10 折交叉验证估算了 kappa 统计量。我们将在这里做同样的事情，使用 kappa 来优化决策树的提升参数。请注意，决策树提升在第五章中已有介绍，并将在本章后面更详细地讨论。

trainControl()函数用于创建一组配置选项，称为控制对象，它指导train()函数的执行。这些选项允许管理模型评估标准，例如重采样策略和用于选择最佳模型的度量标准。尽管此函数可以用来修改调优实验的几乎每个方面，我们将重点关注两个重要的参数：method和selectionFunction。

提示

如果你渴望获取更多细节，可以使用?trainControl命令查看所有参数的列表。

对于trainControl()函数，method参数用于设置重采样方法，如留出法或 k 折交叉验证。下表列出了所有可能的重采样方法类型以及调整样本大小和迭代次数的其他参数。尽管这些重采样方法的默认选项遵循了常见的惯例，但你可以根据数据集的大小和模型的复杂性调整这些选项。

重采样方法	方法名称	额外选项及默认值
留出法采样	LGOCV	p = 0.75（训练数据比例）
k 折交叉验证	cv	number = 10（折数）
重复 k 折交叉验证	repeatedcv	number = 10（折数）repeats = 10（迭代次数）
自助法采样	boot	number = 25（重采样迭代次数）
0.632 自助法	boot632	number = 25（重采样迭代次数）
留一交叉验证	LOOCV	无

selectionFunction 参数用于指定将在各种候选模型中选择最佳模型的函数。包括三种此类函数。best 函数简单地选择在指定性能度量上表现最好的候选模型，这是默认使用的函数。其他两个函数用于选择在性能上接近最佳模型的最简模型。oneSE 函数选择在最佳表现的一个标准误差范围内最简单的候选模型，而 tolerance 函数则在用户指定的百分比范围内选择最简单的候选模型。

提示

caret 包根据简化性对模型进行排名时涉及一些主观性。如需了解如何对模型进行排名，请在 R 命令提示符下输入 ?best，查看选择函数的帮助页面。

若要创建一个名为 ctrl 的控制对象，该对象使用 10 折交叉验证和 oneSE 选择函数，请使用以下命令（请注意，number = 10 仅为清晰起见包含，因为这是 method = "cv" 的默认值，可以省略）：

> ctrl <- trainControl(method = "cv", number = 10,
 selectionFunction = "oneSE")

我们很快将使用这个函数的结果。

与此同时，定义实验的下一步是创建要优化的参数网格。该网格必须包含每个期望模型参数命名的列，并以句点作为前缀。它还必须包含每个期望参数值组合的行。由于我们使用的是 C5.0 决策树，这意味着我们需要名为 .model、.trials 和 .winnow 的列。对于其他机器学习模型，请参考本章前面提供的表格，或者使用 modelLookup() 函数根据之前的描述查找参数。

与其一一填写数据框的每个单元格——如果有许多参数值组合，这将是一项繁琐的任务——我们可以使用 expand.grid() 函数，它可以根据提供的所有值的组合创建数据框。例如，假设我们希望保持 model = "tree" 和 winnow = "FALSE" 不变，同时搜索八种不同的 trials 值。这可以这样创建：

> grid <- expand.grid(.model = "tree",
 .trials = c(1, 5, 10, 15, 20, 25, 30, 35),
 .winnow = "FALSE")

生成的网格数据框包含 1 * 8 * 1 = 8 行：

> grid
 .model .trials .winnow
1   tree       1   FALSE
2   tree       5   FALSE
3   tree      10   FALSE
4   tree      15   FALSE
5   tree      20   FALSE
6   tree      25   FALSE
7   tree      30   FALSE
8   tree      35   FALSE

train() 函数将使用每行的模型参数组合构建一个候选模型进行评估。

给定这个搜索网格和先前创建的控制列表，我们准备运行一个彻底定制的 train() 实验。像之前一样，我们将随机种子设置为任意数字 300，以确保结果可重复。但这次，我们将在传递控制对象和调参网格的同时，添加一个参数 metric = "Kappa"，指示模型评估函数要使用的统计量——在这种情况下为 "oneSE"。完整命令如下：

> set.seed(300)
> m <- train(default ~ ., data = credit, method = "C5.0",
 metric = "Kappa",
 trControl = ctrl,
 tuneGrid = grid)

这将生成一个对象，我们可以通过输入其名称查看：

> m

尽管大部分输出与自动调优后的模型相似，但仍有一些值得注意的差异。由于使用了 10 折交叉验证，构建每个候选模型的样本大小被减少到了 900，而不是自助法中的 1,000。根据我们的要求，测试了八个候选模型。此外，由于model和winnow保持不变，它们的值不再显示在结果中，而是作为脚注列出。

这里的最佳模型与之前的实验差异明显。之前，最佳模型使用了trials = 20，而这里使用的是trials = 1。这个看似奇怪的发现是因为我们使用了oneSE规则，而不是best规则来选择最佳模型。尽管 35 次实验模型根据 kappa 值提供了最好的原始性能，但 1 次实验模型在性能上几乎相同，而且形式更加简单。简单模型不仅计算效率更高，而且还能减少过拟合训练数据的可能性。

通过元学习提高模型性能

作为提高单一模型性能的替代方法，可以将多个模型组合成一个强大的团队。就像最优秀的运动队拥有互补而非重叠技能的球员一样，一些最好的机器学习算法利用互补模型的团队。由于每个模型都会为学习任务带来独特的偏差，它可能很容易学习某个子集的样本，但对另一个子集则可能表现不佳。因此，通过智能地利用多个不同团队成员的优势，可以创建一个由多个弱学习者组成的强大团队。

结合和管理多个模型预测的技术属于元学习方法的一部分，定义了涉及学习如何学习的技术。这包括从简单的算法（通过反复设计决策逐步改进性能——例如，本章早些时候提到的自动参数调优）到使用借鉴自进化生物学和遗传学的概念进行自我修改和适应学习任务的高度复杂算法。

在本章剩余部分，我们将专注于元学习，仅限于建模多个模型预测与期望结果之间关系的内容。本节中介绍的基于团队合作的技术非常强大，并且在构建更有效的分类器时被广泛使用。

理解集成模型

假设你是电视答题节目的一名参赛者，可以选择五个朋友组成团队来帮助你回答最终的百万美元大奖问题。大多数人会试图选取一组多样化的学科专家。一个包含文学、科学、历史、艺术教授以及当代流行文化专家的团队，将是一个均衡的团队。考虑到他们的知识广度，几乎不可能有一个问题能让这个团队感到难倒。

利用类似于创建一个多样化专家团队的原理的元学习方法被称为集成方法。所有的集成方法都基于这样一个理念：通过将多个较弱的学习器组合起来，创造出一个更强的学习器。各种集成方法的区别，主要可以通过以下两个问题的答案来区分：

• 如何选择和/或构建弱学习模型？

• 如何将弱学习器的预测结果组合成一个最终的预测？

在回答这些问题时，想象集成方法的过程图可能会很有帮助；几乎所有的集成方法都遵循这个模式：

首先，使用输入的训练数据来构建多个模型。分配函数决定了每个模型接收多少训练数据。它们是否每个都接收到完整的训练数据集，还是仅仅接收到一个样本？它们是否每个都接收到所有特征，还是仅接收到一部分特征？

虽然理想的集成方法包括多样化的模型集，但分配函数可以通过人为改变输入数据来增加多样性，从而使得生成的学习器产生偏差，即使它们是同一类型的。例如，它可能使用自助抽样（bootstrap sampling）来构建独特的训练数据集，或者将不同的特征或样本子集传递给每个模型。另一方面，如果集成方法已经包含了多种算法——如神经网络、决策树和 k-NN 分类器——那么分配函数可能会将数据传递给每个算法，而数据保持相对不变。

在模型构建完成后，它们可以用于生成一组预测结果，这些预测结果必须以某种方式进行管理。组合函数决定了如何解决预测之间的分歧。例如，集成方法可能会使用多数投票来确定最终的预测结果，或者使用更复杂的策略，比如根据每个模型的历史表现来加权每个模型的投票。

一些集成方法甚至使用另一个模型来学习从各种预测组合中得到一个组合函数。例如，假设当M1和M2都投票“是”时，实际的类别值通常是“否”。在这种情况下，集成方法可以学习忽略M1和M2的投票，当它们一致时。这种使用多个模型的预测结果来训练最终裁定模型的过程被称为堆叠。

使用集成方法的一个好处是，它们可能让你在追求单一最佳模型时花费更少的时间。你可以训练多个合理强大的候选模型并将其结合起来。然而，方便性并不是集成方法在机器学习竞赛中持续获胜的唯一原因；集成方法在多个方面也提供了相对于单一模型的性能优势：

• 更好的泛化能力：由于多个学习者的意见被整合到最终的预测中，因此没有任何单一的偏差能够主导预测结果。这减少了过拟合学习任务的风险。

• 在大规模或极小数据集上提升性能：许多模型在使用非常大规模的特征或样本集时会遇到内存或复杂度的限制，这时训练多个小模型比训练一个完整的模型更高效。相反，集成方法在最小的数据集上也表现良好，因为许多集成设计本身就包含了如自助抽样（bootstrapping）等重采样方法。或许最重要的是，集成方法通常可以通过分布式计算方法并行训练。

• 合成来自不同领域的数据的能力：由于没有一种适用于所有情况的学习算法，集成方法能够结合来自多种学习者的证据，这在复杂现象依赖于来自不同领域的数据时变得越来越重要。

• 对困难学习任务的更细致理解：现实世界中的现象往往非常复杂，包含许多相互作用的细节。将任务划分为更小部分的模型，往往能更准确地捕捉到单一全局模型可能忽略的微妙模式。

如果你无法轻松地在 R 中应用集成方法，那么这些好处将大打折扣，幸运的是，已经有许多包可以用来实现这一点。我们来看一下几种最流行的集成方法，以及它们如何帮助提升我们正在研究的信用模型的性能。

自助法（Bagging）

第一个获得广泛认可的集成方法使用了一种叫做自助聚合（bootstrap aggregating），简称Bagging的方法。正如 Leo Breiman 在 1994 年所描述的那样，自助法通过对原始训练数据进行自助抽样生成多个训练数据集。这些数据集随后被用来生成一组模型，使用同一个学习算法。这些模型的预测结果通过投票（分类）或平均（数值预测）进行结合。

注意事项

关于自助法的更多信息，请参阅 Breiman L. Bagging predictors. Machine Learning. 1996; 24:123-140.

虽然袋装是一种相对简单的集成方法，但只要与相对不稳定的学习器一起使用，它可以表现得相当好。所谓不稳定学习器，就是那些生成的模型在输入数据发生轻微变化时会发生显著变化的模型。为了确保集成的多样性，即使是从 bootstrap 训练数据集之间仅有微小的变化，不稳定的模型是至关重要的。正因如此，袋装方法常常与决策树一起使用，因为决策树在输入数据发生微小变化时往往会发生剧烈变化。

ipred包提供了经典的袋装决策树实现。为了训练模型，bagging()函数的工作方式与之前使用的许多模型类似。nbagg参数用于控制在集成中投票的决策树数量（默认值为25）。根据学习任务的难度和训练数据的数量，增加这个数量可能会提高模型的性能，但也有一个限制。缺点是，这会带来额外的计算开销，因为训练大量的树可能需要一些时间。

安装ipred包后，我们可以按如下方式创建集成。我们将保持默认的 25 棵决策树：

> library(ipred)
> set.seed(300)
> mybag <- bagging(default ~ ., data = credit, nbagg = 25)

生成的模型按预期工作，可以使用predict()函数：

> credit_pred <- predict(mybag, credit)
> table(credit_pred, credit$default)

credit_pred  no yes
 no  699   2
 yes   1 298

根据之前的结果，该模型似乎非常适合训练数据。为了查看这一点如何转化为未来的表现，我们可以使用带有 10 倍交叉验证的袋装树，并使用caret包中的train()函数。请注意，ipred袋装树函数的方法名是treebag：

> library(caret)
> set.seed(300)
> ctrl <- trainControl(method = "cv", number = 10)
> train(default ~ ., data = credit, method = "treebag",
 trControl = ctrl)

Bagged CART 

1000 samples
 16 predictor
 2 classes: 'no', 'yes' 

No pre-processing
Resampling: Cross-Validated (10 fold) 

Summary of sample sizes: 900, 900, 900, 900, 900, 900, ... 

Resampling results

 Accuracy  Kappa      Accuracy SD  Kappa SD 
 0.735     0.3297726  0.03439961   0.08590462

该模型的 kappa 统计量为 0.33，表明集成树模型的表现至少与我们在本章早些时候调优的最佳 C5.0 决策树相当。这说明了集成方法的强大；一组简单的学习器协同工作可以超越非常复杂的模型。

为了超越决策树的袋装，caret包还提供了一个更通用的bag()函数。它原生支持一些模型，尽管通过一些额外的努力，它可以适配到其他类型的模型。bag()函数使用控制对象来配置袋装过程。它需要指定三个函数：一个用于拟合模型，一个用于做出预测，一个用于聚合投票。

例如，假设我们想要创建一个袋装支持向量机模型，可以使用我们在第七章中使用的kernlab包中的ksvm()函数。bag()函数要求我们提供训练 SVM、做出预测和统计投票的功能。

我们无需自己编写这些功能，caret包内置的svmBag列表对象提供了三个我们可以使用的函数：

> str(svmBag)
List of 3
 $ fit      :function (x, y, ...) 
 $ pred     :function (object, x) 
 $ aggregate:function (x, type = "class")

通过查看svmBag$fit函数，我们可以看到它只是调用了kernlab包中的ksvm()函数并返回结果：

> svmBag$fit
function (x, y, ...) 
{
 library(kernlab)
 out <- ksvm(as.matrix(x), y, prob.model = is.factor(y), ...)
 out
}
<environment: namespace:caret>

svmBag的pred和aggregate函数也同样简单。通过研究这些函数并以相同的格式创建自己的函数，可以使用袋装法与任何你想要的机器学习算法。

提示

caret包还包括了朴素贝叶斯模型袋（nbBag）、决策树（ctreeBag）和神经网络（nnetBag）的示例对象。

通过应用svmBag列表中的三个函数，我们可以创建一个袋装控制对象：

> bagctrl <- bagControl(fit = svmBag$fit, 
 predict = svmBag$pred,
 aggregate = svmBag$aggregate)

通过将此与train()函数和之前定义的训练控制对象（ctrl）一起使用，我们可以如下评估袋装 SVM 模型（请注意，kernlab包是必需的，如果尚未安装，需要先安装它）：

> set.seed(300)
> svmbag <- train(default ~ ., data = credit, "bag",
 trControl = ctrl, bagControl = bagctrl)
> svmbag

Bagged Model
1000 samples
 16 predictors
 2 classes: 'no', 'yes' 

No pre-processing
Resampling: Cross-Validation (10 fold) 

Summary of sample sizes: 900, 900, 900, 900, 900, 900, ... 

Resampling results

 Accuracy  Kappa      Accuracy SD  Kappa SD 
 0.728     0.2929505  0.04442222   0.1318101

Tuning parameter 'vars' was held constant at a value of 35

由于 kappa 统计量低于 0.30，似乎袋装 SVM 模型的表现不如袋装决策树模型。值得指出的是，与袋装决策树模型相比，kappa 统计量的标准差相当大。这表明，在交叉验证的各个折叠中，性能变化很大。这种变化可能意味着通过增加集成中的模型数量，性能可能会进一步提高。

提升法

另一种常见的基于集成的方法称为提升法，因为它通过提升弱学习器的表现来达到强学习器的表现。这种方法主要基于 Robert Schapire 和 Yoav Freund 的工作，他们在这个主题上发表了大量的研究。

注意

有关提升法的更多信息，请参考 Schapire RE, Freund Y. Boosting: Foundations and Algorithms。剑桥，马萨诸塞州，麻省理工学院出版社；2012 年。

类似于袋装法，提升法（boosting）使用在重采样数据上训练的模型集和投票来确定最终预测。这里有两个关键的区别。首先，提升法中的重采样数据集是特别构建的，目的是生成互补的学习器。其次，提升法不是给每个学习器平等的投票，而是根据其过去的表现为每个学习器的投票加权。表现较好的模型在集成中的最终预测中有更大的影响。

提升法通常会产生比集成中最好的模型更优的表现，且绝不会差于最强的模型。由于集成中的模型是为了互补性而构建的，假设每个分类器的表现优于随机机会，理论上可以通过增加额外的分类器来随意提高集成的性能。考虑到这一发现的明显实用性，提升法被认为是机器学习领域最重要的发现之一。

提示

尽管提升方法可以创建一个满足任意低误差率的模型，但在实践中这可能并不总是合理的。首先，随着更多学习器的加入，性能提升会越来越小，这使得某些阈值在实际中不可行。此外，追求纯粹的准确度可能会导致模型过拟合训练数据，而无法推广到未见过的数据。

一种名为AdaBoost或自适应增强的提升算法由 Freund 和 Schapire 于 1997 年提出。该算法基于生成弱学习器的思想，通过对频繁被误分类的样本给予更多关注（即，赋予更大的权重），迭代地学习更多难以分类的样本。

从一个未加权的数据集开始，第一个分类器尝试对结果进行建模。分类器正确预测的样本将不太可能出现在下一个分类器的训练数据集中，反之，难以分类的样本将更频繁地出现。随着更多轮弱学习器的加入，它们在越来越难分类的样本上进行训练。该过程持续进行，直到达到预期的整体误差率或性能不再提高为止。此时，每个分类器的投票将根据其在训练数据上的准确性进行加权。

虽然提升原理几乎可以应用于任何类型的模型，但这些原理最常见的应用是与决策树一起使用。我们在第五章中，分治法 – 使用决策树和规则进行分类，已经使用提升方法来提高 C5.0 决策树的性能。

AdaBoost.M1算法提供了另一个基于树的 AdaBoost 实现，用于分类。AdaBoost.M1 算法可以在adabag包中找到。

注意

如需了解更多关于adabag包的信息，请参考 Alfaro E, Gamez M, Garcia N 的文章。adabag – an R package for classification with boosting and bagging。统计软件杂志。2013；54：1-35。

让我们为信用数据创建一个AdaBoost.M1分类器。该算法的一般语法与其他建模技术类似：

> set.seed(300)
> m_adaboost <- boosting(default ~ ., data = credit)

像往常一样，predict()函数应用于结果对象以进行预测：

> p_adaboost <- predict(m_adaboost, credit)

与传统方法不同，该方法并不是返回预测结果的向量，而是返回一个包含模型信息的对象。预测结果存储在名为class的子对象中：

> head(p_adaboost$class)
[1] "no"  "yes" "no"  "no"  "yes" "no"

混淆矩阵可以在confusion子对象中找到：

> p_adaboost$confusion
 Observed Class
Predicted Class  no yes
 no  700   0
 yes   0 300

你注意到 AdaBoost 模型没有犯错吗？在你兴奋之前，记住之前的混淆矩阵是基于模型在训练数据上的表现。由于提升方法允许将错误率降低到任意低的水平，学习器会继续训练直到不再犯错。这很可能导致了在训练数据集上的过拟合。

为了更准确地评估在未见数据上的表现，我们需要使用另一种评估方法。adabag包提供了一个简单的函数来使用 10 折交叉验证（10-fold CV）：

> set.seed(300)
> adaboost_cv <- boosting.cv(default ~ ., data = credit)

根据你计算机的性能，这可能需要一些时间，期间它会在屏幕上记录每次迭代。完成后，我们可以查看一个更合理的混淆矩阵：

> adaboost_cv$confusion
 Observed Class
Predicted Class  no yes
 no  594 151
 yes 106 149

我们可以使用vcd包找到 kappa 统计量，具体方法参见第十章，评估模型性能。

> library(vcd)
> Kappa(adaboost_cv$confusion)
 value       ASE
Unweighted 0.3606965 0.0323002
Weighted   0.3606965 0.0323002

该模型的 kappa 值约为 0.36，这是我们迄今为止表现最好的信用评分模型。让我们看看它与最后一种集成方法的比较。

提示

可以通过在caret中指定method = "AdaBoost.M1"来调整 AdaBoost.M1 算法。

随机森林

另一种基于集成的方法，称为随机森林（或决策树森林），仅专注于决策树的集成。这一方法由 Leo Breiman 和 Adele Cutler 提倡，结合了袋装（bagging）和随机特征选择的基本原则，以增加决策树模型的多样性。在生成了决策树集成（森林）之后，模型通过投票来合并这些树的预测结果。

注意

关于随机森林构建的更多细节，请参阅 Breiman L. Random Forests。机器学习，2001；45：5-32。

随机森林将多样性和强大功能融合为单一的机器学习方法。由于集成只使用了完整特征集中的一小部分随机特征，随机森林可以处理非常大的数据集，在这些数据集中，所谓的“维度灾难”可能会导致其他模型失败。与此同时，它在大多数学习任务上的错误率与几乎所有其他方法相当。

提示

虽然“随机森林”（Random Forests）这个术语是由 Breiman 和 Cutler 注册的商标，但有时人们也用它来泛指任何类型的决策树集成。一个严格的学者会使用更通用的术语“决策树森林”，除非是在指 Breiman 和 Cutler 的具体实现。

值得注意的是，相对于其他基于集成的方法，随机森林非常具有竞争力，并且相较于其他方法具有关键优势。例如，随机森林往往更容易使用，且不容易发生过拟合。下表列出了随机森林模型的一般优缺点：

优点	缺点

|

• 一种在大多数问题上表现良好的通用模型

• 能处理噪声或缺失数据以及类别型或连续型特征

• 仅选择最重要的特征

• 可以用于具有极大量特征或样本的数据

|

• 与决策树不同，该模型不容易解释

• 可能需要一些工作来调整模型以适应数据

|

由于随机森林具有强大的功能、广泛的适应性和易用性，它们正在迅速成为最受欢迎的机器学习方法之一。在本章后面，我们将把随机森林模型与增强版的 C5.0 决策树进行正面比较。

训练随机森林

虽然在 R 中有多个包可以创建随机森林，但randomForest包可能是最忠实于 Breiman 和 Cutler 规范的实现，并且也受到caret包的支持，可以进行自动化调参。训练该模型的语法如下：

默认情况下，randomForest()函数会创建一个包含 500 棵树的集成模型，每棵树在每次分裂时都会考虑sqrt(p)个随机特征，其中p是训练数据集中的特征数量，sqrt()是 R 的平方根函数。是否使用这些默认参数取决于学习任务和训练数据的性质。一般来说，更复杂的学习问题和更大的数据集（无论是更多的特征还是更多的样本）都能通过更多的树来取得更好的效果，但这需要与训练更多树的计算开销进行平衡。

使用大量决策树的目的是训练足够多的树，以便每个特征都有机会出现在多个模型中。这也是mtry参数的默认值sqrt(p)的基础；使用该值限制特征的数量，确保树与树之间有足够的随机变化。例如，由于信用数据有 16 个特征，每棵树在任何时候只能在四个特征上进行分裂。

让我们来看一下默认的randomForest()参数如何与信用数据一起使用。我们将像之前训练其他学习器一样训练模型。再次强调，set.seed()函数确保结果可以复现：

> library(randomForest)
> set.seed(300)
> rf <- randomForest(default ~ ., data = credit)

要查看模型表现的总结，我们只需输入结果对象的名称：

> rf

Call:
 randomForest(formula = default ~ ., data = credit) 
 Type of random forest: classification
 Number of trees: 500
No. of variables tried at each split: 4

 OOB estimate of error rate: 23.8%
Confusion matrix:
 no yes class.error
no  640  60  0.08571429
yes 178 122  0.59333333

输出结果显示随机森林包括了 500 棵树，并且在每次分裂时尝试了四个变量，正如我们所期望的那样。乍一看，你可能会对困惑矩阵中看似不理想的表现感到惊讶——23.8%的错误率远高于到目前为止其他集成方法的重置误差。然而，这个困惑矩阵并没有显示重置误差。相反，它反映的是袋外错误率（在输出中列为OOB estimate of error rate），与重置误差不同，袋外错误率是测试集错误的无偏估计。这意味着它应该是对未来表现的相当合理的估计。

袋外估计在构建随机森林时计算。实际上，任何没有被选为单棵树的自助法样本的例子，都可以用来测试模型在未见过数据上的表现。在森林构建结束时，每个例子在每次被保留时的预测都会被统计，并通过投票来确定该例子的最终预测结果。这样的预测的总误差率即为袋外误差率。

评估随机森林性能

如前所述，randomForest()函数得到了caret的支持，这使我们能够优化模型，同时计算袋外误差率以外的性能度量。为了增加趣味性，我们将比较一个自动调优的随机森林与我们开发的最佳自动调优提升 C5.0 模型。我们将这次实验视为希望找到一个候选模型，提交到机器学习竞赛中。

我们首先需要加载caret并设置训练控制选项。为了对模型性能进行最准确的比较，我们将使用重复的 10 折交叉验证，或者说是重复 10 次的 10 折交叉验证。这意味着模型的构建时间将大大增加，并且评估计算量会更大，但由于这是我们的最终比较，我们必须非常确保我们做出了正确的选择；这场对决的胜者将是我们唯一进入机器学习竞赛的模型。

> library(caret)
> ctrl <- trainControl(method = "repeatedcv",
 number = 10, repeats = 10)

接下来，我们将为随机森林设置调优网格。该模型的唯一调优参数是mtry，它定义了每次分裂时随机选择多少特征。默认情况下，我们知道随机森林将使用sqrt(16)，即每棵树使用四个特征。为了全面起见，我们还将测试该值的一半、两倍以及完整的 16 个特征。因此，我们需要创建一个包含2、4、8和16的网格，如下所示：

> grid_rf <- expand.grid(.mtry = c(2, 4, 8, 16))

提示

一个在每次分裂时考虑完整特征集的随机森林，本质上与一个集成决策树模型相同。

我们可以将生成的网格传递给train()函数，并使用ctrl对象，如下所示。我们将使用 kappa 指标来选择最佳模型：

> set.seed(300)
> m_rf <- train(default ~ ., data = credit, method = "rf",
 metric = "Kappa", trControl = ctrl,
 tuneGrid = grid_rf)

前面的命令可能需要一些时间才能完成，因为它有很多工作要做！完成后，我们将与使用10、20、30和40次迭代的提升树进行比较：

> grid_c50 <- expand.grid(.model = "tree",
 .trials = c(10, 20, 30, 40),
 .winnow = "FALSE")
> set.seed(300)
> m_c50 <- train(default ~ ., data = credit, method = "C5.0",
 metric = "Kappa", trControl = ctrl,
 tuneGrid = grid_c50)

当 C5.0 决策树最终完成时，我们可以将这两种方法进行并排比较。对于随机森林模型，结果是：

> m_rf

Resampling results across tuning parameters:

 mtry  Accuracy  Kappa      Accuracy SD  Kappa SD 
 2    0.7247    0.1284142  0.01690466   0.06364740
 4    0.7499    0.2933332  0.02989865   0.08768815
 8    0.7539    0.3379986  0.03107160   0.08353988
 16    0.7556    0.3613151  0.03379439   0.08891300

对于提升 C5.0 模型，结果是：

> m_c50

Resampling results across tuning parameters:

 trials  Accuracy  Kappa      Accuracy SD  Kappa SD 
 10      0.7325    0.3215655  0.04021093   0.09519817
 20      0.7343    0.3268052  0.04033333   0.09711408
 30      0.7381    0.3343137  0.03672709   0.08942323
 40      0.7388    0.3335082  0.03934514   0.09746073

在mtry = 16时，随机森林模型的 kappa 约为 0.361，是这八个模型中表现最好的。它高于最佳的 C5.0 决策树，其 kappa 约为 0.334，并且略高于 kappa 约为 0.360 的AdaBoost.M1模型。根据这些结果，我们将提交随机森林模型作为最终模型。在没有实际在竞赛数据上评估模型的情况下，我们无法确定它是否最终获胜，但根据我们的性能评估，它是更安全的选择。幸运的话，也许我们能赢得奖项。

总结

阅读本章后，你现在应该知道了在数据挖掘和机器学习竞赛中获胜所需的基本技术。自动调优方法可以帮助从单一模型中挤出每一分性能。另一方面，通过创建多个协同工作的机器学习模型，也可以实现性能的提升。

尽管本章旨在帮助你准备竞赛级别的模型，但请注意，你的竞争对手也可以使用相同的技术。你不能停滞不前；因此，继续将独特的方法加入你的工具箱。也许你能带来独特的专业知识，或者也许你的优势在于数据准备时对细节的关注。无论如何，实践出真知，因此利用开放竞赛来测试、评估和提升你自己的机器学习技能。

在下一章——本书的最后一章——我们将从更高的角度审视如何利用 R 语言将机器学习应用于一些高度专业化且复杂的领域。你将获得将机器学习应用于前沿任务所需的知识。