第3章 统计实验与显著性检验

《Practical Statistics for Data Scientists》书籍英文版
《面向数据科学家的实用统计学》中文版书籍

第 3 章 统计实验与显著性检验

实验设计是统计实践的基石,在几乎所有研究领域都有应用。其目标是设计实验以确认或拒绝某个假设。数据科学家往往需要持续进行实验,尤其是关于用户界面和产品营销方面的实验。本章回顾了传统实验设计,并讨论了数据科学中常见的一些挑战;还介绍了一些统计推断中经常被引用的概念,并解释了它们的含义及其与数据科学的相关性(或不相关性)。

F3.1

当你看到统计显著性、t 检验或 p 值等术语时,通常是在经典统计推断“流水线”的上下文中(见图 3-1)。这个过程从一个假设开始(例如“药物 A 优于现有标准药物”或“价格 A 比现有价格 B 更有利可图”)。然后设计实验(可能是 A/B 测试)以检验这个假设——设计得尽可能能够得出结论性结果。接着收集并分析数据,然后得出结论。术语“推断”体现了这样一种意图:将涉及有限数据集的实验结果,应用到更大的过程或总体上。

A/B 测试

A/B 测试是一种包含两个组的实验,用于确定两种处理、产品、程序等哪一种更优。通常两种处理中的一种是现有的标准处理,或是不做处理。如果使用标准处理(或不处理),则称其为对照组(control)。典型的假设是新处理优于对照组。

A/B 测试关键术语

  • 处理(Treatment):让受试对象接触的某种东西(药物、价格、网页标题等)。
  • 处理组(Treatment group):接触特定处理的一组受试对象。
  • 对照组(Control group):接触标准处理或不处理的一组受试对象。
  • 随机化(Randomization):将受试对象随机分配到不同处理的过程。
  • 受试对象(Subjects):暴露在处理下的个体(网页访问者、患者等)。
  • 检验统计量(Test statistic):用于衡量处理效果的指标。

A/B 测试在网页设计和营销中很常见,因为结果容易衡量。A/B 测试的一些例子包括:

  • 测试两种土壤处理方法,看哪种能让种子发芽更好
  • 测试两种疗法,看哪种对抑制癌症更有效
  • 测试两种定价,看哪种能带来更多净利润
  • 测试两个网页标题,看哪一个能带来更多点击(图 3-2)
  • 测试两个网页广告,看哪一个能带来更多转化

F3.2

一个合格的 A/B 测试要求受试对象可以被分配到某个处理或另一个处理。受试对象可以是人、植物种子、网页访问者;关键是受试对象接触到处理。理想情况下,受试对象要随机分配到处理组。这样你就知道,处理组之间的差异只可能来自两种原因:

  • 不同处理的效果
  • 随机分配带来的偶然性(比如天生表现更好的受试对象恰好集中在 A 组或 B 组)

你还需要注意用于比较 A 组和 B 组的检验统计量或指标。数据科学中最常见的指标是二元变量:点击或不点击、购买或不购买、欺诈或不欺诈等等。这些结果可以汇总在一个 2×2 表中。表 3-1 就是一个实际价格测试的 2×2 表(关于这些结果的进一步讨论见第 103 页“统计显著性与 p 值”)。

T3.1

如果指标是连续变量(购买金额、利润等)或计数(例如住院天数、访问页面数),结果可能以不同方式展示。如果关注的不是转化率而是每次页面浏览的收入,那么表 3-1 中的价格测试结果在典型软件输出中可能看起来是这样的:

  • 价格 A 每次页面浏览收入:平均值 = 3.87,标准差 SD = 51.10
  • 价格 B 每次页面浏览收入:平均值 = 4.11,标准差 SD = 62.98

“SD”指的是各组内值的标准差。

警告:

但要注意:统计软件(包括 R 和 Python)默认生成的输出,并不意味着所有输出都有用或相关。你可以看到,上述标准差并不太有用;表面上它们暗示很多值可能是负的,而负收入在现实中并不可行。这类数据由少量较高值(带转化的页面浏览)和大量 0 值(无转化的页面浏览)组成。用一个数字来概括这种数据的变异性是很困难的,虽然平均绝对偏差**(A 为 7.68,B 为 8.15)比标准差更合理一些。

为什么要有对照组?

为什么不跳过对照组,只对一个组施加感兴趣的处理,然后把结果和以往经验相比呢?

没有对照组,就不能保证“其他条件都相同”,也不能保证任何差异确实是由处理(或随机因素)引起的。当你有一个对照组时,除了研究的处理之外,它与处理组处在同样的条件下。如果只是与“基线”或以往经验比较,处理之外的其他因素也可能不同。

通用注解:

研究中的盲法**:盲法(blind study)是指受试者不知道自己接受的是A处理还是B处理。知道自己接受某种处理可能会影响反应。双盲研究(double-blind study)是指研究者和执行者(如医疗研究中的医生和护士)也不知道哪些受试者接受了哪种处理。当处理本身显而易见时(例如计算机认知疗法与心理学家面对面的认知疗法),就无法进行盲法。

A/B测试在数据科学中通常用于网页场景。处理可能是网页的设计、商品价格、标题的措辞或其他项目。需要一定的思考来保持随机化原则。通常实验中的“受试者”是网站访客,而我们感兴趣的结果是点击、购买、访问时长、访问页面数量、是否访问某一特定页面等在标准A/B实验中,你需要事先决定一个度量指标。可能会收集多个行为指标且都感兴趣,但如果实验预期要在处理A和处理B之间做出决策,那么必须在实验前确定一个单一指标或检验统计量。在实验结束后才选择检验统计量会引入研究者偏倚。

为什么只是A/B?为什么不C、D……?

A/B测试在市场营销和电商领域很流行,但并不是唯一的统计实验类型。可以纳入更多处理。受试者可能会被多次测量。在药物试验中,受试者稀缺、昂贵且需要长期招募,有时会设计多个中途停止实验并得出结论的机会。

传统的统计实验设计侧重于回答一个静态问题

“价格A和价格B之间的差异在统计上显著吗?”

而数据科学家更关心的问题是:

“在多个可能的价格中,哪个最好?”

为此,使用了一种相对较新的实验设计:多臂老虎机算法(multi-arm bandit,见本书第131页)。

警告:

获得许可:**在涉及人类受试者的科学和医学研究中,通常需要获得受试者的许可,并得到机构审查委员会(IRB)的批准。而在商业中作为日常运营一部分进行的实验几乎从不这么做。在大多数情况下(如定价实验、显示哪个标题或提供哪个优惠),这种做法被普遍接受。然而,Facebook在2014年因一次实验引发广泛批评。他们操纵用户新闻推送(newsfeed)中的情绪色彩:用情感分析将帖子分类为正面或负面,然后改变展示给用户的正/负比例。一些随机选择的用户看到更多正面帖子,另一些用户看到更多负面帖子。Facebook发现,看到更正面新闻推送的用户更可能自己发布正面内容,反之亦然。这种效应的幅度虽然很小,但Facebook因为在不告知用户的情况下进行实验而受到大量批评。一些用户甚至推测,若极度抑郁的用户看到负面版本的推送,可能会被“推向崩溃的边缘”。

关键要点

  • 受试者被分配到两个(或更多)组,这些组除了研究的处理不同之外,其他条件完全相同。
  • 理想情况下,受试者被随机分配到各组。

假设检验

假设检验(hypothesis tests),也叫显著性检验(significance tests),在传统的已发表研究的统计分析中无处不在。它们的目的是帮助你判断,随机因素是否可能是观测到的效果的原因。

假设检验的关键术语

  • 原假设(Null hypothesis) 认为随机因素是原因的假设。

  • 备择假设(Alternative hypothesis) 原假设的对立面(你希望证明的东西)。

  • 单尾检验(One-way test) 只在一个方向上计算偶然结果的假设检验。

  • 双尾检验(Two-way test) 在两个方向上都计算偶然结果的假设检验。

A/B 测试(见第88页“A/B Testing”)通常是带着假设来设计的。例如,假设可能是“价格B带来更高利润”。为什么我们需要一个假设?为什么不直接看实验的结果,选那个更好的处理就行?答案在于人类大脑倾向于低估自然随机行为的幅度。这种倾向的一种表现是无法预见极端事件,即所谓的“黑天鹅”(见第73页“长尾分布”)。另一种表现是倾向于将随机事件误解为具有某种有意义的模式。统计假设检验就是为了防止研究人员被随机性所欺骗而发明的。

误解随机性

你可以通过这个小实验观察人类低估随机性的倾向:请几个朋友“编造”一系列 50 次掷硬币的结果,让他们写下一串随机的“H”和“T”。然后再让他们真的掷硬币 50 次,并把结果写下来。让他们把真正的掷硬币结果放在一堆,编造的结果放在另一堆。很容易看出哪些是真实结果:真实结果中“H”或“T”的连续串更长。在 50 次真实掷硬币中,连续出现 5~6 次“H”或“T”一点都不罕见。然而,当我们自己在编造随机掷硬币结果时,如果已经连续 3~4 次“H”,就会告诉自己“为了看起来随机”,最好换成“T”。从另一方面看,当我们在现实世界里看到“连续 6 次H”的等价情况(比如一个标题的点击率比另一个高出10%),我们往往倾向于把它归因于某种真实原因,而不仅仅是偶然。

在一个设计合理的 A/B 测试中,你以这样的方式收集A组和B组的数据,使得任何观测到的差异只能是以下两种原因之一:

  • 分配受试者时的随机因素
  • A和B之间的真实差异

