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

R 机器学习第四版（二）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：概率学习 - 使用朴素贝叶斯进行分类

当气象学家提供天气预报时，降水通常用“70%降雨可能性”这样的短语来描述。这类预报被称为降水概率报告。你有没有考虑过它们是如何计算的？这是一个令人困惑的问题，因为在现实中，要么下雨，要么不下雨，这是绝对确定的。

天气估计基于概率方法，这些方法涉及描述不确定性。它们使用过去事件的数据来预测未来事件。在天气的情况下，降雨的可能性描述了在相似大气条件下发生降水的先前天数所占的比例。70%的降雨可能性意味着在 10 个过去类似条件下，有 7 个地方发生了降水。

本章介绍了朴素贝叶斯算法，它使用概率的方式与天气预报非常相似。在研究这种方法时，你将了解：

• 概率的基本原理

• 使用 R 分析文本数据所需的专业方法和数据结构

• 如何使用朴素贝叶斯构建短信服务（SMS）垃圾信息过滤器

如果你之前上过统计学课程，本章的一些材料可能对你来说是复习。即便如此，刷新你对概率的了解可能也有帮助。你会发现这些原则是朴素贝叶斯获得这样一个奇怪名称的基础。

理解朴素贝叶斯

理解朴素贝叶斯算法所需的基本统计思想已经存在了几个世纪。这项技术源于 18 世纪数学家托马斯·贝叶斯的工作，他开发了描述事件概率及其在额外信息的基础上如何修订的基础原则。这些原则构成了现在被称为贝叶斯方法的基础。

我们将在稍后更详细地介绍这些方法。现在，只需说一个概率是一个介于零和一之间的数字（或从 0 到 100%）即可，它捕捉了在现有证据的基础上事件发生的可能性。概率越低，事件发生的可能性越小。零概率表示事件肯定不会发生，而一概率表示事件将以绝对确定性发生。生活中最有趣的事件往往具有不确定的概率；估计它们发生的可能性有助于我们通过揭示最可能的结果来做出更好的决策。

基于贝叶斯方法的分类器利用训练数据根据特征值提供的证据计算每个结果的概率。当分类器后来应用于未标记的数据时，它使用这些计算出的概率来预测新示例最可能的类别。这是一个简单的想法，但结果是一个可以与更复杂算法相媲美的方法。事实上，贝叶斯分类器已被用于：

• 文本分类，如垃圾邮件（垃圾邮件过滤）

• 计算机网络中的入侵或异常检测

• 根据一组观察到的症状诊断医疗状况

通常，贝叶斯分类器最适合应用于需要同时考虑多个属性信息以估计结果整体概率的问题。虽然许多机器学习算法忽略了具有较弱影响特征，但贝叶斯方法利用所有可用证据微妙地改变预测。这意味着即使大部分特征的影响相对较小，但在贝叶斯模型中它们的综合影响可能相当大。

贝叶斯方法的基本概念

在深入研究朴素贝叶斯算法之前，花些时间定义贝叶斯方法中使用的概念是值得的。用一句话总结，贝叶斯概率理论根植于这样一个观点：估计一个 事件 或潜在结果的似然性应该基于多个 试验 或事件发生机会的证据。

下表展示了几个现实世界结果的事件和试验：

事件	试验
正面朝上	抛硬币
雨天	单日（或另一个时间段）
消息是垃圾邮件	一封 incoming 电子邮件
候选人成为总统	总统选举
死亡率	医院病人
中奖	一张彩票

贝叶斯方法提供了从观察数据中估计这些事件概率的见解。为了了解这一点，我们需要形式化我们对概率的理解。

理解概率

通过将事件发生的试验次数除以总试验次数来估计事件的概率。例如，如果今天有类似条件的 10 天中有 3 天下雨，那么今天下雨的概率可以估计为 3 / 10 = 0.30 或 30%。同样，如果 50 封之前的电子邮件中有 10 封是垃圾邮件，那么任何新收到的邮件是垃圾邮件的概率可以估计为 10 / 50 = 0.20 或 20%。

为了表示这些概率，我们使用形式为 P(A) 的符号，它表示事件 A 的概率。例如，P(rain) = 0.30 表示有 30% 的降雨概率，或 P(spam) = 0.20 描述一个新收到的消息有 20% 的概率是垃圾邮件。

由于试验总是导致某些结果发生，因此试验所有可能结果的概率总和必须始终为 1。因此，如果试验恰好有两个结果且这些结果不能同时发生，那么知道任意一个结果发生的概率就可以揭示另一个结果发生的概率。这种情况适用于许多结果，例如硬币的正反面，或垃圾邮件与合法电子邮件（也称为“ham”），使用这个原理，知道P(spam) = 0.20可以让我们计算出P(ham) = 1 – 0.20 = 0.80。这仅适用于垃圾邮件和 ham 是互斥且穷尽的事件，这意味着它们不能同时发生，并且是唯一的可能结果。

单个事件不能同时发生和未发生。这意味着事件总是与其补集互斥且穷尽，或者包含所有其他结果的补集，其中感兴趣的事件未发生。事件A的补集通常表示为A^c 或A’。

此外，可以使用简写符号P(A^c*)或P(¬A)来表示事件A不发生的概率。例如，符号P(¬spam) = 0.80*表示消息不是垃圾邮件的概率为 80%。

为了说明事件及其补集，想象一个二维空间，该空间被划分为每个事件的概率，通常是有帮助的。在以下图中，矩形代表电子邮件消息的可能结果。圆圈代表消息是垃圾邮件的 20%概率。剩余的 80%代表补集P(¬spam)，或不是垃圾邮件的消息：

图 4.1：所有电子邮件的概率空间可以表示为垃圾邮件和正常邮件的分区

理解联合概率

通常，我们在同一试验中会对几个非互斥事件进行监控。如果某些事件与感兴趣的事件同时发生，我们可能能够利用它们进行预测。例如，考虑一个基于电子邮件消息包含单词Viagra的结果的第二个事件。更新此第二个事件的先前列表可能如下所示：

图 4.2：非互斥事件表示为重叠的分区

注意在图中，Viagra 圆圈与图中的垃圾邮件和 ham 区域重叠，并且垃圾邮件圆圈包括 Viagra 圆圈未覆盖的区域。这表明并非所有垃圾邮件都包含 Viagra 这个词，并且一些包含 Viagra 的消息是 ham。然而，由于这个词在垃圾邮件之外出现得非常少，它在新的传入消息中的出现将是该消息是垃圾邮件的强烈证据。

为了更仔细地观察这两个圆之间的重叠，我们将使用一种称为维恩图的可视化方法。这种图最早在 19 世纪末由数学家约翰·文恩使用，它使用圆来表示项目集合的重叠。与大多数维恩图一样，图中圆的大小和重叠程度没有意义。相反，它被用作提醒，将概率分配给所有事件组合。垃圾邮件和 Viagra 的维恩图可能如下所示：

图 4.3：一个维恩图展示了垃圾邮件和 Viagra 事件的交集

我们知道所有消息中有 20%是垃圾邮件（左边的圆），所有消息中有 5%包含Viagra这个词（右边的圆）。我们希望量化这两个比例之间的重叠程度。换句话说，我们希望估计P(spam)和P(Viagra)同时发生的概率，这可以写成P(spam Viagra)。符号表示两个事件的交集；A B的表示法指的是A和B同时发生的事件。

计算*P(spam Viagra)*取决于两个事件的联合概率，即一个事件的概率如何与另一个事件的概率相关。如果两个事件完全无关，它们被称为独立事件。这并不是说独立事件不能同时发生；事件独立性仅仅意味着知道一个事件的结局不会提供任何关于另一个事件结局的信息。例如，抛硬币得到正面结果的结局与某一天是雨天还是晴天无关。

如果所有事件都是独立的，那么通过观察另一个事件来预测一个事件将是不可能的。换句话说，相关事件是预测建模的基础。就像云的存在预示着雨天一样，Viagra这个词的出现预示着垃圾邮件。

图 4.4：机器学习如何识别有用模式需要相关事件

计算相关事件的概率比独立事件要复杂一些。如果P(spam)和P(Viagra)是独立的，我们可以轻松地计算出P(spam Viagra)，即两个事件同时发生的概率。因为所有消息中有 20%是垃圾邮件，所有邮件中有 5%包含Viagra这个词，我们可以假设所有包含Viagra这个词的消息中有 1%是垃圾邮件。这是因为0.05 * 0.20 = 0.01。更普遍地，对于独立事件A和B，两个事件同时发生的概率可以计算为P(A B) = P(A) * P(B)。

话虽如此，我们知道*P(spam)和P(Viagra)*很可能高度相关，这意味着这个计算是不正确的。为了得到一个更合理的估计，我们需要使用这两个事件之间关系的更谨慎的公式，这个公式基于更先进的贝叶斯方法。

使用贝叶斯定理计算条件概率

可以使用贝叶斯定理来描述相关事件之间的关系，它提供了一种思考如何根据另一个事件提供的证据来修订一个事件概率估计的方法。一种公式如下：

符号P(A|B)读作事件B发生的情况下事件A的概率。这被称为条件概率，因为A的概率依赖于（即条件于）事件B的发生。

贝叶斯定理表明，P(A|B)的最佳估计是在所有发生事件B的试验中，事件A发生的试验比例。这意味着如果每次观察到B时A和B经常一起发生，事件A的概率就会更高。请注意，这个公式调整了P(A B)以反映B发生的概率。如果B非常罕见，P(B)和P(A B)将始终很小；然而，如果A几乎总是与B一起发生，尽管B很罕见，*P(A|B)*仍然会很高。

根据定义，P(A B) = P(A|B) * P(B)，这是一个可以通过对先前公式应用一点代数轻松推导出的事实。利用*P(A B) = P(B A)*的知识重新排列这个公式，我们得出结论，P(A B) = P(B|A) * P(A)，我们可以在贝叶斯定理的以下公式中使用它：

事实上，这是基于我们将它应用于机器学习时将变得清晰的原因的传统贝叶斯定理公式。首先，为了更好地理解贝叶斯定理在实际中的工作原理，让我们回顾一下我们的假设性垃圾邮件过滤器。

在不知道收到的邮件内容的情况下，对其是否为垃圾邮件的最佳估计将是P(spam)，即任何先前邮件是垃圾邮件的概率。这个估计被称为先验概率。我们之前发现这个概率是 20%。

假设你通过更仔细地查看先前收到的邮件集并检查“Viagra”一词出现的频率获得了额外的证据。单词Viagra在先前垃圾邮件中被使用的概率，或P(Viagra|spam)，被称为似然性。Viagra在任何邮件中出现的概率，或P(Viagra)，被称为边缘似然性。

通过将贝叶斯定理应用于这一证据，我们可以计算一个后验概率，该概率衡量一条消息是垃圾邮件的可能性。如果后验概率大于 50%，则该消息更有可能是垃圾邮件而不是正常邮件，可能需要过滤。以下公式显示了贝叶斯定理是如何应用于先前电子邮件消息提供的证据的：

图 4.5：贝叶斯定理作用于先前收到的电子邮件

要计算贝叶斯定理的组成部分，构建一个频率表（如下表中左侧所示）记录Viagra在垃圾邮件和正常邮件中出现的次数是有帮助的。就像一个双向交叉表一样，表的其中一个维度表示类别变量（垃圾邮件或正常邮件）的水平，而另一个维度表示特征（Viagra：是或否）的水平。然后，单元格表示具有指定类别值和特征值的实例数量。

然后，可以使用频率表来构建一个似然表，如下表中右侧所示。似然表的行表示在电子邮件是垃圾邮件或正常邮件的情况下，对于Viagra（是/否）的条件概率。

图 4.6：频率和似然表是计算垃圾邮件后验概率的基础

似然表显示P(Viagra=Yes|spam) = 4 / 20 = 0.20，这表明如果一条消息是垃圾邮件，那么包含术语Viagra的概率是 20%。此外，由于P(A B) = P(B|A) * P(A)，我们可以计算P(spam Viagra)为P(Viagra|spam) * P(spam) = (4 / 20) * (20 / 100) = 0.04。同样的结果可以在频率表中找到，该表指出 100 条消息中有 4 条是垃圾邮件并包含术语Viagra。无论如何，这比我们之前在错误假设独立性下计算的*P(A B) = P(A) * P(B)*的估计值 0.01 高出四倍。这当然说明了贝叶斯定理在估计联合概率中的重要性。

要计算后验概率，P(spam|Viagra)，我们只需取P(Viagra|spam) * P(spam) / P(Viagra)，或者*(4 / 20) * (20 / 100) / (5 / 100) = 0.80*。因此，如果一条消息包含单词Viagra，那么这条消息是垃圾邮件的概率是 80%。鉴于这一发现，任何包含此术语的消息可能应该被过滤。

这正是商业垃圾邮件过滤器的工作方式，尽管在计算频率和似然表时，它们会同时考虑更多的单词。在下一节中，我们将看到如何将这种方法适应涉及额外特征的情况。

简单贝叶斯算法

朴素贝叶斯算法定义了一种简单的方法，将贝叶斯定理应用于分类问题。尽管它不是唯一利用贝叶斯方法的机器学习方法，但它是最常见的。由于它在文本分类中的成功，朴素贝叶斯变得非常流行，一度成为事实上的标准。该算法的优点和缺点如下：

优点	缺点

|

• 简单、快速且非常有效

• 在有噪声和缺失数据以及大量特征的情况下表现良好

• 需要相对较少的训练示例

• 容易获得预测的估计概率

|

• 依赖于一个经常是错误的假设，即特征同等重要且相互独立

• 不适用于具有许多数值特征的集合

• 估计的概率不如预测的类别可靠

|

朴素贝叶斯算法之所以被称为“朴素”，是因为它对数据做出了一些所谓的“朴素”假设。特别是，朴素贝叶斯假设数据集中的所有特征都是同等重要且相互独立的。在大多数实际应用中，这些假设很少是真实的。

例如，当尝试通过监控电子邮件消息来识别垃圾邮件时，几乎可以肯定的是，某些特征将比其他特征更重要。例如，电子邮件发送者可能是比消息文本更重要的垃圾邮件指示器。此外，消息正文中的单词并不是相互独立的，因为某些单词的出现是其他单词也很可能出现的很好指示。包含“Viagra”一词的消息很可能也包含“prescription”或“drugs”一词。

然而，在大多数情况下，即使这些假设被违反，朴素贝叶斯仍然表现出惊人的性能。即使在特征之间存在强依赖性的情况下，也是如此。由于该算法在各种条件下的灵活性和准确性，尤其是在较小的训练数据集上，朴素贝叶斯经常是分类学习任务的合理基线候选者。

尽管朴素贝叶斯存在错误的假设，但它为何表现良好的确切原因一直是许多猜测的对象。一种解释是，只要预测准确，获得概率的精确估计并不重要。例如，如果一个垃圾邮件过滤器正确地识别出垃圾邮件，那么预测垃圾邮件的概率是 51%还是 99%又有什么关系呢？关于这个话题的一个讨论，请参阅*《在零一损失下简单贝叶斯分类器的最优性》，作者：Domingos, P. 和 Pazzani, M.，机器学习，1997 年，第 29 卷，第 103-130 页*。

使用朴素贝叶斯进行分类

让我们通过添加一些额外的监控术语来扩展我们的垃圾邮件过滤器，除了术语Viagra之外，还包括money、groceries和unsubscribe。朴素贝叶斯学习器通过构建这四个词（标记为W¹、W²、W³和W⁴）出现的可能性表来训练，如下图中 100 封电子邮件所示：

图 4.7：扩展的表格增加了垃圾邮件和正常邮件中额外术语的可能性

当收到新消息时，我们需要计算后验概率，以确定它们更有可能是垃圾邮件还是正常邮件，给定在消息文本中找到这些词的可能性。例如，假设一条消息包含术语Viagra和unsubscribe，但不包含money或groceries。

使用贝叶斯定理，我们可以将问题定义为以下公式。这计算了在Viagra = 是、Money = 否、Groceries = 否和Unsubscribe = 是的条件下，一条消息是垃圾邮件的概率：

由于两个原因，这个公式在计算上很难解决。首先，随着附加特征的添加，需要大量的内存来存储所有可能交点事件的概率。想象一下四个词事件 Venn 图的复杂性，更不用说数百个或更多。其次，许多这些潜在的交点在过去的资料中从未被观察到，这会导致联合概率为零，并导致后面会变得明显的问题。

如果我们利用朴素贝叶斯对事件之间独立性的朴素假设，计算将变得更加合理。具体来说，它假设类条件独立性，这意味着只要事件基于相同的类值，它们就是独立的。条件独立性假设允许我们使用独立事件的概率规则，该规则指出 P(A B) = P(A) * P(B)。这通过允许我们乘以单个条件概率而不是计算复杂的条件联合概率来简化分子。

最后，因为分母不依赖于目标类（垃圾邮件或正常邮件），它被视为一个常数，暂时可以忽略。这意味着垃圾邮件的条件概率可以表示为：

并且可以表示消息是正常邮件的概率为：

注意，等号已被比例符号（类似于侧向的、开口的“8”）替换，以表明分母已被省略。

使用可能性表中的值，我们可以开始填写这些方程中的数字。然后，垃圾邮件的整体可能性如下：

(4 / 20) * (10 / 20) * (20 / 20) * (12 / 20) * (20 / 100) = 0.012

而正常邮件的可能性是：

(1 / 80) * (66 / 80) * (71 / 80) * (23 / 80) * (80 / 100) = 0.002

因为0.012 / 0.002 = 6，我们可以说这条消息有 6 倍的可能性是垃圾邮件，而不是 ham。然而，为了将这些数字转换为概率，我们需要最后一步重新引入之前排除的除数。本质上，我们必须通过除以所有可能结果的总似然值来重新缩放每个结果的似然值。

这样，垃圾邮件的概率等于消息是垃圾邮件的可能性除以消息是垃圾邮件或 ham 的可能性：

0.012 / (0.012 + 0.002) = 0.857

同样，ham 的概率等于消息是 ham 的可能性除以消息是垃圾邮件或 ham 的可能性：

0.002 / (0.012 + 0.002) = 0.143

根据这条消息中找到的单词模式，我们预计这条消息有 85.7%的概率是垃圾邮件，有 14.3%的概率是 ham。因为这些是互斥且穷尽的概率事件，所以它们的概率总和为 1。

在前一个例子中使用的朴素贝叶斯分类算法可以用以下公式总结。给定特征[F][1]到[F][n]提供的证据，对于类别C的L级概率等于每个证据在类别级条件下的概率乘积，类别级的先验概率，以及一个缩放因子1 / Z，它将似然值转换为概率。这被公式化为：

尽管这个方程看起来令人畏惧，但正如垃圾邮件过滤示例所示，这一系列步骤相当简单。首先构建一个频率表，然后使用这个表构建一个似然表，并使用朴素假设的独立性乘出条件概率。最后，除以总似然值，将每个类别的似然值转换为概率。尝试手动计算几次后，这将成为第二本能。

拉普拉斯估计器

在我们使用朴素贝叶斯解决更复杂的问题之前，有一些细微之处需要考虑。假设我们收到了另一条消息，这次包含所有四个术语：Viagra，groceries，money，和unsubscribe。使用之前的方法，我们可以计算垃圾邮件的可能性为：

(4 / 20) * (10 / 20) * (0 / 20) * (12 / 20) * (20 / 100) = 0

并且 ham 的似然值为：

(1 / 80) * (14 / 80) * (8 / 80) * (23 / 80) * (80 / 100) = 0.00005

因此，垃圾邮件的概率为：

0 / (0 + 0.00005) = 0

并且 ham 的概率为：

0.00005 / (0 + 0.00005) = 1

这些结果表明，该消息有 0%的概率是垃圾邮件，有 100%的概率是正常邮件。这个预测有道理吗？可能没有。消息中包含了一些通常与垃圾邮件相关的词汇，包括伟哥，这在合法消息中很少使用。因此，该消息被错误分类的可能性非常大。

如果一个事件在类的某个或多个级别上从未发生，那么由此产生的似然值将是零。例如，术语杂货以前从未出现在垃圾邮件中。因此，P(groceries|spam) = 0%。

现在，由于朴素贝叶斯公式中的概率是链式相乘的，这个零值导致垃圾邮件的后验概率为零，使得单词杂货能够有效地抵消并推翻所有其他证据。即使电子邮件在其他方面几乎肯定会被认为是垃圾邮件，垃圾邮件中缺少单词杂货也会始终否决其他证据，导致垃圾邮件的概率为零。

解决这个问题的方法涉及使用一种称为拉普拉斯估计器的东西，该估计器是以法国数学家皮埃尔-西蒙·拉普拉斯的名字命名的。拉普拉斯估计器将一个小数加到频率表中的每个计数上，这确保了每个特征在每种类别中都有非零发生的概率。通常，拉普拉斯估计器设置为 1，这确保每个类别-特征组合至少在数据中出现一次。

拉普拉斯估计器可以设置为任何值，并且不一定需要对每个特征都相同。如果你是一个虔诚的贝叶斯主义者，你可以使用拉普拉斯估计器来反映一个特征与类别相关的假设先验概率。在实践中，给定足够大的训练数据集，这通常是过度的。因此，值 1 几乎总是被使用。

让我们看看这如何影响我们对这条消息的预测。使用拉普拉斯值为 1，我们在似然函数的每个分子上加上 1。然后，我们需要在每个条件概率的分母上加上 4，以补偿分子上增加的 4 个额外值。因此，垃圾邮件的可能性是：

(5 / 24) * (11 / 24) * (1 / 24) * (13 / 24) * (20 / 100) = 0.0004

正常邮件的可能性是：

(2 / 84) * (15 / 84) * (9 / 84) * (24 / 84) * (80 / 100) = 0.0001

通过计算0.0004 / (0.0004 + 0.0001)，我们发现垃圾邮件的概率是 80%，因此正常邮件的概率大约是 20%。这个结果比仅用术语杂货确定结果时计算的P(spam) = 0更合理。

尽管拉普拉斯估计量被添加到了似然函数的分子和分母中，但它并没有被添加到先验概率中——即 20/100 和 80/100 的值。这是因为，根据数据中观察到的结果，我们对于垃圾邮件和正常邮件的整体概率的最佳估计仍然是 20%和 80%。

使用数值特征与简单贝叶斯

简单贝叶斯使用频率表来学习数据，这意味着每个特征必须是分类的，以便创建由类和特征值组成的矩阵的组合。由于数值特征没有值类别，因此前面的算法不能直接处理数值数据。然而，有几种方法可以解决这个问题。

一个简单而有效的方法是将数值特征离散化，这仅仅意味着将数字放入称为分组的类别中。因此，离散化有时也被称为分组。这种方法在有大量训练数据时效果最佳。

有几种不同的方法可以对数值特征进行离散化。可能最常见的方法是探索数据中的自然类别或分布中的切割点。例如，假设你向垃圾邮件数据集中添加了一个特征，记录了邮件发送的白天或夜晚时间，从午夜过后的 0 到 24 小时。使用直方图表示，时间数据可能看起来像以下图表：

图 4.8：可视化接收电子邮件时间分布的直方图

在清晨时分，消息频率较低。在办公时间活动增加，而在晚上逐渐减少。这形成了四个自然的活动分组，如图中虚线所示。这些表示可以将数值数据划分为不同级别以创建新的分类特征的地方，然后可以使用简单贝叶斯。

四个分组的选取是基于数据的自然分布以及对于一整天中垃圾邮件比例可能变化的直觉。我们可能预期垃圾邮件发送者会在深夜活动，或者他们可能在白天人们检查邮件时活动。尽管如此，为了捕捉这些趋势，我们同样可以使用三个或十二个分组。

如果没有明显的切割点，一个选项是使用在第二章中介绍的量数对特征进行离散化。你可以用三分位数将数据分为三个分组，用四分位数分为四个分组，或者用五分位数分为五个分组。

需要注意的是，将数值特征离散化总是会导致信息量的减少，因为特征的原有粒度被减少到更少的类别。重要的是要找到平衡点。分类太少可能导致重要趋势被掩盖。分类太多可能导致朴素贝叶斯频率表中的计数很小，这可能会增加算法对噪声数据的敏感性。

示例 – 使用朴素贝叶斯算法过滤手机垃圾邮件

随着全球手机使用的增长，不道德的营销人员开辟了一条新的电子垃圾邮件途径。这些广告商利用短信文本消息针对潜在消费者发送不受欢迎的广告，即短信垃圾邮件。这种垃圾邮件很麻烦，因为与电子邮件垃圾邮件不同，短信消息由于手机无处不在，特别具有破坏性。开发一个能够过滤短信垃圾邮件的分类算法将为手机服务提供商提供一个有用的工具。

由于朴素贝叶斯分类器在电子邮件垃圾邮件过滤中已被成功应用，因此它似乎也可以应用于短信垃圾邮件。然而，相对于电子邮件垃圾邮件，短信垃圾邮件对自动过滤器提出了额外的挑战。短信消息通常限制在 160 个字符以内，这减少了用于识别消息是否为垃圾信息的文本量。这种限制，加上小型手机键盘，导致许多人采用一种短信缩写语，这进一步模糊了合法消息和垃圾邮件之间的界限。让我们看看简单的朴素贝叶斯分类器如何应对这些挑战。

第 1 步 – 收集数据

为了开发朴素贝叶斯分类器，我们将使用从www.dt.fee.unicamp.br/~tiago/smsspamcollection/的短信垃圾邮件收集中改编的数据。

要了解更多关于短信垃圾邮件收集如何发展的信息，请参阅关于新的短信垃圾邮件收集的有效性，Gómez, J. M.，Almeida, T. A.，和 Yamakami, A.，2012 年第 11 届 IEEE 国际机器学习与应用会议论文集。

此数据集包括短信文本，以及一个标签，表示消息是否不受欢迎。垃圾消息被标记为垃圾邮件，而合法消息被标记为 ham。以下表格显示了垃圾邮件和 ham 的一些示例：

样本短信垃圾信息	样本短信垃圾邮件

|

• 更好。周五我补上了，昨天吃得像猪一样。现在感觉有点糟糕。但至少不是那种让人难以忍受的疼痛。

• 如果他开始寻找，几天内就能找到工作。他有很大的潜力和才能。

• 我又找到一份工作了！是在医院做数据分析之类的，周一开始！不确定我的论文什么时候能完成。

|

• 恭喜你，你获得了 500 元的 CD 代金券或 125 元的礼品保证，并且免费参加每周一次的抽奖活动，发送短信 MUSIC 到 87066。

• 12 月特惠！你的手机使用 11 个月以上了吗？你有权免费升级到最新的彩色摄像头手机！免费拨打 08002986906 至移动更新公司。

• 情人节特别优惠！在我们的问答比赛中赢得 1000 英镑，并带你的伴侣去一次终身难忘的旅行！现在发送 GO 到 83600。每条信息 150 便士。

|

观察前面的消息，你是否注意到垃圾邮件的任何显著特征？一个值得注意的特征是，三个垃圾邮件中有两个使用了单词免费，而这个词并没有出现在任何正常邮件中。另一方面，两条正常邮件提到了具体的星期几，而垃圾邮件中则没有。

我们的朴素贝叶斯分类器将利用单词频率中的这些模式来确定短信消息似乎更适合垃圾邮件还是正常邮件。虽然不能排除单词免费出现在垃圾邮件之外的情况，但合法的消息更有可能提供额外的单词来提供上下文。

例如，一条正常邮件可能会问：“你星期天有空吗？”而一条垃圾邮件可能会使用短语“免费铃声”。分类器将根据消息中所有单词提供的证据计算垃圾邮件和正常邮件的概率。

第 2 步 - 探索和准备数据

构建我们的分类器第一步涉及对原始数据进行处理以进行分析。文本数据准备起来具有挑战性，因为需要将单词和句子转换成计算机可以理解的形式。我们将把我们的短信数据转换成一种称为词袋的表示形式，它提供了一种二元特征，表示每个单词是否出现在给定的示例中，同时忽略单词顺序或单词出现的上下文。尽管这是一种相对简单的表示形式，但正如我们很快将展示的，它对于许多分类任务来说已经足够好了。

这里使用的数据集已经从原始数据集稍作修改，以便更容易在 R 中使用。如果你打算跟随示例，请从 Packt 网站下载sms_spam.csv文件并将其保存到你的 R 工作目录中。

