跳至主要内容

使用 VIM 包可视化缺失数据

学习使用 VIM 包提供的数据可视化工具，快速洞察缺失数据模式。
· 13分钟 阅读

使用 AI 探索

缺失数据是每位数据科学家日常工作中的难题。我们该填补它们吗？如果要填补，哪种方法合适？或者，包含缺失值的观测能否直接删除？要回答这些问题，首先需要了解缺失数据背后的机制。用统计检验来识别机制既复杂又常常只能得到模糊的结论。相比之下，可视化工具易于使用，不仅有助于识别缺失机制，还能帮助了解数据质量的其他方面。在本教程中，我们将介绍 VIM 包中的一组绘图方法，展示它们如何帮助您扎实把握数据缺失的模式。

缺失数据机制

数据集不完整的原因可能有很多。调查可能的缺失原因至关重要，因为这会影响我们处理问题的方式。例如，如果关注的是与工作相关的问卷调查的不响应现象，我们可能会预期最富有和最贫穷的受访者不愿在问卷中披露收入，这意味着缺失值并不是在数据集中均匀分布的。在这种情况下，在进行后续分析前简单删除不完整观测会导致结果偏倚。另一方面，如果数据缺失是由于采集设备故障造成，那么数据集中缺失点的位置也许是完全随机的。

通常有三种不同的模式可以解释数据为何缺失，这些模式被称为缺失数据机制。

  • 完全随机缺失（MCAR）
    在 MCAR 下，数据集中缺失点的位置没有任何系统性模式：它们完全随机出现。形式化地说，某条观测缺失的概率既不取决于其他变量的取值，也不取决于其自身的取值。在这种情况下，删除不完整观测不会给后续分析引入偏倚。
    示例：传感器持续测量温度，将数据通过互联网发送到数据库。由于未知原因，网络连接有时会中断。

  • 随机缺失（MAR）
    在 MAR 下，某条观测缺失的概率仍然不依赖于其自身的取值，但会依赖于其他变量的取值。在这种情况下，删除不完整观测会使样本代表性降低。
    示例：部分夜间时段的数据因传感器维护而缺失，而维护总是在夜间进行。

  • 非随机缺失（MNAR）
    在 MNAR 下，某条观测缺失的概率取决于它自身未被观测到的取值。在这种情况下，同样，删除不完整数据会导致有偏分析。
    示例：传感器在 -20 摄氏度时会冻结，无法测量低于该值的温度。

在实践中，很难判断具体情形对应哪一种缺失机制。为说明这一点，让我们绘制一些 MNAR 数据。我们将使用 ggplot2 和 gridExtra 包，需事先加载。同时还加载另两个稍后会用到的包：用于缺失数据可视化分析的 VIM，以及用于数据预处理的 dplyr。如果这些包尚未安装，请先使用 install.packages() 安装。

library("VIM")
library("dplyr")
library("ggplot2")
library("gridExtra")

我们将创建一个包含两个不相关、服从正态分布的变量 x 和 y 的数据集，在 y 的分布右尾（即最高值部分）设置一些缺失值。注意我们也保留了这些位置的真实值，并在单独的变量中将其标记为缺失。

set.seed(2)
mnar_data <- data.frame(x = rnorm(100), y = rnorm(100)) %>%
    mutate(y_miss = ifelse(y > 1, y, NA),
           y = ifelse(is.na(y_miss), y, NA),
           x_miss = ifelse(is.na(y), x, NA))

现在我们绘制两个 x 对 y 的散点图：一个包含那些我们假定 y 值缺失的点，另一个只包含已观测数据。

