EDA

介绍

本章将向你展示如何使用可视化和转换来系统地探索数据，这是一项统计学家称为探索性数据分析（Exploratory Data Analysis, 简称 EDA）的任务。EDA 是一个迭代循环，其过程包括：

1. 对数据生成问题。

2. 通过可视化、转换和建模来寻找答案。

3. 利用所得的知识来完善你的问题，或生成新的问题。

EDA 并不是一个具有严格规则的正式过程。更确切地说，EDA 是一种思维方式。在 EDA 的初始阶段，你应该自由地探索所有想到的想法。其中有些想法会取得成果，而有些则会走向死胡同。随着探索的深入，你会逐渐聚焦于一些特别有价值的洞察，最终将这些洞察整理成报告并与他人交流。

即使研究的核心问题已经明确地摆在你面前，EDA 仍然是任何数据分析的重要组成部分，因为你始终需要调查数据的质量。数据清洗就是 EDA 的一个应用场景：你需要提出问题，检查数据是否符合你的预期。要进行数据清洗，你需要用到 EDA 的所有工具：可视化、转换和建模。

在进行探索性数据分析（EDA）时，你的目标是深入了解你的数据。最简单的方法是将问题作为工具来引导你的调查。当你提出一个问题时，这个问题会将你的注意力集中在数据集的某个特定部分，并帮助你决定应该绘制哪些图表、构建哪些模型或进行哪些转换。

EDA 本质上是一个创造性的过程。就像大多数创造性过程一样，提出高质量问题的关键在于生成大量的问题。在分析开始时，提出具有洞察力的问题通常很困难，因为你尚不了解数据集中可以发现哪些见解。另一方面，每个新问题都会让你接触到数据的新方面，并增加你发现有价值信息的机会。如果你在每个问题的基础上提出新的问题，你可以快速深入到数据中最有趣的部分，并逐步形成一系列引人深思的问题。

关于如何提出问题来引导研究，并没有固定的规则。然而，有两类问题始终对从数据中发现信息非常有用。可以粗略地将这两类问题表述为：

1. 在我的变量中，存在哪些类型的 变化（variation）？

2. 在我的变量之间，存在哪些类型的 协变化（covariation）？

本章的其余部分将讨论这两个问题。我们将解释什么是变化和协变化，并向你展示多种方法

Variation

变化（Variation） 是指变量的取值在不同测量中发生变化的趋势。在现实生活中，你可以轻松观察到变化的存在；如果你对某个连续变量进行两次测量，结果通常会有所不同。即使是对像光速这样常量的数量进行测量，这种现象依然适用。每次测量都会包含一定量的误差，这些误差会在不同的测量中有所变化。

变量也可能因为测量的对象不同（例如，不同人群的眼睛颜色）或测量的时间不同（例如，不同时刻电子的能量水平）而发生变化。每个变量都有其独特的变化模式，这些模式可以揭示变量在同一观测对象的多次测量中如何变化，以及在不同观测对象之间如何变化的有趣信息。

我们将通过可视化约 54,000 颗钻石的重量（克拉）的分布，开始我们的探索。由于克拉是一个数值变量，我们可以使用直方图来呈现其分布情况：

library(tidyverse)
library(tidymodels)
ggplot(diamonds, aes(x = carat)) +
  geom_histogram(binwidth = 0.5)

现在你已经能够可视化变量的变化，那么在图表中应该关注什么？又应该提出哪些后续问题？以下是我们整理的一份清单，其中列出了图表中最有用的信息类型，以及每种信息对应的一些后续问题。提出好的后续问题的关键在于依靠你的好奇心（你还想了解更多什么？）以及批判性思维（这可能有哪些误导之处？）。

典型值

在柱状图和直方图中，较高的柱子表示变量的常见值，而较矮的柱子表示较少见的值。没有柱子的地方则显示数据中未出现的值。为了将这些信息转化为有用的问题，可以关注任何出乎意料的地方：

• 哪些值是最常见的？为什么？

• 哪些值很少见？为什么？这是否符合你的预期？

• 你能看到任何不寻常的模式吗？可能是什么原因造成的？

现在让我们看看小克拉钻石的重量分布情况。

smaller <- diamonds |> 
  filter(carat < 3)

ggplot(smaller, aes(x = carat)) +
  geom_histogram(binwidth = 0.01)

这个直方图引发了几个有趣的问题：

• 为什么在整数克拉和常见的分数克拉处有更多的钻石？

• 为什么在每个峰值的右侧稍微多一些钻石，而在左侧稍微少一些？

可视化还能够揭示聚类现象，这表明数据中可能存在子群体。为了理解这些子群体，可以提出以下问题：

• 每个子群体中的观测值有哪些相似之处？

• 不同聚类之间的观测值有哪些差异？

• 你如何解释或描述这些聚类？

• 为什么聚类的出现可能具有误导性？