统计假设检验就是对A/B测试或任何随机实验的进一步分析,以评估“随机因素”是否是解释A组和B组差异的合理原因。

原假设

The Null Hypothesis

假设检验的逻辑是这样的:“鉴于人类倾向于对不寻常但随机的行为作出反应,并把它解读为有意义、真实的东西,在我们的实验中,我们将要求证明,组间的差异比偶然因素合理产生的差异更为极端。”这包含一个基线假设:不同处理是等效的,组间的任何差异都源于随机因素。这个基线假设称为原假设(null hypothesis)。我们的希望是能够推翻原假设,证明 A 组与 B 组的结果比随机机会所能产生的差异更大。

一种方法是使用重抽样置换程序:把 A 组和 B 组的结果混合打乱,然后反复按相同的样本量重新分组,再观察得到与观测差异一样极端的情况出现的频率。A 组和 B 组混合后并抽样的结果体现了 A 组与 B 组等价且可互换的原假设,被称为原模型(null model)。更多细节参见第 96 页“重抽样”。

备择假设

Alternative Hypothesis

假设检验本质上不仅涉及原假设,还涉及与之相对的备择假设。例子如下:

  • 原假设:“A 组与 B 组均值无差异”;备择假设:“A 与 B 不同”(可能更大或更小)
  • 原假设:“A ≤ B”;备择假设:“A > B”
  • 原假设:“B 不比 A 高 X%”;备择假设:“B 比 A 高 X%”

原假设和备择假设加在一起必须涵盖所有可能性。原假设的性质决定了假设检验的结构。

单尾检验与双尾检验

One-Way Versus Two-Way Hypothesis Tests

在 A/B 测试中,常见的情形是:你在测试一个新选项(如 B)与一个既定默认选项(A),并假定除非新选项能明确更好,否则你会继续沿用默认选项。在这种情况下,你希望假设检验保护你免于在有利于 B 的方向上被随机性误导。至于相反方向,你并不在意,因为除非 B 被证明确实更好,否则你会坚持 A。因此你需要一个有方向性的备择假设(B 优于 A)。这种情况下使用单尾(单向)假设检验,即只有一个方向上的极端随机结果计入 p 值。

如果你希望假设检验能保护你免于在任意方向上被随机性误导,那么备择假设是双向的(A 与 B 不同,可能更大或更小)。此时应使用双尾(双向)假设检验,即两个方向上的极端随机结果都计入 p 值。

单尾假设检验常常符合 A/B 决策的实际性质:通常一方被赋予“默认”地位,只有另一方更好时才替换。但软件(包括 R 和 Python 的 scipy)通常默认输出双尾检验结果,许多统计学家也倾向于使用更保守的双尾检验以避免争论。单尾与双尾的区别是一个容易混淆的话题,在数据科学中并不太重要,因为 p 值的精确性并不是特别关键。

关键要点

  • 原假设是一种逻辑构造,体现“什么特殊情况都没发生,观测到的效应完全是随机机会造成的”这一概念。
  • 假设检验假定原假设成立,构建“原模型”(一个概率模型),并检验观测到的效应是否是该模型下的合理结果。

重抽样

Resampling

在统计学中,重抽样指的是从已观测到的数据中反复抽取数值,其总体目标是评估某个统计量的随机变动性。它还可以用来评估和提高某些机器学习模型的准确度(例如,把多个自助抽样数据集上建立的决策树模型的预测结果取平均,这一过程称为“袋装法”——参见第 259 页“袋装法与随机森林”)。

重抽样程序主要有两种:自助法(bootstrap)置换检验(permutation test)。自助法用来评估估计值的可靠性,在上一章已讨论过(见第 61 页“自助法”)。置换检验用于检验假设,通常涉及两个或多个组,本节将讨论这种方法。

重抽样的关键术语

置换检验(Permutation test) 把两个或多个样本合并在一起,并随机(或穷尽地)重新分配观测值到新的重抽样中。

同义词 随机化检验、随机置换检验、精确检验

重抽样(Resampling) 从一个已观测到的数据集抽取附加样本(“重抽样”)。

有放回或无放回(With or without replacement) 抽样时,每次抽出一个项目后,是否在下一次抽样前将其放回样本中。

置换检验

Permutation Test

在置换程序中,涉及两个或更多的样本,通常是 A/B 测试或其他假设检验中的各组。Permute 的意思是改变一组数值的顺序。置换检验的第一步是将 A 组和 B 组(如果有 C、D……则包括它们)的结果合并。这体现了原假设的逻辑,即各组接受的处理并无差异。随后,我们通过从这个合并集里随机抽取各组,来检验这一假设并观察它们之间的差异。

置换程序如下:

  1. 将不同组的结果合并为一个数据集。
  2. 将合并后的数据打乱,然后随机抽取(无放回)一个与 A 组同样大小的重抽样(显然其中会包含来自其他组的数据)。
  3. 从剩余数据中随机抽取(无放回)一个与 B 组同样大小的重抽样。
  4. 对 C、D 组等做同样操作。至此,你得到了一组与原始样本大小对应的重抽样。
  5. 对原始样本计算的统计量或估计值(如组间比例差异),在重抽样上重新计算并记录;这构成了一次置换迭代。
  6. 重复上述步骤 R 次,以得到检验统计量的置换分布。

然后回到实际观测到的组间差异,把它与置换得到的差异集合进行比较。如果观测到的差异完全落在置换差异集合之内,我们并不能证明任何东西——观测差异在随机机会可能产生的范围内。但如果观测差异大多落在置换分布之外,则可以认为并非随机机会导致。从技术角度说,这个差异具有统计显著性(参见第 103 页“统计显著性与 p 值”)。

示例:网站黏性

Web Stickiness

一家销售相对高价值服务的公司想要测试两种网页展示哪一种更能促进销售。由于所售服务价值高,销售频率低且周期长,如果直接看销售量,要花很长时间才能积累足够数据判断哪种展示更优。因此,公司决定用一个代理变量来衡量结果,具体是使用描述该服务的详细内部页面。

知识点:

代理变量**(proxy variable)是用来替代真正感兴趣的变量的指标,这个真正的变量可能不可获得、成本太高或耗时太长。例如,在气候研究中,古代冰芯中的氧含量被用作温度的代理指标。最好至少有一些真正感兴趣变量的数据,以便评估它与代理变量之间的关联强度。

对这家公司来说,一个潜在的代理变量是用户点击详细落地页的次数;一个更好的代理变量是用户在页面上停留的时间。合理的假设是,一个能更长时间吸引用户注意的网页展示更有可能带来更多销售。因此,我们的指标是平均会话时长,比较 A 页和 B 页。

由于这是一个内部的、特殊用途的页面,访客数并不多。还要注意,我们用 Google Analytics 来测量会话时长,但它无法测量用户访问的最后一个会话时长。Google Analytics 并不会将这类会话从数据中删除,而是记录为 0,因此数据需要额外处理以剔除这些会话。处理后,两种不同展示共 36 个会话,其中 A 页 21 个,B 页 15 个。

F3.3

我们可以用 ggplot 通过并排箱线图直观比较会话时长:

1
2
ggplot(session_times, aes(x=Page, y=Time)) +
geom_boxplot()

pandas 的 boxplot 命令使用 by 参数来生成图形:

1
2
3
4
ax = session_times.boxplot(by='Page', column='Time')
ax.set_xlabel('')
ax.set_ylabel('Time (in seconds)')
plt.suptitle('')

图 3-3 的箱线图表明,B 页的会话时长比 A 页更长。可以在 R 中这样计算各组均值:

1
2
3
4
mean_a <- mean(session_times[session_times['Page'] == 'Page A', 'Time'])
mean_b <- mean(session_times[session_times['Page'] == 'Page B', 'Time'])
mean_b - mean_a
[1] 35.66667

在 Python 中,我们先按页面筛选 pandas 数据框,再取 Time 列均值:

1
2
3
mean_a = session_times[session_times.Page == 'Page A'].Time.mean()
mean_b = session_times[session_times.Page == 'Page B'].Time.mean()
mean_b - mean_a

结果表明,B 页的平均会话时长比 A 页长 35.67 秒。问题在于,这个差异是否在随机因素可能产生的范围内,即是否具有统计显著性。一种方法是应用置换检验——把所有会话时长合并,然后反复打乱并按 21 个(A 页,nA=21)和 15 个(B 页,nB=15)分组。

为了应用置换检验,我们需要一个函数,把 36 个会话时长随机分配为一组 21 个(A 页)和一组 15 个(B 页)。R 版本的函数是:

1
2
3
4
5
6
7
8
perm_fun <- function(x, nA, nB)
{
n <- nA + nB
idx_b <- sample(1:n, nB)
idx_a <- setdiff(1:n, idx_b)
mean_diff <- mean(x[idx_b]) - mean(x[idx_a])
return(mean_diff)
}

Python 版本的置换检验函数如下:

1
2
3
4
5
def perm_fun(x, nA, nB):
n = nA + nB
idx_B = set(random.sample(range(n), nB))
idx_A = set(range(n)) - idx_B
return x.loc[idx_B].mean() - x.loc[idx_A].mean()

这个函数的工作方式是:不放回地抽取 \(n_B\) 个索引分配给 B 组;剩下的 \(n_A\) 个索引分配给 A 组;返回两组均值之差。 将这个函数调用 R=1,000 次,并指定 \(n_A\)=21、\(n_B\)=15,就可以得到会话时长差异的分布,并画出直方图。 在 R 中用 hist 函数实现如下:

1
2
3
4
5
6
perm_diffs <- rep(0, 1000)
for (i in 1:1000) {
perm_diffs[i] = perm_fun(session_times[, 'Time'], 21, 15)
}
hist(perm_diffs, xlab='Session time differences (in seconds)')
abline(v=mean_b - mean_a)

在 Python 中,我们可以用 matplotlib 画类似的图:

1
2
3
4
5
6
7
perm_diffs = [perm_fun(session_times.Time, nA, nB) for _ in range(1000)]
fig, ax = plt.subplots(figsize=(5, 5))
ax.hist(perm_diffs, bins=11, rwidth=0.9)
ax.axvline(x = mean_b - mean_a, color='black', lw=2)
ax.text(50, 190, 'Observed\ndifference', bbox={'facecolor':'white'})
ax.set_xlabel('Session time differences (in seconds)')
ax.set_ylabel('Frequency')

图 3-4 的直方图显示,随机置换的均值差异经常超过观测到的会话时长差异(竖线所示)。

F3.4

在我们的结果中,这种情况发生在 12.6% 的案例里:

1
2
3
mean(perm_diffs > (mean_b - mean_a))
---
0.126

由于模拟使用了随机数,这个百分比会有所不同。例如,在 Python 版本中,我们得到的是 12.1%:

1
2
3
np.mean(perm_diffs > mean_b - mean_a)
---
0.121

这表明,A 页与 B 页之间观测到的会话时长差异完全在随机波动范围内,因此没有统计显著性

穷举置换检验与自助置换检验

Exhaustive and Bootstrap Permutation Tests

除了前面这种随机洗牌程序(又叫随机置换检验或随机化检验)外,置换检验还有两个变体:

  • 穷举置换检验(exhaustive permutation test)
  • 自助置换检验(bootstrap permutation test)

在穷举置换检验中,我们不是只随机打乱并分组,而是列举出所有可能的分组方式。这只在样本量较小时可行。随机置换检验在重复次数足够多时,其结果会逼近穷举置换检验,并在极限情况下收敛于它。 由于具有保证“原模型”不会以超过检验 α 水平被判定为“显著”的统计性质,穷举置换检验有时也被称为精确检验(exact test)(参见第 103 页“统计显著性与 p 值”)。

在自助置换检验中,随机置换检验步骤 2 和步骤 3 中的抽样改为有放回而非无放回。这样,重抽样程序不仅模拟处理分配给受试者的随机因素,还模拟从总体抽取受试者的随机因素。

这两种程序在统计学中都能遇到,它们之间的区别有些复杂,但在数据科学实践中并不重要。

置换检验:数据科学中的核心要点

Permutation Tests: The Bottom Line for Data Science

置换检验是一种用来探索随机变异作用的有用启发式程序。它们相对容易编写代码、解释和说明,并且提供了一种有别于基于公式统计方法的“形式主义”和“虚假确定性”的有用路径,因为公式给出的“精确”答案往往暗示了不必要的确定性。与公式方法相比,重抽样的一个优点是,它更接近于一种“通用”推断方法。数据可以是数值型的或二值型的;样本容量可以相同也可以不同;不需要假设数据服从正态分布。

关键要点

  • 在置换检验中,先将多个样本合并再打乱顺序。
  • 然后将打乱后的值分成重抽样样本,计算感兴趣的统计量。
  • 重复这一过程,并记录重抽样统计量。
  • 将观察到的统计量与重抽样分布进行比较,可以判断样本间观察到的差异是否可能由随机性产生。

统计显著性与 p 值

Statistical Significance and p-Values

统计显著性是统计学家衡量实验(或对现有数据的研究)所得结果是否比随机性产生的结果更极端的一种方法。如果结果超出了随机变异的范围,就称之为具有统计显著性。

统计显著性与 p 值的关键术语

  • p 值:在一个体现零假设的随机模型下,得到与观测结果一样不寻常或极端的结果的概率。
  • α(显著性水平):判断“异常”程度的概率阈值,实际结果超过该阈值就被视为具有统计显著性。
  • 第一类错误(Type 1 error):错误地认为效应存在(其实是偶然造成的)。
  • 第二类错误(Type 2 error):错误地认为效应是偶然的(其实是真实存在的)。

请看前面网页测试示例的表 3-2:

T3.2

价格 A 的转化率比价格 B 高将近 5%(0.8425% = 200/(23,539+200)*100,而 0.8057% = 182/(22,406+182)*100,差值为 0.0368 个百分点),在高流量业务中已足够有意义。这里我们有 45,000 多个数据点,容易认为这是“大数据”,不需要做统计显著性检验(主要用于处理小样本中的抽样变异)。然而,转化率不到 1%,实际有意义的数值(转化次数)只有几百个,真正决定所需样本量的是这些转化次数。我们可以用重抽样程序来检验价格 A 和 B 之间的转化率差异是否在随机变异的范围内。这里“随机变异”指的是一个体现零假设(即两种价格的转化率没有差别)的概率模型所产生的随机波动。

下述置换程序提出的问题是:“如果两种价格具有相同的转化率,随机波动能否产生像 5% 这样大的差异?”

  1. 把标有 1 和 0 的卡片放入盒中:这代表假定的共享转化率,382 个 1 和 45,945 个 0 = 0.008246 = 0.8246%。
  2. 洗牌并抽取一个大小为 23,739(价格 A 的样本量)的重抽样,记录 1 的数量。
  3. 记录剩下的 22,588(价格 B 的样本量)中 1 的数量。
  4. 记录 1 的比例差异。
  5. 重复步骤 2–4。
  6. 统计差异 ≥ 0.0368 的次数。

在 R 中可以重用第 98 页“示例:网页黏性”里定义的 perm_fun 函数来创建转化率随机置换差异的直方图:

1
2
3
4
5
6
7
8
obs_pct_diff <- 100 * (200 / 23739 - 182 / 22588)
conversion <- c(rep(0, 45945), rep(1, 382))
perm_diffs <- rep(0, 1000)
for (i in 1:1000) {
perm_diffs[i] = 100 * perm_fun(conversion, 23739, 22588)
}
hist(perm_diffs, xlab='Conversion rate (percent)', main='')
abline(v=obs_pct_diff)

对应的 Python 代码为:

1
2
3
4
5
6
7
8
9
10
11
12
13
obs_pct_diff = 100 * (200 / 23739 - 182 / 22588)
print(f'Observed difference: {obs_pct_diff:.4f}%')
conversion = [0] * 45945
conversion.extend([1] * 382)
conversion = pd.Series(conversion)
perm_diffs = [100 * perm_fun(conversion, 23739, 22588)
for _ in range(1000)]
fig, ax = plt.subplots(figsize=(5, 5))
ax.hist(perm_diffs, bins=11, rwidth=0.9)
ax.axvline(x=obs_pct_diff, color='black', lw=2)
ax.text(0.06, 200, 'Observed\ndifference', bbox={'facecolor':'white'})
ax.set_xlabel('Conversion rate (percent)')
ax.set_ylabel('Frequency')

在图 3-5 中展示了 1,000 个重抽样结果的直方图:就本例而言,观察到的 0.0368% 的差异完全处于随机波动范围内。

F3.5

p 值

单纯看图并不是衡量统计显著性的精确方法,因此更有意义的是 p 值。p 值指的是“机会模型”(即零假设模型)产生比观测结果更极端结果的频率。我们可以用置换检验估计 p 值:计算置换检验产生的差异大于或等于观测差异的比例:

1
2
mean(perm_diffs > obs_pct_diff)
[1] 0.308
1
np.mean([diff > obs_pct_diff for diff in perm_diffs])

这里,R 和 Python 都利用了“true 被当作 1、false 被当作 0”这一事实。p 值为 0.308,意味着在随机条件下,我们大约 30% 的时间会得到像这样极端或更极端的结果。

在这种情况下,其实不必用置换检验来得到 p 值。由于我们有二项分布,可以近似计算 p 值。R 代码如下:

1
prop.test(x=c(200, 182), n=c(23739, 22588), alternative='greater')

输出:

1
2
3
4
5
6
7
8
9
2-sample test for equality of proportions with continuity correction
data: c(200, 182) out of c(23739, 22588)
X-squared = 0.14893, df = 1, p-value = 0.3498
alternative hypothesis: greater
95 percent confidence interval:
-0.001057439 1.000000000
sample estimates:
prop 1 prop 2
0.008424955 0.008057376

参数 x 是每组成功次数,参数 n 是每组试验次数。

scipy.stats.chi2_contingency 方法可接受表 3-2 中的数据:

1
2
3
survivors = np.array([[200, 23739 - 200], [182, 22588 - 182]])
chi2, p_value, df, _ = stats.chi2_contingency(survivors)
print(f'p-value for single sided test: {p_value / 2:.4f}')

正态近似得出 p 值为 0.3498,与置换检验的 p 值接近。

α(显著性水平)

统计学家不赞成把“结果是否太不寻常”交给研究者自由裁量,而是提前规定一个阈值,例如“比零假设下 5% 的结果更极端”。这个阈值称为 α。典型的 α 水平有 5% 和 1%。无论选择哪一个水平,本质上都是人为的决策,并不能保证 x% 的时候都作出正确判断。这是因为所回答的概率问题并不是“这个结果发生的概率是多少?”,而是“在机会模型成立的前提下,得到如此极端结果的概率是多少?”。然后我们再反推机会模型是否合适,但这种判断并不具有概率意义。这一点长期以来造成了许多混淆。

