机器学习的数据清理和探索-一-

机器学习的数据清理和探索（一）

原文：annas-archive.org/md5/25ad0fee8d118820a9d79ad1484952bd

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

研究人员为数据分析准备数据的工作——包括提取、转换、清洗和探索——在机器学习工具日益普及的情况下，其本质并没有发生根本变化。30 年前，当我们为多元分析准备数据时，我们同样关注缺失值、异常值、变量分布的形状以及变量之间的相关性，就像我们现在使用机器学习算法时一样。虽然确实，广泛使用相同的机器学习库（如 scikit-learn、TensorFlow、PyTorch 等）确实鼓励了方法上的更大一致性，但良好的数据清洗和探索实践在很大程度上并未改变。

我们谈论机器学习的方式仍然非常侧重于算法；只需选择正确的模型，随之而来的就是组织变革的洞见。但我们必须为过去几十年中我们一直在进行的数据学习留出空间，其中我们从数据中做出的预测、在数据中建模关系以及我们对数据的清洗和探索都是对话的一部分。确保我们的模型正确与从直方图或混淆矩阵中获取尽可能多的信息一样重要，需要仔细调整超参数。

同样，数据分析师和科学家的工作并不从清洗、探索、预处理、建模到评估这一过程有序推进。我们在过程的每一步都怀有潜在模型的想法，并定期更新我们之前的模型。例如，我们最初可能认为我们将使用逻辑回归来建模特定的二元目标，但当我们看到特征的分布时，我们可能至少需要尝试使用随机森林分类。我们将在整篇文章中讨论建模的影响，即使在解释相对常规的数据清洗任务时也是如此。我们还将探讨在早期过程中使用机器学习工具，以帮助我们识别异常、插补值和选择特征。

这指向了数据分析师和科学家在过去十年中工作流程的另一个变化——对“单一模型”的重视减少，以及对模型构建作为迭代过程的接受度提高。一个项目可能需要多个机器学习算法——例如，主成分分析来降低维度（特征的数量），然后是逻辑回归进行分类。

话虽如此，我们在数据清洗、探索和建模方面的方法有一个关键的区别——随着机器学习工具在我们的工作中扮演越来越重要的角色，我们对预测的重视超过了对底层数据的理解。我们更关心我们的特征（也称为自变量、输入或预测因子）如何预测我们的目标（因变量、输出、响应），而不是特征之间的关系以及我们数据的底层结构。我在本书的前两章中指出了这一点如何改变我们的关注点，即使我们在清洗和探索数据时也是如此。

本书面向的对象

在撰写这本书时，我考虑了多个受众，但我最经常想到的是我的一个好朋友，30 年前她买了一本 Transact-SQL 书，立刻对自己的数据库工作有了极大的信心，最终围绕这些技能建立起了自己的职业生涯。如果有人刚开始作为数据科学家或分析师的职业生涯，通过这本书并获得了与我朋友相似的经历，我将感到非常高兴。最重要的是，我希望你通过阅读这本书后感到满意和兴奋，对你可以做到的事情感到自豪。

我还希望这本书对那些已经从事这类工作一段时间的人来说是一本有用的参考书。在这里，我想象有人打开这本书，自己问自己，在我的逻辑回归模型网格搜索中，应该使用哪些好的值？

为了保持本书的实践性质，书中的每一部分输出都可以通过本书中的代码进行重现。我始终坚持一个规则，即使有时会遇到挑战。除了概念性章节外，每个章节都是从原始下载文件中基本未变的数据开始的。你将在每个章节中从数据文件到模型进行转换。如果你忘记了某个特定对象是如何创建的，你只需要翻回一页或两页就能看到。

对于那些对 pandas 和 NumPy 有一定了解的读者来说，在处理一些代码块时会更加得心应手，同样，对 Python 和 scikit-learn 有一定了解的人也是如此。尽管如此，这些都不是必需的。有些部分你可能需要花更多的时间去仔细阅读。如果你需要关于使用 Python 进行数据工作的额外指导，我认为我的 Python 数据清洗食谱 是一本很好的配套书籍。

本书涵盖的内容

第一章，检查特征和目标的分布，探讨了使用常见的 NumPy 和 pandas 技术来更好地了解数据的属性。我们将生成汇总统计量，如均值、最小值和最大值，以及标准差，并计算缺失值的数量。我们还将创建关键特征的可视化，包括直方图和箱线图，以比仅查看汇总统计量更好地了解每个特征的分布。我们将暗示特征分布对数据转换、编码和缩放以及我们在后续章节中用相同数据进行建模的影响。

第二章，检查特征和目标之间的双变量和多变量关系，专注于可能特征和目标变量之间的相关性。我们将使用 pandas 方法进行双变量分析，并使用 Matplotlib 进行可视化。我们将讨论我们发现的特征工程和建模的影响。我们还在本章中使用多元技术来理解特征之间的关系。

第三章，识别和修复缺失值，将介绍识别每个特征或目标缺失值的技术，以及识别大量特征值缺失的观测值。我们将探讨填充值的方法，例如将值设置为整体均值、给定类别的均值或前向填充。我们还将检查用于填充缺失值的多元技术，并讨论它们何时适用。

第四章，编码、转换和缩放特征，涵盖了各种特征工程技术。我们将使用工具删除冗余或高度相关的特征。我们将探索最常见的编码类型——独热编码、有序编码和哈希编码。我们还将使用转换来改善特征的分布。最后，我们将使用常见的分箱和缩放方法来解决偏斜、峰度和异常值，以及调整范围差异较大的特征。

第五章，特征选择将介绍多种特征选择方法，从过滤器到包装器，再到嵌入式方法。我们将探讨它们如何与分类和连续目标一起工作。对于包装器和嵌入式方法，我们将考虑它们与不同算法结合时的效果。

第六章，为模型评估做准备，将展示我们构建第一个完整的流水线，将数据分为测试集和训练集，并学习如何在没有数据泄露的情况下进行预处理。我们将实现 k 折交叉验证，并更深入地研究评估模型性能的方法。

第七章，线性回归模型，是关于使用许多数据科学家喜爱的老方法——线性回归来构建回归模型的几个章节中的第一个。我们将运行一个经典的线性模型，同时考察使特征空间成为线性模型良好候选者的特性。我们将探讨在必要时如何通过正则化和变换来改进线性模型。我们将研究随机梯度下降作为普通最小二乘法（OLS）优化的替代方案。我们还将学习如何使用网格搜索进行超参数调整。

第八章，支持向量回归，讨论了关键的支持向量机概念以及它们如何应用于回归问题。特别是，我们将考察诸如 epsilon 敏感管和软边界等概念如何为我们提供灵活性，以在数据和领域相关挑战下获得最佳拟合。我们还将首次但肯定不是最后一次探索非常实用的核技巧，它允许我们在不进行变换或增加特征数量的情况下建模非线性关系。

第九章，K 近邻、决策树、随机森林和梯度提升回归，探讨了最流行的非参数回归算法中的一些。我们将讨论每个算法的优点，何时可能想要选择一个而不是另一个，以及可能的建模挑战。这些挑战包括如何通过仔细调整超参数来避免欠拟合和过拟合。

第十章，逻辑回归，是关于使用逻辑回归构建分类模型的几个章节中的第一个，逻辑回归是一种高效且偏差低的算法。我们将仔细检查逻辑回归的假设，并讨论数据集和建模问题中使逻辑回归成为良好选择的属性。我们将使用正则化来解决高方差问题或当我们有许多高度相关的预测因子时。我们将通过多项式逻辑回归将算法扩展到多类问题。我们还将讨论如何首次但肯定不是最后一次处理类别不平衡问题。

第十一章, 决策树和随机森林分类，回到在第九章中介绍的决策树和随机森林算法，这次处理分类问题。这为我们提供了另一个学习如何构建和解释决策树的机会。我们将调整包括树深度在内的关键超参数，以避免过拟合。然后我们将探索随机森林和梯度提升决策树作为决策树的优秀、低方差替代方案。

第十二章, 用于分类的 K 近邻，回到k 近邻（KNNs）来处理二元和多类建模问题。我们将讨论并展示 KNN 的优势——构建无装饰模型的简便性以及需要调整的超参数数量有限。到本章结束时，我们将了解两个问题——如何进行 KNN 以及何时应该考虑它来进行我们的建模。

第十三章, 支持向量机分类，探讨了实现支持向量分类（SVC）的不同策略。我们将使用线性 SVC，当我们的类别是线性可分时，它可以表现得非常好。然后我们将考察如何使用核技巧将 SVC 扩展到类别不是线性可分的情况。最后，我们将使用一对一和一对多分类来处理具有两个以上值的标签。

第十四章, 朴素贝叶斯分类，本章讨论了朴素贝叶斯的基本假设以及该算法如何被用来解决我们已探讨的一些建模挑战，以及一些新的挑战，例如文本分类。我们将考虑何时朴素贝叶斯是一个好的选择，何时则不是。我们还将探讨朴素贝叶斯模型的解释。

第十五章, 主成分分析，考察主成分分析（PCA），包括其工作原理以及我们可能想要使用它的时机。我们将学习如何解释 PCA 创建的成分，包括每个特征如何贡献到每个成分以及解释了多少方差。我们将学习如何可视化成分以及如何在后续分析中使用成分。我们还将考察如何使用核函数进行 PCA 以及何时这可能给我们带来更好的结果。

第十六章，K-Means 和 DBSCAN 聚类，探讨了两种流行的聚类技术，k-means 和基于密度的空间聚类应用噪声（DBSCAN）。我们将讨论每种方法的优点，并培养在何时选择一种聚类算法而不是另一种算法的感觉。我们还将学习如何评估我们的聚类以及如何更改超参数以改进我们的模型。

要充分利用本书

要运行本书中的代码，您需要安装一个科学版的 Python，例如 Anaconda。所有代码都使用 scikit-learn 版本 0.24.2 和 1.0.2 进行了测试。

下载示例代码文件

您可以从 GitHub（github.com/PacktPublishing/Data-Cleaning-and-Exploration-with-Machine-Learning）下载本书的示例代码文件。如果代码有更新，它将在 GitHub 仓库中更新。

我们还有其他来自我们丰富的图书和视频目录的代码包，可在github.com/PacktPublishing/找到。查看它们吧！

下载彩色图像

我们还提供了一份包含本书中使用的截图和图表的彩色图像的 PDF 文件。您可以从这里下载：packt.link/aLE6J。

使用的约定

本书使用了多种文本约定。

文本中的代码：表示文本中的代码词汇、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 昵称。以下是一个示例：“为了学习目的，我们在 GitHub 的chapter08文件夹下提供了两个示例mlruns工件和huggingface缓存文件夹。”

代码块设置如下：

client = boto3.client('sagemaker-runtime') 

response = client.invoke_endpoint(

        EndpointName=app_name, 

        ContentType=content_type,

        Accept=accept,

        Body=payload

        )

当我们希望将您的注意力引到代码块的一个特定部分时，相关的行或项目将以粗体显示：

loaded_model = mlflow.pyfunc.spark_udf(

    spark,

    model_uri=logged_model, 

    result_type=StringType())

任何命令行输入或输出都如下所示：

mlflow models serve -m models:/inference_pipeline_model/6

粗体：表示新术语、重要词汇或屏幕上看到的词汇。例如，菜单或对话框中的单词以粗体显示。以下是一个示例：“要执行此单元格中的代码，您只需在右上角的下拉菜单中点击运行单元格。”

小贴士或重要提示

看起来是这样的。

联系我们

欢迎读者反馈。

customercare@packtpub.com并在邮件主题中提及书名。

勘误表：尽管我们已经尽最大努力确保内容的准确性，但错误仍然可能发生。如果您在这本书中发现了错误，如果您能向我们报告，我们将不胜感激。请访问www.packtpub.com/support/err…并填写表格。

copyright@packt.com并附上材料链接。

如果您有兴趣成为作者：如果您在某个领域有专业知识，并且有兴趣撰写或为书籍做出贡献，请访问authors.packtpub.com。

分享您的想法

一旦您阅读了《使用机器学习进行数据清洗和探索》，我们非常想听听您的想法！请点击此处直接进入此书的亚马逊评论页面并分享您的反馈。

您的评论对我们和科技社区都非常重要，并将帮助我们确保我们提供高质量的内容。

第一部分 – 数据清洗和机器学习算法

我尽量避免按顺序思考模型构建过程的各个部分，把自己看作是清理数据，然后是预处理，等等，直到我完成模型验证。我不想把那个过程看作是包含永远结束的阶段。我们在这个节开始时从数据清洗开始，但我希望本节中的章节传达出我们总是在向前看，在清理数据时预测建模挑战；同时，在我们评估模型时，我们通常也会回顾我们所做的数据清洗。

在一定程度上，干净和脏的隐喻隐藏了为后续分析准备数据时的微妙之处。真正的关注点是我们的实例和属性（观察和变量）对感兴趣现象的代表性如何。这总是可以改进的，而且如果不加注意，很容易变得更糟。尽管如此，有一点是可以确定的。在模型构建过程的任何其他部分，我们都无法纠正数据清洗过程中犯下的重要错误。

本书的前三章是关于尽可能准确地获取我们的数据。要做到这一点，我们必须对所有的变量、特征和目标是如何分布的有一个良好的理解。在我们进行任何正式分析之前，我们应该问自己三个问题：1) 我们是否确信我们知道每个感兴趣变量的全部值域和分布形状？2) 我们是否对变量之间的双变量关系有一个良好的了解，即每个变量是如何与其他变量一起变化的？3) 我们尝试解决潜在问题（如异常值和缺失值）的成功率如何？本节中的章节提供了回答这些问题的工具。

本节包括以下章节：

• 第一章，检查特征和目标的分布

• 第二章，检查特征和目标之间的双变量和多变量关系

• 第三章，识别和修复缺失值

第一章：第一章：检查特征和目标的分布

机器学习写作和指导通常是算法导向的。有时，这给人一种印象，我们只需要选择正确的模型，组织变革的见解就会随之而来。但开始机器学习项目的最佳地方是对我们将使用的特征和目标分布的理解。

对于我们作为分析师几十年来一直重视的数据学习，留出空间是非常重要的——研究变量的分布、识别异常值以及检查双变量关系——即使我们越来越关注我们预测的准确性。

我们将在本书的前三章中探索用于此目的的工具，同时考虑对模型构建的影响。

在本章中，我们将使用常见的 NumPy 和 pandas 技术来更好地了解我们数据的属性。在我们进行任何预测分析之前，我们想知道关键特征的分布情况。我们还想知道每个连续特征的集中趋势、形状和分布范围，以及分类特征的每个值的计数。我们将利用非常方便的 NumPy 和 pandas 工具来生成汇总统计信息，例如均值、最小值和最大值，以及标准差。

之后，我们将创建关键特征的可视化，包括直方图和箱线图，以帮助我们更好地了解每个特征的分布，而不仅仅是通过查看汇总统计信息。我们将暗示特征分布对数据转换、编码和缩放以及我们在后续章节中用相同数据进行建模的影响。

具体来说，在本章中，我们将涵盖以下主题：

• 数据子集

• 为分类特征生成频率

• 为连续特征生成汇总统计信息

• 在单变量分析中识别极端值和异常值

• 使用直方图、箱线图和小提琴图来检查连续特征的分布

技术要求

本章将大量依赖 pandas、NumPy 和 Matplotlib 库，但你不需要对这些库有任何先前的知识。如果你从科学发行版，如 Anaconda 或 WinPython 安装了 Python，那么这些库可能已经安装好了。如果你需要安装其中之一来运行本章中的代码，你可以在终端中运行pip install [package name]。

数据子集

几乎我参与的每一个统计建模项目都需要从分析中移除一些数据。这通常是因为缺失值或异常值。有时，我们限制分析数据集的理论原因。例如，我们拥有从 1600 年以来的天气数据，但我们的分析目标仅涉及 1900 年以来的天气变化。幸运的是，pandas 中的子集工具非常强大且灵活。在本节中，我们将使用美国国家青年纵向调查（NLS）的数据。

注意

青年 NLS 是由美国劳工统计局进行的。这项调查始于 1997 年，当时出生在 1980 年至 1985 年之间的一批个人，每年通过 2017 年进行年度跟踪。对于这个配方，我从调查的数百个数据项中提取了关于成绩、就业、收入和对政府态度的 89 个变量。可以从存储库下载 SPSS、Stata 和 SAS 的单独文件。NLS 数据可在www.nlsinfo.org/investigator/pages/search公开使用。

让我们使用 pandas 开始子集数据：

1. 我们将首先加载 NLS 数据。我们还设置了一个索引：

   import pandas as pd
   import numpy as np
   nls97 = pd.read_csv("data/nls97.csv")
   nls97.set_index("personid", inplace=True)
   

2. 让我们从 NLS 数据中选择几个列。以下代码创建了一个新的 DataFrame，其中包含一些人口统计和就业数据。pandas 的一个有用特性是，新的 DataFrame 保留了旧 DataFrame 的索引，如下所示：

   democols = ['gender','birthyear','maritalstatus',
    'weeksworked16','wageincome','highestdegree']
   nls97demo = nls97[democols]
   nls97demo.index.name
   'personid'
   

3. 我们可以使用切片通过位置选择行。nls97demo[1000:1004]选择从冒号左侧的整数（在这种情况下是1000）指示的行开始，直到但不包括冒号右侧的整数（在这种情况下是1004）指示的行。1000行是第 1,001 行，因为索引是从 0 开始的。由于我们将结果 DataFrame 进行了转置，所以每一行都作为输出中的一个列出现：

   nls97demo[1000:1004].T
   personid      195884       195891        195970\
   gender        Male         Male          Female
   birthyear     1981         1980          1982
   maritalstatus NaN          Never-married Never-married
   weeksworked16 NaN          53            53
   wageincome    NaN          14,000        52,000   
   highestdegree 4.Bachelors  2.High School 4.Bachelors
   personid       195996  
   gender         Female  
   birthyear      1980  
   maritalstatus  NaN  
   weeksworked16  NaN  
   wageincome     NaN
   highestdegree  3.Associates    
   

4. 我们还可以通过在第二个冒号之后设置步长值来跳过区间内的行。步长的默认值是 1。以下步长的值是 2，这意味着在1000和1004之间的每隔一行将被选中：

   nls97demo[1000:1004:2].T
   personid        195884       195970
   gender          Male         Female
   birthyear       1981         1982
   maritalstatus   NaN          Never-married
   weeksworked16   NaN          53
   wageincome      NaN          52,000
   highestdegree   4.Bachelors  4\. Bachelors
   

5. 如果我们在冒号左侧不包含值，行选择将从第一行开始。请注意，这返回的 DataFrame 与head方法返回的相同：

   nls97demo[:3].T
   personid       100061         100139          100284
   gender         Female         Male            Male
   birthyear      1980           1983            1984
   maritalstatus  Married        Married         Never-married
   weeksworked16  48             53              47
   wageincome     12,500         120,000         58,000
   highestdegree  2.High School  2\. High School  0.None
   nls97demo.head(3).T
   personid       100061         100139         100284
   gender         Female         Male           Male
   birthyear      1980           1983           1984
   maritalstatus  Married        Married        Never-married
   weeksworked16  48             53             47
   wageincome     12,500         120,000        58,000
   highestdegree  2.High School  2.High School  0\. None
   