这些问题中有一些可以通过数据直接回答，而另一些可能需要依赖对数据的领域知识。此外，许多问题会促使你探索变量之间的关系，例如，查看一个变量的值是否可以解释另一个变量的行为。我们很快就会深入研究这一点。

离群值

离群值是指异常的观测值，即那些似乎不符合数据模式的数据点。有时，离群值可能是数据录入错误，有时它们只是极端值，在数据收集中恰好被观察到了，而有时它们可能暗示着重要的新发现。

当你拥有大量数据时，离群值有时很难在直方图中被观察到。例如，查看 diamonds 数据集中 y 变量的分布时，离群值的唯一迹象可能是 x 轴范围异常宽。

ggplot(diamonds, aes(x = y)) + 
  geom_histogram(binwidth = 0.5)

在常见的区间中有大量观测值，这使得罕见区间的柱子非常矮，难以察觉（尽管如果你仔细盯着 0 的位置，也许会发现一些异常）。为了更容易观察这些不寻常的值，我们需要通过 coord_cartesian() 对 y 轴的小值范围进行缩放：

ggplot(diamonds, aes(x = y)) + 
  geom_histogram(binwidth = 0.5) +
  coord_cartesian(ylim = c(0, 50))

coord_cartesian() 的 xlim() 参数可以用来缩放 x 轴，而 ggplot2 的 xlim() 和 ylim() 函数则稍有不同：它们会直接丢弃超出范围的数据点。

通过缩放，我们可以看到三个不寻常的值：0、大约 30 和大约 60。我们可以使用 dplyr 来提取这些值进行进一步分析。

unusual <- diamonds |> 
  filter(y < 3 | y > 20) |> 
  select(price, x, y, z) |>
  arrange(y)
unusual

y 变量表示钻石的一个维度（以毫米为单位）。我们知道钻石的宽度不可能是 0mm，所以这些值必定是错误的。通过探索性数据分析（EDA），我们发现了用 0 代替的缺失数据。这种错误如果仅仅通过寻找 NA 是不容易发现的。接下来，我们可以选择将这些值重新编码为 NA，从而避免误导性的计算。

此外，我们还可能怀疑 32mm 和 59mm 的测量值是不现实的：这些钻石的长度超过一英寸，但它们的价格却并没有达到数十万美元的水平，这与常识不符。

处理离群值的建议

1. 重复分析：
   在有离群值和没有离群值的情况下分别进行分析。如果离群值对结果的影响很小，并且无法确定它们的来源，那么可以合理地将它们排除并继续分析。

2. 谨慎处理：
   如果离群值对结果有显著影响，你不应该在没有正当理由的情况下删除它们。你需要查明离群值的原因（例如，是否是数据录入错误），并在报告中披露你移除它们的理由。

通过这种方式，你可以确保数据分析的准确性，同时最大程度地避免误导性的结论。

离群值的处理

如果你在数据集中遇到了不寻常的值，并且想直接继续后续的分析，你有两个选项可以处理这些值：

1. 删除包含异常值的整行数据

你可以选择完全删除这些包含异常值的观测行。以下是使用 dplyr 的方法示例：

# 过滤掉 y 值等于 0、30 或 60 的行
diamonds2 <- diamonds |> 
  filter(between(y, 3, 20))

这种方法简单直接，适用于异常值较少且对整体分析影响不大的情况。我们不推荐删除整行数据的选项，因为一个无效值并不意味着该观测的其他值也同样无效。此外，如果你的数据质量较差，那么对每个变量都应用这种方法后，可能会发现几乎没有数据剩下可供分析了！

2. 将异常值替换为缺失值 (NA)

我们推荐将不寻常的值替换为缺失值（NA）。最简单的方法是使用 mutate() 对变量进行修改，生成一个经过处理的副本。可以用 if_else() 函数将不寻常的值替换为 NA。

# 将异常值替换为 NA
diamonds2 <- diamonds |> 
  mutate(y = if_else(y < 3 | y > 20, NA, y))

为什么推荐这种方法？

1. 保留更多数据：
   只有异常值会被替换为 NA，而不是丢弃整行数据，从而最大程度保留其他有效信息。

2. 灵活处理缺失值：
   后续分析中，你可以根据需要选择忽略含有 NA 的值，或选择合适的填充方法对缺失值进行补充。

3. 避免数据偏差：
   如果丢弃整行数据，可能会导致样本量减少，甚至引入偏差，特别是在数据本身质量较低的情况下。

在绘制数据时，缺失值的位置并不显而易见。因此，ggplot2 默认不会在图中显示缺失值，但会发出警告提示这些缺失值已被移除。

例如，当你尝试绘制包含缺失值的数据时，可能会看到类似这样的警告：

ggplot(diamonds2, aes(x = x, y = y)) + 
  geom_point()

要抑制此类警告，可以在绘图函数中设置 na.rm = TRUE，例如：

ggplot(diamonds2, aes(x = x, y = y)) + 
  geom_point(na.rm = TRUE)