p 值争议:p-value controversy

近年来,p 值的使用引发了大量争议。某心理学期刊甚至“禁止”投稿论文使用 p 值,理由是单纯根据 p 值决定是否发表,导致了劣质研究的出版。太多研究者对 p 值的真正含义只有模糊的认识,他们在数据里翻找、在不同假设之间挑选,直到找到一个显著的 p 值,于是就能发表论文。

真正的问题是,人们希望从 p 值中得到超出它本身的信息。我们希望 p 值代表的是:

结果由随机性造成的概率

我们希望这个值越低越好,这样就可以得出“我们证明了某件事”的结论。许多期刊编辑就是这样解读 p 值的。但实际上 p 值真正表示的是:

在一个机会模型成立的条件下,得到像观测结果一样极端结果的概率

两者区别微妙但真实。显著的 p 值并不能像它表面承诺的那样把你推向“证明”的地步。当理解了 p 值的真正含义时,“统计显著”这一结论的逻辑基础其实更弱。

2016 年 3 月,美国统计学会(ASA)经过内部长时间讨论,发布了关于 p 值使用的警示声明,揭示了人们对 p 值的误解。ASA 声明强调了供研究者和期刊编辑参考的六条原则:

  1. p 值可以反映数据与某个特定统计模型的不兼容程度。
  2. p 值不能衡量所研究假设为真的概率,也不能衡量数据纯粹由随机性产生的概率。
  3. 科学结论和商业或政策决策不应仅仅基于 p 值是否超过某个阈值。
  4. 做出恰当推断需要全面报告和透明度。
  5. p 值或统计显著性不能衡量效应大小或结果的重要性。
  6. 单靠 p 值无法提供关于模型或假设的有力证据。

实际显著性:Practical significance

即使一个结果在统计上显著,这也并不意味着它在实际中具有显著性。如果样本量足够大,即便是一个毫无实际意义的微小差异,也可能在统计上显著。大样本能确保即便是小且无意义的效应,也足以排除“随机”作为解释。然而,排除了随机性并不会神奇地使一个本质上无关紧要的结果变得重要。

第一类错误与第二类错误

在评估统计显著性时,可能发生两种错误:

  • 第一类错误(Type 1 error):当效应实际上只是随机造成的,却错误地得出效应真实存在的结论;
  • 第二类错误(Type 2 error):当效应实际上是真实存在的,却错误地得出效应不存在(即归因于随机性)的结论。

实际上,第二类错误与其说是错误,不如说是样本量太小,无法检测到效应。当 p 值未达到统计显著性(例如大于 5%)时,我们真正的意思是“效应未被证明”。可能换成更大的样本就会得到更小的 p 值。

显著性检验(也称假设检验)的基本功能是防止被随机性所欺骗,因此它们通常设计成最小化第一类错误。

数据科学与 p 值

数据科学家的工作通常不是为了在科学期刊上发表,因此围绕 p 值价值的争论更多是学术问题。对于数据科学家来说,p 值在想要知道某个看起来有趣且有用的模型结果是否处于正常随机波动范围内时,是一个有用的度量指标。作为实验中的决策工具,p 值不应被视为“决定性”的,而仅仅是为决策提供信息的一个参考点。例如,在某些统计或机器学习模型中,p 值有时作为中间输入使用——某个特征是否被包含在模型中,可能取决于它的 p 值。

关键要点

  • 显著性检验用于判断观测到的效应是否在零假设模型的随机波动范围内;
  • p 值是在零假设模型成立的前提下,得到像观测结果一样极端结果的概率;
  • α 值是零假设随机模型中“异常程度”的阈值;
  • 显著性检验在正式研究报告中比在数据科学中更相关(但近年即便在研究中也逐渐淡化)。

t 检验

显著性检验有很多种,具体取决于数据是计数数据还是测量数据、样本有多少以及测量的是什么。非常常见的一种是 t 检验,它以 Student’s t 分布命名,最初由 W. S. Gosset 提出,用于近似单一样本均值的分布(参见第 75 页“Student’s t 分布”)。

t 检验关键术语

  • 检验统计量(Test statistic):衡量差异或感兴趣效应的指标;
  • t 统计量(t-statistic):均值等常见检验统计量的标准化版本;
  • t 分布(t-distribution):一个参考分布(此处来源于零假设),用来与观测到的 t 统计量进行比较。

所有显著性检验都要求你指定一个检验统计量,用来衡量你感兴趣的效应,并帮助你判断观测到的效应是否处于正常随机波动的范围内。在重抽样检验中(见第 97 页“置换检验”的讨论),数据的量纲无关紧要。你从数据本身创建参考(零假设)分布,并直接使用该检验统计量。

在 20 世纪 20–30 年代,统计假设检验正在发展之时,随机打乱数据上千次以做重抽样检验在计算上并不可行。统计学家发现,一个对置换(打乱)分布的良好近似是 t 检验,它基于 Gosset 的 t 分布。它被用于非常常见的两样本比较——A/B 测试——当数据是数值型时尤为常见。但要使 t 分布在不考虑量纲的情况下使用,必须采用检验统计量的标准化形式。

一本经典的统计教材在此处会展示各种结合 Gosset 分布的公式,并演示如何将数据标准化以与标准 t 分布进行比较。这里没有列出这些公式,因为所有统计软件(包括 R 和 Python)都内置了实现这些公式的命令。在 R 中,对应的函数是 t.test

1
2
3
4
5
6
7
8
9
10
11
> t.test(Time ~ Page, data=session_times, alternative='less')

Welch Two Sample t-test
data: Time by Page
t = -1.0983, df = 27.693, p-value = 0.1408
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf 19.59674
sample estimates:
mean in group Page A mean in group Page B
126.3333 162.0000

在 Python 中可以使用 scipy.stats.ttest_ind

1
2
3
4
5
res = stats.ttest_ind(
session_times[session_times.Page == 'Page A'].Time,
session_times[session_times.Page == 'Page B'].Time,
equal_var=False)
print(f'p-value for single sided test: {res.pvalue / 2:.4f}')

备择假设是:Page A 的会话时间均值小于 Page B 的会话时间。p 值为 0.1408,与置换检验得到的 p 值 0.121 和 0.126(见第 98 页“示例:网页黏性”)相当接近。

在重抽样模式下,我们根据观测数据和待检验的假设来构造解,不必担心数据是数值型还是二元型,样本量是否均衡、样本方差、以及其他各种因素。而在公式的世界里,存在很多变体,可能令人眼花缭乱。统计学家需要在那个世界中导航并掌握它的地图,但数据科学家不必如此——他们通常无需像准备学术论文的研究者那样在假设检验和置信区间的细节上费尽心思。

关键要点

  • 在计算机出现之前,重抽样检验并不实际,统计学家使用标准参考分布。
  • 然后可以将检验统计量标准化并与参考分布进行比较。
  • 一个广泛使用的标准化统计量就是 t 统计量

多重检验

Multiple Testing

正如我们之前提到过的,统计学里有句俗语:“折磨数据足够久,它终会招供。”这意味着,如果你从足够多的角度看数据、问足够多的问题,几乎总能找到某个“统计显著”的效应。

例如,如果你有 20 个预测变量和 1 个结果变量,且它们都是随机生成的,那么如果你按 α = 0.05 的显著性水平做 20 次显著性检验,至少有一个预测变量(错误地)显得统计显著的概率非常大。正如前文所述,这就是Ⅰ类错误。你可以这样计算这个概率:先求出所有检验都正确地不显著的概率。一个检验正确地不显著的概率是 0.95,那么 20 个都正确地不显著的概率就是0.95 × 0.95 × 0.95… = 0.95²⁰ = 0.36。至少一个预测变量(错误地)显著的概率就是这个概率的反面,即1 – (所有都不显著的概率)= 0.64。这被称为α 膨胀(alpha inflation)

这一问题与数据挖掘中的过拟合有关,也就是“把模型拟合到噪声上”。你加入的变量越多、运行的模型越多,就越可能仅凭偶然性得出“显著”的结果。

多重检验关键术语

  • Ⅰ类错误(Type 1 error) 错误地认为某个效应具有统计显著性。

  • 假发现率(False discovery rate) 在多次检验中犯Ⅰ类错误的比例。

  • α 膨胀(Alpha inflation) 随着进行更多检验,α(Ⅰ类错误概率)不断上升的多重检验现象。

  • p 值调整(Adjustment of p-values) 针对同一数据做多次检验时进行的校正。

  • 过拟合(Overfitting) 拟合到了噪声。

在监督学习任务中,使用一个模型未见过的数据集(保留集/验证集)来评估模型,可以降低这种风险。而在没有标注保留集的统计或机器学习任务中,基于统计噪声得出结论的风险依然存在。

在统计学中,有一些方法在特定情境下用来处理这个问题。例如,如果你在比较多个处理组的结果,你可能会问多个问题。对于 A–C 三个处理,你可能会问:

  • A 是否不同于 B?
  • B 是否不同于 C?
  • A 是否不同于 C?