我们将首先导入 CSV 数据并将其保存到数据框中：

> sms_raw <- read.csv("sms_spam.csv") 

使用str()函数，我们可以看到sms_raw数据框包括 5559 条总短信消息，具有两个特征：type和text。短信类型已被编码为 ham 或 spam。text元素存储完整的原始短信消息文本：

> str(sms_raw) 

'data.frame':    5559 obs. of  2 variables:
 $ type: chr  "ham" "ham" "ham" "spam" ...
 $ text: chr  "Hope you are having a good week. Just checking in"
"K..give back my thanks." "Am also doing in cbe only. But have to
pay." "complimentary 4 STAR Ibiza Holiday or £10,000 cash needs your
URGENT collection. 09066364349 NOW from Landline not to lose out"|
__truncated__ ... 

type元素目前是一个字符向量。由于这是一个分类变量，最好将其转换为因子，如下面的代码所示：

> sms_raw$type <- factor(sms_raw$type) 

使用str()和table()函数检查，我们发现type现在已被适当地重新编码为因子。此外，我们还可以看到在我们的数据中，有 747 条（大约 13%）短信消息被标记为spam，其余的则被标记为ham：

> str(sms_raw$type) 

 Factor w/ 2 levels "ham","spam": 1 1 1 2 2 1 1 1 2 1 ... 

> table(sms_raw$type) 

 ham spam
4812 747 

目前，我们将保持短信文本不变。您将在下一节中了解到，处理原始短信需要使用一套专门为处理文本数据设计的强大工具。

数据准备——清洗和标准化文本数据

短信消息是由单词、空格、数字和标点符号组成的文本字符串。处理这种复杂的数据类型需要大量的思考和努力。需要考虑如何删除数字和标点；处理无趣的词，例如和、但和或；以及将句子分解成单个单词。幸运的是，这种功能已经由 R 社区的成员在一个名为tm的文本挖掘包中提供。

tm包最初是由维也纳经济和商业大学的 Ingo Feinerer 作为学位论文项目创建的。要了解更多信息，请参阅*《R 中的文本挖掘基础设施》，作者：Feinerer, I.，Hornik, K.和 Meyer, D.，统计软件杂志，2008 年，第 25 卷，第 1-54 页*。

可以通过install.packages("tm")命令安装tm包，并通过library(tm)命令加载。即使您已经安装了它，重新安装也可能值得，以确保您的版本是最新的，因为tm包仍在积极开发中。这偶尔会导致其功能发生变化。

本章使用tm版本 0.7-11 进行了测试，截至 2023 年 5 月，这是当前版本。如果您看到输出有差异或代码无法工作，您可能使用的是不同版本。本书的 Packt 支持页面以及其 GitHub 仓库将发布针对未来tm包版本的解决方案，如果发现重大变化。

处理文本数据的第一个步骤是创建一个语料库，它是一组文本文档的集合。文档可以是短或长的，从单个新闻文章、书籍的页面、网页页面，甚至整本书。在我们的情况下，语料库将是一组短信。

要创建一个语料库，我们将使用tm包中的VCorpus()函数，它指的是一个易失性语料库——这里的“易失性”意味着它存储在内存中，而不是存储在磁盘上（使用PCorpus()函数来访问存储在数据库中的永久性语料库）。此函数要求我们指定语料库文档的来源，这可能是计算机的文件系统、数据库、网络或其他地方。

由于我们已经在 R 中加载了短信文本，我们将使用VectorSource()读取器函数从现有的sms_raw$text向量创建一个源对象，然后可以将其如下提供给VCorpus()：

> sms_corpus <- VCorpus(VectorSource(sms_raw$text)) 

生成的语料库对象以sms_corpus命名保存。

通过指定可选的readerControl参数，VCorpus()函数可以用于从 PDF 和 Microsoft Word 文件等来源导入文本。要了解更多信息，请使用vignette("tm")命令检查tm包的数据导入部分。

通过打印语料库，我们看到它包含训练数据中 5,559 条短信消息的文档：

> print(sms_corpus) 

<<VCorpus>>
Metadata:  corpus specific: 0, document level (indexed): 0
Content:  documents: 5559 

现在，因为tm语料库本质上是一个复杂列表，我们可以使用列表操作来选择语料库中的文档。inspect()函数显示了结果的摘要。例如，以下命令将查看语料库中第一条和第二条短信消息的摘要：

> inspect(sms_corpus[1:2]) 

<<VCorpus>>
Metadata:  corpus specific: 0, document level (indexed): 0
Content:  documents: 2
[[1]]
<<PlainTextDocument>>
Metadata:  7
Content:  chars: 49
[[2]]
<<PlainTextDocument>>
Metadata:  7
Content:  chars: 23 

要查看实际的消息文本，必须对所需的消息应用as.character()函数。要查看一条消息，请在单个列表元素上使用as.character()函数，注意需要使用双括号表示法：

> as.character(sms_corpus[[1]]) 

[1] "Hope you are having a good week. Just checking in" 

要查看多个文档，我们需要将as.character()应用于sms_corpus对象中的多个项。为此，我们将使用lapply()函数，它是 R 数据结构中应用程序的一个函数家族的一部分。这些函数，包括apply()和sapply()等，是 R 语言的关键习惯用法之一。经验丰富的 R 程序员使用这些函数的方式类似于在其他编程语言中使用for或while循环，因为它们会产生更易读（有时更高效）的代码。将as.character()应用于语料库元素子集的lapply()函数如下：

> lapply(sms_corpus[1:2], as.character) 

$'1'
[1] "Hope you are having a good week. Just checking in"
$'2'
[1] "K..give back my thanks." 

如前所述，语料库包含 5,559 条文本消息的原始文本。为了进行我们的分析，我们需要将这些消息划分为单个单词。首先，我们需要清理文本以标准化单词并删除会干扰结果的标点符号。例如，我们希望字符串Hello!、HELLO和hello都被计为同一单词的实例。

tm_map()函数提供了一个将转换（也称为映射）应用于tm语料库的方法。我们将使用此函数通过一系列转换清理我们的语料库，并将结果保存在一个名为corpus_clean的新对象中。

我们的第一步转换将标准化消息以仅使用小写字母。为此，R 提供了一个tolower()函数，该函数返回文本字符串的小写版本。为了将此函数应用于语料库，我们需要使用tm包装函数content_transformer()将tolower()视为一个转换函数，该函数可以用于访问语料库。完整的命令如下：

> sms_corpus_clean <- tm_map(sms_corpus,
    content_transformer(tolower)) 

为了检查命令是否按预期工作，让我们检查原始语料库中的第一条消息并将其与转换语料库中的相同消息进行比较：

> as.character(sms_corpus[[1]]) 

[1] "Hope you are having a good week. Just checking in" 

> as.character(sms_corpus_clean[[1]]) 

[1] "hope you are having a good week. just checking in" 

如预期的那样，清理后的语料库中的大写字母已被替换为小写版本。

content_transformer() 函数可用于应用更复杂的文本处理和清理过程，例如 grep 模式匹配和替换。只需编写一个自定义函数，并在应用 tm_map() 函数之前将其包装起来。

让我们继续清理过程，通过从短信中移除数字。尽管一些数字可能提供有用信息，但大多数数字可能仅对个别发送者独特，因此不会在所有消息中提供有用的模式。考虑到这一点，我们将从语料库中移除所有数字，如下所示：

> sms_corpus_clean <- tm_map(sms_corpus_clean, removeNumbers) 

注意，前面的代码没有使用 content_transformer() 函数。这是因为 removeNumbers() 与 tm 包一起提供，还包括其他不需要包装的映射函数。要查看其他内置转换，请简单地输入 getTransformations()。

我们接下来的任务是移除短信中的填充词，如 to、and、but 和 or。这些术语被称为 停用词，通常在文本挖掘之前被移除。这是因为尽管它们出现频率很高，但它们对我们模型的有用信息不多，因为它们不太可能区分垃圾邮件和正常邮件。

我们不会自己定义停用词列表，而是使用 tm 包提供的 stopwords() 函数。此函数允许我们访问来自各种语言的停用词集合。默认情况下，使用的是常见的英语停用词。要在 R 命令提示符中查看默认列表，请输入 stopwords()。要查看其他语言和选项，请输入 ?stopwords 以获取文档页面。

即使在单一语言中，也没有一个单一的、确定的停用词列表。例如，tm 中的默认英语列表包含大约 174 个单词，而另一个选项包含 571 个单词。你甚至可以指定自己的停用词列表。无论你选择哪个列表，都要记住这个转换的目标，即消除无用数据，同时尽可能保留有用信息。

仅定义停用词本身不是一种转换。我们需要的是一种方法来移除任何出现在停用词列表中的单词。解决方案在于 tm 包中包含的 removeWords() 函数。像之前一样，我们将使用 tm_map() 函数应用此映射到数据上，将 stopwords() 函数作为参数提供，以指示我们想要移除的单词。完整的命令如下：

> sms_corpus_clean <- tm_map(sms_corpus_clean,
    removeWords, stopwords()) 

由于 stopwords() 仅返回一个停用词向量，如果我们选择这样做，我们可以用我们自己的单词移除向量替换此函数调用。这样，我们可以根据喜好扩展或减少停用词列表，或者完全移除不同的单词集。

继续我们的清理过程，我们还可以使用内置的 removePunctuation() 转换从文本消息中消除任何标点符号：

> sms_corpus_clean <- tm_map(sms_corpus_clean, removePunctuation) 

removePunctuation() 转换会完全从文本中移除标点符号，这可能会导致意想不到的后果。例如，考虑以下应用方式会发生什么：

> removePunctuation("hello...world") 

[1] "helloworld" 

如所示，省略号后面的空白缺失导致单词 hello 和 world 被合并为一个单词。虽然这目前不是一个重大问题，但值得注意。

为了绕过 removePunctuation() 的默认行为，可以创建一个自定义函数来替换而不是移除标点符号：

> replacePunctuation <- function(x) {
    gsub("[[:punct:]]+", " ", x)
  } 

这使用 R 的 gsub() 函数将 x 中的任何标点符号替换为一个空格。然后可以使用 tm_map() 与其他转换一样使用 replacePunctuation() 函数。这里 gsub() 命令的奇怪语法是由于使用了 正则表达式，它指定了一个匹配文本字符的模式。正则表达式在 第十二章，高级数据准备 中有更详细的介绍。

文本数据的另一个常见标准化涉及通过称为 词干提取 的过程将单词还原到其根形式。词干提取过程将像 learned、learning 和 learns 这样的单词去除后缀，以便将它们转换为基本形式，learn。这允许机器学习算法将相关术语视为单一概念，而不是尝试为每个变体学习一个模式。

tm 包通过集成 SnowballC 包提供词干提取功能。在撰写本文时，SnowballC 并未默认与 tm 一起安装，因此如果您尚未安装，请使用 install.packages("SnowballC") 进行安装。

SnowballC 包由 Milan Bouchet-Valat 维护，并为基于 C 的 libstemmer 库提供 R 接口，该库本身基于 M. F. Porter 的“Snowball”词干提取算法，这是一种广泛使用的开源词干提取方法。有关更多详细信息，请参阅 snowballstem.org。

SnowballC 包提供了一个 wordStem() 函数，对于一个字符向量，它会返回其根形式的相同术语向量。例如，该函数正确地处理了之前描述的单词 learn 的变体：

> library(SnowballC)
> wordStem(c("learn", "learned", "learning", "learns")) 

[1] "learn"   "learn"   "learn"   "learn" 

要将 wordStem() 函数应用于整个文本文档集合，tm 包包括一个 stemDocument() 转换。我们使用 tm_map() 函数以与之前相同的方式应用此转换：

> sms_corpus_clean <- tm_map(sms_corpus_clean, stemDocument) 

如果在应用 stemDocument() 转换时收到错误消息，请确认您已安装 SnowballC 包。

在移除数字、停用词和标点符号后，然后也进行词干提取，文本消息就只剩下之前分隔现在缺失部分的空白空间。因此，我们文本清理过程的最后一步是使用内置的 stripWhitespace() 转换来移除额外的空白：

> sms_corpus_clean <- tm_map(sms_corpus_clean, stripWhitespace) 

下表显示了在清洗过程前后，短信语料库中的前三条消息。消息已被限制为最有趣的单词，并且已经移除了标点和大小写：

清洗前的短信消息	清洗后的短信消息
> as.character(sms_corpus[1:3])``[[1]] Hope you are having a good``week. Just checking in``[[2]] K..give back my thanks.``[[3]] Am also doing in cbe only.``But have to pay.	> as.character(sms_corpus_clean[1:3])``[[1]] hope good week just check``[[2]] kgive back thank``[[3]] also cbe pay

数据准备 – 将文本文档分割成单词

现在数据已经被处理成我们所期望的样子，最后一步是通过一个叫做分词的过程将消息分割成单个术语。术语是文本字符串的单个元素；在这种情况下，术语是单词。

如你所想，tm包提供了对 SMS 消息语料库进行分词的功能。DocumentTermMatrix()函数接受一个语料库并创建一个称为文档-术语矩阵（DTM）的数据结构，其中行表示文档（短信消息），列表示术语（单词）。

tm包还提供了一个术语-文档矩阵（TDM）的数据结构，它实际上是一个转置的文档-术语矩阵（DTM），其中行是术语，列是文档。为什么需要两者？有时，使用其中一个更方便。例如，如果文档数量很少，而单词列表很大，那么使用 TDM 可能是有意义的，因为通常显示多行比显示多列更容易。话虽如此，机器学习算法通常需要 DTM，因为列是特征，行是示例。

矩阵中的每个单元格存储一个数字，表示列中单词在行中文档中出现的次数。以下图仅显示了 SMS 语料库 DTM 的一小部分，因为完整的矩阵有 5,559 行和超过 7,000 列：

图 4.9：SMS 消息的 DTM 主要由零填充

表格中每个单元格为零的事实意味着列顶部的单词列表中的任何一个单词都没有出现在语料库的前五条消息中。这突出了为什么这种数据结构被称为稀疏矩阵的原因；矩阵中的绝大多数单元格都是零。用现实世界的术语来说，尽管每条消息必须至少包含一个单词，但任何单个单词出现在给定消息中的概率很小。

从tm语料库创建 DTM 稀疏矩阵只需要一个命令：

> sms_dtm <- DocumentTermMatrix(sms_corpus_clean) 

这将创建一个sms_dtm对象，它使用默认设置包含分词后的语料库，这些设置应用了最小的额外处理。默认设置是合适的，因为我们已经手动准备了语料库。

另一方面，如果我们还没有执行预处理，我们可以通过提供一个control参数选项的列表来覆盖默认设置。例如，要从原始、未经处理的短信语料库直接创建 DTM，我们可以使用以下命令：

> sms_dtm2 <- DocumentTermMatrix(sms_corpus, control = list(
    tolower = TRUE,
    removeNumbers = TRUE,
    stopwords = TRUE,
    removePunctuation = TRUE,
    stemming = TRUE
)) 

这将按照之前相同的顺序对 SMS 语料库进行相同的预处理步骤。然而，将sms_dtm与sms_dtm2比较，我们发现矩阵中的术语数量略有差异：

> sms_dtm 

<<DocumentTermMatrix (documents: 5559, terms: 6559)>>
Non-/sparse entries: 42147/36419334
Sparsity           : 100%
Maximal term length: 40
Weighting          : term frequency (tf) 

> sms_dtm2 

<<DocumentTermMatrix (documents: 5559, terms: 6961)>>
Non-/sparse entries: 43221/38652978
Sparsity           : 100%
Maximal term length: 40
Weighting          : term frequency (tf) 

这种差异的原因与预处理步骤顺序的微小差异有关。DocumentTermMatrix()函数在将文本字符串分割成单词之后，才对其应用清理函数。因此，它使用了一个稍微不同的停用词去除函数。因此，某些单词在清理之前分割的方式与清理后不同。

要强制两个先前的 DTM（文档-词矩阵）相同，我们可以用我们自己的使用原始替换函数的函数覆盖默认的停用词函数。只需将stopwords = TRUE替换为以下内容：

stopwords = function(x) { removeWords(x, stopwords()) } 

本章的代码文件包括创建一个相同的 DTM 所需的完整步骤，只需一个函数调用即可。

这些差异提出了清理文本数据的一个重要原则：操作的顺序很重要。考虑到这一点，思考过程早期步骤如何影响后续步骤非常重要。这里提供的顺序在许多情况下都会有效，但当过程更细致地针对特定数据集和用例定制时，可能需要重新思考。例如，如果你希望从矩阵中排除某些术语，考虑是否在词干提取之前或之后搜索它们。还要考虑移除标点符号（以及标点符号是否被替换为空白空间）如何影响这些步骤。

数据准备 – 创建训练和测试数据集

在我们的数据准备好分析后，我们现在需要将数据分为训练和测试数据集，以便在构建垃圾邮件分类器之后，它可以在之前未见过的数据上评估。然而，尽管我们需要保持分类器对测试数据集的内容不知情，但在数据清理和加工之后进行分割很重要。我们需要确保训练和测试数据集上发生的准备步骤完全相同。

我们将数据分为两部分：75%用于训练，25%用于测试。由于短信消息是随机排序的，我们可以简单地取前 4,169 条用于训练，剩下的 1,390 条用于测试。幸运的是，DTM 对象非常类似于数据框，可以使用标准的[行，列]操作进行分割。由于我们的 DTM 将短信消息存储为行，将单词存储为列，我们必须为每个请求特定的行范围和所有列：

> sms_dtm_train <- sms_dtm[1:4169, ]
> sms_dtm_test  <- sms_dtm[4170:5559, ] 

为了方便起见，保存一对向量，其中包含训练和测试矩阵中每行的标签也是有帮助的。这些标签没有存储在 DTM 中，因此我们需要从原始 sms_raw 数据框中提取它们：

> sms_train_labels <- sms_raw[1:4169, ]$type
> sms_test_labels  <- sms_raw[4170:5559, ]$type 

为了确认这些子集代表完整的短信数据集，让我们比较训练和测试数据框中垃圾邮件的比例：

> prop.table(table(sms_train_labels)) 

 ham      spam
0.8647158 0.1352842 

> prop.table(table(sms_test_labels)) 

 ham      spam
0.8683453 0.1316547 

训练数据和测试数据中大约有 13%的垃圾邮件。这表明垃圾邮件在这两个数据集中是平均分配的。

可视化文本数据 – 单词云

单词云是一种以视觉方式表示文本数据中单词出现频率的方法。云由围绕图形随机散布的单词组成。在文本中出现频率更高的单词以较大的字体显示，而较少出现的术语则以较小的字体显示。这种类型的图表作为观察社交媒体网站上趋势话题的一种方式而越来越受欢迎。

wordcloud 包提供了一个简单的 R 函数来创建此类图表。我们将使用它来可视化短信中的单词。比较垃圾邮件和正常邮件的云图将帮助我们判断我们的朴素贝叶斯垃圾邮件过滤器是否可能成功。如果您还没有这样做，请在 R 命令行中输入 install.packages("wordcloud") 和 library(wordcloud) 来安装和加载此包。

wordcloud 包是由 Ian Fellows 编写的。有关此包的更多信息，请访问他的博客blog.fellstat.com/?cat=11。

可以直接使用以下语法从 tm 语料库对象创建单词云：

> wordcloud(sms_corpus_clean, min.freq = 50, random.order = FALSE) 

这将根据我们准备好的短信语料库创建单词云。由于我们指定了 random.order = FALSE，云将按非随机顺序排列，高频词将放置在中心附近。如果我们不指定 random.order，则默认随机排列。

min.freq 参数指定单词在语料库中必须出现的次数，才能在云中显示。由于频率为 50 大约是语料库的 1%，这意味着一个单词必须至少出现在 1%的短信中才能包含在云中。

您可能会收到一条警告消息，指出 R 无法将所有单词拟合到图中。如果是这样，请尝试增加 min.freq 以减少云中的单词数量。使用 scale 参数减小字体大小也可能有所帮助。

生成的单词云应类似于以下示例：

图 4.10：展示所有短信中出现的单词的单词云

一个可能更有趣的可视化是将垃圾邮件和正常短信的云进行比较。由于我们没有为垃圾邮件和正常短信构建单独的语料库，这是注意 wordcloud() 函数的一个非常有用的特性的合适时机。给定一个原始文本字符串向量，它将在显示云之前自动应用常见的文本准备过程。

让我们使用 R 的 subset() 函数通过短信类型对 sms_raw 数据进行子集化。首先，我们将创建一个 type 为 spam 的子集：

> spam <- subset(sms_raw, type == "spam") 

接下来，我们将对 ham 子集做同样的事情：

> ham <- subset(sms_raw, type == "ham") 

注意双等号的使用。像许多编程语言一样，R 使用 == 来测试相等性。如果你不小心使用了单个等号，你最终会得到一个比你预期的更大的子集！

现在我们有两个数据框，spam 和 ham，每个数据框都有一个包含短信原始文本字符串的 text 特征。创建词云就像以前一样简单。这次，我们将使用 max.words 参数来查看每个集合中最常见的 40 个单词。scale 参数调整云中单词的最大和最小字体大小。你可以根据需要更改这些参数。以下代码展示了这一点：

> wordcloud(spam$text, max.words = 40, scale = c(3, 0.5))
> wordcloud(ham$text, max.words = 40, scale = c(3, 0.5)) 

注意，当运行此代码时，R 会提供警告信息，指出“转换丢弃文档”。这些警告与 wordcloud() 在接收到原始文本数据而不是术语矩阵时默认执行的 removePunctuation() 和 removeWords() 程序有关。基本上，有一些消息被排除在结果之外，因为在清理后没有剩余的消息文本。例如，带有文本 :) 表示笑脸表情的垃圾邮件消息在清理后被从集合中移除。这对词云没有问题，可以忽略这些警告。

生成的词云应该看起来与下面的类似。你有没有一种预感，哪一个是垃圾邮件云，哪一个是正常短信云？

图 4.11：并排显示的词云，描绘了垃圾邮件和正常短信消息

如你所猜，垃圾邮件云在左边。垃圾邮件包括诸如 call、free、mobile、claim 和 stop 等单词；这些术语根本不会出现在垃圾邮件云中。相反，垃圾邮件使用诸如 can、sorry、love 和 time 等单词。这些明显的差异表明，我们的朴素贝叶斯模型将有一些强有力的关键词来区分这两类。

数据准备 – 为频繁单词创建指示特征

数据准备过程的最后一步是将稀疏矩阵转换成可以用于训练 Naive Bayes 分类器的数据结构。目前，稀疏矩阵包括超过 6,500 个特征；这是每个至少出现在一条短信中的单词的特征。不太可能所有这些都有助于分类。为了减少特征数量，我们将删除在不到 5 条消息或训练数据中不到约 0.1% 的记录中出现的任何单词。

寻找频繁单词需要使用 tm 包中的 findFreqTerms() 函数。此函数接受一个 DTM 并返回一个包含至少出现最小次数的单词的字符向量。例如，以下命令显示在 sms_dtm_train 矩阵中至少出现五次的单词：

> findFreqTerms(sms_dtm_train, 5) 

函数的结果是一个字符向量，所以让我们将我们的频繁单词保存起来以备后用：

> sms_freq_words <- findFreqTerms(sms_dtm_train, 5) 

查看向量的内容显示，至少在 5 条短信中出现的术语有 1,139 个：

> str(sms_freq_words) 

 chr [1:1137] "£wk" "abiola" "abl" "abt" "accept" "access" "account" "across" "act" "activ" ... 

现在我们需要过滤我们的 DTM，只包含频繁单词向量中出现的术语。像之前一样，我们将使用数据框风格的 [row, col] 操作来请求 DTM 的特定部分，注意 DTM 的列名是基于 DTM 包含的单词。我们可以利用这个事实来限制 DTM 到特定的单词。由于我们想要所有行，但只有 sms_freq_words 向量中单词表示的列，我们的命令如下：

> sms_dtm_freq_train <- sms_dtm_train[ , sms_freq_words]
> sms_dtm_freq_test <- sms_dtm_test[ , sms_freq_words] 

训练和测试数据集现在包括 1,137 个特征，这些特征对应于至少在 5 条消息中出现的单词。

Naive Bayes 分类器通常在具有分类特征的数据上训练。这带来一个问题，因为稀疏矩阵中的单元格是数值型的，并衡量一个单词在消息中出现的次数。我们需要将其转换为表示是否出现的分类变量，是或否。

以下定义了一个 convert_counts() 函数，用于将计数转换为 Yes 或 No 字符串：

> convert_counts <- function(x) {
    x <- ifelse(x > 0, "Yes", "No")
} 

到目前为止，前面函数的一些部分应该看起来很熟悉。第一行定义了函数。语句 ifelse(x > 0, "Yes", "No") 将 x 中的值进行转换，如果值大于 0，则将其替换为 "Yes"；否则，将其替换为一个 "No" 字符串。最后，返回新转换的向量 x。

现在我们需要将 convert_counts() 应用到稀疏矩阵的每一列。你可能能够猜到执行此操作的 R 函数的名称。该函数简单地称为 apply()，其用法与之前使用的 lapply() 类似。

apply()函数允许对矩阵中的每一行或每一列使用一个函数。它使用MARGIN参数来指定是行还是列。在这里，我们将使用MARGIN = 2，因为我们感兴趣的是列（MARGIN = 1用于行）。转换训练和测试矩阵的命令如下：

> sms_train <- apply(sms_dtm_freq_train, MARGIN = 2,
    convert_counts)
> sms_test  <- apply(sms_dtm_freq_test, MARGIN = 2,
    convert_counts) 

结果将是两个字符类型的矩阵，每个矩阵的单元格都表示列中单词是否在任何时刻出现在行中代表的邮件中，是“是”还是“否”。

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

现在我们已经将原始短信消息转换成了可以由统计模型表示的格式，是时候应用朴素贝叶斯算法了。该算法将使用单词的存在或不存在来估计给定短信消息是垃圾邮件的概率。

我们将使用的朴素贝叶斯实现是在naivebayes包中。这个包由 Michal Majka 维护，是一个现代且高效的 R 实现。如果您还没有这样做，请确保在继续之前使用install.packages("naivebayes")和library(naivebayes)命令安装并加载该包。

许多机器学习方法在多个 R 包中实现，朴素贝叶斯也不例外。另一个选项是e1071包中的naiveBayes()，这在本书的旧版本中使用过，但在使用上几乎与naive_bayes()相同。本版使用的naivebayes包提供了更好的性能和更高级的功能，详情请访问其网站：majkamichal.github.io/naivebayes/。

与我们在上一章中用于分类的 k-NN 算法不同，朴素贝叶斯学习器的训练和使用分类发生在不同的阶段。尽管如此，如以下表格所示，这些步骤相当直接：

图 4.12：朴素贝叶斯分类语法

使用sms_train矩阵，以下命令训练了一个naive_bayes分类器对象，可以用来进行预测：

> sms_classifier <- naiveBayes(sms_train, sms_train_labels) 

执行前面的命令后，您可能会注意到以下输出：

There were 50 or more warnings (use warnings() to see the first 50) 

目前这没有什么好担心的；输入warnings()命令可以揭示这个问题的原因：

> warnings() 

Warning messages:
1: naive_bayes(): Feature £wk – zero probabilities are present. Consider Laplace smoothing.
2: naive_bayes(): Feature 60 biola – zero probabilities are present. Consider Laplace smoothing.
3: naive_bayes(): Feature abl – zero probabilities are present. Consider Laplace smoothing.
4: naive_bayes(): Feature abt – zero probabilities are present. Consider Laplace smoothing.
5: naive_bayes(): Feature accept – zero probabilities are present. Consider Laplace smoothing. 

这些警告是由那些在零垃圾邮件或零正常邮件中出现的单词引起的，由于它们关联的零概率，它们在分类过程中具有否决权。例如，因为单词接受仅在训练数据中的正常邮件中出现，这并不意味着每个包含这个单词的未来邮件都应该自动被分类为正常邮件。

使用前面描述的拉普拉斯估计器可以轻松解决这个问题，但到目前为止，我们将使用laplace = 0来评估这个模型，这是模型的默认设置。

第 4 步 – 评估模型性能