有时，你可能需要深入了解缺失值的模式，比如分析哪些观测值缺失，或者缺失值的观测与非缺失值的观测有何不同。以飞行数据集 nycflights13::flights 为例，其中 dep_time 列的缺失值表示航班被取消。

有时，你可能希望了解缺失值与非缺失值的观测之间的差异。例如，在 nycflights13::flights 数据集中，dep_time 的缺失值表示航班被取消。可以通过创建一个新变量来标记缺失值，具体方法如下：

nycflights13::flights |> 
  mutate(
    cancelled = is.na(dep_time),
    sched_hour = sched_dep_time %/% 100,
    sched_min = sched_dep_time %% 100,
    sched_dep_time = sched_hour + (sched_min / 60)
  ) |> 
  ggplot(aes(x = sched_dep_time)) + 
  geom_freqpoly(aes(color = cancelled), binwidth = 1/4)

这样就可以区分被取消和未取消的航班，为后续分析提供支持。

然而，这个图表并不理想，因为未取消航班的数量远远多于被取消的航班。在下一节中，我们将探讨一些改进此类比较的技巧。

Covariation

如果说变化（variation）描述的是单个变量的行为，那么协变化（covariation）描述的就是多个变量之间的关联行为。协变化是指两个或多个变量的取值以某种相关方式一起变化的趋势。发现协变化的最佳方法是可视化变量之间的关系。

离散变量与数值变量

例如，我们可以使用 geom_freqpoly() 来探索钻石的价格如何随其质量（通过 cut 来衡量）而变化：

ggplot(diamonds, aes(x = price)) + 
  geom_freqpoly(aes(color = cut), binwidth = 500, linewidth = 0.75)

注意到，ggplot2 为 cut 使用了一个有序的颜色比例，因为在数据中，cut 被定义为有序因子变量（ordered factor）。

然而，geom_freqpoly() 的默认外观在这里并不是很有用，因为高度（由总计数决定）在不同的切工质量之间差异很大，这使得很难看出它们分布形状之间的差异。

为了让比较更容易，我们需要调整 y 轴的显示内容。与其显示计数（count），我们可以显示密度（density），它是经过标准化的计数，使每个频率多边形下的面积为 1。

ggplot(diamonds, aes(x = price, y = after_stat(density))) + 
  geom_freqpoly(aes(color = cut), binwidth = 500, linewidth = 0.75)

注意，我们将密度映射到 y 轴，但由于密度并不是 diamonds 数据集中的一个现有变量，我们需要先计算它。为此，可以使用 after_stat() 函数进行计算。

在这个图中，有一个看起来相当令人惊讶的现象——低质量的切工（Fair）钻石似乎拥有最高的平均价格！ 但这可能是因为频率多边形（frequency polygons）不太容易解释——这个图中包含的信息量太多，导致解读出现困难。

探索这种关系的一种视觉上更简单的方式是使用并排箱线图（side-by-side boxplots）。

ggplot(diamonds, aes(x = cut, y = price)) +
  geom_boxplot()

箱线图虽然提供的分布信息较少，但更加紧凑，因此更容易进行比较（并且可以在一张图中容纳更多内容）。通过箱线图，我们可以支持这个反直觉的发现：质量更好的钻石通常价格更便宜！

cut 是一个有序因子变量：fair 的质量低于 good，而 good 又低于 very good，依此类推。然而，许多分类变量并没有这种内在的顺序，因此你可能需要重新排序这些变量以创建更具信息量的可视化展示。一种实现这种操作的方法是使用 fct_reorder()。

例如，在 mpg 数据集中，如果你想知道高速公路油耗（hwy）如何随汽车类别（class）变化，可以这样做：

ggplot(mpg, aes(x = class, y = hwy)) +
  geom_boxplot()

为了让趋势更清晰，我们可以根据 hwy（高速公路油耗）的中位数值对 class 进行重新排序。这会将车辆类别按高速公路油耗的中位数从低到高进行排列，从而更直观地展示趋势。

ggplot(mpg, aes(x = fct_reorder(class, hwy, median), y = hwy)) +
  geom_boxplot()

如果你的变量名很长，使用 geom_boxplot() 时，将图表旋转 90° 会使其更易于阅读。你可以通过交换 x 和 y 的映射来实现这一点。

ggplot(mpg, aes(x = hwy, y = fct_reorder(class, hwy, median))) +
  geom_boxplot()

练习

1. 根据探索性数据分析（EDA），在 diamonds 数据集中，哪个变量对预测钻石价格最重要？该变量与 cut 的关系如何？为什么这两个关系的结合导致低质量的钻石更贵？

2. 箱线图的一个问题是，它们是在数据集较小的时代开发的，往往会显示大量“异常值”。解决此问题的一种方法是使用字母值图（letter value plot）。安装 lvplot 包，并尝试使用 geom_lv() 来显示价格与 cut 的分布。你学到了什么？如何解释这些图？