grid.arrange(
  # Plot highlighting missing data points
  ggplot(mnar_data, aes(x, y)) +
    geom_point(size = 4, alpha = 0.6) +
    geom_point(aes(x, y_miss), col = "red", size = 4, alpha = 0.6) +
    geom_hline(aes(yintercept = y_miss), col = "red", alpha = 0.6,
               linetype = "dashed") +
    ylim(NA, max(mnar_data$y_miss, na.rm = TRUE)) +
    ggtitle("Where the data points are missing"),
  # What data scientist can see
  ggplot(mnar_data, aes(x, y)) +
    geom_point(size = 4, alpha = 0.6) +
    geom_vline(aes(xintercept = x_miss), col = "red", alpha = 0.6,
               linetype = "dashed") +
    ylim(NA, max(mnar_data$y_miss, na.rm = TRUE)) +
    ggtitle("What the data scientist can see"),
  # Arrange the two plots side-by-side
  ncol = 2
)
scatter plots

从第一个图可以清楚地看出数据确实是 MNAR，但数据科学家所能看到的只有缺失点的 x 坐标。第二个图能得出什么结论呢？数据是 MCAR 吗？也许是，缺失值在 x 轴上的位置似乎分布得相当均匀。或者是 MAR？仔细观察，x 为负时 y 的缺失更多于 x 为正时，因此也可能是这种情况。最后，MNAR 永远无法被排除，因为要判断数据是否 MNAR，从定义上说我们需要观察那些未被观测到的取值。

幸运的是，我们并非完全无从下手。数据可视化工具可以为我们提供一些指引，帮助识别缺失数据中的模式。接下来我们将基于 fivethirtyeight 包中的 biopics 数据集（包含一些传记片的信息）分析若干用于检测这些模式的图。

数据清洗

在进入绘图部分之前，我们需要对数据做一些预处理。下面的代码块将完成这项工作。dplyr 管道开头的 select 提取感兴趣的变量，随后用 mutate 将相关种族合并，以确保每组都有足够的观测。接着调整部分变量的类型：字符串变量需要转换为因子，逻辑变量转换为整数。最后，为提高绘图的可读性，将变量重命名为更短的名称。

data(biopics, package = "fivethirtyeight")

biopics <- biopics %>%
  select(country, year_release, box_office, number_of_subjects,
         type_of_subject, subject_race, person_of_color, subject_sex) %>%
  mutate(subject_race = ifelse(grepl("^Hispanic", subject_race), "Hispanic",
                               subject_race),
         subject_race = ifelse(grepl("^African", subject_race), "African",
                               subject_race),
         subject_race = ifelse(grepl("^Middle", subject_race), "Mid Eastern",
                               subject_race),
         subject_race = ifelse(subject_race %in% c("White", "Asian", "African",
                                                   "Hispanic", "Mid Eastern",
                                                   "Multi racial", NA),
                               subject_race, "other")) %>%
  mutate(country = as.factor(country),
         type_of_subject = as.factor(type_of_subject),
         subject_race = as.factor(subject_race),
         subject_sex = as.factor(subject_sex),
         person_of_color = as.integer(person_of_color)) %>%
      as.data.frame()

colnames(biopics) <- c("country", "year", "earnings", "sub_num",
                       "sub_type", "sub_race", "non_white", "sub_sex")

最终我们得到包含八个变量的数据集：

  • country - 影片的制作国家或地区，
  • year - 上映年份，
  • earnings - 美国票房总收入，
  • sub_num - 影片中出现的人物数量，
  • sub_type - 人物的职业或成名原因，
  • sub_race - 人物的种族，
  • non_white - 表示有色人种的虚拟变量，
  • sub_sex - 人物性别。

聚合图

首先要问的基本问题是：哪些变量有缺失，缺失了多少？聚合图是回答这些问题的有用工具。下面这一行代码就够了。

aggr(biopics, numbers = TRUE, prop = c(TRUE, FALSE))
aggregation plots

我们指定参数 numbers = TRUE 以在柱状图顶部显示数值。参数 prop 为逻辑值，指示是否显示缺失值及其组合的比例而非总量。我们在第一个图中设为 TRUE，在第二个图中设为 FALSE。

很明显，只有两个变量存在缺失：earnings 的缺失超过 40%，sub_race 约为 25%。在右侧的组合图中，网格展示了数据中所有缺失（红）与观测（蓝）取值的组合。共有 317 条完整观测，其中有 77 行这两个变量同时缺失。

