R-应用监督学习-二-R 应用监督学习（二）第三章：监督学习简介学习目标到本章结束时，你将能够：解释监督学习和机

R 应用监督学习（二）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第三章：监督学习简介

学习目标

到本章结束时，你将能够：

解释监督学习和机器学习工作流程
使用和探索北京 PM2.5 数据集
解释连续和分类因变量的区别
在 R 中实现基本的回归和分类算法
区分监督学习与其他类型机器学习的关键差异
与监督学习算法的评估指标一起工作
进行模型诊断，以避免系数估计偏差和大的标准误差

在本章中，我们将介绍监督学习，并通过实际案例展示构建机器学习模型的流程。

简介

在前面的章节中，我们探讨了 R 的一些包，如dplyr、plyr、lubridate和ggplot2，其中我们讨论了在 R 中存储和处理数据的基本方法。后来，这些思想被用于探索性数据分析（EDA），以了解如何将数据分解成更小的部分，从数据中提取洞察力，并探索其他更好地理解数据的方法，在尝试高级建模技术之前。

在本章中，我们将进一步介绍机器学习思想。在广泛地为思考机器学习中的各种算法打下基础的同时，我们将详细讨论监督学习。

监督学习基于由领域专家良好标记的数据。对于从图像中分类猫和狗，算法首先需要看到标记为猫和狗的图像，然后根据标签学习特征。大多数拥有大量历史数据的企业的最大受益者是从这些数据中提取的丰富知识。如果数据干净且标注良好，监督学习可以实现高预测精度，这与其他机器学习算法不同，其他机器学习算法通常在开始时会产生较大的误差。在没有正确标签的情况下，从数据中提取任何意义变得困难，除了能够进行探索性分析和聚类之外。

在解决现实世界问题（如预测贷款违约（是/否）、工厂制造机器故障（是/否）、自动驾驶汽车中的目标检测（道路、汽车、信号）、预测股票市场价格（数值））时，标准组件是一组输入（特征）和一个给定的输出（标签），这通常是从历史数据中获得的。当我们预测定量输出时，我们称之为回归，当我们预测定性输出时，我们称之为分类。

北京 PM2.5 数据集概述

在许多国家的城市和农村地区，主要污染物，细颗粒物，是导致人类许多健康风险的原因，同时也影响气候变化。特别是，PM2.5，定义为直径小于 2.5 µm 的气溶胶颗粒，是大气颗粒物的主要类别。各种研究将 PM2.5 与严重的健康问题联系起来，如心脏病和中风。本节中的表格显示了大气颗粒物的类型及其微米级的尺寸分布。

在本章和后续章节中，我们将使用研究论文《评估北京的 PM2.5 污染：严重程度、天气影响、APEC 和冬季供暖》的作者发布的数据集，其中他们使用位于北京 116.47 E，39.95 N 的美国大使馆的每小时 PM2.5 读数，以及从 weather.nocrew.org 获得的每小时气象测量数据，这些数据是在北京首都国际机场（BCIA）获得的。他们的研究声称是首次在中国 PM2.5 污染中结合了 PM2.5 和气象数据。以下表格描述了数据集中的属性：