library(lvplot)

ggplot(diamonds, aes(x = cut, y = price, fill = cut)) +
  geom_lv(lv = 10, outlier.colour = "red") +
  theme_minimal()

3. 使用 geom_violin() 创建一个钻石价格与分类变量的可视化，然后创建一个分面直方图（geom_histogram()）、一个彩色频率多边形（geom_freqpoly()）以及一个彩色密度图（geom_density()）。比较和对比这四种图。根据分类变量的水平来可视化一个数值变量的分布时，这些方法各有哪些优缺点？

ggplot(diamonds, aes(x = cut, y = price, fill = cut)) +
  geom_violin() +
  theme_minimal()

两个离散变量

对于两个离散变量之间的协变关系的可视化，需要计算每种分类变量水平组合的观测数量。实现这一点的一种方法是使用内置的 geom_count()：

ggplot(diamonds, aes(x = cut, y = color)) +
  geom_count()

图中每个圆的大小表示每种值组合的观测数量。协变关系将表现为特定的 x 值和特定的 y 值之间的强相关性。

另一种探索这些变量关系的方法是使用 dplyr 计算计数：

diamonds |> 
  count(color, cut)

然后使用 geom_tile() 和 fill 进行可视化：

diamonds |> 
  count(color, cut) |>  
  ggplot(aes(x = color, y = cut)) +
  geom_tile(aes(fill = n))

两个连续性变量

你已经见过一种很好的方法来可视化两个数值变量之间的协变关系：使用 geom_point() 绘制散点图。你可以通过点的分布模式观察协变关系。例如，可以看到钻石的克拉大小和价格之间存在正相关关系：克拉数更大的钻石价格更高。这种关系是指数型的。

ggplot(smaller, aes(x = carat, y = price)) +
  geom_point()

随着数据集规模的增长，散点图的实用性会降低，因为点会开始重叠并堆积成均匀的黑色区域，这使得很难判断数据在二维空间中的密度差异，也很难发现趋势。你已经见过一种解决该问题的方法：使用 alpha 属性来增加透明度。

ggplot(smaller, aes(x = carat, y = price)) + 
  geom_point(alpha = 1 / 100)

但是，对于非常大的数据集，使用透明度可能会有一定挑战。另一种解决方案是使用分箱（binning）。此前，你已经使用过 geom_histogram() 和 geom_freqpoly() 对一维数据进行分箱。现在，你将学习如何使用 geom_bin2d() 和 geom_hex() 对二维数据进行分箱。

geom_bin2d() 和 geom_hex() 将坐标平面划分为二维的箱（bins），然后通过填充颜色（fill）来表示每个箱中包含的点的数量。

• geom_bin2d() 创建的是矩形箱。

• geom_hex() 创建的是六边形箱。

要使用 geom_hex()，需要先安装 hexbin 包。

ggplot(smaller, aes(x = carat, y = price)) +
  geom_bin2d()


ggplot(smaller, aes(x = carat, y = price)) +
  geom_hex()

另一种选择是对一个连续变量进行分箱，使其表现得像一个分类变量。然后，你可以使用之前学到的可视化分类变量和连续变量组合的方法之一。例如，你可以对 carat（克拉）进行分箱，然后为每一组绘制一个箱线图（boxplot）：

ggplot(smaller, aes(x = carat, y = price)) + 
  geom_boxplot(aes(group = cut_width(carat, 0.1)))

cut_width(x, width)（如上所用）将变量 x 按指定的宽度 width 分箱。默认情况下，箱线图的外观大致相同（除了离群点的数量），因此很难看出每个箱线图总结的点数是否不同。

Patterns and Models

如果两个变量之间存在系统性关系，它会在数据中表现为一种模式（pattern）。当你发现模式时，可以问自己以下问题：

• 这个模式可能是由于巧合（即随机性）导致的吗？

• 你如何描述模式所暗示的关系？

• 模式所暗示的关系有多强？

• 还有哪些变量可能会影响这种关系？

• 如果你观察数据的个别子群体，这种关系是否会发生变化？

数据中的模式为变量之间的关系提供了线索，即揭示了协变化（covariation）。如果将变化（variation）视为一种导致不确定性的现象，那么协变化则是一种减少不确定性的现象。如果两个变量协变，你可以利用一个变量的值更好地预测另一个变量的值。如果协变化是由于因果关系（一个特殊情况）引起的，那么你可以利用一个变量的值来控制另一个变量的值。

模型 是从数据中提取模式的工具。

模型是从数据中提取模式的工具。例如，考虑钻石数据。由于 cut（切工）和 price（价格），以及 carat（克拉）和 price 之间的紧密关系，很难理解 cut 和 price 的关系。可以使用模型来移除 price 和 carat 之间的强关系，从而探索剩余的细微差异。