6. 如果我们在冒号左侧使用负数-n，将返回 DataFrame 的最后 n 行。这返回的 DataFrame 与tail方法返回的相同：

    nls97demo[-3:].T
   personid       999543          999698        999963
   gender         Female         Female         Female
   birthyear      1984           1983           1982
   maritalstatus  Divorced       Never-married  Married
   weeksworked16  0              0              53
   wageincome     NaN            NaN            50,000
   highestdegree  2.High School  2.High School  4\. Bachelors
    nls97demo.tail(3).T
   personid       999543         999698         999963
   gender         Female         Female         Female
   birthyear      1984           1983           1982
   maritalstatus  Divorced       Never-married  Married
   weeksworked16  0              0              53
   wageincome     NaN            NaN            50,000
   highestdegree  2.High School  2.High School  4\. Bachelors
   

7. 我们可以使用loc访问器通过索引值选择行。回想一下，对于nls97demo DataFrame，索引是personid。我们可以向loc访问器传递一个索引标签列表，例如loc[[195884,195891,195970]]，以获取与这些标签相关的行。我们还可以传递索引标签的下限和上限，例如loc[195884:195970]，以检索指定的行：

    nls97demo.loc[[195884,195891,195970]].T
   personid       195884       195891         195970
   gender         Male         Male           Female
   birthyear      1981         1980           1982
   maritalstatus  NaN          Never-married  Never-married
   weeksworked16  NaN          53             53
   wageincome     NaN          14,000         52,000
   highestdegree  4.Bachelors  2.High School  4.Bachelors
    nls97demo.loc[195884:195970].T
   personid       195884       195891         195970
   gender         Male         Male           Female
   birthyear      1981         1980           1982
   maritalstatus  NaN          Never-married  Never-married
   weeksworked16  NaN          53             53
   wageincome     NaN          14,000         52,000
   highestdegree  4.Bachelors  2.High School  4.Bachelors
   

8. 要按位置选择行，而不是按索引标签，我们可以使用iloc访问器。我们可以传递一个位置数字列表，例如iloc[[0,1,2]]，到访问器以获取那些位置的行。我们可以传递一个范围，例如iloc[0:3]，以获取介于下限和上限之间的行，不包括上限所在的行。我们还可以使用iloc访问器来选择最后 n 行。iloc[-3:]选择最后三行：

    nls97demo.iloc[[0,1,2]].T
   personid       100061         100139         100284
   gender         Female         Male           Male
   birthyear      1980           1983           1984
   maritalstatus  Married        Married        Never-married
   weeksworked16  48             53             47
   wageincome     12,500         120,000        58,000
   highestdegree  2.High School  2.High School  0\. None
    nls97demo.iloc[0:3].T
   personid       100061         100139         100284
   gender         Female         Male           Male
   birthyear      1980           1983           1984
   maritalstatus  Married        Married        Never-married
   weeksworked16  48             53             47
   wageincome     12,500         120,000        58,000
   highestdegree  2.High School  2.High School  0\. None
    nls97demo.iloc[-3:].T
   personid       999543         999698         999963
   gender         Female         Female         Female
   birthyear      1984           1983           1982
   maritalstatus  Divorced       Never-married  Married
   weeksworked16  0              0              53
   wageincome     NaN            NaN            50,000
   highestdegree  2.High School  2.High School  4\. Bachelors
   

通常，我们需要根据列值或几个列的值来选择行。在 pandas 中，我们可以通过布尔索引来完成此操作。这里，我们向loc访问器或方括号运算符传递布尔值向量（可以是序列）。布尔向量需要与 DataFrame 的索引相同。

1. 让我们尝试使用 NLS DataFrame 上的nightlyhrssleep列。我们想要一个布尔序列，对于每晚睡眠 6 小时或更少（第 33 个百分位数）的人为True，如果nightlyhrssleep大于 6 或缺失则为False。sleepcheckbool = nls97.nightlyhrssleep<=lowsleepthreshold创建布尔序列。如果我们显示sleepcheckbool的前几个值，我们将看到我们得到了预期的值。我们还可以确认sleepcheckbool的索引等于nls97的索引：

   nls97.nightlyhrssleep.head()
   personid
   100061     6
   100139     8
   100284     7
   100292     nan
   100583     6
   Name: nightlyhrssleep, dtype: float64
   lowsleepthreshold = nls97.nightlyhrssleep.quantile(0.33)
   lowsleepthreshold
   6.0
   sleepcheckbool = nls97.nightlyhrssleep<=lowsleepthreshold
   sleepcheckbool.head()
   personid
   100061    True
   100139    False
   100284    False
   100292    False
   100583    True
   Name: nightlyhrssleep, dtype: bool
   sleepcheckbool.index.equals(nls97.index)
   True
   

由于sleepcheckbool序列与nls97的索引相同，我们可以直接将其传递给loc访问器以创建一个包含每晚睡眠 6 小时或更少的人的 DataFrame。这里有一些 pandas 的魔法。它为我们处理索引对齐：

lowsleep = nls97.loc[sleepcheckbool]
lowsleep.shape
(3067, 88)

1. 我们本可以在一步中创建数据的lowsleep子集，这是我们通常会做的，除非我们需要布尔序列用于其他目的：

   lowsleep = nls97.loc[nls97.nightlyhrssleep<=lowsleepthreshold]
   lowsleep.shape
   (3067, 88)
   

2. 我们可以向loc访问器传递更复杂的条件并评估多个列的值。例如，我们可以选择nightlyhrssleep小于或等于阈值且childathome（在家居住的儿童数量）大于或等于3的行：

   lowsleep3pluschildren = \
     nls97.loc[(nls97.nightlyhrssleep<=lowsleepthreshold)
       & (nls97.childathome>=3)]
   lowsleep3pluschildren.shape
   (623, 88)
   

nls97.loc[(nls97.nightlyhrssleep<=lowsleepthreshold) & (nls97.childathome>3)]中的每个条件都放在括号内。如果省略括号，将生成错误。&运算符相当于标准 Python 中的and，意味着必须两个条件都为True，行才能被选中。如果我们想选择如果任一条件为True的行，我们可以使用|表示or。

1. 最后，我们可以同时选择行和列。逗号左边的表达式选择行，而逗号右边的列表选择列：

   lowsleep3pluschildren = \
     nls97.loc[(nls97.nightlyhrssleep<=lowsleepthreshold)
       & (nls97.childathome>=3),
       ['nightlyhrssleep','childathome']]
   lowsleep3pluschildren.shape
   (623, 2)
   

在上一两节中，我们使用了三种不同的工具从 pandas DataFrame 中选择列和行：[]括号操作符和两个 pandas 特定访问器loc和iloc。如果你是 pandas 的新手，这可能会有些令人困惑，但经过几个月后，你会清楚地知道在哪种情况下使用哪种工具。如果你带着相当多的 Python 和 NumPy 经验来到 pandas，你可能会发现[]操作符最为熟悉。然而，pandas 文档建议不要在生产代码中使用[]操作符。loc访问器用于通过布尔索引或索引标签选择行，而iloc访问器用于通过行号选择行。

这一部分是关于使用 pandas 选择列和行的简要入门。尽管我们没有对此进行过多详细说明，但涵盖了你需要了解的大部分内容，包括了解本书其余部分中 pandas 特定材料所需的一切。我们将在下一两节中开始应用这些知识，通过为我们的特征创建频率和汇总统计。

生成类别特征的频率

分类别特征可以是名义的或有序的。名义特征，例如性别、物种名称或国家，具有有限的可能值，可以是字符串或数值，但没有内在的数值意义。例如，如果国家用 1 代表阿富汗，2 代表阿尔巴尼亚，以此类推，数据是数值的，但对这些值进行算术运算是没有意义的。

有序特征也有有限的可能值，但与名义特征不同，值的顺序很重要。李克特量表评分（从 1 表示非常不可能到 5 表示非常可能）是一个有序特征的例子。尽管如此，通常不会进行算术运算，因为值之间没有统一和有意义的距离。

在我们开始建模之前，我们希望对可能使用的类别特征的所有可能值进行计数。这通常被称为单向频率分布。幸运的是，pandas 使这变得非常容易。我们可以快速从 pandas DataFrame 中选择列，并使用value_counts方法为每个类别值生成计数：

1. 让我们加载 NLS 数据，创建一个只包含数据前 20 列的 DataFrame，并查看数据类型：

   nls97 = pd.read_csv("data/nls97.csv")
   nls97.set_index("personid", inplace=True)
   nls97abb = nls97.iloc[:,:20]
   nls97abb.dtypes
   loc and iloc accessors. The colon to the left of the comma indicates that we want all the rows, while :20 to the right of the comma gets us the first 20 columns.
   

2. 上一段代码中的所有对象类型列都是类别特征。我们可以使用value_counts来查看maritalstatus每个值的计数。我们还可以使用dropna=False来让value_counts显示缺失值（NaN）：

   nls97abb.maritalstatus.value_counts(dropna=False)
   Married          3066
   Never-married    2766
   NaN              2312
   Divorced         663
   Separated        154
   Widowed          23
   Name: maritalstatus, dtype: int64
   

3. 如果我们只想得到缺失值的数量，我们可以链式调用isnull和sum方法。isnull在maritalstatus缺失时返回包含True值的布尔序列，否则返回False。然后sum计算True值的数量，因为它将True值解释为 1，将False值解释为 0：

   nls97abb.maritalstatus.isnull().sum()
   2312
   

4. 你可能已经注意到，maritalstatus的值默认按频率排序。你可以通过排序索引按值进行字母排序。我们可以利用value_counts返回一个以值为索引的序列这一事实来做到这一点：

   marstatcnt = nls97abb.maritalstatus.value_counts(dropna=False)
   type(marstatcnt)
   <class 'pandas.core.series.Series'>
   marstatcnt.index
   Index(['Married', 'Never-married', nan, 'Divorced', 'Separated', 'Widowed'], dtype='object')
   

5. 要对索引进行排序，我们只需调用sort_index：

   marstatcnt.sort_index()
   Divorced         663
   Married          3066
   Never-married    2766
   Separated        154
   Widowed          23
   NaN              2312
   Name: maritalstatus, dtype: int64
   

6. 当然，我们也可以通过nls97.maritalstatus.value_counts(dropna=False).sort_index()一步得到相同的结果。我们还可以通过将normalize设置为True来显示比率而不是计数。在下面的代码中，我们可以看到 34%的回应是Married（注意我们没有将dropna设置为True，所以缺失值已被排除）：

   nls97.maritalstatus.\
     value_counts(normalize=True, dropna=False).\
        sort_index()
   
   Divorced             0.07
   Married              0.34
   Never-married        0.31
   Separated            0.02
   Widowed              0.00
   NaN                  0.26
   Name: maritalstatus, dtype: float64
   

7. 当一列有有限数量的值时，pandas 的类别数据类型可以比对象数据类型更有效地存储数据。由于我们已经知道所有对象列都包含类别数据，我们应该将这些列转换为类别数据类型。在下面的代码中，我们创建了一个包含对象列名称的列表，catcols。然后，我们遍历这些列，并使用astype将数据类型更改为category：

   catcols = nls97abb.select_dtypes(include=["object"]).columns
   for col in nls97abb[catcols].columns:
   ...      nls97abb[col] = nls97abb[col].astype('category')
   ... 
   nls97abb[catcols].dtypes
   gender                   category
   maritalstatus            category
   weeklyhrscomputer        category
   weeklyhrstv              category
   highestdegree            category
   govprovidejobs           category
   govpricecontrols         category
   dtype: object
   

8. 让我们检查我们的类别特征中的缺失值。gender没有缺失值，highestdegree的缺失值非常少。但govprovidejobs（政府应该提供工作）和govpricecontrols（政府应该控制价格）的绝大多数值都是缺失的。这意味着这些特征可能对大多数建模没有用：

   nls97abb[catcols].isnull().sum()
   gender               0
   maritalstatus        2312
   weeklyhrscomputer    2274
   weeklyhrstv          2273
   highestdegree        31
   govprovidejobs       7151
   govpricecontrols     7125
   dtype: int64
   

9. 我们可以通过将value_counts调用传递给apply来一次生成多个特征的频率。我们可以使用filter来选择我们想要的列——在这种情况下，所有名称中包含*gov*的列。注意，由于我们没有将dropna设置为False，因此已经省略了每个特征的缺失值：

    nls97abb.filter(like="gov").apply(pd.value_counts, normalize=True)
                    govprovidejobs    govpricecontrols
   1\. Definitely              0.25                0.54
   2\. Probably                0.34                0.33
   3\. Probably not            0.25                0.09
   4\. Definitely not          0.16                0.04
   

10. 我们可以在数据的一个子集上使用相同的频率。例如，如果我们只想查看已婚人士对政府角色问题的回答，我们可以通过在filter之前放置nls97abb[nls97abb.maritalstatus=="Married"]来执行这个子集操作：

     nls97abb.loc[nls97abb.maritalstatus=="Married"].\
     filter(like="gov").\
       apply(pd.value_counts, normalize=True)
                     govprovidejobs    govpricecontrols
    1\. Definitely              0.17                0.46
    2\. Probably                0.33                0.38
    3\. Probably not            0.31                0.11
    4\. Definitely not          0.18                0.05
    

11. 在这种情况下，由于只有两个*gov*列，可能更容易执行以下操作：

     nls97abb.loc[nls97abb.maritalstatus=="Married",
       ['govprovidejobs','govpricecontrols']].\
       apply(pd.value_counts, normalize=True)
                      govprovidejobs     govpricecontrols
    1\. Definitely               0.17                 0.46
    2\. Probably                 0.33                 0.38
    3\. Probably not             0.31                 0.11
    4\. Definitely not           0.18                 0.05
    

尽管如此，使用filter通常会更简单，因为在具有相似名称的特征组上执行相同的清理或探索任务并不罕见。

有时候，我们可能希望将连续或离散特征建模为分类特征。NLS DataFrame 包含 highestgradecompleted。从 5 年级到 6 年级的增长可能不如从 11 年级到 12 年级的增长对目标的影响重要。让我们创建一个二进制特征——即当一个人完成了 12 年级或以上时为 1，如果他们完成的少于这个数则为 0，如果 highestgradecompleted 缺失则为缺失。

1. 不过，我们首先需要做一些清理工作。highestgradecompleted 有两个逻辑缺失值——一个 pandas 识别为缺失的实际 NaN 值和一个调查设计者意图让我们在大多数情况下也视为缺失的 95 值。让我们在继续之前使用 replace 来修复它：

   nls97abb.highestgradecompleted.\
     replace(95, np.nan, inplace=True)
   

2. 我们可以使用 NumPy 的 where 函数根据 highestgradecompleted 的值分配 highschoolgrad 的值。如果 highestgradecompleted 为空（NaN），我们将 NaN 分配给我们的新列 highschoolgrad。如果 highestgradecompleted 的值不为空，下一个子句检查值是否小于 12，如果是，则将 highschoolgrad 设置为 0，否则设置为 1。我们可以通过使用 groupby 来获取 highschoolgrad 每个级别的 highestgradecompleted 的最小值和最大值来确认新列 highschoolgrad 包含我们想要的值：

   nls97abb['highschoolgrad'] = \
     np.where(nls97abb.highestgradecompleted.isnull(),np.nan, \
     np.where(nls97abb.highestgradecompleted<12,0,1))
   
   nls97abb.groupby(['highschoolgrad'], dropna=False) \
     ['highestgradecompleted'].agg(['min','max','size'])
                     min       max       size
   highschoolgrad                
   0                   5        11       1231
   1                  12        20       5421
   nan               nan       nan       2332
    nls97abb['highschoolgrad'] = \
   ...  nls97abb['highschoolgrad'].astype('category')
   

虽然 12 作为将我们的新特征 highschoolgrad 分类为类的阈值在概念上是合理的，但如果我们打算将 highschoolgrad 作为目标使用，这可能会带来一些建模挑战。存在相当大的类别不平衡，highschoolgrad 等于 1 的类别是 0 类的四倍多。我们应该探索使用更多组来表示 highestgradecompleted。

1. 使用 pandas 实现这一点的其中一种方法是使用 qcut 函数。我们可以将 qcut 的 q 参数设置为 6 以创建尽可能均匀分布的六个组。现在这些组更接近平衡：

   nls97abb['highgradegroup'] = \
     pd.qcut(nls97abb['highestgradecompleted'], 
      q=6, labels=[1,2,3,4,5,6])
   
   nls97abb.groupby(['highgradegroup'])['highestgradecompleted'].\
       agg(['min','max','size'])
                     min         max      size
   highgradegroup                
   1                   5          11       1231
   2                  12          12       1389
   3                  13          14       1288
   4                  15          16       1413
   5                  17          17        388
   6                  18          20        943
   nls97abb['highgradegroup'] = \
       nls97abb['highgradegroup'].astype('category')
   

2. 最后，我发现通常生成所有分类特征的频率并将其保存下来很有帮助，这样我就可以稍后参考。每当我对数据进行一些可能改变这些频率的更改时，我都会重新运行那段代码。以下代码遍历所有数据类型为分类数据的列，并运行 value_counts：

    freqout = open('views/frequencies.txt', 'w') 
    for col in nls97abb.select_dtypes(include=["category"]):
         print(col, "----------------------",
           "frequencies",
         nls97abb[col].value_counts(dropna=False).sort_index(),
           "percentages",
         nls97abb[col].value_counts(normalize=True).\
           sort_index(),
         sep="\n\n", end="\n\n\n", file=freqout)
   
    freqout.close()
   

这些是在您的数据中生成分类特征的单一频率的关键技术。真正的明星是 value_counts 方法。我们可以一次创建一个 Series 的频率使用 value_counts，也可以使用 apply 对多个列进行操作，或者遍历多个列并在每次迭代中调用 value_counts。我们已经在本节中查看了一些示例。接下来，让我们探讨一些检查连续特征分布的技术。

生成连续和离散特征的摘要统计量

获取连续或离散特征的分布感比分类特征要复杂一些。连续特征可以取无限多个值。一个连续特征的例子是体重，因为某人可能重 70 公斤，或者 70.1，或者 70.01。离散特征具有有限个值，例如看到的鸟的数量，或者购买苹果的数量。思考它们之间差异的一种方式是，离散特征通常是已经被计数的东西，而连续特征通常是通过测量、称重或计时来捕捉的。

连续特征通常会被存储为浮点数，除非它们被限制为整数。在这种情况下，它们可能以整数的形式存储。例如，个人的年龄是连续的，但通常会被截断为整数。

对于大多数建模目的，连续特征和离散特征被同等对待。我们不会将年龄建模为分类特征。我们假设年龄在 25 岁到 26 岁之间的间隔与 35 岁到 36 岁之间的间隔具有大致相同的意义，尽管在极端情况下这种假设会失效。人类 1 岁到 2 岁之间的间隔与 71 岁到 72 岁之间的间隔完全不同。数据分析师和科学家通常对连续特征和目标之间的假设线性关系持怀疑态度，尽管当这种关系成立时建模会更容易。

要了解连续特征（或离散特征）的分布，我们必须检查其中心趋势、形状和分布范围。关键摘要统计量包括均值和中位数用于中心趋势，偏度和峰度用于形状，以及范围、四分位数范围、方差和标准差用于分布范围。在本节中，我们将学习如何使用 pandas，辅以SciPy库，来获取这些统计量。我们还将讨论建模的重要影响。

在本节中，我们将使用 COVID-19 数据。数据集包含每个国家的总病例和死亡数，以及截至 2021 年 6 月的人口统计数据。

注意

我们的世界数据在ourworldindata.org/coronavirus-source-data提供 COVID-19 公共使用数据。本节中使用的数据是在 2021 年 7 月 9 日下载的。数据中的列比我包含的要多。我根据国家创建了地区列。