![图 3.1：北京 PM2.5 数据集的属性。

图 3.1：北京 PM2.5 数据集的属性。

练习 40：探索数据

在这个练习中，我们将通过每个属性的样本值学习数据的结构，并使用summary函数。我们将看到数值变量的五个数字摘要统计量。

执行以下步骤以完成此练习：

首先，使用以下命令将北京 PM2.5 数据集读入 PM25 数据框对象：

PM25 <- read.csv("https://raw.githubusercontent.com/TrainingByPackt/Applied-Supervised-Learning-with-R/master/Lesson03/PRSA_data_2010.1.1-2014.12.31.csv")

接下来，使用str命令打印具有样本值的数据结构：

str(PM25)

上一条命令的输出如下：

'data.frame':	43824 obs. of  13 variables:
 $ No   : int  1 2 3 4 5 6 7 8 9 10 ...
 $ year : int  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
 $ month: int  1 1 1 1 1 1 1 1 1 1 ...
 $ day  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ hour : int  0 1 2 3 4 5 6 7 8 9 ...
 $ pm2.5: int  NA NA NA NA NA NA NA NA NA NA ...
 $ DEWP : int  -21 -21 -21 -21 -20 -19 -19 -19 -19 -20 ...
 $ TEMP : num  -11 -12 -11 -14 -12 -10 -9 -9 -9 -8 ...
 $ PRES : num  1021 1020 1019 1019 1018 ...
 $ cbwd : Factor w/ 4 levels "cv","NE","NW",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ Iws  : num  1.79 4.92 6.71 9.84 12.97 ...
 $ Is   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Ir   : int  0 0 0 0 0 0 0 0 0 0 ...

注意

观察到数据集包含43824个观测值和 13 个属性。观察到数据集包含 2010 年至 2014 年的数据。pm2.5、温度、气压、综合风向、累积风速、累积降雪小时数和累积降雨小时数的值每天每小时都会汇总。

现在，让我们展示数据集的摘要统计量：

summary(PM25)

输出如下：

       No             year          month             day             hour           pm2.5       
 Min.   :    1   Min.   :2010   Min.   : 1.000   Min.   : 1.00   Min.   : 0.00   Min.   :  0.00  
 1st Qu.:10957   1st Qu.:2011   1st Qu.: 4.000   1st Qu.: 8.00   1st Qu.: 5.75   1st Qu.: 29.00  
 Median :21912   Median :2012   Median : 7.000   Median :16.00   Median :11.50   Median : 72.00  
 Mean   :21912   Mean   :2012   Mean   : 6.524   Mean   :15.73   Mean   :11.50   Mean   : 98.61  
 3rd Qu.:32868   3rd Qu.:2013   3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:17.25   3rd Qu.:137.00  
 Max.   :43824   Max.   :2014   Max.   :12.000   Max.   :31.00   Max.   :23.00   Max.   :994.00  
                                                                                 NA's   :2067    
      DEWP              TEMP             PRES      cbwd            Iws               Is          
 Min.   :-40.000   Min.   :-19.00   Min.   : 991   cv: 9387   Min.   :  0.45   Min.   : 0.00000  
 1st Qu.:-10.000   1st Qu.:  2.00   1st Qu.:1008   NE: 4997   1st Qu.:  1.79   1st Qu.: 0.00000  
 Median :  2.000   Median : 14.00   Median :1016   NW:14150   Median :  5.37   Median : 0.00000  
 Mean   :  1.817   Mean   : 12.45   Mean   :1016   SE:15290   Mean   : 23.89   Mean   : 0.05273  
 3rd Qu.: 15.000   3rd Qu.: 23.00   3rd Qu.:1025              3rd Qu.: 21.91   3rd Qu.: 0.00000  
 Max.   : 28.000   Max.   : 42.00   Max.   :1046              Max.   :585.60   Max.   :27.00000  

       Ir         
 Min.   : 0.0000  
 1st Qu.: 0.0000  
 Median : 0.0000  
 Mean   : 0.1949  
 3rd Qu.: 0.0000  
 Max.   :36.0000

以下图像是大气颗粒物尺寸分布（微米）的图形表示：

![图 3.2：大气颗粒物的类型和尺寸分布（微米）。

图 3.2：大气颗粒物的类型和尺寸分布（微米）。

来源：en.wikipedia.org/wiki/File:A…

注意

发表在《胸科疾病杂志》（JTD）上的文章《PM2.5 对人类呼吸系统的影响》的作者讨论了空气污染与呼吸系统疾病之间的关联。他们提供了一种全面的数据驱动方法来解释导致此类呼吸疾病的原因。特别关注北京，研究人员已经广泛研究了 PM2.5 上升的负面影响，并已成为全球各种气候变化论坛的主流讨论点。更多细节可以在文章中找到，链接为 www.ncbi.nlm.nih.gov/pmc/article…

回归问题和分类问题

在我们的日常生活中，我们随处可见分类和回归问题。例如，从 weather.com 获取降雨概率，我们的电子邮件被过滤到垃圾邮件箱和收件箱，我们的个人和家庭贷款被批准或拒绝，决定选择下一个假日目的地，探索购买新房的选项，投资决策以获得短期和长期收益，从亚马逊购买下一本书；这个列表可以一直继续下去。我们周围的世界今天正越来越多地由帮助我们做出选择的算法（这并不总是好事）来运行。

如第二章数据探索分析中所述，我们将使用被称为情境-复杂性-问题（SCQ）的Minto 金字塔原则来定义我们的问题陈述。以下表格展示了针对北京的 PM2.5 问题的 SCQ 方法：

![图 3.3：应用 SCQ 解决北京的 PM2.5 问题。

![图 3.3：应用 SCQ 解决北京的 PM2.5 问题。现在，在前面表格中描述的 SCQ 构造中，我们可以进行简单的相关性分析，以确定影响 PM2.5 水平的因素，或者创建一个预测问题（预测意味着找到一个近似函数，将输入变量映射到输出），使用所有因素估计 PM2.5 水平。为了术语的清晰，我们将把因素称为输入变量。然后，PM2.5 成为因变量（通常被称为输出变量）。因变量可以是分类的或连续的。例如，在电子邮件分类到 SPAM/非 SPAM 的问题中，因变量是分类的。以下表格突出了回归问题和分类问题之间的一些关键区别：![图 3.4：回归问题和分类问题之间的区别。

图 3.4：回归问题和分类问题之间的区别。

机器学习工作流程

为了演示构建预测模型（机器学习或监督学习）的端到端过程，我们创建了一个易于理解的流程。第一步是设计问题，然后是获取和准备数据，这导致为训练和评估编码模型，最后部署模型。在本章的范围内，我们将保持模型解释的简洁，因为它将在第四章和第五章中详细讨论。

下图描述了从准备数据到部署模型所需的整个工作流程：

![图 3.5：机器学习工作流程。图片 C12624_03_05.jpg

图 3.5：机器学习工作流程。

设计问题

一旦确定了工作领域，就会进行问题设计的头脑风暴。首先，将问题定义为回归问题或分类问题。一旦完成，我们选择正确的目标变量，并识别特征。目标变量很重要，因为它决定了训练的方式。监督学习算法将目标变量放在中心，同时试图从给定的特征集中找到模式。

数据来源和准备

数据收集和准备是一项费力的工作，尤其是在数据来源多样且众多的情况下。对于每个数据源，挑战都是不同的，因此处理所需的时间也会有所不同。如果表格数据不包含大量垃圾信息，那么具有表格数据的数据源是最容易处理的，而文本数据由于其自由流动的特性，清理起来最为困难。

编码模型

一旦数据准备就绪，我们就开始选择合适的模型。通常，专家们首先选择一个基线模型，以评估算法使用输入特征和目标变量的可预测性。然后，可以直接尝试最先进的算法，或者决定采用试错法（尝试使用所有可能的模型）。必须理解的是，没有绝对正确或错误的模型，一切取决于数据。在编码过程中，数据被随机分为训练集和测试集。代码被编写来在训练数据集上训练模型，评估则在测试数据上进行，这确保了模型在实际部署时不会表现不佳。

训练和评估

模型评估是模型的重要组成部分，其中决定了其在实际中的可用性。基于给定的一组模型评估指标，我们需要在经过多次尝试和错误后，决定最佳模型。在每次迭代中，计算如 R 平方值、准确率、精确率和 F 分数等指标。通常，整个数据被分为训练数据和测试数据（通常还包括一个用于验证的第三部分）。模型在训练数据上训练，在测试数据上测试。这种分离确保模型不会进行任何机械学习。在更技术性的术语中，模型不会过拟合（关于这一点，请参阅本章的评估指标部分）。通常，在这个工作流程的阶段，一个人可以决定返回并包括更多变量，训练模型，然后重新部署。这个过程会重复进行，直到模型的准确率（或其他重要指标）达到平台期。

我们使用随机数生成函数，如 R 中的sample()函数，将数据随机分割成不同的部分，就像在下一个练习 2 的第 2 步中所做的那样。

练习 41：从北京 PM2.5 数据集中随机生成训练和测试数据集

在这个练习中，我们将从北京 PM2.5 数据集中创建一个随机生成的训练和测试数据集。我们将重用之前练习中创建的PM25对象。

执行以下步骤：

创建一个num_index变量，并将其设置为北京 PM2.5 数据集中观测值的数量：
```
num_index <- nrow(PM25)
```
使用sample()函数，随机选择num_index值的 70%，并将它们存储在train_index中：
```
train_index <- sample(1:num_index, 0.7*nrow(PM25))
```
使用train_index从北京 PM2.5 数据集中选择一个随机子集的行，并将它们存储到一个名为PM25_Train的 DataFrame 中：
```
PM25_Train <- PM25[train_index,]
```
将剩余的观测值存储到一个名为PM25_Test的 DataFrame 中：
```
PM25_Test <- PM25[-train_index,]
```

练习展示了创建训练和测试集的简单示例。随机选择的训练和测试集确保模型没有偏见，并在用于现实世界中的未见数据之前，从所有可能的示例中学习得很好。

部署模型

一旦选定了最佳模型，下一步就是使模型输出能够被业务应用使用。该模型以表示状态转移（REST）API 的形式托管。这些 API 是以端点托管 Web 应用的方式，它监听对模型调用的任何请求，通常返回一个 JSON 对象作为响应。

模型的部署正成为工业界所有机器学习项目的必要部分。一个不可部署的模型对公司来说毫无价值，也许仅仅是为了研发目的。越来越多的专业人士正在专注于模型部署，这有时是一个繁琐且复杂的过程。为了给模型部署应有的重视，我们为其专门设立了一章，即第八章，模型部署。

回归

现在我们已经看到了机器学习的工作流程，我们将探讨两种广泛使用的机器学习算法：回归和分类；两者都采用监督学习来训练模型。本书的整个主题围绕这两种类型的算法展开。北京 PM2.5 数据集将在演示这两种类型时被广泛使用。该数据集有助于理解如何将回归问题转换为分类问题，反之亦然。

简单和多重线性回归

回归是分析学和计量经济学（关注于使用数学方法，尤其是统计学，来描述经济系统的经济学分支）中最有用和最基本的工具之一。在许多方面，现代机器学习的根源在于统计学，这主要归功于弗朗西斯·高尔顿爵士的工作。高尔顿是一位对遗传学、心理学和人类学等领域有深厚兴趣和专长的英国维多利亚时代统计学家和博学家。他是第一个将统计方法应用于研究人类行为和智力的人。值得注意的是，他的出版物《遗传身高回归到中等》基于回归有许多有洞察力的发现。

在本节中，我们将简要分析影响 PM2.5 水平的各种因素，使用北京数据集。特别是，我们将探讨露点、温度、风速和压力等变量对 PM2.5 的影响。

线性回归模型中的假设

由于回归从应用统计学中借鉴了许多概念来建模数据，因此它伴随着许多假设。我们不应将回归算法应用于任何数据集或问题。在我们构建任何模型之前，让我们先检查线性回归的假设。

下表显示了假设以及我们如何从统计上测试线性回归模型是否遵循该假设。该表还显示了一些如果假设被违反的纠正措施。我们将在第四章“回归”中详细讨论这些假设，并执行诊断分析以识别违规情况。

图 3.6：线性回归模型中的各种假设（第一部分）。

图 3.7：线性回归模型中的各种假设（第二部分）。

探索性数据分析（EDA）

构建回归模型需要对目标变量和输入变量之间的模式和关系进行深入分析。北京数据集提供了大量可能影响大气中 PM2.5 水平的环境因素。

练习 42：探索北京 PM2.5 数据集中 PM2.5、DEWP、TEMP 和 PRES 变量的时间序列视图

在这个练习中，我们将可视化pm2.5、DEWP、TEMP和PRES变量在时间序列图中的表现，并观察这些变量在多年中可能出现的任何模式。

执行以下步骤来完成练习：

在系统中导入所有必需的库：

library(dplyr)
library(lubridate)
library(tidyr)
library(grid)
library(ggplot2)

接下来，使用lubridate包的ymd_h函数将年、月和小时转换为日期时间：

PM25$datetime <- with(PM25, ymd_h(sprintf('%04d%02d%02d%02d', year, month, day,hour)))

使用以下命令绘制所有年份的 PM2.5、TEMP、DEWP 和 PRES：

plot_pm25 <- PM25 %>%
  select(datetime, pm2.5) %>%
  na.omit() %>%
  ggplot() + 
  geom_point(aes(x = datetime, y = pm2.5), size = 0.5, alpha = 0.75) +
  ylab("PM2.5")
plot_TEMP <- PM25 %>%
  select(datetime, TEMP) %>%
  na.omit() %>%
  ggplot() + 
  geom_point(aes(x = datetime, y = TEMP), size = 0.5, alpha = 0.75) +
  ylab("TEMP")
plot_DEWP <- PM25 %>%
  select(datetime, DEWP) %>%
  na.omit() %>%
  ggplot() + 
  geom_point(aes(x = datetime, y = DEWP), size = 0.5, alpha = 0.75) +
  ylab("DEWP")
plot_PRES <- PM25 %>%
  select(datetime, PRES) %>%
  na.omit() %>%
  ggplot() + 
  geom_point(aes(x = datetime, y = PRES), size = 0.5, alpha = 0.75) +
  ylab("PRES")

现在，使用以下命令来绘制图表：

grid.newpage()
grid.draw(rbind(ggplotGrob(plot_pm25), ggplotGrob(plot_TEMP),ggplotGrob(plot_DEWP),ggplotGrob(plot_PRES), size = "last"))

图如下：

![图 3.8：散点图显示了环境因素（如温度、露点和压力）的趋势和季节性，以及 2010 年至 2014 年底北京地区的 PM2.5 水平。图 C12624_03_08.jpg

图 3.8：散点图显示了环境因素（如温度、露点和压力）的趋势和季节性，以及 2010 年至 2014 年底北京地区的 PM2.5 水平。

在这个练习中，我们首先展示了数据集中PM2.5、DEWP、TEMP和PRES变量的时间序列视图，并观察其模式。如图图 3.8所示，观察到DEWP、TEMP和PRES具有明显的季节性（每 12 个月重复相同的模式），而 PM2.5 似乎具有随机模式。这是早期迹象，表明我们不太可能看到三个变量对 PM2.5 有任何影响。然而，让我们进一步探究这个假设，使用相关性图并观察变量之间是否存在任何关系。

练习 43：进行相关性分析

在这个练习中，我们将进行相关性分析，研究各种因素之间关系的强度。

执行以下步骤来完成练习：

使用以下命令将corrplot包导入到系统中：
```
library(corrplot)
```

现在，创建一个新的对象并将PM25中的所需值存储在其中：

corr <- cor(PM25[!is.na(PM25$pm2.5),c("pm2.5","DEWP","TEMP","PRES","Iws","Is","Ir")])

使用corrplot包来显示相关性矩阵的图形表示：
```
corrplot(corr)
```
图如下：

![图 3.9：北京数据集中所有变量对的关联性。图 C12624_03_09.jpg

图 3.9：北京数据集中所有变量对的关联性。

首先，我们计算所有变量之间的相关性。生成的相关性图显示，PM2.5 与其他变量之间似乎没有强烈的关联。然而，PM2.5与DEWP、TEMP和Iws显示出一些轻微的相关性，这表明存在某种关系。这并不令人惊讶，因为我们已经在图 3.8中看到，虽然三个变量遵循季节性趋势，但 PM2.5 似乎更随机。请注意，我们对数据集没有进行任何处理或转换；这些发现直接来自我们的第一层分析。我们将在第四章，回归中详细说明。现在，让我们也使用散点图来可视化变量之间的关系。

练习 44：绘制散点图以探索 PM2.5 水平与其他因素之间的关系

在这个练习中，我们将使用散点图来探索 pm2.5 水平与其他因素之间的关系。我们希望看到是否会出现任何有趣的模式或关系。散点图是探索变量之间关系的简单而有效的可视化工具。

执行以下步骤来完成练习：

将 ggplot2 包导入到您的系统中：
```
library(ggplot2)
```
以 month 变量作为颜色绘制 DEWP 和 PM2.5 之间的散点图：
```
ggplot(data = PM25, aes(x = DEWP, y = pm2.5, color = month)) +  geom_point() +  geom_smooth(method='auto',formula=y~x, colour = "red", size =1.5)
```
散点图如下：

![图 3.10：显示 DEWP 和 PM2.5 水平之间关系的散点图。

图 3.10：显示 DEWP 和 PM2.5 水平之间关系的散点图。
以 month 变量作为颜色绘制 TEMP 和 PM2.5 之间的散点图：
```
ggplot(data = PM25, aes(x = TEMP, y = pm2.5, color = month)) +  geom_point() +  geom_smooth(method='auto',formula=y~x, colour = "red", size =1.5)
```
散点图如下：

![图 3.11：显示 TEMP 和 PM2.5 水平之间关系的散点图。

图 3.11：显示 TEMP 和 PM2.5 水平之间关系的散点图。

创建一个以一天中的小时为颜色，按月份分别显示的 DEWP 和 PM2.5 之间的散点图：

ggplot(data = PM25, aes(x = DEWP, y = pm2.5, color = hour)) +  geom_point() +  geom_smooth(method='auto',formula=y~x, colour = "red", size =1) +  facet_wrap(~ month, nrow = 4)

散点图如下：

![图 3.12：按年份月份拆分的 DEWP 和 PM2.5 关系的散点图。 img/C12624_03_12.jpg

图 3.12：按年份月份拆分的 DEWP 和 PM2.5 关系的散点图。

为了评估变量之间的一些关系，我们使用 PM2.5 和 DEWP 之间的散点图并添加了拟合线。观察代码中，我们向 geom_smooth() 传递了一个参数，即 method = "auto"，它会根据数据自动决定使用哪种模型来拟合直线。如图 图 3.10 所示，直线不是线性的。geom_smooth 方法选择了 TEMP 和 PM2.5 的绘图，如图 图 3.11 所示。然而，我们可以更进一步，按月份拆分散点图，如图 图 3.12 所示。这表明存在线性关系，但它高度依赖于季节。例如，在四月（用整数 4 表示），DEWP 和 PM2.5 有一个近乎完美的直线拟合。我们将在 第四章，回归中进一步讨论这个话题。

因此，我们已经看到了一些假设的违反和环境污染因素与 PM2.5 之间缺乏强相关性的情况。然而，似乎还有进一步审查的空间。在本章关于监督学习的介绍中，我们只关注基于我们的机器学习工作流程的方法。

备注

要了解更多关于 GAM 的信息，请参阅此文档：www.stat.cmu.edu/~cshalizi/u…

活动 5：按月份绘制 PRES 和 PM2.5 的散点图

在这个活动中，我们将创建DWEP和PM2.5之间的散点图。通过这个活动，我们将学习使用facet_wrap()函数在ggplot()之上创建一个层，将散点图的可视化分割到每个月，从而帮助观察任何季节性模式。

执行以下步骤来完成活动：

在ggplot中，使用PRES变量分配a()方法的组件。
在geom_smooth()方法的下一层，将colour = "blue"设置为区分。
最后，在facet_wrap()层中，使用month变量为每个月绘制单独的隔离图。

图表如下：

![图 3.13：显示 PRES 和 PM2.5 之间关系的散点图。

图 3.13：显示 PRES 和 PM2.5 之间关系的散点图。

注意

该活动的解决方案可以在第 445 页找到。

模型构建

我们简要探讨了PM2.5与一些因素（如TEMP和DEWP）之间的关系。同样的分析可以应用于其他变量，如PRES、Iwd等。在本节中，让我们创建一个线性模型。（即使我们知道模型的选择不是最好的，我们也不会犹豫去运行模型。机器学习中的试错法总是建立事实的最佳方式。）

通常，线性回归模型输入变量（自变量）和目标变量（因变量或解释变量）之间的线性关系。如果我们有一个解释变量，它被称为简单线性回归；如果有多个解释变量，则称为多元线性回归。以下方程是线性回归或线性预测函数的数学表示，其中包含p个解释变量和n个观测值：

这里，每个包含家具、桌子、座椅、凳子自动生成的描述是一个列值的向量（解释变量）对于包含家具的图片自动生成的描述是未知的参数或系数。包含家具、座椅自动生成的描述，这使得这个方程适合简单的线性回归。有许多算法可以将这个函数拟合到数据上。最流行的一个是普通最小二乘法（OLS）。我们将在下一章关于回归的章节中详细讨论 OLS。

另一种思考方式是包含家具的图片自动生成的描述是它是一个线性预测函数，尽可能地将观测值拟合到维度的空间中，最小化残差平方和（目标值的实际值与预测值之间的差异）。

在接下来的练习中，我们将跳过将数据集分为训练集和测试集的步骤，因为我们仍然处于探索阶段，尚未决定正式进行建模练习。（我们将在下一章中涉及这一点。）我们将使用 R 中的 lm() 方法构建线性模型。同样，关于这一点将在下一章中详细介绍。目前，只需注意 lm() 方法使用一个或多个输入变量将目标变量拟合到一条直线。在简单线性回归中，我们只使用一个变量来拟合直线，而在多重线性回归中，我们可以使用多个变量。

练习 45：探索简单和多重回归模型

在这个练习中，我们将探索简单和多重回归模型。

执行以下步骤以完成练习：

将所需的库和包导入 R-Studio。
接下来，创建一个名为 simple_PM25_linear_model 的 DataFrame 对象，并使用 lm() 方法构建线性模型：
```
simple_PM25_linear_model <- lm(pm2.5 ~ DEWP, data = PM25)
```

使用如下所示的方法，使用 summary 方法打印对象的摘要：

summary(simple_PM25_linear_model)

输出如下：

Call:
lm(formula = pm2.5 ~ DEWP, data = PM25)
Residuals:
    Min      1Q  Median      3Q     Max 
-115.47  -61.26  -28.75   33.83  923.54 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 96.69984    0.44705  216.31   <2e-16 ***
DEWP         1.09325    0.03075   35.55   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 90.69 on 41755 degrees of freedom
  (2067 observations deleted due to missingness)
Multiple R-squared:  0.02939,	Adjusted R-squared:  0.02936 
F-statistic:  1264 on 1 and 41755 DF,  p-value: < 2.2e-16

接下来，创建另一个 DataFrame 对象，并使用 lm() 方法构建线性模型：
```
multiple_PM25_linear_model <- lm(pm2.5 ~ DEWP+TEMP+Iws, data = PM25)
```

使用 summary 函数打印模型对象的摘要：

summary(multiple_PM25_linear_model)

输出如下：

A)
______________________________________________________________
Call:
lm(formula = pm2.5 ~ DEWP + TEMP + Iws, data = PM25)
Residuals:
    Min      1Q  Median      3Q     Max 
-149.02  -53.74  -16.61   34.14  877.82 
______________________________________________________________
B)
______________________________________________________________
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 161.151207   0.768727  209.63   <2e-16 ***
DEWP          4.384196   0.051159   85.70   <2e-16 ***
TEMP         -5.133511   0.058646  -87.53   <2e-16 ***
Iws          -0.274337   0.008532  -32.15   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
______________________________________________________________
C)
______________________________________________________________
Residual standard error: 81.51 on 41753 degrees of freedom
  (2067 observations deleted due to missingness)
Multiple R-squared:  0.216,	Adjusted R-squared:  0.2159 
F-statistic:  3834 on 3 and 41753 DF,  p-value: < 2.2e-16
______________________________________________________________

模型解释

现在，基于简单和多重线性回归模型的先前输出，让我们尝试理解输出中的每一部分代表什么。在本书的这个阶段，了解每一部分的意义就足够了；我们将在 第四章，回归中讨论结果。

使用 lm() 方法与因变量和自变量，用 ~ 符号表示的公式形式表示。这类似于我们的线性预测函数。在简单的回归模型中，只有一个变量——DEWP——而在多重模型中，有 DEWP、TEMP 和 Iws。您还可以看到残差的五个汇总统计量（最小值、第一四分位数、中位数、第三四分位数和最大值）。这表明预测值与实际值之间的差距。
将 X_j 部分纳入我们的预测方程中，我们将得到预测值。名为 Std. Error 的列是估计的标准误差。t 值是通过将 Estimate 和 Std. Error 的比率得到的，而 p 值突出了估计的统计显著性。视觉线索，即 * 和 . 符号是基于 p 值的。小于 0.001 的值得到三个星号，而介于 0.1 和 0.05 之间的值得到一个 .（点）。三个星号意味着最佳情况，即对应于自变量的估计是显著的并且对预测（或解释）因变量是有用的。换句话说，p 值有助于确定回归模型相对于零模型（仅因变量的均值）的显著性。
部分 C：

这一部分展示了模型的功效。最重要的观察值是 R 平方和调整 R 平方值，这些是表示回归模型中由独立变量（s）解释的因变量变化的百分比统计量。

在本章中查看关于评估指标的章节，以了解模型在 R 平方和调整 R 平方指标上的表现解释。

分类

与回归算法类似，分类也通过依赖或目标变量进行学习，并使用所有预测或独立变量来找到正确的模式。主要区别在于，在分类中，目标变量是分类的，而在回归中，它是数字的。在本节中，我们将通过使用北京 PM2.5 数据集来介绍逻辑回归，以展示这一概念。

逻辑回归

逻辑回归是用于二元分类的最受欢迎的透明模型。透明模型定义为提供对预测中整个推理过程的可见性的模型。对于每个做出的预测，我们可以利用模型的数学方程式并解码预测的原因。还有一些完全属于黑盒的分类模型，也就是说，我们无论如何都无法理解模型利用的预测推理。在我们只想关注最终结果的情况下，我们应该更喜欢黑盒模型，因为它们更强大。

简要介绍

虽然名称以回归结尾，但逻辑回归是一种用于预测二元分类结果的技巧，因此是分类问题的良好选择。如前节所述，我们需要一种不同的方法来对分类结果进行建模。这可以通过将结果转换为优势比的对数或事件发生的概率来实现。

让我们将这种方法提炼成更简单的结构。假设一个事件成功的概率为 0.7。那么，同一事件失败的概率将被定义为1 – 0.7 = 0.3。成功的优势被定义为成功概率与失败概率的比率。成功的优势将是0.7/0.3 = 2.33，即成功的优势是 2 比 1。如果成功的概率是 0.5，即 50-50 的机会，那么成功的优势是 1 比 1。逻辑回归模型可以用以下方式数学表示：

![图片 C12624_03_27.jpg]

这里，一个人的插图自动生成的描述是优势比的对数，也称为logit函数。进一步解决数学问题，我们可以推导出结果概率，如下所示：

![图片 C12624_03_29.jpg]

讨论方程的数学背景和推导超出了本章的范围。然而，为了总结，logit 函数（或逻辑函数），作为连接函数，帮助逻辑回归直观地将预测结果（预测结果）重新表述为优势比的对数。这当解决时，有助于我们预测二元因变量的概率。

逻辑回归的机制

就像线性回归一样，其中变量的 beta 系数使用 OLS 方法估计，逻辑回归模型利用最大似然估计（MLE）方法。MLE 函数估计模型参数或 beta 系数的最佳值集，以最大化似然函数，即概率估计。它也可以定义为所选模型与观察数据之间的一致性。当估计出最佳参数值集时，将这些值/β系数插入到前面定义的模型方程中，有助于估计给定样本的输出概率。类似于 OLS，MLE 是一个迭代过程。

模型构建

就像在 R 中使用线性回归构建逻辑回归模型一样，我们有glm()广义线性模型方法来拟合数据和 logit 函数来评分观测值。

使用 glm()函数的语法如下：

glm(Y ~ X1 + X2 + X3, data = <train_data>,family=binomial(link='logit'))

在这里，Y 是我们的因变量，X1、X2 和 X3 是自变量。data 参数将采用训练数据集。family 参数设置为 binomial(link='logit')，这适合逻辑回归模型。

练习 46：将北京 PM2.5 数据集中的滚动 3 小时平均值存储

在这个练习中，我们将创建一个新变量，该变量存储北京 PM2.5 数据集中 PM2.5 变量的滚动 3 小时平均值。滚动平均值将平滑 PM2.5 读数中的任何噪声。

让我们使用zoo包中的rollapply方法来完成练习：

将年、月、日和小时合并成一个名为datetime的新变量：

PM25$datetime <- with(PM25, ymd_h(sprintf('%04d%02d%02d%02d', year, month, day,hour)))

删除 NAs 并查看 PM2.5 数据集中pm2.5变量的前 6 个值：

PM25_subset <- na.omit(PM25[,c("datetime","pm2.5")])
head(PM25_subset$pm2.5)

输出如下：

[1] 129 148 159 181 138 109

将PM25_subset存储到以 datetime 为索引的有序观测值的zoo对象中，并打印前 6 个值：

zoo(PM25_subset$pm2.5,PM25_subset$datetime)

输出如下：

2010-01-02 00:00:00 2010-01-02 01:00:00 2010-01-02 02:00:00 
                129                 148                 159 
2010-01-02 03:00:00 2010-01-02 04:00:00 2010-01-02 05:00:00 
                181                 138                 109

使用rollapply函数创建pm2.5变量的 3 小时滚动平均值，并打印前 6 个值：

PM25_three_hour_pm25_avg <- rollapply(zoo(PM25_subset$pm2.5,PM25_subset$datetime), 3, mean)

输出如下：

2010-01-02 01:00:00 2010-01-02 02:00:00 2010-01-02 03:00:00 
           145.3333            162.6667            159.3333 
2010-01-02 04:00:00 2010-01-02 05:00:00 2010-01-02 06:00:00 
           142.6667            117.3333            112.6667

注意到145.33值是三个小时的pm2.5变量的平均值，如步骤 3 所示（129、148和159）。

活动 6：转换变量并推导新变量以构建模型

在这个活动中，我们将在构建模型之前执行一系列转换并推导新变量。我们需要将pm2.5变量转换为分类变量，以便应用逻辑回归模型。

在我们构建逻辑回归分类模型之前，需要执行以下步骤：

将年、月、日和小时合并成一个名为 datetime 的新变量。
使用 datetime 变量，计算 3 小时窗口内 pm2.5 值的平均值。将这个新变量命名为 PM25_three_hour_pm25_avg。
创建一个名为 pollution_level 的二进制变量。如果 PM25_three_hour_pm25_avg 大于 35，则其值为 1，否则为 0。
以 pollution_level 作为因变量，构建逻辑回归模型。

打印模型的摘要。

最终输出如下：

Call:
glm(formula = pollution_level ~ DEWP + TEMP + Iws, family = binomial(link = "logit"), 
    data = PM25_for_class)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4699  -0.5212   0.4569   0.6508   3.5824  
Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.5240276  0.0273353   92.34   <2e-16 ***
DEWP         0.1231959  0.0016856   73.09   <2e-16 ***
TEMP        -0.1028211  0.0018447  -55.74   <2e-16 ***
Iws         -0.0127037  0.0003535  -35.94   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 49475  on 41754  degrees of freedom
Residual deviance: 37821  on 41751  degrees of freedom
AIC: 37829
Number of Fisher Scoring iterations: 5

注意

这个活动的解决方案可以在第 446 页找到。

解释模型

大部分 glm() 输出看起来与 lm() 方法相似，但有一些新的值，如下所示：

零偏差
残差偏差
赤池信息量准则（AIC）
费舍尔评分

为了避免评分，所有上述措施将在 第五章，分类中详细描述。

参考本章下一节关于 评估指标（基于 混淆矩阵指标 部分）的内容，以找到对模型在 R 平方和调整 R 平方指标上表现如何的解释。

评估指标

在本节中，我们将介绍所有用于评估机器学习模型预测质量的评估措施。根据因变量，我们有几种评估措施的选择。在我们机器学习工作流程的训练和评估步骤中，我们提到，在我们得到期望的结果之前，我们通过添加新变量或更改参数来迭代训练模型。在每次迭代中，我们尝试优化任何一个或两个评估指标。以下表格总结了用于回归、分类和推荐系统的各种类型指标。鉴于本书的范围，我们将更深入地探讨回归和分类算法的细节：

图 3.14：各种机器学习算法的指标

图片 C12624_03_14.jpg

图 3.14：各种机器学习算法的指标。

均方误差 (MAE)

绝对误差是无方向的，这意味着它并不重要，模型在测试数据集上预测的因变量的值是小于还是大于实际值。因此，在我们的北京 PM2.5 数据集示例中，MAE 将给出 PM2.5 预测的平均绝对误差（因变量的预测值和实际值之间的差异），不考虑误差的方向（正或负）：

图片 C12624_03_30.jpg

这里，图片 C12624_03_31.png 是因变量的第 i 个观测值，图片 C12624_03_32.png 是预测值或期望值。

均方根误差 (RMSE)

与 MAE 类似，均方根误差也计算平均预测误差。然而，它基于二次评分，计算平均平方误差的平方根。此外，与 MAE 不同，MAE 取预测值和实际值之间的绝对差，而 RMSE 取平方，这增加了高误差值在开平方前的权重：

在这里，表示第 i 个观测的因变量实际值和估计值之间的差异。

R 平方

R 平方衡量的是线性模型解释响应变量变异的百分比（介于 0 和 1 之间或从 0%到 100%的值）。换句话说，它衡量的是输入特征解释的变异。0%的 R 平方意味着模型的输入特征对响应变量没有任何解释作用。接近 100%意味着该模型是响应变量的良好预测器。例如，如果我们想预测某个地区的房价，特征如卧室数量、面积（平方英尺）以及学校和市场附近的距离决定了房产的价值。然而，仅凭 R 平方不能用于评估模型的好坏。还需要进行关于残差、正态性和异方差性的各种诊断检查。我们将在第四章，回归中详细讨论这一点。

在这里，是因变量实际值和估计值之间平方差的和，而表示因变量实际值和平均值之间平方差的和。

调整后的 R 平方

当我们在回归模型中添加新变量时，随着新变量在解释因变量变化方面的贡献增加，模型的 R 平方值也会提高。（如果新变量设计不当且与解释因变量不相关，则会出现反论点。）因此，为了使评估指标对变量数量保持无偏见，我们在计算中引入了n（观测数）和q（变量数），从而惩罚 R 平方值。这被称为调整后的 R 平方，它同时调整了观测数和变量数。在处理多元线性回归时，查看调整后的 R 平方是一个好的做法。

MSE（均方误差）:

在这里，n是观测数，而q是模型中的系数数。

MST（总均方误差）:

平均倒数排名（MRR）

MRR 常用于评估搜索引擎、推荐算法以及数字空间中的许多其他信息检索算法。MRR 易于解释。通常，它可以用来评估为输入生成响应列表的算法。例如，您在 Google 查询中看到的搜索结果和您在 Amazon 上看到的商品推荐。以下表格显示了计算倒数排名的示例。MRR 的范围从 0 到 1；值越接近 1，表明算法在列表顶部给出了相关结果。

图 3.15：计算倒数排名的示例。

练习 47：寻找评估指标

在这个练习中，我们将找到 MAE、RMSE、R-squared、调整后的 R-squared 和 MRR。

执行以下步骤：

导入所需的库和包。
创建一个名为y_predicted的变量，并将multiple_PM25_linear_model的值赋给它：
```
y_predicted <- predict(multiple_PM25_linear_model, data = PM25)
```
使用以下命令将PM25数据集的值赋给变量：
```
y_actual <- PM25[!is.na(PM25$pm2.5),"pm2.5"]
```

使用均值函数找到 MAE：

MAE <- mean(abs(y_actual - y_predicted))

输出如下：

## [1] 59.82112

接下来，计算 RMSE：

RMSE <- sqrt(mean((y_actual - y_predicted)²))

输出如下：

## [1] 82.09164

现在，使用以下命令计算 R-squared 值：

model_summary <- summary(multiple_PM25_linear_model)
model_summary$r.squared

输出如下：

## [1] 0.216

接下来，使用以下命令找到调整后的 R-squared：
```
model_summary$adj.r.squared
```
输出如下：
```
## [1] 0.2159
```

最后，使用以下命令找到 MRR：

Query_RR_Vector <- c(1/3,1/4,1)
MRR <- sum(Query_RR_Vector)/length(Query_RR_Vector)

输出如下：

## [1] 0.5277778

观察到 MAE 的值为59.82，RMSE 为82.09，这表明误差有很高的方差。换句话说，预测中的观测值有很高的误差（这增加了误差幅度频率分布的方差）；MAE 未能识别误差，而 RMSE 很好地放大了它。如果 MAE 和 RMSE 几乎相等，我们可以推断误差幅度频率分布的方差很低，并且模型对所有观测值的表现都很好。

基于混淆矩阵的指标

基于混淆矩阵的指标用于分类算法。可以从混淆矩阵（也称为A和B）中推导出一系列指标。否则，关于目标变量没有负面或正面的东西。列联表也可以是 NxN，其中N是响应变量中的类别或类别的数量。例如，如果我们想在一个给定的图像中分类 26 个手写英文字符，我们需要一个 26x26 的矩阵：

图 3.16：混淆矩阵的元素。

如果我们将TP、TN、FP和FN排列在一个 2x2 的列联表中，我们就可以得到混淆矩阵，如下表所示：

图 3.17：混淆矩阵。

准确度

准确率衡量模型对正负样本的正确整体分类。矩阵对角线元素（TP 和 TN）之和除以正负样本的总数给出准确率。在现实场景中，准确率并不总是可靠的指标。考虑一下，我们想要区分癌症 CT 扫描和良性 CT 扫描。显然，我们可能会有许多阴性扫描和很少的阳性扫描。这导致了我们所说的不平衡数据集。如果模型主要能够准确预测良性扫描，但在预测癌症 CT 扫描时产生重大错误，准确率可能仍然很高，但模型并不那么有用。

灵敏度

为了解决我们之前讨论的准确率问题，我们可以使用灵敏度（也称为召回率、命中率或真正例率****TPR）和特异性（下一节讨论）的组合。灵敏度给出了模型对阳性案例的预测能力（例如，在 CT 扫描中检测癌症）。我们从所有真正例（TP）案例与阳性（P）案例的数量之比中获得灵敏度。

特异性

特异性提供了对负例正确预测的定量评估（例如，检测良性 CT 扫描）。我们从真阴性案例数与负案例数之比中获得灵敏度。

高灵敏度和特异性值表示模型质量较好。在大多数情况下，我们试图平衡这两个指标以获得最佳模型。

F1 分数

F1 分数通过取两个指标（精确度和灵敏度）的调和平均值（适用于取两个或更多率的平均值）来结合精确度和灵敏度，如下公式所述。阳性预测值（PPV或精确度）衡量的是正确预测的数量与真阳性数和假阳性数之和的比值，即所有阳性病例预测中正确的是多少。

F1 分数比准确率更稳健，但在不平衡类别的情况下仍然存在问题。

对于评估分类模型的好坏，没有好或坏的指标。机器学习从业者通常会查看许多指标的组合来得出模型好坏的结论。这就是为什么了解如何解释上述讨论的每个指标变得很重要。

练习 48：在训练数据上使用模型评估

在这个练习中，我们将使用caret包中的confusionMatrix函数对训练数据进行模型评估。该函数会打印出准确率、灵敏度、特异性和许多其他指标。

执行以下步骤以完成练习：

将所需的库和包导入系统。
创建一个名为predicted的变量，并赋予其值，如图所示：
```
predicted <- ifelse(PM25_logit_model$fitted.values>0.5, 1,0)
```
接下来，创建另一个名为actual的变量，如图所示：
```
actual <- PM25_for_class$pollution_level
```
导入 caret 库：
```
library(caret)
```

最后，使用confusionMatrix方法来描述分类模型的表现：

confusionMatrix(predicted, actual)

输出如下：

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     0     1
##          0  5437  2097
##          1  6232 27989
##                                           
##                Accuracy : 0.8005          
##                  95% CI : (0.7967, 0.8044)
##     No Information Rate : 0.7205          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.4444          
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.4659          
##             Specificity : 0.9303          
##          Pos Pred Value : 0.7217          
##          Neg Pred Value : 0.8179          
##              Prevalence : 0.2795          
##          Detection Rate : 0.1302          
##    Detection Prevalence : 0.1804          
##       Balanced Accuracy : 0.6981          
##                                           
##        'Positive' Class : 0

confusionMatric()输出结果中显示的许多指标在本节中都有描述。然而，在阅读详细内容之前，这里有一个简要的总结。这个逻辑回归模型的准确率为 80%，按照标准来看是好的。这表明我们可以用其他环境因素以 80%的准确率预测正常和高于正常的 PM2.5 值。然而，请注意，准确率是基于整个训练数据集的。我们没有将数据分成两部分来检查过拟合情况，这是一种模型在训练数据上测试时表现非常好，但在测试（或未见）数据上表现较差的情况。

敏感性为 46%，特异性为 93%。这意味着模型在处理负例（高于正常 PM2.5）方面表现良好。通常，这两个指标之间必须有一个权衡。然而，在这种情况下，模型的优先级是尽可能多地预测高于正常的状态。因此，一旦我们有了混淆矩阵，高特异性是可取的；我们可以从中计算出所有指标。

受试者工作特征（ROC）曲线

在分类模型的背景下，预测的输出是一个定量估计，通常是一个概率度量。在二元逻辑回归中，将一个观测值分类为另一个观测值（例如，垃圾邮件与非垃圾邮件）的常用阈值是 0.5。这意味着如果概率大于 0.5，则将其分类为垃圾邮件，否则为非垃圾邮件。现在，根据阈值的不同，你将在我们之前讨论的混淆矩阵中得到不同的 TP、TN、FP 和 FN 值。虽然查看给定阈值（通常是 0.5）下的混淆矩阵是一种标准做法，但它可能无法完全展示模型在现实世界中的表现，这就是为什么阈值的选择至关重要。

ROC 曲线是一种优雅的可视化方式，显示了在每一个可能的阈值下，真正例率（通常由敏感性表示）和真正例率（通常由特异性表示）之间的变化。它帮助我们确定分类的正确阈值。此外，ROC 曲线下的面积（称为 AUC），其值在 0 到 1 之间，告诉我们模型的好坏。越接近 1，意味着模型能够成功地将大多数观测值分类为正类和负类。

在 R 中使用 ROCR 包，我们将使用逻辑回归获取 PM2.5 预测的 ROC 曲线。我们将在下一个练习中观察 AUC。

练习 49：创建 ROC 曲线

在这个练习中，我们将使用 ROCR 包来获取 ROC 曲线。

执行以下步骤：

使用以下命令将 ROCR 包导入系统：
```
library(ROCR)
```

接下来，定义 pred1 和 pref1 对象：

pred1 <- prediction(predict(PM25_logit_model), PM25_for_class$pollution_level)
perf1 <- performance(pred1,"tpr","fpr")

接下来，使用以下命令找到 AUC：

auc <- performance(pred1,"auc")
as.numeric(auc@y.values)

输出如下：

## [1] 0.8077673

使用 plot 命令绘制图形：
```
plot(perf1)
```
![图 3.18：真阳性率（灵敏度）与假阳性率（特异性）之间的 ROC 曲线。

图 3.18：真阳性率（灵敏度）与假阳性率（特异性）之间的 ROC 曲线。

摘要

在本章中，我们从设计问题开始，逐步介绍了构建机器学习工作流程的过程，包括部署模型。我们简要讨论了简单回归、多重回归和逻辑回归，以及所有必要的评估指标，以解释和判断模型的性能。这两个算法分别展示了回归和分类问题的监督学习。

在本章中，我们使用了北京 PM2.5 数据集来构建模型。在这个过程中，我们通过简单地重新设计因变量，将回归问题转换为分类问题。这种重新设计通常用于现实世界问题，以适应特定的用例。

下一章，我们将深入探讨回归算法的细节，并详细阐述除线性回归之外的多种回归算法类型，并讨论何时使用哪种算法。

第四章：回归

学习目标

到本章结束时，你将能够：

构建回归问题
实现各种类型的回归方法及其用例
分析和选择正确的回归方法
通过回归的视角连接统计学和机器学习
深入探讨模型诊断

在本章中，我们将关注各种类型的回归，以及何时使用哪种类型，同时结合 R 语言中的演示。

简介

在上一章中，我们了解了线性回归模型以及输入变量（自变量）和目标变量（因变量或解释变量）之间的线性关系。如果只使用一个变量作为自变量，则定义为简单线性回归。如果使用多个解释变量（自变量），则称为多重线性回归。

回归算法和问题基于预测一个数值目标变量（通常称为因变量），给定所有输入变量（通常称为自变量），例如，根据位置、面积、靠近购物中心等因素预测房价。许多回归的概念都源自统计学。

机器学习的整个领域现在是一个数学、统计学和计算机科学的完美平衡。在本章中，我们将使用回归技术来了解如何建立输入和目标变量之间的关系。我们还将强调模型诊断，因为回归充满了假设，在模型应用于现实世界之前需要检查这些假设。

本质上，所有模型都是错误的，但有些是有用的。 —— 乔治·博克斯

在第三章“监督学习简介”中，我们简要介绍了简单和多重线性回归。在本章中，我们将更专注于模型诊断和其他类型的回归算法，以及它与线性回归的不同之处。

线性回归

让我们回顾第三章“监督学习简介”中的多重线性回归。以下方程是具有p个解释变量和n个观测值的线性方程或线性预测函数的数学表示：

其中每个一张包含家具、桌子的图片自动生成的描述是一个列值的向量（解释变量）和一张包含家具的图片自动生成的描述是未知参数或系数。一张包含家具、座椅的图片自动生成的描述，使这个方程适用于简单线性回归。有许多算法可以将此函数拟合到数据上。最流行的一个是普通最小二乘法（OLS）。

在理解 OLS 的细节之前，首先让我们解释一下在尝试从第三章，监督学习简介中的简单和多元线性回归的模型构建部分拟合北京 PM2.5 数据时得到的方程：

![img/C12624_04_26.jpg]

如果我们代入回归系数的值，一张人脸的插图描述自动生成和![img/C12624_04_28.png]从lm()函数的输出中，我们得到：

![img/C12624_04_29.jpg]

前面的方程试图回答“DEWP、TEMP和Iws这些因素对于预测pm2.5水平是否重要？”的问题。

该模型估计了平均而言，DEWP、TEMP和Iws值如何影响pm2.5水平。例如，DEWP增加一个单位，pm2.5值将增加4.384196。这就是我们通常将这些系数称为权重的原因。需要注意的是，如果R-squared 值低，这些估计的系数是不可靠的。

练习 50：使用`multiple_PM_25_linear_model`对象打印系数和残差值

在这个练习中，我们将使用multiple_PM25_linear_model对象打印系数和残差值。

执行以下步骤以完成练习：

使用multiple_PM25_linear_model对象上的$运算符提取属性系数：

multiple_PM25_linear_model$coefficients

输出如下：

(Intercept)        DEWP        TEMP         Iws 
161.1512066   4.3841960  -5.1335111  -0.2743375

使用multiple_PM25_linear_model对象上的$运算符提取属性残差：
```
multiple_PM25_linear_model$residuals
```

输出如下：

25            26            27            28 
  17.95294914   32.81291348   21.38677872   26.34105878 
           29            30            31            32

活动 7：不使用汇总函数使用模型对象打印各种属性

在第三章的多元线性回归模型部分，监督学习简介中，我们创建了一个多元线性回归模型，并将其存储在模型对象multiple_PM25_linear_model中，使用模型对象。

此活动有助于理解在模型构建后如何提取一些重要的模型属性。在少数情况下，我们将使用$运算符，在其他情况下，我们将进行一些简单的计算。使用multiple_PM25_linear_model对象打印以下模型属性：

残差
拟合值
R-Ssquared 值
F 统计量
系数 p 值

让我们使用模型对象打印这些值：

首先，打印系数值。确保输出类似于使用coefficients选项的summary函数的输出。这些系数是使用 OLS 算法的模型得到的拟合值：
```
(Intercept)        DEWP        TEMP         Iws 
161.1512066   4.3841960  -5.1335111  -0.2743375
```

找到预测值和实际值之间的残差（差异），应尽可能小。残差反映了使用系数拟合的值与实际值之间的距离：

25            26            27            28 
  17.95294914   32.81291348   21.38677872   26.34105878 
           29            30            31            32

接下来，找到拟合值，这些值应更接近实际的 PM2.5 值以获得最佳模型。使用系数，我们可以计算拟合值：

25         26         27         28         29 
111.047051 115.187087 137.613221 154.658941 154.414781 
        30         31         32         33         34

查找 R-Squared 值。它们应该与你在summary函数输出中找到的多个 R-squared 值相同。R-square 有助于评估模型性能。如果值越接近 1，则模型越好：
```
summary(multiple_PM25_linear_model)$r.squared
```
输出如下：
```
[1] 0.2159579
```
查找 F 统计量值。确保输出与你在summary函数旁边找到的 F 统计量文本的输出相同。这将告诉你你的模型是否比仅使用目标变量的均值拟合得更好。在许多实际应用中，F-Statistic 与 p 值一起使用：
```
   value     numdf     dendf 
3833.506     3.000 41753.000 
```
最后，找到系数 p 值，并确保这些值与你在summary函数下系数部分获得的值相同。它将出现在标题为Pr(>|t|):的列下。如果值小于 0.05，则变量在预测目标变量方面具有统计学意义：
```
  (Intercept)          DEWP          TEMP           Iws 
 0.000000e+00  0.000000e+00  0.000000e+00 4.279601e-224
```
注意

本活动的解决方案可在第 449 页找到。

普通最小二乘法 (OLS)

在第三章，监督学习简介中，我们看到了平方残差和![img/C12624_04_30.png]（也称为误差平方和或残差平方和），这是衡量整体模型拟合程度的指标，其表达式如下：

![img/C12624_04_31.jpg]

其中 T 代表矩阵转置，而 X 的行代表与目标变量的特定值相关的所有输入变量的值![img/C12624_04_32.png]。 A picture containing furniture 的值。

自动生成的描述![img/C12624_04_33.png]，该描述最小化了![img/C12624_04_34.png]，被称为β的OLS 估计量。OLS 算法旨在找到 A picture containing furniture 的全局最小值。

自动生成的描述![img/C12624_04_35.png]，该描述将最小化![img/C12624_04_36.png]。

在上一章中，你也了解到，对于北京 PM2.5 数据集上的multiple_PM25_linear_model，其 R-squared 值相当低，以至于这个模型在实用应用中难以发挥作用。解释这些不良结果的一种方式是说预测变量DEWP和TEMP并没有完全解释 PM2.5 的方差，因此它们未能产生良好的结果。

在我们能够跳入此模型的诊断分析之前，让我们看看是否可以用变量month（读数中的月份）来解释 PM2.5 的一些方差。我们还将使用交互变量（更多内容请参阅改进模型部分）DEWPTEMPmonth，在lm()函数中生成DEWP、TEMP和month的所有可能组合。

使用 month 的原因在 第三章 的 监督学习简介 中的 图 3.3 中得到了证明，其中我们看到了 TEMP、DEWP 和 PRES 的值中的季节性影响（显示出漂亮的正弦波模式）。以下练习的输出显示了为解释 PM2.5 数据集而创建的所有交互项。

注意

表达式如 DEWP:TEMP 表示乘法，并且每个 month 的值在 multiple_PM25_linear_model 中都是一个独立的变量，因为我们运行模型之前将 month 转换成了 factor。

练习 51：在 lm()函数中添加交互项 DEWP:TEMP:month

在这个练习中，我们将添加交互项以提高模型性能。

我们将看到添加额外的交互项如何有助于提高模型性能，从 R-squared 值的角度来看。执行以下步骤以完成练习：

使用以下命令读取北京 PM2.5 数据集：

PM25 <- read.csv("PRSA_data_2010.1.1-2014.12.31.csv")

现在，将 month 对象转换为 factor 变量，如下所示：
```
PM25$month <- as.factor(PM25$month)
```
使用具有 DEWP、TEMP 和 month 的交互项的线性模型。观察 DEWP*TEMP*month 项，这将生成 DEWP、TEMP 和 month 变量的所有组合：
```
multiple_PM25_linear_model <- lm(pm2.5 ~ Iws + DEWP*TEMP*month, data = PM25)
```

打印模型的摘要以查看由于交互项而引起的系数和 r-squared 值的变化：

summary(multiple_PM25_linear_model)

输出如下：

## Call:
## lm(formula = pm2.5 ~ Iws + DEWP * TEMP * month, data = PM25)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -298.41  -42.77   -9.35   30.91  967.39 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)
...
## (Intercept)        2.917e+02  4.338e+00  67.257  < 2e-16 ***
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 70.04 on 41708 degrees of freedom
##   (2067 observations deleted due to missingness)
## Multiple R-squared:  0.4217, Adjusted R-squared:  0.4211 
## F-statistic: 633.7 on 48 and 41708 DF,  p-value: < 2.2e-16

注意 R-squared 值从 0.216 到 0.4217 的两次跳跃。然而，这种跳跃是以模型可解释性为代价的。虽然使用单个变量解释模型的可解释性很简单，但它们的乘积会产生难以描述的效果。

在我们北京 PM2.5 的例子中，更合理地认为 DEWP 和 TEMP 与 year 对象的 month 对象的交互，因为这两个都是随季节变化的环境因素。

然而，我们还想进行一些诊断，以全面了解线性回归模型是如何从头到尾研究的，而不仅仅是查看 R-squared 值。

模型诊断

通常，像线性回归和逻辑回归这样的统计模型在接受最终解决方案之前需要验证许多假设。违反假设的模型会导致错误的预测，结果容易产生误解。

以下代码显示了从 lm() 方法的输出中获取诊断图的方法。该图有四个不同的图，查看残差。让我们了解如何解释每个图。所有这些图都是关于拟合如何符合回归假设的。如果有违反，它将在以下代码的图中清楚地显示：

par(mfrow = c(2,2))
plot(multiple_PM25_linear_model)

输出如下：

图 4.1：北京 PM2.5 数据集上线性模型拟合的诊断图

在接下来的四个部分中，我们将使用来自线性方程和二次方程的随机生成数据探索每个图表，稍后回来解释图 4.1中的四个图表与理想情况相比的表现如何。

注意

在二次方程中，被假定为均值为 0，方差为 2 的正态分布。

在下面的练习中，我们将使用线性方程和二次方程生成图表。稍后，我们将深入理解线性模型应遵循的各种假设，这通过使用通过两个方程生成的随机数据进行的模型拟合来实现。

练习 52：使用线性方程和二次方程生成和拟合模型

在这个练习中，我们将了解线性和多项式函数，以及当我们对两者拟合线性模型时会发生什么。

使用线性方程和多项式方程生成随机数，并在两者上拟合线性模型。观察两个图表之间的差异。

执行以下步骤以生成所需的图表：

首先，使用以下代码定义线性函数：

linear_function <- function(x){return (5+(12*x)-(3*x))}

定义如下命令中的二次函数：

quadratic_function <- function(x){return (5+(12*x)-(3*(x²)))}

现在，生成如图所示的均匀随机数（x）：
```
uniform_random_x <- runif(50, min=0, max=15)
```

使用x生成线性值（y），如图所示：

linear_values_y <- linear_function(uniform_random_x) + rnorm(50,mean = 0, sd =sqrt(2))

使用x生成二次值（y）：

quadratic_values_y <- quadractic_function(uniform_random_x) + rnorm(50,mean = 0, sd =sqrt(2))
df <- data.frame(linear_values_y, quadratic_values_y, uniform_random_x)

使用uniform_random_x对linear_values_y拟合线性模型：

model_df_linear <- lm(linear_values_y ~ uniform_random_x, data = df)

绘制线性关系的诊断图：
```
par(mfrow = c(2,2))
plot(model_df_linear)
```
故事情节如下：

图 4.2：使用线性回归的图表

使用uniform_random_x对quadratic_values_y拟合线性模型：

model_df_quad <- lm(quadratic_values_y ~ uniform_random_x, data = df)

生成非线性关系的诊断：
```
par(mfrow = c(2,2))
plot(model_df_quad)
```
输出如下：

图 4.3：使用二次回归的图表

步骤 7 和步骤 9 中的图表差异显示了线性关系的良好和不良拟合。线性模型不能拟合y和x之间的非线性关系。在下一节中，我们将深入理解步骤 7 和步骤 9 生成的图表的四个部分。

在第三章，监督学习简介，图 3.5中，我们讨论了在构建线性回归模型时需要考虑的各种假设。通过本章前面提到的四个图表，我们将检查是否有任何假设被违反。

残差与拟合图

这种图表位于拟合值和残差（lm()方法的差异）之间。如果预测变量和目标变量之间存在非线性关系，该图表将帮助我们识别。

在以下图中，上部分图显示了点均匀地散布，预测变量和目标变量之间的线性关系被清楚地捕捉到。在下部分图中，未解释的非线性关系被留在残差和曲线中，因此底部图清楚地显示它不是线性回归模型的正确拟合，违反了预测变量和目标变量之间的线性关系：

![Figure 4.4: [Top] Residual versus fitted plot of the linear function. [Bottom] Residual versus fitted plot of the quadratic function]

![img/C12624_04_04.jpg]

图 4.4：[上] 线性函数的残差与拟合图。[下] 二次函数的残差与拟合图

正态 Q-Q 图

Q-Q 图，也称为分位数-分位数图，用于检查数据是否可能来自近似的理论分布；在本例中，为正态分布。Q-Q 图是通过绘制两组分位数（数据中低于一定比例的点）相对比而形成的散点图。如果两组分位数都来自相似的分布，我们应看到点形成一条大致的直线。给定一个数据向量，正态 Q-Q 图按顺序绘制数据与标准正态分布的分位数。

线性回归的第二个假设是所有预测变量都是正态分布的。如果这是真的，残差也将是正态分布的。正态 Q-Q 图是标准化残差与理论分位数之间的图。直观上，我们可以检查残差是否遵循直线，如果是正态分布的，或者是否有任何偏差表明违反了假设。

在以下图中，图的上部分展示了线性函数，它与直线对齐，除了 39 号、30 号和 50 号观测值等少数例外。另一方面，图的下半部分展示了二次函数，它意外地与直线有很好的对齐，不像线性函数，因为图中右上角有一些偏差：

![Figure 4.5: [Top] Normal Q-Q plot of the linear function. [Bottom] Normal Q-Q plot of the quadratic function]

![img/C12624_04_05.jpg]

图 4.5：[上] 线性函数的正态 Q-Q 图。[下] 二次函数的正态 Q-Q 图

标度-位置图

标度-位置图显示残差是否均匀地分布在输入变量（预测变量）的范围内。也可以用此图来检查方差相等的假设（同方差性）。如果我们看到一条水平线，点随机分布，这意味着模型是好的。

该图位于拟合值和标准化残差的平方根之间。在下面的图中，顶部的图显示了线性函数，残差沿着水平线随机分布，而在底部的图中，似乎有一个非随机的模式。因此，方差不相等：

![图 4.6：[上] 线性函数的尺度-位置图。[下] 二次函数的尺度-位置图]

![img/C12624_04_06.jpg]

图 4.6：[上] 线性函数的尺度-位置图。[下] 二次函数的尺度-位置图

残差与杠杆

如果数据中存在任何有影响力的点，残差与杠杆图有助于识别它。通常认为所有异常点都是有影响力的，即它决定了回归线的形状。然而，并非所有异常点都是有影响力的。即使一个点在合理的值范围内（不是异常点），它仍然可能是有影响力的点。

在下一个图中，我们将关注右上角或右下角的远离值。这些区域是观察值相对于回归线可能具有影响力的空间。在图 4.7中，红色虚线外的观察值40和39（高 Cook 距离）。请注意，这些观察值在其他三个图中也持续出现，这给了我们一个强有力的理由，如果我们想看到数据中的线性关系，就应该消除这些点。顶部的图似乎没有红色虚线，这证实了一个良好的拟合：

![图 4.7：[上] 线性函数的残差与杠杆作用图。[下] 二次函数的残差与杠杆作用图]

![img/C12624_04_07.jpg]

图 4.7：[上] 线性函数的残差与杠杆作用图。[下] 二次函数的残差与杠杆作用图]

现在，如果我们重新审视图 4.1，即我们从北京 PM2.5 数据集中获得诊断图；看起来这个模型对于实际应用来说拟合并不是最好的。所有四个图都显示出轻微的违反线性、正态性和同方差性假设。

在下一节中，我们列出了一些改进模型的方法，这些方法可能有助于逐步提高 R 平方值并更好地拟合数据。同样，类似于我们刚才讨论的视觉检查方法，许多统计方法，如用于测试正态性的Kolmogorov-Smirnov 检验，用于测试多重共线性的相关系数，用于测试同方差的Goldfeld-Quandt 检验，都可以使用。

改进模型

到目前为止，我们已经看到了数据中的问题，但你可能会问是否可以修复或改进它。让我们讨论一些改进的方法。在本节中，你将学习一些方法，例如变量转换、处理异常点、添加交互作用以及决定采用非线性模型。

转换预测变量或目标变量

提高模型的最常见方法是使用对数函数变换一个或多个变量（也可能是目标变量）。

对数变换可以纠正偏斜分布。它提供了处理数据中偏斜性的能力，同时一旦建立模型，原始值就可以很容易地计算出来。最流行的对数变换是自然对数。关于对数变换的更详细解释，可以在第六章的特征选择和降维部分的对数变换部分找到。

目标是通过变换将数据中的正态分布引入。因此，任何有助于达到这一目标的函数都是一种好的变换。在取对数之后，平方根变换也被广泛使用。查看变换后变量的分布，看看是否得到了对称分布（钟形）；如果是，那么这种变换将是有用的。

选择非线性模型

可能会遇到一种情况，线性模型并不适合，因为预测变量和目标变量之间存在非线性关系，只有非线性函数才能拟合这样的数据。关于此类模型的更多细节，请参阅本章后面的多项式回归部分。

移除异常值或影响点

正如我们在残差与杠杆率部分的诊断图中讨论的那样，我们可能会发现一个异常值或影响点在获得最佳模型时起到了破坏作用。如果你已经正确地识别了它，尝试通过删除观测值来看是否有所改善。

添加交互效应

有时我们可能会看到数据集中两个或更多预测（独立）变量的值以乘法方式影响因变量。一个带有交互项的线性回归方程可能看起来像这样：

可以尝试更高阶的这种交互（例如，使用三个变量）；然而，这些很难解释，通常会被避免。

分位数回归

当数据出现异常值、高偏度和导致异方差性的条件时，我们采用分位数回归进行建模。此外，分位数回归回答的一个关键问题是，线性回归无法回答的：“对于高 PM2.5 和平均 PM2.5，DEWP、TEMP和Iws对 PM2.5 水平的影响是否不同？”

分位数回归与线性回归非常相似；然而，分位数回归参数估计的是由于输入预测变量单位变化而产生的响应变量某个分位数的改变。为了完全理解这个陈述，让我们使用分位数回归（不使用交互项）来拟合我们的北京数据。

我们需要安装quantreg包来将分位数回归拟合到数据中。该包提供了rq()方法，使用tau参数拟合数据，tau是模型参数，指定用于拟合模型到数据中的分位数值。注意，rq()方法的参数的其他部分看起来与lm()类似。

练习 53：在北京 PM2.5 数据集上拟合分位数回归

在这个练习中，我们将观察在不同分位数（尤其是 25 分位数、50 分位数和 75 分位数）处的分位数回归拟合的差异。我们将使用quantreg中的rq()函数来构建模型。在图 4.8中，我们将比较通过lm()函数获得的系数值与通过rq()函数获得的系数值，以比较两种回归类型。

执行以下步骤来完成练习：

使用以下命令读取北京 PM2.5 数据集：

PM25 <- read.csv("PRSA_data_2010.1.1-2014.12.31.csv")

接下来，下一步是安装所需的包。使用以下命令加载quantreg包：
```
library(quantreg)
```
运行分位数回归 tau 值为 0.25、0.5 和 0.75，分别对应 25 分位数、50 分位数和 75 分位数：
```
quantile_regression_PM25_all <- rq(pm2.5 ~ DEWP+TEMP+Iws, data = PM25, tau = seq(0.25,0.99,by = 0.25))
```

打印分位数回归模型的摘要：

summary(quantile_regression_PM25_all)

输出如下：

## tau: [1] 0.25
## 
## Coefficients:
##             Value     Std. Error t value   Pr(>|t|) 
## (Intercept)  63.62367   0.52894  120.28453   0.00000
## DEWP          2.08932   0.01859  112.39914   0.00000
## TEMP         -1.89485   0.02196  -86.27611   0.00000
## Iws          -0.09590   0.00179  -53.59211   0.00000
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value      Std. Error t value    Pr(>|t|)  
## (Intercept)  117.37344    0.73885  158.85921    0.00000
## DEWP           3.43276    0.02835  121.07849    0.00000
## TEMP          -3.37448    0.03225 -104.65011    0.00000
## Iws           -0.16659    0.00202  -82.56604    0.00000
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value      Std. Error t value    Pr(>|t|)  
## (Intercept)  201.16377    1.31859  152.55927    0.00000
## DEWP           5.12661    0.04901  104.59430    0.00000
## TEMP          -5.62333    0.05567 -101.01841    0.00000
## Iws           -0.25807    0.00510  -50.55327    0.00000

以下表格总结了我们在第三章“监督学习简介”的“回归”部分使用lm()获得的线性回归系数值，以及我们使用rq()在三个分位数获得的值。

根据线性回归模型，大气中 PM2.5 的平均水平随着DEWP增加一个单位而增加4.384196。分位数回归的结果如下表所示，并且它表明DEWP对 PM2.5 的高分位数（观察 75 分位数）有更大的负面影响：

图 4.8：25 分位数、50 分位数、75 分位数分位数回归的系数估计以及北京 PM2.5 估计模型的线性回归系数估计

练习 54：使用更细粒度绘制各种分位数

在这个练习中，我们不会使用第 25、50 和 75 分位数，而是会在rq函数中使用更细粒度的 tau 值。这个图表将有助于可视化系数值根据分位数值的变化。使用 R 中的seq()函数，从 0.05 开始设置分位数值，增量是 0.05。

执行以下步骤来完成练习：

创建一个quantile_regression_PM25_granular变量：

quantile_regression_PM25_granular <- rq(pm2.5 ~ DEWP + TEMP + Iws, data = PM25, tau = seq(0.05,0.95,by = 0.05))

现在，使用summary函数存储先前创建的变量的值：
```
plot_granular <- summary(quantile_regression_PM25_granular)
```
让我们使用以下命令来绘制图表。观察不同 tau 值时，Intercept、DEWP、TEMP和Iws的值如何变化：
```
plot(plot_granular)
```
输出图表如下：

![图 4.9：显示了不同分位数下 DEWP、TEMP 和 Iws 的系数的各种值图片 C12624_04_09.jpg

图 4.9：显示了不同分位数下 DEWP、TEMP 和 Iws 的系数的各种值

在这个练习中，我们通过在rq函数中使用更细粒度的 tau 值来探索变量的粒度。上一张图显示了DEWP、TEMP和Iws的系数的各种值。图中 X 轴表示分位数。单条虚线表示分位数回归的估计，灰色区域是置信区间。中间的灰色线是 OLS 系数估计的表示，双虚线显示 OLS 系数的置信区间。观察发现红色和灰色区域没有重叠，这证明了我们使用分位数回归的正确性。如果两条线重叠，那么使用 OLS 和分位数回归的估计就没有差异。

注意

我们并不声称分位数回归比线性回归给出更好的结果。这个模型的调整 R 平方值仍然很低，但在现实世界中它工作得很好。然而，我们声称分位数回归可以帮助估计 PM2.5 的不同水平，而不仅仅是平均值，这为具有异常值、高偏度和异方差性的数据提供了稳健的解释。

多项式回归

在现实世界的数据中，响应变量和预测变量往往没有线性关系，我们可能需要一个非线性多项式函数来拟合数据。各种散点图样的残差与每个预测变量和残差与拟合值的对比揭示了是否存在违反线性关系的情况，这可能会帮助识别在方程中引入二次或三次项的需要。以下是一个通用的多项式方程：

图片 C12624_04_42.jpg

其中k是多项式的次数。对于k=2，*f(X)*被称为二次，h=4被称为三次。请注意，多项式回归仍然被认为是线性回归，因为它在系数一张脸的插图自动生成描述。

在重新审视北京 PM2.5 的例子之前，让我们通过使用我们在线性回归部分介绍的二项式方程的模拟数据来理解多项式回归是如何工作的。

练习 55：使用 runif()函数执行均匀分布

在这个练习中，我们将使用 R 中的runif()函数生成 50 个来自均匀分布的随机数，并将结果存储在uniform_random_x中。我们已经定义了一个使用之前二次方程生成值的函数。请注意，我们将单独添加一个标志的特写自动生成描述到函数返回的值；一个标志的特写自动生成的描述是使用 R 中的rnorm()函数从正态分布生成的。最终值将被存储在quadratic_values_y中：

执行以下步骤以使用runif()函数进行均匀分布：

首先，定义以下命令中所示的二次方程：

quadratic_function <- function(x){return (5+(12*x)-(3*(x²)))}

现在，生成x的均匀随机数：

uniform_random_x <- runif(50, min=0, max=15)

将误差项添加到二次方程中，该误差项是均值为0、方差为2的正态分布（标准差(sd) = 方差的平方根）：
```
quadratic_values_y <- quadratic_function(uniform_random_x) + rnorm(50,mean = 0, sd =sqrt(2))
```
要将数据存储在数据框中，请使用以下命令：
```
df <- data.frame(quadratic_values_y,uniform_random_x)
```
现在，根据二次方程绘制x和y之间的关系图：
```
library(ggplot2)
ggplot(df, aes(x=uniform_random_x,y=quadratic_values_y))+
  geom_point()
```
输出如下：

![图 4.10：使用函数 runif()进行均匀分布

图 4.10：使用函数 runif()进行均匀分布

下一个图清楚地显示了uniform_random_x和quadratic_values_y之间的关系不是线性的，正如预期的那样。现在，如果我们尝试拟合一个线性模型，我们预计会在诊断图中看到一些问题。

图 4.12中的残差与拟合值图显示了曲率，并且它们没有像之前看到的那样表现出均匀的随机性。此外，正态概率图（NPP）似乎偏离了直线，并在远端百分位数处向下弯曲。这些图表明所使用的模型可能存在问题，并表明可能需要一个更高阶的模型。

![图 4.11：该图显示了均匀生成的随机数(x)与二次方程中 x 的值的非线性关系

图 4.11：该图显示了均匀生成的随机数(x)与二次方程中 x 的值的非线性关系
现在，将线性回归模型拟合到多项式（二次）方程，并显示诊断图：
```
par(mfrow = c(2,2))
plot(lm(quadratic_values_y~uniform_random_x,data=df))
```
输出如下：

![图 4.12：使用 lm()方法对二次方程数据拟合的诊断图

图 4.12：使用 lm()方法对二次方程数据拟合的诊断图

现在，让我们看看多项式回归在北京 PM2.5 数据集上的表现。我们引入了一个额外的二次项 DEWP²，这仅仅是 DEWP 的平方。请参考第三章，监督学习导论中的图 3.5，以证明添加这样一个高阶项的合理性。

在北京 PM2.5 数据集上使用多项式回归（二次和三次项）：

multiple_PM25_poly_model <- lm(pm2.5 ~ DEWP² + TEMP + Iws + DEWP*TEMP*month, data = PM25)

要打印模型摘要，请使用以下命令：

summary(multiple_PM25_poly_model)

输出如下：

## Residuals:
##     Min      1Q  Median      3Q     Max 
## -298.41  -42.77   -9.35   30.91  967.39 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.917e+02  4.338e+00  67.257  < 2e-16 ***
## DEWP               1.190e+01  2.539e-01  46.879  < 2e-16 ***
## TEMP              -9.830e+00  8.806e-01 -11.164  < 2e-16 ***
## Iws               -1.388e-01  7.707e-03 -18.009  < 2e-16 ***
## month2            -2.388e+01  5.011e+00  -4.766 1.89e-06 ***
## month3            -1.228e+02  5.165e+00 -23.780  < 2e-16 ***
## DEWP:TEMP:month9   4.455e-01  6.552e-02   6.800 1.06e-11 ***
## DEWP:TEMP:month10  5.066e-01  5.862e-02   8.642  < 2e-16 ***
## DEWP:TEMP:month11  5.111e-02  5.526e-02   0.925  0.35500    
## DEWP:TEMP:month12  1.492e-01  6.599e-02   2.261  0.02375 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 70.04 on 41708 degrees of freedom
##   (2067 observations deleted due to missingness)
## Multiple R-squared:  0.4217, Adjusted R-squared:  0.4211 
## F-statistic: 633.7 on 48 and 41708 DF,  p-value: < 2.2e-16

观察到尽管增加了二次项，但我们并没有比线性模型获得任何更好的 R-squared 值。在这个时候，我们可能得出结论，PM2.5 预测需要一个更好的自变量，它可以解释其中的方差，从而使 R-squared 值达到更高的水平。诊断图似乎也有类似的解释。

使用以下命令绘制诊断图：

par(mfrow = c(2,2))
plot(multiple_PM25_poly_model)

输出如下：

图 4.13：北京 PM2.5 数据集上多项式回归模型拟合的诊断图

Ridge Regression

正如我们在线性回归中看到的，普通最小二乘法（OLS）通过最小化残差平方和来估计包含家具的图片描述自动生成的值。

由于包含家具的图片描述自动生成是从给定样本中计算出的一个估计值，它不是一个真实的人口参数，我们需要注意估计的某些特性。这两个主要特性是偏差和方差。

如果一个标志的特写描述自动生成是在，然后测试数据集上的平均（或期望）可以分解为三个量，方差、平方偏差以及误差项的方差，如以下方程所示：

对于最佳估计，应同时实现低偏差和低方差的一个合适算法，如 OLS。我们通常称之为偏差-方差权衡。以下图中所示的流行牛眼图有助于理解权衡的各种情况：

图 4.14：解释偏差和方差场景的流行牛眼图

牛眼代表 OLS 试图估计的真实人口值，而对其的射击则是我们四个不同估计器的估计值。这些被广泛归类为以下几类：

低偏差和低方差（最理想）
低偏差和高方差
高偏差和低方差
高偏差和高方差（最不理想）

OLS 方法将所有变量视为同等可能，因此具有低偏差（在训练期间导致欠拟合）和高方差（在测试数据中导致预测误差），如图 4.11所示。这种行为对于获得最佳模型复杂度并不理想。解决这个问题的通用方法是牺牲一些偏差来降低方差。这种方法被称为正则化。因此，岭回归可以被视为线性回归的扩展，其中包含一个额外的正则化项。

多元线性回归的一般形式可以表示如下：

图片 C12624_04_52.jpg

其中，argmin表示使函数达到最小值的βs 的最小值。在此上下文中，它找到使 RSS 最小的βs。βs 受到以下约束：

图片 C12624_04_53.jpg

正则化项 – L2 范数

图片 C12624_04_54.jpg 图片 C12624_04_55.jpg

岭回归的惩罚项在 RSS 增加时增加。以下图表显示了模型复杂度（预测变量数量）与误差之间的关系。它表明，当预测变量的数量增加（模型复杂度增加）时，方差上升，而偏差下降。

OLS 估计在右侧找到一个位置，远离最佳权衡点。这种情况需要引入正则化项，因此岭回归成为合适的模型选择：

![图 4.15：偏差与方差 img/C12624_04_15.jpg

图 4.15：偏差与方差

岭回归的 OLS 损失函数可以表示为以下方程：

图片 C12624_04_56.jpg

最小化一张人脸的插图自动生成的描述函数具有提供岭回归估计的正则化项。这个损失函数的有趣特性是，当 img/C12624_04_58.png 变得更大时，方差减少，偏差增加。

练习 56：北京 PM2.5 数据集上的岭回归

这个练习在 PM2.5 数据集上拟合岭回归。我们将使用glmnet库的交叉验证函数cv.glmnet()，参数alpha = 0和变化的 lambda 值。目标是获得一个最佳的 lambda 值，它将被函数输出的lambda.min属性返回。

让我们执行以下步骤来完成练习：

加载glmnet库并对 PM25 DataFrame 进行预处理：

library(glmnet)
PM25 <- na.omit(PM25)
X <- as.matrix(PM25[,c("DEWP","TEMP","Iws")])
Y <- PM25$pm2.5

现在，让我们使用以下代码来设置seed以获得类似的结果：

set.seed(100)
model_ridge = cv.glmnet(X,Y,alpha = 0,lambda = 10^seq(4,-1,-0.1))

要在交叉验证后找到 lambda 的最佳值，请执行以下命令：
```
optimal_lambda <- model_ridge$lambda.min
```

模型拟合的系数值：

ridge_coefficients <- predict(model_ridge, s = optimal_lambda, type = "coefficients")
ridge_coefficients

输出如下：

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) 160.7120263
## DEWP          4.3462480
## TEMP         -5.0902943
## Iws          -0.2756095

使用predict函数再次传递矩阵 X 到newx参数：

ridge_prediction <- predict(model_ridge, s = optimal_lambda, newx = X)
head(ridge_prediction)

输出如下：

         1
25 111.0399
26 115.1408
27 137.3708
28 154.2625
29 154.0172
30 158.8622

我们看到如何使用glmnet库来拟合北京 PM2.5 数据集。

LASSO 回归

最小绝对收缩和选择算子（LASSO）的结构与岭回归相似，除了惩罚项不同，在 LASSO 回归中是 L1（系数估计值的绝对值之和），而在岭回归中是 L2（系数平方之和）：

img/C12624_04_59.jpg

LASSO 回归将一些系数设为零，从而消除了特定变量的影响。这使得它在拟合数据的同时进行特征选择变得高效。

练习 57：LASSO 回归

在这个练习中，我们将对北京 PM2.5 数据集应用 LASSO 回归。我们将使用相同的 cv.glmnet() 函数来找到最优的 lambda 值。

执行以下步骤以完成练习：

首先，让我们设置 seed 以使用以下命令获得相似的结果：

set.seed(100) #Setting the seed to get similar results.
model_LASSO = cv.glmnet(X,Y,alpha = 1,lambda = 10^seq(4,-1,-0.1))

现在，使用以下命令在交叉验证后找到 lambda 的最佳值：
```
optimal_lambda_LASSO <- model_LASSO$lambda.min
```

执行以下命令从模型拟合中找到系数值：

LASSO_coefficients <- predict(model_LASSO, s = optimal_lambda_LASSO, type = "coefficients")
LASSO_coefficients

输出如下：

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) 160.4765008
## DEWP          4.3324461
## TEMP         -5.0725046
## Iws          -0.2739729

使用以下命令从模型中找到预测值：

LASSO_prediction <- predict(model_LASSO, s = optimal_lambda_LASSO, newx = X)
head(LASSO_prediction)

输出如下：

          1
25 110.9570
26 115.0456
27 137.2040
28 154.0434
29 153.7996
30 158.6282

观察岭回归和 LASSO 回归预测的相似性。北京 PM2.5 数据集在这两种方法中显示没有差异。

弹性网络回归

弹性网络结合了岭回归和 LASSO 回归的惩罚项，以避免在变量选择（系数值趋向于零，从而保持高度相关变量受控）上过度依赖数据。弹性网络最小化以下损失函数：

img/C12624_04_60.jpg

其中参数 α 控制岭回归和 LASSO 之间的正确混合。

总结来说，如果一个模型具有许多预测变量或相关变量，引入正则化项有助于在减少方差的同时优化偏差，从而在模型复杂性和误差之间达到正确的平衡。图 4.16 提供了一个流程图，帮助人们选择在多元回归、岭回归、LASSO 回归和弹性网络回归之间的选择：

Figure 4.16: 选择标准以在多元回归、岭回归、LASSO 回归和弹性网络回归之间进行选择

img/C12624_04_16.jpg

图 4.16：选择标准以在多元回归、岭回归、LASSO 回归和弹性网络回归之间进行选择

练习 58：弹性网络回归

在这个练习中，我们将对北京 PM2.5 数据集执行弹性网络回归。

执行以下步骤以完成练习：

让我们先设置 seed 以使用以下命令获得相似的结果：

set.seed(100)
model_elanet = cv.glmnet(X,Y,alpha = 0.5,lambda = 10^seq(4,-1,-0.1))

现在，使用以下命令在交叉验证后找到 lambda 的最佳值：
```
optimal_lambda_elanet <- model_LASSO$lambda.min
```

接下来，执行以下命令从模型拟合中找到系数值：

elanet_coefficients <- predict(model_elanet, s = optimal_lambda_elanet, type = "coefficients")
elanet_coefficients

输出如下：

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) 160.5950551
## DEWP          4.3393969
## TEMP         -5.0814722
## Iws          -0.2747902

使用以下命令从模型中找到预测值：

elanet_prediction <- predict(model_elanet, s = optimal_lambda_elanet, newx = X)

输出如下：

25 110.9987
26 115.0936
27 137.2880
28 154.1538
29 153.9092
30 158.7461

弹性网络回归（Elastic Net Regression）的预测结果与岭回归（ridge regression）和 LASSO 回归大致相同。在下一节中，我们将比较这三种方法。

系数和残差标准误差的比较

以下表格显示了 DEWP、TEMP 和 Iws 的比较，数值之间没有太大差异，这表明岭回归、LASSO 和弹性网络回归（带有正则化项）并不比多重线性回归方法更好。这也表明 DEWP、TEMP 和 Iws 是低相关或无相关的自变量：

图 4.17：线性、岭回归、LASSO 和弹性网络回归的残差标准误差比较

以下图表显示了使用线性、岭回归、LASSO 和弹性网络回归的截距和 DEWP、TEMP 以及 Iws 变量的系数值的比较：

图 4.18：线性、岭回归、LASSO 和弹性网络回归系数值的比较

练习 59：计算线性、岭回归、LASSO 和弹性网络回归的 RSE

在这个练习中，我们将计算线性、岭回归、LASSO 和弹性网络回归的 RSE。

执行以下步骤以完成练习：

使用以下代码使用 Iws、DEWP 和 TEMP 变量拟合线性模型：

multiple_PM25_linear_model <- lm(pm2.5 ~ Iws + DEWP + TEMP, data = PM25)

现在，使用以下命令来查找线性回归的残差标准误差（RSE）：
```
sqrt(sum(multiple_PM25_linear_model$residuals²)/41753)
```
输出如下：
```
## [1] 81.51
```
同样，我们将找到剩余回归的 RSE。

岭回归的 RSE：

sqrt(sum((Y-ridge_prediction)²)/41753)

输出如下：

## [1] 81.51059

LASSO 回归的 RSE：

sqrt(sum((Y-LASSO_prediction)²)/41753)

输出如下：

## [1] 81.51123

弹性网络回归的 RSE：

sqrt(sum((Y-elanet_prediction)²)/41753)

输出如下：

## [1] 81.51087

这表明所有三个的 RSE 并没有显著差异。

泊松回归

在线性回归中，我们看到了以下形式的方程：

在 Y 是一个计数或比率（Y/t），它具有泊松分布，期望（均值）计数为，等于，即方差。

在逻辑回归的情况下，我们会寻找可以最大化对数似然值的值，以获得系数的最大似然估计器（MLEs）。

没有封闭形式的解，因此最大似然估计的估计将使用迭代算法，如牛顿-拉夫森和迭代加权最小二乘法（IRWLS）。

泊松回归适用于计数依赖变量，必须满足以下指南：

它遵循泊松分布
计数不能为负
数值是整数（没有分数）

注意

这里用于演示泊松回归的数据集来自 A. Colin Cameron 和 Per Johansson，"使用级数展开的计数数据回归：应用"，应用计量经济学杂志，第 12 卷，第 3 期，1997 年，第 203-224 页。

以下表格简要描述了变量：

图 4.19：来自澳大利亚健康调查数据集的变量及其描述

注意

博客 www.econ.uiuc.edu/~econ508/Stata/e-ta16_Stata.html 展示了数据集的使用方法。

练习 60：执行泊松回归

在这个练习中，我们将对数据集进行泊松回归。

执行以下步骤以完成练习：

执行泊松回归，加载foreign库以读取dta数据：
```
library(foreign)
```
使用foreign库中的read.data函数读取澳大利亚健康调查数据集：
```
df_health <- read.dta("health.dta")
```

使用glm()函数拟合广义线性模型，其中泊松回归作为家族参数中的值：

poisson_regression_health <- glm(NONDOCCO ~ ., data = df_health, family=poisson(link=log))

打印模型的摘要：

summary(poisson_regression_health)

输出如下：

## Coefficients:##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.116128   0.137763 -22.620  < 2e-16 ***
## SEX          0.336123   0.069605   4.829 1.37e-06 ***
## AGE          0.782335   0.200369   3.904 9.44e-05 ***
## INCOME      -0.123275   0.107720  -1.144 0.252459    
## LEVYPLUS     0.302185   0.097209   3.109 0.001880 ** 
## FREEPOOR     0.009547   0.210991   0.045 0.963910    
## FREEREPA     0.446621   0.114681   3.894 9.84e-05 ***
## ILLNESS      0.058322   0.021474   2.716 0.006610 ** 
## ACTDAYS      0.098894   0.006095  16.226  < 2e-16 ***
## HSCORE       0.041925   0.011613   3.610 0.000306 ***
## CHCOND1      0.496751   0.086645   5.733 9.86e-09 ***
## CHCOND2      1.029310   0.097262  10.583  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 6127.9  on 5189  degrees of freedom
## Residual deviance: 5052.5  on 5178  degrees of freedom
## AIC: 6254.3
## 
## Number of Fisher Scoring iterations: 7

加载ggplot2库：
```
library(ggplot2)
```

将NONDOCCO的实际值与泊松回归拟合的NONDOCCO值合并：

df_pred_actual <- data.frame(cbind(df_health$NONDOCCO,poisson_regression_health$fitted.values))

命名列名：

colnames(df_pred_actual) <- c("actual_NONDOCCO","predicted_NONDOCCO")

绘制NONDOCCO目标变量的实际值与预测值的对比图：

ggplot(df_pred_actual, aes(x=actual_NONDOCCO, y =predicted_NONDOCCO))+
   geom_point()

输出图表如下：

图 4.20：比较实际值和预测值 NONDOCCO

给定残差偏差统计量的值为 5052.5，自由度为 5178，p 值为零，且5052.5/5178 = 0.975小于 1，因此模型在某种程度上是有效的。我们还可以检查过度分散（数据集中比基于给定统计模型预期的更大变异性）。过度分散是通过将sample_variance除以sample_mean来计算的。让我们检查以下练习。

练习 61：计算过度分散

在这个练习中，我们将对数据集进行过度分散的计算。

执行以下步骤以完成练习：

首先，让我们使用以下命令找到样本均值：
```
s_mean <- mean(df_health$NONDOCCO)
s_mean
```
输出如下：
```
## [1] 0.2146435
```

现在，使用以下命令查找样本方差：

s_variance <- var(df_health$NONDOCCO)
s_variance

输出如下：

## [1] 0.931757

同样，可以使用以下命令计算过度分散：
```
s_variance/s_mean
```
输出如下：
```
## [1] 4.34095
```
因此，即使我们尝试添加预测变量到模型拟合中，过度分散开始下降。在我们的例子中，分散度在合理范围内。
现在，让我们使用以下命令计算分散度：
```
summary.glm(poisson_regression_health)$dispersion
```
输出如下：
```
## [1] 1
```

然而，在分散超过限制的情况下，更高阶的泊松回归是一个合适的解决方案。考虑到本书的范围，我们不会在这里详细探讨此类模型。感兴趣的读者可以阅读有关基线密度泊松（p 阶泊松多项式（PPp）模型）的更多内容。

Cox 比例风险回归模型

Cox 回归模型的基础来自生存分析，一组有助于调查事件发生时间的统计方法。以下是一些例子：

将潜在客户转化为销售的时间
从使用开始到产品故障的时间
从保险单开始到死亡的时间
从诊断到死亡的时间
产品保修索赔的时间
从客户注册开始的时间

所有这些例子都是生存分析的一些用例。在大多数生存分析中，有三种广泛使用的方法用于执行此类时间至事件分析：

用于分析不同组的 Kaplan-Meier 生存曲线
用于比较两个或更多生存曲线的对数秩检验
Cox 比例风险回归用于描述变量对生存的影响

考虑到本章和本书的范围，我们将只关注 Cox 比例风险回归。基本思想是前两种方法仅有助于进行单变量分析，换句话说，您只能理解一个因素对时间至事件的影响，而 Cox 回归有助于评估多个因素对生存时间的影响。此外，Cox 回归对分类和数值因素都同样有效，而前两种方法仅适用于分类因素。

Cox 模型由表示死亡风险的危害函数 h(t) 表示，在医学研究中，它主要被使用，表示在时间 t 死亡的风险：

从这个方程中可以观察到以下几点：

t 表示生存时间。
h(t) 表示由 p 个协变量确定的风险函数。在生存分析中，协变量是用来描述预测变量的术语。
系数表示协变量的影响。
术语被称为时间 t 的基线危害。如果所有系数都为零，就成为危害函数的值。

这个函数看起来与逻辑回归（使用指数项）有些相似，这将在 第五章，分类中详细讨论。我们逻辑上把本书中讨论的所有监督学习算法分为 是/否，1/0），但忽略了事件的时间。

如您从危害函数中观察到的，生存模型包括：

表示事件发生时间的连续变量
一个表示事件是否发生的二元变量

NCCTG 肺癌数据

NCCTG 肺癌数据（来自晚期肺癌患者的生存数据）来自 North Central Cancer Treatment Group。这些数据包括一些元数据，例如收集数据的机构、患者的年龄、性别等。该数据集中的性能评分衡量患者进行日常活动的能力。在任何生存分析数据集中，最重要的变量是关于 事件时间 的知识，例如，死亡时间。

生存分析通常定义为一种检查数据的方法，其中结果变量是感兴趣事件发生的时间。

图 4.21：北中癌症治疗组的变量及其描述

在下一个练习中，我们将学习如何使用survival包中的Surv方法创建生存对象。请注意，在添加生存对象后，数据集的摘要中创建了两个额外的变量SurvObject.time和SurvObject.status，它们存储关于事件发生时间（死亡时间）的信息，这随后成为Cox 比例风险回归模型的因变量。

当患者生存时间周围的指示数量稀少时，观察结果会被截尾。通常，最常见的形式是右侧截尾。假设我们正在跟踪一项为期 20 周的研究。在研究期间未经历感兴趣事件的病人可以被称为右侧截尾。该人的生存时间至少是研究持续时间；在这种情况下，20 周。

练习 62：使用 Cox 回归探索 NCCTG 肺癌数据

在这个练习中，我们将使用 Cox 回归探索 NCCTG 肺癌数据。

执行以下步骤以完成练习：

导入survival库：
```
library(survival)
```
导入肺癌数据：
```
df_lung_cancer <- lung
```

使用head函数打印数据集：

head(df_lung_cancer)

输出如下：

##   inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1    3  306      2  74   1       1       90       100     1175      NA
## 2    3  455      2  68   1       0       90        90     1225      15
## 3    3 1010      1  56   1       0       90        90       NA      15
## 4    5  210      2  57   1       1       90        60     1150      11
## 5    1  883      2  60   1       0      100        90       NA       0
## 6   12 1022      1  74   1       1       50        80      513       0

肺癌数据中，status == 2 代表死亡：

df_lung_cancer$SurvObject <- with(df_lung_cancer, Surv(time, status == 2))

找到 Cox 比例风险回归模型：

cox_regression <- coxph(SurvObject ~ age + sex + ph.karno + wt.loss, data =  df_lung_cancer)
cox_regression

输出如下：

## Call:
## coxph(formula = SurvObject ~ age + sex + ph.karno + wt.loss, 
##     data = df_lung_cancer)
## 
## 
##              coef exp(coef) se(coef)     z      p
## age       0.01514   1.01525  0.00984  1.54 0.1238
## sex      -0.51396   0.59813  0.17441 -2.95 0.0032
## ph.karno -0.01287   0.98721  0.00618 -2.08 0.0374
## wt.loss  -0.00225   0.99776  0.00636 -0.35 0.7239
## 
## Likelihood ratio test=18.8  on 4 df, p=0.000844
## n= 214, number of events= 152 
##    (14 observations deleted due to missingness)

本练习演示了使用生存库的 Cox 比例风险回归模型。

摘要

在本章中，我们在上一章简要介绍线性回归之后，对其进行了更详细的讨论。当然，关于线性回归的讨论导致了一系列诊断，为讨论其他类型的回归算法提供了方向。分位数、多项式、岭回归、LASSO 和弹性网络，所有这些都是从线性回归中派生出来的，区别在于线性回归存在一些局限性，而每个算法都帮助克服了这些局限性。泊松回归和 Cox 比例风险回归模型作为回归算法的特殊情况出现，分别适用于计数和事件发生时间依赖变量，而其他算法适用于任何定量依赖变量。

在下一章中，我们将探讨第二常用的机器学习算法及其相关问题。你还将更详细地了解分类。第五章，分类，与本章类似，旨在详细介绍从决策树到深度神经网络的分类算法。

第五章：分类

学习目标

到本章结束时，你将能够：

定义监督机器学习中的二元分类
使用白盒模型进行二元分类：逻辑回归和决策树
评估监督分类模型的表现
使用黑盒集成模型进行二元分类 - 随机森林和 XGBoost
设计和开发用于分类的深度神经网络
为给定的分类用例选择最佳模型

在本章中，我们将专注于解决监督学习的分类用例。我们将使用一个为分类用例设计的数据库，围绕它构建一个业务问题，并探索一些流行的技术来解决该问题。

简介

让我们快速回顾一下我们在第三章中学到的内容，监督学习简介。正如你现在所知道的，监督学习是机器学习和人工智能的一个分支，它帮助机器在没有明确编程的情况下学习。描述监督学习的一种更简单的方式是开发从标记数据中学习的算法。监督学习的广泛类别包括分类和回归，它们的基本区别在于标签的类型，即连续或分类。处理连续变量的算法被称为回归算法，而处理分类变量的算法被称为分类算法。

在分类算法中，我们的目标、依赖或标准变量是一个分类变量。根据类别的数量，我们可以进一步将它们分为以下几组：

二元分类
多项式分类
多标签分类

在本章中，我们将专注于二元分类。讨论多项式分类和多类分类的细节和实际例子超出了本章的范围；然而，在结束本章之前，我们将列出一些高级主题的额外阅读参考。

二元分类算法是机器学习中最受欢迎的算法类别，在商业、研究和学术界有众多应用。简单的模型可以根据学生的过去表现（通过或失败）来预测他们未来考试通过的可能性，预测是否会下雨，预测客户是否会违约，预测患者是否患有癌症等等，这些都是由分类算法解决的常见用例。

在深入研究算法之前，我们将首先从一个小案例开始，这个小案例将帮助我们通过实际练习解决监督学习分类问题。

开始使用用例

在本章中，我们将参考来自澳大利亚联邦气象局并通过 R 提供的weather数据集。该数据集有两个目标变量，RainTomorrow，一个表示明天是否会下雨的标志，以及RISK_MM，它衡量的是下一天降雨量。

简而言之，我们可以使用这个数据集的RainTomorrow进行分类练习。有关数据集的元数据和附加详细信息可在 www.rdocumentation.org/packages/ra… 上探索。由于数据集可以通过 R 直接使用，我们不需要单独下载它；相反，我们可以直接使用rattle库中的 R 函数将数据加载到系统内存中。

关于用例的一些背景信息

记录了温度、方向、气压、云量、湿度和日照等天气参数，持续一年。下一天的降雨量已经在数据集中作为目标变量RainTomorrow进行工程化。我们可以利用这些数据定义一个机器学习模型，该模型从当天的天气参数中学习，并预测下一天的降雨概率。

雨量预测对许多行业至关重要。火车和长途汽车的长途旅行通常关注天气模式的变化，主要是降雨，以估计到达时间和旅行长度。同样，大多数实体店、小型餐馆和食品摊位等也受到降雨的严重影响。了解明天的天气条件可以帮助企业从多个方面更好地准备，以应对业务损失，在某些情况下，甚至可以最大化业务成果。

为了在问题解决练习中建立良好的直觉，让我们使用数据集构建一个业务问题，并为用例制定问题陈述。由于数据是关于降雨预测的，我们将选择当今超本地食品配送服务面临的流行业务问题。DoorDash、Skip the Dishes、FoodPanda、Swiggy、Foodora 等初创公司为不同国家的客户提供超本地食品配送服务。在大多数国家，一个普遍的趋势是随着雨季的到来，食品配送订单量增加。一般来说，大多数配送公司预计在给定的一天内总配送量会增加 30%-40%。由于配送代理人数有限，雨天订单的增加对配送时间影响巨大。为了保持成本最优，这些公司不可能增加全职代理人数；因此，一个常见的策略是在预计服务需求高的日子里动态雇佣更多代理。为了更好地规划，了解下一天的降雨预测至关重要。

定义问题陈述

在设置好问题背景后，让我们尝试为一家超本地食品配送服务公司定义我们的问题陈述，以预测明天的降雨量。为了保持简单和一致性，让我们使用之前研究的框架，即第二章，数据探索分析来构建问题陈述。这将帮助我们以业务优先的方法提炼出我们想要解决的最终目标，同时将机器学习视角放在首位。

下图展示了之前定义的使用案例的简单视觉框架——情境-复杂性-问题（SCQ）框架：

图 5.1：分类使用案例的 SCQ

我们可以清楚地从 SCQ 中回答问题：我们需要一个预测模型来预测明天的降雨概率，作为解决问题的解决方案。让我们继续下一步——收集数据以构建一个预测模型，这将帮助我们解决业务问题。

数据收集

rattle.data包为我们提供了使用案例的数据，可以使用 R 的内部数据集方法访问。如果你还没有安装这些包，你可以使用install.packages("rattle.data")命令轻松安装它们。

练习 63：探索使用案例的数据

在这个练习中，我们将对为使用案例收集的数据集进行初步探索。我们将探索数据的形状，即行数和列数，并研究每个列中的内容。

为了探索数据的形状（行数 x 列数）和内容，执行以下步骤：

首先，使用以下命令加载rattle包：
```
library(rattle.data)
```
加载我们使用案例的数据，这些数据可以从rattle包中获取：
```
data(weatherAUS)
```
注意

weatherAUS数据集是一个 DataFrame，包含来自 45 个以上澳大利亚气象站的超过 140,000 条每日观测数据。
现在，将天气数据直接加载到名为df的 DataFrame 中：
```
df <- weatherAUS
```
使用str命令探索 DataFrame 的内容：
```
str(df)
```
输出如下：

图 5.2：最终输出

我们有近 150,000 行数据，24 个变量。我们需要删除RISK_MM变量，因为它将是回归使用案例（即预测明天下雨量）的目标变量。因此，我们剩下 22 个独立变量和 1 个因变量RainTomorrow，用于我们的使用案例。我们还可以看到连续变量和分类变量的良好混合。Location、WindDir、RainToday等多个变量是分类的，其余的是连续的。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/2Vwgu8Q。

在下一个练习中，我们将计算每个列中缺失值的总百分比。

练习 64：计算所有列的缺失值百分比

我们在练习 1，探索数据集以用于用例中探索的数据集有相当多的空值。在这个练习中，我们将编写一个脚本来计算每个列中空值的百分比。

我们可以看到几个变量中存在空值。让我们检查df数据集中每个列的空值百分比。

执行以下步骤以计算数据集中每列的空值百分比：

首先，移除名为RISK_MM的列，因为它打算用作回归用途的目标变量。（将其添加到我们的模型会导致数据泄露）：
```
df$RISK_MM <- NULL
```

创建一个temp_df DataFrame 对象并将其值存储在其中：

temp_df<-as.data.frame(
  sort(
  round(
  sapply(df, function(y) sum(length(which(is.na(y)))))/dim(df)[1],2)
  )
)
colnames(temp_df) <- "NullPerc"

现在，使用print函数显示每列的空值百分比，使用以下命令：

print(temp_df)

输出如下：

              NullPerc
Date              0.00
Location          0.00
MinTemp           0.01
MaxTemp           0.01
WindSpeed9am      0.01
Temp9am           0.01
Rainfall          0.02
WindSpeed3pm      0.02
Humidity9am       0.02
Temp3pm           0.02
RainToday         0.02
RainTomorrow      0.02
WindDir3pm        0.03
Humidity3pm       0.03
WindGustDir       0.07
WindGustSpeed     0.07
WindDir9am        0.07
Pressure9am       0.10
Pressure3pm       0.10
Cloud9am          0.38
Cloud3pm          0.41
Evaporation       0.43
Sunshine          0.48

我们可以看到最后四个变量有超过30%的缺失或空值。这是一个相当大的下降。最好从我们的分析中删除这些变量。此外，我们还可以看到一些其他变量大约有1%到2%，在某些情况下，高达10%*的缺失或空值。我们可以使用各种缺失值处理技术来处理这些变量，例如用均值或众数替换它们。在某些重要情况下，我们还可以使用基于聚类的均值和众数替换等额外技术来提高处理效果。此外，在非常关键的场景中，我们可以使用回归模型来估计剩余缺失值的剩余部分，通过定义一个模型，其中所需的缺失值列被视为剩余变量的函数。

注意

您可以在 GitHub 上找到完整的代码：bit.ly/2ViZEp1。

在以下练习中，我们将移除空值。如果没有合适的模型，我们将重新审视数据。

练习 65：从数据集中移除空值

约翰正在处理新创建的数据集，在进行分析时，他发现数据集中存在显著的空值。为了使数据集对进一步分析有用，他必须从其中移除空值。

执行以下步骤以从df数据集中移除空值：

首先，使用以下命令选择最后四列，这些列的空值超过30%：
```
cols_to_drop <-tail(rownames(temp_df),4)
```
使用na.omit命令从 DataFrame 中移除所有包含一个或多个空值列的行，该命令会从 DataFrame 中移除所有空行：
```
df_new<- na.omit(df[,!names(df) %in% cols_to_drop])
```

现在，使用以下print命令打印新格式化的数据：

print("Shape of data after dropping columns:")
print(dim(df_new))

输出如下：

Shape of data after dropping columns:
112925     19

使用以下命令，验证新创建的数据集中是否存在空值：

temp_df<-as.data.frame(sort(round(sapply(df_new, function(y) sum(length(which(is.na(y)))))/dim(df)[1],2)))
colnames(temp_df) <- "NullPerc"

现在，使用以下print命令打印数据集：

print(temp_df)

输出如下：

              NullPerc
Date                 0
Location             0
MinTemp              0
MaxTemp              0
Rainfall             0
WindGustDir          0
WindGustSpeed        0
WindDir9am           0
WindDir3pm           0
WindSpeed9am         0
WindSpeed3pm         0
Humidity9am          0
Humidity3pm          0
Pressure9am          0
Pressure3pm          0
Temp9am              0
Temp3pm              0
RainToday            0
RainTomorrow         0

我们现在可以再次检查，看看新的数据集没有更多的缺失值，数据集的总行数也减少到 112,000 行，这大约是训练数据的*20%损失。我们应该使用替换缺失值（如平均值、众数或中位数）等技术来对抗由于缺失值的省略而导致的高损失。一个经验法则是安全地忽略小于5%*的损失。由于我们有超过 100,000 条记录（对于简单用例来说是一个相当高的记录数），我们正在忽略这个经验法则。

注意

您可以在 GitHub 上找到完整的代码：bit.ly/2Q3HIgT。

此外，我们还可以使用Date列来构建日期和时间相关的特征。以下练习创建了诸如日、月、星期和季度等数值特征作为额外的时相关特征，并删除了原始的Date变量。

我们将使用 R 中的lubridate库来处理日期和时间相关特征。它为我们提供了执行日期和时间操作的极其易于使用的函数。如果您尚未安装此包，请使用install.packages('lubridate')命令安装库。

练习 66：从日期变量中构建基于时间的特征

时间和日期相关的属性不能直接用于监督分类模型。为了从日期和时间相关的变量中提取有意义的属性，通常的做法是从日期中创建月份、年份、星期和季度作为特征。

执行以下步骤以在 R 中处理数据和时间函数：

使用以下命令将lubridate库导入 RStudio：
```
library(lubridate)
```
注意

lubridate库提供了方便的日期和时间相关函数。

使用lubridate函数从Date变量中提取day、month、dayofweek和quarter作为新特征：

df_new$day <- day(df_new$Date)
df_new$month <- month(df_new$Date)
df_new$dayofweek <- wday(df_new$Date)
df_new$quarter <- quarter(df_new$Date)

检查新创建的变量：

str(df_new[,c("day","month","dayofweek","quarter")])

现在我们已经创建了所有日期和时间相关的特征，我们不再需要实际的Date变量。因此，删除之前的Date列：

df_new$Date <- NULL

输出如下：

'data.frame':	112925 obs. of  4 variables:
 $ day      : int  1 2 3 4 5 6 7 8 9 10 ...
 $ month    : num  12 12 12 12 12 12 12 12 12 12 ...
 $ dayofweek: num  2 3 4 5 6 7 1 2 3 4 ...
 $ quarter  : int  4 4 4 4 4 4 4 4 4 4 ...

在这个练习中，我们从日期和时间相关的属性中提取了有意义的特征，并删除了实际的日期相关列。

注意

您可以在 GitHub 上找到完整的代码：bit.ly/2E4hOEU。

接下来，我们需要处理或清理 DataFrame 中的另一个特征：location。

练习 67：探索位置频率

Location变量定义了在指定时间实际捕获天气数据的实际位置。让我们快速检查这个变量中捕获的不同值的数量，看看是否有任何可能重要的有趣模式。

在以下练习中，我们将使用Location变量来定义在指定时间实际捕获天气数据的实际位置。

执行以下步骤：

使用dplyr包中的分组函数计算每个位置的降雨频率：

location_dist <- df_new %>%    group_by(Location) %>%     summarise(Rain  = sum(ifelse(RainTomorrow =="Yes",1,0)), cnt=n()) %>%    mutate(pct = Rain/cnt) %>%    arrange(desc(pct))

检查不同位置的数量以确保正确：

print(paste("#Distinct locations:",dim(location_dist)[1]))

输出如下：

"#Distinct locations: 44"

打印summary以检查执行的聚合：

print(summary(location_dist))

输出如下：

    Location        Rain             cnt            pct         
 Adelaide     : 1   Min.   : 102.0   Min.   : 670   Min.   :0.06687  
 Albury       : 1   1st Qu.: 427.8   1st Qu.:2330   1st Qu.:0.18380  
 AliceSprings : 1   Median : 563.5   Median :2742   Median :0.21833  
 BadgerysCreek: 1   Mean   : 568.6   Mean   :2566   Mean   :0.21896  
 Ballarat     : 1   3rd Qu.: 740.5   3rd Qu.:2884   3rd Qu.:0.26107  
 Bendigo      : 1   Max.   :1031.0   Max.   :3117   Max.   :0.36560  
 (Other)      :38

我们可以看到数据中有 44 个不同的位置。定义每个位置记录数（在之前转换的数据中）的cnt变量，平均有 2,566 条记录。第一四分位数、中位数和第三四分位数的相似数量分布表明，位置在数据中分布均匀。然而，如果我们调查记录降雨的记录百分比（pct），我们会看到一个有趣的趋势。在这里，我们有一些位置大约有 6%的降雨概率，还有一些位置大约有 36%的降雨概率。根据位置的不同，降雨的可能性有很大的差异。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/30aKUMx。

由于我们有大约 44 个不同的位置，直接将此变量作为分类特征使用是有困难的。在 R 中，大多数监督学习算法内部将分类列转换为模型可以解释的数值形式。然而，随着分类变量中类别的数量增加，模型的复杂性也随之增加，但没有额外的价值。为了保持简单，我们可以将Location变量转换为一个具有较少级别的新的变量。我们将选择降雨概率最高的五个和最低的五个位置，并将所有其他位置标记为Others。这将减少变量中不同级别的数量为 10+1，这将更适合模型。

练习 68：创建具有较少级别的新的位置

location变量有太多的不同值（44 个位置），通常机器学习模型在具有高频率不同类别的分类变量上表现不佳。因此，我们需要通过减少变量中不同类别的数量来修剪变量。我们将选择降雨概率最高的五个和最低的五个位置，并将所有其他位置标记为Others。这将减少变量中不同级别的数量为 10+1，这将更适合模型。

执行以下步骤以创建一个具有较少不同级别的位置新变量：

将location变量从因子转换为字符：

location_dist$Location <- as.character(location_dist$Location)

创建一个包含降雨概率最高和最低的五个位置的列表。我们可以通过在按升序排序的 DataFrame 中使用head命令获取前五个位置，以及使用tail命令获取后五个位置来实现这一点：
```
location_list <- c(head(location_dist$Location,5),tail(location_dist$Location,5))
```

打印列表以确认我们已经正确存储了位置：

print("Final list of locations - ")
print(location_list)

输出如下：

[1] "Final list of locations - "
 [1] "Portland"      "Walpole"       "Dartmoor"      "Cairns"       
 [5] "NorfolkIsland" "Moree"         "Mildura"       "AliceSprings" 
 [9] "Uluru"         "Woomera"

将主df_new DataFrame 中的Location变量转换为字符：
```
df_new$Location <- as.character(df_new$Location)
```
减少变量中不同位置的数量。这可以通过将所有不属于location_list列表的位置标记为Others来实现：
```
df_new$new_location <- factor(ifelse(df_new$Location %in% location_list,df_new$Location,"Others"))
```
使用以下命令删除旧的Location变量：
```
df_new$Location <- NULL
```

为了确保第五步正确执行，我们可以创建一个临时 DataFrame，并总结记录频率与我们所创建的新location变量之间的对比：

temp <- df_new %>% mutate(loc = as.character(new_location)) %>%    group_by(as.character(loc)) %>%    summarise(Rain  = sum(ifelse(RainTomorrow =="Yes",1,0)), cnt=n()) %>%    mutate(pct = Rain/cnt) %>%    arrange(desc(pct))

打印临时测试 DataFrame 并观察结果。我们应该只看到 11 个不同的位置值：

print(temp)

输出如下：

# A tibble: 11 x 4
   `as.character(loc)`  Rain   cnt    pct
   <chr>               <dbl> <int>  <dbl>
  Portland             1031  2820 0.366 
  Walpole               864  2502 0.345 
  Dartmoor              770  2294 0.336 
  Cairns                910  2899 0.314 
  NorfolkIsland         883  2864 0.308 
  Others              19380 86944 0.223 
  Moree                 336  2629 0.128 
  Mildura               315  2897 0.109 
  AliceSprings          227  2744 0.0827
  Uluru                 110  1446 0.0761
  Woomera               193  2886 0.0669

我们首先将Location变量从因子转换为字符，以简化字符串操作任务。DataFrame 根据降雨概率的百分比降序排序。head和tail命令用于提取列表中的前五个和后五个位置。然后，这个列表被用作参考检查，以减少新特征中的级别数量。最后，在工程化减少级别的新的特征后，我们进行简单的检查以确保我们的特征已经按照预期的方式工程化。

注意

您可以在 GitHub 上找到完整的代码：bit.ly/30fnR31。

让我们现在进入本章最有趣的主题，并探讨监督学习的分类技术。

监督学习的分类技术

要接近监督分类算法，我们首先需要理解算法的基本功能，以抽象的方式探索一点数学，然后使用 R 中现成的包来开发算法。我们将介绍几个基本算法，例如透明算法如逻辑回归和决策树，然后我们将转向高级建模技术，例如黑盒模型如随机森林、XGBoost 和神经网络。我们计划涵盖的算法列表并不全面，但这五个算法将帮助您对主题有一个广泛的理解。

逻辑回归

逻辑回归是用于二元分类的最受欢迎的透明模型。透明模型定义为我们可以看到用于预测的整个推理过程的模型。对于每个做出的预测，我们可以利用模型的数学方程式来解码预测的原因。也存在一组完全黑盒的分类模型，也就是说，我们根本无法理解模型利用的预测推理。在我们只想关注最终结果的情况下，我们应该选择黑盒模型，因为它们更强大。

尽管名称以回归结尾，但逻辑回归是一种用于预测二元分类结果的技巧。我们需要不同的方法来对分类结果进行建模。这可以通过将结果转换为优势比的对数或事件发生的概率来实现。

让我们将这种方法提炼成更简单的结构。假设一个事件成功的概率为 0.8。那么，同一事件失败的概率将被定义为*(1-0.8) = 0.2*。成功的优势被定义为成功概率与失败概率的比率。

在以下示例中，成功的优势将是*(0.8/0.2) = 4*。这意味着成功的优势是四比一。如果成功的概率是 0.5，即 50-50 的机会，那么成功的优势是 0.5 比 1。逻辑回归模型可以用以下方式数学表示：

其中，是优势比的对数。

进一步解决数学问题，我们可以推导出以下结果：

讨论方程的数学背景和推导超出了本章的范围。为了总结，logit函数，即连接函数，帮助逻辑回归直观地将问题（预测结果）重新表述为优势比的对数。当求解时，它帮助我们预测二元因变量的概率。

逻辑回归是如何工作的？

就像线性回归中，变量的贝塔系数是通过普通最小二乘法（OLS）来估计的，逻辑回归模型利用最大似然估计（MLE）。MLE 函数估计模型参数或贝塔系数的最佳值集，以最大化似然函数，即概率估计，这也可以定义为所选模型与观察数据的一致性。当最佳参数值集被估计出来后，将这些值或贝塔系数按先前定义的方式插入到模型方程中，将有助于估计给定样本的输出概率。类似于 OLS，MLE 也是一个迭代过程。

让我们看看逻辑回归模型在我们数据集上的实际应用。为了开始，我们将只使用模型的一小部分变量。理想情况下，建议根据 EDA 练习开始于最重要的变量，然后逐步添加剩余变量。现在，我们将从一个与最高和最低温度值相关的变量开始，一个与风速相关的变量，下午 3 点的气压和湿度，以及当天的降雨量。

我们将整个数据集分为训练集（70%）和测试集（30%）。在将数据拟合到模型时，我们只将使用训练集，稍后将在训练集以及未见过的测试数据上评估模型的表现。这种方法将帮助我们了解我们的模型是否过拟合，并在未见过的数据上提供更现实的模型性能。

练习 69：构建逻辑回归模型

我们将使用逻辑回归和我们在练习 1-6 中探索的数据集构建一个二元分类模型。我们将数据分为训练集和测试集（分别为 70%和 30%），利用训练数据来拟合模型，并使用测试数据来评估模型在未见数据上的性能。

执行以下步骤来完成练习：

首先，使用以下命令设置 seed 以确保可重复性：
```
set.seed(2019)
```

接下来，为训练数据集（70%）创建一个索引列表：

train_index <- sample(seq_len(nrow(df_new)),floor(0.7 * nrow(df_new)))

现在，使用以下命令将数据分割为测试集和训练集：
```
train <- df_new[train_index,]
test <- df_new[-train_index,]
```
使用 RainTomorrow 作为因变量和几个自变量（我们选择了 MinTemp、Rainfall、WindGustSpeed、WindSpeed3pm、Humidity3pm、Pressure3pm、RainToday、Temp3pm 和 Temp9am）构建逻辑回归模型。我们也可以在 DataFrame 中添加所有可用的自变量：
```
model <- glm(RainTomorrow ~ MinTemp + Rainfall + WindGustSpeed +         WindSpeed3pm +Humidity3pm + Pressure3pm +        RainToday +  Temp3pm + Temp9am,         data=train,        family=binomial(link='logit'))
```

使用 summary 函数打印数据集的摘要：

summary(model)

输出如下：

Call:
glm(formula = RainTomorrow ~ MinTemp + Rainfall + WindGustSpeed + 
    WindSpeed3pm + Humidity3pm + Pressure3pm + RainToday + Temp3pm + 
    Temp9am, family = binomial(link = "logit"), data = train)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9323  -0.5528  -0.3235  -0.1412   3.2047  
Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    6.543e+01  1.876e+00  34.878  < 2e-16 ***
MinTemp        9.369e-05  5.056e-03   0.019    0.985    
Rainfall       7.496e-03  1.404e-03   5.337 9.44e-08 ***
WindGustSpeed  5.817e-02  1.153e-03  50.434  < 2e-16 ***
WindSpeed3pm  -4.331e-02  1.651e-03 -26.234  < 2e-16 ***
Humidity3pm    7.363e-02  9.868e-04  74.614  < 2e-16 ***
Pressure3pm   -7.162e-02  1.821e-03 -39.321  < 2e-16 ***
RainTodayYes   4.243e-01  2.751e-02  15.425  < 2e-16 ***
Temp3pm        3.930e-02  5.171e-03   7.599 2.98e-14 ***
Temp9am       -4.605e-02  6.270e-03  -7.344 2.07e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 83718 on 79046 degrees of freedom
Residual deviance: 56557 on 79037 degrees of freedom
AIC: 56577
Number of Fisher Scoring iterations: 5

set.seed 命令确保用于训练和测试数据集分割的随机选择可以重现。我们将数据分为 70%的训练集和 30%的测试集。设置种子函数确保，对于相同的种子，每次都能得到相同的分割。在 R 中，glm 函数用于构建广义线性模型。在模型中使用 family 参数值设置为 binomial(link ='logit') 定义了逻辑回归。glm 函数还可以用于构建其他几种模型（如 gamma、Poisson 和 binomial）。公式定义了因变量，以及一组自变量。它采用一般形式 Var1 ~ Var2 + Var3 + …，表示 Var1 为因变量或目标变量，其余为自变量。如果我们想使用 DataFrame 中的所有变量作为自变量，我们可以使用 formula = Var1 ~ .，这表示 Var1 是因变量，其余都是自变量。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/2HwwUUX。

逻辑回归结果解释

我们在 第二章，数据探索分析 中曾经简要地了解过逻辑回归，但并未深入探讨模型结果的细节。前一个输出片段中展示的结果将类似于你在线性回归中观察到的结果，但也有一些不同之处。让我们逐部分探索和解释这些结果。

首先，我们有的glm函数计算两种类型的残差，即零偏差和残差偏差。两者的区别在于，一个报告了只使用截距（即没有因变量）时的拟合优度，而另一个报告了使用所有提供的自变量时的拟合优度。零偏差和残差偏差之间的偏差减少有助于我们理解独立变量在定义方差或预测正确性时增加的量化值。偏差残差的分布紧随公式之后报告。

接下来，我们有beta 系数和相关的标准误差，z 值和p 值，即显著性概率。对于每个提供的变量，R 内部计算系数，并连同参数值一起报告额外的测试结果，帮助我们解释这些系数的有效性。系数的绝对值是理解该变量对最终预测能力重要性的简单方法，即该变量在确定预测结果最终结果中的影响程度。我们可以看到所有变量的系数值都很低。

接下来，标准误差帮助我们量化值的稳定性。标准误差的值越低，表明 beta 系数的值越一致或稳定。在我们练习中的所有变量的标准误差都很低。z 值和显著性概率共同帮助我们判断结果是否具有统计学意义，或者只是由于随机机会而看似如此。这个想法遵循我们在第二章，数据探索分析中学到的相同原理，类似于我们在第四章，回归中学到的线性回归参数显著性。

解释显著性的最简单方法就是研究每个自变量旁边的星号，即*。星号的数量由实际概率值定义，如下所示。在我们的练习中，注意MinTemp变量不具有统计学意义，即p 值 > 0.05。其余的都是具有统计意义的变量。

赤池信息量准则（AIC）是 R 报告的另一个指标，用于评估模型的拟合度或模型的质量。这个数字在比较同一用例的不同模型时非常有用。比如说，你使用一组独立变量但相同的因变量拟合了几个模型，AIC 可以通过简单比较所有模型中的值来研究最佳模型。该指标的计算来源于模型预测与实际标签之间的偏差，但考虑了没有增加任何价值的变量。因此，类似于线性回归中的R 平方和调整 R 平方，AIC 帮助我们避免构建复杂的模型。要从候选模型列表中选择最佳模型，我们应该选择 AIC 最低的模型。

在上一个输出的末尾，我们可以看到费舍尔评分算法的结果，这是一种牛顿法求解最大似然问题的衍生方法。我们看到它需要五次迭代来将数据拟合到模型中，但除此之外，这些信息对我们来说价值不大。这仅仅是我们得出模型已经收敛的简单指示。

我们现在已经理解了逻辑回归的工作原理，并解释了模型在 R 中报告的结果。然而，我们仍然需要使用我们的训练集和测试集来评估模型结果，并确保模型在未见过的数据上表现良好。为了研究分类模型的表现，我们需要利用各种指标，例如准确率、精确率和召回率。尽管我们已经在第四章，回归中探讨了它们，但现在让我们更详细地研究它们。

评估分类模型

分类模型需要一系列不同的指标来彻底评估，这与回归模型不同。在这里，我们没有像R 平方这样直观的东西。此外，性能要求完全基于特定的用例。让我们简要地看看我们在第三章，监督学习简介中已经研究过的各种指标，用于分类。

混淆矩阵及其衍生指标

研究分类算法模型性能的第一个基础是从混淆矩阵开始。混淆矩阵是每个类别实际值中预测分布的简单表示：

图 5.3：混淆矩阵

之前的表格是混淆矩阵的简单表示。在这里，我们假设1和0；当结果被正确预测时，我们将其正确预测的1分配为1，依此类推，对于剩余的结果也是如此。

基于混淆矩阵及其定义的值，我们可以进一步定义一些指标，帮助我们更好地理解模型的表现。从现在开始，我们将使用缩写 TP 代表 True Positive（真正），FP 代表 False Positive（假正），TN 代表 True Negative（真负），FN 代表 False Negative（假负）：

整体准确率：整体准确率定义为整个测试样本中正确预测总数与预测总数的比率。因此，这将是真正和真负的总和除以混淆矩阵中的所有指标：

精确度或阳性预测值（PPV）：精确度定义为正确预测的正标签数与所有预测为正标签的总数的比率：

召回率或灵敏度：召回率通过表示正确预测的正标签数与实际正标签总数的比率来衡量你的模型灵敏度：

特异性或真正负率（TNR）：特异性定义为正确预测的负标签数与实际负标签总数的比率：

F1 分数：F1 分数是精确度和召回率的调和平均值。对于大多数情况来说，它是一个比整体准确率更好的指标：

你应该选择哪个指标？

在认真考虑的情况下，另一个重要的方面是我们在评估模型时应考虑哪个指标。没有直接的答案，因为最佳指标组合完全取决于我们处理的分类用例类型。在分类用例中，常见的一种情况是类别不平衡。我们并不总是需要在数据中保持正负标签的平等分布。事实上，在大多数情况下，我们会遇到正类别的数据少于 30% 的情况。在这种情况下，整体准确率并不是一个理想的指标来考虑。

让我们通过一个简单的例子来更好地理解这一点。考虑预测信用卡交易中的欺诈的例子。在现实场景中，对于每 100 笔交易，可能只有一两次欺诈交易。现在，如果我们只用整体准确率作为评估模型的唯一指标，即使我们将所有标签预测为否，即非欺诈，我们会有大约 99% 的准确率，0% 的精确度和 0% 的召回率。99% 的准确率可能看起来是一个很好的模型性能指标；然而，在这种情况下，它并不是一个理想的评估指标。

为了处理这种情况，通常需要额外的业务背景来做出有形的决策，但在大多数情况下（对于此类场景），业务希望召回率更高，同时整体准确率和精度有所妥协。使用高召回率作为模型评估指标的合理性在于，即使交易是真实的，预测交易为欺诈仍然是可以接受的；然而，将欺诈交易预测为真实将是一个错误；业务损失将是巨大的。

通常，模型的评估是根据业务需求结合多种指标进行的。最大的决策者将是精度和召回率之间的权衡。如混淆矩阵所示，每当试图提高精度时，都会损害召回率，反之亦然。

这里有一些我们优先考虑不同指标的业务场景：

预测具有灾难性后果的罕见事件：当预测患者是否有癌症或交易是否为欺诈等情况时，预测一个没有癌症的人患有癌症是可以接受的，但反过来预测会导致生命损失。这些场景需要通过妥协精度和整体准确率来保证高召回率。
预测具有非灾难性后果的罕见事件：当预测客户是否会流失或客户是否会积极回应营销活动时，错误的预测不会危及业务结果，而是会危及活动本身。在这种情况下，根据具体情况，拥有高精度并稍微妥协召回率是有意义的。
预测具有非灾难性后果的常规（非罕见）事件：这将处理大多数分类用例，其中正确预测一个类别的成本几乎等于错误预测该类别的成本。在这种情况下，我们可以使用 F1 分数，它代表精度和召回率之间的调和平均值。使用整体准确率与 F1 分数结合使用将是理想的，因为准确率更容易解释。

评估逻辑回归

现在我们来评估我们之前构建的逻辑回归模型。

练习 70：评估逻辑回归模型

在训练数据集上拟合的机器学习模型不能使用同一数据集进行评估。我们需要利用一个单独的测试数据集，并比较模型在训练集和测试集上的性能。caret包有一些方便的函数来计算之前讨论过的模型评估指标。

执行以下步骤以评估我们在练习 7，构建逻辑回归模型中构建的逻辑回归模型：

计算 DataFrame df_new中RainTomorrow目标变量的记录分布：

print("Distribution of labels in the data-")
print(table(df_new$RainTomorrow)/dim(df_new)[1])

输出如下：

"Distribution of labels in the data-"
       No       Yes 
0.7784459 0.2215541

使用 predict 函数在训练数据上预测 RainTomorrow 目标变量，并将概率值大于（>0.5）的观测值转换为 Yes，否则为 No：

print("Training data results -")
pred_train <-factor(ifelse(predict(model,
                           newdata=train, type="response")> 0.5,"Yes","No"))

为训练数据创建混淆矩阵并打印结果：

train_metrics <- confusionMatrix(pred_train,  
                                     train$RainTomorrow,positive="Yes")
print(train_metrics)

输出如下：

[1] "Training data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  58233  8850
       Yes  3258  8706

               Accuracy : 0.8468          
                 95% CI : (0.8443, 0.8493)
    No Information Rate : 0.7779          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.4998          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.4959          
            Specificity : 0.9470          
         Pos Pred Value : 0.7277          
         Neg Pred Value : 0.8681          
             Prevalence : 0.2221          
         Detection Rate : 0.1101          
   Detection Prevalence : 0.1514          
      Balanced Accuracy : 0.7215          

       'Positive' Class : Yes

与第二步类似，在测试数据上预测结果：

print("Test data results -")
pred_test <-factor(ifelse(predict(model,newdata=test,type = "response") > 0.5,"Yes","No"))

为测试数据预测创建混淆矩阵并打印结果：

test_metrics <- confusionMatrix(pred_test, 
                                test$RainTomorrow,positive="Yes")
print(test_metrics)

输出如下：

[1] "Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  25066  3754
       Yes  1349  3709

               Accuracy : 0.8494          
                 95% CI : (0.8455, 0.8532)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.5042          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.4970          
            Specificity : 0.9489          
         Pos Pred Value : 0.7333          
         Neg Pred Value : 0.8697          
             Prevalence : 0.2203          
         Detection Rate : 0.1095          
   Detection Prevalence : 0.1493          
      Balanced Accuracy : 0.7230          

       'Positive' Class : Yes

我们首先加载必要的 caret 库，该库将提供计算所需指标的功能，如前所述。然后我们使用 R 中的 predict 函数，使用先前拟合的模型在训练数据以及测试数据（分别）上预测结果。predict 函数对于逻辑回归默认返回 link 函数的值。使用 type= 'response' 参数，我们可以覆盖该函数以返回目标的概率。为了简单起见，我们在预测中使用 0.5 作为阈值。因此，任何大于 0.5 的值都会被 confusionMatrix 函数从 caret 库提供给我们一个简单的方式来构建混淆矩阵并计算详尽的指标列表。我们需要将实际标签以及预测标签传递给该函数。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/2Q6mYW0。

目标标签的分布不平衡：77% 为无，23% 为有。在这种情况下，我们不能仅仅依靠整体准确率作为评估模型性能的指标。此外，与之前章节中所示的混淆矩阵相比，如图 3 和 5 的输出所示，混淆矩阵是倒置的。我们有预测作为行，实际值作为列。然而，解释和结果将保持不变。下一组输出报告了感兴趣的指标，以及我们尚未探索的一些其他指标。我们已经涵盖了最重要的指标（敏感性、精确度，即阳性预测值）；然而，建议探索剩余的指标，例如阴性预测值、患病率和检测率。我们可以看到，我们得到的精确度约为 73%，召回率约为 *50%，整体准确率为 85%。在训练和测试数据集上，结果相似；因此，我们可以得出结论，该模型没有过拟合。

注意

总体而言，结果还不错。请不要对低召回率感到惊讶；在我们有不平衡数据集的场景中，用于评估模型性能的指标是业务驱动的。

我们可以得出结论，每次有下雨的可能性时，我们至少可以正确预测一半的时间，并且每次我们预测时，我们都是 73% 正确的。从业务角度来看，如果我们试图考虑是否应该努力提高召回率或精确度，我们需要估计误分类的成本。

在我们的用例中，每当预测到第二天有降雨时，运营管理团队就会准备更多的代理人数以更快地交付。由于没有现成的技术来应对与降雨相关的问题，即使我们只有 50%的召回率，我们也有机会覆盖。在这个问题中，由于错误预测降雨的成本对业务来说会更昂贵，也就是说，如果预测降雨的机会，团队将投资于更多的代理人数以交付，这会带来额外的成本。因此，我们希望有更高的精确度，同时我们可以接受在召回率上的妥协。

注意

理想的情况是具有较高的精确度和召回率。然而，在两者之间总是存在权衡。在大多数现实生活中的机器学习用例中，由业务驱动的决策最终确定优先选择精确度或召回率。

在练习 8，评估逻辑回归模型中开发的先前模型仅使用了df_new数据集中可用的一些变量。让我们构建一个包含数据集中所有可用变量的改进模型，并检查在测试数据集上的性能。

为了迭代改进模型，最佳方式是进行特征选择和超参数调整。特征选择涉及通过各种验证方法从可用列表中选择最佳特征集，并最终确定一个具有最佳性能和最少特征数量的模型。超参数调整涉及构建泛化模型，这些模型不会过拟合，即模型在训练和未见过的测试数据上都能表现良好。这些主题将在第六章，特征选择和降维和第七章，模型改进中详细讨论。现在，本章的范围将仅限于演示模型评估。我们将在后续章节中涉及相同的用例进行超参数调整和特征选择。

练习 71：使用我们用例中所有可用独立变量开发逻辑回归模型

在先前的练习中，我们限制了独立变量的数量，只限于几个。在这个例子中，我们将使用df_new数据集中所有可用的独立变量来创建一个改进的模型。我们再次使用训练数据集来拟合模型，并使用测试数据集来评估模型性能。

执行以下步骤以构建一个包含用例中所有可用独立变量的逻辑回归模型：

使用所有可用独立变量拟合逻辑回归模型：

model <- glm(RainTomorrow~., data=train ,family=binomial(link='logit'))

在训练数据集上预测：

print("Training data results-")
pred_train <-factor(ifelse(predict(model,newdata=train,type = "response") >= 0.5,"Yes","No"))

创建混淆矩阵：

train_metrics <- confusionMatrix(pred_train, train$RainTomorrow,positive="Yes")
print(train_metrics)

输出如下：

"Training data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  58189  8623
       Yes  3302  8933

               Accuracy : 0.8491          
                 95% CI : (0.8466, 0.8516)
    No Information Rate : 0.7779          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.5104          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.5088          
            Specificity : 0.9463          
         Pos Pred Value : 0.7301          
         Neg Pred Value : 0.8709          
             Prevalence : 0.2221          
         Detection Rate : 0.1130          
   Detection Prevalence : 0.1548          
      Balanced Accuracy : 0.7276          

       'Positive' Class : Yes

在测试数据上预测结果：

print("Test data results -")
pred_test <-factor(ifelse(predict(model,newdata=test,type = "response") > 0.5,"Yes","No"))

创建混淆矩阵：

test_metrics <- confusionMatrix(pred_test, test$RainTomorrow,positive="Yes")
print(test_metrics)

输出如下：

"Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  25057  3640
       Yes  1358  3823

               Accuracy : 0.8525          
                 95% CI : (0.8486, 0.8562)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.5176          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.5123          
            Specificity : 0.9486          
         Pos Pred Value : 0.7379          
         Neg Pred Value : 0.8732          
             Prevalence : 0.2203          
         Detection Rate : 0.1128          
   Detection Prevalence : 0.1529          
      Balanced Accuracy : 0.7304          

       'Positive' Class : Yes

我们利用数据集中的所有变量，使用glm函数创建逻辑回归模型。然后我们使用拟合的模型来预测训练集和测试集的结果；类似于之前的练习。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/2HgwjaU。

注意整体准确率、精确率和召回率有所提高（尽管幅度很小）。结果尚可，我们可以通过逻辑回归进一步迭代改进它们。现在，让我们探索一些其他分类技术并研究模型的性能。

注意

在这个练习中，我们没有打印模型的摘要统计信息，类似于第一个模型，只有几个变量。如果打印出来，结果将占用不到两页的章节内容。目前，我们将忽略这一点，因为我们不是在探索 R 报告的模型特征；相反，我们是从训练集和测试集上的准确率、精确率和召回率指标来评估模型的。

获取最佳模型的最理想方式是消除所有统计上不显著的变量，消除多重共线性，处理异常值数据等。鉴于本章的范围，现在我们将忽略所有这些步骤。

活动 8：构建具有额外特征的逻辑回归模型

在练习 8，评估逻辑回归模型中，我们构建了一个具有少量特征的简单模型，然后在练习 9，使用我们用例中所有可用独立变量的逻辑回归模型开发中，我们使用了所有特征。在这个活动中，我们将构建一个逻辑回归模型，该模型具有我们可以使用简单数学变换生成的额外特征。添加额外的数值特征变换，如对数变换、平方和立方幂变换、平方根变换等，是一种良好的实践。

执行以下步骤以开发具有额外特征的逻辑回归模型：

为活动创建df_new数据集的副本到df_copy中，并选择任何三个数值特征（例如，MaxTemp，Rainfall和Humidity3pm）。
为每个选定的特征使用平方和立方幂以及平方根变换来生成新特征。
将df_copy数据集分成 70:30 的训练集和测试集。

使用新的训练数据拟合模型，在测试数据上评估它，最后比较结果。

输出如下：

"Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  25057  3640
       Yes  1358  3823

               Accuracy : 0.8525          
                 95% CI : (0.8486, 0.8562)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.5176          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.5123          
            Specificity : 0.9486          
         Pos Pred Value : 0.7379          
         Neg Pred Value : 0.8732          
             Prevalence : 0.2203          
         Detection Rate : 0.1128          
   Detection Prevalence : 0.1529          
      Balanced Accuracy : 0.7304          

       'Positive' Class : Yes

注意

你可以在第 451 页找到这个活动的解决方案。

决策树

与逻辑回归一样，还有一种流行的分类技术因其简单性和白盒特性而非常受欢迎。决策树是一个以树（倒置树）的形式表示的简单流程图。它从一个根节点开始，分支到几个节点，可以根据决策进行遍历，并以叶节点结束，其中确定 最终结果。决策树可以用于回归，也可以用于分类用例。机器学习中实现了多种决策树的变体。这里列出了几个流行的选择：

迭代二分器 3（ID3）
ID3 的继承者（C4.5）
分类和回归树（CART）
CHi-squared 自动交互检测器（CHAID）
条件推理树（C 树）

前面的列表并不详尽。还有其他替代方案，它们在如何处理分类和数值变量、选择树中根节点和连续节点的方法、分支每个决策节点的规则等方面都有细微的差别。在本章中，我们将限制我们的探索范围到 CART 决策树，这是最广泛使用的。R 提供了一些包含 CART 算法实现的包。在我们深入研究实现之前，让我们在以下几节中探讨决策树的一些重要方面。

决策树是如何工作的？

决策树的每个变体都有略微不同的方法。总的来说，如果我们尝试简化通用决策树的伪代码，可以总结如下：

选择根节点（节点对应一个变量）。
将数据划分为组。
对于上一步中的每个组：

创建一个决策节点或叶节点（基于分割标准）。

重复执行，直到节点大小 <= 阈值或特征 = 空集。

不同形式的树实现之间的差异包括处理分类和数值变量的方式、选择树中根节点和连续节点的方法、分支每个决策节点的规则等。

下面的视觉图是一个示例决策树。根节点和决策节点是我们提供给算法的独立变量。叶节点表示最终结果，而根节点和中间决策节点有助于遍历数据到叶节点。决策树的简单性使其非常有效且易于解释。这有助于轻松识别预测任务中的规则。通常，许多研究和商业倡议利用决策树来设计一套简单的分类系统规则：

图 5.4：示例决策树

在一般意义上，给定一组依赖变量和多个独立变量，决策树算法计算一个指标，该指标表示依赖目标变量与所有独立变量之间的拟合优度。对于分类用例，熵和信息增益是 CART 决策树中常用的指标。对于该指标的最好拟合变量被选为根节点，下一个最佳拟合变量被用作决策节点，按拟合优度降序排列。节点根据定义的阈值终止为叶节点。树继续增长，直到耗尽决策节点的变量数量或达到预定义的节点数量阈值。

为了提高树性能并减少过拟合，一些策略，如限制树的深度或宽度，或为叶节点或决策节点添加额外的规则，有助于对预测进行泛化。

让我们使用 R 中的 CART 决策树实现相同的用例。CART 模型通过 R 中的rpart包提供。该算法由 Leo Breiman、Jerome Friedman、Richard Olshen 和 Charles Stone 在 1984 年开发，并在业界得到广泛应用。

练习 72：在 R 中创建决策树模型

在这个练习中，我们将使用与我们在 练习 9 中使用的相同数据和用例，在 R 中创建决策树模型，该用例是 使用我们用例中所有可用独立变量开发逻辑回归模型。我们将尝试研究决策树模型与逻辑回归模型在性能上是否存在任何差异。

在 R 中创建决策树模型的以下步骤：

使用以下命令导入rpart和rpart.plot包：
```
library(rpart)
library(rpart.plot)
```

使用所有变量构建 CART 模型：

tree_model <- rpart(RainTomorrow~.,data=train)

绘制成本参数：
```
plotcp(tree_model)
```
输出如下：

图 5.5：决策树模型
使用以下命令绘制树：
```
rpart.plot(tree_model,uniform=TRUE, main="Predicting RainFall")
```
输出如下：

图 5.6：预测降雨

在训练数据上做出预测：

print("Training data results -")
pred_train <- predict(tree_model,newdata = train,type = "class")
confusionMatrix(pred_train, train$RainTomorrow,positive="Yes")

输出如下：

"Training data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  59667 11215
       Yes  1824  6341

               Accuracy : 0.835           
                 95% CI : (0.8324, 0.8376)
    No Information Rate : 0.7779          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.4098          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.36119         
            Specificity : 0.97034         
         Pos Pred Value : 0.77661         
         Neg Pred Value : 0.84178         
             Prevalence : 0.22210         
         Detection Rate : 0.08022         
   Detection Prevalence : 0.10329         
      Balanced Accuracy : 0.66576         

       'Positive' Class : Yes

在测试数据上做出预测：

print("Test data results -")
pred_test <- predict(tree_model,newdata = test,type = "class")
confusionMatrix(pred_test, test$RainTomorrow,positive="Yes")

输出如下：

[1] "Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  25634  4787
       Yes   781  2676

               Accuracy : 0.8356          
                 95% CI : (0.8317, 0.8396)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.4075          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.35857         
            Specificity : 0.97043         
         Pos Pred Value : 0.77408         
         Neg Pred Value : 0.84264         
             Prevalence : 0.22029         
         Detection Rate : 0.07899         
   Detection Prevalence : 0.10204         
      Balanced Accuracy : 0.66450         

       'Positive' Class : Yes

rpart库为我们提供了决策树的 CART 实现。还有其他库帮助我们可视化 R 中的决策树。我们在这里使用了rpart.plot。如果包尚未安装，请使用install.packages命令安装。我们使用rpart函数创建树模型，并使用所有可用的独立变量。然后我们使用plotcp函数可视化复杂度参数对应的不同迭代中的验证误差。我们还使用plot.rpart函数绘制决策树。

最后，我们在训练集和测试集上做出预测，并使用confusionMatrix函数分别对训练集和测试集计算感兴趣的指标，构建混淆矩阵。

注意

您可以在 GitHub 上找到完整的代码：bit.ly/2WECLgZ。

R 中实现的 CART 决策树已经内置了一些优化。默认情况下，该函数设置了许多参数以获得最佳结果。在决策树中，有几个参数我们可以手动设置，以根据我们的需求调整性能。然而，R 实现通过默认值设置大量参数，并取得了相对较好的效果。这些额外的设置可以通过control参数添加到rpart树中。

我们可以将以下参数添加到树模型中：

control = rpart.control(
    minsplit = 20, 
    minbucket = round(minsplit/3), 
    cp = 0.01, 
    maxcompete = 4, 
    maxsurrogate = 5, 
    usesurrogate = 2, xval = 10, 
    surrogatestyle = 0, 
    maxdepth = 30
)

一个有趣的参数是0.01。我们可以将其进一步更改为一个更小的数字，这将使树生长得更深，变得更加复杂。plotcp函数可视化了不同cp值（即复杂度参数）的相对验证误差。在图 5.4中，最理想的cp值是位于虚线以下的最左侧值。在这种情况下（如图所示），最佳值是 0.017。由于这个值与默认值相差不大，我们没有进一步更改它。

下一个图表（图 5.5）帮助我们可视化算法构建的实际决策树。我们可以看到正在使用可用数据构建的简单规则集。正如您所看到的，只有两个独立变量，即Humidity3pm和WindGustSpeed，被选入树中。如果我们将复杂度参数从0.01改为0.001，我们可以看到一个更深层次的树（这可能导致模型过拟合）将被构建。最后，我们可以看到混淆矩阵（步骤 6）的结果，以及训练集和测试集的感兴趣的其他指标。

我们可以看到，训练集和测试集的结果相似。因此，我们可以得出结论，该模型没有过拟合。然而，准确率（83%）和召回率（35%）有显著下降，而精确率则增加到略高的值（77%）。

我们现在已经与几种白盒建模技术合作过。鉴于白盒模型的简单性和易于解释，它们是商业中分类用例中最受欢迎的技术，在这些用例中，推理和驱动分析至关重要。然而，在某些情况下，企业可能对模型的净结果更感兴趣，而不是对结果的整个解释。在这种情况下，最终模型的表现更为重要。在我们的用例中，我们希望实现高精度。让我们探索一些在模型性能上优于（在大多数情况下）白盒模型，并且可以用更少的努力和更多的训练数据实现的黑盒模型。

活动 9：创建具有额外控制参数的决策树模型

我们在练习 10，在 R 中创建决策树模型中创建的决策树模型使用了树的默认控制参数。在这个活动中，我们将覆盖一些控制参数，并研究其对整体树拟合过程的影响。

执行以下步骤以创建具有额外控制参数的决策树模型：

加载rpart库。
为决策树创建控制对象，并使用新值：minsplit =15和cp = 0.00。
使用训练数据拟合树模型，并将控制对象传递给rpart函数。
绘制复杂性参数图，以查看树在不同CP值下的表现。
使用拟合的模型对训练数据进行预测，并创建混淆矩阵。

使用拟合的模型对测试数据进行预测，并创建混淆矩阵。

输出如下：

"Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  25068  3926
       Yes  1347  3537

               Accuracy : 0.8444          
                 95% CI : (0.8404, 0.8482)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.4828          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.4739          
            Specificity : 0.9490          
         Pos Pred Value : 0.7242          
         Neg Pred Value : 0.8646          
             Prevalence : 0.2203          
         Detection Rate : 0.1044          
   Detection Prevalence : 0.1442          
      Balanced Accuracy : 0.7115          

       'Positive' Class : Yes

注意

你可以在第 454 页找到这个活动的解决方案。

集成建模

当需要在大规模训练样本中提高性能时，集成建模是分类和回归建模技术中最常用的方法之一。简单来说，集成建模可以通过分解其名称为个别术语来定义：集成和建模。我们已经在本书中研究了建模；简单来说，集成就是一个组。因此，为同一任务构建多个模型而不是一个模型，然后将结果通过任何方式（如平均或投票）结合成一个单一结果的过程，称为集成建模。

我们可以构建任何模型的集成，例如线性模型或树模型，实际上甚至可以构建集成模型的集成。然而，最流行的方法是使用树模型作为集成的基座。集成模型主要有两种类型：

Bagging: 在这里，每个模型都是并行构建的，并在每个模型内部引入了一些随机化，所有模型的结果通过简单的投票机制进行组合。比如说我们构建了 100 个树模型，其中 60 个模型预测结果为是，40 个模型预测为否。最终结果将是是。
Boosting: 在这里，模型是按顺序构建的，第一个模型的结果被用来调整下一个模型。每个模型迭代地从先前模型犯的错误中学习，并试图在连续迭代中改进。结果是所有单个结果的加权平均值。

在 Bagging 和 Boosting 中都有几种实现方式。Bagging本身是 R 中可用的一种集成模型。迄今为止，最流行的 Bagging 技术是随机森林。与随机森林类似，另一种 Bagging 技术是额外树。同样，一些 Boosting 技术的例子包括 AdaBoost、随机梯度提升、BrownBoost 等。然而，最流行的 Boosting 技术是XGBoost，它源自名称EXtreme Gradient Boosting。在大多数情况下，对于分类以及回归用例，数据科学家更喜欢使用随机森林或 XGBoost 模型。Kaggle（一个在线数据科学社区）的一项最近调查显示，在大多数机器学习竞赛中最常用的技术总是随机森林和 XGBoost。在本章中，我们将更深入地研究这两种模型。

随机森林

随机森林是机器学习中使用的最流行的 Bagging 技术。它是由 CART 的作者 Leo Brieman 开发的。这种简单技术非常有效，以至于在给定的监督用例中，数据科学家几乎总是首先选择算法。随机森林是分类和回归用例的良好选择。它是一种高度有效的方法，可以以最小的努力减少过拟合。让我们更深入地了解随机森林是如何工作的。

如我们所知，随机森林是一种集成建模技术，其中我们构建了多个模型，并使用简单的投票技术结合它们的结果。在随机森林中，我们使用决策树作为基础模型。算法的内部工作原理可以从其名称本身推测出来，即随机（因为它在构建的每个模型中引入了一层随机化）和森林（因为我们要构建多个树模型）。在我们深入了解算法的实际工作原理之前，我们首先需要了解其前身Bagging的故事，并研究为什么我们需要集成。

为什么使用集成模型？

在你的脑海中浮现的第一个问题可能就是，为什么我们一开始就需要为同一个任务构建多个模型？这是必要的吗？嗯，是的！当我们构建集成模型时，我们并不是多次构建完全相同的模型；相反，我们构建的每个模型都会以某种方式与其他模型不同。这个背后的直觉可以通过我们日常生活中的一个简单例子来理解。它是基于这样一个原则：将几个弱学习器结合起来可以构建一个更强、更健壮的模型。

让我们用一个简单的例子来理解这个想法。假设你到达一个新城市，想知道第二天这个城市下雨的概率。假设技术不是一个可用的选项，你能找到的最简单的方法就是询问在这个地方居住了一段时间的邻居。也许答案并不总是正确的；如果有人说第二天有很大的下雨概率，这并不一定意味着一定会下雨。因此，为了做出改进的猜测，你会询问几个邻居。现在，如果你询问的 10 个人中有 7 个人提到第二天有很高的下雨概率，那么第二天几乎肯定会下雨。这个方法之所以有效，是因为你接触到的每个人都会对降雨模式有一定的了解，而且每个人的了解也会有所不同。尽管这些差异并不大，但当这些了解汇总成一个集体答案时，就会产生更好的答案。

Bagging – 随机森林的前身

集成建模遵循相同的原理。在这里，每个模型中都会引入一定程度的随机性。Bagging 算法为每个模型在训练数据上引入这种随机性。Bagging 这个名字来源于自助聚合；这是一个过程，我们用替换数据从可用的数据中抽取三分之二的数据用于训练，其余的用于测试和验证。在这里，每个模型，即决策树模型，在略微不同的数据集上训练，因此对于相同的测试样本可能会有略微不同的结果。Bagging 在某种程度上模仿了我们讨论的现实世界例子，因此将几个弱学习器（决策树模型）组合成一个强学习器。

随机森林是如何工作的？

随机森林基本上是 bagging 的继承者。在这里，除了训练数据的随机性外，随机森林还通过特征集添加了一个额外的随机层。因此，每个决策树不仅具有自助聚合，即三分之二的训练数据有放回地替换，而且还从可用列表中随机选择特征的一个子集。因此，集成中的每个单个决策树都有略微不同的训练数据集和略微不同的特征集进行训练。这额外的随机层在泛化模型方面非常有效，并减少了方差。

练习 73：在 R 中构建随机森林模型

在这个练习中，我们将在练习 8、9 和 10 中使用的相同数据集上构建一个随机森林模型。我们将利用集成建模，并测试整体模型性能是否优于决策树和逻辑回归。

注意

要开始，我们可以快速使用之前使用的相同数据集构建一个随机森林模型。R 中的 randomForest 包提供了模型的实现，以及一些额外的函数来优化模型。

让我们看看一个基本的随机森林模型。执行以下步骤：

首先，使用以下命令导入 randomForest 库：
```
library(randomForest)
```

使用所有可用的独立特征构建一个随机森林模型：

rf_model <- randomForest(RainTomorrow ~ . , data = train, ntree = 100,                                             importance = TRUE, 
                                            maxnodes=60)

在训练数据上评估：

print("Training data results -")
pred_train <- predict(rf_model,newdata = train,type = "class")
confusionMatrix(pred_train, train$RainTomorrow,positive="Yes")

在测试数据上评估：

print("Test data results -")
pred_test <- predict(rf_model,newdata = test,type = "class")
confusionMatrix(pred_test, test$RainTomorrow,positive="Yes")

绘制特征重要性：

varImpPlot(rf_model)

输出如下：

[1] "Training data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  59630 10133
       Yes  1861  7423

               Accuracy : 0.8483          
                 95% CI : (0.8457, 0.8508)
    No Information Rate : 0.7779          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.472           
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.42282         
            Specificity : 0.96974         
         Pos Pred Value : 0.79955         
         Neg Pred Value : 0.85475         
             Prevalence : 0.22210         
         Detection Rate : 0.09391         
   Detection Prevalence : 0.11745         
      Balanced Accuracy : 0.69628         

       'Positive' Class : Yes             

[1] "Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction    No   Yes
       No  25602  4369
       Yes   813  3094

               Accuracy : 0.847           
                 95% CI : (0.8432, 0.8509)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.4629          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.41458         
            Specificity : 0.96922         
         Pos Pred Value : 0.79191         
         Neg Pred Value : 0.85423         
             Prevalence : 0.22029         
         Detection Rate : 0.09133         
   Detection Prevalence : 0.11533         
      Balanced Accuracy : 0.69190         

       'Positive' Class : Yes

图 5.7：随机森林模型

注意

你可以在 GitHub 上找到完整的代码：bit.ly/2Q2xKwd。

活动 10：构建具有更多树的随机森林模型

在 练习 11，在 R 中构建随机森林模型 中，我们创建了一个只有 100 树的随机森林模型；我们可以构建一个具有更多树的更健壮的模型。在这个活动中，我们将创建一个拥有 500 树的随机森林模型，并研究模型只有 100 树的影响。一般来说，我们预计模型的性能会提高（至少随着树的数量增加而略有提高）。这伴随着模型收敛所需更高的计算时间。

执行以下步骤以构建一个拥有 500 树的随机森林模型：

开发一个具有更多树的随机森林模型；比如说，500 树。鼓励读者尝试更高的数字，如 1,000、2,000 等，并研究每个版本的增量改进。
利用拟合的模型在训练和测试数据上预测估计值，并研究与拥有 100 树的模型相比是否有任何改进。

注意

你可以在第 457 页找到这个活动的解决方案。

XGBoost

XGBoost 是近年来最受欢迎的增强技术。尽管有大型公司开发了各种新版本，但 XGBoost 仍然稳居王座。让我们简要回顾一下增强技术的历史。

提升过程是如何工作的？

提升法与装袋法在核心原则上有区别；学习过程实际上是顺序的。在集成中构建的每个模型理想上都是前一个模型的改进版本。用简单的话来说，想象你正在玩一个游戏，你必须记住你只被展示了一次、持续 30 秒的所有放在桌子上的物品。游戏的主持人将大约 50-100 种不同的物品放在桌子上，如蝙蝠、球、钟、骰子、硬币等等，然后用一块大布覆盖它们。当游戏开始时，他从桌子上取下布，给你 30 秒的时间看它们，然后再次拉上布。你现在必须回忆起你能记住的所有物品。能回忆起最多物品的参与者赢得了游戏。

在这个游戏中，让我们增加一个新维度。假设你是一个团队，玩家们一个接一个地轮流宣布他们能回忆起的所有物品，而其他人则聆听。假设有 10 名参与者；每个参与者走上前来，大声宣布他们从桌子上能回忆起的物品。当第二个玩家走上前来时，他们已经听到了第一个玩家宣布的所有物品。他们可能会提到一些第二个玩家可能没有回忆起的物品。为了改进第一个玩家的表现，第二个玩家从第一个玩家那里学习了一些新的物品，将它们添加到自己的列表中，然后大声宣布。当最后一名玩家走上前来时，他们已经学习了其他玩家回忆起的几个物品，而这些物品他们自己未能回忆起来。

将这些放在一起，那个玩家创建了最详尽的列表，并赢得了比赛。每个玩家依次宣布列表的事实有助于下一个玩家从他们的错误中学习并改进。

提升法以相同的方式工作。每个按顺序训练的模型都获得了额外的知识，使得第一个模型的错误在第二个模型中得到了更好的学习。比如说，第一个模型学会了在特定独立变量的大多数情况下进行良好的分类；然而，它未能正确预测一个特定的类别。下一个模型被赋予了不同的训练样本，使得模型在先前模型失败的类别中学习得更好。一个简单的例子是基于感兴趣变量或类别的过采样。提升法有效地减少了偏差，因此提高了模型的表现。

有哪些流行的提升技术？

之前介绍的提升技术并不太受欢迎，因为它们很容易过拟合，并且通常需要相对较多的努力来调整以达到出色的性能。AdaBoost、BrownBoost、梯度提升和随机梯度提升都是长期流行的提升技术。然而，在 2014 年，当 T Chen 和其他人引入 XGBoost（极端梯度提升）时，它为提升性能带来了新的高度。

XGBoost 是如何工作的？

XGBoost 原生引入了正则化，这有助于模型对抗过拟合，从而实现了高性能。与其他当时可用的提升技术相比，XGBoost 显著减少了过拟合问题，并且所需的工作量最少。在 R 或任何其他语言的当前模型实现中，XGBoost 几乎总是使用默认参数设置表现出色。（尽管，这并不总是正确的；在许多情况下，随机森林的表现优于 XGBoost）。XGBoost 一直是数据科学黑客马拉松和企业项目中使用的最流行的算法之一。

简而言之，XGBoost 在目标函数中引入了正则化，当模型在训练迭代中变得更加复杂时，会对模型进行惩罚。讨论 XGBoost 中涉及的数学结构的深度超出了本章的范围。您可以参考 T Chen 的论文以获取更多信息（arxiv.org/abs/1603.02… GBM 和 XGBoost 之间的数学差异：towardsdatascience.com/boosting-al…

在 R 中实现 XGBoost

我们可以利用 XGBoost 包，它提供了算法的整洁实现。在开始之前，我们需要注意一些实现方法上的差异。与其他 R 中算法的实现不同，XGBoost 不处理分类数据（其他实现会将其内部转换为数值数据）。XGBoost 在 R 中的内部工作方式不处理将分类列自动转换为数值列。因此，我们手动将分类列转换为数值或独热编码形式。

独热编码形式基本上是将单个分类列表示为二进制编码形式。比如说我们有一个包含是/否/可能等值的分类列；然后，我们将这个单一变量转换为，其中我们为分类变量的每个值都有一个单独的变量，表示其值为0或1。因此，是、否和可能的值将根据原始值取0和1。

独热编码在以下表格中演示：

图 5.8：独热编码

让我们将数据转换为所需的形式，并在数据集上构建一个 XGBoost 模型。

练习 74：在 R 中构建 XGBoost 模型

正如我们在练习 11，在 R 中构建随机森林模型中做的那样，我们将尝试通过为与练习 11，在 R 中构建随机森林模型相同的用例和数据集构建 XGBoost 模型来提高分类模型的性能。

执行以下步骤以在 R 中构建 XGBoost 模型。

为目标、分类和数值变量创建列表占位符：

target<- "RainTomorrow"
categorical_columns <- c("RainToday","WindGustDir","WindDir9am",
"WindDir3pm", "new_location")
numeric_columns <- setdiff(colnames(train),c(categorical_columns,target))

将分类因子变量转换为字符。这将在将它们转换为独热编码形式时很有用：
```
df_new <- df_new %>% mutate_if(is.factor, as.character)
```

使用caret包中的dummyVars函数将分类变量转换为独热编码形式：

dummies <- dummyVars(~ RainToday + WindGustDir + WindDir9am + 
                     WindDir3pm + new_location, data = df_new)
df_all_ohe <- as.data.frame(predict(dummies, newdata = df_new))

将第三步中的数值变量和独热编码变量合并到一个名为df_final的单个 DataFrame 中：
```
df_final <- cbind(df_new[,numeric_columns],df_all_ohe)
```
将目标变量转换为数值形式，因为 R 中的 XGBoost 实现不接受因子或字符形式：
```
y <- ifelse(df_new[,target] == "Yes",1,0)
```

将df_final数据集分为训练集（70%）和测试集（30%）：

set.seed(2019)
train_index <- sample(seq_len(nrow(df_final)),floor(0.7 * nrow(df_final)))
xgb.train <- df_final[train_index,]
y_train<- y[train_index]
xgb.test <- df_final[-train_index,]
y_test <- y[-train_index]

使用xgboost函数构建 XGBoost 模型。传递训练数据和y_train目标变量，并定义eta = 0.01、max_depth = 6、nrounds = 200和colsample_bytree = 1超参数，定义评估指标为logloss，目标函数为binary:logistic，因为我们处理的是二元分类：

xgb <- xgboost(data = data.matrix(xgb.train), 
               label = y_train, 
               eta = 0.01,
               max_depth = 6, 
               nround=200, 
               subsample = 1,
               colsample_bytree = 1,
               seed = 1,
               eval_metric = "logloss",
               objective = "binary:logistic",
               nthread = 4
)

使用训练数据集上的拟合模型进行预测，并创建混淆矩阵以评估模型在训练数据上的性能：

print("Training data results -")
pred_train <- factor(ifelse(predict(xgb,data.matrix(xgb.train),type="class")>0.5,1,0))
confusionMatrix(pred_train,factor(y_train),positive='1')

输出如下：

"Training data results -"
Confusion Matrix and Statistics
          Reference
Prediction     0     1
         0 58967  8886
         1  2524  8670

               Accuracy : 0.8557          
                 95% CI : (0.8532, 0.8581)
    No Information Rate : 0.7779          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.5201          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.4938          
            Specificity : 0.9590          
         Pos Pred Value : 0.7745          
         Neg Pred Value : 0.8690          
             Prevalence : 0.2221          
         Detection Rate : 0.1097          
   Detection Prevalence : 0.1416          
      Balanced Accuracy : 0.7264          

       'Positive' Class : 1

现在，就像在之前的步骤中一样，使用拟合模型在测试数据集上进行预测，并创建混淆矩阵以评估模型在测试数据上的性能：

print("Test data results -")
pred_test <- factor(ifelse(predict(xgb,data.matrix(xgb.test),
type="class")>0.5,1,0))
confusionMatrix(pred_test,factor(y_test),positive='1')

输出如下：

[1] "Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction     0     1
         0 25261  3884
         1  1154  3579

               Accuracy : 0.8513          
                 95% CI : (0.8475, 0.8551)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.5017          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.4796          
            Specificity : 0.9563          
         Pos Pred Value : 0.7562          
         Neg Pred Value : 0.8667          
             Prevalence : 0.2203          
         Detection Rate : 0.1056          
   Detection Prevalence : 0.1397          
      Balanced Accuracy : 0.7179          

       'Positive' Class : 1

如果我们仔细观察模型的结果，我们可以看到与随机森林模型结果相比，性能略有提升。从0.5提升到0.54，我们可以在保持比随机森林略高的召回率的同时提高精确度（以匹配随机森林）。XGBoost 的召回率提升幅度远大于精确度的下降。概率截止值的阈值不是一个固定的、硬性的截止值。我们可以根据我们的用例调整阈值。最佳数值可以通过经验实验或研究敏感性、特异性分布来研究。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/30gzSW0。

以下练习使用 0.54 而不是 0.5 作为概率截止值，以研究在牺牲召回率的情况下精确度的提升。

练习 75：提高 XGBoost 模型性能

我们可以通过调整输出的阈值来调整二元分类模型的性能。默认情况下，我们选择 0.5 作为默认概率截止值。因此，所有高于 0.5 的响应都被标记为Yes，否则为No。调整阈值可以帮助我们实现更敏感或更精确的模型。

通过调整概率截止阈值来提高 XGBoost 模型的性能，执行以下步骤：

将训练数据集上预测的概率截止阈值从 0.5 提高到 0.53 并打印结果：

print("Training data results -")
pred_train <- factor(ifelse(predict(xgb,data.matrix(xgb.train),
type="class")>0.53,1,0))
confusionMatrix(pred_train,factor(y_train),positive='1')

输出如下：

[1] "Training data results -"
Confusion Matrix and Statistics
          Reference
Prediction     0     1
         0 59626  9635
         1  1865  7921

               Accuracy : 0.8545        
                 95% CI : (0.852, 0.857)
    No Information Rate : 0.7779        
    P-Value [Acc > NIR] : < 2.2e-16     

                  Kappa : 0.4999        
 Mcnemar's Test P-Value : < 2.2e-16     

            Sensitivity : 0.4512        
            Specificity : 0.9697        
         Pos Pred Value : 0.8094        
         Neg Pred Value : 0.8609        
             Prevalence : 0.2221        
         Detection Rate : 0.1002        
   Detection Prevalence : 0.1238        
      Balanced Accuracy : 0.7104        

       'Positive' Class : 1

将测试数据集上预测的概率截止阈值从 0.5 提高到 0.53 并打印结果：

print("Test data results -")
pred_test <- factor(ifelse(predict(xgb,data.matrix(xgb.test),
type="class")>0.53,1,0))
confusionMatrix(pred_test,factor(y_test),positive='1')

输出如下：

1] "Test data results -"
Confusion Matrix and Statistics
          Reference
Prediction     0     1
         0 25551  4210
         1   864  3253

               Accuracy : 0.8502         
                 95% CI : (0.8464, 0.854)
    No Information Rate : 0.7797         
    P-Value [Acc > NIR] : < 2.2e-16      

                  Kappa : 0.4804         
 Mcnemar's Test P-Value : < 2.2e-16      

            Sensitivity : 0.43588        
            Specificity : 0.96729        
         Pos Pred Value : 0.79014        
         Neg Pred Value : 0.85854        
             Prevalence : 0.22029        
         Detection Rate : 0.09602        
   Detection Prevalence : 0.12152        
      Balanced Accuracy : 0.70159        

       'Positive' Class : 1

我们可以看到，在 44% 的召回率下，我们在测试数据集上达到了 80% 的精确率，训练集和测试集之间的性能差异也可以忽略不计。因此，我们可以得出结论，XGBoost 的模型性能略优于随机森林，尽管只是略好。

注意

你可以在 GitHub 上找到完整的代码：bit.ly/30c5DQ9。

在结束本章之前，让我们通过实验最后一种用于分类的监督技术，即深度神经网络。

深度神经网络

在结束本章之前，我们将讨论最后一种技术，即深度神经网络或深度学习。这是一个漫长且复杂的话题，我们绝对无法在章节的简短部分中公正地处理它。一本书可能甚至不足以覆盖这个话题的表面！我们将从高处探讨这个话题，并快速研究 R 中的简单实现。

深度神经网络主要应用于计算机视觉和自然语言处理领域，在机器学习的回归和分类用例中也具有重要意义，特别是在表格横截面数据上。随着大量数据的出现，深度神经网络已被证明在学习潜在模式并因此训练出性能更好的模型方面非常有效。

深度神经网络的深入探讨

深度神经网络受到人脑神经结构的启发。深度学习领域因解决计算机视觉问题而变得流行，即那些人类容易解决但计算机长期难以解决的问题领域。设计类似微型且高度简化的深度神经网络，旨在解决人类特别容易解决的问题。后来，随着深度学习在计算机视觉领域的成功，它被应用于其他几个领域，包括传统的机器学习监督用例。

神经网络组织成一个神经元层次结构，就像人脑中的神经元一样。每个神经元都与其他神经元相连，这使它们之间能够通过信号进行通信，形成一个可以学习并具有反馈机制的复杂网络。

下图展示了一个简单的神经网络：

图 5.9：简单的神经网络

输入数据构成了网络的第 0 层。这一层然后连接到下一层的神经元，下一层是隐藏的。它被称为隐藏的，因为网络可以被看作是一个黑盒，我们向网络提供输入并直接看到输出。中间层是隐藏的。在神经网络中，一层可以有任意数量的神经元，每个网络可以有任意数量的层。层的数量越多，网络就越“深”。因此得名深度学习和深度神经网络。每个隐藏层中的每个神经元都计算一个数学函数，在深度学习中这被称为激活函数。这个函数有助于模拟两个神经元之间的信号。如果函数（激活）计算出的值大于阈值，它就会向下一层直接连接的神经元发送信号。这两个神经元之间的连接由一个权重调节。权重决定了传入神经元的信号对接收神经元的重要性。深度学习模型中的学习方法通过更新神经元之间的权重，使得最终的预测，类似于机器学习模型，是最准确的。

深度学习模型是如何工作的？

要理解神经网络是如何工作以及如何学习在数据上做出预测的，让我们考虑一个对人类来说相对非常简单的任务。考虑通过人脸识别不同人的任务。我们大多数人每天都会遇到几个不同的人；比如说，在工作、学校或街道上。我们遇到的每个人在某些维度上都是不同的。尽管每个人都会有大量的相似特征，比如两只眼睛、两只耳朵、嘴唇、两只手等等，但我们的大脑可以轻易地区分两个人。当我们第二次遇到某个人时，我们很可能会认出他们，并将他们识别为我们之前遇到过的人。考虑到这种情况发生的规模以及我们的大脑能够轻松地解决这个巨大问题的现实，这让我们不禁想知道这究竟是如何发生的。

为了理解这一点并欣赏我们大脑的美丽，我们需要了解大脑是如何从根本上学习的。大脑是一个由相互连接的神经元组成的大型、复杂结构。当神经元感知到某些重要事物时，它会被激活，并将消息或信号传递给它连接的其他神经元。神经元之间的连接通过从反馈中不断学习而得到加强。在这里，当我们看到一张新面孔时，我们的大脑不是学习面孔的结构来识别人，而是学习给定的面孔与通用基准面孔的不同之处。这可以进一步简化为计算重要面部特征（如眼睛形状、鼻子、嘴唇、耳朵和唇部结构）之间的差异，以及皮肤和头发的颜色偏差和其他属性。这些差异由不同的神经元量化，然后以系统化的方式编排，以便大脑能够区分不同的面孔并从记忆中回忆起面孔。整个计算都是在潜意识中发生的，我们几乎意识不到这一点，因为结果对我们来说是一瞬间就能注意到的。

神经网络本质上试图以极其简化的形式模仿大脑的学习功能。神经元以分层的方式相互连接，并使用随机权重进行初始化。网络中的数学计算结合了所有神经元的输入，最终达到最终结果。最终结果的偏差（预测值）被量化为误差，并作为反馈提供给网络。基于误差，网络试图更新连接的权重，并尝试迭代地减少预测中的误差。经过几次迭代，网络以有序的方式更新其权重，从而学会识别模式以做出正确的预测。

我们使用什么框架进行深度学习模型？

目前，我们将使用 Keras for R 进行深度神经网络实验，以用于我们的分类用例。对于深度学习模型开发，我们需要编写大量的代码，这将构成网络的构建块。为了加快我们的进程，我们可以利用 Keras，这是一个提供深度学习组件整洁抽象的深度学习框架。Keras 有一个 R 接口，并且建立在低级深度学习框架之上。

今天的人工智能社区中可用的深度学习框架要么是低级的，要么是高级的。TensorFlow、Theano、PyTorch、PaddlePaddle 和 mxnet 等框架是低级框架，为深度学习模型提供基本构建块。使用低级框架为最终网络设计提供了大量的灵活性和定制性。然而，我们仍然需要编写相当多的代码才能使一个相对较大的网络工作。为了进一步简化，有一些高级框架可用，它们在低级框架之上工作，并在构建深度学习模型的过程中提供第二层抽象。Keras、Gluon 和 Lasagne 是一些利用上述低级框架作为后端并提供新 API 的框架，这使得整体开发过程变得容易得多。与直接使用 TensorFlow 等低级框架相比，这减少了灵活性，并为大多数网络提供了一个稳健的解决方案。对于我们的用例，我们可以直接利用 Keras 的 R 接口。

使用install.packages('keras')命令将安装 Keras 的 R 接口，并会自动安装 TensorFlow 作为 Keras 的低级后端。

在 Keras 中构建深度神经网络

要在 R 中利用 Keras，我们需要对我们的现有训练数据集进行额外的数据增强。在 R 中大多数机器学习函数中，我们可以直接将分类列编码为因子传递。然而，我们注意到 XGBoost 强制要求数据需要被转换成 one-hot 编码形式，因为它不会将数据内部转换为所需的格式。因此，我们使用了 R 中的dummyVars函数将训练和测试数据集转换为 one-hot 编码版本，这样我们数据集中就只有数值数据。在 Keras 中，我们需要将训练数据集作为矩阵而不是 DataFrame 来提供。因此，除了将数据转换为 one-hot 编码形式外，我们还需要将数据集转换为矩阵。

此外，还建议我们标准化、归一化或缩放所有输入维度。归一化过程将数据值重新缩放到 0 到 1 的范围。同样，标准化将数据缩放到均值为（μ）0 和标准差（σ）为 1（单位方差）。这种转换在机器学习中是一个很好的特性，因为某些算法从中受益并更好地学习。然而，在深度学习中，这种转换变得至关重要，因为如果我们提供一个所有维度都在不同范围或尺度上的输入训练数据集，模型的学习过程就会受到影响。这个问题背后的原因是神经元中使用的激活函数类型。

以下代码片段实现了 Keras 中的一个基本神经网络。在这里，我们使用一个具有三个层，每层有 250 个神经元的架构。找到正确的架构是一个经验过程，没有明确的指南。网络越深，拟合数据所需的计算就越多。以下代码片段中使用的数据集与 XGBoost 中使用的相同，并且已经有一元编码的形式。

练习 76：使用 R Keras 在 R 中构建深度神经网络

在这个练习中，我们将利用深度神经网络构建一个分类模型，用于与练习 13，提高 XGBoost 模型性能相同的用例，并尝试提高性能。深度神经网络并不总是比集成模型表现更好。当我们有非常高的训练样本数量时，例如 1000 万，它们通常是一个更好的选择。然而，我们将进行实验并检查我们是否能实现比我们在练习 10-13 中构建的模型更好的性能。

按照以下步骤在 R 中构建深度神经网络。

将输入数据集缩放到 0 到 1 的范围内。我们首先需要在训练数据上初始化一个preProcess对象。这将随后用于缩放训练数据以及测试数据。神经网络在缩放数据上表现更好。仅使用训练数据来创建缩放对象：
```
standardizer <- preProcess(x_train, method='range',rangeBounds=c(0,1))
```

使用之前步骤中创建的standardizer对象来缩放训练和测试数据：

x_train_scaled <- predict(standardizer, newdata=x_train)
x_test_scaled <- predict(standardizer, newdata=x_test)

将预测变量数量存储在一个名为predictors的变量中。我们将使用这些信息来构建网络：
```
predictors <-  dim(x_train_scaled)[2]
```
定义深度神经网络的架构。我们将使用keras_model_sequential方法。我们将创建一个包含三个隐藏层，每个层有 250 个神经元，激活函数为relu的网络。输出层将有一个神经元，激活函数为sigmoid（因为我们正在开发一个二元分类模型）：
```
dl_model <-  keras_model_sequential()  %>% 
  layer_dense(units = 250, activation = 'relu', 
input_shape =c(predictors)) %>% 
  layer_dense(units = 250, activation = 'relu' ) %>% 
  layer_dense(units = 250, activation = 'relu') %>% 
  layer_dense(units = 1, activation = 'sigmoid') 
```

定义模型优化器为adam，损失函数以及模型训练迭代中要捕获的指标：

dl_model %>% compile(
  loss = 'binary_crossentropy',
  optimizer = optimizer_adam(),
  metrics = c('accuracy')
)
summary(dl_model)

输出如下：

_____________________________________________________________
Layer (type)                 Output Shape               Param #
=============================================================
dense_34 (Dense)             (None, 250)                16750
_____________________________________________________________
dense_35 (Dense)             (None, 250)                62750
_____________________________________________________________
dense_36 (Dense)             (None, 250)                62750
_____________________________________________________________
dense_37 (Dense)             (None, 1)                  251  
=============================================================
Total params: 142,501
Trainable params: 142,501
Non-trainable params: 0

使用步骤 4-5 中创建的模型结构，以及步骤 1-2 中的训练和测试数据来拟合模型：

history <- dl_model %>% fit(
  as.matrix(x_train_scaled), as.matrix(y_train), 
  epochs = 10, batch_size = 32, 
  validation_split = 0.2
)

输出如下：

Train on 63237 samples, validate on 15810 samples
Epoch 1/10
63237/63237 [==============================] - 7s 104us/step – 
loss: 0.3723 - acc: 0.8388 - val_loss: 0.3639 - val_acc: 0.8426
Epoch 2/10
63237/63237 [==============================] - 6s 102us/step – 
loss: 0.3498 - acc: 0.8492 - val_loss: 0.3695 - val_acc: 0.8380
Epoch 3/10
63237/63237 [==============================] - 6s 97us/step – 
loss: 0.3434 - acc: 0.8518 - val_loss: 0.3660 - val_acc: 0.8438
Epoch 4/10
63237/63237 [==============================] - 6s 99us/step – 
loss: 0.3390 - acc: 0.8527 - val_loss: 0.3628 - val_acc: 0.8395
Epoch 5/10
63237/63237 [==============================] - 6s 97us/step – 
loss: 0.3340 - acc: 0.8551 - val_loss: 0.3556 - val_acc: 0.8440
Epoch 6/10
63237/63237 [==============================] - 7s 119us/step – 
loss: 0.3311 - acc: 0.8574 - val_loss: 0.3612 - val_acc: 0.8414
Epoch 7/10
63237/63237 [==============================] - 7s 107us/step – 
loss: 0.3266 - acc: 0.8573 - val_loss: 0.3536 - val_acc: 0.8469
Epoch 8/10
63237/63237 [==============================] - 7s 105us/step – 
loss: 0.3224 - acc: 0.8593 - val_loss: 0.3575 - val_acc: 0.8471
Epoch 9/10
63237/63237 [==============================] - 7s 105us/step – 
loss: 0.3181 - acc: 0.8607 - val_loss: 0.3755 - val_acc: 0.8444
Epoch 10/10
63237/63237 [==============================] - 7s 104us/step – 
loss: 0.3133 - acc: 0.8631 - val_loss: 0.3601 - val_acc: 0.8468

使用拟合的模型在训练数据集上预测响应：

print("Training data results - ")
pred_train <- factor(ifelse(predict(dl_model, 
as.matrix(x_train_scaled))>0.5,1,0))
confusionMatrix(pred_train,factor(y_train),positive='1')

输出如下：

"Training data results - "
Confusion Matrix and Statistics
          Reference
Prediction     0     1
         0 59281  8415
         1  2351  9000

               Accuracy : 0.8638          
                 95% CI : (0.8614, 0.8662)
    No Information Rate : 0.7797          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.547           
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.5168          
            Specificity : 0.9619          
         Pos Pred Value : 0.7929          
         Neg Pred Value : 0.8757          
             Prevalence : 0.2203          
         Detection Rate : 0.1139          
   Detection Prevalence : 0.1436          
      Balanced Accuracy : 0.7393          

       'Positive' Class : 1

使用拟合的模型在测试数据集上预测响应：

#Predict on Test Data
pred_test <- factor(ifelse(predict(dl_model,   
                           as.matrix(x_test_scaled))>0.5,1,0))
confusionMatrix(pred_test,factor(y_test),positive='1')

输出如下：

"Test data results - "
Confusion Matrix and Statistics
          Reference
Prediction     0     1
         0 25028  3944
         1  1246  3660

               Accuracy : 0.8468          
                 95% CI : (0.8429, 0.8506)
    No Information Rate : 0.7755          
    P-Value [Acc > NIR] : < 2.2e-16       

                  Kappa : 0.4965          
 Mcnemar's Test P-Value : < 2.2e-16       

            Sensitivity : 0.4813          
            Specificity : 0.9526          
         Pos Pred Value : 0.7460          
         Neg Pred Value : 0.8639          
             Prevalence : 0.2245          
         Detection Rate : 0.1080          
   Detection Prevalence : 0.1448          
      Balanced Accuracy : 0.7170          
       'Positive' Class : 1

注意

你可以在 GitHub 上找到完整的代码：bit.ly/2Vz8Omb。

预处理函数有助于将数据转换成所需的尺度或范围。在这里，我们将数据缩放到 0 到 1 的尺度。我们应仅考虑使用训练数据作为函数生成器的输入，并使用拟合方法来缩放测试数据。这在实际场景中非常重要，因为我们无法实时访问测试数据。一旦preProcess方法拟合完成，我们就用它来转换训练和测试数据。然后我们定义深度神经网络模型的架构。R 提供了易于扩展的管道操作符%>%，它使得操作符的连接变得简单。我们设计了一个包含三个层次，每个层次有 250 个神经元的网络。输入数据将形成第 0 层，最后一层将是预测结果。网络中用于隐藏层的激活函数是relu，这是任何深度学习用例中最推荐的激活函数。最后一层使用sigmoid激活函数，因为我们有一个二元分类用例。在 Keras 中可以选择许多激活函数，例如prelu、tanh、swish等等。一旦模型架构定义完成，我们就定义损失函数binary_crossentropy，它与二元的logloss（类似于 XGBoost）类似，是模型用来学习和反向传播的技术，即优化器。预测中的错误会被反向传播到网络中，以便它可以调整权重，并迭代地减少错误。

这个功能的数学直观性可以采取多种方法。基于低阶矩自适应估计的 Adam 优化是最受欢迎的选择，对于大多数深度学习用例，我们几乎可以盲目地进行实验。其他一些选项包括rmsprop、随机梯度下降和Adagrad。我们还定义了在每个 epoch 后要在验证数据集上计算的指标，即网络对训练样本的一次完整展示。summary函数显示了使用 Keras 构造在前一节定义的架构的结果。summary函数给我们一个关于每个层中参数数量的简要概念，并且以层次结构的形式表示网络，帮助我们可视化模型架构。最后，我们使用fit函数来训练或“拟合”数据到网络中。我们还定义了模型应该迭代的 epoch 数；epoch 数越高，训练过程将需要更长的时间来计算。

批次大小表示网络在更新网络权重之前在一次单次传递中消耗的训练样本数量；批次较小的数字表示更频繁的权重更新，有助于有效地利用 RAM 内存。验证分割定义了每个 epoch 结束时用于验证的训练样本的百分比。最后，我们在训练数据和测试数据上验证模型的表现。

注意

代码片段中的这种解释无论如何都不能成为主题的正当理由。深度神经网络是一个极其庞大且复杂的主题，可能需要一本书来介绍基础知识。我们已将上下文封装成一段简短的段落，以便您了解模型开发过程中使用的结构。探索这个主题的深度将超出本书的范围。

从结果来看，我们可以看到与之前模型相似的结果。结果几乎可以与之前开发的 XGBoost 模型相媲美。在测试数据集上，我们大约有 48%的召回率和 75%的精确率。结果可以进一步调整以降低召回率并提高精确率（如果需要）。

因此，我们可以得出结论，我们从简单的逻辑回归模型、XGBoost 和深度神经网络模型中得到了相当好的结果。这三个模型之间的差异相对较小。这可能会让你产生重要的问题：在相同的使用案例上迭代各种模型是否值得？哪种模型理论上会给出最佳结果？尽管对这些问题的答案并不直接，但我们可以这样说，总的来说，简单模型总是表现良好；集成模型在大量数据的情况下表现更佳；而深度学习模型在大量数据的情况下表现更佳。在本章中，我们实验的使用案例中，通过超参数调整和；最重要的是；特征工程，我们将从所有模型中获得改进的结果。我们将在第七章“模型改进”中探讨超参数调整，在第六章“特征选择和降维”中探讨轻节点上的特征工程。特征工程的过程非常特定于领域，只能在一定程度上进行概括。我们将在下一章中更详细地探讨这一点。本章的主要议程是介绍涵盖该领域大量建模技术的范围，并帮助您为任何用于分类用例的机器学习技术打下基础。

为您的用例选择正确的模型

到目前为止，我们已经探索了一组白盒模型和几个用于相同分类用例的黑盒机器学习模型。我们还用 Keras 扩展了相同的用例，并研究了其性能。从几个模型和多次迭代的结果来看，我们需要决定哪个模型最适合分类用例。对此并没有简单直接的答案。在更广泛的意义上，我们可以认为对于大多数用例，最佳模型将是随机森林或 XGBoost。然而，这并不适用于所有类型的数据。会有许多场景，集成建模可能并不合适，线性模型会优于它，反之亦然。在数据科学家为分类用例进行的多数实验中，方法将是探索性和迭代的。在机器学习中没有一种适合所有情况的模型。设计和训练机器学习模型的过程是艰巨且极其迭代的，并且始终取决于用于训练的数据类型。

在给定构建监督机器学习模型的任务的情况下，最佳的前进方法如下：

步骤 0: 探索性数据分析、数据处理和特征工程：使用可视化技术的组合对数据进行深入研究，然后处理缺失值、去除异常值、构建新特征，并建立训练集和测试集。（如果需要，也可以创建一个验证集。）
步骤 1: 从简单的白盒模型如逻辑回归开始：在建模迭代中，最佳起点是一个简单的白盒模型，它可以帮助我们以易于量化的方式研究每个预测变量对因变量的影响。几次模型迭代将有助于特征选择，并清晰地理解最佳预测变量和模型基准。
步骤 2: 使用决策树模型重复建模实验：利用决策树模型将始终帮助我们获得对模型和特征模式的全新视角。它可能会给出简单的规则，从而为改进模型提供新的特征工程思路。
步骤 3: 如果有足够的数据，尝试集成建模；否则，尝试其他方法，例如支持向量机。

使用随机森林和 XGBoost 进行集成建模几乎总是实验的一个安全选项。但在数据稀缺的情况下进行训练时，集成建模可能不是一种有效的推进方法。在这种情况下，基于黑盒核的模型在学习和数据模式方面可能更有效，从而提高模型性能。鉴于范围，我们没有在本章中涵盖支持向量机（SVM）。然而，考虑到本章涵盖的广泛主题，对于您来说，开始使用 SVM 将是一个简单直接的任务。本博客提供了一个简单易懂的 SVM 指南：eight2late.wordpress.com/2017/02/07/….

此外，为了了解训练样本的数量是更多还是更少，你可以使用一个简单的经验法则。如果数据集中每个特征至少有 100 个训练样本行，那么对于集成模型来说，数据就足够了；如果样本数量低于那个水平，那么集成模型可能并不总是有效的。尽管如此，尝试一下仍然值得。例如，如果有 15 个特征（自变量）和 1 个因变量，那么如果我们有15 x 100 = 1500个训练样本，集成模型在白盒模型上可能会有更好的性能。
步骤 4：如果数据量足够多，尝试深度神经网络。如果数据集中每个特征至少有 10,000 个样本，那么尝试深度神经网络可能是个好主意。神经网络的问题主要是需要大量的训练数据和大量的迭代才能获得良好的性能。在大多数通用情况下，对于使用表格横截面数据进行分类（本书中我们解决的那种用例），深度神经网络与集成模型一样有效，但需要显著更多的努力来训练和调整以达到相同的结果。当有显著大量的样本进行训练时，它们确实优于集成模型。只有在有显著更多的训练样本时，在深度神经网络上的努力才会带来有利的结果。

摘要

在本章中，我们探讨了监督机器学习中不同类型的分类算法。我们利用澳大利亚天气数据，围绕它设计了一个商业问题，并探索了同一用例上的各种机器学习技术。我们研究了如何在 R 中开发这些模型，并深入研究了这些算法的数学抽象功能。我们总结了每种技术的结果，并研究了处理常见分类用例的通用方法。

在下一章中，我们将研究机器学习模型的特征选择、降维和特征工程。

R-应用监督学习-二-