统计思维：程序员数学之概率统计（第2版）：第1章 探索性数据分析

第1章 探索性数据分析

如果能将数据与实际方法相结合，就可以在存在不确定性时解答问题并指导决策，这就是本书的主题。

举个例子。我的妻子在怀第一胎时，我听到了一个问题：第一胎是不是经常晚于预产期出生？下面所给出的案例研究就是由这个问题引出的。

如果用谷歌搜索这个问题，会看到大量的讨论。有人认为第一胎的生产日期确实经常晚于预产期，有人认为这是无稽之谈，还有人认为恰恰相反，第一胎常常会早产。

在很多此类讨论中，人们会提供数据来支持自己的观点。我发现很多论据是下面这样的。

“我有两个朋友最近都刚生了第一个孩子，她们都是超过预产期差不多两周才出现临产征兆或进行催产的。”

“我的第一个孩子是过了预产期两周才出生的，我觉得第二个孩子可能会早产两周！”

“我认为这种说法不对，因为我姐姐是头生子，而且是早产儿。我还有好些表兄妹也是这样。”

这些说法都是基于未公开的数据，通常来自个人经验，因此称为轶事证据（anecdotal evidence）。在闲聊时讲讲轶事当然无可厚非，所以我并不是要批评以上那几个人。

但是，我们可能需要更具说服力的证据以及更可靠的回答。如果按照这个标准进行衡量，轶事证据通常都靠不住，原因有如下几点。

• 观测值数量较小
  如果第一胎的孕期的确偏长，这个时间差与正常的偏差相比可能很小。在这种情况下，我们可能需要比对大量的孕期数据，才能确定这种时间差确实存在。
• 选择数据时存在偏倚
  人们之所以参与这个问题的讨论，有可能是因为自己的第一个孩子出生较晚。这样的话，这个选择数据的过程就会对结果产生影响。
• 确认数据时存在偏倚
  赞同这种说法的人也许更可能提供例子进行佐证。持怀疑态度的人则更可能引用反例。
• 不精确
  轶事通常都是个人经验，经常会记错、误传或者误解等。

那我们该如何更好地回答这个问题呢？

1.1 统计学方法

为了解决轶事证据的局限性，我们将使用以下统计学工具。

• 数据收集
  我们将使用大型的全国性调查数据，这个调查专门设计用于对美国人口进行有效的统计推断。
• 描述性统计
  得出统计量，对数据进行简要的汇总，并评估可视化数据的不同方法。
• 探索性数据分析
  寻找各种模式、差异，以及其他能够解决我们感兴趣的问题的特征，同时还将检查数据的不一致性，发现局限性。
• 估计
  使用样本数据来估计一般总体的统计特征。
• 假设检验
  如果看到明显的效应，例如两个群组之间存在差异，将衡量该效应是否是偶然产生的。

谨慎执行上面的步骤，并避免各种错误，我们就可以获得合理性和准确性更高的结论。

1.2 全国家庭增长调查

从1973年起，美国疾病控制和预防中心（CDC）就开始进行全国家庭增长调查（NSFG，http://cdc.ogv/nchs/nsfg.htm），以收集“与家庭生活、婚姻状况、妊娠情况、生育情况、避孕情况，以及两性健康相关的信息。此项调查的结果用于……进行健康服务和健康教育项目的规划，以及对家庭、生育及健康情况进行统计研究”。

我们将使用这项调查收集到的数据研究第一胎是否出生较晚，并解答一些其他问题。为了有效地使用这些数据，我们必须理解这项研究是如何设计的。

全国家庭增长调查是一项横截面（cross-sectional）研究，也就是说该研究捕获的是一个群组在某一时刻的快照。在横截面研究之外，最常见的是纵向（longitudinal）研究，指在一个时间段内重复观察一个群组。

全国家庭增长调查进行过7次，每一次都称为一个周期（cycle）。我们将使用第6次的数据，其时间段为2002年1月至2003年3月。

这项调查的目的是对一个总体（population）得出结论。全国家庭增长调查的目标总体是居住在美国、年龄在15~44岁的人。理想情况下，调查要收集这个总体中每个成员的数据，但这是不可能实现的。实际上，我们收集了这个总体的一个子集的数据，这个子集称为样本（sample）。参与调查的人称为调查参与者（respondent）。

通常来说，横截面研究应该是有代表性（representative）的，也就是说目标总体中每个成员参与调查的机会均等。这种理想条件在实践中很难实现，但是进行调查的人员会竭尽所能满足这个条件。