或者,在临床试验中,你可能想在多个阶段观察某种疗法的结果。每一次提问,都会增加被偶然性“愚弄”的概率。统计学中的调整方法可以通过比单个假设检验更严格地设定显著性门槛来补偿这种风险。这些调整方法通常涉及根据检验的次数“分割 α”。这会导致每个检验的 α 更小(即显著性标准更严格)。

  • Bonferroni 调整:简单地把 α 除以比较的次数。
  • Tukey 的“诚实显著性差异”(HSD):用于比较多个组均值。这个检验针对组均值之间的最大差异,将其与一个基于 t 分布的基准进行比较(大致相当于把所有值混在一起,重新抽取与原来组大小相同的样本,然后找出这些重抽样组均值之间的最大差异)。

然而,多重比较的问题超出了这些高度结构化的情形,并与反复“淘洗”数据的现象有关,这正是那句“折磨数据”的俗语的由来。换句话说,给定足够复杂的数据,如果你还没发现有趣的东西,那只是因为你还没看得足够久、问得足够多。如今可用的数据比以往任何时候都多,期刊文章的数量在 2002–2010 年间几乎翻了一番。这给了我们大量机会在数据中发现有趣的东西,同时也带来了多重性问题,例如:

  • 检查多个组间的成对差异
  • 观察多个子组的结果(“整体治疗效果不显著,但我们在 30 岁以下未婚女性中发现了效果”)
  • 尝试大量统计模型
  • 在模型中包含大量变量
  • 提出许多不同的问题(即不同的可能结果)

通用注解

假发现率(False Discovery Rate, FDR)**:“假发现率”这一术语最初用来描述一组假设检验中,错误地识别出显著效应的比例。随着基因组学研究的兴起,它变得尤其有用,因为在基因测序项目中可能要进行海量的统计检验。在这些情境中,这个术语适用于整个检验流程,而单个错误“发现”指的是某次假设检验的结果(例如比较两份样本)。研究人员希望设置检验过程的参数,以在指定水平上控制假发现率。该术语还被用在数据挖掘的分类问题中:它表示在类别 1 的预测中出现的误分类率。换句话说,它是“发现”是假的概率(把一条记录标为“1”时,它实际上是假的概率)。这种情形通常对应 0 类丰富、1 类稀少而有趣的情况(参见第 5 章“稀有类别问题”)。

由于包括“多重性”在内的多种原因,更多的研究并不一定意味着更好的研究。例如,制药公司拜耳在 2011 年发现,当它试图重复 67 项科学研究时,只有 14 项能完全复现,近三分之二完全无法复现。

无论如何,高度结构化的统计检验的各种调整程序过于具体和僵化,不适合数据科学的普遍用途。数据科学家在多重性问题上的要点是:

  • 对于预测建模,通过交叉验证(见第 155 页“交叉验证”)和使用保留样本,可以减轻得到一个看似有效、但实际上主要源自随机性的虚幻模型的风险。
  • 对于没有带标签的保留集来检查模型的其他过程,你必须依赖: — 意识到查询和操作数据越多,随机性起作用的可能性就越大; — 使用重采样和模拟的启发式方法提供随机基准,与观测结果进行比较。

关键要点

  • 研究或数据挖掘项目中的多重性(多重比较、许多变量、许多模型等)增加了仅因偶然性而得出显著结论的风险。
  • 在涉及多次统计比较(即多重显著性检验)的情形下,有统计学上的调整程序。
  • 在数据挖掘情境中,使用带标签结果变量的保留样本可以帮助避免误导性结果。

自由度

Degrees of Freedom

在许多统计检验和概率分布的文档及设置中,你会看到“自由度”一词。该概念应用于从样本数据计算得来的统计量,指可自由变化的数值个数。例如,如果你知道一个包含 10 个值的样本的均值,那么自由度就是 9(因为一旦你知道其中 9 个样本值,第 10 个值就可以计算出来,不再自由)。在许多概率分布中,自由度参数会影响分布的形状。

自由度的数量是许多统计检验的输入参数。例如,自由度就是在方差和标准差计算中分母 n – 1 的名称。为什么这很重要?当你用样本估计总体方差时,如果分母用 n,你会得到一个略微向下偏倚的估计;如果用 n – 1,估计就没有这种偏倚。

自由度关键术语

  • n 或样本量:数据中的观测数(也称行数或记录数)。
  • d.f.(自由度)

传统统计课程或教材中很大一部分都在讲各种标准假设检验(t 检验、F 检验等)。当样本统计量被标准化以应用于传统统计公式时,自由度是标准化计算的一部分,以确保你的标准化数据与相应的参考分布(t 分布、F 分布等)相匹配。

对数据科学来说重要吗?其实不太重要,至少在显著性检验的语境中不是这样。一方面,数据科学中正式的统计检验使用得很有限;另一方面,数据规模通常足够大,因此对于数据科学家来说,分母是 n 还是 n – 1 几乎没有实质差别(当 n 变大时,用 n 作为分母所产生的偏倚会消失)。

不过,有一个情境与之相关:在回归(包括逻辑回归)中使用分解变量时。如果存在完全冗余的预测变量,一些回归算法会出错。这种情况最常见于把分类变量分解为二元指示变量(哑变量)。以“星期几”为例。虽然一周有七天,但“星期几”这个变量只有六个自由度。比如,一旦你知道“不是周一到周六”,就可以推知“必然是周日”。因此,如果你包含了周一到周六的指示变量,再加上周日,就会因为多重共线性导致回归失败。

关键要点

  • 自由度(d.f.)是标准化检验统计量以便与参考分布(t 分布、F 分布等)比较的计算的一部分。
  • 自由度概念是把分类变量分解为 n – 1 个指示(哑)变量以进行回归的理论基础(避免多重共线性)。

ANOVA(方差分析)

假设我们不是进行 A/B 测试,而是要比较多个组,比如 A/B/C/D,每个组都有数值数据。用来检验各组之间是否存在统计学显著差异的统计程序称为方差分析(analysis of variance,简称 ANOVA)。

ANOVA 的关键术语

  • 成对比较(Pairwise comparison) 在多个组中对任意两组之间(如均值)进行的假设检验。

  • 总体检验(Omnibus test) 对多个组均值总体方差进行的单一假设检验。

  • 方差分解(Decomposition of variance) 将对单个值有贡献的成分分离出来(例如相对于总体均值、相对于处理均值和相对于残差误差的成分)。

  • F 统计量(F-statistic) 标准化的统计量,衡量组均值之间的差异是否超过随机模型下可能出现的差异程度。

  • SS(平方和,Sum of Squares) 指的是相对于某个平均值的偏差(deviations)平方之和。

T3.3

表 3-3 显示了四个网页的“黏性”(stickiness),定义为访问者停留在页面上的秒数。四个页面是随机分配的,每位访客只会看到其中一个页面。每个页面有 5 位访客,表 3-3 中的每一列都是一组独立的数据。例如,第 1 页的第 1 位访客与第 2 页的第 1 位访客之间没有关联。注意在这种网页测试中,我们无法完全实现经典的随机抽样设计(即每位访客从某个庞大总体中随机选取),只能按访客到来的顺序分配。访客可能因一天中的时间、星期几、季节、网络条件、所用设备等而系统性地不同,在回顾实验结果时要把这些因素视为潜在偏差。

F3.6

现在我们遇到了一个难题(见图 3-6)。当只比较两个组时,我们只需看两组均值的差异就行。但当有四个均值时,可能的成对比较有六种:

  • 第 1 页与第 2 页比较
  • 第 1 页与第 3 页比较
  • 第 1 页与第 4 页比较
  • 第 2 页与第 3 页比较
  • 第 2 页与第 4 页比较
  • 第 3 页与第 4 页比较

成对比较越多,就越有可能因随机性而得出错误的结论(见第 112 页“多重检验”)。与其担心页面间所有可能的比较,不如做一个总体检验,回答这样的问题:“所有页面的真实‘黏性’是否可能相同,观察到的差异是否仅仅是由于相同会话时间随机分配给四个页面造成的?”

用来检验这一点的程序就是 ANOVA。它的原理可以用下面的重抽样程序来说明(这里针对网页黏性的 A/B/C/D 测试):

  1. 将所有数据合并到一个整体中;
  2. 随机打乱后抽取四个包含 5 个值的重样本;
  3. 记录四个组的均值;
  4. 记录四个组均值之间的方差;
  5. 重复步骤 2–4 多次(比如 1,000 次);

重抽样后有多少次方差超过了观察到的方差?这就是 p 值。

这种置换检验比第 97 页“置换检验”所用的要复杂一些。幸运的是,R 中 lmPerm 包的 aovp 函数可以计算这种情况的置换检验:

1
2
3
4
5
6
7
8
9
library(lmPerm)
summary(aovp(Time ~ Page, data=four_sessions))
[1] "Settings: unique SS "
Component 1 :
Df R Sum Sq R Mean Sq Iter Pr(Prob)
Page 3 831.4 277.13 3104 0.09278 .
Residuals 16 1618.4 101.15
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

这里 Pr(Prob) 给出的 p 值是 0.09278。换句话说,假设底层黏性相同,9.3% 的时候四个页面的响应率可能会像实际观察到的那样不同,仅仅是因为随机性。这种程度的不可能性没有达到传统统计学的 5% 阈值,因此我们得出结论:四个页面之间的差异可能是由随机性造成的。

Iter 列出了置换检验的迭代次数。其他列对应于传统 ANOVA 表,下一节会详细介绍。

在 Python 中,我们可以用下面的代码计算这种置换检验:

1
2
3
4
5
6
7
8
9
10
11
observed_variance = four_sessions.groupby('Page').mean().var()[0]
print('Observed means:', four_sessions.groupby('Page').mean().values.ravel())
print('Variance:', observed_variance)

