数据科学学习指南(三)
原文:
zh.annas-archive.org/md5/9950f82165bee426119f99afc3ab612d译者:飞龙
第十章:探索性数据分析
50 多年前,John Tukey 热情推广了一种与置信区间、假设检验和建模的正式世界不同的数据分析方法。今天,Tukey 的探索性数据分析(EDA)被广泛应用。Tukey 描述 EDA为处理数据的哲学方法:
探索性数据分析是积极的,而不是被动的描述,真正强调发现意外的重要性。
作为一名数据科学家,您将希望在数据生命周期的每个阶段中使用 EDA,从检查数据质量到准备形式建模,再到确认您的模型是否合理。确实,在描述的工作中第九章清理和转换数据过程中,EDA 在指导我们的质量检查和转换中起到了重要作用。
在 EDA 中,我们进入一个发现的过程,不断提出问题,深入探讨未知领域来探索想法。我们使用图表来揭示数据的特征,检查值的分布,并揭示不能从简单的数值摘要中检测到的关系。这种探索包括转换、可视化和总结数据,以建立和确认我们的理解,识别和解决数据可能存在的问题,并为随后的分析提供信息。
EDA 很有趣!但需要实践。学习如何进行 EDA 的最佳方法之一是从他人那里学习,因为他们描述他们在探索数据时的思维过程,我们在本书的例子和案例研究中试图揭示 EDA 思维。
EDA(探索性数据分析)可以提供宝贵的见解,但您需要谨慎对待您所得出的结论。重要的是要认识到,EDA 可能会影响您分析的偏见。EDA 是一个筛选过程和决策过程,可能会影响您后续基于模型的发现的可复制性。有足够的数据,如果您仔细观察,通常可以挖掘出一些完全虚假的有趣内容。
EDA 在科学可重复性危机中的角色已被注意到,数据科学家已经警告不要过度使用它。例如,Gelman 和 Loken指出:
即使在已对给定数据进行了单一分析的情况下,多重比较[数据挖掘]的问题也会出现,因为对变量组合、案例包含和排除、变量转换、在没有主效应的情况下进行交互测试以及分析中的许多其他步骤的不同选择都可能与不同的数据发生。
最佳实践是报告并提供您的 EDA 代码,以便他人了解您所做的选择和您在了解数据过程中采取的路径。
关于可视化的主题分布在三章之间。在第九章中,我们使用图表来辅助我们进行数据整理。那里的图表很基础,发现也很直接。我们没有深入探讨解释和选择图表的问题。在本章中,我们将花更多时间学习如何选择合适的图表并进行解释。由于我们的目标是在执行探索性数据分析时快速生成图表,通常采用绘图函数的默认参数设置。在第十一章,我们将提供制作有效和信息丰富的图表的指南,并提供建议,说明如何使我们的视觉论证清晰和令人信服。
根据Tukey,可视化在探索性数据分析中是核心:
数据中最大的收益来自于惊喜……意外的情况最好通过图片来引起我们的注意。
要制作这些图片,我们需要选择适当类型的图表,我们的选择取决于已收集的数据类型。特征类型与图表选择之间的映射是下一节的主题。然后,我们详细描述如何“读取”图表,要查找的内容以及如何解释所见内容。我们首先讨论单特征图表要查找的内容,然后关注两个特征之间的关系,最后描述三个或更多特征的图表。在介绍了探索性数据分析的可视化工具后,我们提供了执行探索性数据分析的指南,并按照这些指南的步骤进行示例演练。
特征类型
在制作探索性图表或任何图表之前,检查特征(或特征),并确定其特征类型是个好主意。(有时我们将特征称为变量,其类型称为变量类型。)尽管有多种分类特征类型的方法,在本书中,我们考虑三种基本类型。有序和名义数据是分类数据的子类型。分类数据的另一个名称是定性。相反,我们还有定量特征:
名义
代表“命名”类别的特征,其中类别没有自然顺序,称为名义。例如政党隶属(民主党、共和党、绿党、其他)、狗的类型(牧群、猎犬、非体育、体育、梗类、玩具、工作类)和计算机操作系统(Windows、macOS、Linux)。
有序
表示有序类别的测量称为有序测量。有序特征的例子包括 T 恤尺寸(小号、中号、大号)、Likert 量表响应(不同意、中立、同意)和教育水平(高中、大学、研究生院)。重要的是要注意,对于有序特征,例如小号和中号之间的差异不一定等同于中号和大号之间的差异。此外,连续类别之间的差异可能甚至无法量化。考虑餐厅评论中的星级数量及一颗星的含义与两颗星之间的含义。
定量
代表数值测量或数量的数据被称为定量数据。例如,以厘米为单位测量的身高,以美元报告的价格和以公里为单位测量的距离。定量特征可以进一步分为离散,意味着特征只能有几个可能值,以及连续,意味着理论上数量可以测量到任意精度。家庭中兄弟姐妹的数量采用离散的值集(例如 0, 1, 2, …, 8)。相比之下,身高理论上可以报告到任意数量的小数位数,因此我们认为它是连续的。确定数量是离散还是连续没有硬性规定。在某些情况下,这可能是一种判断,而在其他情况下,我们可能希望有意将连续特征视为离散特征。
特征类型与数据存储类型不同。pandas 的每一列 DataFrame 都有自己的存储类型。这些类型可以是整数、浮点数、布尔值、日期时间格式、类别和对象(Python 中以指向字符串的指针存储可变长度的字符串)。我们使用术语特征类型来指代信息的概念性概念,使用术语存储类型来指代计算机中信息的表示。
一个以整数存储的特征可以代表名义数据,字符串可以是定量的(比如 "\$100.00"),在实践中,布尔值通常表示只有两种可能值的名义特征。
注意
pandas 将存储类型称为 dtype,这是数据类型的简称。我们在这里避免使用 数据类型 这个术语,因为它可能会与存储类型和特征类型混淆。
为了确定特征类型,我们经常需要查阅数据集的数据字典或代码簿。数据字典是与数据一起包含的文件,描述了数据表中每一列代表的内容。在以下示例中,我们查看了关于各种狗品种的数据框的列的存储和特征类型,并且发现存储类型通常不能很好地指示字段中包含的信息类型。
示例:狗品种
我们使用 美国肯尼尔俱乐部 (AKC) 的注册狗品种数据来介绍与探索数据分析相关的各种概念。成立于 1884 年的非营利组织 AKC 的宗旨是“推动纯种狗的研究、繁殖、展示、运动和维护”。AKC 组织了诸如全国锦标赛、敏捷邀请赛和服从经典等活动,混种狗在大多数活动中也可以参与。信息美化 网站提供了来自 AKC 关于 172 种狗品种的信息数据集。其可视化作品 最佳展示 包含了多种品种的特征,非常有趣。
AKC 数据集包含多种不同类型的特征,我们已提取了一些显示各种信息类型的特征。这些特征包括品种的名称;其寿命、体重和身高;以及其他信息,例如其适合儿童的程度和学习新技巧所需的重复次数。数据集中的每条记录都是一种狗的品种,提供的信息意在典型代表该品种。让我们将数据读入数据框中:
dogs = pd.read_csv('data/akc.csv')
dogs
| 品种 | 组 | 得分 | 寿命 | ... | 大小 | 重量 | 身高 | 重复次数 | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 边境牧羊犬 | 牧羊 | 3.64 | 12.52 | ... | 中型 | NaN | 51.0 | <5 |
| 1 | 边境梗 | 梗类犬 | 3.61 | 14.00 | ... | 小型 | 6.0 | NaN | 15-25 |
| 2 | 布列塔尼犬 | 狩猎 | 3.54 | 12.92 | ... | 中型 | 16.0 | 48.0 | 5-15 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 169 | 电线狐狸梗 | 梗类犬 | NaN | 13.17 | ... | 小型 | 8.0 | 38.0 | 25-40 |
| 170 | 硬毛指示格里芬猎犬 | 狩猎 | NaN | 8.80 | ... | 中型 | NaN | 56.0 | 25-40 |
| 171 | 墨西哥无毛犬 | 非狩猎 | NaN | NaN | ... | 中型 | NaN | 42.0 | NaN |
172 rows × 12 columns
粗略看一下表格,品种、组和大小似乎是字符串,其他列则是数字。这里显示的数据框摘要提供了索引、名称、非空值计数和每列的dtype:
`dogs``.``info``(``)`
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 172 entries, 0 to 171
Data columns (total 12 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 breed 172 non-null object
1 group 172 non-null object
2 score 87 non-null float64
3 longevity 135 non-null float64
4 ailments 148 non-null float64
5 purchase_price 146 non-null float64
6 grooming 112 non-null float64
7 children 112 non-null float64
8 size 172 non-null object
9 weight 86 non-null float64
10 height 159 non-null float64
11 repetition 132 non-null object
dtypes: float64(8), object(4)
memory usage: 16.2+ KB
这个数据框的几列具有数值计算类型,如float64,这意味着该列可以包含除整数以外的数字。我们还确认pandas将字符串列编码为object类型,而不是string类型。请注意,我们错误地猜测重复次数是定量的。仔细观察数据表后,我们发现重复次数包含范围的字符串值,如"<5"、"15-25"和"25-40",因此这个特征是有序的。
注意
在计算机体系结构中,浮点数或简称“float”指的是可以具有小数部分的数。我们不会在本书中深入探讨计算机体系结构,但在影响术语时,我们会指出它,就像在这种情况下一样。dtype float64 表示该列包含的是在计算机内存中存储时每个占用 64 位空间的十进制数。
另外,pandas 使用了优化的存储类型来存储数值数据,如 float64 或 int64。然而,对于像字符串、字典或集合这样的 Python 对象,它没有优化,因此这些数据类型都存储为 object dtype。这意味着存储类型是不明确的,但在大多数情况下,我们知道 object 列包含字符串或其他某种 Python 类型。
在查看列存储类型时,我们可能会猜测 ailments 和 children 是定量特征,因为它们存储为 float64 dtype。但让我们统计它们的唯一值:
`display_df``(``dogs``[``'``ailments``'``]``.``value_counts``(``)``,` `rows``=``8``)`
ailments
0.0 61
1.0 42
2.0 24
4.0 10
3.0 6
5.0 3
9.0 1
8.0 1
Name: count, dtype: int64
`dogs``[``'``children``'``]``.``value_counts``(``)`
children
1.0 67
2.0 35
3.0 10
Name: count, dtype: int64
ailments 和 children 都只接受几个整数值。children 的值为 3.0 或 ailments 的值为 9.0 代表什么意思?我们需要更多信息才能搞清楚。列名及其在数据框中的存储方式不足以解释。因此,我们需要查阅表 10-1 中显示的数据字典。
表 10-1. AKC 狗品种代码手册
| 特征 | 描述 |
|---|---|
breed | 犬品种,例如边境牧羊犬、达尔马提亚犬、威尔士波音犬 |
group | 美国典狗俱乐部分组(牧羊犬、猎犬、非运动犬、运动犬、梗类犬、玩具犬、工作犬) |
score | AKC 评分 |
longevity | 典型寿命(年) |
ailments | 严重遗传疾病的数量 |
purchase_price | 从 puppyfind.com 平均购买价格 |
grooming | 每次需要的美容频率:1 = 天,2 = 周,3 = 几周 |
children | 对儿童的适合度:1 = 高,2 = 中,3 = 低 |
size | 尺寸:小型、中型、大型 |
weight | 典型体重(公斤) |
height | 肩膀高度(厘米) |
repetition | 理解新命令所需的重复次数:<5,5-15,15-25,25-40,40-80,>80 |
尽管数据字典未明确指定特征类型,但描述足以让我们了解到,children 特征代表品种对儿童的适合程度,1.0 的值对应于“高”适合度。我们还发现 ailments 特征表示该品种狗犬通常具有的严重遗传疾病的数量。根据代码手册,我们将 children 视为分类特征,即使它存储为浮点数,由于低 < 中 < 高,该特征是有序的。由于 ailments 是一个计数,我们将其视为定量(数值)类型,并且对于某些分析,我们进一步定义它为离散型,因为 ailments 可能的取值只有几个。
数据代码表还确认了score、longevity、purchase_price、weight和height特征是定量的。这里的想法是,数值特征具有可以通过差异进行比较的值。可以说吉娃娃的寿命通常比腊肠犬长大约四年(16.5 年对 12.6 年)。另一个检查是是否有意义比较值的比率:一只腊肠犬通常比吉娃娃重大约五倍(11 公斤对 2 公斤)。所有这些定量特征都是连续的;只有ailments是离散的。
对于breed、group、size和repetition,数据字典描述表明这些特征是定性的。每个变量具有不同但常见的特征,值得进一步探索。我们通过检查各个特征的唯一值计数来做到这一点。我们从breed开始:
`dogs``[``'``breed``'``]``.``value_counts``(``)`
breed
Border Collie 1
Great Pyrenees 1
English Foxhound 1
..
Saluki 1
Giant Schnauzer 1
Xoloitzcuintli 1
Name: count, Length: 172, dtype: int64
breed特征有 172 个唯一值,与数据框中的记录数相同,因此我们可以将breed视为数据表的主键。按设计,每个犬种有一条记录,而这个breed特征决定了数据集的粒度。虽然breed也被视为名义特征,但分析它没有真正的意义。我们确实希望确认所有值都是唯一且干净的,但除此之外,我们只会用它来标记绘图中的异常值。
接下来,我们检查group特征:
`dogs``[``'``group``'``]``.``value_counts``(``)`
group
terrier 28
sporting 28
working 27
hound 26
herding 25
toy 19
non-sporting 19
Name: count, dtype: int64
该特征有七个唯一值。由于被标记为“运动型”的犬种与被认为是“玩具型”的犬种在多方面有所不同,这些类别不能轻易归纳为一种顺序。因此,我们将group视为名义特征。名义特征甚至不提供差异方向上的含义。
接下来,我们检查size的唯一值及其计数:
`dogs``[``'``size``'``]``.``value_counts``(``)`
size
medium 60
small 58
large 54
Name: count, dtype: int64
size特征具有自然的顺序:小 < 中 < 大,因此它是序数的。我们不知道“小”类别是如何确定的,但我们知道小型品种在某种意义上比中型品种小,而中型品种又比大型品种小。我们有一个顺序,但概念上差异和比率没有意义。
repetition特征是定量变量的一个示例,已经被折叠成类别,变为序数。数据代码表告诉我们,repetition是新命令需要重复几次才能让狗理解的次数:
`dogs``[``'``repetition``'``]``.``value_counts``(``)`
repetition
25-40 39
15-25 29
40-80 22
5-15 21
80-100 11
<5 10
Name: count, dtype: int64
数值值已被汇总为<5、5-15、15-25、25-40、40-80、80-100,注意这些类别的宽度不同。第一个有 5 次重复,而其他的宽度为 10、15 和 40 次重复。顺序清晰,但从一个类别到下一个的间隔不是相同大小。
现在我们已经再次检查了变量中的值与代码手册中描述的值,我们可以扩充数据字典,包括关于特征类型的额外信息。我们修订后的字典出现在表 10-2 中。
表 10-2. 修订后的 AKC 犬种代码手册
| 特征 | 描述 | 特征类型 | 存储类型 |
|---|---|---|---|
breed | 犬种,例如边境牧羊犬、达尔马提亚犬、威尔士激罗犬 | 主键 | 字符串 |
group | AKC 分组(牧羊犬、猎犬、非运动犬、运动犬、梗类犬、玩具犬、工作犬) | 定性 - 名义 | 字符串 |
| score | AKC score | 定量 | 浮点数 |
longevity | 典型寿命(年) | 定量 | 浮点数 |
ailments | 严重遗传疾病数量(0, 1, …, 9) | 定量 - 离散 | 浮点数 |
purchase_price | 从puppyfind.com平均购买价格 | 定量 | 浮点数 |
grooming | 每隔多久梳理一次:1 = 天、2 = 周、3 = 几周 | 定性 - 序数 | 浮点数 |
children | 适合儿童程度:1 = 高、2 = 中、3 = 低 | 定性 - 序数 | 浮点数 |
size | 尺寸:小、中、大 | 定性 - 序数 | 字符串 |
weight | 典型体重(公斤) | 定量 | 浮点数 |
height | 肩膀高度(厘米) | 定量 | 浮点数 |
repetition | 理解新指令所需的重复次数:<5、5–15、15–25、25–40、40–80、80–100 | 定性 - 序数 | 字符串 |
对 AKC 数据特征类型的更清晰理解有助于我们进行质量检查和转换。我们在第九章中讨论了转换,但还有一些未涉及的额外转换。这些与定性特征的类别有关,我们接下来将对其进行描述。
转换定性特征
无论特征是名义还是序数,我们可能会发现重新标记类别使其更具信息性,合并类别以简化可视化,甚至将数值特征转换为序数以便专注于特定的过渡点都是有用的。我们解释了何时进行每种转换,并给出了示例。
重新标记类别
摘要统计数据,如平均值和中位数,对定量数据是有意义的,但通常对定性数据不是。例如,计算玩具品种的平均价格($687)是有意义的,但是关于儿童的“平均”适合性不是。然而,如果我们要求它,pandas 将乐意计算children列中值的平均值:
`# Don't use this value in actual data analysis!`
`dogs``[``"``children``"``]``.``mean``(``)`
1.4910714285714286
相反,我们希望考虑children的 1、2 和 3 的分布。
注意
存储类型和特征类型之间的主要区别在于,存储类型表示我们可以编写代码来计算哪些操作,而特征类型表示对于数据何种操作是有意义的。
我们可以通过用它们的字符串描述替换数字来转换孩子们。将 1、2 和 3 改为高、中和低,使得很容易识别出孩子们是分类的。使用字符串,我们不会被诱惑计算平均值,类别会与它们的含义相联系,并且绘图的标签默认具有合理的值。例如,让我们只关注玩具品种,并制作适合儿童的酒吧图。首先,我们创建一个新的列,用字符串表示适合性的类别:
`kids` `=` `{``1``:``"``high``"``,` `2``:``"``medium``"``,` `3``:``"``low``"``}`
`dogs` `=` `dogs``.``assign``(``kids``=``dogs``[``'``children``'``]``.``replace``(``kids``)``)`
`dogs`
| 品种 | 分组 | 得分 | 寿命 | ... | 重量 | 身高 | 重复 | 孩子们 | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 边境牧羊犬 | 牧羊犬 | 3.64 | 12.52 | ... | NaN | 51.0 | <5 | 低 |
| 1 | 边境泰瑞尔 | 梗犬 | 3.61 | 14.00 | ... | 6.0 | NaN | 15-25 | 高 |
| 2 | 布列塔尼 | 运动型 | 3.54 | 12.92 | ... | 16.0 | 48.0 | 5-15 | 中等 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 169 | 线毛福克斯梗 | 梗犬 | NaN | 13.17 | ... | 8.0 | 38.0 | 25-40 | NaN |
| 170 | 线毛指示格里芬 | 运动型 | NaN | 8.80 | ... | NaN | 56.0 | 25-40 | NaN |
| 171 | 什么狗品种 | 非运动型 | NaN | NaN | ... | NaN | 42.0 | NaN | NaN |
172 rows × 13 columns
然后我们可以制作玩具品种中每个适合性类别的计数条形图:
`toy_dogs` `=` `dogs``.``query``(``'``group ==` `"``toy``"``'``)``.``groupby``(``'``kids``'``)``.``count``(``)``.``reset_index``(``)`
`px``.``bar``(``toy_dogs``,` `x``=``'``kids``'``,` `y``=``'``breed``'``,` `width``=``350``,` `height``=``250``,`
`category_orders``=``{``"``kids``"``:` `[``"``low``"``,` `"``medium``"``,` `"``high``"``]``}``,`
`labels``=``{``"``kids``"``:` `"``Suitability for children``"``,` `"``breed``"``:` `"``count``"``}``)`
我们并不总是希望用字符串表示分类数据。字符串通常需要更多的存储空间,如果数据集包含许多分类特征,可能会大大增加数据集的大小。
有时候,一个定性特征有许多类别,我们更喜欢对数据进行更高级别的查看,因此我们会合并类别。
折叠类别
让我们创建一个名为play的新列,以表示“目的”是玩耍的狗的组(或不是)。 (这是一个虚构的区分,用于演示目的。)此类别由玩具和非运动品种组成。 新特征play是将特征group折叠的转换:将玩具和非运动组合成一个类别,剩下的类别放在第二个非玩耍类别中。 布尔(bool)存储类型对于指示此特征的存在或不存在非常有用:
`with_play` `=` `dogs``.``assign``(``play``=``(``dogs``[``"``group``"``]` `==` `"``toy``"``)` `|`
`(``dogs``[``"``group``"``]` `==` `"``non-sporting``"``)``)`
将一个两类定性特征表示为布尔值有一些优点。 例如,play的平均值是有意义的,因为它返回True值的比例。 当布尔值用于数字计算时,True变为 1,False变为 0:
`with_play``[``'``play``'``]``.``mean``(``)`
0.22093023255813954
此存储类型为我们提供了计算布尔值的计数和平均值的快捷方式。 在第十五章中,我们将看到它对建模也是一个方便的编码。
有时候,例如当一个离散的定量特征有一个长尾时,我们希望截断更高的值,这将定量特征转化为序数。接下来我们将描述这一点。
将定量转换为序数
最后,我们有时会发现另一种有用的转换是将数值转换为类别。例如,我们可以将ailments中的值合并为类别:0, 1, 2, 3, 4+。换句话说,我们将ailments从定量特征转换为序数特征,映射为 0→0, 1→1, 2→2, 3→3,任何值大于等于 4→4+。我们可能希望进行这种转换,因为几乎没有品种有三种以上的遗传疾病。这种简化可以更清晰和足够用于调查。
注意
截至 2022 年底,pandas也实现了一个设计用于定性数据的category dtype。然而,这种存储类型目前在可视化和建模库中尚未广泛采用,限制了其实用性。因此,我们不将定性变量转换为category dtype。我们预计将来的读者可能会希望在更多库支持它之后使用category dtype。
当我们将定量特征转换为序数时,我们会丢失信息。我们无法回到原来的状态。也就是说,如果我们知道某品种的疾病数为四个或更多,我们无法重新创建实际的数值。当我们合并类别时也是同样的情况。因此,保留原始特征是一个良好的实践。如果需要检查我们的工作或更改类别,我们可以记录和重新创建我们的步骤。
一般来说,特征类型帮助我们确定哪种绘图最合适。接下来我们讨论特征类型与绘图之间的映射。
特征类型的重要性
特征类型指导我们进行数据分析。它们有助于指定我们可以对数据应用的操作、可视化和模型。表 10-3 将特征类型匹配到通常适合的各种绘图类型。变量是定量还是定性通常决定了可以制作的可选绘图集,尽管也有例外情况。决策的其他因素包括观察数量以及特征是否仅取几个不同的值。例如,对于离散定量变量,我们可能制作柱状图而不是直方图。
表 10-3. 将特征类型映射到绘图
| 特征类型 | 维度 | 绘图 |
|---|---|---|
| 定量 | 一个特征 | 地毯图,直方图,密度曲线,箱线图,小提琴图 |
| 定性 | 一个特征 | 柱状图,点图,线图,饼图 |
| 定量 | 两个特征 | 散点图,平滑曲线,等高线图,热力图,分位数-分位数图 |
| 定性 | 两个特征 | 并排柱状图,马赛克图,叠加线条 |
| 混合 | 两个特征 | 叠加密度曲线,并排箱线图,叠加平滑曲线,分位数-分位数图 |
特征类型还帮助我们决定计算哪种摘要统计数据。对于定性数据,我们通常不计算均值或标准差,而是计算每个类别中记录的数量、比例或百分比。对于定量特征,我们计算均值或中位数作为中心测量,以及标准差或四分位间距(第 75 百分位至第 25 百分位)作为扩展测量。除了四分位数外,我们可能还会发现其他百分位数有信息意义。
注意
第n百分位数是使得n%的数据值不超过它的值q。值q可能不唯一,有几种方法可以选择可能的唯一值。有足够的数据时,这些定义之间应该几乎没有区别。
在 Python 中计算百分位数时,我们更喜欢使用:
`np``.``percentile``(``data``,` `method``=``'``lower``'``)`
当探索数据时,我们需要知道如何解释我们的图形显示的形状。接下来的三节将指导这种解释。我们还通过示例介绍了表格 10-3 中列出的许多图表类型。其他类型在第十一章中介绍。
如何选择分布特征
特征的视觉展示可以帮助我们看到观察结果中的模式;它们通常比直接检查数字或字符串本身要好得多。简单的毯子图将每个观察结果定位为沿轴上的“纱线”。“毯子”图在我们有少量观察结果时可能很有用,但是当我们有 100 个值时,很快就难以区分高密度(人口最多)区域,比如说。以下图显示了大约 150 个品种寿命值沿直方图顶部的毯子图:
`px``.``histogram``(``dogs``,` `x``=``"``longevity``"``,` `marginal``=``"``rug``"``,` `nbins``=``20``,`
`histnorm``=``'``percent``'``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``'``longevity``'``:``'``Typical lifespan (yr)``'``}``)`
尽管我们可以看到一个异常大的值大于 16 在毯子图中,但很难比较其他区域的纱线密度。相反,直方图为各种寿命值的观察密度提供了更好的感知。类似地,下图中显示的密度曲线描绘了高低密度区域的图像:
在直方图和密度曲线中,我们可以看到寿命分布是不对称的。在大约 12 年处有一个主要模式,9 到 11 年的范围内有一个肩膀,这意味着虽然 12 年是最常见的寿命,但许多品种的寿命比 12 年短一到三年。我们还看到大约 7 岁的小次要模式,以及一些寿命长达 14 到 16 年的品种。
当解读直方图或密度曲线时,我们会检查分布的对称性和偏斜度;高频区域(模式)的数量、位置和大小;尾部长度(通常与钟形曲线进行比较);未观察到值的间隙;以及异常大或异常值。图 10-1 提供了一个具有这些特征的分布的描述。当我们读取分布时,我们将在图中看到的特征与所测量的数量联系起来。
图 10-1. 示例密度图,根据其形状识别分布的特征
作为另一个例子,犬种中遗传病数量的分布如下直方图所示:
`bins` `=` `[``-``0.5``,` `0.5``,` `1.5``,` `2.5``,` `3.5``,` `9.5``]`
`g` `=` `sns``.``histplot``(``data``=``dogs``,` `x``=``"``ailments``"``,` `bins``=``bins``,` `stat``=``"``density``"``)`
`g``.``set``(``xlabel``=``'``Number of ailments``'``,` `ylabel``=``'``density``'``)``;`
当遗传病数量为 0 时,意味着这个品种没有遗传病,当为 1 时对应一个遗传病,以此类推。从直方图中可以看出,遗传病的分布是单峰的,峰值在 0 处。我们还可以看到,分布向右严重倾斜,右尾长,表明很少的品种有四到九种遗传病。虽然是定量的,但遗传病是离散的,因为只有少数整数值是可能的。因此,我们将区间设定在整数上,例如从 1.5 到 2.5 的区间只包含有两种病的品种。我们还扩宽了最右边的区间。我们将所有有四到九种病的品种归为一组。当区间计数较小时,我们使用更宽的区间进一步平滑分布,因为我们不想过多关注小数字的波动。在这种情况下,没有品种有六或七种病,但有些有四、五、八或九种病。
接下来,我们指出直方图和密度曲线的三个关键方面:y 轴应该使用密度刻度,平滑隐藏了不重要的细节,直方图与条形图基本上是不同的。我们依次描述每个方面:
y 轴上的密度
长寿和遗传病直方图中的 y 轴都标记为“密度”。这个标签意味着直方图中条的总面积等于 1。简单来说,我们可以把直方图想象成一个天际线,高楼密集的地方人口更多,我们可以从矩形的面积中找到任意区间的观察比例。例如,在遗传病直方图中,从 3.5 到 9.5 的矩形大约包含了 10%的品种:6(宽度)× 0.017(高度)大约是 0.10。如果所有的区间宽度相同,那么无论 y 轴表示计数还是密度,天际线看起来都会一样。但是在这个直方图中将 y 轴改为计数会给出一个误导性的图片,右尾的一个非常大的矩形。
平滑处理
使用直方图时,我们隐藏了地毯图中个别纱线的细节,以便查看分布的一般特征。平滑是指这个过程,将点集替换为矩形;我们选择不展示数据集中的每个点,以揭示更广泛的趋势。我们可能希望平滑这些点,因为这是一个样本,我们相信观察到的值附近的其他值是合理的,和/或者我们想要关注的是一般结构而不是个别观测值。没有地毯,我们无法确定一个箱子中的点在哪里。平滑的密度曲线,就像我们之前展示的关于寿命的那个,也具有总面积和为 1 的属性。密度曲线使用平滑的核函数来展开个别纱线,有时被称为核密度估计(KDE)。
条形图 ≠ 直方图
对于定性数据,条形图与直方图起着类似的作用。条形图以视觉方式展示不同组的“流行度”或频率。然而,我们不能像直方图那样解释条形图的形状。在此设置中,尾部和对称性没有意义。此外,类别的频率由条的高度表示,宽度不包含信息。接下来的两个条形图显示了关于各个类别品种数量的相同信息;它们唯一的区别在于条的宽度。在极端情况下,最右边的图表完全消除了条,并通过单个点来表示每个计数。在没有连接线的情况下,这种图称为点图。阅读这个线图时,我们可以看到只有少数品种不适合儿童:
`kid_counts` `=` `dogs``.``groupby``(``[``'``kids``'``]``)``.``count``(``)`
`kid_counts` `=` `kid_counts``.``reindex``(``[``"``high``"``,` `"``medium``"``,` `"``low``"``]``)`
现在我们已经讨论了如何检查单个特征的分布情况,接下来我们转向当我们想要查看两个特征及其关系时的情况。
关系中要寻找的内容
当我们研究多个变量时,我们不仅要检查它们的分布,还要检查它们之间的关系。在这一部分中,我们考虑特征对并描述需要寻找的内容。根据特征类型,表 10-3 提供了绘制图表类型的指南。对于两个特征来说,类型的组合(全是定量、全是定性或混合)很重要。我们逐个考虑每种组合。
两个定量特征
如果两个特征都是定量的,那么我们通常用散点图来考察它们的关系。散点图中的每个点表示一个观测值的一对数值的位置。因此,我们可以将散点图视为一个二维的地毯图。
在散点图中,我们寻找线性和简单的非线性关系,并检查这些关系的强度。我们还要看看是否将一个或两个特征进行变换会导致线性关系。
下面的散点图展示了狗品种的体重和身高(两者都是定量的):
`px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``weight``'``,`
`marginal_x``=``"``rug``"``,` `marginal_y``=``"``rug``"``,`
`labels``=``{``'``height``'``:``'``Height (cm)``'``,` `'``weight``'``:``'``Weight (kg)``'``}``,`
`width``=``350``,` `height``=``250``)`
我们观察到,身高高于平均水平的狗往往体重也高于平均水平。这种关系呈非线性:对于较高的狗,体重的变化速度比对于较矮的狗要快。事实上,如果我们将狗看作基本上呈盒状,那么对于类似比例的盒子,盒子内的内容的重量与其长度呈立方关系是有道理的。
需要注意的是,两个单变量图缺少双变量图中的信息——关于两个特征如何一起变化的信息。实际上,两个定量特征的直方图不包含足够的信息来创建特征的散点图。我们必须谨慎行事,不要过多解读一对单变量图。相反,我们需要使用表 10-3 中适当行中列出的图之一(散点图,平滑曲线,等高线图,热图,分位数-分位数图),以了解两个定量特征之间的关系。
当一个特征是数值的而另一个是定性的时候,表 10-3 提出了不同的建议。我们接下来描述它们。
一个定性和一个定量变量
要检验定量和定性特征之间的关系,我们通常使用定性特征将数据分成组,并比较这些组中定量特征的分布。例如,我们可以比较小型、中型和大型狗品种的身高分布,同时绘制三个重叠的密度曲线:
我们看到小型和中型品种的身高分布都呈双峰性,每组左峰较大。此外,小型和中型组的身高范围比大型品种组要大。
并排的箱线图提供了跨组分布的类似比较。箱线图提供了一种简单的方法,可以粗略了解分布情况。同样,小提琴图沿轴为每个组绘制密度曲线。曲线被翻转以创建对称的“小提琴”形状。小提琴图旨在弥合密度曲线和箱线图之间的差距。我们为品种的高度创建了箱线图(左)和小提琴图(右),给出了大小标签:
狗的三个身高箱线图,每种大小的狗一个,清楚地表明大小分类是基于身高的,因为组间的身高范围几乎没有重叠。(由于平滑处理,这在密度曲线中并不明显。)在这些箱线图中我们看不到小型和中型狗群体的双峰性,但我们仍然可以看到大型狗与其他两组相比具有较窄的分布。
箱线图(也称为箱线图)是对分布的几个重要统计数据的视觉总结。箱子代表第 25 百分位数、中位数和第 75 百分位数,箱须显示尾部,还绘制了异常大或小的值。箱线图不能像直方图或密度曲线那样显示形状。它们主要显示对称性和偏斜、长/短尾巴以及异常大/小的值(也称为异常值)。
图 10-2 是对箱线图各部分的视觉解释。从中位数不在箱子中间可以看出不对称性,箱须的长度显示了尾部的大小,超出箱须的点显示了异常值。最大值被视为异常值,因为它出现在右侧箱须之外。
图 10-2. 带有标记的箱线图摘要统计的图示
当我们研究两个定性特征之间的关系时,我们的重点在于比例,接下来我们会详细解释。
两个定性特征
对于两个定性特征,我们经常比较一个特征在另一个特征定义的子组中的分布。实际上,我们保持一个特征不变,并绘制另一个特征的分布。为此,我们可以使用用于显示一个定性特征分布的一些相同图表,例如线图或条形图。举个例子,让我们来研究品种对儿童适宜性和品种大小之间的关系。
要研究这两个定性特征之间的关系,我们计算三组比例(分别对应低、中、高适宜性)。在每个适宜性类别中,我们找出小型、中型和大型狗的比例。这些比例显示在下表中。注意,每列总和为 1(相当于 100%):
`prop_table_t`
| 儿童 | 高 | 中 | 低 |
|---|---|---|---|
| 大小 | |||
| --- | --- | --- | --- |
| 大型 | 0.37 | 0.29 | 0.1 |
| 中等 | 0.36 | 0.34 | 0.2 |
| 小型 | 0.27 | 0.37 | 0.7 |
下面的线图提供了这些比例的可视化。每个适宜性级别都有一条“线”(连接的点集)。连接的点显示了适宜性类别内大小的分布。我们看到,适宜性低的品种主要是小型犬:
`fig` `=` `px``.``line``(``prop_table_t``,` `y``=``prop_table_t``.``columns``,`
`x``=``prop_table_t``.``index``,` `line_dash``=``'``kids``'``,`
`markers``=``True``,` `width``=``500``,` `height``=``250``)`
`fig``.``update_layout``(`
`yaxis_title``=``"``proportion``"``,` `xaxis_title``=``"``Size``"``,`
`legend_title``=``"``Suitability <br>for children``"`
`)`
我们还可以将这些比例呈现为一组并列条形图,如下所示:
`fig` `=` `px``.``bar``(``prop_table_t``,` `y``=``prop_table_t``.``columns``,` `x``=``prop_table_t``.``index``,`
`barmode``=``'``group``'``,` `width``=``500``,` `height``=``250``)`
`fig``.``update_layout``(`
`yaxis_title``=``"``proportion``"``,` `xaxis_title``=``"``Size``"``,`
`legend_title``=``"``Suitability <br>for children``"`
`)`
到目前为止,我们已经涵盖了包含一个或两个特征的可视化。在下一节中,我们将讨论涵盖超过两个特征的可视化。
多元设置中的比较
当我们检查一个分布或者关系时,我们经常希望能够跨数据子群组进行比较。这个在额外因素的条件下进行的过程通常导致涉及三个或更多变量的可视化。在这一节中,我们将解释如何阅读用于可视化多个变量的常用图表。
例如,让我们比较重复类别之间身高和寿命之间的关系。首先,我们将重复(狗学习新命令的典型次数)从六个类别合并为四个:<15、15-25、25-40 和 40+:
`rep_replacements` `=` `{`
`'``80-100``'``:` `'``40+``'``,` `'``40-80``'``:` `'``40+``'``,`
`'``<5``'``:` `'``<15``'``,` `'``5-15``'``:` `'``<15``'``,`
`}`
`dogs` `=` `dogs``.``assign``(`
`repetition``=``dogs``[``'``repetition``'``]``.``replace``(``rep_replacements``)``)`
现在每个组大约有 30 种品种,并且更少的类别使关系更容易解读。这些类别在散点图中由不同形状的符号表示:
`px``.``scatter``(``dogs``.``dropna``(``subset``=``[``'``repetition``'``]``)``,` `x``=``'``height``'``,` `y``=``'``longevity``'``,`
`symbol``=``'``repetition``'``,` `width``=``450``,` `height``=``250``,`
`labels``=``{``'``height``'``:``'``Height (cm)``'``,`
`'``longevity``'``:``'``Typical lifespan (yr)``'``,`
`'``repetition``'``:``'``Repetition``'``}``,`
`)`
如果“重复”特征内部有更多级别,这种图表将很难解释。
分面图提供了显示这三个特征的另一种方法:
`px``.``scatter``(``dogs``.``dropna``(``subset``=``[``'``repetition``'``]``)``,`
`x``=``'``height``'``,` `y``=``'``longevity``'``,` `trendline``=``'``ols``'``,`
`facet_col``=``'``repetition``'``,` `facet_col_wrap``=``2``,`
`labels``=``{``'``height``'``:``'``Height (cm)``'``,`
`'``longevity``'``:``'``Typical lifespan (yr)``'``}``)`
这四个散点图中的每一个展示了不同重复范围内寿命与身高之间的关系。通过分离散点图,我们可以更好地评估两个定量特征之间的关系如何随着子组的变化而变化。我们还可以更轻松地看到每个重复范围内身高和寿命的范围。我们可以看到,较大的品种往往寿命较短。另一个有趣的特征是,这些线的斜率相似,但 40+重复的线大约比其他线低 1.5 年。这些品种的平均寿命比其他重复类别低约 1.5 年,无论其身高如何。
在这里,我们总结了当我们有三个(或更多)特征时进行比较的各种绘图技术:
两个定量和一个定性
我们已经通过散点图演示了这种情况,根据定性特征的类别变化标记,或者通过每个类别的散点图面板。
两个定性和一个定量特征
我们已经看到了根据品种大小的箱线图集合中,我们可以通过并排箱线图比较不同子组的分布的基本形状。当我们有两个或更多定性特征时,我们可以根据其中一个定性特征将箱线图组织成组。
三个定量特征
当我们绘制两个定量特征和一个定性特征时,我们可以使用类似的技术。这次,我们将其中一个定量特征转换为序数特征,其中每个类别通常具有大致相同数量的记录。然后,我们制作其他两个特征的分面散点图。我们再次寻找跨分面的关系相似性。
三个定性特征
当我们检查定性特征之间的关系时,我们会检查在另一个定性特征定义的子组内一个特征的比例。 在上一节中,一个图中的三条线图和并排的条形图都显示了这些比较。 对于三(或更多)个定性特征,我们可以继续根据特征级别的组合细分数据,并使用线图、点图、并排条形图等比较这些比例。 但是这些图往往在进一步细分时变得越来越难以理解。
注意
将可视化分解以查看由定性特征确定的数据子组是否会改变关系是一种良好的做法。 这种技术称为对特征进行“控制”。 当您在散点图中看到一个线性关系有上升趋势,但在散点图的一些或全部方面中却逆转为下降趋势时,您可能会感到惊讶。 这种现象称为“辛普森悖论”。 这种悖论也可能发生在定性特征中。 在伯克利,男性的研究生院录取率高于女性的情况曾是一个著名案例,但在每个项目中单独检查时,录取率更青睐于女性。 问题在于,女性更多地申请了录取率较低的项目。
涉及多个分类变量的比较可能会很快变得复杂,因为类别组合的可能性增多。例如,有 3 × 4 = 12 种大小重复的组合(如果我们保留原始的重复类别,会有 18 种组合)。检查 12 个子组的分布可能会很困难。此外,我们面临的问题是子组中观察数太少。尽管狗数据框中有近 200 行,但一半的大小重复组合观察次数为 10 或更少。 (当一个特征具有缺失值时,会丢失一个观察。)当我们比较与定量数据的关系时,也会出现这种“维度灾难”。仅仅三个定量变量在多面图中的一些散点图很容易观察到有太少的观察以确认两个变量在子组之间的关系形状。
现在我们已经看到了在探索性数据分析中常用的可视化实例,我们继续讨论一些 EDA 的高级指导原则。
探索指南
到目前为止,在本章中,我们介绍了特征类型的概念,看到了特征类型如何帮助确定要制作的图表,并描述了如何在可视化中读取分布和关系。 EDA 依赖于建立这些技能和灵活地发展对数据的理解。
在第九章中,我们开发了数据质量检查和特征转换,以提高它们在数据分析中的实用性,从而展示了 EDA 的实际应用。接下来是一些指导你进行绘图探索数据时的问题:
-
特征 X 的值是如何分布的?
-
特征 X 和特征 Y 之间的关系如何?
-
特征 X 的分布在由特征 Z 定义的子组中是否相同?
-
特征 X 中是否存在任何不寻常的观察?(X,Y)的组合中是否存在?在 Z 的子组中的 X 中是否存在?
当回答每个问题时,重要的是将你的答案与测量的特征和上下文联系起来。采用积极好奇的探索方法也很重要。为了指导你的探索,问自己“接下来怎么办”和“这对于什么有意义”的问题,例如以下问题:
-
你有理由期望一个组/观察结果可能不同吗?
-
为什么你对形状的发现很重要?
-
哪些额外的比较可能为调查增加价值?
-
是否有任何重要的特征可以进行比较/对比?
在这个过程中,重要的是偶尔离开电脑,深思熟虑你的工作。你可能想要阅读关于这个主题的额外文献,或者去找领域内的专家讨论你的发现。例如,对于一个不寻常的观察可能有很好的理由,领域内的人可以帮助澄清并提供更多背景。
我们将在接下来的具体 EDA 示例中将这些准则付诸实践。
示例:房屋销售价格
在最后一节中,我们使用前一节的问题进行探索性分析来引导我们的调查。尽管 EDA 通常从数据整理阶段开始,但出于演示目的,我们处理的数据已经部分清理,以便我们可以专注于探索感兴趣的特征。还要注意,我们没有详细讨论优化可视化的细节;该主题在第十一章中有涵盖。
我们的数据来自旧金山纪事报(SFChron)网站。该数据包括 2003 年 4 月至 2008 年 12 月旧金山地区出售的所有房屋的完整列表。由于我们没有计划将我们的发现推广到时间段和位置之外,我们正在使用的是普查,人口与访问框架相匹配,样本包括整个人口。
至于粒度,每条记录代表了在指定时间段内在旧金山湾区销售的房屋。这意味着如果一套房在这段时间内销售了两次,那么表中将有两条记录。如果旧金山湾区的一套房在此期间没有上市出售,则它不会出现在数据集中。
数据位于数据帧sfh_df中:
`sfh_df`
| 城市 | 邮政编码 | 街道 | 价格 | 卧室数 | 居住面积 | 建筑面积 | 时间戳 | |
|---|---|---|---|---|---|---|---|---|
| 0 | Alameda | 94501.0 | 1001 Post Street | 689000.0 | 4.0 | 4484.0 | 1982.0 | 2004-08-29 |
| 1 | Alameda | 94501.0 | 1001 Santa Clara Avenue | 880000.0 | 7.0 | 5914.0 | 3866.0 | 2005-11-06 |
| 2 | Alameda | 94501.0 | 1001 Shoreline Drive #102 | 393000.0 | 2.0 | 39353.0 | 1360.0 | 2003-09-21 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 521488 | 温莎 | 95492.0 | 9998 Blasi Drive | 392500.0 | NaN | 3111.0 | NaN | 2008-02-17 |
| 521489 | 温莎 | 95492.0 | 9999 Blasi Drive | 414000.0 | NaN | 2915.0 | NaN | 2008-02-17 |
| 521490 | 温莎 | 95492.0 | 999 Gemini Drive | 325000.0 | 3.0 | 7841.0 | 1092.0 | 2003-09-21 |
521491 rows × 8 columns
The dataset does not have an accompanying codebook, but we can determine the features and their storage types by inspection:
`sfh_df``.``info``(``)`
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 521491 entries, 0 to 521490
Data columns (total 8 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 city 521491 non-null object
1 zip 521462 non-null float64
2 street 521479 non-null object
3 price 521491 non-null float64
4 br 421343 non-null float64
5 lsqft 435207 non-null float64
6 bsqft 444465 non-null float64
7 timestamp 521491 non-null datetime64[ns]
dtypes: datetime64ns, float64(5), object(2)
memory usage: 31.8+ MB
Based on the names of the fields, we expect the primary key to consist of some combination of city, zip code, street address, and date.
Sale price is our focus, so let’s begin by exploring its distribution. To develop your intuition about distributions, make a guess about the shape of the distribution before you start reading the next section. Don’t worry about the range of prices, just sketch the general shape.
Understanding Price
It seems that a good guess for the shape of the distribution of sale price might be highly skewed to the right with a few very expensive houses. The following summary statistics confirm this skewness:
`percs` `=` `[``0``,` `25``,` `50``,` `75``,` `100``]`
`prices` `=` `np``.``percentile``(``sfh_df``[``'``price``'``]``,` `percs``,` `method``=``'``lower``'``)`
`pd``.``DataFrame``(``{``'``price``'``:` `prices``}``,` `index``=``percs``)`
| price | |
|---|---|
| 0 | 22000.00 |
| 25 | 410000.00 |
| 50 | 555000.00 |
| 75 | 744000.00 |
| 100 | 20000000.00 |
The median is closer to the lower quartile than the upper quartile. Also, the maximum is 40 times the median! We might wonder whether that $20M sale price is simply an anomalous value or whether there are many houses that sold at such a high price. To find out, we can zoom in on the right tail of the distribution and compute a few high percentiles:
`percs` `=` `[``95``,` `97``,` `98``,` `99``,` `99.5``,` `99.9``]`
`prices` `=` `np``.``percentile``(``sfh_df``[``'``price``'``]``,` `percs``,` `method``=``'``lower``'``)`
`pd``.``DataFrame``(``{``'``price``'``:` `prices``}``,` `index``=``percs``)`
| price | |
|---|---|
| 95.00 | 1295000.00 |
| 97.00 | 1508000.00 |
| 98.00 | 1707000.00 |
| 99.00 | 2110000.00 |
| 99.50 | 2600000.00 |
| 99.90 | 3950000.00 |
We see that 99.9% of the houses sold for under 20M sale is indeed a rarity. Let’s examine the histogram of sale prices below $4M:
`under_4m` `=` `sfh_df``[``sfh_df``[``'``price``'``]` `<` `4_000_000``]``.``copy``(``)`
`px``.``histogram``(``under_4m``,` `x``=``'``price``'``,` `nbins``=``50``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``'``price``'``:``'``Sale price (USD)``'``}``)`
Even without the top 0.1%, the distribution remains highly skewed to the right, with a single mode around $500,000. Let’s plot the histogram of the log-transformed sale price. The logarithm transformation often does a good job at converting a right-skewed distribution into one that is more symmetric:
`under_4m``[``'``log_price``'``]` `=` `np``.``log10``(``under_4m``[``'``price``'``]``)`
`px``.``histogram``(``under_4m``,` `x``=``'``log_price``'``,` `nbins``=``50``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``'``log_price``'``:``'``Sale price (log10 USD)``'``}``)`
We see that the distribution of log-transformed sale price is roughly symmetric. Now that we have an understanding of the distribution of sale price, let’s consider the so-what questions posed in the previous section on EDA guidelines.
What Next?
我们已经描述了销售价格的形状,但我们需要考虑形状的重要性,并寻找可能有所不同的分布的比较组。
形状很重要,因为基于对称分布的模型和统计性质往往比高度偏斜的分布更具有稳健和稳定的特性(我们在第十五章中详细讨论这个问题)。因此,我们主要使用对数转换后的销售价格进行分析。而且,我们可能还会选择限制分析范围在 400 万美元以下的销售价格,因为超级昂贵的房屋可能表现出截然不同的行为。
至于可能进行的比较,我们需要看背景情况。房地产市场在此期间迅速上涨,然后市场崩溃。因此,比如说,2004 年的销售价格分布可能与 2008 年市场崩盘前的情况有很大不同。为了进一步探索这一观点,我们可以分析价格随时间的变化。或者,我们可以固定时间,分析价格与其他感兴趣特征之间的关系。这两种方法都可能是有价值的。
我们将焦点缩小到一年(在第十一章中我们将研究时间维度)。我们将数据限制在 2004 年的销售情况,因此上涨的房价对我们研究的分布和关系的影响应该是有限的。为了限制非常昂贵和大型的房屋的影响,我们还将数据集限制在售价低于 400 万美元和面积小于 12,000 平方英尺的房屋范围内。这个子集仍然包含大型和昂贵的房屋,但不会过于夸张。稍后,我们将进一步限制我们的研究范围到几个感兴趣的城市:
`def` `subset``(``df``)``:`
`return` `df``.``loc``[``(``df``[``'``price``'``]` `<` `4_000_000``)` `&`
`(``df``[``'``bsqft``'``]` `<` `12_000``)` `&`
`(``df``[``'``timestamp``'``]``.``dt``.``year` `==` `2004``)``]`
`sfh` `=` `sfh_df``.``pipe``(``subset``)`
`sfh`
| 城市 | 邮政编码 | 街道 | 价格 | 卧室数量 | 生活面积 | 建筑面积 | 时间戳 | |
|---|---|---|---|---|---|---|---|---|
| 0 | 阿拉米达 | 94501.00 | 1001 Post Street | 689000.00 | 4.00 | 4484.00 | 1982.00 | 2004-08-29 |
| 3 | 阿拉米达 | 94501.00 | 1001 Shoreline Drive #108 | 485000.00 | 2.00 | 39353.00 | 1360.00 | 2004-09-05 |
| 10 | 阿拉米达 | 94501.00 | 1001 Shoreline Drive #306 | 390000.00 | 2.00 | 39353.00 | 1360.00 | 2004-01-25 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 521467 | 温莎 | 95492.00 | 9960 Herb Road | 439000.00 | 3.00 | 9583.00 | 1626.00 | 2004-04-04 |
| 521471 | 温莎 | 95492.00 | 9964 Troon Court | 1200000.00 | 3.00 | 20038.00 | 4281.00 | 2004-10-31 |
| 521478 | 温莎 | 95492.00 | 9980 Brooks Road | 650000.00 | 3.00 | 45738.00 | 1200.00 | 2004-10-24 |
105996 rows × 8 columns
对于这些数据,销售价格分布的形状保持不变——价格仍然高度右偏。我们继续使用这个子集来解决是否有任何重要特征需要与价格一起研究的问题。
研究其他特征
除了销售价格外,这是我们的主要关注点,还有一些可能对我们的调查很重要的其他特征,例如房屋大小、地块(或物业)大小和卧室数量。我们探索这些特征的分布以及它们与销售价格和彼此之间的关系。
由于房屋和地产的大小可能与其价格有关,因此猜测这些特征也可能向右倾斜,因此我们对建筑物大小进行对数转换是合理的:
`sfh` `=` `sfh``.``assign``(``log_bsqft``=``np``.``log10``(``sfh``[``'``bsqft``'``]``)``)`
我们比较了常规和对数比例尺上建筑面积的分布:
`fig` `=` `make_subplots``(``1``,``2``)`
`fig``.``add_trace``(``go``.``Histogram``(``x``=``sfh``[``'``bsqft``'``]``,` `histnorm``=``'``percent``'``,`
`nbinsx``=``60``)``,` `row``=``1``,` `col``=``1``)`
`fig``.``add_trace``(``go``.``Histogram``(``x``=``sfh``[``'``log_bsqft``'``]``,` `histnorm``=``'``percent``'``,`
`nbinsx``=``60``)``,` `row``=``1``,` `col``=``2``)`
`fig``.``update_xaxes``(``title``=``'``Building size (ft²)``'``,` `row``=``1``,` `col``=``1``)`
`fig``.``update_xaxes``(``title``=``'``Building size (ft², log10)``'``,` `row``=``1``,` `col``=``2``)`
`fig``.``update_yaxes``(``title``=``"``percent``"``,` `row``=``1``,` `col``=``1``)`
`fig``.``update_yaxes``(``range``=``[``0``,` `18``]``)`
`fig``.``update_layout``(``width``=``450``,` `height``=``250``,` `showlegend``=``False``)`
`fig`
分布是单峰的,峰值约为 1,500 平方英尺,许多房屋的面积超过 2,500 平方英尺。我们已经确认了我们的直觉:对数转换后的建筑面积几乎是对称的,尽管保持了轻微的倾斜。对于地块面积的分布也是如此。
鉴于房屋和地块面积都呈现偏斜分布,两者的散点图很可能也应采用对数比例尺:
`sfh` `=` `sfh``.``assign``(``log_lsqft``=``np``.``log10``(``sfh``[``'``lsqft``'``]``)``)`
我们比较了有和没有对数转换的图形:
左边的散点图是以原始单位表示的,这使得很难辨别它们之间的关系,因为大多数点都挤在绘图区域的底部。相比之下,右边的散点图显示了一些有趣的特征:沿着散点图底部有一条水平线,看起来许多房屋的地块大小相同,无论建筑物的大小如何;而且似乎地块和建筑物大小之间存在轻微的正对数线性关系。
让我们看一下地块大小的一些较低分位数,试图弄清这个异常值:
`percs` `=` `[``0.5``,` `1``,` `1.5``,` `2``,` `2.5``,` `3``]`
`lots` `=` `np``.``percentile``(``sfh``[``'``lsqft``'``]``.``dropna``(``)``,` `percs``,` `method``=``'``lower``'``)`
`pd``.``DataFrame``(``{``'``lot_size``'``:` `lots``}``,` `index``=``percs``)`
| lot_size | |
|---|---|
| 0.50 | 436.00 |
| 1.00 | 436.00 |
| 1.50 | 436.00 |
| 2.00 | 436.00 |
| 2.50 | 436.00 |
| 3.00 | 782.00 |
我们发现了一些有趣的东西:大约有 2.5%的房屋的地块面积为 436 平方英尺。这太小了,几乎没有意义,因此我们记录下这个异常值以供进一步调查。
另一种衡量房屋大小的方法是卧室数量。由于这是一个离散的定量变量,我们可以将其视为定性特征并制作条形图。
旧金山湾区的房屋往往比较小,所以我们猜测分布会在三处达到峰值,并向右倾斜,有些房屋有五六间卧室。让我们来检查一下:
`br_cat` `=` `sfh``[``'``br``'``]``.``value_counts``(``)``.``reset_index``(``)`
`px``.``bar``(``br_cat``,` `x``=``"``br``"``,` `y``=``"``count``"``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``'``br``'``:``'``Number of bedrooms``'``}``)`
条形图证实了我们一般的想法。然而,我们发现有些房屋有超过 30 间卧室!这有点难以置信,也指向另一个可能的数据质量问题。由于记录包括房屋的地址,我们可以在房地产应用程序上再次检查这些值。
与此同时,让我们将卧室数量转换为有序特征,将所有大于 8 的值重新分配为 8+,并使用转换后的数据重新创建条形图:
`eight_up` `=` `sfh``.``loc``[``sfh``[``'``br``'``]` `>``=` `8``,` `'``br``'``]``.``unique``(``)`
`sfh``[``'``new_br``'``]` `=` `sfh``[``'``br``'``]``.``replace``(``eight_up``,` `8``)`
`br_cat` `=` `sfh``[``'``new_br``'``]``.``value_counts``(``)``.``reset_index``(``)`
`px``.``bar``(``br_cat``,` `x``=``"``new_br``"``,` `y``=``"``count``"``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``'``new_br``'``:``'``Number of bedrooms``'``}``)`
我们可以看到,即使我们将所有具有 8 个或更多卧室的房屋归为一类,它们的数量也不多。分布几乎对称,高峰出现在 3 个卧室,大约有相同比例的房屋有两个或四个卧室,同样有一个或五个卧室的房屋。存在不对称性,有少数房屋拥有六个或更多卧室。
现在我们来研究卧室数量与销售价格之间的关系。在我们继续之前,我们先保存目前已经完成的转换:
`def` `log_vals``(``df``)``:`
`return` `df``.``assign``(``log_price``=``np``.``log10``(``df``[``'``price``'``]``)``,`
`log_bsqft``=``np``.``log10``(``df``[``'``bsqft``'``]``)``,`
`log_lsqft``=``np``.``log10``(``df``[``'``lsqft``'``]``)``)`
`def` `clip_br``(``df``)``:`
`eight_up` `=` `df``.``loc``[``df``[``'``br``'``]` `>``=` `8``,` `'``br``'``]``.``unique``(``)`
`new_br` `=` `df``[``'``br``'``]``.``replace``(``eight_up``,` `8``)`
`return` `df``.``assign``(``new_br``=``new_br``)`
`sfh` `=` `(``sfh_df`
`.``pipe``(``subset``)`
`.``pipe``(``log_vals``)`
`.``pipe``(``clip_br``)`
`)`
现在我们准备考虑卧室数量与其他变量之间的关系。
深入探讨关系
让我们从检查不同卧室数量的房屋价格分布开始。我们可以通过箱线图来完成这个任务:
`px``.``box``(``sfh``,` `x``=``'``new_br``'``,` `y``=``'``price``'``,` `log_y``=``True``,` `width``=``450``,` `height``=``250``,`
`labels``=``{``'``new_br``'``:``'``Number of bedrooms``'``,``'``price``'``:``'``Sale price (USD)``'``}``)`
中位数销售价格随着卧室数量从一增加到五而增加,但对于最大的房屋(超过六个卧室的房屋),对数转换后的销售价格分布几乎相同。
我们预期一居室的房屋比四居室的房屋小。我们还可能猜想,六个或更多卧室的房屋在大小和价格上是类似的。为了深入了解,我们考虑一种将价格除以建筑面积的转换,得到每平方英尺价格的方法。我们想要检查这个特征是否对所有房屋都是恒定的;换句话说,价格是否主要由房屋大小决定。为此,我们查看了大小和价格、每平方英尺价格和大小之间的关系:
`sfh` `=` `sfh``.``assign``(`
`ppsf``=``sfh``[``'``price``'``]` `/` `sfh``[``'``bsqft``'``]``,`
`log_ppsf``=``lambda` `df``:` `np``.``log10``(``df``[``'``ppsf``'``]``)``)`
我们创建了两个散点图。左侧图显示价格与建筑面积(都进行了对数转换),右侧图显示每平方英尺价格(进行了对数转换)与建筑面积之间的关系。此外,每个图中都添加了一个平滑曲线,反映了大致相同大小建筑的局部平均价格或每平方英尺价格:
左侧图表显示了我们的预期—更大的房屋成本更高。我们还看到这些特征之间大致存在对数关联。
此图中的右侧图表非常有趣地呈现了非线性特征。我们看到较小的房屋每平方英尺的成本比较大的房屋更高,而较大房屋的每平方英尺价格相对平稳。这一特征似乎非常有趣,因此我们将每平方英尺价格转换保存为sfh:
def compute_ppsf(df):
return df.assign(
ppsf=df['price'] / df['bsqft'],
log_ppsf=lambda df: np.log10(df['ppsf']))
到目前为止,我们还没有考虑价格与位置之间的关系。这个数据集中来自 150 多个不同城市的房屋销售数据。有些城市只有少数销售记录,而其他城市则有数千条。我们继续缩小数据范围,并在接下来的几个城市中研究关系。
修正位置
你可能听过这样的表达:房地产有三个重要因素—位置、位置、位置。比较不同城市的房价可能会为我们的调查带来额外的见解。
我们检查了旧金山东湾一些城市的数据:里士满、埃尔塞里托、奥尔巴尼、伯克利、核桃溪、拉莫林达(这是拉斐特、莫拉加和奥林达三个相邻的卧室社区的组合),以及皮德蒙特。
让我们开始比较这些城市的销售价格分布:
`cities` `=` `[``'``Richmond``'``,` `'``El Cerrito``'``,` `'``Albany``'``,` `'``Berkeley``'``,`
`'``Walnut Creek``'``,` `'``Lamorinda``'``,` `'``Piedmont``'``]`
`px``.``box``(``sfh``.``query``(``'``city in @cities``'``)``,` `x``=``'``city``'``,` `y``=``'``price``'``,`
`log_y``=``True``,` `width``=``450``,` `height``=``250``,`
`labels``=``{``'``city``'``:``'``'``,` `'``price``'``:``'``Sale price (USD)``'``}``)`
箱线图显示,拉莫林达和皮德蒙特的房屋更昂贵,里士满最便宜,但许多城市的销售价格存在重叠。
接下来,我们将使用分面散点图更仔细地检查每个四个城市的房价与房屋大小之间的关系:
`four_cities` `=` `[``"``Berkeley``"``,` `"``Lamorinda``"``,` `"``Piedmont``"``,` `"``Richmond``"``]`
`fig` `=` `px``.``scatter``(``sfh``.``query``(``"``city in @four_cities``"``)``,`
`x``=``"``bsqft``"``,` `y``=``"``log_ppsf``"``,` `facet_col``=``"``city``"``,` `facet_col_wrap``=``2``,`
`labels``=``{``'``bsqft``'``:``'``Building size (ft²)``'``,`
`'``log_ppsf``'``:` `"``Price per square foot``"``}``,`
`trendline``=``"``ols``"``,` `trendline_color_override``=``"``black``"``,`
`)`
`fig``.``update_layout``(``xaxis_range``=``[``0``,` `5500``]``,` `yaxis_range``=``[``1.5``,` `3.5``]``,`
`width``=``450``,` `height``=``400``)`
`fig``.``show``(``)`
价格每平方英尺与建筑面积的关系大致呈对数线性,对于这四个地点每个都存在负相关。虽然不是平行的,但似乎在房子方面存在地理位置的提升,例如伯克利的房屋每平方英尺比里士满的房屋贵约 250 美元。我们还看到皮德蒙特和拉莫林达是更昂贵的城市,在这两个城市中,与较小房屋相比,较大房屋每平方英尺的价格没有同样的降低。这些图表支持“地段,地段,地段”的格言。
在探索性数据分析(EDA)中,我们经常回顾早期的图表,以检查新发现是否为先前的可视化增添了新的见解。持续盘点我们的发现并利用它们指导我们进一步的探索至关重要。让我们总结一下到目前为止我们的发现。
EDA 发现
我们的 EDA 揭示了几个有趣的现象。简而言之,其中一些最显著的是:
-
销售价格和建筑面积呈右偏分布,有一个主模式。
-
每平方英尺价格随建筑面积的增加呈非线性下降趋势,较小的房屋每平方英尺的成本高于较大的房屋,并且每平方英尺的价格在大房屋中大致保持不变。
-
更理想的地理位置为房屋的销售价格增加了一点,对于不同大小的房屋增加的金额大致相同。
我们可以(也应该)进行许多其他探索,还有几个我们应该进行的检查。这些包括调查占地面积 436 的价值和用在线房地产应用程序交叉检查异常房屋,例如 30 卧室房屋和 2000 万美元的房屋。
我们将我们的调查限制在一年内,后来又缩小到几个城市。这种缩小帮助我们控制可能干扰发现简单关系的特征。例如,由于数据是在几年内收集的,销售日期可能混淆了销售价格和卧室数量之间的关系。在其他时候,我们希望考虑时间对价格的影响。为了检查随时间的价格变化,我们经常制作折线图,并对通货膨胀进行调整。我们在第十一章重新审视这些数据时,考虑数据范围并更仔细地观察时间趋势。
尽管简短,本节传达了 EDA 实践的基本方法。有关不同数据集的扩展案例研究,请参阅第十二章。
总结
在本章中,我们介绍了名义、序数和数值特征类型及其在数据分析中的重要性。当面对数据集时,我们展示了如何查阅数据字典和数据本身,以确定每列的特征类型。我们还解释了存储类型与特征类型不应混淆。由于大部分 EDA 是通过统计图表进行的,我们描述了如何识别和解释出现的形状和模式,以及如何将其与正在绘制的数据联系起来。最后,我们提供了进行 EDA 的指导方针,并提供了一个示例。
有一种方法可能对你在开发关于特征分布和关系直觉很有帮助,那就是在绘制图表之前先对你将看到的内容进行猜测。试着草拟或描述你认为分布形状会是什么样子,然后再绘制图表。例如,具有自然下限/上限值的变量往往在边界的对面有一个长尾。收入分布(下限为 0)往往有一个长尾在右侧,而考试成绩(上限为 100)往往有一个长尾在左侧。你可以对关系的形状做类似的猜测。我们发现价格和房屋大小几乎呈对数-对数线性关系。当你对形状有了直觉后,进行探索性数据分析(EDA)就变得更容易;你可以更容易地识别出图表显示出令人惊讶的形状。
本章的重点是“阅读”可视化结果。在第十一章中,我们提供了如何创建信息丰富、有效和美观的图表的样式指南。这一章中许多思想也在这里得到了遵循,但我们并未特别指出它们。
第十一章:数据可视化
作为数据科学家,我们创建数据可视化是为了理解我们的数据,并向其他人解释我们的分析。图表应该传达一个信息,我们的工作是尽可能清晰地传达这个信息。
在第十章中,我们将统计图的选择与所绘制数据的类型联系起来;我们还介绍了许多标准图,并展示了如何解读它们。在本章中,我们讨论了有效数据可视化的原则,这些原则使得观众更容易理解我们图中的信息。我们讨论了如何选择轴的刻度,如何通过平滑和聚合处理大量数据,如何进行有意义的比较,如何融入研究设计,并添加上下文信息。我们还展示了如何使用plotly这一流行的 Python 绘图包来创建图表。
撰写关于数据可视化的一章的一个棘手之处在于,可视化软件包经常变动,因此我们展示的任何代码可能很快就会过时。由于这个原因,一些书籍完全避免使用代码。我们相反地取得了平衡,覆盖了广泛有用的高级数据可视化原则。然后我们单独包含实际的绘图代码来实现这些原则。当新软件可用时,读者仍然可以使用我们的原则来指导他们的可视化创建。
选择比例以揭示结构
在第十章中,我们探讨了 2003 年至 2009 年间旧金山湾区房屋销售价格。让我们重新看一下这个例子,并看一看销售价格的直方图:
`px``.``histogram``(``sfh``,` `x``=``'``price``'``,` `nbins``=``100``,`
`labels``=``{``'``price``'``:``"``Sale price (USD)``"``}``,` `width``=``350``,` `height``=``250``)`
虽然这个图准确地显示了数据,但大部分可见的箱子都挤在图的左侧。这使得理解价格分布变得困难。
通过数据可视化,我们希望展示数据的重要特征,如分布的形状和两个或更多特征之间的关系。正如这个例子所示,当我们生成初始图后,仍然有其他方面需要考虑。在本节中,我们涵盖了帮助我们决定如何调整轴限制、放置刻度标记和应用变换的比例原则。我们首先检查何时以及如何调整图形以减少空白区域;换句话说,我们试图用数据填充我们图的数据区域。
填充数据区域
正如我们从销售价格直方图中看到的那样,当大部分数据出现在绘图区域的一个小部分时,读取分布就变得困难了。当这种情况发生时,数据的重要特征,如多模式和偏斜,可能会被掩盖。散点图也存在类似问题。当所有点都挤在散点图的一个角落时,很难看到分布的形状,因此也很难从形状中获得任何见解。
当存在少数异常大的观测时会出现这个问题。为了更好地观察数据的主要部分,我们可以通过调整 x 或 y 轴的限制删除这些观测值,或者在绘制前从数据中删除异常值。无论哪种情况,我们都会在标题或图表本身中提到这种排除。
让我们利用这个想法来改善销售价格的直方图。在接下来的并列图中,我们通过改变 x 轴的限制来裁剪数据。在左图中,我们排除了价格超过 $2 百万的房屋。这样做使得大部分房屋分布的形状在图中更加清晰。例如,我们可以更容易地观察到偏斜和较小的次要模式。在右图中,我们分别展示了分布长尾右侧的详细信息(在线查看更大的版本):
左图的 x 轴包含 0,但右图的 x 轴从 $2M 开始。我们考虑下一步在轴上是否包含或排除 0。
包含零
我们通常不需要在轴上包含 0,特别是如果包含它会使填充数据区域变得困难。例如,让我们制作一个反映狗品种平均寿命与平均身高关系的散点图。(此数据集首次在第十章介绍;它包含 172 种品种的多个特征。)
左图的 x 轴从 10 厘米开始,因为所有狗至少有这么高,类似地,y 轴从 6 年开始。右侧的散点图在两个轴上都包含 0。这将数据推到数据区域的顶部,并留下不利于我们看到线性关系的空白空间。
在柱状图中,通常需要包含 0,这样柱子的高度直接与数据值相关联。例如,我们创建了两个比较狗品种寿命的柱状图。左图包含 0,而右图不包含:
从右图可以轻易地错误地得出小品种的寿命是大品种的两倍的结论。
在处理比例时,我们通常也希望包含 0,因为比例范围从 0 到 1。下图显示了每种类型中品种的比例:
在柱状图和散点图中,包含 0 使您能够更准确地比较各类别的相对大小。
早些时候,当我们调整轴时,实际上是从绘图区域中删除了数据。虽然当少数观测值异常大(或小)时这是一个有用的策略,但在偏斜分布中效果较差。在这种情况下,我们通常需要转换数据以更好地了解其形状。
通过变换显示形状
另一种常见的调整比例的方法是转换数据或图的轴。我们使用变换来处理偏斜的数据,以便更容易检查分布。当变换产生对称分布时,对称性带有有用的建模属性(参见第十五章)。
有多种方法可以转换数据,但对数变换往往特别有用。例如,在下面的图表中,我们重新生成了旧金山房屋销售价格的两个直方图。顶部直方图是原始数据。对于下面的直方图,我们在绘制前对价格取了对数(以 10 为底):
对数变换使得价格分布更对称。现在我们可以更轻松地看到分布的重要特征,如大约 10 5.85 处的模式,约为 70 万美元,以及次要模式接近 10 5.55 处,即 35 万美元。
使用对数变换的缺点是实际值不那么直观——在这个例子中,我们需要将值转换回美元才能理解销售价格。因此,我们通常更喜欢将轴变换为对数刻度,而不是数据本身。这样,我们可以在轴上看到原始的值:
具有对数刻度 x 轴的直方图基本上显示了与转换数据的直方图相同的形状。但由于在美元刻度上的箱子宽度相等,但在对数美元刻度上绘制,所以右侧的箱子变窄了。还要注意,y 轴上的 μ 是 10 − 6 。
对数变换还可以显示散点图中的形状。在这里,我们将建筑物大小绘制在 x 轴上,地块大小绘制在 y 轴上。在这个图中很难看出形状,因为许多点都挤在数据区域的底部:
然而,当我们同时使用对数刻度的 x 轴和 y 轴时,关系的形状更容易看出:
`px``.``scatter``(``sfh``,` `x``=``'``bsqft``'``,` `y``=``'``lsqft``'``,`
`log_x``=``True``,` `log_y``=``True``,`
`labels``=``{``"``bsqft``"``:` `"``Building size (sq ft)``"``,`
`"``lsqft``"``:` `"``Lot size (sq ft)``"``}``,`
`width``=``350``,` `height``=``250``)`
通过变换后的坐标轴,我们可以看到在对数尺度上,地块大小随建筑物大小大致呈线性增加。对数变换将大值——比其他值大几个数量级的值——拉向中心。这种变换有助于填充数据区域并揭示隐藏的结构,就像我们在房价分布和房屋尺寸与地块大小之间的关系中所看到的那样。
除了设置轴的限制和变换轴外,我们还要考虑绘图的纵横比—长度与宽度的比例。调整纵横比称为银行业务,在下一节中,我们将展示如何通过银行业务来帮助揭示特征之间的关系。
解读关系
对于散点图,我们尝试选择刻度,使得两个特征之间的关系大致沿着 45 度线。这种缩放称为银行业务向 45 度。这样做是因为我们的眼睛更容易观察到与直线偏离的情况和趋势。例如,我们重现了显示狗品种寿命与身高关系的图:
`px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``longevity``'``,` `width``=``300``,` `height``=``250``,`
`labels``=``{``"``height``"``:` `"``Height (cm)``"``,`
`"``longevity``"``:` `"``Typical lifespan (yr)``"``}``)`
散点图已经银行业务向 45 度,我们更容易看到数据大致沿直线分布以及它们在极端情况下的偏离。
当银行业务向 45 度倾斜有助于我们判断数据是否遵循线性关系时,当存在明显的曲率时,很难弄清楚关系的形式。在这种情况下,我们尝试能够使数据沿直线分布的变换(例如,请参阅图 11-1)。对数变换在揭示曲线关系的一般形式时非常有用。
透过拉直揭示关系
我们经常使用散点图来观察两个特征之间的关系。例如,在这里我们绘制了狗品种的身高与体重的关系:
`px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``weight``'``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``"``height``"``:` `"``Height (cm)``"``,` `"``weight``"``:` `"``Weight (lb)``"``}``)`
我们看到更高的狗体重更重,但这种关系并非线性的。
当看起来两个变量之间存在非线性关系时,尝试应用对数尺度到 x 轴、y 轴或两者都是有用的。让我们在变换轴的散点图中寻找线性关系。这里我们重新绘制了狗品种的体重与身高的图,但这次我们将对 y 轴应用了对数尺度:
`px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``weight``'``,` `log_y``=``True``,`
`labels``=``{``"``height``"``:` `"``Height (cm)``"``,` `"``weight``"``:` `"``Weight (lb)``"``}``,`
`width``=``300``,` `height``=``300``)`
这个图显示了大致线性的关系,在这种情况下,我们称狗的体重和身高之间存在对数-线性关系。
通常,当我们看到变换一个或两个轴后呈现线性关系时,我们可以使用表 11-1 来揭示原始变量的关系(在表中,a 和 b 是常数)。我们进行这些变换是因为这样更容易看出点是否沿着一条线分布,而不是看它们是否遵循幂律或指数律。
表 11-1。应用变换时两个变量之间的关系
| x 轴 | y 轴 | 关系 | 也称为 |
|---|---|---|---|
| 无转换 | 无转换 | 线性:y = a x + b | 线性 |
| 对数尺度 | 无转换 | 对数:y = a log x + b | 线性-对数 |
| 无变换 | 对数比例 | 指数: y = b a x | 对数-线性 |
| 对数比例 | 对数比例 | 功率: y = b x a | 对数-对数 |
如表 11-1 所示,对数变换可以揭示几种常见类型的关系。正因如此,对数变换被认为是变换的利器。作为另一个虽然是人为的例子,图 11-1 中最左边的图显示了x和y之间的曲线关系。中间的图显示了log(y)和x之间不同的曲线关系;该图看起来也是非线性的。进一步的对数变换,最右边显示了*log(y)对log(x)*的图。这个图证实了数据具有对数-对数(或者功率)关系,因为变换后的点沿着一条直线分布。
图 11-1. 散点图展示了对数变换如何“拉直”两个变量之间的曲线关系
调整比例是数据可视化中的重要实践。虽然对数变换很灵活,但并不适用于所有呈现偏斜或曲率的情况。例如,有时值都大致相同数量级,对数变换影响有限。另一个需要考虑的变换是平方根变换,通常适用于计数数据。
在下一节中,我们将探讨平滑的原则,这是在需要可视化大量数据时使用的方法。
平滑和聚合数据
当我们有大量数据时,通常不想绘制所有单个数据点。以下散点图展示了来自樱花的数据,这是每年四月份在华盛顿特区举办的一场 10 英里长的比赛,当时樱花盛开。这些数据来自比赛网站,包括 1999 年至 2012 年所有注册男性跑步者的官方时间和其他信息。我们将跑步者的年龄绘制在 x 轴上,比赛时间绘制在 y 轴上:
这个散点图包含超过 70,000 个数据点。由于数据点过多,很多点会重叠在一起。这是一个常见的问题,称为过度绘制。在这种情况下,过度绘制使我们无法看到时间和年龄之间的关系。在这个图中,我们唯一能看到的是一群非常年轻的跑步者,这可能指向数据质量存在问题。为了解决过度绘制的问题,我们使用平滑技术在绘图前聚合数据。
平滑技术揭示形状
直方图是一种熟悉的平滑绘图类型。直方图通过将点放入 bin 中并绘制每个 bin 的条形来汇总数据值。这里的平滑意味着我们不能区分 bin 中各个点的位置:点被平滑地分配到它们的 bin 中。对于直方图,bin 的面积对应于 bin 中点的百分比(或计数或比例)。(通常 bin 宽度相等,我们采取一种简便的方式将 bin 的高度标记为比例。)
下面的直方图显示了狗品种寿命的分布:
在这个直方图上方是一个折线图,为每个数据值绘制一条线。我们可以看到在最高的 bin 中,即使少量数据也会导致折线图的重叠。通过平滑处理折线图的点,直方图显示了分布的一般形状。在这种情况下,我们看到许多品种的寿命约为 12 年。有关如何阅读和解释直方图的更多信息,请参见第十章。
另一种常见的平滑技术是核密度估计(KDE)。KDE 图使用平滑曲线而不是柱状图显示分布。在下面的图中,我们展示了同一狗寿命直方图,并覆盖了一个 KDE 曲线。KDE 曲线与直方图具有类似的形状:
把直方图看作平滑方法可能会让人惊讶。KDE 和直方图都旨在帮助我们看到值分布中的重要特征。类似的平滑技术也可以用于散点图。这是下一节的主题。
揭示关系和趋势的平滑技术
我们可以通过分 bin 数据来找到散点图的高密度区域,就像直方图一样。下面的图重新绘制了樱花大赛时间与年龄的散点图。这张图使用了六边形 bin 将点聚合在一起,并根据落入其中的点数对六边形进行了阴影处理:
`runners_over_17` `=` `runners``[``runners``[``"``age``"``]` `>` `17``]`
`plt``.``figure``(``figsize``=``(``4``,` `4``)``)`
`plt``.``hexbin``(``data``=``runners_over_17``,` `x``=``'``age``'``,` `y``=``'``time``'``,` `gridsize``=``35``,` `cmap``=``'``Blues``'``)`
`sns``.``despine``(``)`
`plt``.``grid``(``False``)`
`plt``.``xlabel``(``"``Runner age (yr)``"``)`
`plt``.``ylabel``(``"``Race time (sec)``"``)``;`
注意在 25 到 40 岁组中的高密度区域,通过图中的深色区域表示。图表告诉我们,许多这个年龄段的跑者大约在 5,000 秒(约 80 分钟)内完成比赛。(请注意,我们从这个图中删除了年轻跑者的数据。)我们还可以看到 40 到 60 岁组对应区域的向上曲率,这表明这些跑者通常比 25 到 40 岁组的跑者慢。这张图类似于热图,其中高密度区域通过更热或更亮的颜色传达。
核密度估计在二维中也适用。当我们在二维中使用 KDE 时,通常绘制结果三维形状的轮廓线图,并像读取地形图一样解读图形:
`plt``.``figure``(``figsize``=``(``5``,` `3``)``)`
`fig` `=` `sns``.``kdeplot``(``data``=``runners_over_17``,` `x``=``'``age``'``,` `y``=``'``time``'``)`
`plt``.``xlabel``(``"``Runner age (yr)``"``)`
`plt``.``ylabel``(``"``Race time (sec)``"``)``;`
这个二维 KDE 提供了与前一图中阴影区域相似的见解。我们看到 25 到 40 岁年龄组的跑步者高度集中,这些跑步者的时间大约为 5,000 秒左右。平滑使我们在数据量大时能够更清楚地看到整体情况,因为它可以显示高度集中数据值的位置和这些高集中区域的形状。否则这些区域可能无法看到。
另一种经常更具信息性的平滑方法是对具有相似 x 值的点的 y 值进行平滑处理。简单来说,我们将具有相似年龄的跑步者分组,使用五年递增:20–25、25–30、30–35 等。然后,对于每个五年的跑步者组,我们计算其平均比赛时间,绘制每组的平均时间,并连接点以形成一个“曲线”:
`times` `=` `(`
`runners_over_17``.``assign``(``age_5yr``=``runners_over_17``[``'``age``'``]` `/``/` `5` `*` `5``)`
`.``groupby``(``'``age_5yr``'``)``[``'``time``'``]``.``mean``(``)``.``reset_index``(``)`
`)`
`px``.``line``(``times``,` `x``=``'``age_5yr``'``,` `y``=``'``time``'``,`
`labels``=``{``'``time``'``:``"``Average race time (sec)``"``,` `'``age_5yr``'``:``"``Runner age (5-yr)``"``}``,`
`markers``=``True``,`
`width``=``350``,` `height``=``250``)`
这幅图再次显示,年龄在 25 到 40 岁之间的跑步者典型的跑步时间约为 5,400 秒。它还显示,老年跑步者平均完成比赛时间较长(这并不令人意外,但在早期的图表中没有那么明显)。年龄在 20 岁以下的跑步者时间的下降和 80 岁时曲线的平坦可能仅仅是这些组别中跑步者较少且更健壮的结果。另一种平滑技术使用与 KDE 类似的核平滑。我们这里不详细介绍。
分箱和核平滑技术依赖于一个调整参数,该参数指定箱的宽度或核的扩展,并且在制作直方图、KDE 或平滑曲线时通常需要指定此参数。这是下一节的主题。
平滑技术需要调整
现在我们已经看到平滑在绘图中的用处,我们转向调整的问题。对于直方图,箱宽度或者对于等宽箱,箱的数量会影响直方图的外观。这里显示的长寿左直方图有几个宽箱子,而右侧直方图有许多窄箱子(在线上查看更大图像online):
在两个直方图中,很难看出分布的形状。在几个宽箱子(左图)中,我们对分布进行了过度平滑处理,这使得无法分辨模式和尾部。另一方面,太多的箱子(右图)给出的图像与地毯图相差无几。KDE 图中有一个称为带宽的参数,其作用类似于直方图的箱宽度。
大多数直方图和 KDE 软件会自动选择直方图的箱宽度和核的带宽。然而,这些参数通常需要一点调整才能创建最有用的图形。当您创建依赖调整参数的可视化时,重要的是在最终确定一个参数值之前尝试几种不同的值。
数据减少的另一种不同方法是检查分位数。这是下一节的主题。
减少分布到分位数
在第十章中我们发现,虽然箱线图不如直方图信息丰富,但在一次比较多个组的分布时它们是有用的。箱线图基于数据的四分位数将数据减少到几个基本特征。更一般地,分位数(下四分位数、中位数和上四分位数分别是第 25、50 和 75 个分位数)在比较分布时可以提供有用的数据缩减。
当两个分布在形状上大致相似时,用直方图可能很难比较它们。例如,以下直方图展示了旧金山房屋数据中两卧室和四卧室房屋价格的分布。这些分布在形状上看起来大致相似。但是它们的分位数图可以方便地比较分布的中心、传播和尾部(在线上更大地查看在线):
`px``.``histogram``(``sfh``.``query``(``'``br in [2, 4]``'``)``,`
`x``=``'``price``'``,` `log_x``=``True``,` `facet_col``=``'``br``'``,`
`labels``=``{``'``price``'``:``"``Sale price (USD)``"``}``,`
`width``=``700``,` `height``=``250``)`
我们可以用分位数-分位数图来比较,简称为q-q 图。为了制作这种图,我们首先计算价格的两卧室和四卧室分布的百分位数(也称为分位数):
`br2` `=` `sfh``.``query``(``'``br == 2``'``)`
`br4` `=` `sfh``.``query``(``'``br == 4``'``)`
`percs` `=` `np``.``arange``(``1``,` `100``,` `1``)`
`perc2` `=` `np``.``percentile``(``br2``[``'``price``'``]``,` `percs``,` `method``=``'``lower``'``)`
`perc4` `=` `np``.``percentile``(``br4``[``'``price``'``]``,` `percs``,` `method``=``'``lower``'``)`
`perc_sfh` `=` `pd``.``DataFrame``(``{``'``percentile``'``:` `percs``,` `'``br2``'``:` `perc2``,` `'``br4``'``:` `perc4``}``)`
`perc_sfh`
| 百分位数 | br2 | br4 | |
|---|---|---|---|
| --- | --- | --- | --- |
| 0 | 1 | 1.50e+05 | 2.05e+05 |
| 1 | 2 | 1.82e+05 | 2.50e+05 |
| 2 | 3 | 2.03e+05 | 2.75e+05 |
| ... | ... | ... | ... |
| 96 | 97 | 1.04e+06 | 1.75e+06 |
| 97 | 98 | 1.20e+06 | 1.95e+06 |
| 98 | 99 | 1.44e+06 | 2.34e+06 |
99 rows × 3 columns
然后在散点图上绘制匹配的百分位数。通常我们也会显示参考线 y = x 来帮助比较:
`fig` `=` `px``.``scatter``(``perc_sfh``,` `x``=``'``br2``'``,` `y``=``'``br4``'``,` `log_x``=``True``,` `log_y``=``True``,`
`labels``=``{``'``br2``'``:` `'``Price of 2-bedroom house``'``,`
`'``br4``'``:` `'``Price of 4-bedroom house``'``}``,`
`width``=``350``,` `height``=``250``)`
`fig``.``add_trace``(``go``.``Scatter``(``x``=``[``1e5``,` `2e6``]``,` `y``=``[``1e5``,` `2e6``]``,`
`mode``=``'``lines``'``,` `line``=``dict``(``dash``=``'``dash``'``)``)``)`
`fig``.``update_layout``(``showlegend``=``False``)`
`fig`
当分位点沿直线分布时,变量具有类似形状的分布。与参考线平行的线表示中心的差异,斜率不为 1 的线表示传播的差异,而曲线表示形状的差异。从前述的分位数图中,我们看到四卧室房屋价格的分布在形状上与两卧室的分布类似,除了大约 10 万美元的偏移和稍长的右尾(对大值的上升弯曲指示)。阅读分位数图需要练习。一旦你掌握了技巧,它可以是比较分布的一种便捷方式。注意,房屋数据有超过 10 万个观察值,而分位数图将数据减少到 99 个百分位数。这种数据缩减非常有用。然而,我们并不总是想使用平滑处理器。这是下一节的主题。
当不进行平滑处理时
平滑和聚合可以帮助我们看到重要的特征和关系,但是当我们只有少数观测时,平滑技术可能会误导我们。只有少数观测时,我们更喜欢地毯图而不是直方图、箱线图和密度曲线,而是使用散点图而不是平滑曲线和密度轮廓。这似乎是显而易见的,但是当我们有大量数据时,子组中的数据量可能会迅速减少。这种现象是 维度诅咒 的一个例子。
平滑最常见的误用之一发生在箱线图中。例如,这是寿命的七个箱线图的集合,每个箱线图代表一种狗品种:
`px``.``box``(``dogs``,` `x``=``'``group``'``,` `y``=``'``longevity``'``,`
`labels``=``{``'``group``'``:``"``"``,` `'``longevity``'``:``"``Longevity (yr)``"``}``,`
`width``=``500``,` `height``=``250``)`
其中一些箱线图只有两到三个观测值。接下来的条形图是一种更可取的可视化方法:
`px``.``strip``(``dogs``,` `x``=``"``group``"``,` `y``=``"``longevity``"``,`
`labels``=``{``'``group``'``:``"``"``,` `'``longevity``'``:``"``Longevity (yr)``"``}``,`
`width``=``400``,` `height``=``250``)`
在这个图中,我们仍然可以比较各组,但我们也可以看到每组中的确切值。现在我们可以看出非运动组中只有三个品种;基于箱线图的偏斜分布的印象过于看重箱线的形状。
本节介绍了由于大型数据集而产生重叠点的过绘制问题。为了解决这个问题,我们介绍了聚合数据的平滑技术。我们看到了平滑的两个常见示例—分箱和核平滑—并将它们应用于一维和二维设置中。在一维中,这些是直方图和核密度曲线,它们都帮助我们看到分布的形状。在二维中,我们发现在保持 x 值不变的同时平滑 y 值是有用的,以可视化趋势。我们解决了需要调整平滑量以获得更多信息的直方图和密度曲线,并警告不要用太少的数据进行平滑处理。
有许多其他减少散点图中的重叠的方法。例如,我们可以使点部分透明,以便重叠点显示更暗。如果许多观测值具有相同的值(例如当寿命四舍五入到最接近的年份时),那么我们可以向值添加少量随机噪声,以减少重叠的数量。这个过程称为 jittering,它用于寿命的条形图。透明度和 jittering 对于中等大小的数据很方便。但是,它们对于大型数据集效果不佳,因为绘制所有点仍然会压倒可视化。
我们介绍的分位数图为比较具有远少于点数的分布提供了一种方法;另一种方法是使用并排箱线图,还有一种方法是在同一图中叠加 KDE 曲线。我们经常旨在比较数据子集(或组)之间的分布和关系,接下来我们将讨论几种促进各种图类型的有意义比较的设计原则。
促进有意义的比较
同样的数据可以用许多不同的方式进行可视化,决定制作哪种图表可能令人生畏。一般来说,我们的图表应该帮助读者进行有意义的比较。在本节中,我们讨论了几个有用的原则,可以提高图表的清晰度。
强调重要的差异
每当我们制作一个比较不同群体的图表时,我们都会考虑图表是否突出了重要的差异。一般来说,当图表中的对象以使比较更容易阅读的方式对齐时,读者更容易看到差异。让我们看一个例子。
美国劳工统计局发布了有关收入的数据。我们以 2020 年全职等效周收入的中位数为基础,对超过 25 岁的人群按教育水平和性别进行了分组:^(1)
`labels` `=` `{``"``educ``"``:` `"``Education``"``,`
`"``income``"``:` `"``Weekly earnings (USD)``"``,`
`"``gender``"``:` `"``Sex``"``}`
`fig` `=` `px``.``bar``(``earn``,` `x``=``"``educ``"``,` `y``=``"``income``"``,`
`facet_col``=``"``gender``"``,` `labels``=``labels``,`
`width``=``450``,` `height``=``250``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`
这些条形图显示,随着教育程度的提高,收入也增加。但可以说,一个更有趣的比较是相同教育水平下男女之间的比较。我们可以以不同的方式分组条形图,以便集中在这种比较上:
`px``.``bar``(``earn``,` `x``=``'``educ``'``,` `y``=``'``income``'``,` `color``=``'``gender``'``,`
`barmode``=``'``group``'``,` `labels``=``labels``,`
`width``=``450``,` `height``=``250``)`
这个图更好;我们可以更容易地比较每个教育水平男性和女性的收入。然而,我们可以通过垂直对齐进一步明确这种差异。我们不再使用条形,而是在每个教育水平上使用点来代表男性和女性群体:
`fig` `=` `px``.``line``(``earn``,` `x``=``'``educ``'``,` `y``=``'``income``'``,` `symbol``=``'``gender``'``,`
`color``=``'``gender``'``,` `labels``=``labels``,` `width``=``450``,` `height``=``250``)`
`fig``.``update_traces``(``marker_size``=``10``)`
这个图更清晰地显示了一个重要的差异:男性和女性的收入差距随着教育程度的提高而扩大。我们考虑了三个可视化同一数据的图表,但它们在传达信息的方式上有所不同。我们更喜欢最后一个,因为它将收入差异垂直对齐,更易于比较。
注意,在制作所有三个图表时,我们按照教育水平从少到多的顺序对教育类别进行了排序。这种排序是有意义的,因为教育水平是有序的。当我们比较名义类别时,我们使用其他排序方法。
分组顺序
对于有序特征,当我们制作图表时,我们保持类别的自然顺序,但对于名义特征,这个原则不适用。相反,我们选择一种有助于比较的顺序。在条形图中,按照高度排序条形是一个好的做法,而对于箱线图和条带图,通常按照中位数的顺序排列箱子/条带。
后面的两个条形图各自比较了不同狗品种的平均寿命:
左边的图按字母顺序排列了条形,我们更喜欢右边的图,因为它按寿命的顺序排列了条形,这样更容易比较各个类别的寿命。我们不需要来回跳跃或眯起眼睛去猜测牧群犬是否比玩具犬寿命更短。
作为另一个例子,以下两组箱线图比较了旧金山东湾地区不同城市的房屋销售价格分布:
我们更喜欢右侧的图表,因为它按每个城市的中位数价格排序了箱子。同样,这种排序使得更容易比较组(在这种情况下是城市)之间的分布。我们看到 Albany 和 Walnut Creek 的下四分位数和中位数价格大致相同,但是 Walnut Creek 的价格有更大的右偏。
在可能的情况下,按高度排序条形图中的条和箱线图中的箱子,使得我们更容易比较不同组之间的数据。用于呈现分组数据的另一种技术是堆叠。我们将在下一节介绍堆叠,并提供示例以说服您远离这种类型的图表。
避免堆积
接下来的图示是一个堆积条形图,其中每个城市都有一根条形,并按照销售房屋中卧室数从一到八个或更多的比例进行划分。这称为堆积条形图。条形图基于一个交叉制表:
`br_crosstab`
| br | 1.0 | 2.0 | 3.0 | 4.0 | 5.0 | 6.0 | 7.0 | 8.0 |
|---|---|---|---|---|---|---|---|---|
| city | ||||||||
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Albany | 1.21e-01 | 0.56 | 0.25 | 0.05 | 9.12e-03 | 1.01e-03 | 2.03e-03 | 4.05e-03 |
| Berkeley | 6.91e-02 | 0.38 | 0.31 | 0.16 | 4.44e-02 | 1.42e-02 | 6.48e-03 | 7.23e-03 |
| El Cerrito | 1.81e-02 | 0.34 | 0.47 | 0.14 | 2.20e-02 | 6.48e-03 | 0.00e+00 | 6.48e-04 |
| Piedmont | 8.63e-03 | 0.22 | 0.40 | 0.26 | 9.50e-02 | 1.29e-02 | 7.19e-03 | 1.44e-03 |
| Richmond | 3.60e-02 | 0.36 | 0.42 | 0.15 | 2.52e-02 | 7.21e-03 | 7.72e-04 | 7.72e-04 |
| Walnut Creek | 1.16e-01 | 0.35 | 0.30 | 0.18 | 4.37e-02 | 5.08e-03 | 4.12e-04 | 2.75e-04 |
绘图中的每个条形的高度都相同,因为段代表城市中一卧室或更多卧室的房屋比例,因此总和为 1 或 100%(在线查看详细信息):
`fig` `=` `px``.``bar``(``br_crosstab``,` `width``=``450``,` `height``=``300``)`
`fig``.``update_layout``(``yaxis_title``=``None``,` `xaxis_title``=``None``,`
`legend_title``=``"``# Bedrooms``"``)`
`fig``.``show``(``)`
通过简单地扫描每列第一个片段顶部,很容易比较每个城市一卧室房屋的比例。但是,比较四卧室房屋就更加困难。由于段的底部水平不对齐,我们必须用眼睛判断绘图中上下移动的段长度。这种上下移动称为基线摆动。(我们意识到,由于颜色过多,这个图在灰度中显示效果不佳,但我们的目标是引导你远离这样的图表,所以对于在线阅读的读者,我们保留了所有的颜色。)
堆叠线图更难阅读,因为我们必须判断曲线之间的间隙,而它们上下摇晃。下图显示了 1950 年至 2012 年排放量最高的 10 个国家的二氧化碳(CO[2])排放(在线查看颜色online):
由于这些线条是堆叠在一起的,很难看出特定国家的排放量如何变化,也很难比较各个国家。相反,我们可以绘制每个国家的线条,而不是堆叠,如下图所示(在线查看颜色online):
现在查看单个国家的变化以及比较国家变得更加容易,因为我们只需要判断 y 轴位置而不是具有不同基线的短垂直段。我们还在 y 轴上使用了对数刻度。现在我们可以看到,一些国家的二氧化碳排放率增长缓慢,例如美国和日本,而其他国家的增长速度更快,如中国和印度,而德国减缓了其二氧化碳排放。在每个国家的基线在图中摇晃时,几乎不可能检测到这些方面。
在这两个图中,为了更容易区分各个国家,我们使用了不同的线型和颜色。选择颜色以促进比较依赖于许多考虑因素。这是下一节的主题。
选择调色板
选择颜色在数据可视化中也起着重要作用。我们希望避免使用过亮或过暗的颜色,以免给读者的眼睛带来压力。我们还应避免对色盲人士可能困难的调色板——7%到 10%的人(主要是男性)是红绿色盲。
对于分类数据,我们希望使用一种可以清楚区分类别的调色板。图 11-2 中的顶部展示了一个示例(参见 Figure 11-2)。从上到下,这些调色板用于分类数据的定性,用于数值数据的分散调色板,其中希望吸引对大值和小值的注意力,以及用于数值数据的顺序调色板,其中希望强调大或小值。
图 11-2. 来自 ColorBrewer 2.0 的三种适合打印的调色板(在线查看颜色online)
对于数值数据,我们希望使用顺序调色板,强调光谱的一侧而不是另一侧,或者使用分散调色板,两端平衡强调,中间淡化。图 11-2 显示了顺序调色板位于底部,分散调色板位于中间(参见 Figure 11-2):
当我们想要强调低或高值(如癌症率)时,我们选择一个顺序调色板。当我们想要强调两个极端(如两党选举结果)时,我们选择一个分散调色板。
选择感知统一的颜色调色板非常重要。这意味着当数据值加倍时,可视化中的颜色在人眼中看起来会变得两倍鲜艳。我们还希望避免颜色会在我们从图表的一个部分看向另一个部分时产生余像,不同强度的颜色会使一个属性看起来比另一个属性更重要,以及色盲人士难以区分的颜色。我们强烈建议使用专为数据可视化制作的调色板或调色板生成器。
图表应该被长时间检查,因此我们应该选择不妨碍读者仔细研究图表的颜色。更重要的是,颜色的使用不应是多余的——颜色应该代表信息。相关的是,人们通常很难区分超过七种颜色,因此我们限制图表中颜色的数量。最后,颜色在纸上打印时可能与在计算机屏幕上查看时相差很大。在选择颜色时,我们考虑我们的图表将如何显示。
在可视化中进行准确的比较是如此重要,以至于研究人员已经研究了人们感知颜色和其他绘图特征(如角度和长度)差异的程度。这是下一节的主题。
绘图比较的指南
研究人员研究了人们在不同类型的图表中读取信息的准确度。他们发现以下排序,从最准确到最不准确:
-
沿着共同刻度的位置,比如在地毯图、条带图或点图中
-
在相同但非对齐刻度上的位置,比如在条形图中
-
长度,在堆叠条形图中
-
角度和斜率,在饼图中
-
面积,在堆叠线图或气泡图中
-
体积和密度,在三维条形图中
-
颜色的饱和度和色调,比如在使用半透明点进行重叠绘制时
例如,这是一个显示旧金山售出的房屋中拥有从一到八个或更多卧室比例的饼图,以及具有相同比例的条形图:
在饼图中很难判断角度,需要用实际百分比进行注释。我们还失去了卧室数量的自然顺序。条形图没有这些问题。
然而,任何规则都有例外。每个饼图中只有两个或三个片段的多个饼图可以提供有效的可视化效果。例如,一组饼图显示旧金山东湾地区六个城市中每个城市售出的两卧室房屋的比例,按比例排序,可以产生深远的可视化效果。然而,坚持使用条形图通常至少和任何饼图一样清晰。
根据这些准则,我们建议在进行比较时坚持使用位置和长度。读者可以更准确地根据位置或长度来判断比较,而不是角度、面积、体积或颜色。但是,如果我们想向图表添加额外信息,通常会使用颜色、符号和线型,除了位置和长度。我们在本章中展示了几个例子。
接下来我们转向数据设计的话题,讨论如何在可视化中反映数据收集的时间、地点和方式的方面。这是一个微妙但重要的话题。如果忽略数据的范围,我们可能会得到非常误导人的图表。
结合数据设计
当我们创建可视化时,考虑数据的范围尤为重要,特别是数据设计(见第二章)。考虑数据收集方式的问题可以影响我们的图表选择和所描绘的比较。这些考虑包括数据收集的时间和地点,以及用于选择样本的设计。我们看几个例子,展示数据范围如何影响我们所做的可视化。
数据随时间的收集
当数据随时间收集时,我们通常制作一条线图,将时间戳放在 x 轴上,将感兴趣的特征放在 y 轴上,以便观察时间趋势。例如,让我们重新审视旧金山房价数据。这些数据从 2003 年到 2008 年收集,展示了 2008/2009 年美国房地产泡沫的崩溃(美国房地产泡沫)。由于时间是这些数据范围的关键方面,让我们将销售价格可视化为时间序列。我们之前的探索显示销售价格高度偏斜,所以让我们使用百分位数而不是平均数。我们绘制中位数价格(这是我们在本章早些时候看到的一种平滑形式):
此图显示了 2003 年到 2007 年价格的上升和 2008 年的下降。但是,我们可以通过绘制几个额外的百分位数线来展示更多信息,而不仅仅是中位数。让我们为第 10、30、50(中位数)、70 和 90 百分位数的销售价格分别绘制不同的线条。当我们随时间检查价格时,通常需要校正通货膨胀,以便比较处于相同基础上。除了通货膨胀调整外,让我们将价格相对于 2003 年的起始价格进行绘制,每个百分位数的起始值为 y = 1。 (2006 年 90 百分位数的值为 1.5 表明,销售价格是 2003 年 90 百分位数的 1.5 倍)。这种归一化方法使我们能够看到住房市场不同部分的房屋危机如何影响房主:
当我们随时间跟踪第 10 百分位线的折线图时,我们发现它在 2005 年迅速上升,在接下来的几年里相对于 2003 年的值保持较高水平,然后比其他百分位数更早更快地下降。这告诉我们,较便宜的房屋,如起始家庭住房,在房地产市场崩溃中遭受了更大的波动并且损失了更多的价值。相比之下,高端住房受到的冲击较小;在 2008 年底,第 90 百分位数的房价仍高于 2003 年的价格。应用这一领域知识有助于揭示我们可能会忽略的数据趋势,并展示如何利用数据设计来改善可视化。
住房数据是一个关于特定时期内特定地理区域的完全普查的观察数据的例子。接下来,我们考虑另一个观察性研究,其中自我选择和时间段影响了可视化效果。
观察性研究
我们需要特别小心那些不构成人口普查或科学样本的数据。我们也应该注意横断面研究,无论是来自人口普查还是科学样本。在这个例子中,我们重新审视了樱花 10 英里赛的数据。在本章的早些时候,我们制作了一个平滑曲线来研究赛时与年龄之间的关系。我们在这里再次重现这个图表,以突显解释可能存在的潜在陷阱:
诱人的是,看到这张图表可能会得出结论,例如,60 岁的跑步者通常比 40 岁时需要额外花费 600 秒来完成比赛。然而,这是一个横断面研究,而不是纵向研究。该研究并未随时间跟踪人们的变化;相反,它获取了一个人群的横截面。在图表中代表 60 岁跑步者的人与代表 40 岁跑步者的人是不同的。这两组人可能在影响赛时与年龄关系的方式上有所不同。作为一个群体,比赛中的 60 岁老人可能比 40 岁的老人更适合他们的年龄。换句话说,数据设计不允许我们对个体跑步者做出结论。这个可视化并没有错,但我们在从中得出结论时需要小心。
这个设计更加复杂,因为我们有多年的比赛结果。每一年形成一个队列,一个选手组,而从一年到下一年,队列会改变。我们通过比较不同年份的选手来清晰地传达这个信息。在这里,我们分别画出了 1999 年、2005 年和 2010 年的选手线路图:
我们看到 2010 年每个年龄组的中位数比 2005 年的跑步者更高,而 2005 年的跑步者的时间比 1999 年的跑步者更高。有趣的是,随着比赛的增加,尤其是近年来新手跑步者的参与增加,比赛时间逐年放缓。这个例子显示了我们在解释模式时需要注意数据范围。在科学研究中,我们也需要牢记数据范围。这是下一节的主题。
不均匀抽样
在科学研究中,我们必须考虑样本设计,因为它可能影响我们的图表。一些样本以不均等的速率抽取个体,这需要在我们的可视化中进行考虑。我们在第八章和第九章中看到了科学研究的例子:药物滥用警告网络(DAWN)调查。这些数据来自于关于药物相关急诊访问的复杂随机研究,每个记录都带有一个权重,我们必须使用这些权重来准确地代表人口中的急诊访问情况。接下来的两个条形图显示了急诊访问类型的分布。可以在在线链接上查看更大的版本。左侧的图表未使用调查权重,右侧的图表使用了:
在未加权的条形图中,“其他”类别与“不良反应”类别一样频繁。然而,加权后,“其他”类别下降到“不良反应”的约三分之二。忽略抽样权重可能会导致分布的误导性呈现。无论是直方图、条形图、箱线图、二维轮廓还是光滑曲线,我们都需要使用权重来获得代表性的图表。数据范围的另一个影响我们图表选择的方面是数据收集的地点,这是下一节的主题。
地理数据
当我们的数据包含像纬度和经度这样的地理信息时,除了典型的图表外,我们应该考虑制作地图。例如,下图展示了美国空气质量传感器的位置,这是第十二章案例研究的重点:
注意加利福尼亚州和东海岸有更多的数据点。使用所有这些传感器数据的简单空气质量直方图会误代表美国的空气质量分布。为了将空间因素纳入分布中,我们可以使用不同颜色的标记将空气质量测量添加到地图上,并且可以按位置分面显示空气质量的直方图。
除了绘制条形、颜色和线条样式等特征外,我们还可以选择添加文本以提供背景信息,从而使我们的图表更具信息性。这是下一节的主题。
添加背景信息
在本章中,我们在图表中使用文本来提供包括测量单位的有意义的轴标签、用于类别的刻度标签和标题。在广泛分享可视化时,这是一个良好的实践。一个好的目标是在图中包含足够的上下文,使其能够独立存在——读者应该能够理解图的主要内容,而不需要到处寻找解释。尽管如此,统计图的每个元素都应该有一个目的。多余的文本或图表功能,通常称为图表垃圾,应该被消除。在本节中,我们提供了如何通过添加上下文来创建一个出版准备好的图表的简要概述。
文本上下文包括标签和标题。在刻度标记和轴上一贯使用信息丰富的标签是一个良好的做法。例如,轴标签通常受益于包含测量单位。当需要时,图表应包含标题和图例。信息丰富的标签对于其他人会看到和解释的图表尤为重要。然而,即使在进行我们自己的探索性数据分析时,我们通常也希望包含足够的上下文,以便在稍后返回分析时,我们能够轻松地理解我们绘制了什么。
标题起到了几个作用。它们描述了绘制的内容并指导读者。标题还指出了图的重要特征并评论了它们的含义。标题中重复文本中的信息是可以接受的。读者通常会快速浏览出版物并关注章节标题和可视化内容,因此图的标题应该是自包含的。
参考标记为绘图区域带来了额外的背景信息。提供基准、历史值和其他外部信息的参考点和参考线帮助进行比较和解释。例如,我们经常在分位数-分位数图上添加斜率为 1 的参考线。我们还可能在时间序列图上添加垂直线以标记特殊事件,如自然灾害。
下面的例子演示了如何向绘图添加这些上下文元素。
例子:100 米短跑时间
以下图显示了自 1968 年以来男子 100 米短跑的比赛时间。这些数据仅包括在正常风条件下、使用电子计时的户外比赛,并且仅包括在 10 秒内完成比赛的选手的时间。图表是一个基本的散点图,显示了时间与年份的关系。从这张图开始,我们对其进行了增强,以创建一幅出现在FiveThirtyEight 文章中的图表。
当我们想要准备一张图形供他人阅读时,我们考虑的是主要信息。在这种情况下,我们的主要信息是双重的:在过去的 50 年中,最佳选手们的速度越来越快,而乌塞恩·博尔特在 2009 年创下的 9.58 秒的非凡记录至今未被打破(事实上,第二快的比赛时间也属于博尔特)。我们通过添加直接陈述主要信息的标题、y 轴标签中的测量单位以及散点图中关键点的注释来提供这张图形的背景信息。此外,我们添加了一个水平参考线,在 10 秒以下的时间中仅绘制时间,以澄清世界纪录时间使用特殊符号来引起读者的注意。
这些背景信息描述了我们所绘制的内容,帮助读者看到主要信息,并指出数据中几个有趣的特征。现在,这张图可以成为幻灯片、技术报告或社交媒体帖子的有用部分。根据我们的经验,查看我们的数据分析的人记住我们的图形,而不是段落文本或方程式。为他人准备的图形添加背景信息是非常重要的。
在下一节中,我们将详细介绍如何使用 plotly Python 包创建图形。
使用 plotly 创建图形
在本节中,我们介绍了使用 plotly Python 包的基础知识,这是我们在本书中创建图形所使用的主要工具。
plotly 包相比其他绘图库有几个优点。它创建的是交互式图形而不是静态图像。当你在 plotly 中创建一个图形时,你可以平移和缩放,以查看通常太小以至于看不见的图形部分。你还可以悬停在图形元素上,比如散点图中的符号,以查看原始数据值。此外,plotly 可以使用 SVG 文件格式保存图形,这意味着即使放大后图像仍然清晰。如果你在 PDF 或书本的纸质副本中阅读本章节,我们使用了这个功能来渲染图形图像。最后,它有一个简单的“express” API,用于创建基本图形,这在进行探索性分析并想快速创建多个图形时非常有帮助。
在本节中,我们介绍了 plotly 的基础知识。如果你在这里遇到未涵盖的内容,我们建议查阅官方 plotly 文档。
图形和轨迹对象
每个 plotly 中的图形都包装在一个 Figure 对象中。Figure 对象跟踪要绘制的内容。例如,一个单独的 Figure 可以在左侧绘制散点图,在右侧绘制线图。Figure 对象还跟踪图形的布局,包括图形的大小、标题、图例和注释。
plotly.express 模块提供了一个简洁的 API 来制作图形:
`import` `plotly``.``express` `as` `px`
我们在以下代码中使用 plotly.express 制作了一张关于狗品种体重与身高的散点图。请注意,.scatter() 的返回值是一个 Figure 对象:
`fig` `=` `px``.``scatter``(`
`dogs``,` `x``=``"``height``"``,` `y``=``"``weight``"``,`
`labels``=``dict``(``height``=``"``Height (cm)``"``,` `weight``=``"``Weight (kg)``"``)``,`
`width``=``350``,` `height``=``250``,`
`)`
`fig``.``__class__`
plotly.graph_objs._figure.Figure
显示 Figure 对象将其渲染到屏幕上:
`fig`
这个特定的 Figure 包含了一个图表,但是 Figure 对象可以包含任意数量的图表。在这里,我们创建了一个包含三个散点图的分面图:
`# The plot titles are cut off; we'll fix them in the next snippet`
`px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``weight``'``,`
`facet_col``=``'``size``'``,`
`labels``=``dict``(``height``=``"``Height (cm)``"``,` `weight``=``"``Weight (kg)``"``)``,`
`width``=``550``,` `height``=``250``)`
这三个图表存储在 Trace 对象中。但是,我们尽量避免手动操作 Trace 对象。相反,plotly 提供了自动创建分面子图的函数,例如我们这里使用的 px.scatter 函数。现在我们已经看到如何制作简单的图表,接下来我们展示如何修改图表。
修改布局
我们经常需要更改图表的布局。例如,我们可能想调整图表的边距或轴范围。我们可以使用 Figure.update_layout() 方法来实现这一点。在我们制作的分面散点图中,标题被裁剪,因为图表的边距不够大。我们可以通过 Figure.update_layout() 来更正这个问题:
`fig` `=` `px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``weight``'``,`
`facet_col``=``'``size``'``,`
`labels``=``dict``(``height``=``"``Height (cm)``"``,` `weight``=``"``Weight (kg)``"``)``,`
`width``=``550``,` `height``=``250``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``40``)``)`
`fig`
.update_layout() 方法允许我们修改布局的任何属性。这包括图表标题 (title)、边距 (margins 字典) 和是否显示图例 (showlegend)。plotly 文档中列出了完整的 布局属性。
Figure 对象还具有 .update_xaxes() 和 .update_yaxes() 函数,与 .update_layout() 类似。这两个函数允许我们修改轴的属性,如轴限制 (range)、刻度数 (nticks) 和轴标签 (title)。在这里,我们调整了 y 轴的范围,并更改了 x 轴的标题。我们还为图表添加了标题,并更新了布局,以避免标题被裁剪:
`fig` `=` `px``.``scatter``(`
`dogs``,` `x``=``"``weight``"``,` `y``=``"``longevity``"``,`
`title``=``"``Smaller dogs live longer``"``,`
`width``=``350``,` `height``=``250``,`
`)`
`fig``.``update_yaxes``(``range``=``[``5``,` `18``]``,` `title``=``"``Typical lifespan (yr)``"``)`
`fig``.``update_xaxes``(``title``=``"``Average weight (kg)``"``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`
`fig`
plotly 包提供了许多绘图方法;我们在下一节中描述了其中几种。
绘图函数
plotly 方法包括折线图、散点图、条形图、箱线图和直方图。每种类型的绘图方法的 API 都类似。数据帧是第一个参数。然后我们可以使用 x 和 y 关键字参数指定数据帧的列放置在 x 轴和 y 轴上。
我们首先制作了每年樱花赛跑者的中位时间折线图:
`px``.``line``(``medians``,` `x``=``'``year``'``,` `y``=``'``time``'``,` `width``=``350``,` `height``=``250``)`
接下来,我们制作不同大小狗品种的平均寿命条形图:
`lifespans` `=` `dogs``.``groupby``(``'``size``'``)``[``'``longevity``'``]``.``mean``(``)``.``reset_index``(``)`
`px``.``bar``(``lifespans``,` `x``=``'``size``'``,` `y``=``'``longevity``'``,`
`width``=``350``,` `height``=``250``)`
plotly 中的绘图方法还包含用于制作分面图的参数。我们可以根据颜色、绘图符号或线条样式在同一图表中进行分面,或者将其分为多个子图。以下是每种方式的示例。我们首先制作了一个关于狗品种身高和体重的散点图,并使用不同的绘图符号和颜色来进行分面:
`fig` `=` `px``.``scatter``(``dogs``,` `x``=``'``height``'``,` `y``=``'``weight``'``,`
`color``=``'``size``'``,` `symbol``=``'``size``'``,`
`labels``=``dict``(``height``=``"``Height (cm)``"``,`
`weight``=``"``Weight (kg)``"``,` `size``=``"``Size``"``)``,`
`width``=``350``,` `height``=``250``)`
`fig`
下一个图显示了每个品种大小寿命的并列直方图。这里我们按列进行分面:
`fig` `=` `px``.``histogram``(``dogs``,` `x``=``'``longevity``'``,` `facet_col``=``'``size``'``,`
`width``=``550``,` `height``=``250``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`
要获取绘图函数的完整列表,请参阅plotly的主要文档或plotly.express,这是我们在本书中主要使用的plotly子模块。
要为图表添加背景信息,我们使用plotly的注释方法;接下来将对此进行描述。
注释
Figure.add_annotation()方法在plotly图表上放置注释。这些注释是带有文本和可选箭头的线段。箭头的位置使用x和y参数设置,我们可以使用ax和ay参数将文本从其默认位置移动。在这里,我们使用信息注释散点图的一个点:
fig = px.scatter(dogs, x='weight', y='longevity',
labels=dict(weight="Weight (kg)",
longevity="Typical lifespan (yr)"),
width=350, height=250)
fig.add_annotation(text='Chihuahuas live 16.5 years on average!',
x=2, y=16.5,
ax=30, ay=5,
xshift=3,
xanchor='left')
fig
本节介绍了使用plotly Python 包创建图表的基础知识。我们介绍了Figure对象,这是plotly用来存储图表及其布局的对象。我们讨论了plotly提供的基本图表类型,以及通过调整布局和坐标轴以及添加注释来自定义图表的几种方式。在下一节中,我们将简要比较plotly与 Python 中其他常见的可视化工具。
其他可视化工具
有许多软件包和工具可用于创建数据可视化。在本书中,我们主要使用plotly。但值得了解的是其他几种常用工具。在本节中,我们将plotly与matplotlib和图形语法工具进行比较。
matplotlib
matplotlib库是为 Python 创建的最早的数据可视化工具之一。因此,它被广泛使用,并拥有大量的软件包生态系统。值得注意的是,用于pandas DataFrames的内置绘图方法使用matplotlib。在matplotlib之上构建的一个流行包被称为seaborn。与仅使用matplotlib相比,seaborn提供了一个更简单的 API 来创建统计图,例如带有置信区间的点图。事实上,seaborn的 API 启发了plotly的 API。如果你将plotly和seaborn的代码并排比较,你会发现创建基本图表的方法使用了相似的代码。
使用matplotlib的一个优势是其流行性。因为许多现有项目使用它,因此在线获取帮助来创建或调整图表相对容易。对于本书而言,使用plotly的主要优势是我们创建的图表是交互式的。matplotlib中的图表通常是静态图像,不支持平移、缩放或悬停在标记上。尽管如此,我们预计matplotlib仍将继续用于数据分析。事实上,本书中的几个图表是使用seaborn和matplotlib制作的,因为plotly尚未支持我们想要制作的所有图表。
图形语法
图形语法是由李·威尔金森(Lee Wilkinson)为创建数据可视化而开发的理论。其基本思想是使用常见的构建块来制作图表。例如,柱状图和点图几乎是相同的,唯一的区别在于柱状图绘制矩形,而点图绘制点。这个想法被体现在图形语法中,它会说柱状图和点图只在它们的“几何”组件上有所不同。图形语法是一个优雅的系统,我们可以使用它来描述几乎我们想要制作的每一种图表。
这个系统在流行的绘图库ggplot2(用于 R 编程语言)和Vega(用于 JavaScript)中实现。Vega-Altair,一个 Python 包,提供了一种使用 Python 创建Vega图表的方法,我们鼓励感兴趣的读者阅读其文档。
使用像Vega-Altair这样的图形语法工具可以实现可视化的灵活性。像plotly一样,altair也创建交互式可视化。然而,这些工具的 Python API 可能没有plotly的 API 那么直接。在本书中,我们通常不需要plotly无法创建的图表,因此我们选择了plotly更简单的 API。
为了简洁起见,我们略去了 Python 中许多其他绘图工具。但对于本书的目的,依赖于plotly提供了交互性和灵活性的有用平衡。
总结
当我们分析数据集时,我们使用可视化来发现数据中难以以其他方式检测到的模式。数据可视化是一个迭代的过程。我们创建一个图,然后决定是否进行调整或选择一个全新的图表类型。本章介绍了我们用来做出这些决策的原则。
我们首先介绍了规模的原则,并看到通过改变或转换绘图轴来调整比例尺可以揭示数据中隐藏的结构。然后,我们讨论了平滑和聚合技术,这些技术帮助我们处理大型数据集,否则会导致过度绘制。为了进行有意义的比较,我们应用了知觉原则,例如调整基线以使线条、柱状图和点易于比较。我们展示了如何考虑数据设计来改进可视化效果。并且我们看到,为绘图添加上下文有助于读者理解我们的信息。
在这一章之后,你应该能够创建一个图表,并理解哪些调整能使图表更有效。当你学会如何制作信息丰富的可视化时,请保持耐心并进行迭代。我们没有人能在第一次尝试时做出完美的图表,在分析中我们发现了新的东西时,我们会继续改进我们的图表。然后,在向他人展示我们的发现时,我们会筛选出少数几个最能说服我们未来读者分析正确性和重要性的图表。这甚至可能会导致创建一个更好地传达我们发现的新图表,我们会进行迭代开发。
在下一章中,我们将通过一个扩展案例研究来展示我们迄今为止在书中学到的一切。我们希望你会对自己已经能做到的事情感到惊讶。
^(1) 美国政府调查仍然根据性别的二元定义收集数据,但进展正在取得。例如,从 2022 年开始,美国公民可以在他们的护照申请中选择“X”作为他们的性别标记。
第十二章:案例研究:空气质量测量的准确性有多高?
加利福尼亚州容易发生森林大火,以至于该州的居民(如本书的作者们)有时会说加利福尼亚州“总是在火灾中”。2020 年,40 起不同的火灾使得整个州笼罩在烟雾之中,迫使成千上万的人撤离,并造成超过 120 亿美元的损失(图 12-1)。
图 12-1 卫星图像,显示 2020 年 8 月加利福尼亚州被烟雾覆盖的情况(图片来源于Wikipedia,根据 CC BY-SA 3.0 IGO 许可)。
在像加利福尼亚这样的地方,人们利用空气质量测量来了解他们需要采取哪些保护措施。根据情况,人们可能希望戴口罩、使用空气过滤器或完全避免外出。
在美国,一个重要的空气质量信息来源是由美国政府运行的空气质量系统(AQS)。AQS 在美国各地的位置上安装了高质量的传感器,并向公众提供它们的数据。这些传感器经过严格的校准到严格的标准——实际上,AQS 传感器通常被视为准确度的黄金标准。然而,它们也有一些缺点。这些传感器很昂贵:每台通常在 15,000 至 40,000 美元之间。这意味着传感器数量较少,并且它们之间的距离较远。住在传感器远离的人可能无法获取 AQS 数据用于个人使用。此外,AQS 传感器不提供实时数据。由于数据经过广泛的校准,它们仅每小时发布一次,并且有一到两小时的时间滞后。实质上,AQS 传感器精确但不及时。
相比之下,PurpleAir传感器,我们在第三章介绍过,售价约 250 美元,可以轻松安装在家中。由于价格较低,成千上万的美国人购买了这些传感器用于个人使用。这些传感器可以连接到家庭 WiFi 网络,因此可以轻松监测空气质量,并可以向 PurpleAir 报告数据。2020 年,数千名 PurpleAir 传感器的所有者将他们传感器的测量数据公开发布。与 AQS 传感器相比,PurpleAir 传感器更及时。它们每两分钟报告一次测量结果,而不是每小时。由于部署了更多的 PurpleAir 传感器,更多的人住在接近传感器的地方,可以利用这些数据。然而,PurpleAir 传感器的准确性较低。为了使传感器价格更加合理,PurpleAir 使用了更简单的方法来计算空气中的颗粒物。这意味着 PurpleAir 的测量可能会报告空气质量比实际情况更差(见Josh Hug 的博客文章)。实质上,PurpleAir 传感器倾向于及时但准确性较低。
在本章中,我们计划利用 AQS 传感器的测量结果来改进 PurpleAir 的测量。这是一项重大任务,我们首先采用了由卡洛琳·巴克约恩、布雷特·甘特和安德里亚·克莱门斯从美国环境保护署开发的分析方法。巴克约恩及其同事的工作非常成功,以至于截至撰写本文时,类似 AirNow 的官方美国政府地图,如火灾与烟雾地图,都包括 AQS 和 PurpleAir 传感器,并对 PurpleAir 数据应用了巴克约恩的校正。
我们的工作遵循数据科学生命周期,从考虑问题和可用数据的设计和范围开始。我们大部分的工作都花费在清洗和整理数据以进行分析,但我们也进行了探索性数据分析并建立了一个泛化模型。我们首先考虑问题、设计和数据的范围。
问题、设计和范围
理想情况下,空气质量的测量应该既准确又及时。不准确或有偏差的测量可能意味着人们对空气质量的重视程度不够。延迟的警报可能会使人们暴露于有害空气中。在引言中提到廉价空气质量传感器的流行背景让我们对它们的质量和实用性产生了兴趣。
两种不同类型的仪器测量了一个自然现象——空气中颗粒物的数量。AQS 传感器具有测量误差小和偏差可以忽略不计的优势(参见第二章)。另一方面,PurpleAir 仪器的精度较低;测量值具有更大的变异性并且也有偏差。我们最初的问题是:我们能否利用 AQS 测量结果使 PurpleAir 的测量更准确?
我们现在处于一个有大量可用数据的情况中。我们可以访问少量来自 AQS 的高质量测量数据,并且可以从成千上万个 PurpleAir 传感器中获取数据。为了缩小我们问题的焦点,我们考虑如何利用这两个数据源来改进 PurpleAir 的测量。
这两个来源的数据包括传感器的位置。因此,我们可以尝试将它们配对,找到接近每个 AQS 传感器的 PurpleAir 传感器。如果它们很接近,那么这些传感器实际上就是在测量相同的空气。我们可以将 AQS 传感器视为地面真实情况(因为它们非常精确),并研究给定真实空气质量情况下 PurpleAir 测量值的变化。
即使合并的 AQS 和 PurpleAir 传感器对相对较少,将我们发现的任何关系普遍化到其他 PurpleAir 传感器似乎是合理的。如果 AQS 和 PurpleAir 测量之间存在简单的关系,那么我们可以利用这种关系来调整来自任何 PurpleAir 传感器的测量结果,使其更加准确。
我们已经明确了我们的问题:我们能否建立 PurpleAir 传感器读数与相邻 AQS 传感器读数之间的关系模型?如果可以,那么希望可以利用该模型改进 PurpleAir 的读数。剧透警告:确实可以!
本案例研究很好地整合了本书这一部分介绍的概念。它为我们提供了一个机会,看看数据科学家如何在实际环境中处理、探索和可视化数据。特别是,我们看到了如何通过大型、不太精确的数据集来增强小型、精确的数据集的有用性。像这样结合大型和小型数据集对数据科学家来说尤为激动人心,并广泛适用于从社会科学到医学等其他领域。
在接下来的部分中,我们通过寻找彼此靠近的 AQS 和 PurpleAir 传感器配对来进行数据处理。我们特别关注 PM2.5 颗粒的读数,这些颗粒直径小于 2.5 微米。这些颗粒足以吸入到肺部,对健康构成最大风险,并且在木材燃烧中特别常见。
寻找配对传感器
我们的分析始于寻找 AQS 和 PurpleAir 传感器的配对,即安装在基本相邻位置的传感器。这一步骤很重要,因为它使我们能够减少可能导致传感器读数差异的其他变量的影响。想象一下,如果我们比较一个放置在公园中的 AQS 传感器和一个放置在繁忙高速公路旁的 PurpleAir 传感器会发生什么。这两个传感器将有不同的读数,部分原因是传感器暴露在不同的环境中。确保传感器真正配对让我们可以声明传感器读数差异是由传感器构造方式和小范围空气波动引起的,而不是其他潜在混淆变量的影响。
由 EPA 小组进行的 Barkjohn 分析找到了在彼此距离不到 50 米的 AQS 和 PurpleAir 传感器配对。该小组联系了每个 AQS 站点,确认是否也安装了 PurpleAir 传感器。这额外的工作让他们对传感器配对的真实配对性充满信心。
在这一部分中,我们探索和清理来自 AQS 和 PurpleAir 的位置数据。然后,我们进行一种连接操作,构建一个潜在的配对传感器列表。我们不会自己联系 AQS 站点;相反,我们将在后续章节中使用 Barkjohn 确认的配对传感器列表。
我们下载了 AQS 和 PurpleAir 传感器列表,并将数据保存在文件 data/list_of_aqs_sites.csv 和 data/list_of_purpleair_sensors.json 中。让我们开始将这些文件读入 pandas DataFrames。首先,我们检查文件大小,看看它们是否可以合理地加载到内存中:
`!`ls -lLh data/list_of*
-rw-r--r-- 1 sam staff 4.8M Oct 27 16:54 data/list_of_aqs_sites.csv
-rw-r--r-- 1 sam staff 3.8M Oct 22 16:10 data/list_of_purpleair_sensors.json
两个文件都相对较小。让我们从 AQS 站点列表开始。
处理 AQS 站点列表
我们已过滤出仅显示测量 PM2.5 的 AQS 站点的AQS 站点地图,然后使用该地图的 Web 应用程序下载站点列表作为 CSV 文件。现在我们可以将其加载到pandas DataFrame中:
`aqs_sites_full` `=` `pd``.``read_csv``(``'``data/list_of_aqs_sites.csv``'``)`
`aqs_sites_full``.``shape`
(1333, 28)
表中有 28 列。让我们检查列名:
`aqs_sites_full``.``columns`
Index(['AQS_Site_ID', 'POC', 'State', 'City', 'CBSA', 'Local_Site_Name',
'Address', 'Datum', 'Latitude', 'Longitude', 'LatLon_Accuracy_meters',
'Elevation_meters_MSL', 'Monitor_Start_Date', 'Last_Sample_Date',
'Active', 'Measurement_Scale', 'Measurement_Scale_Definition',
'Sample_Duration', 'Sample_Collection_Frequency',
'Sample_Collection_Method', 'Sample_Analysis_Method',
'Method_Reference_ID', 'FRMFEM', 'Monitor_Type', 'Reporting_Agency',
'Parameter_Name', 'Annual_URLs', 'Daily_URLs'],
dtype='object')
为了找出对我们最有用的哪些列,我们参考了 AQS 在其网站上提供的数据字典。在那里,我们确认数据表包含有关 AQS 站点的信息。因此,我们可能期望粒度对应于 AQS 站点,即每行表示一个单个站点,标记为AQS_Site_ID的列是主键。我们可以通过每个 ID 的记录计数来确认这一点:
`aqs_sites_full``[``'``AQS_Site_ID``'``]``.``value_counts``(``)`
06-071-0306 4
19-163-0015 4
39-061-0014 4
..
46-103-0020 1
19-177-0006 1
51-680-0015 1
Name: AQS_Site_ID, Length: 921, dtype: int64
看起来有些站点在这个数据框中出现了多次。不幸的是,这意味着粒度比单个站点级别更细。为了弄清楚站点重复的原因,让我们更仔细地查看一个重复站点的行:
`dup_site` `=` `aqs_sites_full``.``query``(``"``AQS_Site_ID ==` `'``19-163-0015``'``"``)`
根据列名选择几列进行检查——那些听起来可能会揭示重复原因的列:
`some_cols` `=` `[``'``POC``'``,` `'``Monitor_Start_Date``'``,`
`'``Last_Sample_Date``'``,` `'``Sample_Collection_Method``'``]`
`dup_site``[``some_cols``]`
| POC | Monitor_Start_Date | Last_Sample_Date | Sample_Collection_Method | |
|---|---|---|---|---|
| 458 | 1 | 1/27/1999 | 8/31/2021 | R & P Model 2025 PM-2.5 Sequential Air Sampler... |
| 459 | 2 | 2/9/2013 | 8/26/2021 | R & P Model 2025 PM-2.5 Sequential Air Sampler... |
| 460 | 3 | 1/1/2019 | 9/30/2021 | Teledyne T640 at 5.0 LPM |
| 461 | 4 | 1/1/2019 | 9/30/2021 | Teledyne T640 at 5.0 LPM |
POC列看起来对于区分表中的行是有用的。数据字典对该列有以下说明:
这是用于区分在同一站点上测量相同参数的不同仪器的“参数发生代码”。
因此,站点19-163-0015有四个仪器都测量 PM2.5。数据框的粒度是单个仪器级别。
由于我们的目标是匹配 AQS 和 PurpleAir 传感器,我们可以通过选择每个 AQS 站点的一个仪器来调整粒度。为此,我们根据站点 ID 对行进行分组,然后在每个组中取第一行:
`def` `rollup_dup_sites``(``df``)``:`
`return` `(`
`df``.``groupby``(``'``AQS_Site_ID``'``)`
`.``first``(``)`
`.``reset_index``(``)`
`)`
`aqs_sites` `=` `(``aqs_sites_full`
`.``pipe``(``rollup_dup_sites``)``)`
`aqs_sites``.``shape`
(921, 28)
现在行数与唯一 ID 的数量匹配。
要与 PurpleAir 传感器匹配 AQS 站点,我们只需要站点 ID、纬度和经度。因此,我们进一步调整结构,仅保留这些列:
`def` `cols_aqs``(``df``)``:`
`subset` `=` `df``[``[``'``AQS_Site_ID``'``,` `'``Latitude``'``,` `'``Longitude``'``]``]`
`subset``.``columns` `=` `[``'``site_id``'``,` `'``lat``'``,` `'``lon``'``]`
`return` `subset`
`aqs_sites` `=` `(``aqs_sites_full`
`.``pipe``(``rollup_dup_sites``)`
`.``pipe``(``cols_aqs``)``)`
现在aqs_sites数据框已准备就绪,我们转向 PurpleAir 站点。
整理 PurpleAir 站点列表
不同于 AQS 站点,包含 PurpleAir 传感器数据的文件是以 JSON 格式提供的。我们将在第十四章中更详细地讨论这种格式。现在,我们使用 shell 工具(参见第八章)来查看文件内容:
`!`head data/list_of_purpleair_sensors.json `|` cut -c `1`-60
{"version":"7.0.30",
"fields":
["ID","pm","pm_cf_1","pm_atm","age","pm_0","pm_1","pm_2","pm
"data":[
[20,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,97,0.0,0.0,0.0
[47,null,null,null,4951,null,null,null,null,null,null,null,9
[53,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,1.2,5.2,6.0,97,0.0,0.5,702
[74,0.0,0.0,0.0,0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,97,0.0,0.0,0.0
[77,9.8,9.8,9.8,1,9.8,10.7,11.0,11.2,13.8,15.1,15.5,97,9.7,9
[81,6.5,6.5,6.5,0,6.5,6.1,6.1,6.6,8.1,8.3,9.7,97,5.9,6.8,405
从文件的前几行可以猜测数据存储在"data"键中,列标签存储在"fields"键中。我们可以使用 Python 的json库将文件读取为 Python 的dict:
`import` `json`
`with` `open``(``'``data/list_of_purpleair_sensors.json``'``)` `as` `f``:`
`pa_json` `=` `json``.``load``(``f``)`
`list``(``pa_json``.``keys``(``)``)`
['version', 'fields', 'data', 'count']
我们可以从data中的值创建一个数据框,并使用fields的内容标记列:
`pa_sites_full` `=` `pd``.``DataFrame``(``pa_json``[``'``data``'``]``,` `columns``=``pa_json``[``'``fields``'``]``)`
`pa_sites_full``.``head``(``)`
| ID | pm | pm_cf_1 | pm_atm | ... | Voc | Ozone1 | Adc | CH | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 20 | 0.0 | 0.0 | 0.0 | ... | NaN | NaN | 0.01 | 1 |
| 1 | 47 | NaN | NaN | NaN | ... | NaN | 0.72 | 0.72 | 0 |
| 2 | 53 | 0.0 | 0.0 | 0.0 | ... | NaN | NaN | 0.00 | 1 |
| 3 | 74 | 0.0 | 0.0 | 0.0 | ... | NaN | NaN | 0.05 | 1 |
| 4 | 77 | 9.8 | 9.8 | 9.8 | ... | NaN | NaN | 0.01 | 1 |
5 rows × 36 columns
与 AQS 数据类似,此数据框中的列比我们需要的多得多:
`pa_sites_full``.``columns`
Index(['ID', 'pm', 'pm_cf_1', 'pm_atm', 'age', 'pm_0', 'pm_1', 'pm_2', 'pm_3',
'pm_4', 'pm_5', 'pm_6', 'conf', 'pm1', 'pm_10', 'p1', 'p2', 'p3', 'p4',
'p5', 'p6', 'Humidity', 'Temperature', 'Pressure', 'Elevation', 'Type',
'Label', 'Lat', 'Lon', 'Icon', 'isOwner', 'Flags', 'Voc', 'Ozone1',
'Adc', 'CH'],
dtype='object')
在这种情况下,我们可以猜测我们最感兴趣的列是传感器 ID (ID)、传感器标签 (Label)、纬度 (Lat)和经度 (Lon)。但我们确实在 PurpleAir 网站的数据字典中进行了查阅以进行双重检查。
现在让我们检查ID列是否存在重复,就像我们对 AQS 数据所做的那样:
`pa_sites_full``[``'``ID``'``]``.``value_counts``(``)``[``:``3``]`
85829 1
117575 1
118195 1
Name: ID, dtype: int64
由于 value_counts() 方法按降序列出计数,我们可以看到每个 ID 仅包含一次。因此,我们已验证粒度处于单个传感器级别。接下来,我们保留仅需要的列以匹配来自两个源的传感器位置:
`def` `cols_pa``(``df``)``:`
`subset` `=` `df``[``[``'``ID``'``,` `'``Label``'``,` `'``Lat``'``,` `'``Lon``'``]``]`
`subset``.``columns` `=` `[``'``id``'``,` `'``label``'``,` `'``lat``'``,` `'``lon``'``]`
`return` `subset`
`pa_sites` `=` `(``pa_sites_full`
`.``pipe``(``cols_pa``)``)`
`pa_sites``.``shape`
(23138, 4)
注意,PurpleAir 传感器比 AQS 传感器多数万台。我们的下一个任务是找到每个 AQS 传感器附近的 PurpleAir 传感器。
匹配 AQS 和 PurpleAir 传感器
我们的目标是通过找到每个 AQS 仪器附近的 PurpleAir 传感器来匹配两个数据框中的传感器。我们认为附近意味着在 50 米范围内。这种匹配比我们到目前为止见过的连接更具挑战性。例如,使用pandas的merge方法的天真方法会失败:
`aqs_sites``.``merge``(``pa_sites``,` `left_on``=``[``'``lat``'``,` `'``lon``'``]``,` `right_on``=``[``'``lat``'``,` `'``lon``'``]``)`
| site_id | lat | lon | id | label | |
|---|---|---|---|---|---|
| 0 | 06-111-1004 | 34.45 | -119.23 | 48393 | VCAPCD OJ |
我们不能简单地匹配具有完全相同纬度和经度的仪器;我们需要找到足够接近 AQS 仪器的 PurpleAir 站点。
要找出两个位置之间有多远,我们使用一个基本的近似方法:在南北方向上,111,111米大致相当于一度纬度,在东西方向上,111,111 * cos(纬度)相当于一度经度。^(1) 因此,我们可以找到对应于每个点周围 50 米×50 米矩形的纬度和经度范围:
`magic_meters_per_lat` `=` `111_111`
`offset_in_m` `=` `25`
`offset_in_lat` `=` `offset_in_m` `/` `magic_meters_per_lat`
`offset_in_lat`
0.000225000225000225
为了进一步简化,我们使用 AQS 站点的中位纬度:
`median_latitude` `=` `aqs_sites``[``'``lat``'``]``.``median``(``)`
`magic_meters_per_lon` `=` `111_111` `*` `np``.``cos``(``np``.``radians``(``median_latitude``)``)`
`offset_in_lon` `=` `offset_in_m` `/` `magic_meters_per_lon`
`offset_in_lon`
0.000291515219937587
现在我们可以将坐标匹配到offset_in_lat和offset_in_lon之内。在 SQL 中执行此操作比在pandas中要容易得多,因此我们将表推入临时 SQLite 数据库,然后运行查询以将表读回到数据框中:
`import` `sqlalchemy`
`db` `=` `sqlalchemy``.``create_engine``(``'``sqlite://``'``)`
`aqs_sites``.``to_sql``(``name``=``'``aqs``'``,` `con``=``db``,` `index``=``False``)`
`pa_sites``.``to_sql``(``name``=``'``pa``'``,` `con``=``db``,` `index``=``False``)`
`query` `=` `f``'''`
`SELECT`
`aqs.site_id AS aqs_id,`
`pa.id AS pa_id,`
`pa.label AS pa_label,`
`aqs.lat AS aqs_lat,`
`aqs.lon AS aqs_lon,`
`pa.lat AS pa_lat,`
`pa.lon AS pa_lon`
`FROM aqs JOIN pa`
`ON pa.lat -` `{``offset_in_lat``}` `<= aqs.lat`
`AND aqs.lat <= pa.lat +` `{``offset_in_lat``}`
`AND pa.lon -` `{``offset_in_lon``}` `<= aqs.lon`
`AND aqs.lon <= pa.lon +` `{``offset_in_lon``}`
`'''`
`matched` `=` `pd``.``read_sql``(``query``,` `db``)`
`matched`
| aqs_id | pa_id | pa_label | aqs_lat | aqs_lon | pa_lat | pa_lon | |
|---|---|---|---|---|---|---|---|
| 0 | 06-019-0011 | 6568 | IMPROVE_FRES2 | 36.79 | -119.77 | 36.79 | -119.77 |
| 1 | 06-019-0011 | 13485 | AMTS_Fresno | 36.79 | -119.77 | 36.79 | -119.77 |
| 2 | 06-019-0011 | 44427 | Fresno CARB CCAC | 36.79 | -119.77 | 36.79 | -119.77 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 146 | 53-061-1007 | 3659 | Marysville 7th | 48.05 | -122.17 | 48.05 | -122.17 |
| 147 | 53-063-0021 | 54603 | Augusta 1 SRCAA | 47.67 | -117.36 | 47.67 | -117.36 |
| 148 | 56-021-0100 | 50045 | WDEQ-AQD Cheyenne NCore | 41.18 | -104.78 | 41.18 | -104.78 |
149 rows × 7 columns
我们已经达到了我们的目标——我们匹配了 149 个 AQS 站点与 PurpleAir 传感器。我们对位置的整理完成了,现在转向整理和清理传感器测量数据的任务。我们从 AQS 站点获取的测量值开始。
整理和清理 AQS 传感器数据
现在我们已经找到彼此附近的传感器,我们准备整理和清理这些站点的测量数据文件。我们演示了一个 AQS 仪器及其匹配的 PurpleAir 传感器涉及的任务。我们选择了位于加利福尼亚州萨克拉门托的一对。AQS 传感器 ID 是 06-067-0010,PurpleAir 传感器名称是 AMTS_TESTINGA。
AQS 提供了一个网站和 API 用于下载传感器数据。我们下载了从 2018 年 5 月 20 日到 2019 年 12 月 29 日的每日测量数据,保存在 data/aqs_06-067-0010.csv 文件中。让我们从加载这个文件到数据框架开始:
`aqs_full` `=` `pd``.``read_csv``(``'``data/aqs_06-067-0010.csv``'``)`
`aqs_full``.``shape`
(2268, 31)
从 数据字典 中我们发现,arithmetic_mean 列对应于实际的 PM2.5 测量。一些 AQS 传感器每小时进行一次测量。对于我们的分析,我们下载了每日平均值(算术平均值)的 24 小时平均值的传感器测量。
让我们进行一些质量检查和必要的数据清理。我们关注与值的范围和质量相关的检查:
-
检查和修正数据的粒度。
-
删除不需要的列。
-
检查
date_local列中的值。 -
检查
arithmetic_mean列中的值。
为简洁起见,我们选择了一些重要的质量检查,特别是加强了我们在数据整理、探索性数据分析和可视化中涵盖的思想。
检查粒度
我们希望我们数据的每一行对应于单个日期,具有该日期的平均 PM2.5 读数。正如我们之前看到的,一个简单的检查方法是查看 date_local 列中是否有重复值:
`aqs_full``[``'``date_local``'``]``.``value_counts``(``)`
date_local
2019-01-03 12
2018-12-31 12
2018-12-28 12
..
2018-11-28 12
2018-11-25 12
2018-11-22 12
Name: count, Length: 189, dtype: int64
实际上,每个日期有 12 行数据,因此粒度并不是在个体日期级别上。
从数据字典中我们得知,有多种标准用于从原始传感器数据计算最终测量结果。pollutant_standard 列包含每个标准的名称。event_type 列标记了是否包括“异常事件”期间测得的数据。让我们通过计算这 12 个测量的范围来检查这些平均值有多么不同:
`(``aqs_full`
`.``groupby``(``'``date_local``'``)`
`[``'``arithmetic_mean``'``]`
`.``agg``(``np``.``ptp``)` `# np.ptp computes max() - min()`
`.``value_counts``(``)`
`)`
arithmetic_mean
0.0 189
Name: count, dtype: int64
对于所有 189 个日期,最大 PM2.5–最小 PM2.5 为 0。这意味着我们只需取每个日期的第一个 PM2.5 测量值:
`def` `rollup_dates``(``df``)``:`
`return` `(`
`df``.``groupby``(``'``date_local``'``)`
`.``first``(``)`
`.``reset_index``(``)`
`)`
`aqs` `=` `(``aqs_full`
`.``pipe``(``rollup_dates``)``)`
`aqs``.``shape`
(189, 31)
此数据清理步骤使我们获得了所需的细粒度:每一行代表一个单独的日期,具有该日期的平均 PM2.5 测量值。接下来,我们进一步修改数据框的结构并删除不需要的列。
移除不必要的列
我们计划将 AQS 数据框中的 PM2.5 测量与 PurpleAir 的 PM2.5 测量进行匹配,以每日为单位。为简化结构,我们可以舍弃除日期和 PM2.5 列外的所有列。我们还将 PM2.5 列重命名,以便更容易理解:
`def` `drop_cols``(``df``)``:`
`subset` `=` `df``[``[``'``date_local``'``,` `'``arithmetic_mean``'``]``]`
`return` `subset``.``rename``(``columns``=``{``'``arithmetic_mean``'``:` `'``pm25``'``}``)`
`aqs` `=` `(``aqs_full`
`.``pipe``(``rollup_dates``)`
`.``pipe``(``drop_cols``)``)`
`aqs``.``head``(``)`
| date_local | pm25 | |
|---|---|---|
| 0 | 2018-05-20 | 6.5 |
| 1 | 2018-05-23 | 2.3 |
| 2 | 2018-05-29 | 11.8 |
| 3 | 2018-06-01 | 6.0 |
| 4 | 2018-06-04 | 8.0 |
现在我们已经为数据表获得了所需的形状,我们转而检查数据值。
检查日期的有效性
让我们仔细查看日期。我们已经看到当没有 PM2.5 读数时存在间隙,因此我们预计存在缺失日期。让我们将日期解析为时间戳对象,以便更容易确定缺失的日期。与我们在 第九章 中所做的一样,我们检查格式:
`aqs``[``'``date_local``'``]``.``iloc``[``:``3``]`
0 2018-05-20
1 2018-05-23
2 2018-05-29
Name: date_local, dtype: object
日期以 YYYY-MM-DD 表示,因此我们在 Python 表示中描述为 '%Y-%m-%d'。为了解析日期,我们使用 pd.to_datetime() 函数,并将 date_local 列重新分配为 pd.TimeStamp:
`def` `parse_dates``(``df``)``:`
`date_format` `=` `'``%Y``-``%m``-``%d``'`
`timestamps` `=` `pd``.``to_datetime``(``df``[``'``date_local``'``]``,` `format``=``date_format``)`
`return` `df``.``assign``(``date_local``=``timestamps``)`
`aqs` `=` `(``aqs_full`
`.``pipe``(``rollup_dates``)`
`.``pipe``(``drop_cols``)`
`.``pipe``(``parse_dates``)``)`
该方法运行无误,表明所有字符串都匹配了格式。
注意
仅仅因为日期可以解析,这并不意味着日期立即可以用于进一步的分析。例如,字符串 9999-01-31 可以解析为 pd.TimeStamp,但该日期无效。
现在日期已转换为时间戳,我们可以计算缺失的日期数。我们找到最早日期和最晚日期之间的天数——这对应于我们可能记录的最大测量次数:
`date_range` `=` `aqs``[``'``date_local``'``]``.``max``(``)` `-` `aqs``[``'``date_local``'``]``.``min``(``)`
`date_range``.``days`
588
减去时间戳会得到 Timedelta 对象,正如我们所看到的,它们具有一些有用的属性。数据中有很多日期缺失。然而,当我们将此传感器的这些数据与其他传感器的数据结合起来时,我们希望有足够的数据来拟合一个模型。
我们的最终整理步骤是检查 PM2.5 测量的质量。
检查 PM2.5 测量质量
颗粒物浓度以每立方米空气中的微克数(µg/m³)来衡量。(1 克中有 1 百万微克,1 磅约等于 450 克。)EPA 制定了标准,即 PM2.5 的日均值为 35 µg/m³,年均值为 12 µg/m³。我们可以利用这些信息对 PM2.5 测量进行几项基本检查。首先,PM2.5 不能低于 0。其次,我们可以寻找异常高的 PM2.5 值,并查看它们是否对应于重大事件,如野火。
进行这些检查的一种视觉方法是将 PM2.5 测量值绘制成日期图:
`px``.``scatter``(``aqs``,` `x``=``'``date_local``'``,` `y``=``'``pm25``'``,`
`labels``=``{``'``date_local``'``:``'``Date``'``,` `'``pm25``'``:``'``AQS daily avg PM2.5``'``}``,`
`width``=``500``,` `height``=``250``)`
我们发现 PM2.5 的测量数值不会低于 0,通常低于 EPA 的标准。我们还发现 2018 年 11 月中旬 PM2.5 出现了大幅上升。该传感器位于萨克拉门托,因此我们可以检查该地区是否发生了火灾。
实际上,2018 年 11 月 8 日标志着“加州历史上最致命和破坏性的野火”——“大火营”(见由美国人口普查局管理的大火营页面)。火灾发生在距离萨克拉门托仅 80 英里的地方,因此这个 AQS 传感器捕捉到了 PM2.5 的剧烈增长。
我们已经清理和探索了一个 AQS 传感器的数据。在下一节中,我们将对其附近的 PurpleAir 传感器进行相同的操作。
整理 PurpleAir 传感器数据
在上一节中,我们分析了 AQS 站点 06-067-0010 的数据。匹配的 PurpleAir 传感器命名为 AMTS_TESTINGA,我们已经使用 PurpleAir 网站将此传感器的数据下载到 data/purpleair_AMTS 文件夹中:
`!`ls -alh data/purpleair_AMTS/* `|` cut -c `1`-72
-rw-r--r-- 1 nolan staff 50M Jan 25 16:35 data/purpleair_AMTS/AMTS_
-rw-r--r-- 1 nolan staff 50M Jan 25 16:35 data/purpleair_AMTS/AMTS_
-rw-r--r-- 1 nolan staff 48M Jan 25 16:35 data/purpleair_AMTS/AMTS_
-rw-r--r-- 1 nolan staff 50M Jan 25 16:35 data/purpleair_AMTS/AMTS_
有四个 CSV 文件。它们的名称相当长,并且每个的开头都相同。PurpleAir 数据的数据字典表明,每个传感器都有两个单独的仪器 A 和 B,每个都记录数据。请注意,我们用于收集这些数据和配套数据字典的 PurpleAir 网站已降级。数据现在通过 REST API 可用。记录 API 的站点还包含有关字段的信息。(REST 的主题在第十四章中涵盖。)让我们检查文件名的后部分:
`!`ls -alh data/purpleair_AMTS/* `|` cut -c `73`-140
TESTING (outside) (38.568404 -121.493163) Primary Real Time 05_20_20
TESTING (outside) (38.568404 -121.493163) Secondary Real Time 05_20_
TESTING B (undefined) (38.568404 -121.493163) Primary Real Time 05_2
TESTING B (undefined) (38.568404 -121.493163) Secondary Real Time 05
我们可以看到前两个 CSV 文件对应于仪器 A,最后两个对应于 B。拥有两个仪器对数据清理很有用;如果 A 和 B 对某个测量结果有异议,我们可能会质疑测量结果的完整性,并决定删除它。
数据字典还提到,每个仪器记录主要和次要数据。主要数据包含我们感兴趣的字段:PM2.5、温度和湿度。次要数据包含其他粒径的数据,如 PM1.0 和 PM10。因此我们只使用主文件。
我们的任务与上一节类似,只是增加了处理两台仪器读数的内容。
我们首先加载数据。当 CSV 文件名很长时,我们可以将文件名分配给 Python 变量,以更轻松地加载文件:
`from` `pathlib` `import` `Path`
`data_folder` `=` `Path``(``'``data/purpleair_AMTS``'``)`
`pa_csvs` `=` `sorted``(``data_folder``.``glob``(``'``*.csv``'``)``)`
`pa_csvs``[``0``]`
PosixPath('data/purpleair_AMTS/AMTS_TESTING (outside) (38.568404 -121.493163) Primary Real Time 05_20_2018 12_29_2019.csv')
`pa_full` `=` `pd``.``read_csv``(``pa_csvs``[``0``]``)`
`pa_full``.``shape`
(672755, 11)
让我们看一下列,确定我们需要哪些列:
`pa_full``.``columns`
Index(['created_at', 'entry_id', 'PM1.0_CF1_ug/m3', 'PM2.5_CF1_ug/m3',
'PM10.0_CF1_ug/m3', 'UptimeMinutes', 'RSSI_dbm', 'Temperature_F',
'Humidity_%', 'PM2.5_ATM_ug/m3', 'Unnamed: 10'],
dtype='object')
虽然我们对 PM2.5 感兴趣,但似乎有两列包含 PM2.5 数据:PM2.5_CF1_ug/m3 和 PM2.5_ATM_ug/m3。我们调查了这两列之间的差异,发现 PurpleAir 传感器使用两种不同的方法将原始激光记录转换为 PM2.5 数值。这两种计算对应于 CF1 和 ATM 列。Barkjohn 发现使用 CF1 比 ATM 产生了更好的结果,因此我们保留该列,以及日期、温度和相对湿度:
`def` `drop_and_rename_cols``(``df``)``:`
`df` `=` `df``[``[``'``created_at``'``,` `'``PM2.5_CF1_ug/m3``'``,` `'``Temperature_F``'``,` `'``Humidity_``%``'``]``]`
`df``.``columns` `=` `[``'``timestamp``'``,` `'``PM25cf1``'``,` `'``TempF``'``,` `'``RH``'``]`
`return` `df`
`pa` `=` `(``pa_full`
`.``pipe``(``drop_and_rename_cols``)``)`
`pa``.``head``(``)`
| 时间戳 | PM25cf1 | TempF | RH | |
|---|---|---|---|---|
| 0 | 2018-05-20 00:00:35 UTC | 1.23 | 83.0 | 32.0 |
| 1 | 2018-05-20 00:01:55 UTC | 1.94 | 83.0 | 32.0 |
| 2 | 2018-05-20 00:03:15 UTC | 1.80 | 83.0 | 32.0 |
| 3 | 2018-05-20 00:04:35 UTC | 1.64 | 83.0 | 32.0 |
| 4 | 2018-05-20 00:05:55 UTC | 1.33 | 83.0 | 32.0 |
接下来我们检查粒度。
检查粒度
为了使这些测量的粒度与 AQS 数据匹配,我们希望每个日期(24 小时)有一个平均 PM2.5。PurpleAir 表示传感器每两分钟进行一次测量。在我们将其聚合到 24 小时期间之前,让我们再次检查原始测量的粒度。
要实现这一点,我们将包含日期信息的列从字符串转换为 pd.TimeStamp 对象。日期的格式与 AQS 格式不同,我们描述为 '%Y-%m-%d %X %Z'。正如我们所见,pandas 对于具有时间戳索引的数据框架有特殊支持:
`def` `parse_timestamps``(``df``)``:`
`date_format` `=` `'``%Y``-``%m``-``%d` `%X` `%``Z``'`
`times` `=` `pd``.``to_datetime``(``df``[``'``timestamp``'``]``,` `format``=``date_format``)`
`return` `(``df``.``assign``(``timestamp``=``times``)`
`.``set_index``(``'``timestamp``'``)``)`
`pa` `=` `(``pa_full`
`.``pipe``(``drop_and_rename_cols``)`
`.``pipe``(``parse_timestamps``)``)`
`pa``.``head``(``2``)`
| PM25cf1 | TempF | RH | |
|---|---|---|---|
| 时间戳 | |||
| --- | --- | --- | --- |
| 2018-05-20 00:00:35+00:00 | 1.23 | 83.0 | 32.0 |
| 2018-05-20 00:01:55+00:00 | 1.94 | 83.0 | 32.0 |
时间戳很棘手 — 注意原始时间戳是以 UTC 时区给出的。然而,AQS 数据根据 加利福尼亚州当地时间 平均,这比 UTC 时间晚七或八小时,这取决于是否实行夏令时。这意味着我们需要更改 PurpleAir 时间戳的时区以匹配当地时区。df.tz_convert() 方法作用于数据框架的索引,这也是我们将 pa 的索引设置为时间戳的一个原因:
`def` `convert_tz``(``pa``)``:`
`return` `pa``.``tz_convert``(``'``US/Pacific``'``)`
`pa` `=` `(``pa_full`
`.``pipe``(``drop_and_rename_cols``)`
`.``pipe``(``parse_timestamps``)`
`.``pipe``(``convert_tz``)``)`
`pa``.``head``(``2``)`
| PM25cf1 | TempF | RH | |
|---|---|---|---|
| 时间戳 | |||
| --- | --- | --- | --- |
| 2018-05-19 17:00:35-07:00 | 1.23 | 83.0 | 32.0 |
| 2018-05-19 17:01:55-07:00 | 1.94 | 83.0 | 32.0 |
如果我们将此版本的数据框架的前两行与上一个进行比较,我们会看到时间已更改以表明与 UTC 相差七个小时。
可视化时间戳可以帮助我们检查数据的粒度。
可视化时间戳
可视化时间戳的一种方法是计算每个 24 小时内出现的次数,然后绘制这些计数随时间的变化。要在 pandas 中对时间序列数据进行分组,可以使用 df.resample() 方法。此方法适用于具有时间戳索引的数据框。它的行为类似于 df.groupby(),但我们可以指定希望如何分组时间戳——可以分组为日期、周、月等多种选项(D 参数告诉 resample 将时间戳聚合成单独的日期):
`per_day` `=` `(``pa``.``resample``(``'``D``'``)`
`.``size``(``)`
`.``rename``(``'``records_per_day``'``)`
`.``to_frame``(``)`
`)`
`percs` `=` `[``10``,` `25``,` `50``,` `75``,` `100``]`
`np``.``percentile``(``per_day``[``'``records_per_day``'``]``,` `percs``,` `method``=``'``lower``'``)`
array([ 293, 720, 1075, 1440, 2250])
我们可以看到每天的测量次数变化很大。这些计数的折线图能更好地展示这些变化:
`px``.``line``(``per_day``,` `x``=``per_day``.``index``,` `y``=``'``records_per_day``'``,`
`labels``=``{``'``timestamp``'``:``'``Date``'``,` `'``records_per_day``'``:``'``Records per day``'``}``,`
`width``=``550``,` `height``=``250``,``)`
这是一个引人入胜的图表。我们可以看到数据中存在明显的缺失测量间隙。在 2018 年 7 月和 2019 年 9 月,数据的大部分似乎都丢失了。即使传感器运行正常,每天的测量次数也略有不同。例如,在 2018 年 8 月至 10 月之间的折线图上“崎岖不平”,日期上的测量次数有所变化。我们需要决定如何处理缺失数据。但也许更紧迫的是:折线图上出现了奇怪的“阶梯”。一些日期大约有 1,000 个读数,一些大约有 2,000 个,一些大约有 700 个,一些大约有 1,400 个。如果传感器每两分钟进行一次测量,则每天的最大测量次数应为 720 次。对于完美的传感器,折线图应显示为 720 次的平直线。显然情况并非如此。让我们来调查一下。
检查采样率
进一步的挖掘显示,尽管 PurpleAir 传感器当前每 120 秒记录一次数据,但以前并非如此。2019 年 5 月 30 日之前,传感器每 80 秒记录一次数据,即每天 1,080 个数据点。采样率的变化确实解释了 2019 年 5 月 30 日的数据下降。接下来我们看看时间段内存在远高于预期点数的情况。这可能意味着数据中存在重复的测量。我们可以通过查看一天的测量数据来验证这一点,例如 2019 年 1 月 1 日。我们通过向 .loc 传入字符串来筛选该日期的时间戳:
`len``(``pa``.``loc``[``'``2019-01-01``'``]``)`
2154
几乎有 1,080 个预期读数的两倍。让我们检查是否存在重复读数:
`pa``.``loc``[``'``2019-01-01``'``]``.``index``.``value_counts``(``)`
2019-01-01 13:52:30-08:00 2
2019-01-01 12:02:21-08:00 2
2019-01-01 11:49:01-08:00 2
..
2019-01-01 21:34:10-08:00 2
2019-01-01 11:03:41-08:00 2
2019-01-01 04:05:38-08:00 2
Name: timestamp, Length: 1077, dtype: int64
每个时间戳正好出现两次,我们可以验证所有重复的日期包含相同的 PM2.5 读数。由于温度和湿度也是如此,我们从数据框中删除重复行:
`def` `drop_duplicate_rows``(``df``)``:`
`return` `df``[``~``df``.``index``.``duplicated``(``)``]`
`pa` `=` `(``pa_full`
`.``pipe``(``drop_and_rename_cols``)`
`.``pipe``(``parse_timestamps``)`
`.``pipe``(``convert_tz``)`
`.``pipe``(``drop_duplicate_rows``)``)`
`pa``.``shape`
(502628, 3)
为了检查,我们重新制作了每天记录数的折线图,这次我们会在预期包含的区域内填充颜色:
`per_day` `=` `(``pa``.``resample``(``'``D``'``)`
`.``size``(``)``.``rename``(``'``records_per_day``'``)`
`.``to_frame``(``)`
`)`
`fig` `=` `px``.``line``(``per_day``,` `x``=``per_day``.``index``,` `y``=``'``records_per_day``'``,`
`labels``=``{``'``timestamp``'``:``'``Date``'``,` `'``records_per_day``'``:``'``Records per day``'``}``,`
`width``=``550``,` `height``=``250``)`
`fig``.``add_annotation``(``x``=``'``2019-07-24``'``,` `y``=``720``,`
`text``=``"``720``"``,` `showarrow``=``False``,` `yshift``=``10``)`
`fig``.``add_annotation``(``x``=``'``2019-07-24``'``,` `y``=``1080``,`
`text``=``"``1080``"``,` `showarrow``=``False``,` `yshift``=``10``)`
`fig``.``add_hline``(``y``=``1080``,` `line_width``=``3``,` `line_dash``=``"``dot``"``,` `opacity``=``0.6``)`
`fig``.``add_hline``(``y``=``720``,` `line_width``=``3``,` `line_dash``=``"``dot``"``,` `opacity``=``0.6``)`
`fig``.``add_vline``(``x``=``"``2019-05-30``"``,` `line_width``=``3``,` `line_dash``=``"``dash``"``,` `opacity``=``0.6``)`
`fig`
去除重复日期后,每天的测量数据图与我们预期的计数更一致。细心的读者会发现每年 11 月左右,在取消夏令时后,有两次高于最大测量值的峰值。当时钟被倒退一小时时,那一天有 25 小时而不是通常的 24 小时。时间戳确实很棘手!
但仍然有缺失的测量值,我们需要决定如何处理。
处理缺失值
计划是创建测量值的 24 小时平均,但我们不想使用测量不足的日子。我们遵循 Barkjohn 的分析,仅在该日有至少 90%可能数据点时保留 24 小时平均。请记住,在 2019 年 5 月 30 日之前,一天有 1080 个可能数据点,之后为 720 个可能数据点。我们计算了每天需要保留的最少测量数:
`needed_measurements_80s` `=` `0.9` `*` `1080`
`needed_measurements_120s` `=` `0.9` `*` `720`
现在我们可以确定哪些日子有足够的测量数据可以保留:
`cutoff_date` `=` `pd``.``Timestamp``(``'``2019-05-30``'``,` `tz``=``'``US/Pacific``'``)`
`def` `has_enough_readings``(``one_day``)``:`
`[``n``]` `=` `one_day`
`date` `=` `one_day``.``name`
`return` `(``n` `>``=` `needed_measurements_80s`
`if` `date` `<``=` `cutoff_date`
`else` `n` `>``=` `needed_measurements_120s``)`
`should_keep` `=` `per_day``.``apply``(``has_enough_readings``,` `axis``=``'``columns``'``)`
`should_keep``.``head``(``)`
timestamp
2018-05-19 00:00:00-07:00 False
2018-05-20 00:00:00-07:00 True
2018-05-21 00:00:00-07:00 True
2018-05-22 00:00:00-07:00 True
2018-05-23 00:00:00-07:00 True
Freq: D, dtype: bool
我们已经准备好将每天的读数进行平均,并删除不足的天数:
`def` `compute_daily_avgs``(``pa``)``:`
`should_keep` `=` `(``pa``.``resample``(``'``D``'``)`
`[``'``PM25cf1``'``]`
`.``size``(``)`
`.``to_frame``(``)`
`.``apply``(``has_enough_readings``,` `axis``=``'``columns``'``)``)`
`return` `(``pa``.``resample``(``'``D``'``)`
`.``mean``(``)`
`.``loc``[``should_keep``]``)`
`pa` `=` `(``pa_full`
`.``pipe``(``drop_and_rename_cols``)`
`.``pipe``(``parse_timestamps``)`
`.``pipe``(``convert_tz``)`
`.``pipe``(``drop_duplicate_rows``)`
`.``pipe``(``compute_daily_avgs``)``)`
`pa``.``head``(``2``)`
| PM25cf1 | TempF | RH | |
|---|---|---|---|
| timestamp | |||
| --- | --- | --- | --- |
| 2018-05-20 00:00:00-07:00 | 2.48 | 83.35 | 28.72 |
| 2018-05-21 00:00:00-07:00 | 3.00 | 83.25 | 29.91 |
现在我们有了仪器 A 的平均每日 PM2.5 读数,需要在仪器 B 上重复刚才在仪器 A 上执行的数据整理工作。幸运的是,我们可以重复使用同样的流程。为了简洁起见,我们不在此处包含该数据整理步骤。但我们需要决定如果 PM2.5 的平均值不同要怎么办。如果 A 和 B 的 PM2.5 值相差超过 61%,或者超过 5 µg m⁻³,Barkjohn 会删除 12 行中的一对传感器中的此行。
如您所见,准备和清理这些数据需要大量工作:我们处理了缺失数据,对每个仪器的读数进行了聚合,将两个仪器的读数平均在一起,并删除了它们不一致的行。这项工作使我们对 PM2.5 读数更加自信。我们知道最终数据框中的每个 PM2.5 值都是来自生成一致且完整读数的两个独立仪器的日均值。
要完全复制 Barkjohn 的分析,我们需要在所有 PurpleAir 传感器上重复此过程。然后我们会在所有 AQS 传感器上重复 AQS 清洁程序。最后,我们将 PurpleAir 和 AQS 数据合并在一起。这个过程产生了每对同位置传感器的日均读数。为了简洁起见,我们略过此代码。而是继续使用小组数据集进行分析的最后步骤。我们从 EDA 开始,着眼于建模。
探索 PurpleAir 和 AQS 的测量数据
让我们探索匹配的 AQS 和 PurpleAir PM2.5 读数的清理数据集,并寻找可能帮助我们建模的见解。我们主要关注两种空气质量测量数据之间的关系。但我们要记住数据的范围,比如这些数据在时间和地点上的分布。从数据清理中我们了解到,我们处理的是几年内的 PM2.5 日均值数据,并且我们有来自美国数十个地点的数据。
首先,我们回顾整个清理过的数据框:
`full_df`
| date | id | region | pm25aqs | pm25pa | temp | rh | dew | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2019-05-17 | AK1 | 阿拉斯加 | 6.7 | 8.62 | 18.03 | 38.56 | 3.63 |
| 1 | 2019-05-18 | AK1 | 阿拉斯加 | 3.8 | 3.49 | 16.12 | 49.40 | 5.44 |
| 2 | 2019-05-21 | AK1 | 阿拉斯加 | 4.0 | 3.80 | 19.90 | 29.97 | 1.73 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 12427 | 2019-02-20 | WI6 | 北部 | 15.6 | 25.30 | 1.71 | 65.78 | -4.08 |
| 12428 | 2019-03-04 | WI6 | 北部 | 14.0 | 8.21 | -14.38 | 48.21 | -23.02 |
| 12429 | 2019-03-22 | WI6 | 北部 | 5.8 | 9.44 | 5.08 | 52.20 | -4.02 |
12246 rows × 8 columns
我们在下表中为数据框中的每一列提供解释:
| 列名 | 描述 |
|---|---|
| date | 观测日期 |
| id | 一个唯一标识符,格式为带有数字的美国州缩写(我们为站点 ID CA1 执行了数据清理) |
| region | 地区的名称,对应一组站点(CA1站点位于西部地区) |
| pm25aqs | AQS 传感器测量的 PM2.5 |
| pm25pa | 紫外传感器测量的 PM2.5 |
| temp | 温度,摄氏度 |
| rh | 相对湿度,范围从 0%到 100% |
| dew | 露点(较高的露点意味着空气中有更多的湿气) |
让我们从制作几个简单的可视化图开始获取见解。由于范围涉及特定位置的时间测量,我们可以选择一个有很多测量记录的位置,并制作每周平均空气质量的折线图。首先,让我们找出记录较多的站点:
`full_df``[``'``id``'``]``.``value_counts``(``)``[``:``3``]`
id
IA3 830
NC4 699
CA2 659
Name: count, dtype: int64
标记为NC4的位置有近 700 次观测记录。为了稍微平滑折线图,让我们制作每周平均值图:
`nc4` `=` `full_df``.``loc``[``full_df``[``'``id``'``]` `==``'``NC4``'``]`
`ts_nc4` `=` `(``nc4``.``set_index``(``'``date``'``)`
`.``resample``(``'``W``'``)`
`[``'``pm25aqs``'``,` `'``pm25pa``'``]`
`.``mean``(``)`
`.``reset_index``(``)`
`)`
`fig` `=` `px``.``line``(``ts_nc4``,` `x``=``'``date``'``,` `y``=``'``pm25aqs``'``,`
`labels``=``{``'``date``'``:``'``'``,` `'``pm25aqs``'``:``'``PM2.5 weekly avg``'``}``,`
`width``=``500``,` `height``=``250``)`
`fig``.``add_trace``(``go``.``Scatter``(``x``=``ts_nc4``[``'``date``'``]``,` `y``=``ts_nc4``[``'``pm25pa``'``]``,`
`line``=``dict``(``color``=``'``black``'``,` `dash``=``'``dot``'``)``)``)`
`fig``.``update_yaxes``(``range``=``[``0``,``30``]``)`
`fig``.``update_layout``(``showlegend``=``False``)`
`fig``.``show``(``)`
我们看到大多数 AQS 传感器(实线)测得的 PM2.5 值介于 5.0 到 15.0 µg m⁻³之间。PurpleAir 传感器跟随 AQS 传感器的上下波动模式,这让人感到安心。但测量值始终高于 AQS,并且在某些情况下相当高,这表明可能需要进行修正。
接下来,让我们考虑两个传感器 PM2.5 读数的分布情况:
`left` `=` `px``.``histogram``(``nc4``,` `x``=``'``pm25aqs``'``,` `histnorm``=``'``percent``'``)`
`right` `=` `px``.``histogram``(``nc4``,` `x``=``'``pm25pa``'``,` `histnorm``=``'``percent``'``)`
`fig` `=` `left_right``(``left``,` `right``,` `width``=``600``,` `height``=``250``)`
`fig``.``update_xaxes``(``title``=``'``AQS readings``'``,` `col``=``1``,` `row``=``1``)`
`fig``.``update_xaxes``(``title``=``'``PurpleAir readings``'``,` `col``=``2``,` `row``=``1``)`
`fig``.``update_yaxes``(``title``=``'``percent``'``,` `col``=``1``,` `row``=``1``)`
`fig``.``show``(``)`
这两个分布都是右偏的,这在值有下限时经常发生(在这种情况下是 0)。比较这两个分布的更好方法是使用量化-量化图(见 第十章)。使用 q-q 图可以更容易地比较均值、扩展和尾部。
`percs` `=` `np``.``arange``(``1``,` `100``,` `1``)`
`aqs_qs` `=` `np``.``percentile``(``nc4``[``'``pm25aqs``'``]``,` `percs``,` `interpolation``=``'``lower``'``)`
`pa_qs` `=` `np``.``percentile``(``nc4``[``'``pm25pa``'``]``,` `percs``,` `interpolation``=``'``lower``'``)`
`perc_df` `=` `pd``.``DataFrame``(``{``'``percentile``'``:` `percs``,` `'``aqs_qs``'``:``aqs_qs``,` `'``pa_qs``'``:``pa_qs``}``)`
`fig` `=` `px``.``scatter``(``perc_df``,` `x``=``'``aqs_qs``'``,` `y``=``'``pa_qs``'``,`
`labels``=``{``'``aqs_qs``'``:` `'``AQS quantiles``'``,`
`'``pa_qs``'``:` `'``PurpleAir quantiles``'``}``,`
`width``=``350``,` `height``=``250``)`
`fig``.``add_trace``(``go``.``Scatter``(``x``=``[``2``,` `13``]``,` `y``=``[``1``,` `25``]``,`
`mode``=``'``lines``'``,` `line``=``dict``(``dash``=``'``dash``'``,` `width``=``4``)``)``)`
`fig``.``update_layout``(``showlegend``=``False``)`
`fig`
量化-量化图大致是线性的。我们叠加了一个斜率为 2.2 的虚线,它很好地对齐了分位数,这表明 PurpleAir 测量的分布约是 AQS 的两倍。
在 q-q 图或并排直方图中看不到的是传感器读数如何一起变化。让我们来看看这个。首先,我们看一下两个读数之间的差异分布:
`diffs` `=` `(``nc4``[``'``pm25pa``'``]` `-` `nc4``[``'``pm25aqs``'``]``)`
`fig` `=` `px``.``histogram``(``diffs``,` `histnorm``=``'``percent``'``,`
`width``=``350``,` `height``=``250``)`
`fig``.``update_xaxes``(``range``=``[``-``10``,``30``]``,` `title``=``"``Difference: PA–AQS reading``"``)`
`fig``.``update_traces``(``xbins``=``dict``(`
`start``=``-``10.0``,` `end``=``30.0``,` `size``=``2``)``)`
`fig``.``update_layout``(``showlegend``=``False``)`
`fig``.``show``(``)`
如果仪器完全一致,我们会在 0 处看到一个尖峰。如果仪器一致且没有偏差的测量误差,我们期望看到以 0 为中心的分布。然而,我们看到 90% 的时间,PurpleAir 测量大于 AQS 的 24 小时平均值,大约 25% 的时间高出 10 µg/m³,考虑到 AQS 的平均值通常在 5 µg/m³ 到 10 µg/m³ 之间,这是很多的。
散点图可以让我们进一步了解这两台仪器的测量之间的关系。由于我们有兴趣找到一个总体关系,不论时间和地点,我们在图中包含所有的平均读数:
`px``.``scatter``(``full_df``,` `x``=``'``pm25aqs``'``,` `y``=``'``pm25pa``'``,` `width``=``350``,` `height``=``250``,`
`labels``=``{``'``pm25aqs``'``:``'``AQS PM2.5``'``,` `'``pm25pa``'``:``'``PurpleAir PM2.5``'``}``)`
虽然关系看起来是线性的,但除了少数几个读数外,其余都集中在图的左下角。让我们重新制作散点图并放大数据的大部分,以便更好地观察。我们还向图中添加了一个平滑曲线,以帮助更好地看到关系:
`full_df` `=` `full_df``.``loc``[``(``full_df``[``'``pm25aqs``'``]` `<` `50``)``]`
`px``.``scatter``(``full_df``,` `x``=``'``pm25aqs``'``,` `y``=``'``pm25pa``'``,`
`trendline``=``'``lowess``'``,` `trendline_color_override``=``"``orange``"``,`
`labels``=``{``'``pm25aqs``'``:``'``AQS PM2.5``'``,` `'``pm25pa``'``:``'``PurpleAir PM2.5``'``}``,`
`width``=``350``,` `height``=``250``)`
这种关系看起来大致是线性的,但在 AQS 的小值区间中曲线略微弯曲。当空气非常清洁时,PurpleAir 传感器不会检测到太多颗粒物质,因此更准确。此外,我们可以看到曲线应通过点 (0, 0)。尽管关系中有轻微弯曲,但这两个测量之间的线性关联(相关性)很高:
`np``.``corrcoef``(``full_df``[``'``pm25aqs``'``]``,` `full_df``[``'``pm25pa``'``]``)`
array([[1\. , 0.88],
[0.88, 1\. ]])
在开始这项分析之前,我们预期 PurpleAir 测量通常会高估 PM2.5。事实上,在散点图中反映了这一点,但我们也看到这两台仪器的测量之间似乎有强烈的线性关系,这将有助于校准 PurpleAir 传感器。
创建校正 PurpleAir 测量的模型
现在我们已经探讨了 AQS 和 PurpleAir 传感器之间 PM2.5 读数的关系,我们准备进行分析的最后一步:创建一个修正 PurpleAir 测量的模型。Barkjohn 的原始分析对数据拟合了许多模型,以找到最合适的模型。在本节中,我们使用来自第四章的技术来拟合一个简单的线性模型。我们还简要描述了 Barkjohn 为实际应用选择的最终模型。由于这些模型使用了本书后面介绍的方法,我们不会在这里深入解释技术细节。相反,我们鼓励您在阅读第十五章后重新访问本节。
首先,让我们总结一下我们的建模目标。我们希望创建一个尽可能准确预测 PM2.5 的模型。为此,我们建立了一个根据 AQS 测量调整 PurpleAir 测量的模型。我们将 AQS 测量视为真实的 PM2.5 值,因为它们来自精确校准的仪器,并且被美国政府积极用于决策制定。因此,我们有理由信任 AQS PM2.5 值的精度和接近真实的特性。
在我们建立了通过 AQS 调整 PurpleAir 测量的模型之后,我们将模型反转并使用它来预测未来的真实空气质量,这是一个校准场景。由于 AQS 测量接近真实值,我们将更为变化的 PurpleAir 测量与之对齐;这就是校准过程。然后我们使用校准曲线来纠正未来的 PurpleAir 测量。这个两步过程包含在即将介绍的简单线性模型及其反转形式中。
首先,我们拟合一条线来预测从真实数据(由 AQS 仪器记录)中记录的 PurpleAir(PA)测量:
PA ≈ b + m AQS
接下来,我们将线条反转,使用 PA 测量来预测空气质量:
True air quality ≈ − b / m + 1 / m PA
在我们的探索性数据分析期间制作的散点图和直方图表明 PurpleAir 测量更为变化,这支持校准方法。我们发现 PurpleAir 测量值大约是 AQS 测量值的两倍,这表明m可能接近 2,而1 / m可能接近1 / 2。
现在让我们拟合模型。根据第四章的概念,我们选择一个损失函数并最小化平均误差。回想一下,损失函数衡量我们的模型与实际数据的差距。我们使用平方损失,在这种情况下是[ P A − ( b + m A Q S ) ] 2 。为了使模型与我们的数据拟合,我们最小化数据上的平均平方损失:
1 n ∑ i = 1 n [ P A i − ( b + m A Q S i ] 2
我们使用scikit-learn提供的线性建模功能来做到这一点(现在不必担心这些细节):
`from` `sklearn``.``linear_model` `import` `LinearRegression`
`AQS``,` `PA` `=` `full_df``[``[``'``pm25aqs``'``]``]``,` `full_df``[``'``pm25pa``'``]`
`model` `=` `LinearRegression``(``)``.``fit``(``AQS``,` `PA``)`
`m``,` `b` `=` `model``.``coef_``[``0``]``,` `model``.``intercept_`
通过反转线路,我们得到估计:
`print``(``f``"``True air quality estimate =` `{``-``b``/``m``:``.2``}` `+` `{``1``/``m``:``.2``}``PA``"``)`
True air quality estimate = 1.4 + 0.53PA
这接近我们的预期。对 PurpleAir 测量的调整约为1 / 2 。
Barkjohn 确定的模型包含了相对湿度:
PA ≈ b + m 1 AQS + m 2 RH
这是一个多变量线性回归模型的示例—它使用多个变量进行预测。我们可以通过最小化数据上的平均平方误差来拟合它:
1 n ∑ i = 1 n [ P A i − ( b + m 1 A Q S i + m 2 R H i ] 2
接着我们反转校准,使用以下方程找到预测模型:
True air quality ≈ − b m 1 + 1 m 1 PA − m 2 m 1 RH
我们拟合这个模型并检查系数:
`AQS_RH``,` `PA` `=` `full_df``[``[``'``pm25aqs``'``,` `'``rh``'``]``]``,` `full_df``[``'``pm25pa``'``]`
`model_h` `=` `LinearRegression``(``)``.``fit``(``AQS_RH``,` `PA``)`
`[``m1``,` `m2``]``,` `b` `=` `model_h``.``coef_``,` `model_h``.``intercept_`
`print``(``f``"``True Air Quality Estimate =` `{``-``b``/``m``:``1.2``}` `+` `{``1``/``m1``:``.2``}``PA +` `{``-``m2``/``m1``:``.2``}``RH``"``)`
True Air Quality Estimate = 5.7 + 0.53PA + -0.088RH
在第 15 和 16 章中,我们将学习如何通过检查预测误差的大小和模式等因素来比较这两个模型。目前,我们注意到包含相对湿度的模型表现最佳。
摘要
在本章中,我们复制了 Barkjohn 的分析。我们创建了一个模型,校正 PurpleAir 测量,使其与 AQS 测量接近。这个模型的准确性使得 PurpleAir 传感器可以包含在官方的美国政府地图上,比如 AirNow Fire and Smoke 地图。重要的是,这个模型提供了及时的和准确的空气质量测量。
我们看到众包开放数据如何通过来自精确、严格维护、政府监控设备的数据进行改进。在这个过程中,我们专注于清理和合并来自多个来源的数据,但我们也适合模型以调整和改进空气质量测量。
对于这个案例研究,我们应用了本书的这一部分涵盖的许多概念。正如您所见,整理文件和数据表以便分析是数据科学中的一个重要部分。我们使用文件整理和来自第八章的粒度概念来准备两个来源以进行合并。我们将它们组织成结构,以便匹配相邻的空气质量传感器。数据科学中的这个“肮脏”部分对于通过增加众包开放数据来扩展严格维护的精确政府监控设备的数据的覆盖范围至关重要。
这个准备过程涉及了对数据的深入、仔细的检查、清理和改进,以确保它们在两个来源之间的兼容性和在我们的分析中的可信度。来自第九章的概念帮助我们有效地处理时间数据,并找到和修正了许多问题,如缺失数据点甚至是重复的数据值。
文件和数据整理、探索性数据分析以及可视化是许多分析工作的重要组成部分。虽然拟合模型可能看起来是数据科学中最激动人心的部分,但了解和信任数据是至关重要的,通常会在建模阶段带来重要的见解。与建模相关的主题占据了本书其余大部分内容。然而,在我们开始之前,我们会涵盖与数据整理相关的另外两个主题。在下一章中,我们将展示如何从文本创建可分析的数据,而在接下来的一章中,我们将研究我们在第八章中提到的其他源文件格式。
在你进入下一章之前,回顾一下你到目前为止学到的内容。对自己的成就感到自豪——你已经走了很长的一段路!我们在这里涵盖的原则和技术对几乎每一种数据分析都很有用,你可以立即开始将它们应用到你自己的分析中。
^(1) 这种估算是基于假设地球是完美的球体。然后,一度纬度是地球的半径(以米为单位)。插入地球的平均半径后,得到每度纬度 111,111 米。经度也是一样的,但是每个“环”围绕地球的半径随着接近极点而减小,所以我们通过一个cos ( lat ) 因子进行调整。事实证明地球并不是完美的球体,所以这些估算不能用于像着陆火箭这样精确计算的目的。但对于我们的目的来说,它们表现得相当好。