按照以下步骤生成摘要统计量：

1. 让我们将 COVID .csv文件加载到 pandas 中，设置索引，并查看数据。有 221 行和 16 列。我们设置的索引iso_code为每一行包含一个唯一值。我们使用sample来随机查看两个国家，而不是前两个（我们为random_state设置了一个值，以便每次运行代码时都能得到相同的结果）：

   import pandas as pd
   import numpy as np
   import scipy.stats as scistat
   covidtotals = pd.read_csv("data/covidtotals.csv",
       parse_dates=['lastdate'])
   covidtotals.set_index("iso_code", inplace=True)
   covidtotals.shape
   (221, 16)
   covidtotals.index.nunique()
   221
   covidtotals.sample(2, random_state=6).T
   iso_code                         ISL               CZE
   lastdate                  2021-07-07        2021-07-07
   location                     Iceland           Czechia
   total_cases                    6,555         1,668,277
   total_deaths                      29            30,311
   total_cases_mill              19,209           155,783
   total_deaths_mill                 85             2,830
   population                   341,250        10,708,982
   population_density                 3               137
   median_age                        37                43
   gdp_per_capita                46,483            32,606
   aged_65_older                     14                19
   total_tests_thous                NaN               NaN
   life_expectancy                   83                79
   hospital_beds_thous                3                 7
   diabetes_prevalence                5                 7
   region                Western Europe    Western Europe
   

只需看这两行，我们就可以看到冰岛和捷克在病例和死亡人数方面的显著差异，即使在人口规模方面也是如此。（total_cases_mill 和 total_deaths_mill 分别表示每百万人口中的病例和死亡人数。）数据分析师非常习惯于思考数据中是否还有其他因素可以解释捷克比冰岛病例和死亡人数显著更高的原因。从某种意义上说，我们总是在进行特征选择。

1. 让我们来看看每列的数据类型和空值数量。几乎所有列都是连续的或离散的。我们有关于病例和死亡的数据，分别有 192 行和 185 行。我们必须要做的一个重要数据清洗任务是确定我们能对那些对于我们的目标有缺失值的国家的数据做些什么。我们将在稍后讨论如何处理缺失值：

   covidtotals.info()
   <class 'pandas.core.frame.DataFrame'>
   Index: 221 entries, AFG to ZWE
   Data columns (total 16 columns):
    #   Column             Non-Null Count         Dtype 
   ---  -------            --------------  --------------
    0   lastdate             221 non-null  datetime64[ns]
    1   location             221 non-null          object
    2   total_cases          192 non-null         float64
    3   total_deaths         185 non-null         float64
    4   total_cases_mill     192 non-null         float64
    5   total_deaths_mill    185 non-null         float64
    6   population           221 non-null         float64
    7   population_density   206 non-null         float64
    8   median_age           190 non-null         float64
    9   gdp_per_capita       193 non-null         float64
    10  aged_65_older        188 non-null         float64
    11  total_tests_thous     13 non-null         float64
    12  life_expectancy      217 non-null         float64
    13  hospital_beds_thous  170 non-null         float64
    14  diabetes_prevalence  200 non-null         float64
    15  region               221 non-null          object
   dtypes: datetime64ns, float64(13), object(2)
   memory usage: 29.4+ KB
   

2. 现在，我们已经准备好检查一些特征的分部情况。我们可以通过使用 describe 方法来获取我们想要的绝大部分的摘要统计信息。平均值和中位数（50%）是分布中心的良好指标，各有其优势。注意到平均值和中位数之间的显著差异，作为偏斜的指示也是很好的。例如，我们可以看到每百万病例的平均值几乎是中位数的两倍，分别是 36.7 千和 19.5 千。这是一个明显的正偏斜指标。死亡人数每百万的情况也是如此。

病例和死亡人数的四分位距也相当大，两种情况下的第 75 百分位数值大约是第 25 百分位数值的 25 倍。我们可以将这一点与 65 岁及以上人口比例和糖尿病患病率进行比较，其中第 75 百分位数分别是第 25 百分位数的四倍或两倍。我们可以立即得出结论，这两个可能的特征（aged_65_older 和 diabetes_prevalence）必须做大量工作来解释我们目标中的巨大方差：

 keyvars = ['location','total_cases_mill','total_deaths_mill',
...  'aged_65_older','diabetes_prevalence']
 covidkeys = covidtotals[keyvars]
 covidkeys.describe()
total_cases_mill total_deaths_mill aged_65_older diabetes_prevalence
count        192.00       185.00    188.00     200.00
mean      36,649.37       683.14      8.61       8.44
std       41,403.98       861.73      6.12       4.89
min            8.52         0.35      1.14       0.99
25%        2,499.75        43.99      3.50       5.34
50%       19,525.73       293.50      6.22       7.20
75%       64,834.62     1,087.89     13.92      10.61
max      181,466.38     5,876.01     27.05      30.53

1. 我有时发现查看十分位数值有助于更好地了解分布。quantile 方法可以接受一个百分位数值，例如 quantile(0.25) 表示第 25 百分位数，或者一个列表或元组，例如 quantile((0.25,0.5)) 表示第 25 和第 50 百分位数。在下面的代码中，我们使用 NumPy 的 arange 函数（np.arange(0.0, 1.1, 0.1)) 生成一个从 0.0 到 1.0，增量是 0.1 的数组。如果我们使用 covidkeys.quantile([0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0])，我们会得到相同的结果：

    covidkeys.quantile(np.arange(0.0, 1.1, 0.1))
         total_cases_mill  total_deaths_mill  aged_65_older  diabetes_prevalence
   0.00         8.52       0.35       1.14     0.99
   0.10       682.13      10.68       2.80     3.30
   0.20     1,717.39      30.22       3.16     4.79
   0.30     3,241.84      66.27       3.86     5.74
   0.40     9,403.58      145.06      4.69     6.70
   0.50     19,525.73     293.50      6.22     7.20
   0.60     33,636.47     556.43      7.93     8.32
   0.70     55,801.33     949.71     11.19    10.08
   0.80     74,017.81    1,333.79    14.92    11.62
   0.90     94,072.18    1,868.89    18.85    13.75
   1.00    181,466.38    5,876.01    27.05    30.53
   

对于病例、死亡和糖尿病患病率，大部分的范围（最小值和最大值之间的距离）位于分布的最后 10%。这一点对于死亡尤其如此。这暗示了可能存在建模问题，并邀请我们仔细查看异常值，我们将在下一节中这样做。

1. 一些机器学习算法假设我们的特征具有正态分布（也称为高斯分布），它们分布对称（具有低偏度），并且它们的尾部相对正常（既不过度高也不过度低的高斯峰）。我们迄今为止看到的统计数据已经表明，我们的两个可能的目标——即总人口中的总病例数和死亡人数——具有很高的正偏度。让我们通过计算一些特征的偏度和峰度来对此进行更细致的分析。对于高斯分布，我们期望偏度为 0 附近，峰度为 3。total_deaths_mill的偏度和峰度值值得关注，而total_cases_mill和aged_65_older特征具有过低的峰度（瘦尾巴）：

   covidkeys.skew()
   total_cases_mill        1.21
   total_deaths_mill       2.00
   aged_65_older           0.84
   diabetes_prevalence     1.52
   dtype: float64
    covidkeys.kurtosis()
   total_cases_mill        0.91
   total_deaths_mill       6.58
   aged_65_older          -0.56
   diabetes_prevalence     3.31
   dtype: float64
   

2. 我们还可以通过遍历keyvars列表中的特征并运行scistat.shapiro(covidkeys[var].dropna())来显式测试每个分布的正态性。请注意，我们需要使用dropna删除缺失值以运行测试。p 值小于 0.05 表明我们可以拒绝正态性的零假设，这对于四个特征中的每一个都是如此：

   for var in keyvars[1:]:
         stat, p = scistat.shapiro(covidkeys[var].dropna())
         print("feature=", var, "     p-value=", '{:.6f}'.format(p))
   
   feature= total_cases_mill       p-value= 0.000000
   feature= total_deaths_mill      p-value= 0.000000
   feature= aged_65_older          p-value= 0.000000
   feature= diabetes_prevalence    p-value= 0.000000
   

这些结果应该让我们在考虑参数模型，如线性回归时停下来思考。没有任何分布近似正态分布。然而，这并不是决定性的。这并不像在具有正态分布特征时使用某些模型（比如 k-最近邻）而在不具有时使用非参数模型那么简单。

在做出任何建模决策之前，我们希望进行额外的数据清洗。例如，我们可能决定删除异常值或确定数据转换是合适的。我们将在本书的几个章节中探索转换，例如对数和多项式转换。

本节向您展示了如何使用 pandas 和 SciPy 来了解连续和离散特征的分布情况，包括它们的中心趋势、形状和分布范围。为任何可能包含在我们建模中的特征或目标生成这些统计数据是有意义的。这也指引我们进一步的工作方向，以准备我们的数据进行分析。我们需要识别缺失值和异常值，并找出我们将如何处理它们。我们还应该可视化连续特征的分布。这很少不会带来额外的见解。我们将在下一节学习如何识别异常值，并在下一节创建可视化。

在单变量分析中识别极端值和异常值

异常值可以被视为具有特征值或特征值之间关系非常不寻常的观测值，以至于它们无法解释数据中其他关系。这对于建模很重要，因为我们不能假设异常值将对我们的参数估计产生中性影响。有时，我们的模型会努力构建参数估计，以解释异常观测值中的模式，这可能会损害模型对所有其他观测值的解释或预测能力。如果你曾经花了好几天时间尝试解释一个模型，最终发现一旦移除几个异常值，你的系数和预测就完全改变了，请举手。

我应该快速补充一下，目前还没有关于异常值的统一定义。我提供前面的定义是为了在本书中使用，因为它有助于我们区分异常值，正如我描述的那样，以及极端值。这两个概念之间有相当大的重叠，但许多极端值并不是异常值。这是因为这样的值反映了特征中的自然和可解释的趋势，或者因为它们反映了与数据中观察到的特征之间的相同关系。反之亦然。有些异常值不是极端值。例如，一个目标值可能正好位于分布的中间，但具有相当意外的预测值。

因此，对于我们的建模来说，很难在没有参考多元关系的情况下说一个特定的特征或目标值是异常值。但是，当我们在单变量分析中看到位于中心左侧或右侧的值时，这至少应该引起我们的注意。这应该促使我们进一步调查该值，包括检查其他特征值。我们将在下一章中更详细地探讨多元关系。在这里，以及在下节关于可视化的内容中，我们将专注于在查看单个变量时识别极端值和异常值。

识别极端值的一个良好起点是查看其与分布中间的距离。实现这一目标的一种常见方法是将每个值与四分位数范围（IQR）的距离进行计算，即第一个四分位数值与第三个四分位数值之间的距离。我们通常将任何高于第三四分位数 1.5 倍或低于第一四分位数的值标记出来。我们可以使用这种方法来识别 COVID-19 数据中的异常值。

让我们开始吧：

1. 首先，让我们导入我们将需要的库。除了 pandas 和 NumPy，我们还将使用 Matplotlib 和 statsmodels 来创建图表。我们还将加载 COVID 数据并选择所需的变量：

   import pandas as pd
   import numpy as np
   import matplotlib.pyplot as plt
   import statsmodels.api as sm
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   keyvars = ['location','total_cases_mill','total_deaths_mill',
     'aged_65_older','diabetes_prevalence','gdp_per_capita']
   covidkeys = covidtotals[keyvars]
   

2. 让我们来看看total_cases_mill。我们得到第一和第三四分位数，并计算四分位数范围，1.5*(thirdq-firstq)。然后，我们计算它们的阈值以确定高和低极端值，分别是interquartilerange+thirdq和firstq-interquartilerange（如果你熟悉箱线图，你会注意到这是用于箱线图须的相同计算；我们将在下一节中介绍箱线图）：

   thirdq, firstq = covidkeys.total_cases_mill.quantile(0.75), covidkeys.total_cases_mill.quantile(0.25)
   interquartilerange = 1.5*(thirdq-firstq)
   extvalhigh, extvallow = interquartilerange+thirdq, firstq-interquartilerange
   print(extvallow, extvalhigh, sep=" <--> ")
   -91002.564625 <--> 158336.930375
   

3. 这个计算表明，任何total_cases_mill的值，如果高于 158,337，都可以被认为是极端的。我们可以忽略低端的极端值，因为它们将是负数：

   covidtotals.loc[covidtotals.total_cases_mill>extvalhigh].T
   iso_code                   AND         MNE         SYC
   lastdate            2021-07-07  2021-07-07  2021-07-07
   location            Andorra     Montenegro  Seychelles
   total_cases          14,021        100,392      16,304
   total_deaths            127          1,619          71
   total_cases_mill    181,466        159,844     165,792
   total_deaths_mill     1,644          2,578         722
   population           77,265        628,062      98,340
   population_density      164             46         208
   median_age              NaN             39          36
   gdp_per_capita          NaN         16,409      26,382
   aged_65_older           NaN             15           9
   total_tests_thous       NaN            NaN         NaN
   life_expectancy          84             77          73
   hospital_beds_thous     NaN              4           4
   diabetes_prevalence       8             10          11
   region              Western        Eastern        East 
                        Europe         Europe      Africa
   

4. 安道尔、黑山和塞舌尔的所有total_cases_mill都超过了阈值。这促使我们探索这些国家可能异常的其他方式，以及我们的特征是否能够捕捉到这一点。我们不会在这里深入多元分析，因为我们将那部分内容放在下一章中，但开始思考为什么这些极端值可能或可能没有意义是个好主意。在整个数据集上有一个平均值可能有助于我们：

   covidtotals.mean()
   total_cases               963,933
   total_deaths              21,631
   total_cases_mill          36,649
   total_deaths_mill         683
   population                35,134,957
   population_density        453
   median_age                30
   gdp_per_capita            19,141
   aged_65_older             9
   total_tests_thous         535
   life_expectancy           73
   hospital_beds_thous       3
   diabetes_prevalence       8 
   

这三个国家与其他国家的主要区别在于它们的人口非常少。令人惊讶的是，每个国家的人口密度都比平均水平低得多。这与你预期的相反，值得我们在这本书的整个分析中进一步考虑。

IQR 计算的替代方法

使用四分位数范围来识别极端值的另一种方法是使用与平均值几个标准差的位置，比如说 3 个。这种方法的一个缺点是，它比使用四分位数范围更容易受到极端值的影响。

我认为为我的数据中的所有关键目标和特征生成这种分析是有帮助的，因此让我们自动化识别极端值的方法。我们还应该将结果输出到文件中，这样我们就可以在我们需要时使用它们：

1. 让我们定义一个函数，getextremevalues，它遍历我们的 DataFrame（除了包含位置列的第一个列之外的所有列），计算该列的四分位数范围，选择所有高于该列高阈值或低于低阈值的观测值，然后将结果追加到一个新的 DataFrame（dfout）中：

   def getextremevalues(dfin):
         dfout = pd.DataFrame(columns=dfin.columns, 
                              data=None)
         for col in dfin.columns[1:]:
           thirdq, firstq = dfin[col].quantile(0.75), \
             dfin[col].quantile(0.25)
           interquartilerange = 1.5*(thirdq-firstq)
           extvalhigh, extvallow = \
             interquartilerange+thirdq, \
             firstq-interquartilerange
           df = dfin.loc[(dfin[col]>extvalhigh) | (dfin[col]<extvallow)]
           df = df.assign(varname = col,
             threshlow = extvallow,
             threshhigh = extvalhigh)
           dfout = pd.concat([dfout, df])
       return dfout
   

2. 现在，我们可以将我们的covidkeys DataFrame 传递给getextremevalues函数，以获取包含每个列的极端值的 DataFrame。然后，我们可以显示每个列的极端值数量，这告诉我们人口中每百万人的总死亡数（total_deaths_mill）有四个极端值。现在，我们可以将数据输出到 Excel 文件中：

   extremevalues = getextremevalues(covidkeys)
   extremevalues.varname.value_counts()
   gdp_per_capita         9
   diabetes_prevalence    8
   total_deaths_mill      4
   total_cases_mill       3
   Name: varname, dtype: int64
   extremevalues.to_excel("views/extremevaluescases.xlsx")
   

3. 让我们更仔细地看看每百万极端死亡值。我们可以查询我们刚刚创建的 DataFrame，以获取total_deaths_mill的threshhigh值，该值为2654。我们还可以获取具有极端值的那些国家的其他关键特征，因为我们已经将那些数据包含在新的 DataFrame 中：

   extremevalues.loc[extremevalues.varname=="total_deaths_mill",
       'threshhigh'][0]
   2653.752
   extremevalues.loc[extremevalues.varname=="total_deaths_mill",
         keyvars].sort_values(['total_deaths_mill'], ascending=False)
         location              total_cases_mill  total_deaths_mill
   PER   Peru                    62,830            5,876
   HUN   Hungary                 83,676            3,105
   BIH   Bosnia and Herzegovina  62,499            2,947
   CZE   Czechia                 155,783           2,830
        _65_older  diabetes_prevalence  gdp_per_capita  
   PER      7              6                12,237
   HUN     19              8                26,778
   BIH     17             10                11,714
   CZE     19              7                32,606
   

秘鲁、匈牙利、波斯尼亚和黑塞哥维那以及捷克共和国的total_deaths_mill值超过了极端值阈值。这三个国家中的一个突出特点是，它们 65 岁或以上人口百分比远高于平均水平（该特征的均值为 9，正如我们在先前的表中显示的那样）。尽管这些是死亡极端值，但老年人口百分比与 COVID 死亡之间的关系可能解释了其中很大一部分，并且可以在不过度拟合这些极端案例的情况下做到这一点。我们将在下一章中介绍一些提取这些策略的方法。

到目前为止，我们讨论了异常值和极端值，但没有提及分布形状。我们迄今为止所暗示的是，极端值是一个罕见值——比分布中心附近的值要罕见得多。但是，当特征分布接近正态分布时，这最有意义。另一方面，如果一个特征具有均匀分布，那么一个非常高的值并不比任何其他值罕见。

在实践中，我们会考虑相对于特征分布的极端值或异常值。分位数-分位数（Q-Q）图可以通过允许我们相对于理论分布（正态分布、均匀分布、对数分布或其他）图形化地查看它来提高我们对该分布的感觉：让我们看一下：

1. 让我们创建一个与正态分布相关的每百万总病例的 Q-Q 图。我们可以使用statsmodels库来完成这个任务：

   sm.qqplot(covidtotals[['total_cases_mill']]. \
     sort_values(['total_cases_mill']).dropna(),line='s')
   )
   plt.show()
   

这会产生以下图表：

图 1.1 – 每百万总病例的 Q-Q 图

这个 Q-Q 图清楚地表明，各国总病例的分布不是正态分布。我们可以通过数据点与红线偏离的程度来看到这一点。这是一个我们预期来自具有某些正偏斜的分布的 Q-Q 图。这与我们已为总病例特征计算出的摘要统计一致。它进一步强化了我们正在形成的认识，即我们需要对参数模型保持谨慎，我们可能需要考虑的不仅仅是单个或两个异常观测值。

让我们看看一个分布略接近正态分布的特征的 Q-Q 图。在 COVID 数据中没有很好的候选者，所以我们将使用美国国家海洋和大气管理局 2019 年陆地温度的数据。

数据备注

陆地温度数据集包含了来自全球超过 12,000 个站点的 2019 年平均温度读数（摄氏度），尽管大多数站点位于美国。原始数据是从全球历史气候学网络综合数据库中检索的。美国国家海洋和大气管理局已将其公开供公众使用，网址为 www.ncdc.noaa.gov/data-access/land-based-station-data/land-based-datasets/global-historical-climatology-network-monthly-version-2。