全国家庭增长调查不具有代表性，而是特意进行过度抽样（oversample）。这项研究的设计者招募了拉美裔美国人、非洲裔美国人和青少年3个群组的参与者，每个群组的招募比例都超过其在美国人口中所占的比例，以确保各群组的参与者数量足够多，从而进行有效的统计推断。

当然，过度抽样也有缺点，那就是不容易从调查的统计数据中得出关于总体的结论。我们稍后会对此进行讨论。

在使用这种调查数据时，我们必须熟悉代码本（codebook），这一点非常重要。代码本记录了一项研究的设计、使用的调查问题，以及调查中响应变量的编码。你可以从美国疾病控制和预防中心的网站（http://www.cdc.gov/nchs/nsfg/nsfg_cycle6.htm）下载全国家庭增长调查数据的代码本和使用手册。

1.3 数据导入

本书所用的代码和数据都可以通过GitHub（https://github.com/AllenDowney/ThinkStats2）获取。前言中介绍了如何下载和使用这些代码。

下载代码后，你会得到一个名为ThinkStats2/code的文件夹，其中包含一个名为nsfg.py的文件。运行nsfg.py会读取数据文件，执行测试，然后打印出一条消息，例如“All test passed”。

让我们看看这个文件所执行的工作。第6次全国家庭增长调查的妊娠数据保存在名为2002FemPreg.dat.gz的文件中，这是一个纯文本（ASCII码）形式的gzip压缩文件，有固定宽度的列。这个文件中的每一行都是一个记录（record），包含一次妊娠的数据。

2002FemPreg.dct是一个Stata字典文件，记录了数据文件的格式。Stata是一个统计软件。Stata“字典”是由变量名、变量类型及标识变量位置的索引值组成的列表。

下面几行摘自2002FemPreg.dct：

infile dictionary { 
  _column(1) str12 caseid %12s "RESPONDENT ID NUMBER" 
  _column(13) byte pregordr %2f "PREGNANCY ORDER (NUMBER)" 
}

这个字典描述了两个变量：caseid是一个长度为12的字符串，代表调查参与者的ID；pregorder是一个单字节整数，说明这条记录描述的是这位调查参与者的第几次妊娠。

下载的代码包含一个thinkstats2.py文件，这是一个Python模块，包含了本书中用到的很多类和函数，其中有读取Stats字典和全国家庭增长调查数据文件的函数。这两个函数在nsfg.py中的用法如下：

def ReadFemPreg(dct_file='2002FemPreg.dct', 
                dat_file='2002FemPreg.dat.gz'): 
    dct = thinkstats2.ReadStataDct(dct_file) 
    df = dct.ReadFixedWidth(dat_file, compression='gzip') 
    CleanFemPreg(df) 
    return df

ReadStataDct的参数是字典文件名，返回值dct是一个FixedWidthVariables对象，其中包含从字典文件中得到的信息。dct对象提供ReadFixdWidth方法进行数据文件的读取。

1.4 DataFrame

ReadFixedWidth方法返回一个DataFrame对象。DataFrame是pandas提供的基础数据结构。pandas是一个Python数据和统计包，它的使用会贯穿本书。在DataFrame中，每个记录为一行（在我们的例子中就是每个妊娠数据为一行），每个变量为一列。

除了数据，DataFrame还包含变量名和变量类型信息，并提供访问和修改数据的方法。

如果打印df对象，你会看到其中行列的部分数据和DataFrame的大小：13 593行/记录，244列/变量。

>>> import nsfg
>>> df = nsfg.ReadFemPreg()
>>> df
...
[13593 rows x 244 columns]

df的columns属性将列名返回为一列Unicode字符串。

>>> df.columns
Index([u'caseid', u'pregordr', u'howpreg_n', u'howpreg_p', ... ])

df.columns的结果是一个Index对象，Index也是一个pandas数据结构。我们稍后会详细介绍Index，现在可以暂时将其视为一个列表。

>>> df.columns[1]
'pregordr'

要访问DataFrame中的一列，你可以将列名作为键值。

>>> pregordr = df['pregordr']
>>> type(pregordr)
<class 'pandas.core.series.Series'>

其结果是一个Series对象，这又是一个pandas数据结构。Series与Python列表类似，还能提供一些附加功能。打印一个Series对象会得到索引和对应的数值。

>>> pregordr
0     1
1     2
2     1
3     2
...
13590    3
13591    4
13592    5
Name: pregordr, Length: 13593, dtype: int64