以下代码拟合了一个从 carat 预测 price 的模型，然后计算残差（预测值与实际值之间的差异）。残差让我们能够观察钻石的价格，而不受 carat 影响的影响。请注意，在使用 price 和 carat 的原始值之前，我们先对它们进行对数变换（log transform），然后对对数变换后的值拟合模型。最后，我们对残差进行指数化操作（exponentiate），以将其恢复到原始价.

library(tidymodels)

diamonds <- diamonds |>
  mutate(
    log_price = log(price),
    log_carat = log(carat)
  )

diamonds_fit <- linear_reg() |>
  fit(log_price ~ log_carat, data = diamonds)

diamonds_aug <- augment(diamonds_fit, new_data = diamonds) |>
  mutate(.resid = exp(.resid))

ggplot(diamonds_aug, aes(x = carat, y = .resid)) + 
  geom_point()

一旦移除了 carat（克拉）和 price（价格）之间的强关系，你就可以观察到 cut（切工）和 price 之间的预期关系：相对于其大小，高质量的钻石更昂贵。

ggplot(diamonds_aug, aes(x = cut, y = .resid)) + 
  geom_boxplot()

练习

基于某公司应收账款数据预测坏账准备。导入real_world_credit_sales.csv表格

library(tidyverse)
credit_sales_journal <- read_csv('AR_case/credit_sales_journal.csv')

在应收账款的数据中，我们最关注的变量是collection delay，也就是账龄，一个应收账款要多久才能被收回。

计算出collection-delay

credit_sales_journal <- credit_sales_journal |> 
  mutate(
    collection_delay = ymd(collection_date) -
          ymd(invoice_date)
  )

这段代码在计算收款延迟天数：

ymd() 是 lubridate 包的函数，用于将日期字符串转换为日期类型（Date）

• ymd(collection_date) 将收款日期转换为日期类型

• ymd(invoice_date) 将发票日期转换为日期类型

两个日期相减 - 会得到它们之间的天数差（difftime 对象）

• ymd(collection_date) - ymd(invoice_date) 计算收款日期与发票日期之间相差的天数

2. 探索并展示collection_delay的分布，这个分布有哪些特点，有哪些值得关注的异常值。

• 查看分布

ggplot(credit_sales_journal, aes(x = collection_delay)) +
  geom_histogram(binwidth = 10)

• 可以发现，有一些潜在的，收款期很长的应收款。

ggplot(credit_sales_journal, aes(x = collection_delay)) +
  geom_histogram(binwidth = 10)+
  coord_cartesian(ylim = c(0, 20))

• 把超过200天的视为异常值

credit_sales_journal <- credit_sales_journal |> 
  mutate(unusual = collection_delay >=200)

• 或者根据3-sigma rule来筛选，计算均值和标准差，并筛选出异常值

credit_sales_journal <- credit_sales_journal |> 
  mutate(
    mean_value = mean(collection_delay, na.rm = TRUE),          # 计算均值
    sd_value = sd(collection_delay, na.rm = TRUE),              # 计算标准差
    z_score = (collection_delay - mean_value) / sd_value,       # 计算 Z-score
    is_outlier = abs(z_score) > 3                    # 判断是否为异常值
  )

# 看一下异常值的分布

credit_sales_journal |> 
  filter(unusual == TRUE) |> 
  ggplot(aes(x = collection_delay)) +
  geom_histogram(binwidth = 10)


credit_sales_journal |> 
  filter(is_outlier == TRUE) |> 
  ggplot(aes(x = collection_delay)) +
  geom_histogram(binwidth = 10)

• 可以看到，两者筛选出来的异常值基本一致，因此，后续采用is_outlier来作为异常值

  探索collection_delay和数据中其他变量的协变关系

• 不同客户的还款分布

ggplot(credit_sales_journal, aes(x = collection_delay)) + 
  geom_freqpoly(aes(color = customer_no), binwidth = 10, linewidth = 0.75)

# 不同客户的异常还款

credit_sales_journal |> 
  filter(is_outlier == TRUE) |> 
  ggplot(aes(x = collection_delay)) + 
  geom_freqpoly(aes(color = customer_no), binwidth = 10, linewidth = 0.75)

大致还是比较均匀的，也就是说，每个客户直接还款市场、异常还款的分布情况都差不多。

# 在sales_return的分别情况
ggplot(credit_sales_journal, aes(x = sales_return, y = collection_delay)) +
  geom_boxplot(aes(group = sales_return))

• 还款期限是否有季节性呢？

ggplot(credit_sales_journal, aes(x = month(invoice_date), 
                                 y = collection_delay)) +
  geom_boxplot(aes(group = month(invoice_date)))

月份的波动似乎也很正常

• 和一些连续变量的关系

# 和单价sales_unit_price的关系

ggplot(credit_sales_journal, aes(x = sales_unit_price, y = collection_delay)) +
  geom_point(aes(color = is_outlier))+
  geom_smooth(method = "lm", color = "red")

单价在80以上的，超长未还款的情况似乎要少一些

# 和数量sales_count