1. 首先，让我们将数据加载到 pandas DataFrame 中，并对温度特征 avgtemp 运行一些描述性统计。我们必须向正常的 describe 输出中添加一些百分位数统计，以更好地了解值的范围。avgtemp 是每个 12,095 个气象站一年的平均温度。所有站点的平均温度为 11.2 摄氏度。中位数为 10.4。然而，有一些非常低的值，包括 14 个平均温度低于 -25 的气象站。这导致了一中度负偏斜，尽管偏斜和峰度都更接近正态分布的预期：

   landtemps = pd.read_csv("data/landtemps2019avgs.csv")
   landtemps.avgtemp.describe(percentiles=[0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95])
   count        12,095.0
   mean             11.2
   std               8.6
   min             -60.8
   5%               -0.7
   10%               1.7
   25%               5.4
   50%              10.4
   75%              16.9
   90%              23.1
   95%              27.0
   max              33.9
   Name: avgtemp, dtype: float64
    landtemps.loc[landtemps.avgtemp<-25,'avgtemp'].count()
   14
   landtemps.avgtemp.skew()
   -0.2678382583481768
   landtemps.avgtemp.kurtosis()
   2.1698313707061073
   

2. 现在，让我们看一下平均温度的 Q-Q 图：

   sm.qqplot(landtemps.avgtemp.sort_values().dropna(), line='s')
   plt.tight_layout()
   plt.show()
   

这产生了以下图表：

图 1.2 – 平均温度的 Q-Q 图

在大多数范围内，平均温度的分布看起来非常接近正态分布。例外的是极低温度，这导致了一小部分负偏斜。在高端也有一些偏离正态分布的情况，尽管这并不是一个大问题（你可能已经注意到，具有负偏斜的特征的 Q-Q 图具有雨伞状形状，而具有正偏斜的特征，如总病例数，则具有更多碗状形状）。

我们在理解可能特征和目标分布以及识别极端值和异常值的相关努力中取得了良好的开端。当我们构建、改进和解释模型时，这些重要信息对我们来说至关重要。但我们还可以做更多的事情来提高我们对数据的直觉。一个好的下一步是构建关键特征的可视化。

使用直方图、箱线图和提琴图来检查特征分布

我们已经生成了构成直方图或箱线图数据点的许多数字。但当我们看到数据以图形方式表示时，我们往往能更好地理解数据。我们看到观测值围绕着平均值聚集，我们注意到尾部的尺寸，我们看到似乎极端的值。

使用直方图

按以下步骤创建直方图：

1. 在本节中，我们将同时处理 COVID 数据和温度数据。除了我们迄今为止使用的库之外，我们还必须导入 Seaborn，以便比在 Matplotlib 中更容易地创建一些图表：

   import pandas as pd
   import numpy as np
   import matplotlib.pyplot as plt
   import seaborn as sns
   landtemps = pd.read_csv("data/landtemps2019avgs.csv")
   covidtotals = pd.read_csv("data/covidtotals.csv", parse_dates=["lastdate"])
   covidtotals.set_index("iso_code", inplace=True)
   

2. 现在，让我们创建一个简单的直方图。我们可以使用 Matplotlib 的hist方法创建每百万人口中总病例的直方图。我们还将绘制均值和中位数的线条：

   plt.hist(covidtotals['total_cases_mill'], bins=7)
   plt.axvline(covidtotals.total_cases_mill.mean(), color='red',
      linestyle='dashed', linewidth=1, label='mean')
   plt.axvline(covidtotals.total_cases_mill.median(), color='black',
      linestyle='dashed', linewidth=1, label='median')
   plt.title("Total COVID Cases")
   plt.xlabel('Cases per Million')
   plt.ylabel("Number of Countries")
   plt.legend()
   plt.show()
   

这产生了以下图表：

图 1.3 – 总 COVID 病例

这个直方图突出显示的总体分布的一个方面是，大多数国家（192 个中的 100 多个）都在第一个区间内，即每百万人口 0 例到 25,000 例之间。在这里，我们可以看到正偏斜，由于极端高值，均值被拉向右边。这与我们在上一节中使用 Q-Q 图发现的情况一致。

1. 让我们创建一个来自陆地温度数据集的平均温度直方图：

   plt.hist(landtemps['avgtemp'])
   plt.axvline(landtemps.avgtemp.mean(), color='red', linestyle='dashed', linewidth=1, label='mean')
   plt.axvline(landtemps.avgtemp.median(), color='black', linestyle='dashed', linewidth=1, label='median')
   plt.title("Average Land Temperatures")
   plt.xlabel('Average Temperature')
   plt.ylabel("Number of Weather Stations")
   plt.legend()
   plt.show()
   

这产生了以下图表：

图 1.4 – 平均陆地温度

来自陆地温度数据集的平均陆地温度直方图看起来相当不同。除了少数几个高度负值外，这种分布看起来更接近正态分布。在这里，我们可以看到均值和众数非常接近，分布看起来相当对称。

1. 我们应该看一下分布最左端的观测值。它们都在南极洲或加拿大最北部。在这里，我们必须怀疑是否应该将具有如此极端值的观测值包含在我们构建的模型中。然而，仅基于这些结果就做出这样的决定还为时过早。我们将在下一章中回到这个问题，当我们检查用于识别异常值的多变量技术时：

    landtemps.loc[landtemps.avgtemp<-25,['station','country','avgtemp']].\
   ...  sort_values(['avgtemp'], ascending=True)
         station               country          avgtemp
   827   DOME_PLATEAU_DOME_A   Antarctica      -60.8
   830   VOSTOK                Antarctica      -54.5
   837   DOME_FUJI             Antarctica      -53.4
   844   DOME_C_II             Antarctica      -50.5
   853   AMUNDSEN_SCOTT        Antarctica      -48.4
   842   NICO                  Antarctica      -48.4
   804   HENRY                 Antarctica      -47.3
   838   RELAY_STAT            Antarctica      -46.1
   828   DOME_PLATEAU_EAGLE    Antarctica      -43.0
   819   KOHNENEP9             Antarctica      -42.4
   1299  FORT_ROSS             Canada          -30.3
   1300  GATESHEAD_ISLAND      Canada          -28.7
   811   BYRD_STATION          Antarctica      -25.8
   816   GILL                  Antarctica      -25.5
   

同时可视化集中趋势、离散度和异常值的一个极好方法是使用箱线图。

使用箱线图

箱线图显示了四分位距，其中触须代表四分位距的 1.5 倍，以及超出该范围的数据点，这些数据点可以被认为是极端值。如果这个计算看起来很熟悉，那是因为我们之前在本章中用来识别极端值的就是这个计算！让我们开始吧：

1. 我们可以使用 Matplotlib 的boxplot方法创建每百万人口中总病例的箱线图。我们可以画箭头来表示四分位距（第一四分位数、中位数和第三四分位数）以及极端值阈值。阈值上方的三个圆圈可以被认为是极端值。从四分位距到极端值阈值的线通常被称为触须。通常在四分位距上方和下方都有触须，但在这个情况下，低于第一四分位数阈值的值将是负数：

   plt.boxplot(covidtotals.total_cases_mill.dropna(), labels=['Total Cases per Million'])
   plt.annotate('Extreme Value Threshold', xy=(1.05,157000), xytext=(1.15,157000), size=7, arrowprops=dict(facecolor='black', headwidth=2, width=0.5, shrink=0.02))
   plt.annotate('3rd quartile', xy=(1.08,64800), xytext=(1.15,64800), size=7, arrowprops=dict(facecolor='black', headwidth=2, width=0.5, shrink=0.02))
   plt.annotate('Median', xy=(1.08,19500), xytext=(1.15,19500), size=7, arrowprops=dict(facecolor='black', headwidth=2, width=0.5, shrink=0.02))
   plt.annotate('1st quartile', xy=(1.08,2500), xytext=(1.15,2500), size=7, arrowprops=dict(facecolor='black', headwidth=2, width=0.5, shrink=0.02))
   plt.title("Boxplot of Total Cases")
   plt.show()
   

这产生了以下图表：

图 1.5 – 总病例数的箱线图

仔细观察四分位数范围是有帮助的，特别是中位数在范围内的位置。对于这个箱线图，中位数位于范围的低端。这就是我们在正偏分布中看到的情况。

1. 现在，让我们为平均温度创建一个箱线图。所有的极端值现在都位于分布的低端。不出所料，鉴于我们已经看到的平均温度特征，中位数线比我们之前的箱线图更接近四分位数范围的中心（这次我们不会注释图表——我们上次这样做是为了解释目的）：

    plt.boxplot(landtemps.avgtemp.dropna(), labels=['Boxplot of Average Temperature'])
    plt.title("Average Temperature")
    plt.show()
   

这产生了以下图表：

图 1.6 – 平均温度的箱线图

直方图帮助我们了解分布的分布情况，而箱线图则使识别异常值变得容易。我们可以通过小提琴图在一个图表中很好地了解分布的分布情况和异常值。

使用小提琴图

小提琴图结合了直方图和箱线图，在一个图表中显示。它们显示了四分位数范围、中位数和触须，以及所有值范围内的观测频率。

让我们开始吧：

1. 我们可以使用 Seaborn 为每百万个 COVID 病例和平均温度特征创建 violin 图。我在这里使用 Seaborn 而不是 Matplotlib，因为我更喜欢其默认的 violin 图选项：

   import seaborn as sns
   fig = plt.figure()
   fig.suptitle("Violin Plots of COVID Cases and Land Temperatures")
   ax1 = plt.subplot(2,1,1)
   ax1.set_xlabel("Cases per Million")
   sns.violinplot(data=covidtotals.total_cases_mill, color="lightblue", orient="h")
   ax1.set_yticklabels([])
   ax2 = plt.subplot(2,1,2)
   ax2.set_xlabel("Average Temperature")
   sns.violinplot(data=landtemps.avgtemp, color="wheat",  orient="h")
   ax2.set_yticklabels([])
   plt.tight_layout()
   plt.show()
   

这产生了以下图表：

图 1.7 – COVID 病例和陆地温度的小提琴图

中间带有白色圆点的黑色条形表示四分位数范围，而白色圆点代表中位数。当小提琴图水平时，每个点的高度（即当小提琴图水平时）给出了相对频率。四分位数范围右侧的细黑线表示每百万个案例，平均温度的左右两侧是触须。极端值显示在触须之外的分布部分。

如果我要为数值特征创建一个图表，我会创建一个小提琴图。小提琴图让我能够在一个图表中看到集中趋势、形状和分布。空间不允许在这里展示，但我通常喜欢创建所有连续特征的 violin 图，并将它们保存到 PDF 文件中以供以后参考。

摘要

在本章中，我们探讨了探索数据的一些常用技术。我们学习了如何在需要时检索数据的子集以供我们的分析使用。我们还使用了 pandas 方法来生成关于均值、四分位数范围和偏度等特征的关键统计数据。这让我们对每个特征的分布的中心趋势、分布范围和形状有了更好的理解。这也使我们能够更好地识别异常值。最后，我们使用了 Matplotlib 和 Seaborn 库来创建直方图、箱线图和小提琴图。这为我们提供了关于特征分布的额外见解，例如尾部长度和与正态分布的偏离程度。

可视化是本章讨论的单变量分析工具的绝佳补充。直方图、箱线图和小提琴图显示了每个特征分布的形状和分布范围。从图形上看，它们显示了通过检查一些汇总统计数据可能遗漏的内容，例如分布中的隆起（或隆起）位置以及极端值所在的位置。当我们探索特征与目标之间的双变量和多变量关系时，这些可视化将同样有用，我们将在第二章“检查特征与目标之间的双变量和多变量关系”中进行探讨。

第二章：第二章：检查特征与目标变量之间的双变量和多变量关系

在本章中，我们将探讨可能特征与目标变量之间的相关性。使用交叉表（双向频率）、相关性、散点图和分组箱线图的双变量探索性分析可以揭示建模的关键问题。常见问题包括特征之间高度相关以及特征与目标变量之间的非线性关系。在本章中，我们将使用 pandas 方法进行双变量分析，并使用 Matplotlib 进行可视化。我们还将讨论我们在特征工程和建模方面发现的影响。

我们还将使用多元技术来理解特征之间的关系。这包括依赖一些机器学习算法来识别可能存在问题的观测值。之后，我们将就消除某些观测值以及转换关键特征提供初步建议。

在本章中，我们将涵盖以下主题：

• 在双变量关系中识别异常值和极端值

• 使用散点图查看连续特征之间的双变量关系

• 使用分组箱线图查看连续和分类特征之间的双变量关系

• 使用线性回归识别具有显著影响的数据点

• 使用 K 近邻算法寻找异常值

• 使用隔离森林算法寻找异常值

技术要求

本章将大量依赖 pandas 和 Matplotlib 库，但你不需要对这些库有任何先前的知识。如果你是从科学发行版安装的 Python，例如 Anaconda 或 WinPython，那么这些库可能已经安装好了。我们还将使用 Seaborn 进行一些图形，并使用 statsmodels 库进行一些汇总统计。如果你需要安装任何包，可以从终端窗口或 Windows PowerShell 运行pip install [package name]。本章的代码可以在本书的 GitHub 仓库github.com/PacktPublishing/Data-Cleaning-and-Exploration-with-Machine-Learning中找到。

在双变量关系中识别异常值和极端值

没有一个好的双变量关系感，很难开发出一个可靠的模型。我们不仅关心特定特征与目标变量之间的关系，还关心特征如何一起移动。如果特征高度相关，那么建模它们的独立效应可能变得棘手或没有必要。即使特征只是在某个值域内高度相关，这也可能是一个挑战。

对双变量关系的良好理解对于识别异常值也很重要。一个值可能出乎意料，即使它不是一个极端值。这是因为当第二个特征具有某些值时，某些特征的值可能是不寻常的。当其中一个特征是分类的，而另一个是连续的时，这一点很容易说明。

下面的图表说明了多年来每天鸟类观测的数量，但显示了两个地点不同的分布。一个地点每天的平均观测数为 33，而另一个地点为 52。（这是一个从我的Python 数据清洗食谱中提取的虚构示例。）整体平均数（未显示）为 42。对于每天 58 次观测的值我们应该如何看待？它是异常值吗？这取决于观察的是哪个地点。如果地点 A 一天有 58 次观测，那么 58 将是一个非常高的数字。然而，对于地点 B 来说，58 次观测并不会与该地点的平均值有很大不同：

图 2.1 – 每日鸟类观测

这暗示了一个有用的经验法则：当感兴趣的特征与另一个特征相关时，我们在尝试识别异常值（或任何与该特征相关的建模）时应该考虑这种关系。更精确地表述这一点并将其扩展到两个特征都是连续的情况是有帮助的。如果我们假设特征x和特征y之间存在线性关系，我们可以用熟悉的y = mx + b方程来描述这种关系，其中m是斜率，b是y轴截距。然后，我们可以预期y的值将接近x乘以估计的斜率，加上y轴截距。意外值是那些与这种关系有较大偏差的值，其中y的值远高于或低于根据x的值预测的值。这可以扩展到多个x，或预测变量。

在本节中，我们将学习如何通过考察一个特征与另一个特征之间的关系来识别异常值和意外值。在本章的后续部分，我们将使用多元技术来进一步提高我们的异常值检测。

在本节中，我们将基于各国 COVID-19 病例的数据进行工作。数据集包含每百万人口中的病例和死亡人数。我们将把这两个列都视为可能的靶标。它还包含每个国家的人口统计数据，例如人均 GDP、中位数年龄和糖尿病患病率。让我们开始吧：

注意

Our World in Data 在ourworldindata.org/coronavirus-source-data提供 COVID-19 公共使用数据。本节所使用的数据集是在 2021 年 7 月 9 日下载的。数据中包含的列比我包含的要多。我根据国家创建了region列。

1. 让我们从加载 COVID-19 数据集并查看其结构开始。我们还将导入 Matplotlib 和 Seaborn 库，因为我们将会进行一些可视化：

   import pandas as pd
   import matplotlib.pyplot as plt
   import seaborn as sns
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   covidtotals.info()
   <class 'pandas.core.frame.DataFrame'>
   Index: 221 entries, AFG to ZWE
   Data columns (total 16 columns):
    #  Column                  Non-Null Count     Dtype
   --  --------                ---------------    -------
    0  lastdate                221    non-null    object 
    1  location                221    non-null    object 
    2  total_cases             192    non-null    float64
    3  total_deaths            185    non-null    float64
    4  total_cases_mill        192    non-null    float64
    5  total_deaths_mill       185    non-null    float64
    6  population              221    non-null    float64
    7  population_density      206    non-null    float64
    8  median_age              190    non-null    float64
    9  gdp_per_capita          193    non-null    float64
   10  aged_65_older           188    non-null    float64
   11  total_tests_thous        13    non-null    float64
   12  life_expectancy         217    non-null    float64
   13  hospital_beds_thous     170    non-null    float64
   14  diabetes_prevalence     200    non-null    float64
   15  region                  221    non-null    object
   dtypes: float64(13), object(3)
   memory usage: 29.4+ KB
   

2. 在我们检查双变量关系的过程中，从相关性开始是一个很好的起点。首先，让我们创建一个包含一些关键特征的 DataFrame：

   totvars = ['location','total_cases_mill', 
     'total_deaths_mill']
   demovars = ['population_density','aged_65_older', 
     'gdp_per_capita','life_expectancy', 
     'diabetes_prevalence']
   covidkeys = covidtotals.loc[:, totvars + demovars]
   

3. 现在，我们可以获取这些特征的皮尔逊相关矩阵。病例和每百万死亡之间的正相关系数为 0.71。65 岁或以上的人口百分比与病例和死亡都呈正相关，两者均为 0.53。预期寿命也与每百万病例高度相关。似乎至少有一些与国内生产总值（GDP）每人相关的病例：

   corrmatrix = covidkeys.corr(method="pearson")
   corrmatrix
                     total_cases_mill  total_deaths_mill\
   total_cases_mill              1.00               0.71
   total_deaths_mill             0.71               1.00
   population_density            0.04              -0.03
   aged_65_older                 0.53               0.53
   gdp_per_capita                0.46               0.22
   life_expectancy               0.57               0.46
   diabetes_prevalence           0.02              -0.01
                   population_density  aged_65_older  gdp_per_capita\
   total_cases_mill         0.04          0.53       0.46
   total_deaths_mill       -0.03          0.53       0.22
   population_density       1.00          0.06       0.41
   aged_65_older            0.06          1.00       0.49
   gdp_per_capita           0.41          0.49       1.00
   life_expectancy          0.23          0.73       0.68
   diabetes_prevalence      0.01         -0.06       0.12
                     life_expectancy  diabetes_prevalence
   total_cases_mill        0.57           0.02  
   total_deaths_mill       0.46          -0.01  
   population_density      0.23           0.01  
   aged_65_older           0.73          -0.06  
   gdp_per_capita          0.68           0.12  
   life_expectancy         1.00           0.19  
   diabetes_prevalence     0.19           1.00
   

值得注意的是，可能特征之间的相关性，例如，预期寿命和人均 GDP（0.68）以及预期寿命和 65 岁或以上人群（0.73）之间的相关性。

1. 将相关矩阵作为热图查看可能会有所帮助。这可以通过将相关矩阵传递给 Seaborn 的heatmap方法来完成：

   sns.heatmap(corrmatrix, xticklabels =
     corrmatrix.columns, yticklabels=corrmatrix.columns, 
     cmap="coolwarm")
   plt.title('Heat Map of Correlation Matrix')
   plt.tight_layout()
   plt.show()
   