为了评估 SMS 分类器，我们需要在测试数据中未见过的新消息上测试其预测结果。回想一下，未见过的新消息特征存储在一个名为 sms_test 的矩阵中，而类别标签（垃圾邮件或正常邮件）存储在一个名为 sms_test_labels 的向量中。我们训练的分类器已被命名为 sms_classifier。我们将使用这个分类器来生成预测，然后比较预测值和真实值。

使用 predict() 函数进行预测。我们将这些预测存储在一个名为 sms_test_pred 的向量中。我们只需向这个函数提供我们分类器和测试数据集的名称，如下所示：

> sms_test_pred <- predict(sms_classifier, sms_test) 

为了将预测值与真实值进行比较，我们将使用 gmodels 包中的 CrossTable() 函数，我们在前面的章节中使用过它。这次，我们将添加一些额外的参数来消除不必要的单元格比例，并使用 dnn 参数（维度名称）来重新标记行和列，如下面的代码所示：

> library(gmodels)
> CrossTable(sms_test_pred, sms_test_labels,
    prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE,
    dnn = c('predicted', 'actual')) 

这产生了以下表格：

Total Observations in Table:  1390 

             | actual 
   predicted |       ham |      spam | Row Total |
-------------|-----------|-----------|-----------|
         ham |      1201 |        30 |      1231 |
             |     0.864 |     0.022 |           |
-------------|-----------|-----------|-----------|
        spam |         6 |       153 |       159 |
             |     0.004 |     0.110 |           |
-------------|-----------|-----------|-----------|
Column Total |      1207 |       183 |      1390 |
-------------|-----------|-----------|-----------| 

观察表格，我们可以看到总共只有 6 + 30 = 36 条 1,390 条短信被错误分类（2.6%）。在这些错误中，有 6 条在 1,207 条正常邮件中被错误地识别为垃圾邮件，以及 30 条在 183 条垃圾邮件中被错误地标记为正常邮件。考虑到我们在项目上投入的少量努力，这种性能水平似乎相当令人印象深刻。这个案例研究说明了为什么朴素贝叶斯在文本分类中如此经常被使用：直接使用，它表现出令人惊讶的良好性能。

另一方面，被错误分类为垃圾邮件的六条合法消息可能会对我们过滤算法的部署造成重大问题，因为过滤器可能会使一个人错过一条重要的短信。我们应该尝试看看我们是否可以稍微调整模型以获得更好的性能。

第 5 步 – 提高模型性能

你可能还记得，我们在训练模型时没有为拉普拉斯估计器设置值；事实上，很难错过 R 警告我们超过 50 个特征概率为零的信息！为了解决这个问题，我们将像以前一样构建一个朴素贝叶斯模型，但这次将 laplace = 1 设置：

> sms_classifier2 <- naiveBayes(sms_train, sms_train_labels,
    laplace = 1) 

接下来，我们将像以前一样进行预测：

> sms_test_pred2 <- predict(sms_classifier2, sms_test) 

最后，我们将使用交叉表比较预测类别与实际分类：

> CrossTable(sms_test_pred2, sms_test_labels,
    prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE,
    dnn = c('predicted', 'actual')) 

这产生了以下表格：

Total Observations in Table:  1390 
             | actual 
   predicted |       ham |      spam | Row Total |
-------------|-----------|-----------|-----------|
         ham |      1202 |        28 |      1230 |
             |     0.865 |     0.020 |           |
-------------|-----------|-----------|-----------|
        spam |         5 |       155 |       160 |
             |     0.004 |     0.112 |           |
-------------|-----------|-----------|-----------|
Column Total |      1207 |       183 |      1390 |
-------------|-----------|-----------|-----------| 

通过设置laplace = 1添加拉普拉斯估计器，将误报（将正常邮件错误分类为垃圾邮件）的数量从 6 减少到 5，将漏报（将垃圾邮件错误分类为正常邮件）的数量从 30 减少到 28。尽管这似乎是一个微小的变化，但考虑到模型的准确率已经相当令人印象深刻，这是一个重大的改进。在进一步调整模型之前，我们需要谨慎，因为在过滤垃圾邮件时，保持过于激进和过于被动之间的平衡是很重要的。用户更希望少量垃圾邮件能够通过过滤器，而不是另一种情况，即正常邮件被过度过滤。

摘要

在本章中，我们学习了使用天真贝叶斯进行分类。该算法构建用于估计新示例属于各种类别的概率表。这些概率是通过称为贝叶斯定理的公式计算的，该公式指定了依赖事件之间的关系。尽管贝叶斯定理在计算上可能很昂贵，但一个简化的版本，即所谓的“天真”假设关于特征的独立性，能够处理更大的数据集。

天真贝叶斯分类器常用于文本分类。为了说明其有效性，我们在涉及垃圾短信的分类任务中使用了天真贝叶斯。准备文本数据进行分析需要使用专门的 R 包进行文本处理和可视化。最终，该模型能够正确地将 97%以上的短信分类为垃圾邮件或正常邮件。

在下一章中，我们将探讨两种更多的机器学习方法。每种方法通过将数据划分为相似值组来进行分类。正如你很快会发现的那样，这些方法本身非常有用。然而，展望未来，这些基本算法也成为了当今一些最强大的机器学习方法的重要基础，这些方法将在第十四章“构建更好的学习者”中介绍。

加入我们书籍的 Discord 空间

加入我们的 Discord 社区，与志同道合的人相聚，并和超过 4000 人一起学习：

packt.link/r

第五章：分而治之 - 使用决策树和规则进行分类

在决定是否接受工作邀请时，许多人首先列出利弊清单，然后使用简单的规则消除选项。例如，他们可能会决定，“如果我必须通勤超过一小时，我会不开心，”或者“如果我收入低于 5 万美元，我就无法养家糊口。”通过这种方式，预测个人未来职业幸福感的复杂决策可以简化为一系列简单决策。

本章介绍了决策树和规则学习器——两种机器学习方法，它们也能从简单选择集合中做出复杂决策。这些方法以逻辑结构的形式呈现其知识，无需统计知识即可理解。这一特性使得这些模型在商业策略和流程改进方面特别有用。

到本章结束时，你将学到：

• 树和规则如何“贪婪”地将数据分割成有趣的片段

• 最常见的决策树和分类规则学习器，包括 C5.0、1R 和 RIPPER 算法

• 如何使用这些算法执行现实世界的分类任务，例如识别有风险的银行贷款和有毒蘑菇

我们将首先检查决策树，然后看看分类规则。然后，我们将通过预览后续章节来总结我们所学的知识，这些章节讨论了使用树和规则作为更高级机器学习技术基础的方法。

理解决策树

决策树学习器是强大的分类器，它们利用树结构来模拟特征与潜在结果之间的关系。如图所示，这种结构之所以得名，是因为它反映了真实树木的生长方式，从底部宽大的树干开始，随着向上生长，逐渐分裂成越来越窄的树枝。同样地，决策树分类器使用分支决策的结构来引导示例进入最终的预测类别值。

为了更好地理解这在实践中是如何工作的，让我们考虑以下树，它预测是否应该接受工作邀请。一个正在考虑的工作邀请从根节点开始，然后通过决策节点，这些节点需要根据工作的属性做出选择。这些选择将数据分割成表示决策潜在结果的分支。在这里，它们被描绘为是或否的结果，但在其他情况下，可能存在超过两种可能性。

如果可以做出最终决策，树将终止在叶节点（也称为终端节点），这些节点表示一系列决策的结果所采取的行动。在预测模型的情况下，叶节点提供了给定树中的事件序列的预期结果。

图 5.1：一个表示确定是否接受新工作邀请过程的决策树

决策树算法的一个巨大优势是，类似于流程图的树结构不仅用于机器的内部使用。在模型创建之后，许多决策树算法会将结果结构以人类可读的格式输出。这为理解模型如何以及为什么在特定任务中表现良好或不好提供了洞察。这也使得决策树特别适合于需要出于法律原因使分类机制透明或需要与他人共享结果以告知未来业务实践的应用。考虑到这一点，一些潜在用途包括以下内容：

• 需要明确记录并消除偏见的信用评分模型，其中导致申请人被拒绝的标准需要被清楚地记录

• 市场研究客户行为，如满意度或流失，这些将与管理层或广告机构共享

• 基于实验室测量、症状或疾病进展率的医疗状况诊断

尽管前面的应用说明了树在告知决策过程的价值，但这并不意味着它们的效用到此为止。实际上，决策树是单一最广泛使用的机器学习技术之一，可以应用于几乎任何类型的数据——通常具有出色的即插即用性能。

话虽如此，尽管它们的适用性很广，但值得注意的是，在某些情况下，树可能不是理想的选择。这包括数据具有许多名义特征和多个级别或大量数值特征的任务。这些情况可能导致非常多的决策和过于复杂的树。它们还可能加剧决策树过度拟合数据的倾向，尽管正如我们很快就会看到的，通过调整一些简单的参数，甚至这种弱点也可以克服。

分而治之

决策树是使用一种称为递归划分的启发式方法构建的。这种方法也通常被称为分而治之，因为它将数据分割成子集，然后这些子集被反复分割成更小的子集，如此类推，直到算法确定子集中的数据足够同质，或者满足另一个停止标准。

要了解如何分割数据集可以创建决策树，想象一个根节点，它将成长为一棵成熟的树。最初，根节点代表整个数据集，因为还没有发生分割。在这里，决策树算法必须选择一个特征进行分割；理想情况下，它选择对目标类最具有预测性的特征。

然后将示例根据该特征的独特值进行分组，并形成第一组树分支。

沿着每个分支向下工作，算法继续分割和征服数据，每次选择最佳候选特征来创建另一个决策节点，直到达到停止标准。如果以下条件之一成立，分割和征服可能会在节点处停止：

• 节点上的所有（或几乎所有）示例都属于同一类别

• 没有剩余的特征来区分示例

• 树已经增长到预定义的大小限制

为了说明树构建过程，让我们考虑一个简单的例子。想象一下，你为一家好莱坞工作室工作，你的角色是决定工作室是否应该继续制作由有潜力的新作者提出的剧本。度假回来后，你的桌子堆满了提案。没有时间逐个阅读每个提案，你决定开发一个决策树算法来预测潜在电影是否会落入以下三个类别之一：评论家成功、主流热门或票房失败。

为了创建决策树模型而获取数据，你转向电影制片厂的档案来检查导致公司最近 30 部发行作品成功或失败的因素。你很快注意到电影的预估拍摄预算、主演角色的 A 名单明星数量以及电影的成功程度之间存在关联。对这个发现感到兴奋，你制作了一个散点图来展示这种模式：

图 5.2：散点图展示了电影预算与名人数量之间的关系

使用分割和征服策略，你可以从这个数据中构建一个简单的决策树。首先，为了创建树的根节点，你分割了表示名人数量的特征，将电影分为有显著数量的 A 名单明星和无显著数量的 A 名单明星的组：

图 5.3：决策树的第一次分割将数据分为高名人数量和低名人数量电影

接下来，在拥有更多名人的电影组中，你在有高预算和无高预算的电影之间进行另一个分割：

图 5.4：决策树的第二次分割进一步将高名人数量电影分为低预算和高预算两类

到目前为止，你已经将数据分为三个组。图的最左上角组完全由获得评论家赞誉的电影组成。这个组的特点是名人数量多而预算相对较低。在右上角，几乎所有电影都是票房大赢家，预算高且名人众多。最后一个组，虽然星光不足，但预算从小到大不等，包含了失败的作品。

如果需要，你可以通过在预算和名人数量越来越具体的范围内分割数据来继续细分和征服数据，直到当前所有错误分类的值都正确分类在其自己的小分区中。然而，以这种方式过度拟合决策树是不明智的。尽管算法可以无限分割数据，但过于具体的决策并不总是能更广泛地推广。因此，你选择在这里停止算法，因为每个组中超过 80%的例子都来自单个类别。这是决策树模型的停止标准。

你可能已经注意到，对角线可能更干净地分割了数据。这是决策树知识表示的一个局限性，它使用轴平行分割。每个分割只考虑一个特征的事实阻止了决策树形成更复杂的决策边界。例如，可以通过一个询问“名人的数量是否大于估计的预算？”的决策来创建一条对角线。如果是这样，那么“它将是一个关键的成功。”

预测电影未来成功的模型可以用一个简单的树来表示，如下面的图所示。树中的每一步都显示了落入每个类别的示例比例，这显示了随着分支接近叶子，数据如何变得更加同质化。为了评估一个新电影剧本，沿着每个决策分支前进，直到预测出剧本的成功或失败。使用这种方法，你将能够快速识别剧本库中最有希望的选项，并回到更重要的工作，例如撰写奥斯卡颁奖典礼的获奖感言！

图 5.5：基于历史电影数据的决策树可以预测未来电影的性能

由于现实世界的数据包含超过两个特征，决策树很快就会比这复杂得多，具有更多的节点、分支和叶子。在下一节中，你将了解一个流行的自动构建决策树模型的算法。

C5.0 决策树算法

决策树有多种实现方式，但其中最著名的一种是 C5.0 算法。该算法由计算机科学家 J. Ross Quinlan 开发，是他先前算法 C4.5 的改进版本，而 C4.5 本身又是他迭代二分器 3（ID3）算法的改进。尽管 Quinlan 将 C5.0 推广给商业客户（详情请见www.rulequest.com/），但该算法的单线程版本源代码已被公开，因此被整合到 R 等程序中。

为了进一步混淆问题，一个流行的基于 Java 的开源 C4.5 替代品，名为J48，包含在 R 的RWeka包中（在本章后面介绍）。由于 C5.0、C4.5 和 J48 之间的差异很小，本章中的原则适用于这三种方法中的任何一种，并且算法应被视为同义的。

C5.0 算法已成为生成决策树的行业标准，因为它可以直接解决大多数类型的问题。与其他高级机器学习模型相比，例如在第七章中描述的黑盒方法 - 神经网络和支持向量机，C5.0 构建的决策树通常表现几乎一样好，但更容易理解和部署。此外，如以下表格所示，该算法的缺点相对较小，并且可以很大程度上避免。

优点	缺点

|

• 一种通用的分类器，在许多类型的问题上表现良好

• 高度自动化的学习过程，可以处理数值或名义特征，以及缺失数据

• 排除不重要的特征

• 可用于小型和大型数据集

• 生成的模型可以在没有数学背景的情况下进行解释（对于相对较小的树）

• 比其他复杂模型更高效

|

• 决策树模型往往偏向于具有大量级别的特征的分割

• 容易过拟合或欠拟合模型

• 由于依赖于轴平行分割，可能难以模拟某些关系

• 训练数据的小幅变化可能导致决策逻辑发生大幅变化

• 大型树可能难以解释，它们做出的决策可能看起来不符合直觉

|

为了保持简单，我们之前的决策树示例忽略了机器如何采用分而治之策略涉及的数学。让我们更详细地探讨这一点，以检查这种启发式方法在实际中的工作方式。

选择最佳分割点

决策树将面临的第一挑战是确定要分割哪个特征。在前面的例子中，我们寻找一种分割数据的方法，使得结果分区主要包含单个类别的示例。

一个示例子集只包含单个类别的程度被称为纯度，而只由单个类别组成的任何子集都称为纯集。

有各种纯度度量可以用来识别最佳的决策树分割候选者。C5.0 使用熵，这是从信息理论中借用的一个概念，它量化了类别值集合中的随机性或无序性。具有高熵的集合非常多样化，并且关于可能也属于该集合的其他项目提供的信息很少，因为没有明显的共同点。决策树希望找到可以减少熵的分割点，从而最终增加组内的同质性。

通常，熵以 比特 为单位测量。如果只有两个可能的类别，熵值可以从 0 到 1 变化。对于 n 个类别，熵的范围从 0 到 log2。在每种情况下，最小值表示样本完全同质，而最大值表示数据尽可能多样化，并且没有组有哪怕是很小的多数。

在数学概念中，熵被定义为：

在这个公式中，对于给定数据段 (S)，术语 c 指的是类别级别数，而 p[i] 指的是落在第 i 个类别级别中的值的比例。

例如，假设我们有一个包含两个类别的数据分区：红色（60%）和白色（40%）。我们可以计算熵如下：

> -0.60 * log2(0.60) - 0.40 * log2(0.40) 

[1] 0.9709506 

我们可以可视化所有可能的二分类排列的熵。如果我们知道一个类别的示例比例是 x，那么另一个类别的比例是 (1 – x)。然后，使用 curve() 函数，我们可以绘制所有可能的 x 值的熵：

> curve(-x * log2(x) - (1 - x) * log2(1 - x),
        col = "red", xlab = "x", ylab = "Entropy", lwd = 4) 

这导致了以下图表：

图 5.6：在二分类结果中，随着一个类别的比例变化，总熵

如 x = 0.50 处的峰值所示，50-50 的分裂会产生最大熵。当一个类别越来越多地主导另一个类别时，熵减少到零。

要使用熵来确定最优的分裂特征，算法会计算在每种可能的特征上分裂所导致的同质性的变化，这种度量称为 信息增益。特征 F 的信息增益是分裂前段 (S[1]) 的熵与分裂后的分区 (S[2]) 的熵之间的差值：

InfoGain(F) = Entropy(S[1]) – Entropy(S[2])

一个复杂的问题是，在分裂之后，数据被分成多个分区。因此，计算 Entropy(S[2]) 的函数需要考虑分裂后所有分区的总熵。它是通过根据所有记录落在该分区中的比例来加权每个分区的熵来做到这一点的。这可以用以下公式表示：

简而言之，分裂产生的总熵是每个 n 个分区的熵的总和，每个分区的熵都按落在该分区中的示例比例 (w[i]) 加权。

信息增益越高，特征在分裂该特征后创建同质组的性能就越好。如果信息增益为零，则在该特征上分裂不会减少熵。另一方面，最大信息增益等于分裂前的熵。这表明分裂后的熵为零，这意味着分裂产生了完全同质的组。

之前的公式假设名义特征，但决策树在分割数值特征时也使用信息增益。为此，一种常见的做法是测试各种分割，将值分为大于或小于阈值的组。这把数值特征简化为两个级别的分类特征，从而允许像往常一样计算信息增益。产生最大信息增益的数值切割点被选为分割点。

虽然 C5.0 使用了信息增益，但它不是构建决策树时可以使用的唯一分割标准。其他常用的标准包括基尼指数、卡方统计量和增益率。关于这些（以及许多其他）标准的综述，请参阅Mingers, J, Machine Learning, 1989, Vol. 3, pp. 319-342。

剪枝决策树

如前所述，决策树可以无限增长，选择分割特征并将数据分成越来越小的分区，直到每个示例都被完美分类或算法耗尽可用于分割的特征。然而，如果树长得过于庞大，它所做的许多决策将过于具体，模型将过度拟合训练数据。剪枝决策树的过程涉及减小其大小，以便更好地泛化到未见数据。

解决这个问题的方法之一是在树达到一定数量的决策或决策节点只包含少量示例时停止树的生长。这被称为早期停止或预剪枝决策树。由于树避免了做无谓的工作，这是一种吸引人的策略。然而，这种方法的一个缺点是，无法知道树是否会错过它如果长得更大可能会学习到的微妙但重要的模式。

另一种称为后剪枝的方法涉及有意生长一个过大的树，并剪枝叶节点以将树的大小减小到更合适的水平。这通常比预剪枝更有效，因为在不首先生长树的情况下很难确定决策树的最佳深度。在稍后剪枝树允许算法确信所有重要的数据结构都已发现。

剪枝操作的实现细节非常技术性，超出了本书的范围。关于一些可用方法的比较，请参阅Esposito, F, Malerba, D, Semeraro, G, IEEE Transactions on Pattern Analysis and Machine Intelligence, 1997, Vol. 19, pp. 476-491。

C5.0 算法的一个优点是它在剪枝方面有明确的观点——它使用合理的默认值自动处理许多决策。其整体策略是在剪枝后对树进行修剪。它首先生长出一个大的树，该树过度拟合了训练数据。随后，移除对分类错误影响较小的节点和分支。在某些情况下，整个分支会被移动到树的更高位置或被更简单的决策所取代。这些移除分支的过程分别被称为子树提升和子树替换。

在过度拟合和欠拟合之间取得正确的平衡有点像一门艺术，但如果模型准确性至关重要，那么花些时间尝试不同的剪枝选项以查看是否可以提高测试数据集的性能可能是值得的。正如你很快就会看到的，C5.0 算法的一个优点是它非常容易调整训练选项。

示例 – 使用 C5.0 决策树识别风险银行贷款

2007-2008 年的全球金融危机突出了银行实践中透明度和严谨性的重要性。由于信贷的可用性有限，银行收紧了其贷款系统，并转向机器学习以更准确地识别风险贷款。

由于决策树具有高准确性和用普通语言制定统计模型的能力，因此在银行业中得到广泛应用。由于许多国家的政府仔细监控贷款实践的公平性，因此高管必须能够解释为什么某个申请人被拒绝贷款而另一个被批准。这些信息对希望确定其信用评级为何不满意的客户也很有用。

自动化信用评分模型可能用于信用卡邮寄和即时在线批准流程。在本节中，我们将使用 C5.0 决策树开发一个简单的信用批准模型。我们还将看到如何调整模型结果以最小化导致财务损失的错误。

第 1 步 – 收集数据

我们信用模型的动力在于识别与贷款违约高风险相关的因素。为此，我们必须获取关于过去银行贷款的数据以及当时在信用申请时可能可用的贷款申请人信息。

具有这些特征的数据集由汉堡大学的 Hans Hofmann 捐赠给 UCI 机器学习仓库 (archive.ics.uci.edu/ml)。该数据集包含来自德国一家信用机构的贷款信息。

本章中展示的数据集已经对原始数据进行了轻微修改，以消除一些预处理步骤。为了跟随示例，请从本章的 Packt Publishing GitHub 仓库下载 credit.csv 文件，并将其保存到您的 R 工作目录中。

信用数据集包括 1,000 个贷款示例，以及一组表示贷款和贷款申请人特征的数值和名义特征。一个类别变量表示贷款是否违约。让我们看看我们是否能识别出任何预测这种结果的模式。

第 2 步 – 探索和准备数据

如我们之前所做的那样，我们将使用read.csv()函数导入数据。现在，因为字符数据完全是分类的，我们可以将stringsAsFactors = TRUE设置为自动将所有字符类型列转换为结果数据框中的因子：

> credit <- read.csv("credit.csv", stringsAsFactors = TRUE) 

我们可以通过检查str()函数输出的前几行来检查结果对象：

> str(credit) 

'data.frame':1000 obs. of  17 variables:
 $ checking_balance : Factor w/ 4 levels "< 0 DM","> 200 DM",..
 $ months_loan_duration: int  6 48 12 42 24 36 24 36 12 30 ...
 $ credit_history : Factor w/ 5 levels "critical","good",..
 $ purpose : Factor w/ 6 levels "business","car",..
 $ amount : int  1169 5951 2096 7882 4870 9055 2835 6948 ... 

我们看到预期的 1,000 个观测值和 17 个特征，这些特征是因子和整型数据类型的组合。

让我们看看几个似乎可以预测违约的贷款特征的table()输出。申请人的检查和储蓄账户余额被记录为分类变量：

> table(credit$checking_balance) 

 < 0 DM   > 200 DM 1 - 200 DM    unknown
       274         63        269        394 

> table(credit$savings_balance) 

 < 100 DM > 1000 DM  100 - 500 DM 500 - 1000 DM   unknown
          603        48           103            63       183 

检查和储蓄账户余额可能是贷款违约状态的重要预测指标。请注意，由于贷款数据来自德国，这些值使用的是德国马克（DM），这是欧元采用之前德国使用的货币。

一些贷款的特征是数值型的，例如其期限和请求的信贷金额：

> summary(credit$months_loan_duration) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    4.0    12.0    18.0    20.9    24.0    72.0 

> summary(credit$amount) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    250    1366    2320    3271    3972   18424 

贷款金额从 250 DM 到 18,420 DM 不等，期限为 4 到 72 个月。它们的中间值为 2,320 DM，中间期限为 18 个月。

default向量表示贷款申请人是否能够满足约定的还款条款，或者他们是否违约。在这个数据集中，总共有 30%的贷款违约：

> table(credit$default) 

no yes
700 300 

对于银行来说，高违约率是不可取的，因为这意味着银行不太可能完全收回其投资。如果我们成功，我们的模型将识别出高违约风险的申请人，使银行能够在发放资金之前拒绝信贷申请。

数据准备 – 创建随机训练和测试数据集

如我们在前面的章节中所做的那样，我们将把我们的数据分成两部分：一个用于构建决策树的训练数据集和一个用于评估其在新数据上性能的测试数据集。我们将使用 90%的数据进行训练，10%的数据进行测试，这将为我们提供 100 条记录来模拟新申请人。在这里使用 90-10 的分割而不是更常见的 75-25 分割，是因为信用数据集相对较小；鉴于预测贷款违约是一个具有挑战性的学习任务，我们需要尽可能多的训练数据，同时还要保留足够的测试样本。

在第十章“评估模型性能”中介绍了用于训练和评估相对较小数据集的更复杂方法。

由于前几章使用的数据是随机排序的，我们只需通过取记录的第一个子集用于训练，剩余的子集用于测试，就可以将数据集分为两部分。相比之下，信用数据集并不是随机排序的，这使得先前的做法不明智。假设银行是按照贷款金额对数据进行排序的，最大的贷款位于文件末尾。如果我们用前 90%的数据用于训练，剩下的 10%用于测试，那么我们将在只有小额贷款上训练模型，并在大额贷款上测试模型。显然，这可能会出现问题。

我们将通过在信用数据的随机样本上训练模型来解决此问题。随机样本简单来说就是一个随机选择记录子集的过程。在 R 中，sample()函数用于执行随机抽样。然而，在付诸实践之前，一个常见的做法是设置一个种子值，这会导致随机化过程遵循一个可以后来复制的序列。这似乎违背了生成随机数的初衷，但这样做有很好的理由。通过set.seed()函数提供种子值可以确保，如果分析在未来重复进行，将获得相同的结果。

你可能会想知道一个所谓的随机过程如何被播种以产生相同的结果。这是因为计算机使用一种称为伪随机数生成器的数学函数来创建看似非常随机的随机数序列，但实际上，如果知道序列中的前一个值，它们是非常可预测的。在实践中，现代伪随机数序列几乎与真正的随机序列无法区分，但它们的好处是计算机可以快速轻松地生成它们。

以下命令使用带有种子值的sample()。请注意，set.seed()函数使用任意值9829。省略此种子值将导致您的训练和测试分割与本章其余部分所示的不同。以下命令从 1 到 1,000 的整数序列中随机选择 900 个值：

> set.seed(9829)
> train_sample <- sample(1000, 900) 

如预期的那样，生成的train_sample对象是一个包含 900 个随机整数的向量：

> str(train_sample) 

int [1:900] 653 866 119 152 6 617 250 343 367 138 ... 

通过使用这个向量从信用数据中选择行，我们可以将其分为我们想要的 90%训练数据和 10%测试数据集。回想一下，在测试记录的选择中使用的否定运算符（-字符）告诉 R 选择不在指定行中的记录；换句话说，测试数据只包括不在训练样本中的行：

> credit_train <- credit[train_sample, ]
> credit_test  <- credit[-train_sample, ] 

如果随机化做得正确，我们应在每个数据集中大约有 30%的贷款出现违约：

> prop.table(table(credit_train$default)) 

 no       yes 
0.7055556 0.2944444 

> prop.table(table(credit_test$default)) 

 no  yes 
0.65 0.35 

训练集和测试集在贷款违约的分布上大致相似，因此我们现在可以构建我们的决策树。如果比例差异很大，我们可能会决定重新采样数据集，或者尝试更复杂的采样方法，例如在第十章“评估模型性能”中介绍的方法。

如果您的结果不完全匹配，请确保在创建train_sample向量之前立即运行了set.seed(9829)命令。请注意，R 的默认随机数生成器在 R 版本 3.6.0 中发生了变化，如果在此代码在早期版本上运行，则结果将不同。这也意味着这里的结果与本书先前版本中的结果略有不同。

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

我们将使用C50包中的 C5.0 算法来训练我们的决策树模型。如果您尚未安装，请使用install.packages("C50")安装该包，并使用library(C50)将其加载到 R 会话中。