ggplot(credit_sales_journal, aes(x = sales_count, y = collection_delay)) +
  geom_point(aes(color = is_outlier))+
  geom_smooth(method = "lm", color = "red")

# 和总价sales_extended
ggplot(credit_sales_journal, aes(x = sales_extended, y = collection_delay)) +
  geom_point(aes(color = is_outlier))+
  geom_smooth(method = "lm", color = "red")

ggplot(credit_sales_journal, aes(x = sales_extended, y = collection_delay)) +
  geom_point(aes(color = is_outlier))+
  geom_smooth(method = "lm", color = "red")

基本没有太大的相关性。

构建模型

delay_fit <- linear_reg() |>
  fit(as.numeric(collection_delay) ~ 
        customer_no + 
        sales_unit_price + 
        sales_extended+
        sales_count + 
        as.factor(sales_return) +
        as.factor(is_outlier), data = credit_sales_journal)
# 这个模型加入了之前考虑的一些变量

tidy(delay_fit) |> 
  mutate(across(where(is.numeric), ~ round(.x, 4)))

结果告诉我们，除了sales_return和is_outlier之外，其他变量和因变量的关系并不显著。这里，is_outlier本身就是y是否比较大，其实是一种机械相关，我们不能用y来预测y。

delay_fit <- linear_reg() |>
  fit(as.numeric(collection_delay) ~ 
        customer_no + 
        sales_unit_price + 
        sales_extended+
        sales_count + 
        as.factor(sales_return), data = credit_sales_journal)

tidy(delay_fit) |> 
  mutate(across(where(is.numeric), ~ round(.x, 4)))

去掉is_outlier之后，解释性依然很差。

# 用augment函数提取出残差和预测值
delay_aug <- augment(delay_fit, new_data = credit_sales_journal)

#画出残差图
ggplot(delay_aug, aes(x = collection_delay, y = .resid)) + 
  geom_point(aes(color = sales_return))

可以看到，很明显，residual，他也就是预测误差会随着collection_delay的上升而上升，这里大致可以说明，我们拟合的模型，并不能很好的预测collection delay，我们需要一个更好的模型来预测collection delay。

我们可以和之前预测frieght的模型进行对比：

vendor_invoices_dec <-  read_csv('case/InvoicePurchases12312016.csv')|> 
  janitor::clean_names()

freight_fit <- linear_reg() |>
  fit(freight ~ dollars + quantity + approval + approval*dollars + approval * quantity, data = vendor_invoices_dec)

tidy(freight_fit) |> 
  mutate(across(where(is.numeric), ~ round(.x, 4)))

# 用augment函数提取出残差和预测值
freight_aug <- augment(freight_fit, new_data = vendor_invoices_dec)

#画出残差图
ggplot(freight_aug, aes(x = freight, y = .resid)) + 
  geom_point(aes(color = approval))

可以发现，除了有一小簇数据有一个直线向上的趋势，其他大部分残差是在0周围分布的。

回归模型的诊断

接下来，我们介绍如何对回归模型进行诊断。这里，我们开始用baseR来拟合回归函数，因为tidymodels主要功能是拟合并预测，对于诊断的支持并不好。

之前的大部分内容中我们都是通过OLS回归法通过一系列的自变量来预测因变量。

多元回归模型的数学表达形式为：

\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots +\beta_p x_p + \epsilon\] 其中：

\(y\) 是因变量（被预测的变量）。 \(\beta_0\) 是截距（常数项）。

\((\beta_1, \beta_2, \dots, \beta_p)\) 是自变量（预测变量）对应的回归系数。

\(( x_1, x_2, \dots, x_p)\) 是自变量。

\(\epsilon\) 是误差项，表示模型无法解释的随机误差。

如果用矩阵形式表示：

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} \]

其中： \[ \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}\ \]

是因变量的向量。

\(\mathbf{X} = \begin{bmatrix} 1 & x_{11} & x_{12} & \dots & x_{1p} \\ 1 & x_{21} & x_{22} & \dots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \dots & x_{np} \end{bmatrix}\) 是包含自变量的矩阵。

\(\boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_p \end{bmatrix}\)是回归系数的向量。 \(\boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix}\)是误差项的向量。

最小二乘法的目标就是通过最小化误差平方和（RSS, Residual Sum of Squares）来求解回归系数：