def perm_test(df):
df = df.copy()
df['Time'] = np.random.permutation(df['Time'].values)
return df.groupby('Page').mean().var()[0]

perm_variance = [perm_test(four_sessions) for _ in range(3000)]
print('Pr(Prob)', np.mean([var > observed_variance for var in perm_variance]))

F 统计量

就像 t 检验可以用来替代置换检验以比较两个组的均值一样,基于 F 统计量也有一个用于方差分析(ANOVA)的统计检验。F 统计量是“组均值之间的方差(the variance across group means)”(即处理效应)与“残差误差的方差(the variance due to residual error)”的比值。这个比值越大,结果就越显著。如果数据服从正态分布,统计理论规定该统计量应当服从某种特定分布。基于此,我们就可以计算 p 值。

在 R 中,可以用 aov 函数计算 ANOVA 表:

1
2
3
4
5
6
> summary(aov(Time ~ Page, data=four_sessions))
Df Sum Sq Mean Sq F value Pr(>F)
Page 3 831.4 277.1 2.74 0.0776 .
Residuals 16 1618.4 101.2
---
Signif. codes: 0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

statsmodels 包在 Python 中也提供了 ANOVA 的实现:

1
2
3
model = smf.ols('Time ~ Page', data=four_sessions).fit()
aov_table = sm.stats.anova_lm(model)
aov_table

Python 输出与 R 的输出几乎相同。

Df 表示“自由度(degrees of freedom)”,Sum Sq 表示“平方和(sum of squares)”,Mean Sq 表示“均方(mean squares,均方偏差)”,F value 表示 F 统计量。对于总平均值,平方和是总平均值偏离 0 的平方乘以 20(观测数)。总平均值的自由度按定义是 1。

对于处理均值,自由度是 3(当三个值确定后,再加上总平均值就固定了,另一个处理均值不能再变)。处理均值的平方和是各处理均值与总平均值之间偏差的平方和。

对于残差,自由度是 20(所有观测可以变动),平方和是各观测值与处理均值之间差异的平方和。均方(MS)就是平方和除以自由度。F 统计量是 MS(处理)/MS(误差)。F 值因此只依赖这个比率,并可以与标准 F 分布比较,以判断处理均值之间的差异是否超出了随机波动的期望范围。

通用注解

方差分解**:数据集中每个观测值都可以视为不同组成部分的和。对于数据集中的任意观测值,我们可以将其分解为总平均值、处理效应和残差误差。这就是所谓的“方差分解”:

  1. 从总平均值开始(网页粘性数据为 173.75)。
  2. 加上处理效应,可能是负的(自变量 = 网页)。
  3. 加上残差误差(residual error,),可能是负的。

因此,A/B/C/D 测试表左上角值的方差分解如下:

  1. 从总平均值开始:173.75。
  2. 加上处理(组)效应:–1.75(172 – 173.75)。
  3. 加上残差:–8(164 – 172)。
  4. 得到:164

双因素 ANOVA

Two-Way ANOVA

刚才描述的 A/B/C/D 测试是一种“单因素”ANOVA,即只有一个因子(组)在变化。如果有第二个因子——例如“周末 vs 工作日”——并且在每种组合(A 组周末、A 组工作日、B 组周末等)上收集数据,这就是“双因素 ANOVA”。处理方法与单因素 ANOVA 类似,但需要识别“交互效应”。在确定总平均效应和处理效应后,我们把每组的周末和工作日观测分开,计算这些子集平均值与处理平均值的差异。

可以看出,ANOVA 和双因素 ANOVA 是迈向完整统计建模(如回归和逻辑回归)的第一步,在完整模型中可以同时刻画多个因素及其效应(见第 4 章)。

关键概念

  • ANOVA 是一种分析多个组实验结果的统计方法。
  • 它是 A/B 测试等类似程序的推广,用于评估组间总体差异是否在随机波动范围内。
  • ANOVA 的一个有用结果是识别与组处理、交互效应和误差相关的方差组成部分。

卡方检验

Chi-Square Test

在网络测试中,人们经常会超越简单的 A/B 测试,同时检验多种处理方案。卡方检验用于计数型数据,以检验其与某个期望分布的拟合程度。在统计实践中,卡方统计量最常见的用途是用于 \(r \times c\) 列联表,以评估变量之间独立性的原假设是否合理(参见第 80 页“卡方分布”)。

卡方检验最早由 Karl Pearson 于 1900 年提出。“chi”这个词来源于 Pearson 在文章中使用的希腊字母 Χ。

卡方检验关键术语

  • 卡方统计量(Chi-square statistic) 衡量一组观测数据偏离期望程度的指标。

  • 期望值(Expectation or expected) 在某个假设(通常是原假设)下,我们期望数据呈现的结果。

通用注解 \(r \times c\) 表示“行 × 列”。例如,一个 \(2 \times 3\) 的表格表示两行三列。

卡方检验:一种重抽样方法

假设你在测试三条不同的标题——A、B 和 C——并分别在 1,000 位访客上运行测试,结果如表 3-4 所示。

T3.4

标题之间显然存在差异。标题 A 的点击率几乎是 B 的两倍。但实际数字很小。我们可以通过重抽样方法来检验点击率的差异是否超出了偶然性造成的范围。对于这个检验,我们需要“期望”的点击分布。在原假设下,假定三条标题的点击率相同,总体点击率为 34/3,000。在这种假设下,我们的列联表如表 3-5 所示。

T3.5

Pearson 残差(Pearson residual) 定义为:

\[ R = \frac{\text{Observed} - \text{Expected}}{ \sqrt{\text{Expected}}} \]

\(R\) 衡量实际计数与期望计数的差异程度(见表 3-6)。

T3.6

卡方统计量(Chi-square statistic) 定义为所有 Pearson 残差平方和:

\[ \chi^2 = \sum_{i}^{r}\sum_{j}^{c}R_{ij}^2 \]

其中 \(r\)\(c\) 分别为行数和列数。本例中的卡方统计量为 1.666。这是否超过了偶然模型中合理出现的范围呢?

我们可以用以下重抽样算法进行检验:

  1. 构建一个包含 34 个“1”(点击)和 2,966 个“0”(未点击)的总体。
  2. 打乱顺序,分别抽取三个 1,000 人样本,并统计每个样本的点击数。
  3. 计算抽样计数与期望计数之间的平方差并求和。
  4. 重复步骤 2 和 3,例如 1,000 次。
  5. 观察重抽样平方偏差和超过实际观测值的频率。这就是 p 值。

chisq.test 函数可用于在 R 中计算重抽样的卡方统计量。对于点击数据,卡方检验如下:

1
2
3
4
> chisq.test(clicks, simulate.p.value=TRUE)
Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)
data: clicks
X-squared = 1.6659, df = NA, p-value = 0.4853

该检验结果表明,这个观测值很可能只是随机性所致。

在 Python 中运行置换检验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
box = [1] * 34
box.extend([0] * 2966)
random.shuffle(box)

def chi2(observed, expected):
pearson_residuals = []
for row, expect in zip(observed, expected):
pearson_residuals.append([(observe - expect) ** 2 / expect
for observe in row])
# 返回平方和
return np.sum(pearson_residuals)

expected_clicks = 34 / 3
expected_noclicks = 1000 - expected_clicks
expected = [34 / 3, 1000 - 34 / 3]
chi2observed = chi2(clicks.values, expected)

def perm_fun(box):
sample_clicks = [sum(random.sample(box, 1000)),
sum(random.sample(box, 1000)),
sum(random.sample(box, 1000))]
sample_noclicks = [1000 - n for n in sample_clicks]
return chi2([sample_clicks, sample_noclicks], expected)

perm_chi2 = [perm_fun(box) for _ in range(2000)]

resampled_p_value = sum(perm_chi2 > chi2observed) / len(perm_chi2)
print(f'Observed chi2: {chi2observed:.4f}')
print(f'Resampled p-value: {resampled_p_value:.4f}')

卡方检验:统计理论

Chi-Square Test: Statistical Theory

渐近统计理论表明,卡方统计量的分布可以用卡方分布来近似(参见第 80 页“卡方分布”)。合适的标准卡方分布由自由度决定(参见第 116 页“自由度”)。对于一个列联表,自由度与行数 \(r\) 和列数 \(c\) 的关系如下:

\[ \text{degrees of freedom} = (r-1)\times(c-1) \]

卡方分布通常是偏斜的,具有右长尾;见图 3-7(自由度分别为 1、2、5 和 20 的分布)。

F3.7

观测统计量在卡方分布上的位置越靠右,p 值就越小。

chisq.test 函数可使用卡方分布作为参考计算 p 值:

1
2
3
4
> chisq.test(clicks, simulate.p.value=FALSE)
Pearson's Chi-squared test
data: clicks
X-squared = 1.6659, df = 2, p-value = 0.4348

在 Python 中,使用 scipy.stats.chi2_contingency 函数:

1
2
3
chisq, pvalue, df, expected = stats.chi2_contingency(clicks)
print(f'Observed chi2: {chi2observed:.4f}')
print(f'p-value: {pvalue:.4f}')

这个 p 值比重抽样的 p 值稍小;这是因为卡方分布只是对统计量实际分布的一种近似。

费舍尔确切检验

Fisher’s Exact Test