同时看两个图可能略显混乱，因此我们用下一行代码将同样的信息压缩到一个可视化中。

aggr(biopics, combined = TRUE, numbers = TRUE)

combined 参数允许将两个图合并为一个。此处，网格右侧的水平条显示相应组合的频数，而顶部的垂直条展示各变量的缺失比例。将 earnings 的所有组合值相加，即 0.32 和 0.10，可见我们先前估计为 40% 以上的总体缺失比例实际上是 42%。

旋转直方图与脊柱图

到目前为止，我们对缺失数据有了宏观认识。现在是时候更深入地观察具体变量之间的相互作用了。旋转直方图（spinogram） 和 脊柱图（spineplot） 用于研究一个变量在另一个变量不同取值下的缺失比例。如果后者是数值型变量，则使用旋转直方图；若是分类型变量，则称为脊柱图。

两者都可以使用 spineMiss() 函数生成，输入为两列的数据框。第一个指定的变量用于划分数据，映射到图的水平轴。第二个变量是我们关注其缺失模式的变量。来看一个脊柱图的示例。我们先指定分类变量 sub_race，然后是 earnings，也就是：在 sub_race 的各类别下，earnings 的缺失比例是多少？

spineMiss(biopics[, c("sub_race", "earnings")])
Spinogram and Spineplot

sub_race 各类别条形的相对宽度反映了该类别在数据集中的频率：例如，绝大多数影片的主角是白人。每个条形内部展示了 earnings 的缺失比例，而右侧阴影条展示了整个数据集的该比例。看起来当主角为非洲裔时，我们最有可能获得完整的票房信息。

通过像下面这样交换变量顺序，我们可以生成一个旋转直方图，回答另一个方向的问题：在 earnings 的不同取值下，sub_race 的缺失比例是多少？

spineMiss(biopics[, c("earnings", "sub_race")])
Spinogram and Spineplot 2

由于 earnings 是数值型，其取值被分箱，箱宽对应变量自身的分布。该旋转直方图显示 earnings 高度左偏：只有少数影片获得了最高利润。有趣的是，对于这些票房大片，我们最可能缺失主角种族信息，这体现在 earnings 最高值对应的红色条形最高。

马赛克图

上文的旋转直方图和脊柱图用于研究两个变量之间的交互。将这一思路推广到更多变量，就得到马赛克图。该图由许多小块组成，每个小块对应两个或多个变量的特定类别（对分类型变量）或分箱（对数值型变量）的组合。在每个小块内，展示另一个变量的缺失比例。原则上，马赛克图可用于任意数量的变量，但如果太多就会显得杂乱。并且，当用于划分的数据为因子且水平数不多时，图形更易读。

让我们看看按 sub_sex 和 US_movie 划分的 earnings 缺失比例。后者在下面的代码块中创建，是一个逻辑变量，表示美国是否参与了该片的制作。我们将包含三个变量的数据框作为 mosaicMiss() 的第一个参数。将 plotvars 设为包含 1 和 2 的向量，表示我们希望按照该数据框的前两列划分数据。将 highlight 设为 3，表示在小块中显示第三列变量的缺失比例。

biopics <- biopics %>%
    mutate(US_movie = ifelse(grepl("US", country), TRUE, FALSE))
mosaicMiss(biopics[, c("sub_sex", "US_movie", "earnings")], highlight = 3,
           plotvars = 1:2, miss.labels = FALSE)
Mosaic Plot

同样，小块的大小对应该组合在数据集中的出现频率。例如，右下角最大的方块意味着大多数电影的主角为男性，且至少部分由美国参与制作。

对于主角为男性的电影，非美国电影更可能缺失票房信息。另一方面，如果主角为女性，美国电影在票房上略多一些缺失。不过，这些差异似乎过小，难言显著。

并行箱线图

另一种可视化是并行箱线图。其思路是将数据集拆分为两个子集：一个只包含某个不完整变量的观测值，另一个只包含其缺失值。对这两个子集各画一个选定数值变量的箱线图。这样您可以检查该变量的分布是否受到拆分变量缺失情况的影响。