以下语法框列出了在构建决策树时使用的一些最常见参数。与之前使用的机器学习方法相比，C5.0 算法提供了许多更多的方式来定制模型以适应特定的学习问题。

图 5.7：C5.0 决策树语法

C5.0()函数使用一种称为R 公式接口的新语法来指定要训练的模型。公式语法使用~运算符（称为波浪号）来表示目标变量与其预测变量之间的关系。要学习的类别变量放在波浪号的左侧，预测特征写在右侧，由+运算符分隔。

如果您想建模目标y与预测变量x1和x2之间的关系，您将公式写成y ~ x1 + x2。要包含模型中的所有变量，使用点字符。例如，y ~ .指定了y与数据集中所有其他特征之间的关系。

R 公式接口在许多 R 函数中使用，并提供了一些强大的功能来描述预测变量之间的关系。我们将在后面的章节中探讨一些这些功能。然而，如果您急于预览，请随时使用?formula命令阅读文档。

对于信用审批模型的第一次迭代，我们将使用默认的 C5.0 设置，如下面的代码所示。目标类别命名为default，所以我们将其放在波浪号~的左侧，后面跟着一个点表示credit_train数据框中的所有其他列都将用作预测变量：

> credit_model <- C5.0(default ~ ., data = credit_train) 

credit_model对象现在包含一个 C5.0 决策树。我们可以通过输入其名称来查看一些关于树的基本数据：

> credit_model 

Call:
C5.0.formula(formula = default ~ ., data = credit_train)
Classification Tree
Number of samples: 900 
Number of predictors: 16 
Tree size: 67 
Non-standard options: attempt to group attributes 

输出显示了关于树的一些简单事实，包括生成它的函数调用、用于构建树的特性数量（标记为predictors）和示例（标记为samples）。还包括树的大小为67，这表明树有 67 个决策深度——比我们迄今为止考虑的示例树大得多！

要查看树的决策，我们可以在模型上调用summary()函数：

> summary(credit_model) 

这会导致以下输出，其中已截断以仅显示前几行：

> summary(credit_model) 