卡方分布是前面描述的洗牌重抽样检验(shuffled resampling test)的一个很好的近似,但在计数值非常小(个位数,尤其是 5 或更少)时就不再适用。在这种情况下,重抽样程序会给出更精确的 p 值。事实上,大多数统计软件都有一个程序,能够枚举所有可能发生的重新排列(置换),统计它们的频率,并精确地判断观察到的结果有多极端。这种方法被称为 费舍尔确切检验,以纪念著名统计学家 R. A. Fisher。R 语言中调用费舍尔确切检验的代码非常简单:

1
2
3
4
5
6
> fisher.test(clicks)

Fisher's Exact Test for Count Data
data: clicks
p-value = 0.4824
alternative hypothesis: two.sided

得到的 p 值与使用重抽样方法得到的 0.4853 非常接近。

在某些情况下,如果部分计数非常小而其他计数很大(例如转化率中的分母很大),则可能需要执行 洗牌置换检验(a shuffled permutation test ),而不是完整的确切检验,因为计算所有可能的置换会非常困难。上述 R 函数提供了几个参数来控制是否使用近似方法(simulate.p.value = TRUE 或 FALSE)、迭代次数(B = ...),以及一个计算约束参数(workspace = ...),用于限制精确结果的计算规模。

Python 中并没有易用的费舍尔确切检验实现。

科学造假的检测

一个有趣的案例是塔夫茨大学(Tufts University)的研究员 Thereza Imanishi-Kari,她在 1991 年被指控在研究中捏造数据。美国国会议员 John Dingell 介入此事,最终导致她的同事、当时任洛克菲勒大学校长的 David Baltimore 辞职。

案件中的一个关键环节涉及其实验数据中数字的预期分布。她的每条观测数据包含多个数字,调查人员将注意力集中在中间的数字(忽略首位和末位数字)。按统计规律,这些中间数字应该遵循均匀随机分布,即每个数字出现的概率相等(首位数字可能偏向某个值,而末位数字可能受到四舍五入的影响)。表 3-7 给出了实际数据中各中间数字的频数。

T3.7

图 3-8 显示了 315 个数字的分布,看起来明显不像是随机的。调查人员计算了与期望值的偏差(理论上严格均匀分布下每个数字应出现 31.5 次),并使用卡方检验(同样也可以使用重抽样方法)证明实际分布远远超出正常的随机波动范围,表明数据可能被捏造。

(注:Imanishi-Kari 最终在经历漫长调查后被澄清无罪。)

F3.8

个人注:辩护方的观点: 伊马尼西-卡里及其支持者(包括大卫·巴尔的摩)坚称,数据中的不一致性是由于实验过程中的疏忽、不规范的记录习惯,以及当时生物学实验固有的可变性所致。他们认为,指控方过于依赖统计学分析,而忽略了生物学实验的复杂性和现实情况。此外,他们还批评了调查程序缺乏正当性,对伊马尼西-卡里不公平

与数据科学的相关性

Relevance for Data Science

当你想知道某个效应是否真实存在,还是只是偶然产生时,就会用到 卡方检验(chi-square test)费舍尔确切检验(Fisher’s exact test)。在大多数经典统计学的应用中,卡方检验的作用是建立统计显著性,通常这是研究或实验能够发表前所必须具备的条件。 但对于数据科学家来说,这一点并不那么重要。在大多数数据科学实验中(无论是 A/B 还是 A/B/C……),目标不仅仅是建立统计显著性,而是找到最佳处理方案。为此,多臂老虎机算法(见第 131 页 “Multi-Arm Bandit Algorithm”)提供了一个更完整的解决方案。

卡方检验(特别是费舍尔确切检验)在数据科学中的一个应用,是用来确定网络实验的适当样本量。这些实验往往点击率极低,即使有成千上万的曝光,计数率可能仍然太小,无法在实验中得出确定结论。在这种情况下,费舍尔确切检验、卡方检验以及其他检验可以作为功效(power)和样本量计算(见第 135 页 “Power and Sample Size”)的一部分发挥作用。

研究人员广泛使用卡方检验,以寻找那个“难以捉摸的统计显著性 p 值”来支撑论文发表。而在数据科学应用中,卡方检验或类似的重抽样模拟更多地被用作筛选工具,用来判断某个效应或特征是否值得进一步研究,而不是作为形式化的显著性检验。

例如,它们可以用在空间统计和制图中,以判断空间数据是否符合某个指定的零假设分布(比如:犯罪是否集中在某个区域,超出了随机机会所允许的程度?)。它们也可以用在机器学习中的自动化特征选择,评估不同特征上的类别分布,从而识别某些特征在某个类别上出现率异常高或异常低的情况,这种偏差无法用随机波动解释。

关键要点(Key Ideas)

  • 统计学中的一个常见做法是检验观测到的数据计数是否与独立性假设一致(例如,购买某种商品的倾向是否与性别无关)。
  • 卡方分布是一个参照分布(体现了独立性假设),用来与观测数据计算出的卡方统计量进行比较。

多臂老虎机算法

Multi-Arm Bandit Algorithm

多臂老虎机为测试(尤其是网络测试)提供了一种方法,它能比传统的实验设计与统计方法实现显式优化更快速的决策

多臂老虎机关键术语

  • 多臂老虎机(Multi-arm bandit) 一个假想的老虎机,顾客可以选择多个拉杆,每个拉杆的收益不同,这里用作多处理(多方案)实验的类比。

  • 拉杆(Arm) 实验中的一种处理(例如:“网页测试中的标题 A”)。

  • 中奖(Win) 在实验中对应老虎机中奖的事件(例如:“顾客点击了链接”)。

传统的 A/B 测试是在实验中按预先设定的设计收集数据,用来回答一个特定的问题,例如:“哪一个更好,处理 A 还是处理 B?” 这种方法默认一旦我们得到问题的答案,实验就结束,然后基于结果采取行动。

但这种方法存在几个问题:

  1. 结果可能不明确: “效果未被证明。” 换句话说,实验结果可能暗示有某种效应,但如果真的存在效应,我们的样本量可能不够大,无法按照传统统计标准“证明”它。那么我们该如何决策?

  2. 想在实验结束前利用已有结果: 我们可能希望在实验完全结束前就根据已经获得的数据做出决策。

  3. 想在实验后期调整或改变方案: 根据实验结束后新出现的数据,我们可能希望更换策略或做其他尝试。

传统的实验与假设检验方法源自 20 世纪 20 年代,相当不灵活。而计算机算力和软件的发展使得更强大、更灵活的方法成为可能。此外,数据科学(以及商业领域)并不太在乎统计显著性,而更关注优化整体努力与结果

老虎机算法(Bandit algorithms)在网络测试中非常流行,它允许你同时测试多种处理方案,并比传统统计设计更快得出结论。它的名字来自赌博中的老虎机(slot machine,又称“单臂强盗”one-armed bandit,因为它们“稳定地”从赌徒那里赢钱)。如果你想象一台有多个拉杆的老虎机,每个拉杆的支付率不同,你就得到了一个“多臂老虎机”,也就是该算法的全称。

你的目标是赢取尽可能多的钱,更具体地说,尽早识别并选定获胜的拉杆。困难在于你并不知道每个拉杆的总体中奖率——你只知道每次拉动时的个别结果。假设每次中奖的奖金都是一样的(无论拉哪个杆),不同之处在于中奖的概率。再假设你最初每个拉杆都尝试 50 次,得到如下结果:

  • 拉杆 A:50 次中赢了 10 次
  • 拉杆 B:50 次中赢了 2 次
  • 拉杆 C:50 次中赢了 4 次

一种极端做法是说:“看起来 A 拉杆是赢家——那我们就不再尝试其他拉杆,专注于 A 吧。”这样可以充分利用初始试验的信息。如果 A 真的更好,我们可以早早享受到这一优势。但另一方面,如果 B 或 C 实际上更好,我们就失去了发现它们的机会。

另一种极端做法是说:“这些结果都可能是偶然的——那就继续平等地拉所有拉杆吧。”这样可以最大限度地给 A 之外的备选方案展示自己的机会。然而,在此过程中,我们也在持续投放看起来较差的处理。那我们要容忍多久?

多臂老虎机算法采取折中方法: 我们开始更频繁地拉 A,以利用它的明显优势,但并不放弃 B 和 C,只是减少拉它们的次数。如果 A 持续优于其他,我们就继续把资源(拉动次数)从 B、C 转移到 A;如果相反,C 开始表现更好而 A 变差,我们可以把拉动次数从 A 转回 C。如果其中某个在初始试验中因偶然性被“埋没”的处理其实比 A 优秀,那么通过后续测试它就有机会显现出来。

现在把这个思路应用到网络测试: 不再是多个老虎机拉杆,而是多个优惠、标题、颜色等在网站上进行测试。顾客要么点击(商家“中奖”),要么不点击。最初,所有方案都是随机且均等展示。然而,如果某个方案开始优于其他方案,它就可以被更频繁地展示(“拉”)。但是,修改展示比例的算法参数应该怎么设定?应该改成什么比例?什么时候改?

一个简单的算法:A/B 测试中的 ε-贪婪算法(epsilon-greedy)

  1. 生成一个 0 到 1 之间均匀分布的随机数。

  2. 如果这个数在 0 和 ε 之间(ε 是 0 到 1 之间的数,通常比较小),掷一枚公平的硬币(50/50 概率):

    • 如果正面,展示方案 A;
    • 如果反面,展示方案 B。
  3. 如果这个数 ≥ ε,就展示目前响应率最高的方案。