我们可以使用 pbox() 函数绘制并行箱线图。只需传入一个两列数据框。第一列的变量将用于绘制箱线图，第二列用于拆分。在下面的示例中，我们对高度偏斜的 earnings 取对数，否则箱线图会近似为一条线。

biopics <- biopics %>%
  mutate(log_earnings = log(earnings))
pbox(biopics[, c("log_earnings", "sub_race")])
Parallel Boxplot

左侧白色箱线展示了 log_earnings 的总体分布，蓝色和红色箱线分别展示了 sub_race 为已观测和缺失子集中的分布。箱线的相对宽度反映了对应子集的大小：蓝箱更宽意味着 sub_race 中观测值多于缺失值。除此之外，两个箱线彼此相似，也与整体白箱相似。这表明缺失的种族信息对票房分布没有影响。

平行坐标图

我们已经观察了单个变量及它们的交互。现在来一次性分析数据集中的所有变量。允许这样做的可视化是平行坐标图。在该图中，每个变量被转换到相同刻度，并由一条平行坐标轴表示。对于分类变量，坐标轴的刻度被分为等距的若干点，每个类别一个。缺失值被放在纵轴上方、图形区域之外。图上的每条线对应数据中的一条观测，线的颜色表示所选变量是否缺失。

我们可以使用 parcoordMiss() 函数生成平行坐标图。默认情况下，它使用作为第一个参数传入的数据框中的所有变量。将 highlight 设为 earnings 会用不同颜色标记那些 earnings 缺失的观测对应的线。将 alpha 混合设为 0.6 可使线条半透明，提升可读性。

# Remove variables created for the previous plots
biopics <- biopics %>%
  select(- US_movie, - log_earnings)

parcoordMiss(biopics, highlight = 'earnings', alpha = 0.6)
Parallel Coordinate Plot

在上图中，深红色线表示 earnings 缺失的观测。这些高亮观测似乎与其余数据的行为不同。尤其是，只有少量它们对应于 country 的次高水平，而该水平原本包含大量观测！结果发现这是美/英联合制作的电影（要验证，只需运行 levels(biopics$country) 并查看倒数第二个值）。

此外，快速查看 year 轴可以发现，在其大约三分之二高度处有一段区间，几乎没有深红色线穿过。似乎在不久前有一段时间票房数据更为完整。这表明 country 和 year 都可能有助于解释票房缺失值的分布！

矩阵图

我们将要看的最后一种可视化是矩阵图。它以矩形的形式展示数据矩阵的所有单元格。观测数据使用连续的灰—黑配色（颜色越深，数值越大），缺失值以红色高亮。一个好习惯是按照某个不完整变量对数据排序——这会让图更易解读。下面的代码块按 earnings 进行排序。

matrixplot(biopics, sortby = c('earnings'))
Matrix Plot

该图验证了之前的一些发现：earnings 缺失的观测往往来自较早时期（year 取值较小），且在 country 上也偏低——这体现在那些票房为红色的行里这两个变量较浅的颜色。此外，似乎当 earnings 缺失时，主角为非白人的电影更少。

结论

总而言之，biopics 数据集中的缺失数据似乎不是 MCAR。earnings 与 sub_race 中缺失点的位置可以相互解释，这已被脊柱图和旋转直方图证明。此外，矩阵图和平行坐标图表明，country、year 和 non_white 也可能有助于解释 earnings 的缺失分布。因此，不宜删除不完整观测，因为这很可能会给推断带来偏倚。

最后的想法

您已经走了很远！现在您知道了缺失数据机制是什么、它们有何不同，更重要的是，您也掌握了如何用多种可视化工具对其进行推断。干得漂亮！您现在已完全具备分析自己不完整数据集的能力！

如果您想进一步学习 R，请参加 DataCamp 的 Data Visualization with ggplot2 (Part 1) 课程，并查看我们的 R Formula Tutorial。

主题

数据可视化数据科学