这个示例中的索引是从0到13 592的整数，但通常索引可以使用任何可排序的数据类型。这个示例中的元素也是整数，但元素可以是任何类型的。

示例中的最后一行列出了变量名、Series长度和数据类型。int64是NumPy提供的类型之一。如果在32位机器上运行这个示例，得到的数据类型可能是int32。

你可以使用整数的index和slice值访问Series中的元素。

>>> pregordr[0]
1
>>> pregordr[2:5]
2    1
3    2 
4    3
Name: pregordr, dtype: int64

index操作符的结果是int64，slice的结果还是一个Series。

你也可以使用点标记法来访问DataFrame中的列。

>>> pregordr = df.pregordr

只有当列名为合法的Python标识符时（即以字母开头，不包含空格等），才能使用这种写法。

1.5 变量

我们已经使用了全国家庭增长调查数据集中的两个变量——caseid和pregordr，还看到数据集中共有244个变量。本书的探索性分析用到如下变量。

• caseid：调查参与者的整数ID。
• prglength：妊娠周数，是一个整数。
• outcome：怀孕结果的整数代码。1代表成功生产。
• pregordr：妊娠的顺序号。例如，一位调查参与者的第一次妊娠为1，第二次为2，以此类推。
• birthord：成功生产的顺序号，一位调查参与者的第一个孩子代码为1，以此类推。对没有成功生产的其他妊娠结果，此字段为空。
• birthwgt_lb和birthwgt_oz：新生儿体重的磅部分数值和盎司部分数值。
• agepreg：妊娠结束时母亲的年龄。
• finalwgt：调查参与者的统计权重。这是一个浮点数，表示这位调查参与者在全美人口中代表的人数。

如果你仔细阅读了代码本，就会发现这些变量中很多都是重编码（recode），也就是说这些不是调查收集的原始数据（raw data），而是使用原始数据计算得到的。

例如，如果成功生产，prglngth的值就与原始变量wksgest（妊娠周数）相等；否则，prglngth的值估算为mosgest * 4.33（妊娠月数乘以一个月的平均周数）。

重编码通常都基于一定的逻辑，这种逻辑用于检查数据的一致性和准确性。一般情况下，如果数据中存在重编码，我们就直接使用，除非有特殊的原因需要自己处理原始数据。

1.6 数据变换

导入调查数据时，经常需要检查数据中是否存在错误，处理特殊值，将数据转换为不同的格式并进行计算。这些操作都称为数据清洗（data cleaning）。

nsfg.py包含一个CleanFemPreg函数，用于清洗计划使用的变量。

def CleanFemPreg(df): 
    df.agepreg /= 100.0

    na_vals = [97, 98, 99]
    df.birthwgt_lb.replace(na_vals, np.nan, inplace=True)
    df.birthwgt_oz.replace(na_vals, np.nan, inplace=True)

    df['totalwgt_lb'] = df.birthwgt_lb + df.birthwgt_oz / 16.0

agepreg包含母亲在妊娠结束时的年龄。在数据文件中，agepreg是以百分之一年为单位的整数值。因此CleanFemPreg的第一行将每个agepreg除以100，从而获得以年为单位的浮点数值。

birthwgt_lb和birthwgt_oz包含成功生产时的新生儿体重，分别是磅和盎司的部分。这两个变量还使用几个特殊的代码。

97 NOT ASCERTAINED 
98 REFUSED 
99 DON'T KNOW

用数字编码特殊值是一种危险的做法，因为如果没有进行正确的处理，这些数字可能产生虚假结果，例如，99磅重的新生儿。replace方法可以将这些值替换为np.nan，这是一个特殊的浮点数值，表示“不是数字”。replace方法使用inplace标识，说明直接修改现有的Series对象，而不是创建新对象。

IEEE浮点数表示法标准中规定，在任何算术运算中，如果有参数为nan，结果都返回nan。

>>> import numpy as np
>>> np.nan / 100.0
nan

因此使用nan进行计算会得到正确的结果，而且大部分的pandas函数都能恰当地处理nan。但我们经常需要处理数据缺失的问题。

CleanFemPreg函数的最后一行生成一个新列totalwgt_lb，将磅和盎司值结合在一起，得到一个以磅为单位的值。

需要注意的是，向DataFrame添加新列时，必须使用如下字典语法：

# 正确 
df['totalwgt_lb'] = df.birthwgt_lb + df.birthwgt_oz / 16.0 

而不是使用点标记：