ε 是唯一控制这个算法的参数:

  • 如果 ε = 1,我们得到的是标准的简单 A/B 实验(每个用户在 A 和 B 间随机分配)。
  • 如果 ε = 0,我们得到的是纯贪婪算法——始终选择当前表现最好的方案(局部最优),不再继续探索,仅把用户(网站访问者)分配到当前最佳处理。

更高级的算法:Thompson 采样(Thompson’s sampling)

这个方法在每一步“抽样”(拉动拉杆),以最大化选择最佳拉杆的概率。当然你并不知道哪一个是最好的——这正是问题所在!—— 但随着每次抽样后你观察到收益,你获得越来越多的信息。Thompson 采样采用贝叶斯方法:最初假设奖励的先验分布(使用 Beta 分布,这是在贝叶斯问题中指定先验信息的常见机制)。随着每次抽样积累信息,这些信息被更新,从而使下一次抽样在选择正确拉杆时更优化。

多臂老虎机算法可以高效地处理 3 个以上处理方案,并趋向于最优的“最佳”选择。而在传统统计检验程序中,处理 3 个以上方案的决策复杂度远远超过传统 A/B 测试,这使得多臂老虎机算法的优势更大。

关键要点

  • 传统的 A/B 测试设想的是随机抽样过程,这可能导致对劣质处理的过多暴露。
  • 多臂老虎机算法则相反,它在实验过程中利用获得的信息,动态调整抽样过程,减少对劣质处理的分配频率。
  • 它们还能高效处理多于两种处理方案的情况。
  • 有多种算法可以将抽样概率从较差处理方案转向(假定的)较优方案。

检验效能与样本量

Power and Sample Size

如果你要运行一次网页测试,如何决定测试要运行多久(即每个处理需要多少展示次数)?尽管许多网页测试指南都有所提及,但实际上没有什么普遍适用的“标准答案”——主要取决于你所期望达成的目标出现的频率。

检验效能与样本量的关键术语

  • 效应量(Effect size) 你希望在统计检验中能够检测到的最小效应,例如“点击率提高 20%”。

  • 检验效能(Power) 在给定样本量下检测到指定效应量的概率。

  • 显著性水平(Significance level) 检验所采用的统计显著性水平。

在进行样本量的统计计算时,其中一步要问:“假设检验能否真正揭示处理 A 和 B 之间的差异?”假设检验的结果(即 p 值)取决于处理 A 和 B 之间真实的差异,也取决于抽样的运气——即实验中谁被分配到哪个组。但直觉上,A 和 B 之间真实差异越大,我们越有可能在实验中发现它;差异越小,就需要更多数据来检测它。比如,要区分棒球打击率 0.350 和 0.200 的打者,不需要太多打席就可以;而要区分 0.300 和 0.280 的打者,就需要更多打席。

检验效能是指在给定样本特征(样本量与变异度)下,检测到指定效应量的概率。例如,我们可以(假设性地)说:在 25 次打席中区分打击率 0.330 和 0.200 的打者的概率是 0.75。这里的效应量是 0.130,而“检测到”意味着假设检验能够拒绝“无差异”的原假设,并得出“存在真实效应”的结论。所以这个包含两位打者、每人 25 次打席(n=25)、效应量 0.130 的实验,其(假设的)检验效能为 0.75,即 75%。

你可以看到这里有好几个变量,而且很容易在众多统计假设和公式中打结(需要指定样本变异、效应量、样本量、假设检验的 α 水平等,才能计算效能)。实际上,有专门的统计软件可以用来计算效能。大多数数据科学家并不需要完全按照发表论文那样严格地报告检验效能。不过,他们可能会遇到这样一种情况:要为一次 A/B 测试收集一些数据,而收集或处理数据是有成本的。这时,大致知道需要多少数据有助于避免“辛辛苦苦收集数据却得到无结论结果”的情况。下面是一种比较直观的替代方法:

  1. 用一些假设数据作为最初的“最佳猜测”,代表预期的实验结果(或许基于以往数据)——例如,用一个包含 20 个“1”和 80 个“0”的盒子来代表打击率 0.200 的打者,或用一个盒子存放一些“在网站停留时间”的观测值。
  2. 通过在第一个样本的基础上加上期望的效应量来创建第二个样本——例如,第二个盒子中包含 33 个“1”和 67 个“0”,或者在每条最初“网站停留时间”上加 25 秒。
  3. 从每个盒子中各抽取一个大小为 n 的自助(bootstrap)样本。
  4. 对这两个自助样本进行一次置换检验(或基于公式的假设检验),并记录它们之间的差异是否具有统计显著性。
  5. 多次重复前两步,并计算差异显著的比例——这就是估计的检验效能。

样本量

Sample Size

检验效能计算最常见的用途就是估算需要多大的样本量

例如,假设你在观察点击率(点击次数占展示次数的百分比),并测试一个新广告与现有广告。那在研究中你需要积累多少点击数呢? 如果你只对巨大的差异(比如 50% 的差异)感兴趣,相对较小的样本量可能就足够了;反之,如果你连很小的差异也想检测出来,那么就需要大得多的样本量。一种常见做法是先制定一个政策:新广告必须比现有广告好一定百分比,比如 10%;否则现有广告继续使用。这个目标——即“效应量”——就决定了所需的样本量。

例如,假设当前点击率约为 1.1%,你希望获得 10% 的提升,即 1.21%。那么我们有两个“盒子”:盒子 A 含有 1.1% 的“1”(比如 110 个“1”和 9,890 个“0”),盒子 B 含有 1.21% 的“1”(比如 121 个“1”和 9,879 个“0”)。起步时,我们从每个盒子各抽取 300 次(就像每个广告 300 次展示)。 假设第一次抽取得到如下结果:

  • 盒子 A:3 个“1”
  • 盒子 B:5 个“1”

显然,任何假设检验都能看出这种差异(5 比 3)完全在随机波动范围之内。这个样本量(每组 n=300)和效应量(10% 差异)的组合太小,无法让任何假设检验可靠地显示差异。

所以我们可以尝试增加样本量(比如 2,000 次展示),并要求更大的提升(50% 而不是 10%)。 例如,假设当前点击率仍为 1.1%,但我们现在寻求 50% 的提升至 1.65%。于是两个盒子:盒子 A 仍为 1.1% 的“1”(比如 110 个“1”和 9,890 个“0”),盒子 B 为 1.65% 的“1”(比如 165 个“1”和 9,868 个“0”)。 我们从每个盒子各抽取 2,000 次。假设第一次抽取得到:

  • 盒子 A:19 个“1”
  • 盒子 B:34 个“1”

对这个差异(34–19)进行显著性检验,仍显示为“非显著”(但比之前的 5–3 更接近显著)。要计算检验效能,我们需要多次重复前面的步骤,或使用可以计算效能的统计软件,但我们初步抽取就已经提示:即便是检测 50% 的提升,也需要几千次广告展示。

总结:在计算检验效能或所需样本量时,有四个“活动部件”:

  • 样本量
  • 希望检测的效应量
  • 检验的显著性水平(α)
  • 检验效能(Power)

指定其中三个,就可以计算出第四个。最常见的需求是计算样本量,因此必须指定其他三个条件。 在 R 和 Python 中,你还需要指定备择假设是“greater”还是“larger”以得到单尾检验;关于单尾与双尾检验的讨论见本书第 95 页“一尾与双尾假设检验”。

下面是一个 R 代码示例,计算两个比例的检验,其中两个样本大小相同(使用 pwr 包):

1
2
effect_size = ES.h(p1=0.0121, p2=0.011)
pwr.2p.test(h=effect_size, sig.level=0.05, power=0.8, alternative='greater')

输出:

1
2
3
4
5
6
Difference of proportion power calculation for binomial distribution (arcsine transformation)
h = 0.01029785
n = 116601.7
sig.level = 0.05
alternative = greater
NOTE: same sample sizes

函数 ES.h 用来计算效应量。我们看到,如果我们要 80% 的检验效能,需要的样本量接近 120,000 次展示。如果我们寻求的是 50% 的提升(p1=0.0165),所需样本量则减少到 5,500 次展示。

statsmodels 包中包含若干效能计算方法。这里用 proportion_effectsize 计算效应量,用 TTestIndPower 求样本量:

1
2
3
4
5
effect_size = sm.stats.proportion_effectsize(0.0121, 0.011)
analysis = sm.stats.TTestIndPower()
result = analysis.solve_power(effect_size=effect_size,
alpha=0.05, power=0.8, alternative='larger')
print('Sample Size: %.3f' % result)

输出:

1
Sample Size: 116602.393

关键要点

  • 确定所需样本量时,要提前考虑将要进行的统计检验。
  • 必须指定你希望检测的最小效应量。
  • 必须指定检测该效应量的概率(检验效能)。
  • 最后,必须指定检验的显著性水平(α)。

总结

实验设计的原则——将受试者随机分配到接受不同处理的两个或多个组——使我们能够对处理效果得出有效的结论。最好包含一个“无变化”的对照处理。正式统计推断这一主题——假设检验、p 值、t 检验及更多类似内容——在传统统计课程或教材中占据大量时间和篇幅,但从数据科学的角度来看,这种形式化在多数情况下并非必需。然而,认识到随机变动在误导人类大脑方面可能发挥的作用仍然很重要。直观的重抽样程序(置换和自助法)让数据科学家能够评估随机变动在其数据分析中可能起到的作用程度。