Call:
C5.0.formula(formula = default ~ ., data = credit_train)
C5.0 [Release 2.07 GPL Edition]
-------------------------------
Class specified by attribute `outcome'
Read 900 cases (17 attributes) from undefined.data
Decision tree:
checking_balance in {> 200 DM,unknown}: no (415/55)
checking_balance in {< 0 DM,1 - 200 DM}:
:...credit_history in {perfect,very good}: yes (59/16)
    credit_history in {critical,good,poor}:
    :...months_loan_duration > 27:
        :...dependents > 1:
        :   :...age <= 45: no (12/2)
        :   :   age > 45: yes (2) 

上述输出显示了决策树的一些最初分支。前三行可以用普通语言表示为：

1. 如果支票账户余额未知或超过 200 DM，则将其分类为“不太可能违约”

2. 否则，如果支票账户余额小于零 DM 或介于 1 到 200 DM 之间...

3. …如果信用记录是完美或非常好的，则将其分类为“可能违约”

括号中的数字表示满足该决策标准的示例数量和被该决策错误分类的示例数量。例如，在第一行中，415/55表示，在达到决策的 415 个示例中，有 55 个被错误地分类为“不太可能违约”。换句话说，415 个申请人中有 55 人实际上违约了，尽管模型预测他们不会违约。

有时候，一棵树会导致一些逻辑上不太合理的决策。例如，为什么信用记录完美或非常好的申请人可能会违约，而那些支票账户余额未知的人却不太可能违约？这样的矛盾规则有时会出现。它们可能反映了数据中的真实模式，或者可能是统计异常。在两种情况下，调查这些奇怪的决策以查看树的逻辑是否适合商业用途都是非常重要的。

在树之后，summary(credit_model)输出显示了一个混淆矩阵，它是一个交叉表，指示模型在训练数据中的错误分类记录：

Evaluation on training data (900 cases):
        Decision Tree   
      ----------------  
      Size      Errors  
        66  118(13.1%)   <<
       (a)   (b)    <-classified as
      ----  ----
       604    31    (a): class no
        87   178    (b): class yes 

Errors标题显示，该模型正确分类了 900 个训练实例中的 878 个，错误率为 13.1%。共有 31 个实际值为no的实例被错误地分类为yes（假阳性），而 87 个值为yes的实例被错误地分类为no（假阴性）。鉴于决策树倾向于过度拟合训练数据，这里报告的错误率，即基于训练数据性能的错误率，可能过于乐观。因此，将决策树应用于未见过的测试数据集尤为重要，我们将在稍后进行。

输出还包括一个标记为Attribute usage的部分，它提供了关于决策树模型中使用的重要预测因子的一般感觉。以下是一些输出内容的前几行：

Attribute usage:
    100.00%       checking_balance
    53.89%        credit_history
    47.33%        months_loan_duration
    26.11%        purpose
    24.33%        savings_balance
    18.22%        job
    12.56%        dependents
    12.11%        age 

决策树输出中的属性使用统计信息指的是在训练数据中使用列出的特征进行最终预测的行百分比。例如，100%的行需要checking_balance特征，因为检查账户余额在树的第一个分裂点被使用。第二个分裂使用credit_history，但 46.11%的行已经根据检查账户余额被分类为非违约。这留下了只有 53.89%的行需要考虑申请人的信用历史。在这个列表的底部，只有 12.11%的例子需要申请人的年龄来做出预测，这表明申请人的年龄不如他们的检查账户余额或信用历史重要。

这条信息，连同树结构本身，提供了对模型工作方式的洞察。这两者都易于理解，即使没有统计学背景也是如此。当然，即使模型易于理解，如果它不能做出准确的预测，那么它也是无用的，因此我们现在将对其进行更正式的性能评估。

C5.0 决策树模型可以使用plot()函数进行可视化，该函数依赖于partykit包中的功能。不幸的是，这仅适用于相对较小的决策树。例如，我们可以通过输入plot(credit_model)来可视化我们的决策树，但除非你有非常大的显示屏，否则由于树中节点和分裂的数量很大，生成的图表可能会显得杂乱无章。

第 4 步 – 评估模型性能

要将我们的决策树应用于测试数据集，我们使用predict()函数，如下面的代码行所示：

> credit_pred <- predict(credit_model, credit_test) 

这创建了一个预测类值的向量，我们可以使用gmodels包中的CrossTable()函数将其与实际类值进行比较。将prop.c和prop.r参数设置为FALSE将从表中删除列和行百分比。剩余的百分比（prop.t）表示单元格中记录数占总记录数的比例：

> library(gmodels)
> CrossTable(credit_test$default, credit_pred,
             prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE,
             dnn = c('actual default', 'predicted default')) 

这导致了以下表格：

 | predicted default 
actual default |        no |       yes | Row Total |
---------------|-----------|-----------|-----------|
            no |        56 |         9 |        65 |
               |     0.560 |     0.090 |           |
---------------|-----------|-----------|-----------|
           yes |        24 |        11 |        35 |
               |     0.240 |     0.110 |           |
---------------|-----------|-----------|-----------|
  Column Total |        80 |        20 |       100 |
---------------|-----------|-----------|-----------| 

在测试集的 100 个贷款申请中，我们的模型正确预测了 56 个没有违约，11 个违约，从而实现了 67%的准确率和 33%的错误率。这比它在训练数据上的表现略差，但考虑到模型在未见数据上的表现通常更差，这是可以预料的。此外，请注意，该模型在测试数据中只正确预测了 35 个实际贷款违约中的 11 个，即 31%。不幸的是，这种类型的错误可能是一个代价非常高的错误，因为银行在每次违约中都会损失金钱。让我们看看我们是否可以通过更多的努力来提高结果。

第 5 步 – 提高模型性能

我们模型的错误率可能太高，无法将其部署到实时信用评分应用中。事实上，如果模型对每个测试案例都预测“无违约”，那么它将有 65%的时间是正确的——这个结果几乎与我们的模型一样，但需要付出多得多的努力！仅使用 900 个训练示例来预测贷款违约似乎是一个具有挑战性的问题。

更糟糕的是，我们的模型在识别那些违约的申请人方面表现特别差。幸运的是，有几个简单的方法可以调整 C5.0 算法，这可能会帮助提高模型的整体性能，以及对于更昂贵的错误类型的性能。

提升决策树的准确性

C5.0 算法在 C4.5 算法的基础上进行改进的一种方式是通过添加自适应提升法。这是一个构建许多决策树的过程，这些树对每个示例的最佳类别进行投票。

提升法的想法主要基于 Rob Schapire 和 Yoav Freund 的研究。有关更多信息，请尝试在网络上搜索他们的出版物或他们的教科书 Boosting: Foundations and Algorithms, Cambridge, MA, The MIT Press, 2012。

由于提升法可以更普遍地应用于任何机器学习算法，它将在本书的 第十四章，构建更好的学习者 中更详细地介绍。现在，只需说提升法基于这样的观点：通过结合几个表现不佳的学习者，可以创建一个比任何单个学习者都强大的团队。每个模型都有其独特的一组优势和劣势，可能在某些问题上表现更好或更差。因此，使用具有互补优势和劣势的几个学习者的组合可以显著提高分类器的准确性。

C5.0() 函数使得将提升法添加到我们的决策树变得简单。我们只需添加一个额外的 trials 参数，表示在提升团队中使用的独立决策树的数量。trials 参数设置了一个上限；如果算法认识到额外的试验似乎没有提高准确性，它将停止添加树。我们将从 10 次试验开始，这个数字已经成为事实上的标准，因为研究表明这可以减少测试数据上的错误率大约 25%。

除了新的参数之外，命令与之前相似：

> credit_boost10 <- C5.0(default ~ ., data = credit_train,
                         trials = 10) 

在检查生成的模型时，我们可以看到输出现在指示了提升法的添加：

> credit_boost10 

Number of boosting iterations: 10 
Average tree size: 57.3 

新的输出显示，在 10 次迭代中，我们的树大小缩小了。如果您愿意，可以在命令提示符中键入summary(credit_boost10)来查看所有 10 棵树。请注意，其中一些树，包括为第一次试验构建的树，包含一个或多个子树，如下面的输出摘录中所示，其中一个子树被标记为 [S1]：

dependents > 1: yes (8.8/0.8)
dependents <= 1:
:...years_at_residence <= 1: no (13.4/1.6)
    years_at_residence > 1:
:...age <= 23: yes (11.9/1.6)
    age > 23: [S1] 

注意到那条说如果 age > 23，结果将是 [S1] 的行。为了确定这意味着什么，我们必须将 [S1] 与输出中稍低位置的对应子树匹配，在那里我们看到最终的决策需要几个额外的步骤：

SubTree [S1]
employment_duration in {< 1 year,> 7 years,4 - 7 years,
:                       unemployed}: no (27.7/6.3)
employment_duration = 1 - 4 years:
:...months_loan_duration > 30: yes (7.2)
    months_loan_duration <= 30:
    :...other_credit = bank: yes (2.4)
        other_credit in {none,store}: no (16.6/5.6) 

这样的子树是后剪枝选项（如子树提升和子树替换）的结果，如本章前面提到的。

树的 summary() 输出还显示了该树在训练数据上的性能：

> summary(credit_boost10) 

 (a)   (b)    <-classified as
      ----  ----
       633     2    (a): class no
        17   248    (b): class yes 

分类器在 900 个训练示例上犯了 19 个错误，错误率为 2.1%。这比我们之前提升前的 13.1% 训练错误率有了相当大的改进！然而，我们还需要看看测试数据上是否会有类似的改进。让我们看一下：

> credit_boost_pred10 <- predict(credit_boost10, credit_test)
> CrossTable(credit_test$default, credit_boost_pred10,
             prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE,
             dnn = c('actual default', 'predicted default')) 

结果表格如下：

 | predicted default 
actual default |        no |       yes | Row Total |
---------------|-----------|-----------|-----------|
            no |        58 |         7 |        65 |
               |     0.580 |     0.070 |           |
---------------|-----------|-----------|-----------|
           yes |        19 |        16 |        35 |
               |     0.190 |     0.160 |           |
---------------|-----------|-----------|-----------|
  Column Total |        77 |        23 |       100 |
---------------|-----------|-----------|-----------| 

在这里，我们将提升前的总错误率从 33% 降低到提升模型中的 26%。这可能看起来改进不大，但它离我们预期的 25% 减少并不太远。话虽如此，如果提升可以这么容易地添加，为什么不将其默认应用于每个决策树呢？原因有两个。首先，如果构建一个决策树需要大量的计算时间，构建多个树可能在计算上不切实际。其次，如果训练数据非常嘈杂，那么提升可能根本不会带来改进。尽管如此，如果需要更高的准确性，尝试提升是值得的。

另一方面，模型在识别真实违约方面仍然表现不佳，正确预测的只有 46%（35 个中的 16 个），而简单模型中的正确率是 31%（35 个中的 11 个）。让我们再调查一个选项，看看我们是否可以减少这些代价高昂的错误。

一些错误比其他错误代价更高

给一个可能违约的申请人贷款可能是一个代价高昂的错误。减少错误阴性的一个解决方案可能是拒绝更多的边缘申请人，假设银行从风险贷款中获得的利息远远低于如果钱全部不还所造成的巨大损失。

C5.0 算法允许我们为不同类型的错误分配惩罚，以阻止树犯更多代价高昂的错误。这些惩罚在 成本矩阵 中指定，该矩阵指定每个错误相对于任何其他错误的代价高多少倍。

要开始构建成本矩阵，我们需要首先指定维度。由于预测值和实际值都可以取两个值，是 或 否，我们需要使用两个值组成的两个向量的列表来描述一个 2x2 的矩阵。同时，我们还将命名矩阵维度，以避免以后混淆：

> matrix_dimensions <- list(c("no", "yes"), c("no", "yes"))
> names(matrix_dimensions) <- c("predicted", "actual") 

检查新对象显示我们的维度已经设置正确：

> matrix_dimensions 

$predicted
[1] "no"  "yes"
$actual
[1] "no"  "yes" 

接下来，我们需要通过提供四个值来填充矩阵，以分配各种类型错误的惩罚。由于 R 是从上到下逐列填充矩阵，因此我们需要按照特定顺序提供这些值：

1. 预测为“否”，实际为“否”

2. 预测为“是”，实际为“否”

3. 预测为“否”，实际为“是”

4. 预测为“是”，实际为“是”

假设我们认为贷款违约的成本是银行错过机会的四倍。那么我们的惩罚值可以这样定义：

> error_cost <- matrix(c(0, 1, 4, 0), nrow = 2,
                  dimnames = matrix_dimensions) 

这就产生了以下矩阵：

> error_cost 

 actual
predicted no yes
      no   0   4
      yes  1   0 

如此矩阵定义，当算法正确地将“否”或“是”分类时，没有分配成本，但假阴性有 4 的成本，而假阳性有 1 的成本。为了了解这如何影响分类，让我们使用C5.0()函数的costs参数将其应用于我们的决策树。我们还将使用之前相同的步骤：

> credit_cost <- C5.0(default ~ ., data = credit_train,
                        costs = error_cost)
> credit_cost_pred <- predict(credit_cost, credit_test)
> CrossTable(credit_test$default, credit_cost_pred,
             prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE,
             dnn = c('actual default', 'predicted default')) 

这产生了以下混淆矩阵：

 | predicted default 
actual default |        no |       yes | Row Total |
---------------|-----------|-----------|-----------|
            no |        34 |        31 |        65 |
               |     0.340 |     0.310 |           |
---------------|-----------|-----------|-----------|
           yes |         5 |        30 |        35 |
               |     0.050 |     0.300 |           |
---------------|-----------|-----------|-----------|
  Column Total |        39 |        61 |       100 |
---------------|-----------|-----------|-----------| 

与我们的提升模型相比，这个版本总体上犯的错误更多：这里的错误率为 36%，而在提升模型中为 26%。然而，错误类型非常不同。在先前的模型中，只有 31%和 46%的违约被正确分类，而在本模型中，*30/35=86%*的实际违约被正确预测为违约。这种权衡在减少假阴性错误的同时增加了假阳性错误，如果我们的成本估计准确，这可能是可以接受的。

理解分类规则

分类规则以逻辑 if-else 语句的形式表示知识，为未标记的示例分配一个类别。它们由前件和后件指定，形成一个陈述，即“如果这个发生，那么那个就会发生。”前件包括某些特征值的组合，而后件指定如果规则的条件得到满足，则分配的类别值。一个简单的规则可能声明，“如果计算机正在发出点击声，那么它即将失败。”

规则学习者是决策树学习者的紧密相关的兄弟，通常用于类似类型的任务。像决策树一样，它们可以用于生成未来行动知识的应用，例如：

• 识别导致机械设备硬件故障的条件

• 描述人群群体的关键特征以进行客户细分

• 寻找导致股票市场价格大幅下跌或上涨的条件

规则学习者在与决策树相比有一些明显的差异。与必须通过一系列分支决策的树不同，规则是类似于独立事实陈述的命题。此外，由于以下将要讨论的原因，规则学习者的结果可能比基于相同数据的决策树更简单、更直接、更容易理解。

你可能已经意识到决策树的分支几乎与规则学习算法的 if-else 语句相同，实际上，规则可以从树中生成。那么，为什么还要有一个单独的规则学习算法组呢？继续阅读以发现区分两种方法的细微差别。

规则学习器通常应用于特征主要是或完全是名义性的问题。它们在识别罕见事件方面做得很好，即使罕见事件只发生在特征值之间非常具体的交互中。

分而治之

规则学习分类算法利用一种称为分而治之的启发式方法。这个过程包括识别一个覆盖训练数据中子集的规则，然后从这个分区中分离出剩余的数据。随着规则的添加，更多的数据子集被分离，直到整个数据集都被覆盖，没有更多的例子剩下。尽管分而治之在很多方面与之前提到的分而治之启发式方法相似，但它以微妙的方式有所不同，很快就会变得清楚。

想象分而治之的规则学习过程的一种方式是通过创建越来越具体的规则来深入数据。假设你被要求创建规则来识别动物是否是哺乳动物。你可以将所有动物集合表示为一个大的空间，如下面的图表所示：

图 5.8：一个规则学习算法可能有助于将动物分为哺乳动物和非哺乳动物组

规则学习器首先使用可用的特征来寻找同质群体。例如，使用一个表示物种是否通过陆地、海洋或空中旅行的特征，第一条规则可能建议任何陆生动物都是哺乳动物：

图 5.9：一个潜在的规则将陆地上旅行的动物视为哺乳动物

你注意到这个规则有任何问题吗？如果你是一个动物爱好者，你可能已经意识到青蛙是两栖动物，而不是哺乳动物。因此，我们的规则需要更加具体。让我们进一步深入，建议哺乳动物必须在陆地上行走并且有尾巴：

图 5.10：一个更具体的规则表明，在陆地上行走并且有尾巴的动物是哺乳动物

如前图所示，新规则导致了一组完全由哺乳动物组成的动物子集。因此，哺乳动物的子集可以与其他数据分开，青蛙被返回到剩余动物池中——无意中开玩笑！

可以定义一个额外的规则来区分蝙蝠，这是唯一剩下的哺乳动物。一个可能区分蝙蝠和其他动物的特征是毛皮的存在。使用围绕这个特征构建的规则，我们正确地识别了所有动物：

图 5.11：一条规则表明有毛皮的动物完美地分类了剩余的动物

在这一点上，由于所有训练实例都已分类，规则学习过程将停止。我们总共学习了三条规则：

• 在陆地上行走并带有尾巴的动物是哺乳动物

• 如果动物没有毛皮，它就不是哺乳动物

• 否则，该动物是哺乳动物

之前的例子说明了规则如何逐渐消耗越来越大的数据段，最终对所有实例进行分类。由于规则似乎覆盖了数据的一部分，因此分离和征服算法也被称为覆盖算法，而由此产生的规则被称为覆盖规则。在下一节中，我们将通过检查一个简单的规则学习算法来了解覆盖规则在实际中的应用。然后我们将检查一个更复杂的规则学习器，并将这两个算法应用于一个现实世界的问题。

1R 算法

假设一个电视游戏节目有一个动物隐藏在大型帘子后面。你被要求猜测它是否是哺乳动物，如果猜对了，你将赢得一大笔现金奖金。你没有得到任何关于动物特征的线索，但你知道世界上非常少的动物是哺乳动物。因此，你猜测“非哺乳动物”。你认为你赢得的机会有多大？

当然，选择这个选项最大化了你赢得奖金的机会，因为在假设动物是随机选择的情况下，这是最可能的结果。显然，这个游戏节目有点荒谬，但它展示了最简单的分类器，ZeroR，这是一个不考虑任何特征的规则学习器，实际上它没有学习任何规则（因此得名）。对于每个未标记的例子，无论其特征值如何，它都预测最常见的类别。这个算法在现实世界中的实用性非常有限，除了它为比较其他更复杂的规则学习器提供了一个简单的基线。

1R 算法（也称为One Rule或OneR），通过选择一条规则来改进 ZeroR。尽管这可能看起来过于简单，但它往往比你预期的表现要好。正如实证研究所证明的，这个算法的准确性可以接近许多现实世界任务中更复杂算法的准确性。

想要深入了解 1R 令人惊讶的性能，请参阅非常简单的分类规则在大多数常用数据集上表现良好，Holte, RC, 机器学习，1993，第 11 卷，第 63-91 页。

1R 算法的优缺点如下表所示：

优点	缺点

|

• 生成一个单一、易于理解、人类可读的规则

• 通常表现惊人地好

• 可以作为更复杂算法的基准

|

• 只使用单个特征

• 可能过于简单化

|

该算法的工作方式很简单。对于每个特征，1R 将数据划分为具有相似特征值的组。然后，对于每个段，算法预测多数类别。基于每个特征的规则错误率被计算，选择错误最少的规则作为单一规则。

以下表格显示了这对于我们之前查看的动物数据是如何工作的：

图 5.12：1R 算法选择具有最低误分类率的单一规则

对于通过特征，数据集被划分为三个组：空中、陆地和海洋。空中和海洋组的动物被预测为非哺乳动物，而陆地组的动物被预测为哺乳动物。这导致了两个错误：蝙蝠和青蛙。

有毛皮特征将动物分为两组。有毛皮的动物被预测为哺乳动物，而无毛皮的动物则不是。共有三个错误：猪、大象和犀牛。由于通过特征导致的错误较少，1R 算法会返回以下结果：

• 如果动物通过空中旅行，它就不是哺乳动物

• 如果动物通过陆地旅行，它就是哺乳动物

• 如果动物通过海洋旅行，它就不是哺乳动物

算法在这里停止，因为它找到了最重要的单一规则。

显然，这种规则学习算法可能对于某些任务来说过于基础。您希望医疗诊断系统只考虑单一症状，或者自动驾驶系统只基于单一因素来停止或加速您的汽车吗？对于这些类型的任务，一个更复杂的规则学习者可能更有用。我们将在下一节中了解一个。

RIPPER 算法

早期的规则学习算法存在几个问题。首先，它们以速度慢而闻名，这使得它们对于日益增多的大型数据集无效。其次，它们通常在噪声数据上容易不准确。

1994 年，Johannes Furnkranz 和 Gerhard Widmer 提出了解决这些问题的第一步。他们的增量减少错误剪枝（IREP）算法结合了预剪枝和后剪枝方法，这些方法可以生成非常复杂的规则，并在从完整数据集中分离实例之前对其进行剪枝。尽管这种策略有助于规则学习者的性能，但决策树通常仍然表现更好。

关于 IREP 的更多信息，请参阅 Incremental Reduced Error Pruning, Furnkranz, J and Widmer, G, Proceedings of the 11th International Conference on Machine Learning, 1994, pp. 70-77.

1995 年，威廉·W·科恩引入了重复增量剪枝以产生误差减少（RIPPER）算法，该算法在 IREP 的基础上进行了改进，以生成与决策树性能匹配或超过性能的规则，这使得规则学习者在 1995 年又迈出了新的一步。

关于 RIPPER 的更多详细信息，请参阅Cohen, WW，第 12 届国际机器学习会议论文集，1995 年，第 115-123 页《快速有效的规则归纳》。

如下表概述，RIPPER 的优点和缺点通常与决策树相当。主要好处是它们可能导致一个略微更节俭的模型。

优点	缺点

|

• 生成易于理解、人类可读的规则

• 在大型和噪声数据集上效率高

• 通常，产生的模型比可比较的决策树更简单

|

• 可能会导致看似违反常识或专家知识的规则

• 不适合处理数值数据

• 可能不如更复杂的模型表现得好

|

RIPPER 算法是从几个规则学习算法的迭代中演变而来的，它是一系列用于规则学习的有效启发式方法的拼凑。由于其复杂性，对实现细节的讨论超出了本书的范围。然而，可以将其大致理解为三个步骤的过程：

1. 生长

2. 剪枝

3. 优化

生长阶段使用分离和征服技术贪婪地向规则添加条件，直到它完美地分类数据子集或耗尽用于分割的属性。像决策树一样，使用信息增益标准来识别下一个分割属性。当增加规则的具体性不再减少熵时，规则立即被剪枝。第一步和第二步重复进行，直到达到停止标准，此时使用各种启发式方法对整个规则集进行优化。

与 1R 算法相比，RIPPER 算法可以创建更复杂的规则，因为它可以考虑多个特征。这意味着它可以创建具有多个前件的规则，例如“如果一个动物会飞并且有毛皮，那么它是一种哺乳动物。”这提高了算法建模复杂数据的能力，但就像决策树一样，这也意味着规则可能很快变得难以理解。

分类规则学习者的进化并没有随着 RIPPER 而停止。新的规则学习算法正在迅速提出。文献综述显示，有 IREP++、SLIPPER、TRIPPER 等许多其他算法。

决策树中的规则

分类规则也可以直接从决策树中获得。从叶节点开始，沿着分支回溯到根节点，你可以得到一系列决策。这些决策可以组合成一条规则。以下图显示了如何从预测电影成功的决策树构建规则：

图 5.13：可以通过从根节点到每个叶节点的路径生成规则

沿着从根节点到每个叶节点的路径，规则如下：

1. 如果名人数量少，那么这部电影将是一个 票房炸弹

2. 如果名人数量多且预算高，那么这部电影将是一个 主流热门

3. 如果名人数量多且预算低，那么这部电影将是一个 成功的关键

在下一节中将要解释的原因是，使用决策树生成规则的缺点是生成的规则通常比规则学习算法学习的规则更复杂。决策树使用的划分和征服策略与规则学习器的结果偏差不同。另一方面，有时从树中生成规则在计算上更有效。

在 C50 包中，使用 C5.0() 函数生成模型时，如果您在训练模型时指定 rules = TRUE，则会使用分类规则。

什么使得树和规则是贪婪的？

决策树和规则学习器被称为 贪婪学习器，因为它们基于先到先得的原则使用数据。决策树使用的划分和征服启发式方法以及规则学习器使用的单独和征服启发式方法都试图一次划分一个部分，首先找到最同质的划分，然后是下一个最好的，以此类推，直到所有示例都被分类。

贪婪方法的缺点是贪婪算法不能保证为特定数据集生成最优、最准确或最少数量的规则。通过早期获取低垂的果实，贪婪学习器可能会迅速找到一个适用于数据子集的准确规则；然而，在这样做的同时，学习器可能会错过开发一个更细致的规则集的机会，该规则集在整个数据集上具有更好的整体准确性。然而，如果不使用贪婪方法进行规则学习，对于除了最小的数据集之外的所有数据集，规则学习可能从计算上不可行。

图 5.14：决策树和分类规则学习器都是贪婪算法

虽然树和规则都采用贪婪学习启发式方法，但它们构建规则的方式存在细微差别。也许最好的区分方式是注意，一旦在特征上进行了划分和征服的分割，分割创建的分区可能不会被重新征服，而只会进一步细分。这样，树就永久性地受限于其过去的决策历史。相比之下，一旦单独和征服找到一个规则，任何未覆盖到该规则所有条件的示例可能会被重新征服。

为了说明这种对比，考虑我们之前构建规则学习器来确定动物是否是哺乳动物的情况。规则学习器确定了三条完美分类示例动物的规则：

1. 在陆地上行走并且有尾巴的动物是哺乳动物（熊、猫、狗、大象、猪、兔子、老鼠和犀牛）

2. 如果动物没有皮毛，它就不是哺乳动物（鸟类、鳗鱼、鱼类、青蛙、昆虫和鲨鱼）

3. 否则，动物是哺乳动物（蝙蝠）

相比之下，基于相同数据的决策树可能提出了四条规则来实现相同的完美分类：

1. 如果动物在陆地上行走并且有尾巴，那么它是哺乳动物（熊、猫、狗、大象、猪、兔子、老鼠和犀牛）

2. 如果动物在陆地上行走但没有尾巴，那么它不是哺乳动物（青蛙）

3. 如果动物不在陆地上行走并且有皮毛，那么它就是哺乳动物（蝙蝠）

4. 如果动物不在陆地上行走并且没有皮毛，那么它不是哺乳动物（鸟类、昆虫、鲨鱼、鱼类和鳗鱼）

这两种方法的不同结果与青蛙在“在陆地上行走”决策后发生的情况有关。规则学习器允许青蛙被“没有皮毛”的决策重新征服，而决策树不能修改现有的分区，因此必须将青蛙放入它自己的规则中。

图 5.15：青蛙的处理区分了分而治之与分离和征服启发式方法。后者允许青蛙被后来的规则重新征服。

一方面，因为规则学习器可以重新审视那些被视为但最终未作为先前规则一部分的案件，规则学习器通常比决策树生成的规则更简洁。另一方面，这种数据重用意味着规则学习器的计算成本可能比决策树高一些。

示例 - 使用规则学习器识别有毒蘑菇

每年，许多人因食用有毒野生蘑菇而生病，甚至有些人因此死亡。由于许多蘑菇在外观上非常相似，有时即使是经验丰富的蘑菇采集者也会中毒。

与识别有毒植物，如毒橡树或毒常春藤不同，没有像“三叶草，让它去吧”这样的明确规则来识别野生蘑菇是有毒的还是可食用的。

事情变得更加复杂，许多传统规则，例如“有毒蘑菇颜色鲜艳”，提供了危险或误导性的信息。如果能够有简单、明确和一致的规则来识别有毒蘑菇，它们就能拯救采集者的生命。

规则学习算法的一个优点是它们生成的规则易于理解，这使得它们似乎非常适合这个分类任务。然而，这些规则的有效性将取决于它们的准确性。

第 1 步 – 收集数据

为了识别区分有毒蘑菇的规则，我们将使用卡内基梅隆大学的杰夫·施利默（Jeff Schlimmer）提供的蘑菇数据集。原始数据集可以从 UCI 机器学习仓库免费获取（archive.ics.uci.edu/ml）。

该数据集包括来自 23 种有菌盖蘑菇的 8,124 个蘑菇样本的信息，这些蘑菇在《奥杜邦北美蘑菇野外指南》（1981 年）中列出。在野外指南中，每种蘑菇物种都被标识为“肯定可食用”、“肯定有毒”或“可能有毒，不建议食用”。为了本数据集的目的，后者组被合并到“肯定有毒”组中，形成两个类别：有毒和非有毒。UCI 网站上可用的数据字典描述了蘑菇样本的 22 个特征，包括如菌盖形状、菌盖颜色、气味、菌褶大小和颜色、菌柄形状和栖息地等特征。

本章使用蘑菇数据的略微修改版本。如果你打算跟随示例进行操作，请从本章的 Packt Publishing GitHub 仓库下载mushrooms.csv文件，并将其保存到你的 R 工作目录中。

第 2 步 – 探索和准备数据

我们首先使用read.csv()函数导入分析所需的数据。由于所有 22 个特征和目标类别都是名义的，我们将stringsAsFactors = TRUE设置为利用自动因子转换的优势：

> mushrooms <- read.csv("mushrooms.csv", stringsAsFactors = TRUE) 

str(mushrooms)命令的输出指出，数据包含 8,124 个观测值和 23 个变量，正如数据字典所描述的那样。虽然大多数str()输出并不引人注目，但有一个特征值得提及。你注意到以下行中的veil_type变量有什么特别之处吗？

$ veil_type : Factor w/ 1 level "partial": 1 1 1 1 1 1 ... 

如果你认为一个因子只有一个级别很奇怪，你是正确的。数据字典列出了该特征的两个级别：partial和universal；然而，我们数据中的所有示例都被分类为partial。很可能这个数据元素被错误地编码了。无论如何，由于菌幕类型在样本之间没有变化，它不会为预测提供任何有用的信息。我们将使用以下命令从分析中删除这个变量：

> mushrooms$veil_type <- NULL 

通过将NULL赋值给veil_type向量，R 将消除蘑菇数据框中的该特征。

在深入探讨之前，我们应该快速查看我们数据集中蘑菇类型变量的分布情况：

> table(mushrooms$type) 

 edible poisonous
     4208      3916 

大约 52%的蘑菇样本是可食用的，而 48%是有毒的。

对于这个实验，我们将蘑菇数据中的 8,214 个样本视为所有可能野生蘑菇的完整集合。这是一个重要的假设，因为它意味着我们不需要为测试目的从训练数据中保留一些样本。我们不是试图开发覆盖未预见蘑菇类型的规则；我们只是试图找到准确描述已知蘑菇类型完整集合的规则。因此，我们可以在相同的数据上构建和测试模型。

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

如果我们在这个数据上训练一个假设的 ZeroR 分类器，它会预测什么？由于 ZeroR 忽略了所有特征，只是简单地预测目标的众数，用简单的话说，它的规则会声明“所有蘑菇都是可食用的”。显然，这不是一个非常有用的分类器，因为它会让蘑菇采集者在近一半的蘑菇样本中生病或死亡！我们的规则需要比这个基准做得更好，才能提供可以发布的安全建议。同时，我们需要简单易记的规则。

由于简单的规则仍然可能是有用的，让我们看看一个非常简单的规则学习者在蘑菇数据上的表现。为此，我们将应用 1R 分类器，它将识别对目标类别最有预测性的单个特征，并使用这个特征来构建规则。

我们将使用阿沙芬堡应用科学大学霍格·冯·约恩-迪德里希（Holger von Jouanne-Diedrich）在OneR包中找到的 1R 实现。这是一个相对较新的包，它使用原生 R 代码实现了 1R，以提高速度和易用性。如果您还没有这个包，可以使用install.packages("OneR")命令进行安装，并通过输入library(OneR)来加载。

图 5.16：1R 分类规则语法

与 C5.0 一样，OneR()函数使用 R 公式语法来指定要训练的模型。使用type ~ .公式与OneR()一起允许我们的第一个规则学习者在预测蘑菇类型时考虑蘑菇数据中的所有可能特征：

> mushroom_1R <- OneR(type ~ ., data = mushrooms) 

要检查它创建的规则，我们可以输入分类器对象的名称：

> mushroom_1R 

Call:
OneR.formula(formula = type ~ ., data = mushrooms)
Rules:
If odor = almond   then type = edible
If odor = anise    then type = edible
If odor = creosote then type = poisonous
If odor = fishy    then type = poisonous
If odor = foul     then type = poisonous
If odor = musty    then type = poisonous
If odor = none     then type = edible
If odor = pungent  then type = poisonous
If odor = spicy    then type = poisonous
Accuracy:
8004 of 8124 instances classified correctly (98.52%) 

检查输出，我们看到odor特征被用于规则生成。odor的类别，如almond、anise等，指定了蘑菇是否可能为edible或poisonous的规则。例如，如果蘑菇闻起来有fishy、foul、musty、pungent、spicy或像creosote，那么蘑菇很可能是有毒的。另一方面，有更愉快气味的蘑菇，如almond和anise，以及完全没有气味的蘑菇，被预测为edible。对于蘑菇采集野外指南的目的，这些规则可以总结为一条简单的经验法则：“如果蘑菇闻起来令人不快，那么它很可能是有毒的。”

第 4 步 – 评估模型性能

输出的最后一行指出，规则正确预测了 8,124 个蘑菇样本中的 8,004 个的可食用性，即近 99%。然而，任何不够完美的预测都存在风险，如果模型将有毒蘑菇错误地分类为可食用，可能会导致中毒。

为了确定是否发生了这种情况，让我们检查预测值与实际值之间的混淆矩阵。这需要我们首先生成 1R 模型的预测，然后将预测与实际值进行比较：

> mushroom_1R_pred <- predict(mushroom_1R, mushrooms)
> table(actual = mushrooms$type, predicted = mushroom_1R_pred) 

 predicted
actual      edible poisonous
  edible      4208         0
  poisonous    120      3796 

在这里，我们可以看到我们的规则出错的地方。表格的列表示预测的蘑菇可食用性，而表格的行将 4,208 个实际可食用的蘑菇和 3,916 个实际有毒的蘑菇分开。查看表格，我们可以看到尽管 1R 分类器没有将任何可食用的蘑菇分类为有毒，但它将 120 个有毒蘑菇错误地分类为可食用，这造成了极其危险的错误！

考虑到学习器只利用了一个特征，它表现得相当不错；如果你在寻找蘑菇时避免不愉快的气味，你几乎总是能避免去医院。然而，当涉及到生命安全时，即使是接近完美也不够，更不用说当读者生病时，指南出版商可能不会对诉讼的前景感到高兴。让我们看看我们能否添加一些更多的规则，并开发出一个更好的分类器。

第 5 步 - 提高模型性能

对于一个更复杂的规则学习器，我们将使用JRip()，这是基于 Java 实现的 RIPPER 算法。JRip()函数包含在RWeka包中，它通过 Ian H. Witten 和 Eibe Frank 开发的基于 Java 的 Weka 软件应用，为 R 提供了访问这些机器学习算法的途径。

Weka 是一个流行的开源和功能齐全的图形应用程序，用于执行数据挖掘和机器学习任务——这类工具中最早的之一。有关 Weka 的更多信息，请参阅www.cs.waikato.ac.nz/~ml/weka/。

RWeka包依赖于rJava包，而rJava包本身要求在主机计算机上安装Java 开发工具包（JDK）才能安装。这可以从https://www.java.com/下载，并使用针对您平台的特定说明进行安装。安装 Java 后，使用install.packages("RWeka")命令安装RWeka及其依赖项，然后使用library(RWeka)命令加载RWeka包。

Java 是一套无需付费的编程工具，它允许开发和使用跨平台应用程序，如 Weka。尽管它曾经被许多计算机默认包含，但这已经不再是这样了。不幸的是，它可能很难安装，尤其是在苹果电脑上。如果你遇到麻烦，确保你有最新的 Java 版本。此外，在 Microsoft Windows 上，你可能需要正确设置你的环境变量，如 JAVA_HOME，并检查你的 PATH 设置（在网上搜索详细信息）。在 macOS 或 Linux 电脑上，你也可以尝试从终端窗口执行 R CMD javareconf，然后使用 R 命令 install.packages("rJava", type = "source") 从源安装 R 包 rJava。如果所有其他方法都失败了，你可以尝试一个免费的 Posit Cloud 账户 (posit.cloud/)，它提供了一个已经安装了 Java 的 RStudio 环境。

在安装了 rJava 和 RWeka 之后，训练 JRip() 模型的过程与训练 OneR() 模型的过程非常相似，如下面的语法框所示。这是 R 公式接口的一个令人愉悦的好处：语法在算法之间是一致的，这使得比较各种模型变得简单。

图 5.17：RIPPER 分类规则语法

让我们像训练 OneR() 一样训练 JRip() 规则学习器，允许它在所有可用特征中找到规则：

> mushroom_JRip <- JRip(type ~ ., data = mushrooms) 

要检查规则，请输入分类器的名称：

> mushroom_JRip 

JRIP rules:
===========
(odor = foul) => type=poisonous (2160.0/0.0)
(gill_size = narrow) and (gill_color = buff)
  => type=poisonous (1152.0/0.0)
(gill_size = narrow) and (odor = pungent)
  => type=poisonous (256.0/0.0)
(odor = creosote) => type=poisonous (192.0/0.0)
(spore_print_color = green) => type=poisonous (72.0/0.0)
(stalk_surface_below_ring = scaly)
  and (stalk_surface_above_ring = silky)
    => type=poisonous (68.0/0.0)
(habitat = leaves) and (gill_attachment = free)
  and (population = clustered)
  => type=poisonous (16.0/0.0)
=> type=edible (4208.0/0.0)
Number of Rules : 8 

JRip() 分类器从蘑菇数据中学习了总共八条规则。

读取这些规则的一个简单方法是将它们视为一系列的 if-else 语句，类似于编程逻辑。前三个规则可以表示为：

• 如果气味难闻，那么蘑菇类型是有毒的

• 如果菌盖大小狭窄且菌盖颜色为米色，那么蘑菇类型是有毒的

• 如果菌盖大小狭窄且气味刺鼻，那么蘑菇类型是有毒的

最后，第八条规则意味着任何前七条规则未涵盖的蘑菇样本都是可食用的。按照我们的编程逻辑，这可以读作：

• 否则，蘑菇是可食用的

每条规则旁边的数字表示该规则覆盖的实例数量和误分类实例的数量。值得注意的是，使用这八条规则没有误分类的蘑菇样本。因此，最后一条规则覆盖的实例数量正好等于数据中可食用蘑菇的数量（N = 4,208）。

下图提供了一个如何将规则应用于蘑菇数据的粗略说明。如果你想象这个大椭圆形包含所有蘑菇种类，规则学习器会识别出将同质段与较大群体区分开来的特征，或特征集。首先，算法发现了一群独特的有毒蘑菇，它们以其恶臭而闻名。接下来，它发现了更小、更具体的有毒蘑菇群体。通过为每种有毒蘑菇识别覆盖规则，所有剩余的蘑菇都是可食用的。

感谢大自然母亲，每种蘑菇的独特性足够让分类器达到 100% 的准确性。

图 5.18：一个复杂的规则学习算法识别了规则，完美地覆盖了所有类型的有毒蘑菇

摘要

本章介绍了两种使用所谓“贪婪”算法根据特征值对数据进行划分的分类方法。决策树使用分而治之的策略来创建类似流程图的结构，而规则学习器则分离并征服数据以识别逻辑的 if-else 规则。这两种方法产生的模型可以在没有统计背景的情况下进行解释。

一种流行且高度可配置的决策树算法是 C5.0。我们使用 C5.0 算法创建了一棵树，用于预测贷款申请人是否会违约。通过使用提升和成本敏感的错误选项，我们提高了准确性，并避免了可能给银行带来更多金钱损失的风险贷款。

我们还使用了两个规则学习器，1R 和 RIPPER，来开发识别有毒蘑菇的规则。1R 算法使用单个特征实现了 99% 的准确性，用于识别可能致命的蘑菇样本。另一方面，由更复杂的 RIPPER 算法生成的八条规则正确地识别了每种蘑菇的可食用性。

这仅仅触及了如何使用树和规则的一角。下一章，第六章，预测数值数据 – 回归方法，描述了被称为回归树和模型树的技巧，这些技巧使用决策树进行数值预测而不是分类。在第八章，寻找模式 – 使用关联规则进行市场篮子分析中，我们将看到关联规则——分类规则的近亲——如何被用来识别交易数据中的物品组。最后，在第十四章，构建更好的学习者中，我们将发现如何通过将决策树组合在一个称为随机森林的模型中来提高决策树的表现，此外还有其他依赖于决策树的先进建模技术。

加入我们书籍的 Discord 空间

加入我们的 Discord 社区，与志同道合的人相聚，并和超过 4000 人在以下链接处一起学习：

packt.link/r

第六章：预测数值数据 – 回归方法

数学关系帮助我们理解日常生活的许多方面。例如，体重是个人摄入卡路里的函数；收入通常与教育和工作经验相关；而民意调查数字有助于估计总统候选人连任的机会。

当用数字来表述这些模式时，我们能够获得额外的清晰度。例如，每天额外摄入 250 千卡路里可能会导致每月增加近 1 公斤的体重；每增加一年的工作经验可能会使年薪增加 1000 美元；当经济强劲时，总统更有可能获得连任。显然，这些方程并不完美地适用于每一种情况，但我们预期它们在大多数时候是相当准确的。

本章通过超越之前介绍的分类方法，并引入用于估计数值数据中关系的技巧，扩展了我们的机器学习工具箱。在考察几个现实世界的数值预测任务时，你将学习：

• 回归中使用的基本统计原理，这是一种模拟数值关系大小和强度的技术

• 如何为回归分析准备数据，估计和解释回归模型，以及应用回归变体，如广义线性模型

• 一对称为回归树和模型树的混合技术，这些技术将决策树分类器适应于数值预测任务

基于统计学领域的丰富研究成果，本章所使用的方法在数学方面比之前介绍的方法更为复杂，但请放心！即使你的代数技能有些生疏，R 语言会帮你处理繁重的工作。

理解回归

回归涉及指定单个数值因变量（要预测的值）与一个或多个数值自变量（预测因子）之间的关系。正如其名称所暗示的，因变量依赖于自变量或变量的值。回归的最简单形式假设自变量和因变量之间的关系遵循一条直线。

术语“回归”用来描述将线拟合到数据的过程，其起源可以追溯到 19 世纪末弗朗西斯·高尔顿爵士对遗传学的研究。他发现，身高极矮或极高的父亲往往会有身高更接近平均值的儿子。他把这种现象称为“回归到平均值”。

你可能还记得从基础代数中，直线可以用类似于 y = a + bx 的斜率截距形式来定义。在这种形式中，字母 y 表示因变量，x 表示自变量。斜率项 b 指定了直线在 x 增加时上升的量。正值定义了向上倾斜的直线，而负值定义了向下倾斜的直线。项 a 被称为截距，因为它指定了直线与垂直 y 轴相交或截取的点。它表示当 x = 0 时 y 的值。

图 6.1：具有各种斜率和截距的直线示例

回归方程使用类似的斜率截距格式来模拟数据。机器的任务是确定* a* 和 b 的值，使得指定的直线最能将提供的 x 值与 y 的值联系起来。

可能并不总是存在一组完美的 a 和 b 参数来完美地关联这些值，因此机器还必须有一些方法来量化误差范围并选择最佳拟合。我们将在稍后深入讨论这个问题。

回归分析被用于各种任务——它几乎肯定是应用最广泛的机器学习方法。它可以用于解释过去并预测未来，并且可以应用于几乎任何任务。一些具体的用例包括：

• 检验人口和个体通过其测量的特征如何变化，这在经济学、社会学、心理学、物理学和生态学等科学研究中

• 量化事件与其响应之间的因果关系，例如在临床试验、工程安全测试或市场研究中

• 识别可以使用已知标准来预测未来行为的模式，例如用于预测保险索赔、自然灾害损失、选举结果和犯罪率

回归方法也用于统计假设检验，这决定了在观察数据的基础上，一个前提是否可能是真实的或错误的。回归模型对关系强度和一致性的估计提供了可用于评估观察结果是否仅由偶然性引起的信息。

假设检验非常微妙，超出了机器学习的范围。如果你对这个主题感兴趣，一本入门统计学教科书是一个好的起点，例如，直观入门统计学，沃尔夫，D. A. 和施奈德，G.，斯普林格，2017。

回归分析并不等同于单个算法。相反，它是一个涵盖许多方法的术语，这些方法可以适应几乎任何机器学习任务。如果你只能选择一个机器学习方法来研究，回归将是一个不错的选择。一个人可能可以将整个职业生涯都投入到这个领域，也许仍然有很多东西要学习。

在本章中，我们将从最基本的线性回归模型开始——那些使用直线的模型。只有一个自变量的情况被称为简单线性回归。有两个或更多自变量的情况被称为多元线性回归，或简称多元回归。这两种技术都假设有一个单一的因变量，它在连续尺度上被测量。

回归还可以用于其他类型的因变量，甚至可以用于某些分类任务。例如，逻辑回归用于建模二元分类结果，而泊松回归（以法国数学家西莫恩·泊松的名字命名）用于建模整数计数数据。被称为多项式逻辑回归的方法用于建模分类结果，因此可以用于分类。

这些专门的回归方法属于广义线性模型（GLMs）类别，它们将传统回归模型的直线调整为允许建模其他形式的数据。这些将在本章后面进行描述。

由于类似的统计原理适用于所有回归方法，一旦你理解了线性情况，了解其他变体就变得简单直接。我们将从简单线性回归的基本情况开始。尽管名称上看似简单，但这种方法并不简单到无法解决复杂问题。在下一节中，我们将看到简单线性回归模型的使用如何可能避免一场悲剧性的工程灾难。

简单线性回归

1986 年 1 月 28 日，美国航天飞机“挑战者”号的七名机组人员在火箭助推器故障导致灾难性解体时丧生。在事故发生后，专家们迅速将发射温度视为潜在的罪魁祸首。负责密封火箭接头的橡胶 O 形圈从未在 40°F（4°C）以下进行过测试！F (4C)，而发射当天的天气异常寒冷，低于冰点。

借助事后诸葛的优势，这次事故已成为数据分析可视化的重要性案例研究。尽管不清楚火箭工程师和决策者在发射前可以获得哪些信息，但不可否认的是，更好的数据，如果被谨慎使用，很可能避免这场灾难。

本节的分析基于*《航天飞机风险分析：挑战者号之前的故障预测》，Dalal, S. R., Fowlkes, E. B., 和 Hoadley, B.，美国统计学会杂志，1989 年，第 84 卷，第 945-957 页中呈现的数据。关于数据如何可能改变结果的一个观点，请参阅《视觉解释：图像与数量，证据与叙事》，Tufte, E. R.，Cheshire, C. T.：Graphics Press，1997 年*。对于相反的观点，请参阅*《表现与误表现：图夫特与挑战者号上的莫顿·索尔科工程师》，Robison, W.，Boisjoly, R.，Hoeker, D.，和 Young, S.，《科学和工程伦理》，2002 年，第 8 卷，第 59-81 页*。

火箭工程师几乎肯定知道低温会使组件变得更加脆弱，密封性能降低，这会导致危险燃料泄漏的概率增加。然而，鉴于继续发射的政治压力，他们需要数据来支持这一假设。一个展示温度与 O 形圈故障之间联系的回归模型，并且能够根据预期的发射温度预测故障概率，可能会非常有帮助。

为了构建回归模型，科学家可能使用了在 23 次之前成功的航天飞机发射中记录的发射温度和组件故障的数据。组件故障表明两种类型的问题之一。第一种问题，称为侵蚀，发生在过度的热量烧毁 O 形圈时。第二种问题，称为吹过，发生在热气体通过或“吹过”密封不良的 O 形圈时。由于航天飞机共有六个主要 O 形圈，每次飞行可能发生多达六个故障。尽管火箭可以生存一个或多个故障事件或被一个故障摧毁，但每个额外的故障都会增加灾难性故障的概率。以下散点图显示了之前 23 次发射检测到的初级 O 形圈故障，与发射温度的对比：

图 6.2：航天飞机 O 形圈故障与发射温度的可视化

检查图表，存在一个明显的趋势：在较高温度下进行的发射往往有较少的 O 形圈故障事件。此外，最冷的发射（53°F）发生了两个故障事件，这个水平只在另一次发射中达到。考虑到这些信息，挑战者号计划在比这低 20 多度的条件下发射似乎令人担忧。但他们应该有多担心？为了回答这个问题，我们可以转向简单线性回归。

简单线性回归模型定义了因变量与单个自变量预测变量之间的关系，使用以下形式的方程定义的线：

除了希腊字母外，这个方程几乎与之前描述的斜截式相同。截距 （alpha）描述了直线与 y 轴的交点，而斜率 （beta）描述了在 x 增加时 y 的变化。对于航天飞机发射数据，斜率将告诉我们发射温度每增加一度，O 形圈故障预期的变化。

在统计学领域，希腊字母通常用于表示统计函数的参数变量。因此，进行回归分析涉及找到 和 的参数估计。alpha 和 beta 的参数估计通常用 a 和 b 表示，尽管你可能发现一些术语和符号被交替使用。

假设我们知道航天飞机发射数据方程中估计的回归参数是 a = 3.70 和 b = -0.048。因此，完整的线性方程是 y = 3.70 – 0.048x。暂时忽略这些数字是如何得到的，我们可以像这样在散点图上绘制这条线：

图 6.3：一个回归线，用于模拟压力事件与发射温度之间的关系

如线所示，在 60 华氏度时，我们预测不到一个 O 形圈的故障事件。在 50 华氏度时，我们预计大约有 1.3 次故障。如果我们使用该模型外推到 31 度——挑战者号预测的温度——我们预计大约会有 3.70 - 0.048 * 31 = 2.21 次 O 形圈的故障事件。

假设每个 O 形圈故障同样可能引起灾难性的燃料泄漏，这意味着挑战者号在 31 度发射的风险几乎是 60 度典型发射的三倍，比 70 度发射的风险超过八倍。

注意到直线并没有精确地穿过每个数据点。相反，它在大约均匀地穿过数据，一些预测值低于或高于直线。在下一节中，我们将学习为什么选择这条特定的线。

普通最小二乘估计

为了确定 和 的最佳估计值，使用了一种称为普通最小二乘法（OLS）的估计方法。在 OLS 回归中，斜率和截距被选择以最小化平方误差和（SSE）。误差，也称为残差，是预测的 y 值与实际 y 值之间的垂直距离。由于误差可能是高估或低估，它们可以是正数或负数；平方它们使得误差无论方向如何都是正数。以下图表展示了几个点的残差：

图 6.4：回归线的预测值与实际值之间的差异由残差量决定

用数学术语来说，OLS 回归的目标可以表达为最小化以下方程的任务：

用简单的话说，这个方程定义e（误差）为实际y值与预测y值之间的差异。误差值被平方以消除负值，并在数据中的所有点上进行求和。

在y项上方的撇号（^）是统计符号中常用的一个特征。它表示该项是对真实y值的估计。这被称为y的估计值。

a的解取决于b的值。可以使用以下公式获得：

要理解这些方程，你需要了解另一部分统计符号。出现在x和y项上方的水平横线表示x或y的均值值。这被称为x横或y横。

虽然证明超出了本书的范围，但可以使用微积分证明，导致最小平方误差的b的值是：

如果我们将这个方程分解为其组成部分，我们可以稍微简化它。b的分母应该看起来很熟悉；它与x的方差非常相似，表示为 Var(x)。正如我们在第二章，管理和理解数据中学到的，方差涉及找到x的均值与均值的平均平方偏差。这可以表示为：

分子涉及将每个数据点的偏差（与均值x值的偏差）乘以该点偏离均值y值的偏差之和。这与x和y的协方差函数类似，表示为 Cov(x, y)。协方差公式是：

如果我们将协方差函数除以方差函数，分子和分母中的n项会相互抵消，我们可以将b的公式重写为：

给出这种重述后，使用内置的 R 函数很容易计算出b的值。让我们将它们应用于航天飞机发射数据来估计回归线。

如果你想跟随这些示例，请从 Packt Publishing 网站下载challenger.csv文件，并使用launch <- read.csv("challenger.csv")命令将其加载到数据框中。

如果航天飞机数据存储在名为launch的数据框中，自变量x命名为temperature，因变量y命名为distress_ct，则可以使用 R 函数cov()和var()来估计b：

> b <- cov(launch$temperature, launch$distress_ct) /
         var(launch$temperature)
> b 

[1] -0.04753968 

然后，我们可以通过使用计算出的b值并应用mean()函数来估计a：

> a <- mean(launch$distress_ct) - b * mean(launch$temperature)
> a 

[1] 3.698413 

手动估计回归方程显然不是最佳选择，因此 R 预期地提供了一个用于自动拟合回归模型的功能。我们很快就会使用这个功能。在此之前，通过首先学习一种衡量线性关系强度的方法，来扩展你对回归模型拟合的理解是很重要的。此外，你很快就会学习如何将多元线性回归应用于具有多个自变量的问题。

相关性

两个变量之间的相关性是一个数字，表示它们之间的关系有多接近一条直线。如果没有额外的限定，相关性通常指的是皮尔逊相关系数，这是 20 世纪数学家卡尔·皮尔逊开发的。相关性的范围在 -1 到 +1 之间。最大值和最小值表示完全的线性关系，而接近零的相关性表示没有线性关系。

以下公式定义了皮尔逊相关系数：

这里引入了更多的希腊符号：第一个符号（看起来像小写的 p）是 rho，它用来表示皮尔逊相关统计量。看起来像 q 字符逆时针旋转的符号是希腊字母小写 sigma，它们表示 x 或 y 的标准差。

使用这个公式，我们可以计算发射温度和 O 型圈故障事件数之间的相关性。回想一下，协方差函数是 cov()，标准差函数是 sd()。我们将结果存储在 r 中，这是一个常用来表示估计相关性的字母：

> r <- cov(launch$temperature, launch$distress_ct) /
         (sd(launch$temperature) * sd(launch$distress_ct))
> r 

[1] -0.5111264 

或者，我们可以使用 cor() 相关函数得到相同的结果：

> cor(launch$temperature, launch$distress_ct) 

[1] -0.5111264 

温度和故障 O 型圈数量之间的相关系数是 -0.51。负相关性表明温度的增加与故障 O 型圈数量的减少有关。对于研究 O 型圈数据的 NASA 工程师来说，这将是一个非常明确的指标，表明低温发射可能存在问题。相关性还告诉我们温度和 O 型圈故障之间的相对强度。因为 -0.51 是最大负相关系数 -1 的一半，这意味着存在中等强度的负线性关联。

在解释相关性强度方面，有许多经验法则。一种方法是将 0.1 到 0.3 之间的值标记为“弱”；0.3 到 0.5 之间的范围标记为“中等”；而高于 0.5 的值标记为“强”（这些也适用于类似范围的负相关性）。然而，这些阈值对于某些目的可能过于严格或过于宽松。通常，相关性必须在特定背景下进行解释。

对于涉及人类的数据，0.5 的相关性可能被认为非常高；对于由机械过程生成的数据，0.5 的相关性可能非常弱。

你可能听说过“相关性不等于因果关系”这个说法。这源于这样一个事实，即相关性仅描述了一对变量之间的关联，但可能存在其他未考虑的解释，这些解释可能是观察到的关系的责任。例如，寿命与每天观看电影的时间之间可能存在强烈的关联，但在医生建议我们所有人都看更多电影之前，我们需要排除另一个解释：年轻人看更多的电影，而年轻人（总的来说）不太可能死亡。

测量两个变量之间的相关性为我们提供了一种快速检查独立变量和因变量之间线性关系的方法。随着我们开始定义具有更多预测因子的回归模型，这一点将变得越来越重要。

多元线性回归

大多数现实世界的分析都有多个自变量。因此，你很可能会在大多数数值预测任务中使用多元线性回归。多元线性回归的优缺点如下表所示：

优点	缺点

|

• 到目前为止，建模数值数据最常见的方法

• 可以适应几乎任何建模任务

• 提供了特征和结果之间关系的大小和强度的估计

|

• 对数据有强烈的假设

• 模型的形式必须由用户事先指定

• 无法处理缺失数据

• 只适用于数值特征，因此分类数据需要额外的准备

• 需要一些统计学知识来理解模型

|

我们可以将多元回归视为简单线性回归的扩展。在两种情况下，目标都是相似的——找到斜率系数的值，以最小化线性方程的预测误差。关键的区别是，对于额外的自变量，有额外的项。

多元回归模型具有以下方程的形式。因变量 y 被指定为截距项 与每个 i 个特征对应的估计值 与 x 变量的乘积之和。这里增加了一个误差项 （用希腊字母 epsilon 表示），作为提醒预测并不完美。这代表了之前提到的残差项：

让我们暂时考虑一下估计的回归参数的解释。你会注意到，在前面的方程中，为每个特征提供了一个系数。这允许每个特征对 y 的值有单独的估计影响。换句话说，当特征 ![img/B17290_06_023.png] 每增加一个单位时，y 的变化量为 ![img/B17290_06_022.png]。因此，截距 ![img/B17290_06_005.png] 是当独立变量都为零时 y 的期望值。

由于截距项 ![img/B17290_06_005.png] 真的与其他任何回归参数没有区别，它有时也用 ![img/B17290_06_024.png] （发音为 beta naught）表示，如下面的公式所示：

![img/B17290_06_025.png]

就像之前一样，截距与任何独立变量 x 都无关。然而，由于以下原因将在短时间内变得清楚，想象 ![img/B17290_06_024.png] 被乘以一个项 x[0] 有助于理解。我们将 x[0] 分配为一个常量，其值为 1：

![img/B17290_06_027.png]

为了估计回归参数，必须将因变量 y 的每个观测值与独立变量 x 的观测值通过前面形式的回归方程联系起来。以下图是多个回归任务设置的图形表示：

图 6.5：多重回归试图找到![img/B17290_06_006.png]值，这些值将 X 值与 Y 值相关联，同时最小化![img/B17290_06_020.png]

上述图中所展示的许多行和列的数据可以用加粗的矩阵符号来压缩表示，以表明每个项代表多个值。以这种方式简化后，公式如下：

![img/B17290_06_030.png]

在矩阵表示法中，因变量是一个向量，Y，其中每一行代表一个示例。独立变量被组合成一个矩阵，X，其中每一列代表一个特征，还有一个额外的截距列，值为 1。每一列都有每一行的示例。回归系数 ![img/B17290_06_031.png] 和残差误差 ![img/_eqn_032.png] 现在也是向量。

目标现在是要求解![img/B17290_06_031.png]，这是使预测值和实际Y值之间平方误差和最小的回归系数向量。找到最优解需要使用矩阵代数；因此，推导需要比本文中提供的更仔细的注意。然而，如果你愿意相信他人的工作，向量![img/B17290_06_031.png]的最佳估计可以计算如下：

![img/B17290_06_034.png]

这个解决方案使用了一对矩阵运算：T 表示矩阵 X 的转置，而负指数表示 矩阵逆。利用 R 的内置矩阵运算，我们可以实现一个简单的多元回归学习器。让我们将这个公式应用到挑战者号发射数据上。

如果你对前面的矩阵运算不熟悉，Wolfram MathWorld 的转置页面 (mathworld.wolfram.com/Transpose.html) 和矩阵逆页面 (mathworld.wolfram.com/MatrixInverse.html) 提供了全面的介绍，即使没有高级数学背景也能理解。

使用以下代码，我们可以创建一个名为 reg() 的基本回归函数，它接受参数 y 和 x，并返回一个估计的贝塔系数向量：

> reg <- function(y, x) {
    x <- as.matrix(x)
    x <- cbind(Intercept = 1, x)
    b <- solve(t(x) %*% x) %*% t(x) %*% y
    colnames(b) <- "estimate"
    print(b)
  } 

这里创建的 reg() 函数使用了几个我们之前没有用过的 R 命令。首先，由于我们将使用该函数与数据框的列集一起使用，as.matrix() 函数将数据框转换为矩阵形式。

接下来，cbind() 函数将一个额外的列绑定到 x 矩阵上；命令 Intercept = 1 指示 R 将新列命名为 Intercept 并用重复的 1 值填充该列。然后，对 x 和 y 对象执行一系列矩阵运算：

• solve() 取矩阵的逆

• t() 用于矩阵的转置

• %*% 用于两个矩阵的乘法

通过将这些 R 函数组合起来，我们的函数将返回一个向量 b，其中包含与 x 到 y 相关的线性模型的估计参数。函数中的最后两行给 b 向量命名并在屏幕上打印结果。

让我们将这个函数应用到航天飞机发射数据上。如下面的代码所示，数据集包括三个特征和灾情计数（distress_ct），这是我们感兴趣的结果：

> str(launch) 

'data.frame':    23 obs. of  4 variables:
 $ distress_ct         : int  0 1 0 0 0 0 0 0 1 1 ...
 $ temperature         : int  66 70 69 68 67 72 73 70 57 63 ...
 $ field_check_pressure: int  50 50 50 50 50 50 100 100 200 200 ...
 $ flight_num          : int  1 2 3 4 5 6 7 8 9 10 ... 

我们可以通过比较其对于之前找到的 O 型环故障与温度的简单线性回归模型的结果来确认我们的函数是否正确工作，该模型具有参数 a = 3.70 和 b = -0.048。由于温度位于发射数据的第二列，我们可以如下运行 reg() 函数：

> reg(y = launch$distress_ct, x = launch[2]) 

 estimate
Intercept    3.69841270
temperature -0.04753968 

这些值与我们之前的结果完全匹配，因此让我们使用这个函数来构建一个多元回归模型。我们将像之前一样应用它，但这次我们将指定 x 参数的第二到第四列以添加两个额外的预测因子：

> reg(y = launch$distress_ct, x = launch[2:4]) 

 estimate
Intercept             3.527093383
temperature          -0.051385940
field_check_pressure  0.001757009
flight_num            0.014292843 

该模型使用温度、现场检查压力和发射 ID 号来预测 O 形圈故障事件的数量。值得注意的是，加入这两个新的预测因子并没有改变我们从简单线性回归模型中得到的结果。正如之前一样，温度变量的系数为负，这表明随着温度的升高，预期的 O 形圈事件数量会减少。这种影响的大小也大致相同：发射温度每增加一度，预计的故障事件将减少大约 0.05 次。

这两个新的预测因子也对预测的故障事件有贡献。现场检查压力是指在发射前测试中对 O 形圈施加的压力量。虽然检查压力最初是 50 磅力，但在某些发射中提高到 100 和 200 磅力，这导致一些人认为这可能是 O 形圈侵蚀的原因。系数为正，但很小，至少为这个假设提供了一些证据。飞行次数代表了航天飞机的年龄。每次飞行，它都会变老，部件可能会更脆弱或更容易损坏。飞行次数与故障计数之间的小的正相关可能反映了这一事实。

总体而言，我们对航天飞机数据的回顾性分析表明，考虑到天气条件，挑战者号的发射风险很高。也许如果工程师们事先应用线性回归，灾难可能就可以避免。当然，当时的情况现实，以及涉及的政治影响，肯定没有现在回望时看起来那么简单。

广义线性模型和逻辑回归

如同在挑战者号航天飞机发射数据分析中所展示的那样，标准线性回归是用于建模数值结果与一个或多个预测因子之间关系的有用方法。回归能够经受时间的考验，这并不奇怪。即使在一百多年后，它仍然是我们工具箱中最重要的技术之一，尽管它并不比找到最佳直线来拟合数据更复杂。

然而，并非每个问题都适合用一条线来建模，而且，回归模型在许多实际任务中违反了统计假设。即使是挑战者号的数据，对于线性回归来说也不够理想，因为它违反了回归假设，即目标变量是在连续尺度上测量的。由于 O 形圈的故障次数只能取可数值，因此模型预测出恰好 2.21 次故障事件，而不是两个或三个，这是没有意义的。

对于建模计数值、分类或二元结果，以及其他目标变量不是正态分布的连续变量的情况，标准线性回归并不是最佳工具——尽管许多人仍然将这些类型的问题应用于它，并且它通常表现得相当出色。

为了解决这些不足，可以使用名为 GLM 的适当名称来适应其他用例，该模型最早由统计学家 John Nelder 和 Robert Wedderburn 于 1972 年描述。GLM 放宽了传统回归模型的两个假设。首先，它允许目标变量是非正态分布的、非连续变量。其次，它允许目标变量的方差与其均值相关。前者为建模分类数据或计数数据，甚至预测值范围有限的情况（例如，落在 0 到 1 之间的概率值）打开了大门。后者允许模型更好地拟合预测变量与预测以非线性方式相关的情况，例如指数增长，其中时间的增加导致结果的增加越来越大。

要阅读关于广义线性模型（GLM）的原始出版物，请参阅Nelder, J. A. 和 Wedderburn, T. W. M.，皇家统计学会杂志，1972 年，第 135 卷，第 370-384 页。对于更温和、非数学性的介绍，请参阅Dunteman, G. H. 和 Ho, M. H. R.，社会科学定量应用，2006 年，第 145 卷。

这两种对线性回归的推广反映在任何 GLM 的两个关键组件中：

1. 家族指的是目标特征的分布，必须从指数分布族的成员中选择，该族包括正态高斯分布以及其他如泊松、二项和伽马分布。所选分布可以是离散的或连续的，并且可以跨越不同的值范围，例如仅正数或仅介于零和一之间的值。

2. 链接函数将预测变量与目标变量之间的关系转换为可以使用线性方程进行建模的形式，尽管原始关系是非线性的。始终存在一个标准链接函数，它由所选的家族确定并默认使用，但在某些情况下，可以选择不同的链接来改变模型的理解方式或获得更好的模型拟合。

调整家族和链接函数使广义线性模型（GLM）方法具有极大的灵活性，以适应许多不同的实际应用场景，并符合目标变量的自然分布。了解使用哪种组合取决于模型的应用方式以及目标变量的理论分布。详细了解这些因素需要了解指数族中的各种分布以及统计学理论背景。幸运的是，大多数广义线性模型（GLM）的应用符合一些常见的家族和链接组合，这些组合在下面的表中列出：

家族	规范链接函数	目标范围	注意事项和应用
高斯（正态）	恒等	-![img/B17290_06_036.png]到![img/B17290_06_036.png]	用于线性响应建模；将广义线性模型（GLM）简化为标准线性回归。
泊松	对数	整数 0 到![img/B17290_06_036.png]	称为泊松回归；通过估计事件发生的频率来建模事件发生的次数（如 O 型圈故障的总数）。
二项式	对数几率	0 到 1	称为逻辑回归；通过估计结果发生的概率来建模二元结果（如是否有 O 型圈故障）。
伽马	负逆	0 到![img/B17290_06_036.png]	建模右偏斜数据的一种可能性；可用于建模事件发生的时间（如 O 型圈故障的秒数）或成本数据（如汽车事故的保险索赔成本）。
多项式	对数几率	K个类别中的 1 个	称为多项式逻辑回归；通过估计示例落在每个类别中的概率来建模分类结果（如成功、失败或中止的航天飞机发射）。通常使用专门的软件包而不是广义线性模型（GLM）函数来帮助解释。

由于解释广义线性模型（GLM）的细微差别，要熟练应用一个模型需要大量的实践和仔细的研究，而且很少有人能声称自己是所有这些模型的专家。整本教科书都致力于每种广义线性模型（GLM）变体。幸运的是，在机器学习的领域，解释和理解的重要性不如能够将正确的广义线性模型（GLM）形式应用于实际问题并产生有用的预测。虽然本章不能涵盖列出的每种方法，但关键细节的介绍将允许你后来追求与你自己的工作最相关的广义线性模型（GLM）变体。

从表中列出的最简单变体开始，标准线性回归可以被视为一种特殊的 GLM，它使用高斯族和恒等链接函数。恒等链接意味着目标 y 与预测变量 x[i] 之间的关系没有进行任何转换。因此，与标准回归一样，一个估计的回归参数 可以相当简单地解释为在 x[i] 增加 1 个单位时 y 的增加量，假设所有其他因素都保持不变。

使用其他链接函数的 GLM 形式并不容易解释，要完全理解单个预测变量的影响需要更加仔细的分析。这是因为回归参数必须被解释为在 x 增加 1 个单位时 y 的增加量，但这是在通过链接函数转换之后。

例如，使用 对数链接 函数来模拟事件预期计数的泊松族，通过自然对数将 y 与预测变量 x[i] 相关联；因此， 对 y 的加性效应在响应变量的原始尺度上变成了乘性效应。这是因为利用对数的性质，我们知道 ，在指数化以消除对数后，这变成了 。

由于这种乘性影响，参数估计被理解为相对增加率，而不是像线性回归中那样是 y 的绝对增加量。

为了在实践中看到这一点，假设我们构建了一个关于 O 型圈故障数与发射温度的泊松回归模型。如果 x[1] 是温度，并且估计的 = -0.103，那么我们可以确定，在发射时每增加 1 度温度，平均 O 型圈故障数会减少约 9.8%。这是因为 exp(-0.103) = 0.902，或者说每度故障的 90.2%，这意味着我们预计每增加 1 度温度，故障数会减少 9.8%。将此应用于挑战者号发射时的 36 华氏度温度，我们可以外推，如果发射温度再高 17 度（53 华氏度是之前的最低发射温度），那么预期的故障数将大约是 (0.902)¹⁷ = 17.2%，相当于减少 82.8%。

在航天飞机数据的背景下，假设我们为预测发射过程中是否会发生一个或多个 O 形圈故障的二元分类任务构建了一个逻辑回归模型。一个不会改变 O 形圈故障概率的因素将保持赔率为 1:1（50-50 的概率），这转化为对数赔率 log(0.5 / (1 - 0.5)) = 0 以及对于这个特征的估计回归系数 = 0。找到赔率比 exp(0) = 1 表明，无论这个因素的价值如何，赔率保持不变。现在，假设像温度这样的因素降低了结果发生的可能性，并且在这个以 x[1] 作为温度的逻辑回归模型中，那么估计的 = -0.232。通过指数化这个值，我们找到赔率比 exp(-0.232) = 0.793，这意味着在保持其他条件不变的情况下，温度每上升一度，故障的赔率降低约 20%。需要注意的是，这并不意味着每次温度上升一度，故障的概率就会降低 20%。

因为赔率和概率之间的关系是非线性的，温度变化对失效概率的影响取决于温度变化发生的具体环境！

概率和几率通过 logit 和逻辑函数之间的逆关系相关联。逻辑函数有一个方便的性质，即对于任何输入x值，输出都在 0 到 1 的范围内——正好与概率的范围相同。此外，逻辑函数在绘制时创建一个 S 形曲线，如图图 6.6所示，该图显示了 O 形圈故障概率与发射温度的假设逻辑回归模型。在温度范围的中间，失败概率在y轴上变化最为强烈；在温度极端情况下，每增加或减少一度温度，预测的失败概率变化很小。

图 6.6：表示航天飞机发射数据的假设逻辑回归曲线

拟合的逻辑回归模型在 0 到 1 的范围内创建了一条曲线，表示对连续尺度上的概率估计，尽管目标结果（如图中圆圈所示）只取值为y = 0 或y = 1。为了获得二元预测，只需定义一个概率阈值，目标结果将据此进行预测。例如，如果预测的 O 形圈故障概率大于 0.50，则预测“故障”，否则预测“无故障”。使用 50%的阈值是常见的，但可以使用更高的或更低的阈值来调整模型对成本的敏感度。

检视图 6.6中的逻辑曲线会引出另一个问题：建模算法是如何确定最适合数据的曲线的呢？毕竟，鉴于这并非一条直线，标准线性回归中使用的 OLS 算法似乎不再适用。

事实上，广义线性模型使用一种称为最大似然估计（MLE）的不同技术，该技术找到指定分布中最有可能生成观察数据的参数值。

由于 OLS 估计是最大似然估计的特殊情况，只要满足 OLS 建模的假设，使用 OLS 或 MLE 对线性模型进行建模就没有区别。对于线性建模之外的应用，MLE 技术会产生不同的结果，并且必须使用 MLE 而不是 OLS。MLE 技术内置在 GLM 建模软件中，通常通过在数据上反复迭代以识别最优模型参数，而不是直接找到正确解。幸运的是，正如我们很快将看到的，在 R 中构建 GLM 几乎不比训练一个更简单的线性模型更具挑战性。

本介绍仅触及了线性回归和 GLM 可能性的表面。尽管理论和像挑战者数据集这样的简单例子有助于理解回归模型的工作原理，但构建一个有用的模型所涉及的不仅仅是我们所看到的。R 内置的回归函数包括拟合更复杂模型所需的功能，同时提供额外的诊断输出，以帮助模型解释和评估拟合度。让我们应用这些函数，通过尝试一个现实世界的学习任务来扩展我们对回归的了解。

示例 – 使用线性回归预测汽车保险索赔成本

对于一家汽车保险公司来说，为了盈利，它需要收取的会员费比支付给其受益人的车辆盗窃、损坏或事故中生命损失索赔要多。因此，保险公司投入时间和金钱来开发模型，以准确预测受保人群的索赔成本。这是被称为精算科学的领域，它使用复杂的统计技术来估计受保人群的风险。

由于事故，尤其是致命事故相对罕见——在美国，每行驶 1 亿英里车辆中略超过 1 起死亡事故——因此，个人保险费用难以准确预测。然而，当事故发生时，它们代价极高。此外，导致任何特定事故的具体条件基于难以衡量的因素，它们看起来似乎是随机的。一个有良好驾驶记录的优秀驾驶员可能会遭遇不幸，被醉酒驾驶员撞到，而另一个人可能会因分心驾驶手机，由于好运，从未造成事故。

由于预测个人费用的几乎不可能性，保险公司应用平均法则，计算具有相似风险特征的人群保险段的平均成本。如果每个风险段的费用估计是正确的，保险公司可以为风险较低的段定价较低的保险费，并可能从竞争的保险公司吸引新的低风险客户。在接下来的分析中，我们将模拟这种场景。

第 1 步 – 收集数据

本例的数据集是基于美国政府的人口统计和交通统计数据创建的模拟数据。它的目的是近似美国密歇根州汽车保险公司的现实世界条件，该州约有 1000 万居民和 700 万持牌驾驶员。

如果你想交互式地跟随，请从本书的 Packt Publishing GitHub 仓库下载autoinsurance.csv文件，并将其保存到你的 R 工作文件夹中。

保险数据集包括 20,000 个受益者参加的假设汽车车辆保险计划的示例。这比实际精算师使用的典型数据集要小得多，特别是对于非常罕见的结果，但规模已经缩小，以便即使在内存有限的计算机上也能进行分析。每个示例代表一个被保险个人的特征和计划在日历年度中收取的总保险索赔成本（费用）。在登记时可用的是以下特征：

• age: 驾驶者的年龄，从 16 岁到 89 岁

• geo_area: 车主主要居住地的地理区域，以及车辆最常使用的地方；邮编被归入城市、郊区和农村类别

• est_value: 根据年龄和折旧估计的车辆（们）的市场价值；上限为 125,000 美元——允许的最大保险价值

• vehicle_type: 乘客车辆的类型，可以是汽车、卡车、微型货车或运动型多功能车（SUV）

• miles_driven: 在日历年度内驾驶的距离（英里）

• college_grad_ind: 如果受益者拥有大学学历或更高，则该二进制指示器设置为1

• speeding_ticket_ind: 如果在过去五年内收到过超速罚单或违规，则该二进制指示器设置为1

• clean_driving_ind: 如果在过去五年内没有支付过责任保险索赔，则该二进制指示器设置为1

在这个示例场景中，有 20,000 名受益者参加了“安全驾驶折扣”计划，该计划要求使用设备或手机应用程序进行位置跟踪，以监控全年安全驾驶条件。这有助于验证miles_driven的准确性，并创建了以下两个额外的预测因子，旨在反映更危险的驾驶行为：

• hard_braking_ind: 如果车辆经常应用“硬刹车”（例如突然停车的情况），则该二进制指示器设置为1

• late_driving_ind: 如果车辆在午夜后经常驾驶，则该二进制指示器设置为1

考虑这些变量如何与保险费用相关是很重要的——有些方式比其他方式更明显。例如，我们显然预期经常驾驶的汽车比那些停在车库里的汽车更容易发生事故。另一方面，城市、农村或郊区的驾驶员哪个风险更高并不那么明显；农村驾驶员可能驾驶的距离更远，但城市驾驶涉及更多的交通，可能更容易发生车辆盗窃。回归模型将帮助我们解开这些关系，但需要我们自己指定特征之间的联系，而不是自动检测，这与许多其他机器学习方法不同。我们将在下一节中探讨一些潜在的关系。

考虑哪些可能的有用预测因子没有被包含在训练数据集中也许也很有趣。性别常被用于汽车保险定价（并且男女成本不同！）但密歇根州在 2020 年禁止了为此目的使用性别和信用评分。这些特征可能具有高度预测性，但可能导致对受保护群体的系统性偏见，如第一章介绍机器学习中所述。

第 2 步 – 探索和准备数据

如我们之前所做的那样，我们将使用read.csv()函数来加载数据进行分析。我们可以安全地使用stringsAsFactors = TRUE，因为将三个名义变量转换为因子是合适的：

> insurance <- read.csv("insurance.csv", stringsAsFactors = TRUE) 

str()函数确认数据格式正如我们所期望的：

> str(insurance) 

'data.frame':    20000 obs. of  11 variables:
 $ age                : int  19 30 39 64 33 27 62 39 67 38 ...
 $ geo_area           : Factor w/ 3 levels "rural","suburban", ...
 $ vehicle_type       : Factor w/ 4 levels "car","minivan", ...
 $ est_value          : int  28811 52603 113870 35228 ...
 $ miles_driven       : int  11700 12811 9784 17400 ...
 $ college_grad_ind   : int  0 1 1 0 0 1 1 0 1 1 ...
 $ speeding_ticket_ind: int  1 0 0 0 0 0 0 0 0 0 ...
 $ hard_braking_ind   : int  1 0 0 0 0 0 0 0 0 0 ...
 $ late_driving_ind   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ clean_driving_ind  : int  0 1 0 1 1 0 1 1 0 1 ...
 $ expenses           : num  0 6311 49684 0 0 ... 

我们的模型因变量是expenses，它衡量每个人在保险计划下一年内所声称的损失或损害。在构建线性回归模型之前，检查正态性通常很有帮助。尽管没有正态分布的因变量，线性回归也不会失败，但模型在这一点上通常拟合得更好。让我们看看总结统计：

> summary(insurance$expenses) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0    1709       0  232797 

最小值、第一四分位数、中位数和第三四分位数都是零，这意味着至少 75%的受益人在日历年度内没有费用。平均值大于中位数的事实让我们感觉到保险费用的分布是右偏的，但偏斜可能非常极端，因为平均费用是 1,709 美元，而最大值是 232,797 美元。我们可以使用直方图来直观地确认这一点：

> hist(insurance$expenses) 

图 6.7：年度保险索赔成本的分布

如预期的那样，该图显示了一个右偏分布，在零处有一个巨大的峰值，反映了只有一小部分（大约 8.6%）提出了保险索赔。在那些确实提出了车辆损失或损害索赔的人中，分布的尾部延伸到右侧，超过了 200,000 美元的昂贵伤害费用。尽管这种分布不适合线性回归，但提前知道这种弱点可能有助于我们稍后设计一个更好的拟合模型。现在，仅使用expenses的分布，我们可以这样说，平均受益人应该被收取每年 1,709 美元的保险费，这样保险公司才能收支平衡，或者每月每名订阅者约 150 美元的轻微利润。当然，这假设风险和成本是平均分担的。一个改进的保险模型会将更大的成本转嫁给风险较高的驾驶员，并为安全驾驶员提供经济上的节省。

在添加额外的预测变量之前，重要的是要注意回归模型要求每个特征都是数值型的，而我们的数据框中有两个因素类型的特征。例如，geo_area变量被分为urban、suburban和rural等级，而vehicle_type有car、truck、suv和minivan等类别。

让我们更仔细地看看它们的分布情况：

> table(insurance$geo_area) 

 rural suburban    urban 
    3622     8727     7651 

> table(insurance$vehicle_type) 

 car minivan     suv   truck 
   5801     726    9838    3635 

在这里，我们看到数据几乎均匀地分布在城市和郊区之间，但农村数据占比较小。此外，SUV 是最受欢迎的车型，其次是汽车和卡车，微型货车位于遥远的第四位。我们将很快看到 R 的线性回归函数如何处理这些因素变量。

探索特征之间的关系——相关矩阵

在将回归模型拟合到数据之前，确定自变量如何与因变量以及彼此相关可能很有用。相关矩阵提供了这些关系的快速概述。给定一组变量，它为每对关系提供相关性。

要为保险数据框中的四个数值型、非二进制变量创建相关矩阵，请使用cor()命令：

> cor(insurance[c("age", "est_value", "miles_driven", "expenses")]) 

 age   est_value miles_driven     expenses
age           1.000000000 -0.05990552   0.04812638 -0.009121269
est_value    -0.059905524  1.00000000  -0.01804807  0.088100468
miles_driven  0.048126376 -0.01804807   1.00000000  0.062146507
expenses     -0.009121269  0.08810047   0.06214651  1.000000000 

在每一行和每一列的交叉处，列出该行和列所指示变量的相关性。对角线始终是1.0000000，因为变量与其自身之间总是存在完美的相关性。对角线上方和下方的值是相同的，因为相关性是对称的。换句话说，cor(x, y)等于cor(y, x)。

矩阵中的相关性都不强，但关联与常识相符。例如，age和expenses似乎存在弱负相关性，这意味着随着年龄的增长，预期的保险费用会略有下降——这可能是由于更丰富的驾驶经验。还有est_value和expenses以及miles_driven和expenses之间的正相关性，这表明更昂贵的汽车和更广泛的驾驶会导致更高的费用。当我们构建最终的回归模型时，我们将尝试更清晰地揭示这些类型的关系。

可视化特征之间的关系——散点图矩阵

使用散点图可视化数值特征之间的关系也可能很有帮助。尽管我们可以为每种可能的关系创建一个散点图，但对于大量特征来说，这样做很快就会变得繁琐。

另一种选择是创建一个散点图矩阵（有时简称为SPLOM），它只是按网格排列的散点图的集合。它用于检测三个或更多变量之间的模式。散点图矩阵不是真正的多维可视化，因为一次只检查两个特征。尽管如此，它提供了数据可能相互关联的一般感觉。

我们可以使用 R 的图形功能来创建四个非二进制数值特征（age、est_value、miles_driven和expenses）的散点图矩阵。pairs()函数是 R 默认安装的一部分，提供了生成散点图矩阵的基本功能。要调用此函数，只需提供要绘制的数据框的子集。鉴于我们的insurance数据集相对较大，我们将设置绘图字符参数pch = "."为点，以便使可视化更容易阅读，然后限制列为我们感兴趣的四个变量：

> pairs(insurance[c("age", "est_value", "miles_driven",
                  "expenses")], pch = ".") 

这会产生以下散点图矩阵：

图 6.8：保险数据集中数值特征的散点图矩阵

在散点图矩阵中，每一行和每一列的交叉处持有由行和列对指示的变量的散点图。对角线以上和以下的图是转置的，因为x轴和y轴已经交换了。您在这些图中注意到任何模式吗？尽管它们大多看起来像随机的点云，但其中一些似乎显示出一些微妙趋势。est_value和miles_driven与expenses的关系似乎显示出轻微的上升趋势，这从相关矩阵中我们已经学到的内容得到了视觉上的证实。

通过向图中添加更多信息，可以使它变得更加有用。可以使用psych包中的pairs.panels()函数创建增强的散点图矩阵。如果您尚未安装此包，请在系统上输入install.packages("psych")进行安装，并使用library(psych)命令加载它。然后，我们可以使用pch参数设置绘图字符，就像我们之前做的那样来创建散点图矩阵：

> library(psych)
> pairs.panels(insurance[c("age", "est_value", "miles_driven",
                           "expenses")], pch = ".") 

这会产生一个更信息丰富的散点图矩阵，如下所示：

图 6.9：pairs.panels()函数为散点图矩阵添加了细节

在pairs.panels()输出中，对角线以上的散点图被相关矩阵所取代。对角线现在包含描述每个特征值分布的直方图。最后，对角线以下的散点图展示了额外的视觉信息。

每个散点图上的椭圆形物体（由于大量黑色点可能难以在打印中看到，但在电脑屏幕上更容易看到）是一个相关椭圆。它提供了一个简单的视觉指示，表示相关强度。在这个数据集中，没有强烈的关联，所以椭圆大多是平的；有更强的关联时，椭圆会向上或向下倾斜，以表示正相关或负相关。椭圆中心的小点是一个反映x轴和y轴变量平均值的点。

投影在散点图上的线（在计算机屏幕上显示为红色）称为局部加权回归曲线。它表示x轴和y轴变量之间的一般关系。最好通过例子来理解。尽管图表的小尺寸使得这种趋势难以看到，但age和miles_driven的曲线在达到中年之前略微上升，然后趋于平稳。这意味着驾驶往往随着年龄的增长而增加，直到它随时间大致保持恒定。

虽然在此处未观察到，但局部加权回归曲线有时可以非常显著，具有 V 形或 U 形曲线以及阶梯状模式。识别这些模式可以帮助后来开发更好的拟合回归模型。

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

要使用 R 拟合线性回归模型，可以使用lm()函数。这是stats包的一部分，应该默认包含并加载到您的 R 安装中。lm()的语法如下：

图 6.10：多重回归语法

以下命令拟合了一个线性回归模型，该模型将十个独立变量与总保险费用相关联。R 公式语法使用波浪线字符（~）来描述模型；因变量expenses写在波浪线的左侧，而独立变量则写在右侧，由加号（+）分隔。

没有必要指定回归模型的截距项，因为它默认包含：

> ins_model <- lm(expenses ~ age + geo_area + vehicle_type +
                    est_value + miles_driven +
                    college_grad_ind + speeding_ticket_ind +
                    hard_braking_ind + late_driving_ind +
                    clean_driving_ind,
                  data = insurance) 

因为点号字符（.）可以用来指定所有特征（不包括公式中已指定的那些），所以以下命令与之前的命令等价：

> ins_model <- lm(expenses ~ ., data = insurance) 

在构建模型后，只需输入模型对象的名称即可查看估计的贝塔系数。请注意，options(scipen = 999)命令关闭了科学记数法，以便更容易阅读输出：

> options(scipen = 999)
> ins_model 

Call:
lm(formula = expenses ~ ., data = insurance)
Coefficients:
        (Intercept)                  age     geo_areasuburban  
        -1154.91486             -1.88603            191.07895  
      geo_areaurban  vehicle_typeminivan      vehicle_typesuv  
          169.11426            115.27862            -19.69500  
  vehicle_typetruck            est_value         miles_driven  
           21.56836              0.03115              0.11899  
   college_grad_ind  speeding_ticket_ind     hard_braking_ind  
          -25.04030            155.82410             11.84522  
   late_driving_ind    clean_driving_ind  
          362.48550           -239.04740 

理解线性回归模型的回归系数相对简单。截距是当独立变量等于零时expenses的预测值。然而，在许多情况下，截距本身具有很小的解释价值，因为通常不可能所有特征都具有零值。

这就是这种情况，因为没有任何投保人可以有零岁或没有行驶英里，因此截距没有实际世界的解释。因此，在实践中，截距通常被忽略。

贝塔系数表示每个特征增加一个单位时，保险索赔成本的估计增加量，假设其他所有值保持不变。例如，对于每增加一岁，我们预计平均保险索赔成本将降低 1.89 美元，假设其他所有条件保持不变。同样，我们预计每增加一英里行驶将增加 0.12 美元的索赔，每增加一美元的保险价值将增加 0.03 美元，其他条件相同。

您可能会注意到，尽管我们在模型公式中只指定了 10 个特征，但除了截距之外，还报告了 13 个系数。这是因为 lm() 函数自动对模型中包含的每个因子类型变量应用虚拟编码。

如 第三章 所述，懒惰学习 – 使用最近邻进行分类，虚拟编码允许将名义特征通过为特征的每个类别（除了一个参考类别外）创建一个二元变量来处理为数值。每个虚拟变量在观测值属于指定类别时设置为 1，否则设置为 0。例如，geo_area 特征有三个类别：urban、suburban 和 rural。因此，使用了两个虚拟变量，分别命名为 geo_areaurban 和 geo_areasuburban。对于 geo_area = "rural" 的观测值，geo_areaurban 和 geo_areasuburban 都将设置为零。同样，对于四类 vehicle_type 特征，R 创建了三个虚拟变量，分别命名为 vehicle_typeminivan、vehicle_typesuv 和 vehicle_typetruck。这使 vehicle_type = "car" 在三个虚拟变量都为零时作为参考类别。

当在回归模型中使用虚拟编码特征时，回归系数的解释是相对于省略的类别。在我们的模型中，R 自动保留了 geo_arearural 和 vehicle_typecar 变量，使农村汽车所有者成为参考组。因此，与农村地区相比，城市居民每年的索赔成本高出 $169.11，而卡车每年比汽车使保险公司多花费$21.57。为了清楚起见，这些差异假设所有其他特征都保持相等，因此它们独立于农村驾驶员可能行驶更多里程或拥有更便宜车辆的事实。我们预计两个在其他方面完全相同的人，除了一个住在农村地区，一个住在城市地区，平均差异约为 $170。

默认情况下，R 使用因子变量的第一级作为参考。如果您希望使用其他级别，可以使用 relevel() 函数手动指定参考组。在 R 中使用 ?relevel 命令获取更多信息。

通常，线性回归模型的结果具有逻辑意义；然而，我们目前还没有关于模型如何拟合数据的良好感觉。我们将在下一节回答这个问题。

第 4 步 – 评估模型性能

通过输入 ins_model 获得的参数估计告诉我们独立变量与因变量之间的关系，但它们没有告诉我们模型如何拟合我们的数据。为了评估模型性能，我们可以使用存储的模型上的 summary() 命令：

> summary(ins_model) 

这产生了以下输出，其中已添加注释以供说明：

图 6.11：回归模型的总结输出可以分为三个主要部分，如图中所示

summary()输出的内容一开始可能看起来令人不知所措，但基本内容很容易掌握。如前述输出中的编号标签所示，评估我们模型性能或拟合度主要有三种方式：

1. 残差部分提供了预测误差的汇总统计，其中一些误差显然相当大。由于残差等于真实值减去预测值，最大误差231252表明模型至少对一个观测值预测费用低于 230,000 美元以上。另一方面，大多数误差是相对较小的负值，这意味着我们对大多数受保人的费用估计过高。这正是保险公司能够承担昂贵事故费用的原因。

2. 对于每个估计的回归系数，p 值（用Pr(>|t|)表示），提供了在给定估计值的情况下，真实系数为零的概率估计。小的 p 值表明真实系数极不可能为零，这意味着该特征极不可能与因变量没有关系。请注意，一些 p 值有星号（***），这对应于指定估计所达到的显著性水平的脚注。这个水平是一个阈值，在构建模型之前选择，将用于指示“真实”发现，而不是仅由偶然引起的那些；小于显著性水平的 p 值被认为是统计显著的。如果模型中这样的项很少，这可能是一个值得关注的问题，因为这表明所使用的特征对结果不是很有预测性。在这里，我们的模型有几个高度显著的自变量，并且它们似乎以预期的方式与结果相关。

3. 多重 R 平方值（也称为确定系数）提供了衡量我们的模型整体解释因变量值好坏的指标。它与相关系数类似，即值越接近 1.0，模型对数据的解释就越完美。由于 R 平方值为 0.01241，我们知道该模型解释了因变量变化的约 1.2%。由于具有更多特征的模型总是解释更多的变化，调整 R 平方值通过惩罚具有大量独立变量的模型来纠正 R 平方。这对于比较具有不同数量解释变量的模型性能是有用的。

根据前面的三个性能指标，我们的模型表现足够好。一些误差的大小有点令人担忧，但考虑到保险费用数据的性质，这并不令人惊讶。

此外，现实世界数据的回归模型具有低 R-squared 值并不罕见。尽管 0.01241 的值特别低，但它反映了我们没有汽车事故的近因预测因子；要真正预测事故，我们需要实时驾驶数据，或者至少需要一些衡量真实驾驶技能的指标。话虽如此，正如我们将在下一节中看到的，我们仍然可以通过以略微不同的方式指定模型来提高模型的表现。

第 5 步 – 提高模型性能

如前所述，回归建模与其他机器学习方法的显著区别在于，回归通常将特征选择和模型指定留给用户。因此，如果我们对某个特征与结果之间的关系有专业知识，我们可以利用这些信息来指导模型指定，并可能提高模型的表现。

模型指定 – 添加非线性关系

在线性回归中，假设自变量和因变量之间的关系是线性的，但这并不一定正确。例如，年龄对保险支出的影响可能不会在所有年龄值上保持恒定；对于最年轻和最年长的群体，治疗可能变得不成比例地昂贵——如果将费用与年龄绘制成曲线，则呈现 U 形曲线。

如果你还记得，典型的回归方程遵循类似以下的形式：

为了考虑非线性关系，我们可以在回归方程中添加一个更高阶的项，将模型视为多项式。实际上，我们将模拟如下关系：

这两个模型之间的区别在于，将估计一个额外的回归参数，目的是捕捉x²项的影响。这允许将年龄的影响作为年龄和年龄平方的函数来衡量。

要将非线性年龄添加到模型中，我们只需创建一个新的变量：

> insurance$age2 <- insurance$age² 

然后，当我们生成改进的模型时，我们将使用形式expenses ~ age + age2将age和age2都添加到lm()公式中。这将允许模型分离年龄对医疗费用的线性和非线性影响。

模型指定 – 添加交互效应

到目前为止，我们只考虑了每个特征对结果的单个贡献。如果某些特征对因变量有联合影响怎么办？例如，硬刹车和晚开车的不良习惯可能分别有有害的影响，但可以合理地假设它们的联合影响可能比单独每个的影响更糟。

当两个特征有联合效应时，这被称为交互作用。如果我们怀疑两个变量之间存在交互作用，我们可以通过将它们的交互作用添加到模型中来测试这个假设。交互作用使用 R 公式语法来指定。为了将硬刹车指标（hard_braking_ind）与晚驾指标（late_driving_ind）进行交互，我们将编写一个形式为expenses ~ hard_braking_ind*late_driving_ind的公式。

*运算符是一个简写，指示 R 对expenses ~ hard_braking_ind + late_driving_ind + hard_braking_ind:late_driving_ind进行建模。在展开形式中，冒号运算符（:）表示hard_braking_ind:late_driving_ind是两个变量之间的交互作用。请注意，展开形式自动还包括了单个hard_braking_ind和late_driving_ind变量以及它们的交互作用。

如果你在决定是否包含一个变量时遇到困难，一个常见的做法是先包含它，然后检查 p 值。如果变量在统计上不显著，你就有了一个合理的理由将其排除在模型之外。

将所有内容整合在一起——一个改进的回归模型

基于一些关于保险成本可能与报名者特征相关的主观知识，我们开发了一个我们认为更精确指定的回归公式。为了总结改进之处，我们：

• 为年龄添加了一个非线性项

• 指定硬刹车和晚驾之间的交互作用

我们将使用与之前相同的lm()函数来训练模型，但这次我们将添加交互项以及age2，它将自动包含：

> ins_model2 <- lm(expenses ~ . + hard_braking_ind:late_driving_ind,
                 data = insurance) 

接下来，我们总结一下结果：

> summary(ins_model2) 

输出如下：

Call:
lm(formula = expenses ~ . hard_barking_ind:late_driving_ind, 
    data = insurance)
Residuals: 
    Min       1Q   Median       3Q      Max  
  -6618    -1996    -1491    -1044   231358  
Coefficients:
                      Estimate   Std. Error  t value  Pr(>|z|)   
(Intercept)        -535.038171   457.146614   -1.170   0.2419 
age                 -33.142400   15.366892    -2.157   0.0310 *
geo_areasuburban    178.825158   143.305863    1.248   0.2121 
geo_areaurban       132.463265   158.726709    0.835   0.4040 
vehicle_typeminivan 178.825158   143.305863    1.248   0.2121 
vehicle_typesuv      -8.006108   118.116633   -0.068   0.9460 
vehicle_typetruck    26.426396   153.650455    0.172   0.8634 
est_value             0.031179   0.002496    12.489   <0.000000002 ***
miles_driven          0.118748   0.014327     8.289   <0.000000002 ***
college_grad_ind     17.248581   117.398583    0.147   0.8832 
speeding_ticket_ind 155.061583   140.143658    1.107   0.2658 
hard_braking_ind    -12.442358   109.794208   -0.113   0.9098 
late_driving_ind    183.329848   284.218859    0.645   0.5189 
clean driving_ind    -232.843170   111.106714    -2.096   0.0361 
age2                    0.343165     0.165340     2.076   0.0380 
hard_braking_ind:     469.079140   461.685886     1.016   0.3096 
 late_driving_ind
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6995 on 19984 degrees of freedom
Multiple R-squared: 0.01267,  Adjusted R-squared: 0.01193
F-statistic: 17.1 on 15 and 19984 DF,
p-value: <0.00000000000000022 

尽管与之前的模型相比，R 平方和调整 R 平方值变化不大，但新特征提供了一些有趣的见解。特别是，age的估计值相对较大且为负（支出较低），而age2的估计值相对较小且为正（支出较高）。然而，由于年龄平方的增长速度快于年龄，对于年龄非常高的群体，支出将会开始上升。整体效果是一个 U 形的支出曲线，预测最年轻和最年长的报名者将有更高的支出。hard_braking_ind和late_driving_ind的交互作用也很有趣，因为它是一个相对较大的正值。尽管交互作用在统计上不显著，但效应的方向暗示了如果你是那种已经危险驾驶的司机，那么在晚上晚些时候驾驶尤其危险。

严格来说，回归建模对数据做出了一些强烈的假设。这些假设对于数值预测并不那么重要，因为模型的价值并不在于它是否真正捕捉到了潜在的过程——我们只关心其预测的准确性。然而，如果你想要从回归模型系数中得出明确的推断，就必须运行诊断测试以确保回归假设没有被违反。关于这个主题的优秀介绍，请参阅 《多元回归：入门》，Allison, P. D.，Pine Forge Press，1998 年。

使用回归模型进行预测

在检查了估计的回归系数和拟合统计量之后，我们还可以使用该模型来预测保险计划未来参保人的支出。为了说明预测过程，让我们首先使用 predict() 函数将模型应用于原始训练数据，如下所示：

> insurance$pred <- predict(ins_model2, insurance) 

这将预测结果保存为名为 pred 的新向量，在保险数据框中。然后我们可以计算预测的保险成本与实际成本之间的相关性：

> cor(insurance$pred, insurance$expenses) 

[1] 0.1125714 

0.11 的相关性表明预测值和实际值之间存在着相对较弱的线性关系，这在某种程度上令人失望，但鉴于交通事故看似随机的性质，这并不太令人惊讶。将这一发现作为散点图来考察也是有用的。以下 R 命令绘制了这种关系，并添加了一条截距为零、斜率为一的识别线。col、lwd 和 lty 参数分别影响线的颜色、宽度和类型：

> plot(insurance$pred, insurance$expenses)
> abline(a = 0, b = 1, col = "red", lwd = 3, lty = 2) 

图 6.12：在这个散点图中，落在或接近对角虚线（y = x）上的点表示预测值与实际值非常接近

位于线上方的非对角点表示实际支出大于预期的案例，而位于线下方的案例表示支出小于预期的案例。我们可以看到，少数支出远大于预期的个人被大量支出略小于预期的个人所平衡。

现在，假设你想预测保险计划中潜在新参保人的支出。为此，你必须向 predict() 函数提供一个包含潜在驾驶员数据的数据框。对于许多驾驶员，你可能考虑创建一个 CSV 电子表格文件以在 R 中加载，或者对于较少的驾驶员，你可以在 predict() 函数内部直接创建一个数据框。例如，为了估计一个 30 岁、居住在农村、驾驶价值为 25,000 美元的卡车、每年行驶约 14,000 英里且驾驶记录良好的驾驶员的保险费用：

> predict(ins_model2,
          data.frame(age = 30, age2 = 30², geo_area = "rural", 
                     vehicle_type = "truck", est_value = 25000,
                     miles_driven = 14000, college_grad_ind = 0,
                     speeding_ticket_ind = 0, hard_braking_ind = 0,
                     late_driving_ind = 0, clean_driving_ind = 1)) 

 1 
1015.059 

使用这个值，保险公司需要每年收取约 1015 美元才能为这个人口群体实现盈亏平衡。要比较一个在其他方面都与上述情况相似，但最近有过事故记录的人的费率，可以使用predict()函数以类似的方式：

> predict(ins_model2,
          data.frame(age = 30, age2 = 30², geo_area = "rural", 
                     vehicle_type = "truck", est_value = 25000,
                     miles_driven = 14000, college_grad_ind = 0,
                     speeding_ticket_ind = 0, hard_braking_ind = 0,
                     late_driving_ind = 0, clean_driving_ind = 0)) 

 1 
1247.903 

注意，这两个值之间的差异，1015.059 – 1247.903 = -232.844，与估计的回归模型系数clean_driving_ind相同。平均而言，拥有良好驾驶记录的驾驶员预计每年在计划上的支出将少约 232.84 美元，其他条件相同。

这说明了更普遍的事实，即预测的支出是每个回归系数与其预测数据框中相应值的乘积之和。例如，使用模型对行驶里程的回归系数 0.118748，我们可以预测增加 10,000 英里将导致支出增加10,000 * 0.118748 = 1187.48，如下所示：

> predict(ins_model2,
          data.frame(age = 30, age2 = 30², geo_area = "rural", 
                     vehicle_type = "truck", est_value = 25000,
                     miles_driven = 14000, college_grad_ind = 0,
                     speeding_ticket_ind = 0, hard_braking_ind = 0,
                     late_driving_ind = 0, clean_driving_ind = 0)) 

 1 
1247.903 

> 2435.384 - 1247.903 

[1] 1187.481 

通过对多个额外的客户风险细分进行类似的步骤，保险公司能够开发出一个定价结构，该结构根据驾驶员的估计风险水平公平地设定成本，同时在整个细分市场保持一致的利润。

导出模型的回归系数允许你构建自己的预测函数。这样做的一个潜在用例是在客户数据库中实现回归模型，以进行实时预测。

进一步——使用逻辑回归预测保险保单持有人流失率

在保险公司内部，机器学习的潜在应用不仅限于对索赔成本的精算估计。营销和客户保留团队很可能对预测流失率，即选择不续保保险计划后离开公司的客户感兴趣。在许多业务中，防止流失率被高度重视，因为流失的客户不仅会减少一个企业的收入流，而且通常还会增加直接竞争对手的收入流。此外，营销团队知道，获取新客户的成本通常远高于保留现有客户的成本。

因此，提前知道哪些客户最有可能流失可以帮助将保留资源用于干预，防止流失发生。

历史上，营销团队使用一个名为最近购买频率、货币价值（RFM）的简单模型来识别高价值客户以及最有可能流失的客户。RFM 分析考虑了每位客户的三个特征：

• 他们最近一次购买是什么时候？一段时间没有购买的客户可能价值较低，也更可能永远不会回来。

• 他们购买频率如何？他们是一年复一年地回来，还是购买行为中存在不规则的间隔？表现出忠诚度的客户可能更有价值，也更可能回头。

• 他们在购买时花费多少钱？他们是否比平均客户花费更多或升级到高级产品？这些客户在财务上更有价值，同时也表现出对品牌的喜爱。

收集历史客户购买数据是为了开发这三个因素中每个因素的度量标准。然后，将这些度量标准转换为每个三个领域的标准尺度（例如，从零到十的尺度），并将它们相加以创建每个客户的最终 RFM 分数。一个最近且频繁购买且平均花费的客户可能的总分为 10 + 10 + 5 = 25，而一个很久以前只购买过一次的客户可能得分为 2 + 1 + 4 = 7，这在 RFM 尺度上要低得多。

这种分析是一种粗略但有用的工具，用于理解一组客户并帮助识别可能对预测流失有用的数据类型。然而，RFM 分析并不特别科学，也不提供正式的流失概率估计或增加其可能性的因素。相比之下，预测二元流失结果的逻辑回归模型为每个客户提供流失的估计概率，以及每个预测变量的影响。

就像保险索赔成本示例一样，我们将使用为本书创建的模拟数据集构建一个流失模型，该数据集旨在模拟汽车保险公司的客户行为。

如果您想交互式地跟随学习，请从本书的 Packt Publishing GitHub 仓库下载 insurance_churn.csv 文件，并将其保存到您的 R 工作文件夹中。

流失数据集包括 5,000 个当前和以前的受益人，他们参加了假设的汽车保险计划。每个示例包括衡量计划年度客户行为的特征，以及一个二进制指标（churn），表示他们是否在年底未续约而退出计划。可用的特征包括：

• member_id：一个随机分配的客户识别号。

• loyalty_years：连续参加保险计划的年数。

• vehicles_covered：保险计划覆盖的车辆数量。

• premium_plan_ind：一个二进制指标，表示成员支付了包含额外福利的昂贵计划的高级版本。

• mobile_app_user：一个二进制指标，表示成员使用手机伴侣应用程序。

• home_auto_bundle：一个二进制指标，表示成员还持有由同一家公司提供的房主保险计划。

• auto_pay_ind：一个二进制指标，表示成员已开启自动支付。

• recent_rate_increase：一个二进制指标，表示成员的价格最近有所提高。

注意，许多这些因素与 RFM 的三个组成部分相关，因为它们是忠诚度和货币价值的衡量标准。因此，即使后来构建了一个更复杂的模型，仍然可以在处理过程的早期步骤中进行 RFM 分析，这仍然是有帮助的。

要将此数据集读入 R，请输入：

> churn_data <- read.csv("insurance_churn.csv") 

使用 table() 和 prop.table() 函数，我们可以看到整体流失率仅为 15%以上：

> prop.table(table(churn_data$churn)) 

 0      1 
0.8492 0.1508 

在更正式的分析中，在进行进一步分析之前进行更多数据探索是明智的。在这里，我们将跳过创建逻辑回归模型来预测这些流失客户。

glm() 函数是 R 内置的 stats 包的一部分，用于拟合 GLM 模型，如逻辑回归以及本章前面描述的其他变体，如泊松回归。逻辑回归的语法如下所示：

图 6.13：逻辑回归语法

注意 glm() 函数语法与之前用于标准线性回归的 lm() 函数之间的许多相似之处。除了指定家族和链接函数外，拟合模型并不更困难。主要差异主要在于如何解释生成的模型。

注意，R 的 glm() 函数默认使用高斯分布和恒等链接，因此在需要其他 GLM 形式时，很容易意外地执行标准线性回归！因此，在 R 中构建 GLM 时，始终指定家族和链接是一个明智的习惯。

要拟合逻辑回归流失模型，我们指定 binomial 家族和 logit 链接函数。在这里，我们将流失建模为数据集中除 member_id 之外所有其他特征的函数，member_id 对每个成员来说是唯一的，因此对预测无用：

> churn_model <- glm(churn ~ . -member_id, data = churn_data,
                     family = binomial(link = "logit")) 

使用 summary() 对生成的 churn_model 对象进行操作，可以显示估计的回归参数：

> summary(churn_model) 

Call:
glm(formula = churn ~ . - member_id,
      family = binomial(link = "logit"), data = ins_churn)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1244  -0.6152  -0.5033  -0.3950   2.4995  
Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -0.488893   0.141666  -3.451 0.000558 ***
loyalty_years        -0.072284   0.007193 -10.050  < 2e-16 ***
vehicles_covered     -0.212980   0.055237  -3.856 0.000115 ***
premium_plan_ind     -0.370574   0.148937  -2.488 0.012842 *  
mobile_app_user      -0.292273   0.080651  -3.624 0.000290 ***
home_auto_bundle     -0.267032   0.093932  -2.843 0.004472 ** 
auto_pay_ind         -0.075698   0.106130  -0.713 0.475687    
recent_rate_increase  0.648100   0.102596   6.317 2.67e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 4240.9  on 4999  degrees of freedom
Residual deviance: 4059.2  on 4992  degrees of freedom
AIC: 4075.2
Number of Fisher Scoring iterations: 5 

从高层次来看，逻辑回归的输出与线性回归的输出相当相似。p 值（标记为 Pr(>|z|)) 和显著性代码（由 * 字符表示）表明变量是否具有统计学意义。除了 auto_pay_ind 之外的所有特征在 0.05 水平或更好上都是显著的。通过查看 Estimate 值之前的关系（正或负）也可以简单地理解预测变量与目标结果之间的关系。几乎所有估计值都是负的，这意味着这些特征会减少流失，除了 recent_rate_increase 是正的，因此会增加流失。这些联系是有意义的；预计保险计划的价格上涨会增加流失，而忠诚度多年或支付高级计划特征的成员不太可能离开。

解释特定特征对流失的影响比线性回归更困难，因为估计值以对数优势形式显示。假设我们想知道在保险计划价格最近增加后，流失的可能性增加了多少。由于recent_rate_increase的估计值为 0.6481，这意味着当增加指标为1时，流失的对数优势增加 0.6481，而当它为0时。对此进行指数化以消除对数并找到优势比，我们发现exp(0.6481) = 1.911905，这意味着在增加费率后，流失的可能性几乎是两倍（或增加了 91.2%）。

在相反的方向上，使用移动应用（mobile_app_user）的成员与未使用应用的成员相比，估计的对数优势差异为-0.292273。将优势比计算为exp(-0.292273) = 0.7465647表明，应用用户的流失率大约是未使用应用用户的 75%，或者应用用户减少了大约 25%。同样，我们可以发现，对于忠诚度每增加一年，流失率会减少大约 7%，如exp(-0.072284) = 0.9302667。对于模型中的所有其他预测因子，包括截距（代表所有预测因子为零时的流失优势），都可以进行类似的计算。

要使用此模型来防止流失，我们可以在当前计划成员的数据库上进行预测。让我们从这个章节可用的test数据集开始，加载包含 1000 个订阅者的数据集：

> churn_test <- read.csv("insurance_churn_test.csv") 

然后，我们将使用带有predict()函数的逻辑回归模型对象来向这个数据框添加一个新列，其中包含每个成员的预测值：

> churn_test$churn_prob <- predict(churn_model, churn_test,
                                   type = "response") 

注意，type = "response"参数被设置为使预测以概率形式呈现，而不是默认的type = "link"设置，后者以对数优势值的形式产生预测。

总结这些预测概率，我们看到平均流失概率大约为 15%，但一些用户预测的流失率非常低，而其他用户的流失概率高达 41%：

> summary(churn_test$churn_prob) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02922 0.09349 0.13489 0.14767 0.18452 0.41604 

假设客户保留团队有资源在有限数量的案例中干预。通过排序成员以识别那些预测流失可能性最高的成员，我们可以为团队提供最有可能产生最大影响的指导方向。

首先，使用order()函数获取一个按流失概率降序排列的行号向量：

> churn_order <- order(churn_test$churn_prob, decreasing = TRUE) 

接下来，在根据churn_order向量对churn_test数据框进行排序并取两个感兴趣的列之后，使用head()函数取前n行；在这种情况下，我们将n = 5设置为限制流失可能性最高的五个成员：

> head(churn_test[churn_order,
        c("member_id", "churn_prob")], n = 5) 

 member_id churn_prob
406  29603520  0.4160438
742  12588881  0.4160438
390  23228258  0.3985958
541  86406649  0.3985958
614  49806111  0.3985958 

在将结果保存到具有n设置为更高数值的电子表格后，就可以为保留团队提供一份最有可能流失的保险计划成员名单。将保留努力集中在这些成员上，可能比随机针对成员更有利于营销预算的利用，因为大多数成员的流失概率非常低。通过这种方式，机器学习可以在最小投资的情况下提供实质性的回报，其影响可以通过比较干预前后的流失率来轻松量化。

可以通过关于将响应保留努力的客户比例的简单假设来获得防止客户流失所保留的收益的估计。例如，如果我们假设流失模型针对的N名成员将被保留，这将导致N倍的*$X*保留收入，其中*$X*是平均客户消费。将这个数字带给利益相关者有助于为实施机器学习项目提供正当理由。

这个例子只是冰山一角，因为通过额外的努力，客户流失建模可以变得更加复杂。例如，我们不仅可以将目标对准流失概率最高的客户，还可以考虑如果客户流失将失去的收入；即使比低价值客户流失概率低，优先考虑高价值客户也可能是值得的。此外，由于一些客户无论是否干预都肯定会流失，而另一些客户可能更愿意留下，因此除了建模流失概率外，还可以建模保留的可能性。在任何情况下，即使以简单形式，客户流失建模对大多数企业来说都是低垂的果实，并且是实施机器学习的绝佳起点。

理解回归树和模型树

如果你还记得第五章，分而治之 – 使用决策树和规则进行分类，决策树构建了一个模型，就像流程图一样，其中决策节点、叶节点和分支定义了一系列决策，这些决策用于分类示例。这样的树也可以通过仅对树生长算法进行微小调整来用于数值预测。在本节中，我们将考虑用于数值预测的树与用于分类的树的不同之处。

用于数值预测的树分为两类。第一类，称为回归树，在 20 世纪 80 年代作为分类与回归树（CART）算法的一部分被引入。尽管名称如此，回归树并不使用本章前面描述的线性回归方法；相反，它们根据达到叶节点的示例的平均值进行预测。

CART 算法在*《分类与回归树》，Breiman, L.，Friedman, J. H.，Stone, C. J.，Olshen, R. A.，Chapman and Hall，1984*中有详细描述。

用于数值预测的第二种树被称为模型树。比回归树晚几年引入，它们不太为人所知，但可能更强大。模型树的生长方式与回归树非常相似，但在每个叶节点，都会从到达该节点的示例中构建一个多元线性回归模型。根据叶节点的数量，模型树可能会构建数十甚至数百个这样的模型。

这使得模型树比等效的回归树更难以理解，但它们可能产生更精确的模型。

最早的模型树算法M5在使用连续类别的学习，Quinlan, J. R.，第五次澳大利亚联合人工智能会议论文集，1992 年，第 343-348 页中进行了描述。

将回归添加到树中

能够执行数值预测的树提供了一个引人注目但往往被忽视的回归建模替代方案。以下表格列出了相对于更常见的回归方法，回归树和模型树的优缺点：

优点	缺点

|

• 结合了决策树的优势以及建模数值数据的能力

• 不需要用户事先指定模型

• 使用自动特征选择，这使得方法可以与大量特征一起使用

• 可能比线性回归更适合某些类型的数据

• 解释模型不需要了解统计学知识

|

• 不如线性回归知名

• 需要大量的训练数据

• 难以确定单个特征对结果的整体净效应

• 大树可能比回归模型更难以解释

|

尽管传统的回归方法通常是数值预测任务的第一个选择，但在某些情况下，数值决策树提供了独特的优势。例如，决策树可能更适合具有许多特征或特征与结果之间存在许多复杂、非线性关系的任务；这些情况对回归构成了挑战。回归建模也假设数据具有某些属性，但在现实世界的数据中这些属性通常被违反；而对于树来说并非如此。

用于数值预测的树与用于分类的树构建方式非常相似。从根节点开始，数据根据将导致分割后结果同质性最大增加的特征，使用分而治之的策略进行分区。

在分类树中，你可能记得同质性是通过熵来衡量的。对于数值数据，这将是未定义的。相反，对于数值决策树，同质性是通过诸如方差、标准差或平均值的绝对偏差等统计量来衡量的。

一种常见的分割标准被称为标准差减少（SDR）。它由以下公式定义：

在这个公式中，sd(T)函数指的是集合T中值的标准差，而T[1]，T[2]，...，T[n]是特征分割后得到的值集。|T|项指的是集合T中的观测值数量。本质上，该公式通过比较分割前后的加权标准差来衡量标准差的减少。

例如，考虑以下情况，其中树正在决定是否在二元特征 A 上进行分割或在二元特征 B 上进行分割：

图 6.14：算法考虑在特征 A 和 B 上进行分割，这会创建不同的 T[1]和 T[2]组

使用由提议的分割产生的组，我们可以计算 A 和 B 的 SDR 如下。这里使用的length()函数返回向量的元素数量。请注意，整体组 T 被命名为tee，以避免覆盖 R 的内置T()和t()函数。

> tee <- c(1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7)
> at1 <- c(1, 1, 1, 2, 2, 3, 4, 5, 5)
> at2 <- c(6, 6, 7, 7, 7, 7)
> bt1 <- c(1, 1, 1, 2, 2, 3, 4)
> bt2 <- c(5, 5, 6, 6, 7, 7, 7, 7)
> sdr_a <- sd(tee) - (length(at1) / length(tee) * sd(at1) +
             length(at2) / length(tee) * sd(at2))
> sdr_b <- sd(tee) - (length(bt1) / length(tee) * sd(bt1) +
             length(bt2) / length(tee) * sd(bt2)) 

让我们比较 A 的 SDR 与 B 的 SDR：

> sdr_a 

[1] 1.202815 

> sdr_b 

[1] 1.392751 

在特征 A 上的分割 SDR 约为 1.2，而在特征 B 上的分割 SDR 为 1.4。由于 B 的分割降低了更多的标准差，决策树将首先使用 B。这比 A 产生了稍微更均匀的集合。

假设树在这里停止生长，使用这个唯一的一次分割。回归树的工作就完成了。它可以根据示例在特征 B 上的值，将示例放入组T[1]或T[2]，从而对新示例进行预测。如果示例最终落在T[1]中，模型将预测mean(bt1) = 2，否则它将预测mean(bt2) = 6.25。

相比之下，模型树会更进一步。使用落在T[1]组的七个训练示例和T[2]组的八个示例，模型树可以构建一个关于特征 A 的线性回归模型。请注意，特征 B 在构建回归模型时没有帮助，因为所有叶节点上的 B 值都相同——它们根据 B 的值被放入T[1]或T[2]。然后模型树可以使用这两个线性模型中的任何一个对新示例进行预测。

为了进一步说明这两种方法之间的差异，让我们通过一个现实世界的例子来分析。

示例 - 使用回归树和模型树估计葡萄酒的质量

葡萄酒酿造是一项具有挑战性和竞争性的业务，具有巨大的盈利潜力。然而，有许多因素会影响酒庄的盈利能力。作为一个农产品，从天气到生长环境等各种各样的变量都会影响品种的质量。装瓶和制造也可能影响口味，无论是好是坏。甚至产品的营销方式，从瓶子设计到价格点，都可能影响顾客对味道的认知。

因此，酿酒业在数据收集和可能协助酿酒决策科学的机器学习方法上投入了大量资金。例如，机器学习已被用于发现不同地区葡萄酒化学成分的关键差异，以及识别导致葡萄酒味道更甜的化学因素。

最近，机器学习已被用于协助评估葡萄酒的质量——这是一个众所周知困难的任务。一篇由著名葡萄酒评论家撰写的评论通常决定了产品最终会放在货架的顶部还是底部，尽管即使是专家评委在盲品测试中对葡萄酒的评分也不一致。

在本案例研究中，我们将使用回归树和模型树创建一个能够模仿专家葡萄酒评分的系统。由于树可以产生易于理解的模型，这可能会让酿酒师识别出有助于获得更高评分的关键因素。也许更重要的是，该系统不受品尝过程中的人类因素影响，如评分者的情绪或味蕾疲劳。因此，计算机辅助的葡萄酒测试可能因此产生更好的产品，以及更客观、一致和公平的评分。

第一步 – 收集数据

为了开发葡萄酒评分模型，我们将使用 P. Cortez、A. Cerdeira、F. Almeida、T. Matos 和 J. Reis 捐赠给 UCI 机器学习仓库（archive.ics.uci.edu/ml）的数据。他们的数据集包括来自葡萄牙的红色和白色 Vinho Verde 葡萄酒的示例——葡萄牙是世界上领先的葡萄酒生产国之一。由于影响高度评分葡萄酒的因素可能在红葡萄酒和白葡萄酒品种之间有所不同，因此，在本分析中，我们将仅检查更受欢迎的白葡萄酒。

为了跟随这个示例，请从本书的 Packt Publishing GitHub 仓库下载whitewines.csv文件，并将其保存到您的 R 工作目录中。如果您想自己探索这些数据，redwines.csv文件也是可用的。

白葡萄酒数据包括 4,898 个酒样 11 种化学特性的信息。对于每一款酒，实验室分析测量了如酸度、糖含量、氯化物、硫、酒精、pH 值和密度等特征。然后，由不少于三位评委在从 0（非常差）到 10（优秀）的质量等级上进行盲品评分。如果评委对评分有分歧，则使用中位数值。

Cortez 的研究评估了三种机器学习方法建模葡萄酒数据的能力：多元回归、人工神经网络和支持向量机。我们在本章前面介绍了多元回归，我们将在第七章，黑盒方法 – 神经网络和支持向量机中学习神经网络和支持向量机。该研究发现，支持向量机模型比线性回归模型提供了显著更好的结果。然而，与回归不同，支持向量机模型难以解释。使用回归树和模型树，我们可能能够在保持模型相对容易理解的同时提高回归结果。

要了解更多关于这里描述的葡萄酒研究的信息，请参阅Cortez, P.，Cerdeira, A.，Almeida, F.，Matos, T.，和 Reis, J.，通过数据挖掘物理化学特性建模葡萄酒偏好，决策支持系统，2009 年，第 47 卷，第 547-553 页。

第 2 步 – 探索和准备数据

如往常一样，我们将使用read.csv()函数将数据加载到 R 中。由于所有特征都是数值型的，我们可以安全地忽略stringsAsFactors参数：

> wine <- read.csv("whitewines.csv") 

wine数据包括 11 个特征和品质结果，如下所示：

> str(wine) 

'data.frame':	4898 obs. of  12 variables:
 $ fixed.acidity       : num  6.7 5.7 5.9 5.3 6.4 7 7.9 ...
 $ volatile.acidity    : num  0.62 0.22 0.19 0.47 0.29 ...
 $ citric.acid         : num  0.24 0.2 0.26 0.1 0.21 0.41 ...
 $ residual.sugar      : num  1.1 16 7.4 1.3 9.65 0.9 ...
 $ chlorides           : num  0.039 0.044 0.034 0.036 0.041 ...
 $ free.sulfur.dioxide : num  6 41 33 11 36 22 33 17 34 40 ...
 $ total.sulfur.dioxide: num  62 113 123 74 119 95 152 ...
 $ density             : num  0.993 0.999 0.995 0.991 0.993 ...
 $ pH                  : num  3.41 3.22 3.49 3.48 2.99 3.25 ...
 $ sulphates           : num  0.32 0.46 0.42 0.54 0.34 0.43 ...
 $ alcohol             : num  10.4 8.9 10.1 11.2 10.9 ...
 $ quality             : int  5 6 6 4 6 6 6 6 6 7 ... 

与其他类型的机器学习模型相比，树模型的一个优点是它们可以处理许多类型的数据而无需预处理。这意味着我们不需要对特征进行归一化或标准化。

然而，为了了解模型性能的评估，我们需要对结果变量的分布进行一些检查。例如，假设葡萄酒之间的质量变化非常小，或者葡萄酒落入双峰分布：要么非常好，要么非常差。这可能会影响我们设计模型的方式。为了检查这种极端情况，我们可以使用直方图检查葡萄酒质量的分布：

> hist(wine$quality) 

这产生了以下图示：

图 6.15：白葡萄酒质量评分的分布

葡萄酒质量值似乎遵循一个大致正常的、钟形的分布，中心值为六。从直观上看，这是有道理的，因为大多数葡萄酒的平均质量；很少有特别差或好的。尽管这里没有显示结果，但检查summary(wine)输出以查找异常值或其他潜在的数据问题也是有用的。尽管树模型对杂乱的数据相当稳健，但始终检查严重问题总是谨慎的。目前，我们假设数据是可靠的。

因此，我们的最后一步是将数据集划分为训练集和测试集。由于wine数据集已经被随机排序，我们可以将其划分为两个连续行集，如下所示：

> wine_train <- wine[1:3750, ]
> wine_test <- wine[3751:4898, ] 

为了与 Cortez 使用的条件相匹配，我们分别使用了 75% 和 25% 的数据集用于训练和测试。我们将评估基于树的模型在测试数据上的性能，以查看我们是否能获得与先前研究相似的结果。

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

我们将首先训练一个回归树模型。尽管几乎任何决策树的实现都可以用于回归树建模，但 rpart（递归分区）包提供了最忠实的回归树实现，正如 CART 团队所描述的那样。作为经典的 R 实现 CART，rpart 包也具有良好的文档支持，并提供了用于可视化和评估 rpart 模型的函数。

使用 install.packages("rpart") 命令安装 rpart 包。然后可以使用 library(rpart) 语句将其加载到您的 R 会话中。以下语法将使用默认设置训练一个树，这些设置在大多数情况下都很好，但并不总是如此。如果您需要更精细的设置，请参阅 ?rpart.control 命令的文档以了解控制参数。

图 6.16：回归树语法

使用 R 公式接口，我们可以将 quality 指定为结果变量，并使用点符号允许 wine_train 数据框中的所有其他列用作预测变量。生成的回归树模型对象命名为 m.rpart，以区分我们稍后将要训练的模型树：

> m.rpart <- rpart(quality ~ ., data = wine_train) 

要获取关于树的基本信息，只需输入模型对象的名称：

> m.rpart 

n= 3750 
node), split, n, deviance, yval
      * denotes terminal node
 1) root 3750 2945.53200 5.870933  
   2) alcohol< 10.85 2372 1418.86100 5.604975  
     4) volatile.acidity>=0.2275 1611  821.30730 5.432030  
       8) volatile.acidity>=0.3025 688  278.97670 5.255814 *
       9) volatile.acidity< 0.3025 923  505.04230 5.563380 *
     5) volatile.acidity< 0.2275 761  447.36400 5.971091 *
   3) alcohol>=10.85 1378 1070.08200 6.328737  
     6) free.sulfur.dioxide< 10.5 84   95.55952 5.369048 *
     7) free.sulfur.dioxide>=10.5 1294  892.13600 6.391036  
      14) alcohol< 11.76667 629  430.11130 6.173291  
        28) volatile.acidity>=0.465 11   10.72727 4.545455 *
        29) volatile.acidity< 0.465 618  389.71680 6.202265 *
      15) alcohol>=11.76667 665  403.99400 6.596992 * 

对于树中的每个节点，都会列出到达决策点的示例数量。例如，所有 3,750 个示例都从根节点开始，其中 2,372 个具有 alcohol < 10.85，1,378 个具有 alcohol >= 10.85。因为 alcohol 在树中首先被使用，所以它是葡萄酒质量的最重要预测变量。

由 * 标记的节点是终端或叶节点，这意味着它们会导致一个预测（在此处列为 yval）。例如，节点 5 的 yval 值为 5.971091。当使用此树进行预测时，任何酒精含量小于 10.85 且挥发性酸度小于 0.2275 的葡萄酒样品都会被预测为具有 5.97 的质量值。

使用 summary(m.rpart) 命令可以获得关于树拟合的更详细摘要，包括每个节点的均方误差以及特征重要性的总体度量。

可视化决策树

虽然可以使用前面的输出理解树，但使用可视化通常更容易理解。Stephen Milborrow 的 rpart.plot 包提供了一个易于使用的函数，可以生成高质量的决策树。

关于rpart.plot的更多信息，包括该函数可以生成的决策树图类型的额外示例，请参阅作者的网站www.milbo.org/rpart-plot/。

使用install.packages("rpart.plot")命令安装包后，rpart.plot()函数可以从任何rpart模型对象生成树形图。以下命令绘制了我们之前构建的回归树：

> library(rpart.plot)
> rpart.plot(m.rpart, digits = 3) 

生成的树形图如下：

图 6.17：葡萄酒质量回归树模型的可视化

除了控制图中包含的数字位数的digits参数外，许多其他可视化方面都可以进行调整。以下命令显示了其中一些有用的选项：

> rpart.plot(m.rpart, digits = 4, fallen.leaves = TRUE,
               type = 3, extra = 101) 

fallen.leaves参数强制叶节点对齐在图的底部，而type和extra参数影响决策和节点标签的方式。数字3和101指的是特定的样式格式，可以在命令的文档中找到，或者通过尝试不同的数字进行实验。

这些更改的结果是一个看起来非常不同的树形图：

图 6.18：更改绘图函数参数允许自定义树形可视化

这种类型的可视化有助于传播回归树结果，因为即使没有数学背景也能轻松理解。在两种情况下，叶节点中显示的数字都是到达该节点的示例的预测值。因此，向葡萄酒生产商展示此图可能有助于识别预测高分葡萄酒的关键因素。

第 4 步 – 评估模型性能

要使用回归树模型对测试数据进行预测，我们使用predict()函数。默认情况下，这返回了结果变量的估计数值，我们将将其保存在名为p.rpart的向量中：

> p.rpart <- predict(m.rpart, wine_test) 

快速查看预测的摘要统计数据显示了一个潜在问题：预测值落在比真实值更窄的范围内：

> summary(p.rpart) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  4.545   5.563   5.971   5.893   6.202   6.597 

> summary(wine_test$quality) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  3.000   5.000   6.000   5.901   6.000   9.000 

这个发现表明，模型没有正确识别极端情况，特别是最好的和最差的葡萄酒。另一方面，在第一和第三四分位数之间，我们可能做得很好。

预测值和实际质量值之间的相关性提供了一个简单的方法来衡量模型的表现。回想一下，cor()函数可以用来测量两个等长向量之间的关系。我们将使用它来比较预测值与真实值之间的对应程度：

> cor(p.rpart, wine_test$quality) 

[1] 0.5369525 

0.54的相关性当然是可以接受的。然而，相关性仅衡量预测与真实值之间的相关性强度；它不是衡量预测偏离真实值程度的指标。

使用平均绝对误差衡量性能

另一种思考模型性能的方法是考虑其预测值平均偏离真实值的程度。这种测量称为平均绝对误差（MAE）。

MAE 的公式如下，其中n表示预测的数量，而e[i]表示预测i的错误：

如其名所示，此公式计算误差绝对值的平均值。由于误差仅仅是预测值与实际值之间的差异，我们可以创建一个简单的MAE()函数，如下所示：

> MAE <- function(actual, predicted) {
    mean(abs(actual - predicted))
} 

我们的预测的 MAE 如下：

> MAE(p.rpart, wine_test$quality) 

[1] 0.5872652 

这意味着，平均而言，我们的模型预测与真实质量评分之间的差异约为0.59。在 0 到 10 的质量尺度上，这似乎表明我们的模型表现相当不错。

另一方面，回想一下，大多数酒既不是非常好也不是非常差；典型的质量评分在 5 到 6 之间。因此，仅预测平均值的分类器根据这个指标也可能表现相当不错。

训练数据中的平均质量评分如下：

> mean(wine_train$quality) 

[1] 5.870933 

如果我们预测每个酒样值为5.87，我们的 MAE 将只有大约0.67：

> MAE(5.87, wine_test$quality) 

[1] 0.6722474 

我们的回归树（MAE = 0.59）平均来说比插补的均值（MAE = 0.67）更接近真实质量评分，但差距不大。相比之下，Cortez 报告了神经网络模型的 MAE 为 0.58，支持向量机的 MAE 为0.45。这表明还有改进的空间。

第 5 步 - 提高模型性能

为了提高学习器的性能，让我们应用模型树算法，这是一种更复杂的将树应用于数值预测的应用。回想一下，模型树通过用回归模型替换叶节点来扩展回归树。这通常比只使用单个数值值进行叶节点预测的回归树产生更准确的结果。

在模型树领域，当前最先进的技术是Cubist算法，该算法本身是对 M5 模型树算法的改进——这两者都是由 J. R. Quinlan 在 20 世纪 90 年代初发表的。尽管实现细节超出了本书的范围，但 Cubist 算法涉及构建决策树，根据树的分支创建决策规则，并在每个叶节点构建回归模型。使用额外的启发式方法，如剪枝和提升，以提高预测质量并平滑预测值的范围。

关于 Cubist 和 M5 算法的更多背景信息，请参阅*《使用连续类学习》，Quinlan, J. R.，1992 年第五次澳大利亚联合人工智能会议论文集，第 343-348 页*。此外，请参阅*《结合实例学习和基于模型的学习》，Quinlan, J. R.，1993 年第十次国际机器学习会议论文集，第 236-243 页*。

Cubist 算法在 R 中通过Cubist包和相关的cubist()函数可用。此函数的语法在以下表格中显示：

图 6.19：模型树语法

我们将使用与回归树不同的语法来拟合 Cubist 模型树，因为cubist()函数不接受 R 公式语法。相反，我们必须指定用于x独立变量和y因变量的数据帧列。由于要预测的葡萄酒质量位于第 12 列，并使用所有其他列作为预测变量，完整的命令如下：

> library(Cubist)
> m.cubist <- cubist(x = wine_train[-12], y = wine_train$quality) 

模型树的基本信息可以通过输入其名称来检查：

> m.cubist 

Call:
cubist.default(x = wine_train[-12], y = wine_train$quality)
Number of samples: 3750
Number of predictors: 11
Number of committees: 1
Number of rules: 25 

在这个输出中，我们看到算法生成了 25 条规则来模拟葡萄酒质量。要检查这些规则中的几个，我们可以对模型对象应用summary()函数。由于完整的树非常大，这里只包括描述第一个决策规则的输出前几行：

> summary(m.cubist) 

 Rule 1: [21 cases, mean 5.0, range 4 to 6, est err 0.5]
   if
        free.sulfur.dioxide > 30
        total.sulfur.dioxide > 195
        total.sulfur.dioxide <= 235
        sulphates > 0.64
        alcohol > 9.1
   then
        outcome = 573.6 + 0.0478 total.sulfur.dioxide
                  - 573 density - 0.788 alcohol
                  + 0.186 residual.sugar - 4.73 volatile.acidity 

你会注意到输出中的if部分与我们在早期构建的回归树有些相似。一系列基于二氧化硫、硫酸盐和酒精的葡萄酒特性的决策创建了一个以最终预测为终点的规则。然而，这个模型树输出与早期回归树输出的一个关键区别是，这里的节点不是以数值预测结束，而是一个线性模型。

该规则的线性模型显示在outcome =语句之后的then输出中。这些数字可以与我们在本章早期构建的多元回归模型完全相同地解释。每个值都是相关特征的估计影响，即该特征单位增加对预测葡萄酒质量的净效应。例如，残留糖的系数为 0.186 意味着残留糖增加一个单位时，预计葡萄酒质量评分将增加 0.186。

重要的是要注意，此模型估计的回归效应仅适用于达到此节点的葡萄酒样本；检查 Cubist 输出的全部内容揭示，在这个模型树中总共构建了 25 个线性模型，每个决策规则一个，每个模型对残留糖和 10 个其他特征的参数估计都不同。

为了检验这个模型的表现，我们将查看它在未见过的测试数据上的表现如何。predict()函数为我们提供了一个预测值的向量：

> p.cubist <- predict(m.cubist, wine_test) 

模型树似乎在预测值范围上比回归树更广：

> summary(p.cubist) 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  3.677   5.416   5.906   5.848   6.238   7.393 

相关性似乎也显著更高：

> cor(p.cubist, wine_test$quality) 

[1] 0.6201015 

此外，该模型略微降低了 MAE：

> MAE(wine_test$quality, p.cubist) 

[1] 0.5339725 

尽管我们没有在回归树模型上取得很大的改进，但我们超越了 Cortez 发布的神经网络模型的性能，并且在使用一个简单得多的学习方法的同时，我们正接近支持向量机模型发布的 MAE 值 0.45。

毫不奇怪，我们证实了预测葡萄酒质量是一个困难的问题；毕竟，品酒本质上是一种主观行为。如果您想进行额外的练习，您可以在阅读第十四章，构建更好的学习者之后，再次尝试这个问题，该章节涵盖了可能带来更好结果的技术。

摘要

在本章中，我们研究了两种建模数值数据的方法。第一种方法，线性回归，涉及将直线拟合到数据上，但一种称为广义线性建模的技术可以使回归适应其他环境。第二种方法使用决策树进行数值预测。后者有两种形式：回归树，它使用叶节点上示例的平均值进行数值预测，以及模型树，它以混合方法在每个叶节点构建回归模型，这种混合方法在某些方面是两者的最佳结合。

我们通过使用回归模型来调查挑战者号航天飞机灾难的原因，开始理解回归模型的有用性。然后，我们使用线性回归模型来计算各种汽车驾驶者的预期保险索赔成本。

由于估计的回归模型很好地记录了特征与目标变量之间的关系，我们能够识别出某些人口统计特征，例如高里程和深夜驾驶者，他们可能需要支付更高的保险费来覆盖他们高于平均的索赔成本。然后，我们将逻辑回归，一种用于二元分类的回归变体，应用于建模保险客户保留的任务。这些例子展示了回归如何灵活地适应许多类型的现实世界问题。

在机器学习的某种不太商业化的应用中，回归树和模型树被用来根据可测量的特征建模葡萄酒的主观质量。在这个过程中，我们了解到回归树提供了一种简单的方法来解释特征与数值结果之间的关系，但更复杂的模型树可能更准确。在这个过程中，我们还学习了新的方法来评估数值模型的表现。

与本章所涵盖的、导致对输入和输出之间关系有清晰理解的机器学习方法截然不同，下一章涵盖了导致几乎无法理解模型的算法。优点是它们是极其强大的技术——属于最强大的股票分类器之一——可以应用于分类和数值预测问题。

加入我们书籍的 Discord 空间

加入我们的 Discord 社区，与志同道合的人相聚，并与其他 4000 多人一起学习：

packt.link/r