# 错误！
df.totalwgt_lb = df.birthwgt_lb + df.birthwgt_oz / 16.0

使用点标记的写法会给DataFrame对象添加一个新属性，而不是创建一个新列。

1.7 数据验证

当数据从一个软件环境导出，再导入另一个环境时，可能会产生错误。如果不熟悉新数据集，可能会对数据进行不正确的解释，或者引入其他的误解。如果能抽出一些时间进行数据验证，就可以节省后续可能花费的时间，避免可能出现的错误。

验证数据的一种方法是计算基本的统计量，并与已发布的结果进行比较。例如，全国家庭增长调查的代码本为每个变量提供了概要表。outcome变量对每个妊娠结果进行了编码，其概要表如下：

value label       Total 
1 LIVE BIRTH          9148 
2 INDUCED ABORTION    1862 
3 STILLBIRTH           120 
4 MISCARRIAGE         1921 
5 ECTOPIC PREGNANCY    190 
6 CURRENT PREGNANCY    352

Series类提供了一个value_counts方法，可用于计算每个值出现的次数。如果得到DataFrame中的outcome Series，我们可以使用value_counts方法，将结果与已发布的数据进行比较。

>>> df.outcome.value_counts().sort_index()
1    9148
2    1862
3     120
4    1921 
5     190 
6     352

value_counts返回的结果是一个Series对象。sort_index方法将Series对象按索引排序，使结果按序显示。

我们将得到的结果与官方发布的表格进行对比，outcome变量的值似乎没有问题。类似地，已发布的关于birthwgt_lb的概要表如下：

value label                Total 
. INAPPLICABLE          4449 
0-5 UNDER 6 POUNDS        1125 
6 6 POUNDS              2223 
7 7 POUNDS              3049 
8 8 POUNDS              1889 
9-95 9 POUNDS OR MORE       799

birthwgt_lb的value_counts结果如下：

>>> df.birthwgt_lb.value_counts(sort=False) 
0        8 
1       40 
2       53 
3       98 
4      229 
5      697 
6     2223 
7     3049 
8     1889 
9      623 
10     132 
11      26 
12      10 
13       3 
14       3 
15       1 
51       1

数值6、7、8的出现次数是正确的。如果计算出0~5和9~95的次数，结果也是正确的。但是，如果再看仔细些，你会发现有一个数值肯定是错的——一个51磅的新生儿！

为了处理这个错误，可以在CleanFemPreg中加入一行代码。

df.birthwgt_lb[df.birthwgt_lb > 20] = np.nan

这行代码将非法值替换为np.nan。方括号中的表达式产生一个bool类型的Series对象，值为True表示满足该条件。当一个布尔Series用作索引时，它只选择满足该条件的元素。

1.8 解释数据

要想有效使用数据，就必须同时在两个层面上思考问题：统计学层面和上下文层面。

例如，让我们看一看几位调查参与者的outcome序列。由于数据文件的组织方式，我们必须进行一些处理才能得到每位调查参与者的妊娠数据。以下函数实现了我们需要的处理：

def MakePregMap(df): 
    d = defaultdict(list) 
    for index, caseid in df.caseid.iteritems(): 
        d[caseid].append(index) 
    return d

df是包含妊娠数据的DataFrame对象。iteritems方法遍历所有妊娠记录的索引（行号）和caseid。

d是将每个caseID映射到一列索引的字典。如果你不熟悉defaultdict，可以到Python的collections模块中查看其定义。使用d，我们可以查找一位调查参与者，获得其妊娠数据的索引。

下面的示例就查找了一位调查参与者，并打印出其妊娠结果列表：

>>> caseid = 10229
>>> indices = preg_map[caseid]
>>> df.outcome[indices].values
[4 4 4 4 4 4 1]

indices是调查参与者10229的妊娠记录索引列表。

以这个列表为索引可以访问df.outcome中指定的行，获得一个Series。上面的示例没有打印整个Series对象，而是选择输出values属性，这个属性是一个NumPy数组。

输出结果中的代码1表示成功分娩。代码4表示流产，即自发终止的妊娠，终止原因通常未知。

从统计学上看，这位调查参与者并无异常。流产并不少见，其他一些调查参与者的流产次数相同或者更多。

但是考虑到上下文，这个数据说明一位妇女怀孕6次，每次都以流产告终。她第7次也是最近一次怀孕成功产下了孩子。如果我们抱着同情心看待这些数据，就很容易被数据背后的故事感动。