这会创建以下图表：

图 2.2 – COVID 数据的热图，最强的相关性用红色和桃色表示

我们需要关注用较暖色调显示的细胞——在这种情况下，主要是桃色。我发现使用热图有助于我在建模时记住相关性。

注意

本书包含的所有彩色图像都可以下载。请查看本书的前言以获取相应的链接。

1. 让我们更仔细地看看每百万总病例数和每百万死亡数之间的关系。比仅仅通过相关系数更好地理解这一点的一种方法是比较每个的高值和低值，看看它们是如何一起变化的。在下面的代码中，我们使用qcut方法为病例创建一个具有五个值的分类特征，这些值从非常低到非常高分布得相对均匀：我们为死亡也做了同样的事情：

   covidkeys['total_cases_q'] = \
     pd.qcut(covidkeys['total_cases_mill'],
     labels=['very low','low','medium','high',
     'very high'], q=5, precision=0)
   covidkeys['total_deaths_q'] = \
     pd.qcut(covidkeys['total_deaths_mill'],
     labels=['very low','low','medium','high',
     'very high'], q=5, precision=0)
   

2. 我们可以使用crosstab函数查看每个病例五分位数和死亡五分位数的国家数量。正如我们所预期的，大多数国家都在对角线上。有 27 个国家病例和死亡都非常低，有 25 个国家病例和死亡都非常高。有趣的是那些不在对角线上的计数，例如，四个病例非常高但死亡中等的国家，以及一个病例中等但死亡非常高的国家。让我们也看看我们特征的均值，这样我们以后可以参考它们：

   pd.crosstab(covidkeys.total_cases_q, 
     covidkeys.total_deaths_q)
   total_deaths_q  very low  low  medium  high  very high
   total_cases_q                                         
   very low              27    7       0     0          0
   low                    9   24       4     0          0
   medium                 1    6      23     6          1
   high                   0    0       6    21         11
   very high              0    0       4    10         25
   covidkeys.mean()
   total_cases_mill         36,649
   total_deaths_mill           683
   population_density          453
   aged_65_older                 9
   gdp_per_capita           19,141
   life_expectancy              73
   diabetes_prevalence           8
   

3. 让我们仔细看看远离对角线的国家。四个国家——塞浦路斯、科威特、马尔代夫和卡塔尔——的每百万死亡数低于平均水平，但每百万病例数则远高于平均水平。有趣的是，这四个国家在人口规模上都非常小；其中三个国家的人口密度远低于平均的 453；再次，其中三个国家的 65 岁或以上人口比例远低于平均水平：

   covidtotals.loc[(covidkeys.total_cases_q=="very high")
     & (covidkeys.total_deaths_q=="medium")].T
   iso_code            CYP        KWT        MDV    QAT
   lastdate            2021-07-07 2021-07-07 2021-07-07  2021-07-07
   location            Cyprus    Kuwait   Maldives  Qatar
   total_cases         80,588    369,227  74,724  222,918
   total_deaths        380       2,059    213     596
   total_cases_mill    90,752    86,459   138,239 77,374
   total_deaths_mill   428       482      394     207
   population         888,005 4,270,563 540,542 2,881,060
   population_density  128       232      1,454    227
   median_age          37        34       31       32
   gdp_per_capita      32,415    65,531   15,184   116,936
   aged_65_older       13        2        4        1
   total_tests_thous   NaN       NaN      NaN     NaN
   life_expectancy     81        75       79      80
   hospital_beds_thous 3         2        NaN     1
   diabetes_prevalence 9         16       9       17
   region              Eastern   West     South   West  
                       Europe    Asia     Asia    Asia
   

4. 让我们仔细看看那些根据病例数预期死亡数较多的国家。对于墨西哥来说，每百万病例数远低于平均水平，而每百万死亡数则相当高于平均水平：

   covidtotals.loc[(covidkeys. total_cases_q=="medium")
     & (covidkeys.total_deaths_q=="very high")].T
   iso_code                         MEX
   lastdate                  2021-07-07
   location                      Mexico
   total_cases                2,558,369
   total_deaths                 234,192
   total_cases_mill              19,843
   total_deaths_mill              1,816
   population               128,932,753
   population_density                66
   median_age                        29
   gdp_per_capita                17,336
   aged_65_older                      7
   total_tests_thous                NaN
   life_expectancy                   75
   hospital_beds_thous                1
   diabetes_prevalence               13
   region                 North America
   

当我们想要了解数据集中双变量关系时，相关系数和热图是一个好的起点。然而，仅用相关系数来可视化连续变量之间的关系可能很困难。这尤其适用于关系不是线性的情况——也就是说，当它基于特征的取值范围而变化时。我们通常可以通过散点图来提高我们对两个特征之间关系的理解。我们将在下一节中这样做。

使用散点图查看连续特征之间的双变量关系

在本节中，我们将学习如何获取数据的散点图。

我们可以使用散点图来获取比仅通过相关系数所能检测到的更完整的两个特征之间的关系图景。这在数据关系随数据范围变化而变化时尤其有用。在本节中，我们将创建一些与上一节中考察的相同特征的散点图。让我们开始吧：