\[ RSS = (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^\top (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \]

求导并令其等于 0：

\[ \frac{\partial RSS}{\partial \boldsymbol{\beta}} = -2 \mathbf{X}^\top (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) = 0 \]

化简后得到：

\[ \mathbf{X}^\top \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y} \]

这是一组线性方程组，称为 正则方程。

如果\(\mathbf{X}^\top \mathbf{X}\) 是非奇异矩阵（可逆），则回归系数的解为：

\[ \boldsymbol{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \]

为了能够恰当的解释OLS回归模型的系数，数据必须满足以下统计假设

正态性：因变量成正态分布

独立性：\(Y_i\) 值之间相互独立

线性：因变量与自变量直接线性相关

同方差性：因变量的方差不随自变量的取值不同二变化。

如果违背以上假设，那么统计显著性检验结果和执行区间很可能不精确。

刚才的练习中，我们可以用lm函数（基于baseR）建立一个多元回归模型,再用summary函数查看结果。

delay_fit <- lm(as.numeric(collection_delay) ~ 
            customer_no + 
            sales_unit_price + 
            sales_extended+
            sales_count + 
            as.factor(sales_return) , data = credit_sales_journal)
summary(delay_fit)

相比于tidymodels的输出，用标准包里的summary函数可以看到更多诊断信息，比如说F-statistics， R-squared等。

标准评估方法

R的基础包中提供了大量检验回归分析中统计假设的方法。最常见的就是对lm函数返回的对象使用plot函数。

library(car)
par(mfrow = c(2,2))
plot(delay_fit)
par(mfrow = c(1,1))

为了理解这些图形，我们来回顾一下OLS回归的统计假设：

1. 正态性：当自变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。“正态Q-Q图”(Normal Q-Q，右上)是在正态分布对应的值下准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上，若不是如此，那么就违反了正态性的假设。
2. 独立性：我们无法从这些图中分辨出因变量值是否相互独立，只能从数据的收集方式来验证。上面的例子中，我们需要判断一个应收款的延迟，是否会影响到另一个应收款的延迟。
3. 线性：若因变量与自变量线性相关，那么残差值与预测(拟合)值就没何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“图一拟合图”(Residualsvs Fitted，左上)中并不能看出残差和拟合值有很明显的关系。
4. 同方差性：若满足同方差性假设，那么在“位置-尺度关系图”(Scale- Location，左下)中，水平线周围的点应该随机分布。该图似乎满足此假设。

最后一幅“残差-杠杆关系图”(Residuals vs Leverage，右下)提供了我们可能关注的单个观测点的信息。从图形可以识别出离群点、高杠杆值点和强影响点。

• 一个观测点是离群点，表明拟合回归模型对其预测效果不佳(产生了巨大的或正或负的残差)。

• 一个观测点有很高的杠杆值，表明它是一个异常的自变量值的组合。也就是说，在自变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。

• 一个观测点是强影响点(influential observation)，表明它对模型参数的估计产生的影响过大，非常不成比例。强影响点可以通过Cook距离即Cook’sD统计量来识别。

改进的方法

1. 正态性

car包中的qqPlot函数提供了更精确的正态假设检验方法。

library(car)

qqPlot(delay_fit,id = list(method ='identity'),simulate = TRUE, main = 'Q-Q Plot')

2. 误差独立性

car包中提供了一个可以做Durbin-Watson检验的函数，能够检验误差的序列相关性。

durbinWatsonTest(delay_fit)

p值略微显著(0.092)，说明误差直接并不完全独立，lag = 1 表示数据集每个数据都是和其后面一个数据进行比较的。

3.线性

通过成分残差图，或者叫偏残差图，可以看看因变量和自变量之间是否存在非线性关系。

crPlots(delay_fit)

4. 同方差性

car包还提供两个有用的函数，可以判断误差方差是否恒定不变

ncvTest(delay_fit)

ncvTest函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平变化而变化。

这里p < 0.01，说明存在异方差性。

spreadLevelPlot(delay_fit)

这里可以看到，建议幂次转换（power transformation）接近0.5，也就是如果用\(\sqrt y\) 代替\(y\) 可能能够解决异方差性。

delay_fit2 <- lm(sqrt(as.numeric(collection_delay)) ~ 
            customer_no + 
            sales_unit_price + 
            sales_extended+
            sales_count + 
            as.factor(sales_return) , data = credit_sales_journal)
summary(delay_fit2)

ncvTest(delay_fit2)

可以看到，转换之后ncvtest不再显著，说明异方差性确实得到了缓解。

多重共线性

当回归结果都不显著时，有一种可能性是自变量直接相关性太高，导致了模型参数的置信区间过大，从而形成多重共线性问题。

多重共线性可以用统计量VIF（variance inflation factor，方差膨胀因子）来检验。一般原则下vif大于10就表明可能存在多重共线性问题。

vif(delay_fit2)

可以看到sales_unit_price，sales_extended的gvif都大于10，说明存在多重共线性问题。仔细思考一下，单价、总价和数量这三个变量一定是高度相关的，因此，可以尝试在模型中删除。

delay_fit3 <- lm(sqrt(as.numeric(collection_delay)) ~ 
            customer_no + 
            sales_count+
            as.factor(sales_return) , data = credit_sales_journal)
summary(delay_fit3)
vif(delay_fit3)

异常观测值

离群点、高杠杆值点和强影响点

离群点指的是模型效果预测不佳的观测点，通常有很大的残差。

高杠杆值观测点，即由许多异常自变量值组合起来的离群点，和因变量没有关系。

强影响点，对模型参数估计值的影响有些失衡的点，例如，如果移除某一个观察点，会对整个模型发出巨大的改变，那么我们就需要检验一下。

car包的influencePlot可以整合这三种异常观测点。

influencePlot(delay_fit3)

纵轴是residuals，residuals大于2或小于-2可以被认为是离群点，这里是270 和 630

横轴是hat-values，纵向的虚线标注了超过hat-values均值2倍和3倍，超过这两条线很多的会被标记出来为高杠杆值点，hat-value的均值为p/n，其中p是模型估计的参数数目（包括截距项），n是样本量。超过hat均值2倍或者3倍，可以判定为高杠杆值点。

颜色是cook’s D统计量，用来识别强影响点，这里颜色越深表示影响越大，一般来说，我们可以用D = 1 作为判别标准。

因此，我们可以先计算这些参考值，在考虑是否要删除这些离群点。

p <-  length(coefficients(delay_fit3))
n <- length(fitted(delay_fit3))
hat_mean <- p/n
print(hat_mean*2)
print(hat_mean*3)

# 提取影响点的行号
influential_points <- as.numeric(rownames(influencePlot(delay_fit3)))

# 打印这些行的invoice_no
credit_sales_journal[influential_points, "invoice_no"]

# 删除影响点
cleaned_credit <- credit_sales_journal[-influential_points, ]

delay_fit4 <- lm(sqrt(as.numeric(collection_delay)) ~ 
            customer_no + 
            sales_count+
            as.factor(sales_return) , data = cleaned_credit)
summary(delay_fit4)

删除了这些异常观测值之后，可以看到，sales_return的影响开始变得显著了，模型的拟合得到了一定的提升。

找到“最佳”的模型

baseR中的anova函数可以比较两个嵌套模型的拟合优度，所谓嵌套模型，就是回归方程的所有项完全包含在另一个模型中。

delay_fit4 <- lm(sqrt(as.numeric(collection_delay)) ~ 
            customer_no + 
            sales_count+
            as.factor(sales_return) , data = cleaned_credit)

delay_fit4_reduced <- lm(sqrt(as.numeric(collection_delay)) ~ 
            sales_count+
            as.factor(sales_return) , data = cleaned_credit)

anova(delay_fit4,delay_fit4_reduced)

p = 0.4804，可以得出结论，reduced的模型和原模型没有差别，也就不需要添加customer_no进入模型。

另一种方式是比较AIC，AIC越小的模型要有限选择。用AIC比较时，没有嵌套模型的限制。

AIC(delay_fit4, delay_fit4_reduced)

广义线性模型

logistic回归

把collection_delay转为二值型因子

credit_sales_journal <- credit_sales_journal |>
  mutate(
    collection_delay = as.numeric(collection_date - invoice_date),
    long_ar = factor(if_else(collection_delay > 200, 1, 0),
                     levels = c(0,1),
                     labels = c('no','yes'))
  )

table(credit_sales_journal$long_ar)

这个因子变量作为因变量

fit <-glm(long_ar ~ customer_no + sales_count + sales_unit_price + sales_extended
          + as.factor(sales_return), data = credit_sales_journal,
          family = binomial())
summary(fit)

fit_reduced <-glm(long_ar ~  sales_extended
          , data = credit_sales_journal,
          family = binomial())
summary(fit_reduced)

尝试缩小判定。

credit_sales_journal <- credit_sales_journal |>
  mutate(
    collection_delay = as.numeric(collection_date - invoice_date),
    long_ar = factor(if_else(collection_delay > 100, 1, 0),
                     levels = c(0,1),
                     labels = c('no','yes'))
  )

table(credit_sales_journal$long_ar)

fit <-glm(long_ar ~ customer_no + sales_count + sales_unit_price + sales_extended
          + as.factor(sales_return), data = credit_sales_journal,
          family = binomial())
summary(fit)

fit_reduced <-glm(long_ar ~ as.factor(sales_return), data = credit_sales_journal,
          family = binomial())
summary(fit_reduced)

可以发现，sales_return的影响是显著的。

anova(fit,fit_reduced)
AIC(fit,fit_reduced)

poisson回归

fit <- glm(collection_delay ~ customer_no + sales_count + sales_unit_price + sales_extended
           + as.factor(sales_return), data = credit_sales_journal,
           family = poisson())
summary(fit)

fit_reduced <- glm(collection_delay ~ customer_no + sales_unit_price
           + as.factor(sales_return), data = credit_sales_journal,
           family = poisson())
summary(fit_reduced)

判断是否存在过度离趋

deviance(fit)/df.residual(fit)
deviance(fit_reduced)/df.residual(fit_reduced)

远大于1，说明有可能高估了自变量的影响。

解决方法：

fit_od<- glm(collection_delay ~ customer_no + sales_unit_price
                   + as.factor(sales_return), data = credit_sales_journal,
                   family = quasipoisson())
summary(fit_od)