全国家庭增长调查数据集中的每一条记录都代表一位参与者，这些参与者诚实地回答了很多非常私密而且难以回答的问题。我们可以使用这些数据解答与家庭生活、生育和健康相关的统计学问题。同时，我们有义务思及这些数据所代表的参与者，对他们心存敬意和感谢。

1.9 练习

• 练习1.1

  你下载的代码中应该有一个名为chap01ex.ipynb的文件，这是一个IPython记事本。你可以用如下命令从命令行启动IPython记事本：

  $ ipython notebook &

  如果系统安装了IPython，会启动一个在后台运行的服务器，并打开一个浏览器查看记事本。如果你不熟悉IPython，我建议你从IPython网站（http://ipython.org/ipython-doc/stable/notebook/notebook.html）开始学习。

  你可以添加一个命令行选项，使图片在“行内”（即在记事本中）显示，而非弹出窗口：

  $ ipython notebook --pylab=inline &

  打开chap01ex.ipynb。记事本中一些单元已经填好了代码，可以直接执行。其他单元列出了你应该尝试的练习。

  本练习的参考答案在chap01soln.ipynb中。

• 练习1.2

  创建一个名为chp01ex.py的文件，编写代码，读取参与者文件2001FemResp.dat.gz。你可以复制nsfg.py文件并对其进行修改。

  变量pregnum是一个重编码，用于说明每位调查参与者有过多少次妊娠经历。打印这个变量中不同值的出现次数，将结果与全国家庭增长调查代码本中发布的结果进行比较。

  你也可以将每位调查参与者的pregnum值与妊娠文件中的记录数进行比较，对调查参与者文件和妊娠文件进行交叉验证。

  你可以使用nsfg.MakePregMap生成一个字典，将每个caseid映射到妊娠DataFrame的索引列表。

  本练习的参考答案在chp01soln.py中。

• 练习1.3

  学习统计学的最好方法是使用一个你感兴趣的项目。你想研究“第一胎是否都会晚出生”这样的问题吗？

  请思考一些你个人感兴趣的问题，可以是传统观点、争议话题或影响政局的问题，看是否可以构想出一个能以统计调查进行验证的问题。

  寻找能帮助你回答这个问题的数据。公共研究的数据经常可以免费获取，因此政府网站是很好的数据来源，如http://www.data.gov/和http://www.science.gov/。如果想获得英国的数据，可以访问http://data.gov.uk/。

  我个人最喜爱的两个数据集是General Social Survey（http://www3.norc.org/gss+website/）和European Social Survey（http://www.europeansocialsurvey.org）。

  如果有人看似已经解答了你的问题，那么仔细检查该回答是否合理。数据和分析中可能存在的缺陷都会使结论不可靠。如果发现别人的解答存在问题，你可以对同样的数据进行不同的分析，或者寻找更好的数据来源。

  如果有一篇论文解答了你的问题，那么你应该能够获得论文使用的原始数据。很多论文作者会把数据放在网上供大家使用，但如果涉及敏感信息，你可能需要向作者写信索要，提供你计划如何使用这些数据的信息，或者同意某些使用条款。坚持就是胜利！

1.10 术语

• 轶事证据（anecdotal evidence）

  随意收集，而非通过精心设计的研究获得的证据，通常是个人证据。

• 总体（population）

  在研究中，我们感兴趣的群组。“总体”经常指一组人，但这个词也可以用于其他对象。

• 横截面研究（cross-sectional study）

  收集一个总体在某个特定时间点的数据的研究。

• 周期（cycle）

  在重复进行的横截面研究中，每次研究称为一个周期。

• 纵向研究（longtitudinal study）

  在一段时间内跟踪一个总体的研究，从同一个群体重复收集数据。

• 记录（record）

  在数据集中，关于单个人或其他对象的信息集合。

• 调查参与者（respondent）

  参与调查的人。

• 样本（sample）

  总体中用于数据收集的一个子集。

• 有代表性（representative）

  如果总体中的每个成员被选入样本的机会都均等，那么这个样本就是有代表性的。

• 过度抽样（oversampling）

  一种通过增加一个子总体的样本数来避免因样本规模过小产生错误的技术。

• 原始数据（raw data）

  没有经过或只经过少许检查、计算或解释，直接收集和记录的值。

• 重编码（recode）

  通过计算和应用于原始数据的其他逻辑生成的值。

• 数据清洗（data cleaning）

  数据处理过程，包括数据验证、错误检查，以及数据类型和表示的转换等。

from: http://www.ituring.com.cn/tupubarticle/3914