1. 通过数据点绘制回归线是有帮助的。我们可以使用 Seaborn 的regplot方法来做这件事。让我们再次加载 COVID-19 数据，以及 Matplotlib 和 Seaborn 库，并生成total_cases_mill与total_deaths_mill的散点图：

   import pandas as pd
   import numpy as np
   import matplotlib.pyplot as plt
   import seaborn as sns
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   ax = sns.regplot(x="total_cases_mill",
     y="total_deaths_mill", data=covidtotals)
   ax.set(xlabel="Cases Per Million", ylabel="Deaths Per 
     Million", title="Total COVID Cases and Deaths by 
     Country")
   plt.show()
   

这会产生以下图表：

图 2.3 – 按国家划分的 COVID 病例和死亡总数

回归线是对每百万病例和每百万死亡之间关系的估计。线的斜率表明，我们可以预期每百万死亡数随着每百万病例数增加 1 个单位而增加多少。那些在散点图上显著高于回归线的点应该被更仔细地检查。

1. 死亡率每百万接近 6,000 例，病例率每百万低于 75,000 例的国家显然是异常值。让我们仔细看看：

   covidtotals.loc[(covidtotals.total_cases_mill<75000) \
     & (covidtotals.total_deaths_mill>5500)].T
   iso_code                           PER
   lastdate                    2021-07-07
   location                          Peru
   total_cases                  2,071,637
   total_deaths                   193,743
   total_cases_mill                62,830
   total_deaths_mill                5,876
   population                  32,971,846
   population_density                  25
   median_age                          29
   gdp_per_capita                  12,237
   aged_65_older                        7
   total_tests_thous                  NaN
   life_expectancy                     77
   hospital_beds_thous                  2
   diabetes_prevalence                  6
   region                   South America
   

在这里，我们可以看到异常值国家是秘鲁。秘鲁的确每百万人口病例数高于平均水平，但其每百万死亡人数仍然远高于按病例数预期的数值。如果我们画一条垂直于x轴的线，在 62,830 处，我们可以看到它在大约每百万 1,000 人死亡处与回归线相交，这比秘鲁的 5,876 要少得多。在秘鲁的数据中，除了病例数之外，还有两个非常不同于数据集平均值的数值也引人注目，那就是人口密度和人均 GDP，这两个数值都显著低于平均水平。在这里，我们没有任何特征能帮助我们解释秘鲁的高死亡率。

注意

当创建散点图时，通常会将一个特征或预测变量放在x轴上，将目标变量放在y轴上。如果画了一条回归线，那么它代表的是预测变量增加 1 个单位时目标变量的增加。但散点图也可以用来检查两个预测变量或两个可能的目标变量之间的关系。

回顾我们在第一章中定义的异常值，即检查特征和目标变量的分布，可以认为秘鲁是一个异常值。但我们还需要做更多的工作才能得出这个结论。秘鲁并不是唯一一个在散点图上点远高于或低于回归线的国家。通常来说，调查这些点中的许多点是件好事。让我们来看看：

1. 创建包含大多数关键连续特征的散点图可以帮助我们识别其他可能的异常值，并更好地可视化我们在本章第一部分观察到的相关性。让我们创建 65 岁及以上人口和人均 GDP 与每百万总病例数的散点图：

   fig, axes = plt.subplots(1,2, sharey=True)
   sns.regplot(x=covidtotals.aged_65_older, 
     y=covidtotals.total_cases_mill, ax=axes[0])
   sns.regplot(x=covidtotals.gdp_per_capita, 
     y=covidtotals.total_cases_mill, ax=axes[1])
   axes[0].set_xlabel("Aged 65 or Older")
   axes[0].set_ylabel("Cases Per Million")
   axes[1].set_xlabel("GDP Per Capita")
   axes[1].set_ylabel("")
   plt.suptitle("Age 65 Plus and GDP with Cases Per 
     Million")
   plt.tight_layout()
   fig.subplots_adjust(top=0.92)
   plt.show()
   

这会产生以下图表：

图 2.4 – 65 岁及以上人口和人均 GDP 与每百万总病例数

这些散点图显示，一些每百万病例数非常高的国家，其数值接近我们根据人口年龄或 GDP 预期的数值。这些是极端值，但并不一定是我们定义的异常值。

使用散点图可以展示两个特征与目标之间的关系，所有这些都在一个图形中。让我们回到上一章中我们处理过的陆地温度数据，来探讨这一点。

数据注意

陆地温度数据集包含了 2019 年从全球超过 12,000 个气象站读取的平均温度（以摄氏度为单位），尽管大多数站点位于美国。该数据集是从全球历史气候学网络综合数据库中检索的。它已由美国国家海洋和大气管理局在www.ncdc.noaa.gov/data-access/land-based-station-data/land-based-datasets/global-historical-climatology-network-monthly-version-4上提供给公众使用。

1. 我们预计气象站的平均温度会受到纬度和海拔的影响。假设我们之前的分析表明，海拔对温度的影响直到大约 1,000 米才开始显著。我们可以将landtemps数据框分为低海拔和高海拔站，以 1,000 米为阈值。在下面的代码中，我们可以看到这给我们带来了 9,538 个低海拔站，平均温度为 12.16 摄氏度，以及 2,557 个高海拔站，平均温度为 7.58：

   landtemps = pd.read_csv("data/landtemps2019avgs.csv")
   low, high = landtemps.loc[landtemps.elevation<=1000],
     landtemps.loc[landtemps.elevation>1000]
   low.shape[0], low.avgtemp.mean()
   (9538, 12.161417937651676)
   high.shape[0], high.avgtemp.mean()
   (2557, 7.58321486951755)
   

2. 现在，我们可以在一个散点图中可视化海拔和纬度与温度之间的关系：

   plt.scatter(x="latabs", y="avgtemp", c="blue",
     data=low)
   plt.scatter(x="latabs", y="avgtemp", c="red", 
     data=high)
   plt.legend(('low elevation', 'high elevation'))
   plt.xlabel("Latitude (N or S)")
   plt.ylabel("Average Temperature (Celsius)")
   plt.title("Latitude and Average Temperature in 2019")
   plt.show()
   

这产生了以下散点图：

图 2.5 – 2019 年纬度和平均温度

在这里，我们可以看到随着赤道距离（以纬度衡量）的增加，温度逐渐降低。我们还可以看到高海拔气象站（那些带有红色圆点的）通常位于低海拔站下方——也就是说，在相似的纬度上，它们的温度较低。

1. 似乎高海拔和低海拔站点之间的斜率也有至少一些差异。随着纬度的增加，高海拔站点的温度似乎下降得更快。我们可以在散点图上绘制两条回归线——一条用于高海拔站点，一条用于低海拔站点——以获得更清晰的图像。为了简化代码，让我们为低海拔和高海拔站点创建一个分类特征，elevation_group：

   landtemps['elevation_group'] = 
     np.where(landtemps.elevation<=1000,'low','high')
   sns.lmplot(x="latabs", y="avgtemp", 
     hue="elevation_group", palette=dict(low="blue", 
     high="red"), legend_out=False, data=landtemps)
   plt.xlabel("Latitude (N or S)")
   plt.ylabel("Average Temperature")
   plt.legend(('low elevation', 'high elevation'), 
     loc='lower left')
   plt.yticks(np.arange(-60, 40, step=20))
   plt.title("Latitude and Average Temperature in 2019")
   plt.tight_layout()
   plt.show()
   

这产生了以下图形：

图 2.6 – 2019 年纬度和平均温度及回归线

在这里，我们可以看到高海拔站点的更陡峭的负斜率。

1. 如果我们想看到一个包含两个连续特征和一个连续目标的散点图，而不是像上一个例子中那样将一个特征强制二分，我们可以利用 Matplotlib 的 3D 功能：

   fig = plt.figure()
   plt.suptitle("Latitude, Temperature, and Elevation in 
     2019")
   ax = plt.axes(projection='3d')
   ax.set_xlabel("Elevation")
   ax.set_ylabel("Latitude")
   ax.set_zlabel("Avg Temp")
   ax.scatter3D(landtemps.elevation, landtemps.latabs, 
     landtemps.avgtemp)
   plt.show()
   

这产生了以下三维散点图：

图 2.7 – 2019 年的纬度、温度和海拔

散点图是揭示连续特征之间关系的首选可视化工具。我们比仅通过相关系数更能感受到这些关系。然而，如果我们正在检查连续特征和分类特征之间的关系，我们需要一个非常不同的可视化。分组箱线图在这种情况下非常有用。在下一节中，我们将学习如何使用 Matplotlib 创建分组箱线图。

使用分组箱线图查看连续和分类特征之间的双变量关系

分组箱线图是一种被低估的视觉化工具。当我们检查连续和分类特征之间的关系时，它们非常有帮助，因为它们显示了连续特征的分布如何因分类特征的不同值而变化。

我们可以通过回到上一章中我们使用过的国家纵向调查（NLS）数据来探索这一点。NLS 对每位调查受访者有一个观测值，但收集有关教育和就业的年度数据（每年数据被捕获在不同的列中）。

数据备注

如第一章中所述，检查特征和目标分布，青年国家纵向调查由美国劳工统计局进行。可以从相应的存储库下载 SPSS、Stata 和 SAS 的单独文件。NLS 数据可以从www.nlsinfo.org/investigator/pages/search下载。

按照以下步骤创建分组箱线图：

1. 在 NLS DataFrame 的许多列中，有highestdegree和weeksworked17，分别代表受访者获得的最高学位和 2017 年该人工作的周数。让我们看看获得不同学位的人的工作周数的分布。首先，我们必须定义一个函数gettots来获取我们想要的描述性统计。然后，我们必须使用apply将一个groupby系列对象groupby(['highestdegree'])['weeksworked17']传递给该函数：

   nls97 = pd.read_csv("data/nls97.csv")
   nls97.set_index("personid", inplace=True)
   def gettots(x):
     out = {}
     out['min'] = x.min()
     out['qr1'] = x.quantile(0.25)
     out['med'] = x.median()
     out['qr3'] = x.quantile(0.75)
     out['max'] = x.max()
     out['count'] = x.count()
     return pd.Series(out)
   nls97.groupby(['highestdegree'])['weeksworked17'].\
     apply(gettots).unstack()
                      min   qr1   med   qr3   max   count
   highestdegree                                  
   0\. None              0     0    40    52    52     510
   1\. GED               0     8    47    52    52     848
   2\. High School       0    31    49    52    52   2,665
   3\. Associates        0    42    49    52    52     593
   4\. Bachelors         0    45    50    52    52   1,342
   5\. Masters           0    46    50    52    52     538
   6\. PhD               0    46    50    52    52      51
   7\. Professional      0    47    50    52    52      97
   

在这里，我们可以看到没有高中学位的人的工作周数的分布与拥有学士学位或更高学位的人的分布有多么不同。对于没有学位的人，超过 25%的人在一年中工作了 0 周。对于拥有学士学位的人，即使在第 25 百分位的人，一年中也工作了 45 周。四分位距覆盖了没有学位的个人的整个分布（0 到 52），但对于拥有学士学位的个人的分布（45 到 52）只覆盖了一小部分。

我们还应该注意highestdegree的类别不平衡。在硕士学位之后，计数变得相当小，而高中学位的计数几乎是下一个最大组的两倍。在我们使用这些数据进行任何建模之前，我们可能需要合并一些类别。

1. 分组箱线图使分布差异更加明显。让我们用相同的数据创建一些。我们将使用 Seaborn 来绘制此图：

   import seaborn as sns
   myplt = sns.boxplot(x='highestdegree', 
     y= 'weeksworked17' , data=nls97,
     order=sorted(nls97.highestdegree.dropna().unique()))
   myplt.set_title("Boxplots of Weeks Worked by Highest 
     Degree")
   myplt.set_xlabel('Highest Degree Attained')
   myplt.set_ylabel('Weeks Worked 2017')
   myplt.set_xticklabels(myplt.get_xticklabels(),
     rotation=60, horizontalalignment='right')
   plt.tight_layout()
   plt.show()
   

这产生了以下图表：

图 2.8 – 按最高学位划分的每周工作箱线图

分组箱线图说明了按学位获得的每周工作时间之间的四分位距的显著差异。在副学士学位水平（美国两年制学院学位）或以上，有低于须股的值，用点表示。在副学士学位以下，箱线图没有识别出任何异常值或极端值。例如，对于没有学位的人来说，0 周工作值不是极端值，但对于拥有副学士学位或更高学位的人来说，则是。

1. 我们还可以使用分组箱线图来展示 COVID-19 病例的分布如何因地区而异。让我们也添加一个 swarmplot 来查看数据点，因为它们并不多：

   sns.boxplot(x='total_cases_mill', y='region',
     data=covidtotals)
   sns.swarmplot(y="region", x="total_cases_mill",
     data=covidtotals, size=1.5, color=".3", linewidth=0)
   plt.title("Boxplots of Total Cases Per Million by 
     Region")
   plt.xlabel("Cases Per Million")
   plt.ylabel("Region")
   plt.tight_layout()
   plt.show()
   

这产生了以下图表：

图 2.9 – 按地区划分的每百万总病例箱线图

这些分组箱线图显示了每百万个案例的中位数如何因地区而异，从东非和东亚的低端到东欧和西欧的高端。东亚的极端高值低于西欧的第一四分位数。鉴于大多数地区的计数（国家数量）相当小，我们可能应该避免从那以后得出太多结论。

到目前为止，在本章中，我们主要关注特征之间的双变量关系，以及特征与目标之间的那些关系。我们生成的统计和可视化将指导我们将要进行的建模。我们已经开始对可能的特征、它们对目标的影响以及某些特征的分布如何随着另一个特征的值而变化有所了解。

我们将在本章剩余部分探索多元关系。在我们开始建模之前，我们希望对多个特征如何一起移动有所了解。一旦包含其他特征，某些特征是否不再重要？哪些观测值对我们的参数估计的影响大于其他观测值，这对模型拟合有什么影响？同样，哪些观测值与其他观测值不同，因为它们要么具有无效值，要么似乎在捕捉与其他观测值完全不同的现象？我们将在接下来的三个部分中开始回答这些问题。尽管我们不会在构建模型之前得到任何明确的答案，但我们可以通过预测来开始做出困难的建模决策。

使用线性回归来识别具有显著影响的数据点

发现一些观测值对我们的模型、参数估计和预测有出奇高的影响力并不罕见。这可能是或可能不是所希望的。如果这些观测值反映了与数据中的其他部分不同的社会或自然过程，那么具有显著影响力的观测值可能是有害的。例如，假设我们有一个关于飞行动物的飞行距离的数据集，这些动物几乎都是鸟类，除了关于帝王蝶的数据。如果我们使用翅膀结构作为预测迁移距离的因素，那么帝王蝶的数据可能应该被移除。

我们应该回到第一部分中提到的极端值和异常值之间的区别。我们提到，异常值可以被视为具有特征值或特征值之间关系异常的观测值，这些关系在数据中的其他部分无法解释。另一方面，极端值可能反映了特征中的自然且可解释的趋势，或者在整个数据中观察到的特征之间的相同关系。

在具有高影响力的观测值中区分异常值和极端值最为重要。回归分析中影响的一个标准度量是库克距离（Cook's D）。这给出了如果从数据中删除一个观测值，我们的预测将如何变化的度量。

让我们在本节中构建一个相对简单的多元回归模型，使用我们一直在使用的 COVID-19 数据，并为每个观测值生成一个 Cook's D 值：

1. 让我们加载 COVID-19 数据和 Matplotlib 以及 statsmodels 库：

   import pandas as pd
   import matplotlib.pyplot as plt
   import statsmodels.api as sm
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   

2. 现在，让我们来看看每百万人口中的总病例分布以及一些可能的预测因素：

   xvars = ['population_density','aged_65_older',
     'gdp_per_capita','diabetes_prevalence']
   covidtotals[['total_cases_mill'] + xvars].\
     quantile(np.arange(0.0,1.05,0.25)) 
     total_cases_mill  population_density  aged_65_older\
   0.00             8.52                0.14         1.14
   0.25         2,499.75               36.52         3.50 
   0.50        19,525.73               87.25         6.22
   0.75        64,834.62              213.54        13.92
   1.00       181,466.38           20,546.77        27.05
            gdp_per_capita  diabetes_prevalence  
   0.00             661.24                 0.99  
   0.25           3,823.19                 5.34  
   0.50          12,236.71                 7.20  
   0.75          27,216.44                10.61  
   1.00         116,935.60                30.53
   

3. 接下来，让我们定义一个函数getlm，该函数使用 statsmodels 运行线性回归模型并生成影响统计量，包括 Cook's D。此函数接受一个 DataFrame，目标列的名称，以及特征列的名称（通常将目标称为y，将特征称为X）。

我们将使用dropna函数删除任何特征值缺失的观测值。该函数返回估计的系数（包括pvalues），每个观测值的影响度量，以及完整的回归结果（lm）：

def getlm(df, ycolname, xcolnames):
  df = df[[ycolname] + xcolnames].dropna()
  y = df[ycolname]
  X = df[xcolnames]
  X = sm.add_constant(X)
  lm = sm.OLS(y, X).fit()
  influence = lm.get_influence().summary_frame()
  coefficients = pd.DataFrame(zip(['constant'] + 
    xcolnames, lm.params, lm.pvalues),
    columns=['features','params','pvalues'])
  return coefficients, influence, lm 

1. 现在，我们可以调用getlm函数，同时指定每百万人口中的总病例数作为目标，人口密度（每平方英里的人数），65 岁及以上的人口百分比，人均 GDP 和糖尿病患病率作为预测因素。然后，我们可以打印参数估计值。通常，我们希望查看模型的完整摘要，这可以通过lm.summary()生成。这里我们将跳过这一步，以便于理解：

   coefficients, influence, lm = getlm(covidtotals,
     'total_cases_mill', xvars)
   coefficients
   features                        params       pvalues
   0  constant                 -1,076.471         0.870
   1  population_density           -6.906         0.030
   2  aged_65_older             2,713.918         0.000
   3  gdp_per_capita                0.532         0.001
   4  diabetes_prevalence         736.809         0.241
   

人口密度、65 岁及以上人口和 GDP 的系数在 95%的水平上都是显著的（p 值小于 0.05）。人口密度的结果很有趣，因为我们的双变量分析没有揭示人口密度与每百万病例数之间的关系。系数表明，每平方英里人口增加 1 人，每百万病例数减少 6.9 点。更广泛地说，一旦我们控制了 65 岁或以上人口的比例和人均 GDP，人口密度更高的国家每百万人口病例数会更少。这可能是一个偶然的关系，也可能是一个只能通过多元分析检测到的关系。（也可能是因为人口密度与一个对每百万病例数有更大影响的特征高度相关，但这个特征被遗漏在模型之外。这将给我们一个关于人口密度的有偏系数估计。）

1. 我们可以使用我们在调用getlm时创建的影响 DataFrame 来更仔细地查看那些 Cook's D 值高的观测值。定义高 Cook's D 的一种方法是用所有观测值的 Cook's D 平均值的三倍。让我们创建一个包含所有高于该阈值的值的covidtotalsoutliers DataFrame。

有 13 个国家的 Cook's D 值超过了阈值。让我们按 Cook's D 值降序打印出前五个国家。巴林和马尔代夫在病例分布的前四分之一（参见本节之前我们打印的描述性统计）中。它们的人口密度也较高，65 岁或以上的人口比例较低。在其他条件相同的情况下，根据我们的模型关于人口密度与病例之间关系的说法，我们预计这两个国家的每百万病例数会较低。然而，巴林的人均 GDP 非常高，我们的模型告诉我们这与病例数的高发有关。

新加坡和香港的人口密度极高，每百万病例数低于平均水平，尤其是香港。这两个地方单独可能就解释了人口密度系数的方向。它们的人均 GDP 也非常高，这可能会对该系数产生拖累。可能只是我们的模型不应该包括城邦这样的地区：

influencethreshold = 3*influence.cooks_d.mean()
covidtotals = covidtotals.join(influence[['cooks_d']])
covidtotalsoutliers = \
  covidtotals.loc[covidtotals.cooks_d >
  influencethreshold]
covidtotalsoutliers.shape
(13, 17)
covidtotalsoutliers[['location','total_cases_mill', 
  'cooks_d'] + xvars].sort_values(['cooks_d'],
  ascending=False).head()
     location  total_cases_mill   cooks_d    population_density\
iso_code               
BHR  Bahrain        156,793.409     0.230    1,935.907
SGP  Singapore       10,709.116     0.200    7,915.731
HKG  Hong Kong        1,593.307     0.181    7,039.714
JPN  Japan            6,420.871     0.095    347.778
MDV  Maldives       138,239.027     0.069    1,454.433
         aged_65_older  gdp_per_capita  diabetes_prevalence  
iso_code                                                      
BHR              2.372      43,290.705          16.520
SGP             12.922      85,535.383          10.990
HKG             16.303      56,054.920           8.330
JPN             27.049      39,002.223           5.720
MDV              4.120      15,183.616           9.190

1. 那么，让我们看看如果我们移除香港和新加坡，我们的回归模型估计会有什么变化：

   coefficients, influence, lm2 = \
     getlm(covidtotals.drop(['HKG','SGP']),
     'total_cases_mill', xvars)
   coefficients
      features                  params     pvalues
   0  constant              -2,864.219       0.653
   1  population_density        26.989       0.005
   2  aged_65_older          2,669.281       0.000
   3  gdp_per_capita             0.553       0.000
   4  diabetes_prevalence      319.262       0.605
   

模型中的重大变化是人口密度系数现在已改变方向。这表明人口密度估计对异常观测值的敏感性，这些观测值的特征和目标值可能无法推广到其他数据。在这种情况下，这可能适用于像香港和新加坡这样的城邦。

使用线性回归生成影响度量是一种非常有用的技术，并且它有一个优点，那就是它相对容易解释，正如我们所看到的。然而，它确实有一个重要的缺点：它假设特征之间存在线性关系，并且特征是正态分布的。这通常并不成立。我们还需要足够理解数据中的关系，以创建标签，将每百万总病例数识别为目标。这也不总是可能的。在接下来的两个部分中，我们将探讨不做出这些假设的异常值检测机器学习算法。

使用 K 近邻算法寻找异常值

当我们拥有未标记数据时，机器学习工具可以帮助我们识别与其他观测值不同的观测值——也就是说，当没有目标或因变量时。即使选择目标和特征相对简单，也可能有助于在不假设特征之间的关系或特征分布的情况下识别异常值。

尽管我们通常使用K 近邻算法（KNN）来处理标记数据，用于分类或回归问题，我们也可以用它来识别异常观测值。这些观测值的特点是它们的价值与最近邻的价值之间差异最大。KNN 是一个非常受欢迎的算法，因为它直观、对数据的结构假设很少，并且相当灵活。KNN 的主要缺点是它不如许多其他方法高效，尤其是像线性回归这样的参数技术。我们将在第九章，K 近邻、决策树、随机森林和梯度提升回归，以及第十二章，K 近邻在分类中的应用中更详细地讨论这些优点。

我们将使用PyOD，即Python 异常值检测，来识别 COVID-19 数据中与其他国家显著不同的国家。PyOD 可以使用多种算法来识别异常值，包括 KNN。让我们开始吧：

1. 首先，我们需要从 PyOD 库中导入 KNN 模块，以及从sklearn的预处理实用函数中导入StandardScaler。我们还加载了 COVID-19 数据：

   import pandas as pd
   from pyod.models.knn import KNN
   from sklearn.preprocessing import StandardScaler
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   

接下来，我们需要标准化数据，这在我们的特征范围差异很大时很重要，例如从每百万总病例数和人均 GDP 超过 100,000 到糖尿病患病率和 65 岁及以上人口数不到 20。我们可以使用 scikit-learn 的标准缩放器，它将每个特征值转换为 z 分数，如下所示：

这里，是第 j 个特征的第 i 个观测值的值，是特征的均值，而是该特征的标准差。

1. 我们可以使用缩放器仅针对我们将在模型中包含的特征，然后删除一个或多个特征缺失值的所有观测值：

   standardizer = StandardScaler()
   analysisvars =['location', 'total_cases_mill', 
     'total_deaths_mill','population_density',
     'diabetes_prevalence', 'aged_65_older', 
     'gdp_per_capita']
   covidanalysis = 
     covidtotals.loc[:,analysisvars].dropna()
   covidanalysisstand = 
     standardizer.fit_transform(covidanalysis.iloc[:,1:])
   

2. 现在，我们可以运行模型并生成预测和异常分数。首先，我们必须将contamination设置为0.1，以表示我们希望 10%的观测值被识别为异常值。这相当随意，但不是一个坏的开始。在用fit方法运行 KNN 算法后，我们得到预测（异常值为 1，内点值为 0）和异常分数，这是预测的基础（在这种情况下，异常分数最高的前 10%将被预测为 1）：

   clf_name = 'KNN'
   clf = KNN(contamination=0.1)
   clf.fit(covidanalysisstand)
   y_pred = clf.labels_
   y_scores = clf.decision_scores_
   

3. 我们可以将预测和异常分数的两个 NumPy 数组（分别命名为y_pred和y_scores）合并，并将它们转换为 DataFrame 的列。这使得查看异常分数的范围及其相关的预测变得更容易。有 18 个国家被识别为异常值（这是将contamination设置为0.1的结果）。异常值的异常分数在 1.77 到 9.34 之间，而内点值的分数在 0.11 到 1.74 之间：

   pred = pd.DataFrame(zip(y_pred, y_scores), 
     columns=['outlier','scores'], 
     index=covidanalysis.index)
   pred.outlier.value_counts()
   0    156
   1     18
   pred.groupby(['outlier'])[['scores']].\
     agg(['min','median','max'])
                  scores            
              min      median    max
   outlier                   
   0         0.11        0.84   1.74
   1         1.77        2.48   9.34
   

4. 让我们更仔细地看看异常分数最高的国家：

   covidanalysis = covidanalysis.join(pred).\
     loc[:,analysisvars + ['scores']].\
     sort_values(['scores'], ascending=False)
   covidanalysis.head(10)
         location           total_cases_mill   total_deaths_mill\
   iso_code                                                     …  
   SGP   Singapore                 10,709.12         6.15
   HKG   Hong Kong                  1,593.31        28.28
   PER   Peru                      62,830.48     5,876.01
   QAT   Qatar                     77,373.61       206.87
   BHR   Bahrain                  156,793.41       803.37
   LUX   Luxembourg               114,617.81     1,308.36
   BRN   Brunei                       608.02         6.86
   KWT   Kuwait                    86,458.62       482.14
   MDV   Maldives                 138,239.03       394.05
   ARE   United Arab Emirates      65,125.17       186.75
             aged_65_older    gdp_per_capita    scores
   iso_code                                         
   SGP               12.92         85,535.38      9.34
   HKG               16.30         56,054.92      8.03
   PER                7.15         12,236.71      4.37
   QAT                1.31        116,935.60      4.23
   BHR                2.37         43,290.71      3.51
   LUX               14.31         94,277.96      2.73
   BRN                4.59         71,809.25      2.60
   KWT                2.35         65,530.54      2.52
   MDV                4.12         15,183.62      2.51
   ARE                1.14         67,293.48      2.45
   

在上一节中我们确定的具有高影响力的几个地点，包括新加坡、香港、巴林和马尔代夫，它们的异常分数都很高。这进一步证明我们需要对这些国家的数据进行更仔细的审查。可能存在无效数据，或者有理论上的原因导致它们与其他数据差异很大。

与上一节中的线性模型不同，这里没有定义的目标。在这种情况下，我们包括每百万总病例数和每百万总死亡数。秘鲁在这里被识别为异常值，尽管在线性模型中并不是。这部分是因为秘鲁每百万死亡人数非常高，是数据集中的最高值（我们没有在线性回归模型中使用每百万死亡人数）。

1. 注意到日本并不在这个异常值列表中。让我们看看它的异常分数：

   covidanalysis.loc['JPN','scores']
   2.03
   

异常分数是数据集中第 15 高的。与上一节中日本第 4 高的 Cook's D 分数进行比较。

将这些结果与我们可以用 Isolation Forest 进行的类似分析进行比较是很有趣的。我们将在下一节中这样做。

注意

这只是一个简化的例子，展示了我们通常在机器学习项目中采取的方法。这里最重要的遗漏是我们对整个数据集进行分析。正如我们将在第四章“编码、转换和特征缩放”的开始部分讨论的那样，我们希望在早期就将数据分成训练集和测试集。我们将在本书的剩余章节中学习如何在机器学习管道中集成异常值检测。

使用 Isolation Forest 寻找异常值

隔离森林是一种相对较新的机器学习技术，用于识别异常值。它迅速变得流行，部分原因在于其算法被优化以寻找异常值，而不是正常值。它通过连续划分数据来找到异常值，直到数据点被孤立。需要较少划分才能孤立的数据点会获得更高的异常分数。这个过程在系统资源上表现得相当容易。在本节中，我们将学习如何使用它来检测异常的 COVID-19 病例和死亡。

1. 我们可以使用与上一节类似的分析方法，用隔离森林（Isolation Forest）而不是 KNN。让我们首先加载 scikit-learn 的StandardScaler和IsolationForest模块，以及 COVID-19 数据：

   import pandas as pd
   import matplotlib.pyplot as plt
   from sklearn.preprocessing import StandardScaler
   from sklearn.ensemble import IsolationForest
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   

2. 接下来，我们必须标准化数据：

   analysisvars = ['location','total_cases_mill','total_deaths_mill',
     'population_density','aged_65_older','gdp_per_capita']
   standardizer = StandardScaler()
   covidanalysis = covidtotals.loc[:, analysisvars].dropna()
   covidanalysisstand =
     standardizer.fit_transform(covidanalysis.iloc[:, 1:])
   

3. 现在，我们已经准备好运行我们的异常检测模型。n_estimators参数表示要构建多少棵树。将max_features设置为1.0将使用我们所有的特征。predict方法为我们提供异常预测，对于异常值为-1。这是基于异常分数，我们可以通过decision_function获取：

   clf=IsolationForest(n_estimators=50, 
     max_samples='auto', contamination=.1, 
     max_features=1.0)
   clf.fit(covidanalysisstand)
   covidanalysis['anomaly'] = 
     clf.predict(covidanalysisstand)
   covidanalysis['scores'] = 
     clf.decision_function(covidanalysisstand)
   covidanalysis.anomaly.value_counts()
    1    156
   -1     18
   Name: anomaly, dtype: int64
   

   inlier, outlier = 
     covidanalysis.loc[covidanalysis.anomaly==1],\
     covidanalysis.loc[covidanalysis.anomaly==-1]
   outlier[['location','total_cases_mill',
     'total_deaths_mill',
     'scores']].sort_values(['scores']).head(10)
        location   total_cases_mill  total_deaths_mill   scores
   iso_code                                              
   SGP  Singapore      10,709.12          6.15     -0.20
   HKG  Hong Kong       1,593.31         28.28     -0.16
   BHR  Bahrain       156,793.41        803.37     -0.14
   QAT  Qatar          77,373.61        206.87     -0.13
   PER  Peru           62,830.48      5,876.01     -0.12
   LUX  Luxembourg    114,617.81      1,308.36     -0.09
   JPN  Japan           6,420.87        117.40     -0.08
   MDV  Maldives      138,239.03        394.05     -0.07
   CZE  Czechia       155,782.97      2,830.43     -0.06
   MNE  Montenegro    159,844.09      2,577.77     -0.03
   

4. 查看异常值和内属值的可视化很有帮助：

   fig = plt.figure()
   ax = plt.axes(projection='3d')
   ax.set_title('Isolation Forest Anomaly Detection')
   ax.set_zlabel("Cases Per Million (thous.)")
   ax.set_xlabel("GDP Per Capita (thous.)")
   ax.set_ylabel("Aged 65 Plus %")
   ax.scatter3D(inlier.gdp_per_capita/1000,
     inlier.aged_65_older, inlier.total_cases_mill/1000, 
     label="inliers", c="blue")
   ax.scatter3D(outlier.gdp_per_capita/1000,
     outlier.aged_65_older, 
     outlier.total_cases_mill/1000, label="outliers", 
     c="red")
   ax.legend()
   plt.show()
   

这会产生以下图表：

图 2.10 – 隔离森林异常检测 – 人均 GDP 和每百万人口病例数

尽管我们只能通过这种可视化看到三个维度，但图表确实展示了使异常值成为异常值的一些原因。我们预计随着人均 GDP 和 65 岁及以上人口比例的增加，案例数量也会增加。我们可以看到，异常值偏离了预期的模式，其每百万人口病例数明显高于或低于具有相似 GDP 和 65 岁及以上人口值的其他国家。

摘要

在本章中，我们使用了双变量和多变量统计技术和可视化方法，以更好地理解特征之间的双变量关系。我们研究了常见的统计量，例如皮尔逊相关系数。我们还通过可视化来考察双变量关系，当两个特征都是连续变量时使用散点图，当一个特征是分类变量时使用分组箱线图。本章的最后三节探讨了用于考察关系和识别异常值的多变量技术，包括 KNN 和隔离森林等机器学习算法。

现在我们已经对数据的分布有了很好的了解，我们准备开始构建我们的特征，包括填充缺失值和编码、转换以及缩放我们的变量。这将是下一章的重点。

第三章：第三章：识别和修复缺失值

当我说，很少有看似微小且微不足道的事情像缺失值那样具有如此重大的影响时，我想我代表了许多数据科学家。我们花费大量时间担心缺失值，因为它们可能会对我们的分析产生戏剧性和出人意料的效应。这种情况最有可能发生在缺失值不是随机的情况下——也就是说，当它们与一个特征或目标相关时。例如，假设我们正在进行一项关于收入的纵向研究，但受教育程度较低的人更有可能在每年跳过收入问题。这有可能导致我们对教育参数估计的偏差。

当然，识别缺失值甚至不是战斗的一半。接下来，我们需要决定如何处理它们。我们是否要删除具有一个或多个特征缺失值的任何观测值？我们是否基于样本的统计量，如平均值，来插补一个值？或者我们是否基于更具体的统计量，如某个特定类别的平均值，来分配一个值？我们是否认为对于时间序列或纵向数据，最近的时序值可能最有意义？或者我们应该使用更复杂的多元技术来插补值，比如基于线性回归或k 近邻（KNN）？

所有这些问题的答案都是是的。在某个时候，我们将想要使用这些技术中的每一个。当我们做出关于缺失值插补的最终选择时，我们希望能够回答为什么或为什么不选择所有这些可能性。根据具体情况，每个答案都有意义。

在本章中，我们将探讨识别每个特征或目标缺失值的技巧，以及对于大量特征值缺失的观测值的策略。然后，我们将探讨插补值的策略，例如将值设置为整体平均值、给定类别的平均值或前向填充。我们还将检查用于插补缺失值的多元技术，并讨论它们何时适用。

特别地，在本章中，我们将涵盖以下主题：

• 识别缺失值

• 清理缺失值

• 使用回归进行值插补

• 使用 KNN 插补

• 使用随机森林进行插补

技术要求

本章将大量依赖 pandas 和 NumPy 库，但你不需要对这些库有任何先前的知识。如果你从科学发行版，如 Anaconda 或 WinPython 安装了 Python，这些库可能已经安装好了。我们还将使用statsmodels库进行线性回归，以及来自sklearn和missingpy的机器学习算法。如果你需要安装这些包中的任何一个，你可以通过在终端窗口或 Windows PowerShell 中运行pip install [package name]来安装。

识别缺失值

由于识别缺失值是分析师工作流程中如此重要的一个部分，我们使用的任何工具都需要使其能够轻松地定期检查这些值。幸运的是，pandas 使得识别缺失值变得相当简单。

我们将处理 weeksworked16 和 weeksworked17，分别代表 2016 年和 2017 年的工作周数。

注意

我们还将再次处理 COVID-19 数据。这个数据集为每个国家提供了一个观测值，指定了总 COVID-19 病例和死亡人数，以及每个国家的某些人口统计数据。

按照以下步骤来识别我们的缺失值：

1. 让我们从加载 NLS 和 COVID-19 数据开始：

   import pandas as pd
   import numpy as np
   nls97 = pd.read_csv("data/nls97b.csv")
   nls97.set_index("personid", inplace=True)
   covidtotals = pd.read_csv("data/covidtotals.csv")
   covidtotals.set_index("iso_code", inplace=True)
   

2. 接下来，我们计算可能用作特征的列的缺失值数量。我们可以使用 isnull 方法来测试每个特征值是否缺失。如果值缺失，则返回 True，如果不缺失则返回 False。然后，我们可以使用 sum 来计算 True 值的数量，因为 sum 将每个 True 值视为 1，每个 False 值视为 0。我们使用 axis=0 来对每个列的行进行求和：

   covidtotals.shape
   (221, 16)
   demovars = ['population_density','aged_65_older',
     'gdp_per_capita', 'life_expectancy', 
     'diabetes_prevalence']
   covidtotals[demovars].isnull().sum(axis=0)
   population_density        15
   aged_65_older             33
   gdp_per_capita            28
   life_expectancy            4
   diabetes_prevalence       21
   

如我们所见，221 个国家中有 33 个国家的 aged_65_older 有空值。我们几乎对所有国家的 life_expectancy 都有数据。

1. 如果我们想要每行的缺失值数量，我们可以在求和时指定 axis=1。以下代码创建了一个 Series，demovarsmisscnt，包含每个国家人口统计特征的缺失值数量。181 个国家所有特征都有值，11 个国家缺失五个特征中的四个，三个国家缺失所有特征：

   demovarsmisscnt = covidtotals[demovars].isnull().sum(axis=1)
   demovarsmisscnt.value_counts().sort_index()
   0        181
   1        15
   2         6
   3         5
   4        11
   5         3
   dtype: int64
   

2. 让我们看看有四个或更多缺失值的几个国家。这些国家的人口统计数据非常少：

   covidtotals.loc[demovarsmisscnt > = 4, ['location'] +
     demovars].sample(6, random_state=1).T
   iso_code                         FLK   NIU        MSR\
   location            Falkland Islands  Niue  Montserrat
   population_density               NaN   NaN         NaN
   aged_65_older                    NaN   NaN         NaN
   gdp_per_capita                   NaN   NaN         NaN
   life_expectancy                   81    74          74
   diabetes_prevalence              NaN   NaN         NaN
   iso_code                         COK    SYR        GGY
   location                Cook Islands  Syria   Guernsey
   population_density               NaN    NaN        NaN
   aged_65_older                    NaN    NaN        NaN
   gdp_per_capita                   NaN    NaN        NaN
   life_expectancy                   76      7        NaN
   diabetes_prevalence              NaN    NaN        NaN
   

3. 让我们也检查一下总病例和死亡病例的缺失值。29 个国家在每百万人口中的病例数有缺失值，36 个国家每百万死亡人数有缺失值：

   totvars = 
     ['location','total_cases_mill','total_deaths_mill']
   covidtotals[totvars].isnull().sum(axis=0)
   location                0
   total_cases_mill       29
   total_deaths_mill      36
   dtype: int64
   

4. 我们还应该了解哪些国家同时缺失这两个数据。29 个国家同时缺失病例和死亡数据，而我们只有 185 个国家同时拥有这两个数据：

   totvarsmisscnt = 
     covidtotals[totvars].isnull().sum(axis=1)
   totvarsmisscnt.value_counts().sort_index()
   0        185
   1        7
   2        29
   dtype: int64
   

有时，我们会有逻辑缺失值，需要将其转换为实际缺失值。这发生在数据集设计者使用有效值作为缺失值的代码时。这些值通常是 9、99 或 999 等数值，基于变量的允许数字位数。或者可能是一个更复杂的编码方案，其中存在不同原因导致的缺失值的代码。例如，在 NLS 数据集中，代码揭示了受访者为什么没有回答某个问题的原因：-3 是无效跳过，-4 是有效跳过，而 -5 是非访谈。

1. NLS DataFrame 的最后四列包含受访者母亲和父亲完成的最高学历、家庭收入以及受访者出生时母亲的年龄的数据。让我们从受访者母亲完成的最高的学历开始，检查这些列的逻辑缺失值：

   nlsparents = nls97.iloc[:,-4:]
   nlsparents.shape
   (8984, 4)
   nlsparents.loc[nlsparents.motherhighgrade.between(-5, 
     -1), 'motherhighgrade'].value_counts()
   -3        523
   -4        165
   Name: motherhighgrade, dtype: int64
   

2. 有 523 个无效跳过和 165 个有效跳过。让我们看看这四个特征中至少有一个非响应值的几个个体：

   nlsparents.loc[nlsparents.apply(lambda x: x.between(
     -5,-1)).any(axis=1)]
           motherage  parentincome  fatherhighgrade  motherhighgrade
   personid  
   100284    22            50000            12         -3
   100931    23            60200            -3         13
   101122    25               -4            -3         -3
   101414    27            24656            10         -3
   101526    -3            79500            -4         -4
            ...              ...            ...       ...
   999087    -3           121000            -4         16
   999103    -3            73180            12         -4
   999406    19               -4            17         15
   999698    -3            13000            -4         -4
   999963    29               -4            12         13
   [3831 rows x 4 columns]
   

3. 对于我们的分析，非响应的原因并不重要。让我们只计算每个特征的非响应数量，无论非响应的原因是什么：

   nlsparents.apply(lambda x: x.between(-5,-1).sum())
   motherage              608
   parentincome          2396
   fatherhighgrade       1856
   motherhighgrade        688
   dtype: int64
   

4. 在使用这些特征进行分析之前，我们应该将这些值设置为missing。我们可以使用replace将-5 到-1 之间的所有值设置为missing。当我们检查实际缺失值时，我们得到预期的计数：

   nlsparents.replace(list(range(-5,0)), 
     np.nan, inplace=True)
   nlsparents.isnull().sum()
   motherage            608
   parentincome        2396
   fatherhighgrade     1856
   motherhighgrade      688
   dtype: int64
   

本节展示了识别每个特征的缺失值数量以及具有大量缺失值的观测值的一些非常有用的 pandas 技术。我们还学习了如何查找逻辑缺失值并将它们转换为实际缺失值。接下来，我们将首次探讨清理缺失值。

清理缺失值

在本节中，我们将介绍一些处理缺失值的最直接方法。这包括删除存在缺失值的观测值；将样本的汇总统计量，如平均值，分配给缺失值；以及根据数据适当子集的平均值分配值：

1. 让我们加载 NLS 数据并选择一些教育数据：

   import pandas as pd
   nls97 = pd.read_csv("data/nls97b.csv")
   nls97.set_index("personid", inplace=True)
   schoolrecordlist = 
     ['satverbal','satmath','gpaoverall','gpaenglish',
     'gpamath','gpascience','highestdegree',
     'highestgradecompleted']
   schoolrecord = nls97[schoolrecordlist]
   schoolrecord.shape
   (8984, 8)
   

2. 我们可以使用在前一节中探讨的技术来识别缺失值。schoolrecord.isnull().sum(axis=0)为我们提供了每个特征的缺失值数量。绝大多数观测值在satverbal上存在缺失值，共有 7,578 个，占 8,984 个观测值中的大部分。只有 31 个观测值在highestdegree上存在缺失值：

   schoolrecord.isnull().sum(axis=0)
   satverbal                        7578
   satmath                          7577
   gpaoverall                       2980
   gpaenglish                       3186
   gpamath                          3218
   gpascience                       3300
   highestdegree                      31
   highestgradecompleted            2321
   dtype: int64
   

3. 我们可以创建一个 Series，misscnt，它指定了每个观测值的缺失特征数量，misscnt = schoolrecord.isnull().sum(axis=1)。946 个观测值在教育数据上有七个缺失值，而 11 个观测值的所有八个特征都缺失：

   misscnt = schoolrecord.isnull().sum(axis=1)
   misscnt.value_counts().sort_index()
   0         1087
   1          312
   2         3210
   3         1102
   4          176
   5          101
   6         2039
   7          946
   8           11
   dtype: int64
   

4. 让我们再看看一些具有七个或更多缺失值的观测值。看起来highestdegree通常是唯一存在的特征，这并不令人惊讶，因为我们已经发现highestdegree很少缺失：

   schoolrecord.loc[misscnt>=7].head(4).T
   personid              101705  102061  102648  104627
   satverbal                NaN     NaN     NaN     NaN
   satmath                  NaN     NaN     NaN     NaN
   gpaoverall               NaN     NaN     NaN     NaN
   gpaenglish               NaN     NaN     NaN     NaN
   gpamath                  NaN     NaN     NaN     NaN
   gpascience               NaN     NaN     NaN     NaN
   highestdegree          1.GED  0.None   1.GED  0.None
   highestgradecompleted    NaN     NaN     NaN     NaN
   

5. 让我们删除在八个特征中至少有七个缺失值的观测值。我们可以通过将dropna的thresh参数设置为2来实现这一点。这将删除具有少于两个非缺失值的观测值；也就是说，0 个或 1 个非缺失值。使用dropna后，我们得到预期的观测值数量；即，8,984 - 946 - 11 = 8,027：

   schoolrecord = schoolrecord.dropna(thresh=2)
   schoolrecord.shape
   (8027, 8)
   schoolrecord.isnull().sum(axis=1).value_counts().sort_index()
   0      1087
   1       312
   2      3210
   3      1102
   4       176
   5       101
   6      2039
   dtype: int64
   

gpaoverall有相当数量的缺失值——即 2,980 个——尽管我们有三分之二的观测值是有效的（(8,984 - 2,980)/8,984）。如果我们能很好地填充缺失值，我们可能能够将其作为特征挽救。这比仅仅删除这些观测值更可取。如果我们能避免，我们不想失去这些数据，尤其是如果缺失gpaoverall的个体在其他方面与我们预测相关的方面有所不同。

1. 最直接的方法是将gpaoverall的整体平均值分配给缺失值。以下代码使用 pandas Series 的fillna方法将gpaoverall的所有缺失值分配给 Series 的平均值。fillna的第一个参数是你想要为所有缺失值设置的值——在这种情况下，schoolrecord.gpaoverall.mean()。请注意，我们需要记住将inplace参数设置为True以覆盖现有值：

   schoolrecord.gpaoverall.agg(['mean','std','count'])
   mean         281.84
   std           61.64
   count      6,004.00
   Name: gpaoverall, dtype: float64
   schoolrecord.gpaoverall.fillna(
     schoolrecord.gpaoverall.mean(), inplace=True)
   schoolrecord.gpaoverall.isnull().sum()
   0
   schoolrecord.gpaoverall.agg(['mean','std','count'])
   mean      281.84
   std        53.30
   count   8,027.00
   Name: gpaoverall, dtype: float64
   

平均值没有变化。然而，标准差有显著下降，从 61.6 下降到 53.3。这是使用数据集的平均值填充所有缺失值的一个缺点。

1. NLS 数据中wageincome也有相当数量的缺失值。以下代码显示有 3,893 个观测值有缺失值：

   wageincome = nls97.wageincome.copy(deep=True)
   wageincome.isnull().sum()
   copy method, setting deep to True. We wouldn't normally do this but, in this case, we don't want to change the values of wageincome in the underlying DataFrame. We have avoided this here because we will demonstrate a different method of imputing values in the next couple of code blocks.
   

2. 我们与其将wageincome的平均值分配给缺失值，不如使用另一种常见的值填充技术：我们可以将前一个观测值中的最近非缺失值分配给缺失值。fillna的ffill选项会为我们完成这项工作：

   wageincome.fillna(method='ffill', inplace=True)
   wageincome.head().T
   personid
   100061       12,500
   100139      120,000
   100284       58,000
   100292       58,000
   100583       30,000
   Name: wageincome, dtype: float64
   wageincome.isnull().sum()
   0
   wageincome.agg(['mean','std','count'])
   mean      49,549.33
   std       40,014.34
   count      8,984.00
   Name: wageincome, dtype: float64
   

3. 我们可以通过将fillna的method参数设置为bfill来执行向后填充。这会将缺失值设置为最近的后续值。这会产生以下输出：

   wageincome = nls97.wageincome.copy(deep=True)
   wageincome.std()
   40677.69679818673
   wageincome.fillna(method='bfill', inplace=True)
   wageincome.head().T
   personid
   100061       12,500
   100139      120,000
   100284       58,000
   100292       30,000
   100583       30,000
   Name: wageincome, dtype: float64
   wageincome.agg(['mean','std','count'])
   mean    49,419.05
   std     41,111.54
   count    8,984.00
   Name: wageincome, dtype: float64
   

如果缺失值是随机分布的，那么向前填充或向后填充与使用平均值相比有一个优点：它更有可能近似该特征的非缺失值的分布。请注意，标准差并没有大幅下降。

有时候，基于相似观测值的平均值或中位数来填充值是有意义的；比如说，那些对于相关特征具有相同值的观测值。如果我们正在为特征 X1 填充值，而 X1 与 X2 相关，我们可以利用 X1 和 X2 之间的关系来为 X1 填充一个可能比数据集的平均值更有意义的值。当 X2 是分类变量时，这通常很简单。在这种情况下，我们可以为 X2 的关联值填充 X1 的平均值。

1. 在 NLS DataFrame 中，2017 年的工作周数与获得的最高学位相关。以下代码显示了工作周数的平均值如何随着学位获得而变化。工作周数的平均值是 39，但没有学位的人（28.72）要低得多，而有专业学位的人（47.20）要高得多。在这种情况下，将 28.72 分配给未获得学位的个人缺失的工作周数，而不是 39，可能是一个更好的选择：

   nls97.weeksworked17.mean()
   39.01664167916042
   nls97.groupby(['highestdegree'])['weeksworked17'
     ].mean()
   highestdegree
   0\. None                  28.72
   1\. GED                   34.59
   2\. High School           38.15
   3\. Associates            40.44
   4\. Bachelors             43.57
   5\. Masters               45.14
   6\. PhD                   44.31
   7\. Professional          47.20
   Name: weeksworked17, dtype: float64
   

2. 以下代码将缺失工作周数的观测值的平均工作周数分配给具有相同学位获得水平的观测值。我们通过使用groupby创建一个按highestdegree分组的 DataFrame，即groupby(['highestdegree'])['weeksworked17']来实现这一点。然后，我们在apply中使用fillna来填充这些缺失值，用最高学位组的平均值来填充。请注意，我们确保只对最高学位不缺失的观测值进行插补，~nls97.highestdegree.isnull()。对于既缺失最高学位又缺失工作周数的观测值，我们仍然会有缺失值：

   nls97.loc[~nls97.highestdegree.isnull(),
     'weeksworked17imp'] = 
     nls97.loc[ ~nls97.highestdegree.isnull() ].
     groupby(['highestdegree'])['weeksworked17'].
     apply(lambda group: group.fillna(np.mean(group)))
   nls97[['weeksworked17imp','weeksworked17',
     'highestdegree']].head(10)
          weeksworked17imp  weeksworked17   highestdegree
   personid                                              
   100061            48.00         48.00   2\. High School
   100139            52.00         52.00   2\. High School
   100284             0.00          0.00          0\. None
   100292            43.57           NaN     4\. Bachelors
   100583            52.00         52.00   2\. High School
   100833            47.00         47.00   2\. High School
   100931            52.00         52.00    3\. Associates
   101089            52.00         52.00   2\. High School
   101122            38.15           NaN   2\. High School
   101132            44.00         44.00          0\. None
   nls97[['weeksworked17imp','weeksworked17']].\
     agg(['mean','count'])
          weeksworked17imp  weeksworked17
   mean          38.52         39.02
   count      8,953.00      6,670.00
   

这些插补策略——删除缺失值的观测值、分配数据集的平均值或中位数、使用前向或后向填充，或使用相关特征的组均值——对于许多预测分析项目来说都是可行的。当缺失值与目标不相关时，它们效果最好。当这一点成立时，插补值允许我们保留这些观测值的其他信息，而不会对我们的估计产生偏差。

然而，有时情况并非如此，需要更复杂的插补策略。接下来的几节将探讨清理缺失数据的多变量技术。

使用回归插补值

我们在上一个部分结束时，将组均值分配给缺失值，而不是整体样本均值。正如我们讨论的那样，当确定组的特征与具有缺失值的特征相关时，这很有用。使用回归来插补值在概念上与此相似，但我们通常在插补将基于两个或更多特征时使用它。

回归插补用相关特征的回归模型预测的值来替换一个特征的缺失值。这种特定类型的插补被称为确定性回归插补，因为插补值都位于回归线上，没有引入错误或随机性。

这种方法的潜在缺点是它可能会显著降低具有缺失值的特征的方差。我们可以使用随机回归插补来解决这个问题。在本节中，我们将探讨这两种方法。

NLS 数据集中的wageincome特征有几个缺失值。我们可以使用线性回归来插补值。工资收入值是 2016 年报告的 earnings：

1. 让我们从再次加载 NLS 数据开始，并检查 wageincome 以及可能与 wageincome 相关的特征是否存在缺失值。我们还会加载 statsmodels 库。

info 方法告诉我们，对于近 3,000 个观测值，我们缺少 wageincome 的值。其他特征的缺失值较少：

import pandas as pd
import numpy as np
import statsmodels.api as sm
nls97 = pd.read_csv("data/nls97b.csv")
nls97.set_index("personid", inplace=True)
nls97[['wageincome','highestdegree','weeksworked16',
  'parentincome']].info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 8984 entries, 100061 to 999963
Data columns (total 4 columns):
 #  Column               Non-Null Count       Dtype
--  -------               --------------      -----  
 0wageincome            5091 non-null       float64
 1  highestdegree         8953 non-null       object 
 2  weeksworked16         7068 non-null       float64
 3  parentincome8984 non-null       int64
dtypes: float64(2), int64(1), object(1)
memory usage: 350.9+ KB

1. 让我们将 highestdegree 特征转换为数值。这将使我们在本节剩余部分进行的分析更容易：

   nls97['hdegnum'] =
     nls97.highestdegree.str[0:1].astype('float')
   nls97.groupby(['highestdegree','hdegnum']).size() 
   highestdegree    hdegnum
   0\. None                0            953
   1\. GED                 1           1146
   2\. High School         2           3667
   3\. Associates          3            737
   4\. Bachelors           4           1673
   5\. Masters             5            603
   6\. PhD                 6             54
   7\. Professional        7            120
   

2. 正如我们已经发现的，我们需要将 parentincome 的逻辑缺失值替换为实际缺失值。之后，我们可以进行一些相关性分析。每个特征都与 wageincome 有关联，特别是 hdegnum：

   nls97.parentincome.replace(list(range(-5,0)), np.nan,
     inplace=True)
   nls97[['wageincome','hdegnum','weeksworked16', 
     'parentincome']].corr()
               wageincome  hdegnum  weeksworked16  parentincome
   wageincome     1.00      0.40         0.18        0.27
   hdegnum        0.40      1.00         0.24        0.33
   weeksworked16  0.18      0.24         1.00        0.10
   parentincome   0.27      0.33         0.10        1.00
   

3. 我们应该检查对于工资收入存在缺失值的观测值是否在某些重要方面与那些没有缺失值的观测值不同。以下代码显示，这些观测值的学位获得水平、父母收入和工作周数显著较低。这是一个整体均值分配不是最佳选择的情况：

   nls97['missingwageincome'] =
     np.where(nls97.wageincome.isnull(),1,0)
   nls97.groupby(['missingwageincome'])[['hdegnum', 
     'parentincome', 'weeksworked16']].agg(['mean', 
     'count'])
                    hdegnum    parentincome weeksworked16
                    mean count mean   count  mean   count
   missingwageincome                                    
   0                2.76 5072  48,409.13 3803 48.21  5052
   1                1.95 3881  43,565.87 2785 16.36  2016
   

4. 让我们尝试回归插补。首先，让我们进一步清理数据。我们可以用平均值替换缺失的 weeksworked16 和 parentincome 值。我们还应该将 hdegnum 合并到那些获得低于大学学位、拥有大学学位和拥有研究生学位的人群中。我们可以将它们设置为虚拟变量，当它们为 False 或 True 时分别具有 0 或 1 的值。这是在回归分析中处理分类数据的一个经过验证的方法，因为它允许我们根据组别估计不同的 y 截距：

   nls97.weeksworked16.fillna(nls97.weeksworked16.mean(),
     inplace=True)
   nls97.parentincome.fillna(nls97.parentincome.mean(),
     inplace=True)
   nls97['degltcol'] = np.where(nls97.hdegnum<=2,1,0)
   nls97['degcol'] = np.where(nls97.hdegnum.between(3,4),
     1,0)
   nls97['degadv'] = np.where(nls97.hdegnum>4,1,0)
   

   注意

   scikit-learn 具有预处理功能，可以帮助我们完成这类任务。我们将在下一章中介绍其中的一些。

5. 接下来，我们定义一个函数 getlm，用于使用 statsmodels 模块运行线性模型。此函数具有目标或因变量名称的参数 ycolname 和特征或自变量名称的参数 xcolnames。大部分工作是由 statsmodels 的 fit 方法完成的；即 OLS(y, X).fit()：

   def getlm(df, ycolname, xcolnames):
     df = df[[ycolname] + xcolnames].dropna()
     y = df[ycolname]
     X = df[xcolnames]
     X = sm.add_constant(X)
     lm = sm.OLS(y, X).fit()
     coefficients = pd.DataFrame(zip(['constant'] +
       xcolnames,lm.params, lm.pvalues), columns = [
       'features' , 'params','pvalues'])
     return coefficients, lm
   

6. 现在，我们可以使用 getlm 函数来获取参数估计和模型摘要。所有系数在 95% 的置信水平上都是正值且显著的，因为它们的 pvalues 小于 0.05。正如预期的那样，工资收入随着工作周数和父母收入的增加而增加。拥有大学学位与没有大学学位相比，可以增加近 16K 的收入。研究生学位甚至能将收入预测提升更多——比那些低于大学学位的人多近 37K：

   xvars = ['weeksworked16', 'parentincome', 'degcol', 
     'degadv']
   coefficients, lm = getlm(nls97, 'wageincome', xvars)
   coefficients
        features           params           pvalues
   0    constant           7,389.37         0.00
   1    weeksworked16      494.07           0.00
   2    parentincome       0.18             0.00
   3    degcol             15,770.07        0.00
   4    degadv             36,737.84        0.00
   

7. 我们可以使用这个模型来插补工资收入缺失处的值。由于我们的模型包含一个常数，我们需要为预测添加一个常数。我们可以将预测转换为 DataFrame，然后将其与 NLS 数据的其余部分合并。然后，我们可以创建一个新的工资收入特征，wageincomeimp，当工资收入缺失时获取预测值，否则获取原始工资收入值。让我们也看看一些预测，看看它们是否有意义：

   pred = lm.predict(sm.add_constant(nls97[xvars])).
     to_frame().rename(columns= {0: 'pred'})
   nls97 = nls97.join(pred)
   nls97['wageincomeimp'] = 
     np.where(nls97.wageincome.isnull(),
     nls97.pred, nls97.wageincome)
   pd.options.display.float_format = '{:,.0f}'.format
   nls97[['wageincomeimp','wageincome'] + xvars].head(10)
   wageincomeimp  wageincome  weeksworked16  parentincome  degcol  degadv
   personid                                          
   100061     12,500     12,500    48     7,400    0    0
   100139    120,000    120,000    53    57,000    0    0
   100284     58,000     58,000    47    50,000    0    0
   100292     36,547        NaN     4    62,760    1    0
   100583     30,000     30,000    53    18,500    0    0
   100833     39,000     39,000    45    37,000    0    0
   100931     56,000     56,000    53    60,200    1    0
   101089     36,000     36,000    53    32,307    0    0
   101122     35,151        NaN    39    46,362    0    0
   101132          0          0    22     2,470    0    0
   

8. 我们应该查看我们预测的一些汇总统计信息，并将其与实际工资收入值进行比较。插补的工资收入特征的均值低于原始工资收入均值。这并不奇怪，因为我们已经看到，工资收入缺失的个体在正相关特征上的值较低。令人惊讶的是标准差的急剧下降。这是确定性回归插补的一个缺点：

   nls97[['wageincomeimp','wageincome']].
     agg(['count','mean','std'])
          wageincomeimp  wageincome
   count        8,984        5,091
   mean        42,559       49,477
   std         33,406       40,678
   

9. 随机回归插补基于我们模型的残差向预测中添加一个正态分布的错误。我们希望这个错误具有 0 的均值和与残差相同的方差。我们可以使用 NumPy 的正常函数来实现这一点，np.random.normal(0, lm.resid.std(), nls97.shape[0])。lm.resid.std()参数获取我们模型残差的方差。最后一个参数值，nls97.shape[0]，表示要创建多少个值；在这种情况下，我们希望为数据中的每一行创建一个值。

我们可以将这些值与我们的数据合并，然后向我们的预测中添加错误，randomadd：

randomadd = np.random.normal(0, lm.resid.std(),
  nls97.shape[0])
randomadddf = pd.DataFrame(randomadd, 
  columns=['randomadd'], index=nls97.index)
nls97 = nls97.join(randomadddf)
nls97['stochasticpred'] = nls97.pred + nls97.randomadd
nls97['wageincomeimpstoc'] =
  np.where(nls97.wageincome.isnull(),
  nls97.stochasticpred, nls97.wageincome)

1. 这应该会增加方差，但不会对均值产生太大影响。让我们来确认这一点：

   nls97[['wageincomeimpstoc','wageincome']].agg([
     'count','mean','std'])
   
          wageincomeimpstoc  wageincome
   count        8,984           5,091
   mean        42,517          49,477
   std         41,381          40,678
   

这似乎已经起作用了。我们的随机预测与原始工资收入特征的方差几乎相同。

回归插补是利用我们拥有的所有数据为特征插补值的好方法。它通常优于我们在上一节中检查的插补方法，尤其是在缺失值不是随机的情况下。如果我们使用随机回归插补，我们不会人为地降低我们的方差。

在我们开始使用机器学习进行这项工作之前，这是我们用于插补的多变量方法的首选。现在我们有选择使用 KNN 等算法进行这项任务，在某些情况下，这种方法比回归插补具有优势。与回归插补不同，KNN 插补不假设特征之间存在线性关系，或者那些特征是正态分布的。我们将在下一节中探讨 KNN 插补。

使用 KNN 插补

KNN 是一种流行的机器学习技术，因为它直观、易于运行，并且在特征和观测数不是很多时能产生良好的结果。出于同样的原因，它通常用于插补缺失值。正如其名称所暗示的，KNN 识别出与每个观测值特征最相似的 k 个观测值。当它用于插补缺失值时，KNN 使用最近邻来确定要使用的填充值。

我们可以使用 KNN 插补来完成与上一节回归插补相同的插补：

1. 让我们从导入 scikit-learn 的KNNImputer并再次加载 NLS 数据开始：

   import pandas as pd
   import numpy as np
   from sklearn.impute import KNNImputer
   nls97 = pd.read_csv("data/nls97b.csv")
   nls97.set_index("personid", inplace=True)
   

2. 接下来，我们必须准备特征。我们将学位获得合并为三个类别 – 大专以下、大学和大学后学位 – 每个类别由不同的虚拟变量表示。我们还必须将父母收入的逻辑缺失值转换为实际缺失值：

   nls97['hdegnum'] =
     nls97.highestdegree.str[0:1].astype('float')
   nls97['degltcol'] = np.where(nls97.hdegnum<=2,1,0)
   nls97['degcol'] = 
     np.where(nls97.hdegnum.between(3,4),1,0)
   nls97['degadv'] = np.where(nls97.hdegnum>4,1,0)
   nls97.parentincome.replace(list(range(-5,0)), np.nan, 
     inplace=True)
   

3. 让我们创建一个只包含工资收入和一些相关特征的 DataFrame：

   wagedatalist = ['wageincome','weeksworked16', 
     'parentincome','degltcol','degcol','degadv']
   wagedata = nls97[wagedatalist]
   

4. 现在，我们已准备好使用 KNN 插补器的fit_transform方法来获取传递的 DataFrame wagedata中所有缺失值的值。fit_transform返回一个包含wagedata中所有非缺失值以及插补值的 NumPy 数组。我们可以使用与wagedata相同的索引将此数组转换为 DataFrame。这将使在下一步中合并数据变得容易。

   注意

   我们将在整本书中使用这种技术，当我们使用 scikit-learn 的transform和fit_transform方法处理 NumPy 数组时。

我们需要指定用于最近邻数量的值，即 k。我们使用一个经验法则来确定 k – 将观测数的平方根除以 2 (sqrt(N)/2)。在这种情况下，k 为 47：

impKNN = KNNImputer(n_neighbors=47)
newvalues = impKNN.fit_transform(wagedata)
wagedatalistimp = ['wageincomeimp','weeksworked16imp',
  'parentincomeimp','degltcol','degcol','degadv']
wagedataimp = pd.DataFrame(newvalues,
  columns=wagedatalistimp, index=wagedata.index)

1. 现在，我们必须将插补的工资收入和工作周数列与原始 NLS 工资数据合并，并做出一些观察。请注意，使用 KNN 插补时，我们不需要对相关特征的缺失值进行任何预插补（使用回归插补时，我们将工作周数和父母收入设置为数据集的平均值）。但这确实意味着，即使没有很多信息，KNN 插补也会返回一个插补值，例如以下代码块中的personid的101122：

   wagedata = wagedata.join(wagedataimp[['wageincomeimp', 
     'weeksworked16imp']])
   wagedata[['wageincome','weeksworked16','parentincome',
     'degcol','degadv','wageincomeimp']].head(10)
   wageincome  weeksworked16  parentincome degcol  degadv wageincomeimp
   personid         
   100061     12,500    48     7,400    0    0     12,500
   100139    120,000    53    57,000    0    0    120,000
   100284     58,000    47    50,000    0    0     58,000
   100292        NaN     4    62,760    1    0     28,029
   100583     30,000    53    18,500    0    0     30,000
   100833     39,000    45    37,000    0    0     39,000
   100931     56,000    53    60,200    1    0     56,000
   101089     36,000    53    32,307    0    0     36,000
   101122        NaN   NaN       NaN    0    0     33,977
   101132          0     22    2,470    0    0          0
   

2. 让我们来看看原始特征和插补特征的汇总统计。不出所料，插补工资收入的平均值低于原始平均值。正如我们在上一节中发现的，缺失工资收入的观测值在学位获得、工作周数和父母收入方面较低。我们还在工资收入中失去了一些方差：

   wagedata[['wageincome','wageincomeimp']].agg(['count',
     'mean','std'])
   wageincome  wageincomeimp 
   count          5,091        8,984
   mean          49,477        44,781
   std           40,678        32,034
   

KNN 填充时不假设底层数据的分布。在回归填充中，线性回归的标准假设适用——也就是说，特征之间存在线性关系，并且它们是正态分布的。如果情况不是这样，KNN 可能是更好的填充方法。

尽管有这些优点，KNN 填充确实存在局限性。首先，我们必须根据对 k 的一个良好初始假设来调整模型，有时这个假设仅基于我们对数据集大小的了解。KNN 计算成本高，可能不适合非常大的数据集。最后，当要填充的特征与预测特征之间的相关性较弱时，KNN 填充可能表现不佳。作为 KNN 填充的替代方案，随机森林填充可以帮助我们避免 KNN 和回归填充的缺点。我们将在下一节中探讨随机森林填充。

使用随机森林进行填充

随机森林是一种集成学习方法。它使用自助聚合，也称为 bagging，来提高模型精度。它通过重复多次取多棵树的平均值来进行预测，从而得到越来越好的估计。在本节中，我们将使用 MissForest 算法，这是随机森林算法的一个应用，用于寻找缺失值填充。

MissForest 首先填充缺失值的均值或众数（对于连续或分类特征分别适用），然后使用随机森林来预测值。使用这个转换后的数据集，将缺失值替换为初始预测，MissForest 生成新的预测，可能用更好的预测来替换初始预测。MissForest 通常会经历至少四次这样的迭代过程。

运行 MissForest 甚至比使用我们在上一节中使用的 KNN 填充器还要简单。我们将为之前处理过的相同工资收入数据填充值：

1. 让我们先导入 MissForest 模块并加载 NLS 数据：

   import pandas as pd
   import numpy as np
   import sys
   import sklearn.neighbors._base
   sys.modules['sklearn.neighbors.base'] =
     sklearn.neighbors._base
   from missingpy import MissForest
   nls97 = pd.read_csv("data/nls97b.csv")
   nls97.set_index("personid", inplace=True)
   

   注意

   我们需要解决 sklearn.neighbors._base 名称冲突的问题，它可以是 sklearn.neighbors._base 或 sklearn.neighbors.base，具体取决于您使用的 scikit-learn 版本。在撰写本文时，MissForest 使用的是旧名称。

2. 让我们进行与上一节相同的数据清理：

   nls97['hdegnum'] = 
     nls97.highestdegree.str[0:1].astype('float')
   nls97.parentincome.replace(list(range(-5,0)), np.nan,
     inplace=True)
   nls97['degltcol'] = np.where(nls97.hdegnum<=2,1,0)
   nls97['degcol'] = np.where(nls97.hdegnum.between(3,4), 
     1,0)
   nls97['degadv'] = np.where(nls97.hdegnum>4,1,0)
   wagedatalist = ['wageincome','weeksworked16',
     'parentincome','degltcol','degcol','degadv']
   wagedata = nls97[wagedatalist]
   

3. 现在，我们已经准备好运行 MissForest。请注意，这个过程与我们使用 KNN 填充器的过程非常相似：

   imputer = MissForest()
   newvalues = imputer.fit_transform(wagedata)
   wagedatalistimp = ['wageincomeimp','weeksworked16imp', 
     'parentincomeimp','degltcol','degcol','degadv']
   wagedataimp = pd.DataFrame(newvalues, 
     columns=wagedatalistimp , index=wagedata.index)
   Iteration: 0
   Iteration: 1
   Iteration: 2
   Iteration: 3
   

4. 让我们查看一些填充值和一些汇总统计信息。填充值的均值较低，这是预料之中的，因为我们已经了解到缺失值不是随机分布的，拥有较低学历和较少工作周数的人更有可能缺少工资收入的数据：

   wagedata = wagedata.join(wagedataimp[['wageincomeimp', 
     'weeksworked16imp']])
   wagedata[['wageincome','weeksworked16','parentincome',
     'degcol','degadv','wageincomeimp']].head(10)
        wageincome  weeksworked16  parentincome  degcol  degadv  wageincomeimp
   personid                                         
   100061     12,500    48     7,400    0    0     12,500
   100139    120,000    53    57,000    0    0    120,000
   100284     58,000    47    50,000    0    0     58,000
   100292        NaN     4    62,760    1    0     42,065
   100583     30,000    53    18,500    0    0     30,000
   100833     39,000    45    37,000    0    0     39,000
   100931     56,000    53    60,200    1    0     56,000
   101089     36,000     5    32,307    0    0     36,000
   101122        NaN   NaN       NaN    0    0     32,384
   101132          0    22     2,470    0    0          0
   wagedata[['wageincome','wageincomeimp', 
     'weeksworked16','weeksworked16imp']].agg(['count', 
     'mean','std'])
       wageincome  wageincomeimp  weeksworked16  weeksworked16imp
   count    5,091        8,984        7,068         8,984
   mean    49,477       43,140           39            37
   std     40,678       34,725           21            21
   

MissForest使用随机森林算法生成高度准确的预测。与 KNN 不同，它不需要用 k 的初始值进行调优。它也比 KNN 计算成本更低。也许最重要的是，随机森林插补对特征之间低或非常高的相关性不太敏感，尽管在本例中这不是一个问题。

摘要

在本章中，我们探讨了缺失值插补最流行的方法，并讨论了每种方法的优缺点。分配一个总体样本均值通常不是一个好的方法，尤其是在缺失值的观测与其他观测在重要方面不同时。我们还可以显著减少我们的方差。前向填充或后向填充允许我们保持数据中的方差，但它在观测的邻近性有意义时效果最好，例如时间序列或纵向数据。在大多数非平凡情况下，我们将希望使用多元技术，例如回归、KNN 或随机森林插补。

到目前为止，我们还没有涉及到数据泄露的重要问题以及如何创建独立的训练和测试数据集。为了避免数据泄露，我们需要在开始特征工程时独立于测试数据工作训练数据。我们将在下一章更详细地研究特征工程。在那里，我们将编码、转换和缩放特征，同时也要小心地将训练数据和测试数据分开。