既然你已经了解了因果关系的基本知识,现在该讲讲因果推断中的 “推断” 部分了。本章首先会在随机实验的背景下,回顾上一章的一些概念。随机实验是因果推断的黄金标准,所以理解是什么让它们如此特殊非常重要。即便无法进行随机化,将其作为努力追求的理想目标,在思考因果关系时也会大有裨益。
接下来,我会借助随机实验回顾一些重要的统计概念和工具,比如误差、置信区间、假设检验、检验效能和样本量计算。如果你已经了解这些内容,我会明确标出回顾部分从哪里开始,这样你就可以跳过。
在上一章中,你了解了相关性与因果关系存在差异的原因及方式。你也了解了要让相关性等同于因果关系所需满足的条件:
\(E[Y \vert T = 1] - E[Y \vert T = 0] = E[Y_1 - Y_0 \vert T = 1] + \left\{ E[Y_0 \vert T = 1] - E[Y_0 \vert T = 0] \right\}\)
其中,\(E[Y_1 - Y_0 \vert T = 1]\) 为平均处理效应\((ATT, Average \space Treatment \space Effect \space on \space the \space Treated)\),\(\left\{ E[Y_0 \vert T = 1] - E[Y_0 \vert T = 0] \right\}\) 为偏差\((BIAS)\)。
回顾一下,相关性可被描述为两个部分的总和:对处理组的平均处理效应和偏差。只有当偏差部分为零时,测得的相关性才完全归因于因果关系。若\(E[Y_t \vert T = 0] = E[Y_t \vert T = 1]\),则不存在偏差。换句话说,若处理组和对照组除了所接受的处理外,是相同或可比的,那么相关性就是因果关系。或者,用更专业一点的术语来说,当处理组的潜在结果与未处理组的潜在结果至少在期望上相等时,就满足条件。要记住,潜在结果 \(Y_{ti}\) 是指若单元 \(i\) 接受处理\(t\),你会观测到的结果 。
在第 \(1\) 章中,我还简要提及了在潜在结果与处理相互独立的情况下,如何将相关性等同于因果关系,即:
\((Y_0, Y_1) \perp T\)
需要重点注意的是,我所说的并非处理和结果之间的独立性。要是那样的话,处理对结果就不会产生任何可供你衡量的影响。举个例子,假设处理是你公司应用程序中的一项新功能,结果是用户在该应用程序上花费的时间。若说 “功能\((Feature)\)与使用时长\((TimeSpent)\)相互独立”,那就意味着在处理组(使用新功能的用户)和未处理组(未使用新功能的用户)中,用户在应用程序上花费的时间是相同的。换句话说,这项新功能完全没有效果。
相反,你希望的是潜在结果与处理相互独立。这里有一个重要的区别。说 \(Y_1 \perp T\) 意味着,若受试者接受处理,原本会观测到的结果,与他们实际上是否真的接受了处理无关。类似地,\(Y_0 \perp T\) 表明,若受试者未接受处理,原本会观测到的结果,并不依赖于实际的处理分配情况。总而言之,实际观测到的结果 \(Y\) 仍然取决于实际分配的处理。
换一种更简单的说法,独立性假设意味着处理组和对照组具有可比性。或者说,知道处理分配情况并不会为我提供任何关于基准潜在结果 \(Y_0\) 的信息。因此,\((Y_0, Y_1) \perp T\) 意味着处理是导致处理组和对照组结果存在差异的唯一因素,即:
\(E[Y_0 \vert T = 0] = E[Y_0 \vert T = 1] = E[Y_0]\)
并且
\(E[Y_1 \vert T = 0] = E[Y_1 \vert T = 1] = E[Y_1]\)
如你所见,这使得通过对处理组和对照组的均值进行简单比较,就能识别出平均处理效应\(ATE\):
\(E[Y \vert T = 1] - E[Y \vert T = 0] = E[Y_1 - Y_0] = ATE\)
尽管独立性仅仅是一种假设,但如果你对处理 \(T\) 进行随机化,就能让这一假设变得更有说服力。通过随机化,你是将处理分配与抛硬币(一种我们完全知晓的随机机制)挂钩。这枚硬币不一定是公平的。你可以只将处理分配给 \(10\%\)、\(1\%\),甚至更少比例的受试者。只要分配机制是随机的,你就能获得识别处理效应所需的恰当条件。
通过对处理进行随机化,你能确保处理组和对照组大致(在期望上)具有可比性。它们之间唯一的系统性差异就是处理本身,这使得你可以将结果中的任何差异归因于该处理。从本质上讲,随机化以强制的方式让你朝着处理与潜在结果相互独立的方向推进。
现在,让我们结合这些数学知识,通过一个例子来演示,这样你会发现其实它相当简单。在下一部分内容中,我会使用随机对照试验(RCT)来理解交叉销售邮件的影响。
企业常用的一种策略是推出一款廉价甚至免费、本身不盈利的产品,将其作为吸引新客户的入口。一旦企业吸引到这些客户,就可以交叉销售其他利润更高的产品。假设你在一家咖啡配送公司工作。公司的主要产品是一项低成本的月度订阅服务,客户订阅后,每周能收到优质且经过精心挑选的咖啡。除了这项基础且低成本的订阅服务外,公司还提供更高端的订阅服务,包含冲泡福利以及世界上最优质的咖啡(比如来自巴西迪维诺兰迪亚小镇当地生产商的咖啡)。这是目前公司利润最高的服务,因此你的目标是向已订阅该低成本入门产品的用户推广,增加该高端服务的销量。为此,公司的营销团队尝试向客户推销高端咖啡配送订阅服务,主要通过交叉销售邮件来进行。作为因果推断专家,你的目标是了解这些邮件的效果如何。
当你查阅现有数据(非随机化数据)来回答这个问题时,你能清楚地看到,收到邮件的客户更有可能购买高端订阅服务。用专业术语来说,当客户购买了你试图推销的产品时,你可以说他们完成了转化。所以,你可以得出结论:收到邮件的客户转化率更高,即:
\(E[\text{Conversion} \vert \text{Email} = 1] > E[\text{Conversion} \vert \text{Email} = 0]\)
遗憾的是,你还发现,营销团队往往会优先给他们认为更有可能转化的客户发送邮件。目前还不完全清楚他们是如何判断的,可能是寻找与公司互动最多的客户,或者是在满意度调查中给出积极反馈的客户。不管怎样,这都强有力地证明了:
\(E[\text{Conversion}_0 \vert \text{Email} = 1] > E[\text{Conversion}_0 \vert \text{Email} = 0]\)
换句话说,即便根本没发送邮件,实际收到邮件的客户转化数量也会比其他客户更多。因此,对均值进行简单比较,会得到交叉销售邮件真实因果效应的有偏估计。要解决这个问题,你需要让处理组(收到邮件的客户 )和未处理组(未收到邮件的客户 )具有可比性,即满足 \(E[Y_0 \vert T = 1] = E[Y_0 \vert T = 0]\),而随机分配邮件发送可以实现这一点。要是你成功做到了,那么除了所接受的处理不同外,处理组和未处理组的平均转化率会是相同的。
那么,假设你真的这么做了。你从客户群中选取了三个随机样本。其中一个样本,你不发送任何邮件;另一个样本,你发送一封内容丰富、撰写精美的关于高端订阅服务的邮件;对于最后一个样本,你发送一封简短、直切主题的关于高端订阅服务的邮件。在收集了一段时间的数据后,你得到了如下类似的数据:
import pandas as pd # for data manipulation
import numpy as np # for numerical computation
# pd.set_option('display.max_rows', 5)
data = pd.read_csv("./data/cross_sell_email.csv")
data
## gender cross_sell_email age conversion
## 0 0 short 15 0
## 1 1 short 27 0
## 2 1 long 17 0
## 3 1 long 34 0
## 4 1 no_email 14 0
## .. ... ... ... ...
## 318 0 long 18 0
## 319 1 no_email 16 0
## 320 0 no_email 15 0
## 321 1 no_email 16 0
## 322 1 long 24 1
##
## [323 rows x 4 columns]
您可以看到,您有 \(323\) 个观测值。这算不上严格意义上的大数据,但也足以开展分析工作了。
模拟数据与真实世界数据
在教授因果推断相关内容时,使用模拟数据会非常有帮助。首先,因果推断往往伴随着对数据生成方式的说明。通过模拟,我能够确切无疑地讲解这种分配机制。其次,因果推断涉及反事实量,我可以有选择地展示这些反事实量,以便更好地阐释实际发生的情况。不过,为了避免数据显得过于虚假,我常常会选取真实世界的数据,并对其进行转换,使其契合我想要讲解的示例。例如,本示例所用数据源自威廉・T・阿尔珀特(William T. Alpert)等人 \(2016\) 年发表的论文《在线学习的随机评估》(A Randomized Assessment of Online Learning ),并经过转换,使其看起来像是交叉销售邮件数据。
要估算因果效应,您只需计算每个处理组的平均转化率即可:
data.groupby(["cross_sell_email"]).mean()
## gender age conversion
## cross_sell_email
## long 0.550459 21.752294 0.055046
## no_email 0.542553 20.489362 0.042553
## short 0.633333 20.991667 0.125000
没错,就是这么简单。可以看到,未发送邮件组的转化率为 \(4.2\%\),而发送长邮件组和短邮件组的转化率分别为 \(5.5\%\) 和高达 \(12.5\%\) 。因此,平均处理效应\((ATEs)\)通过对比每个处理组与对照组的差异来衡量,即\(ATE = E[Y|T = t] - E[Y|T = 0]\) ,长邮件和短邮件组对应的转化率分别提升了 \(1.3\) 和 \(8.3\) 个百分点 。有意思的是,发送简洁明了的短邮件似乎比内容详尽的长邮件效果更好。
随机对照试验\((RCTs)\)的妙处在于,您无需再担心营销团队是否有意瞄准了本就容易转化的客户,也不必担心不同处理组的客户在系统层面存在除所接受处理之外的其他差异 。从设计角度而言,随机试验就是为消除这些差异而设置的,理论上至少能实现 \((Y_0, Y_1) \perp T\) (即潜在结果与处理分配相互独立 )。
在实际操作中,要检验随机化是否合理(或者确认自己查看的是正确数据 ),一个不错的合理性检验方法是查看处理组和对照组在预处理变量上是否一致 。比如,您掌握了性别和年龄数据,就可以看看这些特征在不同处理组间是否均衡分布 。
当观察年龄这一变量时,各处理组(不同邮件策略组 )在年龄分布上似乎非常相似,但在性别(\(女性=1,男性=0\))分布上存在差异。具体而言,收到短邮件的组中男性占比达 \(63\%\) ,而对照组(未发邮件组 )男性占比为 \(54\%\),收到长邮件组的男性占比为 \(55\%\) 。
这一情况多少有些令人不安,因为您发现影响效果最显著的处理组(短邮件组 ),在性别分布上也与其他组不同。所以,即便理论上随机对照试验中(潜在结果与处理分配的)独立性应成立,但在实际操作中未必如此。有可能您观察到的短邮件带来的显著效果,是因为无论出于何种原因,存在 \(E [Y_0|man] > E [Y_0|woman]\)(即男性的潜在对照结果期望高于女性 )这一情况 。
对于如何评估(组间特征的)平衡性,目前尚未形成清晰的共识,但有一个非常简单的建议是:检查各处理组之间的标准化差异,计算公式如下:
\(\frac{\hat{\mu}_{tr} - \hat{\mu}_{co}}{\sqrt{(\hat{\sigma}_{tr}^{2} + \hat{\sigma}_{co}^{2}) / 2}}\)
其中,\(\hat{\mu}\) 、\(\hat{\sigma}^{2}\) 分别为样本均值和样本方差 。由于在你的示例中存在三个处理组,你只需针对对照组来计算这种差异即可 :
X = ["gender", "age"]
mu = data.groupby("cross_sell_email")[X].mean()
var = data.groupby("cross_sell_email")[X].var()
norm_diff = ((mu - mu.loc["no_email"])/ np.sqrt((var + var.loc["no_email"])/2))
norm_diff
## gender age
## cross_sell_email
## long 0.015802 0.221423
## no_email 0.000000 0.000000
## short 0.184341 0.087370
若此差异过小或过大,都应引起担忧。遗憾的是,对于差异达到何种程度算 “过大”,目前尚未有明确阈值,但 \(0.5\) 似乎是个实用的经验法则。在本示例中,虽未出现差异高达该数值的情况,但确实能发现:收到短邮件的组在性别方面存在较大差异,而收到长邮件的组在年龄方面存在较大差异 。
另见
若想深入探讨此主题,可查阅圭多・W・因本斯(Guido W. Imbens)与唐纳德・B・鲁宾(Donald B. Rubin)所著《因果推断:统计学、社会科学与生物医学科学导论》(Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction ,剑桥大学出版社出版 )一书的 14.2 节 。
要是此刻你觉得前面的公式有点像变魔术(让人摸不着头脑),别担心。等你看完本章的统计复习部分,它就会更清晰易懂了。眼下,我只想让你留意在小数据集里会出现的情况。即便经过随机化处理,也有可能因为偶然因素,某个组和另一个组存在差异。在大样本中,这种差异往往会消失。这也引出了一个问题,即多大的差异才足以让你断定处理(干预)确实有效,而非仅仅是偶然因素导致的,这是我很快会探讨的内容。
随机实验,亦即随机对照试验,是获取因果效应最为可靠的途径。这是一项简洁明了的技术,且具有极强的说服力。其效力之强大,致使大多数国家将其作为证明新药有效性的一项必备要求。不妨这般思考:倘若条件允许,随机对照试验会成为你为揭示因果关系所开展的唯一研究手段。一个设计精良的随机对照试验,是任何科学家与决策者的理想追求。
遗憾的是,这类试验往往要么成本极其高昂 —— 不仅耗费金钱,更关键的是耗费时间 —— 要么完全有悖伦理道德。在某些情况下,你根本无法对分配机制加以控制。设想你是一名医生,试图评估孕期吸烟对新生儿出生体重的影响。你无法简单地迫使随机选取的一部分孕妇在孕期吸烟。或者,假设你就职于一家大型银行,需要评估信用额度对客户流失的影响。为客户随机分配信用额度的成本会高到难以承受。再或者,你希望了解提高最低工资对失业率的影响。你无法随意给不同的国家指定不同的最低工资标准。此外,正如你将在第 \(3\) 章中看到的,在某些情形下(存在选择偏差的情形),即便采用随机对照试验也无法解决问题。
不过,我希望你能跳出 “随机实验仅作为揭示因果效应的工具” 这一思维。相反,此处的目标是将其用作一种基准。每当你在不借助随机对照试验进行因果推断时,都应始终问问自己:为回答你的问题,怎样的完美实验才是恰当的。即便那个理想实验不可行,它也能充当有价值的基准。它常常能为你指明方向,让你知晓即便没有这样的实验,该如何去发现因果效应 。
既然你已理解实验的价值,接下来就该审视 “没有无限数据” 意味着什么。因果推断是一个两步走的过程。随机对照试验在助力”识别”因果关系方面有着不可估量的价值,但要是实验的样本量较小,你在第二步(即 “推断” 环节)就会遇到难题。为理解这一点,回顾一些统计概念和工具是值得的。要是你已经熟悉这些内容,大可直接跳到下一章。
在 \(2007\)年的一篇著名文章中,霍华德・韦纳\((Howard \space Wainer)\)论述了一些极具危险性的方程:
“有些方程,知晓它们会带来危险,而另一些方程,不知晓它们才会带来危险。第一类方程可能带来危险,是因为其边界内的奥秘会开启通往可怕危险的大门。这类方程中最典型的当属爱因斯坦标志性的方程 \(E = MC^2\),因为它给出了普通物质中隐藏的巨大能量的度量方式…… 但我感兴趣的是另一类方程,它们的危险并非在我们了解它们时释放,而是在我们不了解它们时显现。这些方程若随时可用,能让我们清晰地理解事物,可一旦缺失,就会使我们陷入危险的无知状态。”
他所谈及的方程是德・莫伊弗\((Moivre)\)方程:
\(SE = \frac{\sigma}{\sqrt{n}}\)
其中,\(SE\) 是均值的标准误,\(\sigma\) 是标准差,\(n\) 是样本量。这个数学公式是你绝对应该掌握的内容,那咱们就开始学习吧。
要弄清楚为何不了解这个方程十分危险,让我们来看一些教育数据。我整理了不同学校\(3\)年内的 \(ENEM\) 成绩数据(巴西标准化高中学业成绩,类似于美国学术能力评估测试)。我也对数据进行了清理,只保留了本节中与你相关的信息。
如果你查看表现最优异的学校,会发现一个引人注意的点 —— 这些学校的学生数量都相当少:
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
from scipy import stats
import seaborn as sns
from matplotlib import pyplot as plt
from cycler import cycler
import matplotlib
default_cycler = (cycler(color=['0.1', '0.5', '1.0']))
color=['0.3', '0.5', '0.7', '0.9']
linestyle=['-', '--', ':', '-.']
marker=['o', 'v', 'd', 'p']
plt.rc('axes', prop_cycle=default_cycler)
matplotlib.rcParams.update({'font.size': 10})
df = pd.read_csv("data/enem_scores.csv")
df.sort_values(by="avg_score", ascending=False).head(10)
## year school_id number_of_students avg_score
## 16670 2007 33062633 68 82.97
## 16796 2007 33065403 172 82.04
## 16668 2005 33062633 59 81.89
## 16794 2005 33065403 177 81.66
## 10043 2007 29342880 43 80.32
## 18121 2007 33152314 14 79.82
## 16781 2007 33065250 80 79.67
## 3026 2007 22025740 144 79.52
## 14636 2007 31311723 222 79.41
## 17318 2007 33087679 210 79.38
换个角度看,你可以只筛选出排名前 \(1\%\) 的顶尖学校并对其展开研究。这些学校是什么样的呢?或许你能从最优秀的学校身上学到一些经验,并将其复制到其他地方。而事实也的确如此,如果你去看排名前 \(1\%\) 的学校,就会发现它们的平均学生数量更少:
plot_data = (df
.assign(top_school = df["avg_score"] >= np.quantile(df["avg_score"], .99))[["top_school", "number_of_students"]]
.query(f"number_of_students<{np.quantile(df['number_of_students'], .98)}")) # remove outliers
plt.figure(figsize=(8,4))
ax = sns.boxplot(x="top_school", y="number_of_students", data=plot_data)
plt.title("Number of Students of 1% Top Schools (Right)")
一个自然而然的结论是,小规模学校能带来更高的学业表现。这在直观上说得通,因为我们认为每个老师对应的学生数量越少,老师就越能专注关注每个学生。但这和德・莫伊弗的方程有什么关系呢?又为什么说它危险呢?
嗯,一旦人们开始依据这些信息做出重要且代价高昂的决策,危险就来了。在他的文章中,霍华德接着写道:
\(20\) 世纪 \(90\) 年代,倡导缩小学校规模开始流行起来。众多慈善组织和政府机构为拆分大型学校提供资金,因为小规模学校的学生在考试高分群体中占比过高。
人们忘了做的是,也去看看排名垫底的 \(1\%\) 的学校:它们的学生数量也非常少!
你在 \(图2-1\) 中看到的情况,完全符合德・莫伊弗方程的预期。随着学生数量增多,平均成绩会变得越来越精确。学生数量极少的学校(样本量小),仅因偶然因素,就可能出现极高或极低的成绩。这种情况在大规模学校中不太可能发生。德・莫伊弗方程阐述了一个与以数据形式呈现的信息和记录的实际情况相关的基本事实:数据信息始终是不精确的。接下来的问题就变成:不精确到什么程度?以及你能采取什么措施将这些不准确性考虑进去?
q_99 = np.quantile(df["avg_score"], .99)
q_01 = np.quantile(df["avg_score"], .01)
plot_data = (df
.sample(10000)
.assign(Group = lambda d: np.select([(d["avg_score"] > q_99) | (d["avg_score"] < q_01)], ["Top and Bottom"], "Middle")))
plt.figure(figsize=(10,5))
sns.scatterplot(y="avg_score", x="number_of_students", data=plot_data.query("Group=='Middle'"), label="Middle")
ax = sns.scatterplot(y="avg_score", x="number_of_students", data=plot_data.query("Group!='Middle'"), color="0.7", label="Top and Bottom")
plt.title("ENEM Score by Number of Students in the School")
量化我们不确定性的一种方法是利用估计值的方差。方差能告诉你观测值与其中心(期望)值的偏离程度。正如德・莫伊弗方程所表明的,这种不确定性会随着你观测的数据量增加而减小。这说得通,对吧?要是你看到一所学校里很多学生表现优异,你就能更有把握地认为这确实是一所好学校。然而,要是你看到一所只有 \(10\) 名学生的学校,其中 \(8\) 名学生表现良好,你就需要多持些怀疑态度。有可能只是碰巧这所学校招到了一些成绩高于平均水平的学生。
您在\(图2-1\)中看到的漂亮三角图 精准地讲述了这个情况。它向您展示,当样本量较小时,您对学校表现的估计值会有极大的方差。它也表明,随着样本量增加,方差会缩小。这一点对于学校的平均成绩而言是成立的,但对于您可能用到的任何描述性统计量(包括您常想要估计的\(ATE\))同样成立。回到我们的交叉销售邮件应用场景,如果每个处理组中有数千名客户,而非几百名,您就会更有把握认为,处理组和对照组之间看到的转化率差异不单纯是偶然因素导致的。
随机误差与系统误差
思考数据中这种不确定性的另一种方式,是对比系统误差和随机误差。系统误差是始终存在的偏差,会以相同方式影响所有测量结果;而随机误差是数据中因偶然因素产生的不可预测的波动。系统误差(或偏差)不会随着您收集更多数据而减小,因为它会把所有测量结果推向同一方向,偏离您想要估计的量。相反,如德・莫伊弗方程所示,随机误差会随着样本量的增加而减小。统计学是一门研究由随机误差导致的这些不精确性的科学,这样这些误差就不会让您猝不及防。它是一种将不确定性纳入考量的方法。
鉴于这只是统计学知识回顾,我就冒昧加快讲解节奏了。要是你不熟悉分布、方差和标准误这些概念,请继续阅读,但要记住你可能需要一些额外资料。我建议你在谷歌上搜索麻省理工学院的统计学入门课程,这类课程通常质量很高,而且你能在 \(YouTube\) 上免费观看。
在”最危险的方程”部分,你将平均处理效应 \(E[Y_1 - Y_0]\) 估计为处理组和未处理组均值的差值,即\(E[Y|T = 1] - E[Y|T = 0]\) 。具体而言,你算出了两种交叉销售邮件对转化率的平均处理效应 。随后你发现,短邮件带来了非常可观的提升,超过 \(8\) 个百分点,而长邮件的影响较小,仅使转化率提升 \(1.3\) 个百分点 。但仍有一个悬而未决的问题:这些效应是否足够大,大到你可以确信它们并非由偶然因素导致?用专业术语来说,你是否知道它们在统计上是否显著?
要解答这个问题,你首先需要依据我之前给出的方程估计标准误\((SE)\) 。样本量 \(n\) 很容易获取,你只需知道每个处理组的样本数量。或者,你可以使用 \(pandas\) 的 \(groupby\) 方法,接着进行 \(size\) 聚合操作,如下所示:
data = pd.read_csv("./data/cross_sell_email.csv")
short_email = data.query("cross_sell_email=='short'")["conversion"]
long_email = data.query("cross_sell_email=='long'")["conversion"]
email = data.query("cross_sell_email!='no_email'")["conversion"]
no_email = data.query("cross_sell_email=='no_email'")["conversion"]
data.groupby("cross_sell_email").size()
## cross_sell_email
## long 109
## no_email 94
## short 120
## dtype: int64
要得到标准差的估计值,你可以运用以下公式:
\(\hat{\sigma} = \sqrt{\frac{1}{N - 1} \sum_{i = 0}^{N} (x - \bar{x})^2}\)
其中,\(\bar{x}\) 是 \(x\) 的均值。
帽子符号
在本书中,我会用帽子符号(^)表示参数和预测结果的样本估计值 。
对你而言幸运的是,大多数编程软件已实现了这一功能。在 \(pandas\) 中,你可以使用 \(std\) 方法。把这些内容整合起来,你就能得到用于计算标准误的如下函数:
def se(y: pd.Series):
return y.std() / np.sqrt(len(y))
print("SE for Long Email:", se(long_email))
## SE for Long Email: 0.02194602460918551
print("SE for Short Email:", se(short_email))
## SE for Short Email: 0.030316953129541618
知晓这个公式极其有用(相信我,我们会多次用到它 ),不过要知道,\(pandas\) 也有一个内置方法用于计算标准误,即 \(.sem()\)(对应均值的标准误 ) 。
print("SE for Long Email:", long_email.sem())
## SE for Long Email: 0.02194602460918551
print("SE for Short Email:", short_email.sem())
## SE for Short Email: 0.030316953129541618
你的估计值的标准误是一种置信度量。若要精准理解其确切含义,你需涉足统计学中充满争议且复杂的领域。从频率主义统计学观点来看,我们会认为,我们的数据不过是潜在数据生成过程的一种体现。这一过程抽象且理想化,由恒定却对我们而言未知的真实参数所支配。就交叉销售邮件的情境而言,倘若你能开展多次实验并计算每次实验的转化率,这些转化率会围绕真实的潜在转化率分布,尽管不会与之完全等同。这与柏拉图有关 “理型论\((Forms)\)” 的论述极为相似:
每一种本质形式都会在各种各样的组合、行动、物质事物以及彼此之间展现自身,且每一种形式看似都呈现出多样的形态。
为理解这一点,假设你掌握了短交叉销售邮件转化率的真实抽象分布。由于转化率要么是 \(0\) 要么是 \(1\),它遵循伯努利分布,且假设该分布中的成功概率为 \(0.08\) 。也就是说,每当客户收到短邮件时,有 \(8\%\) 的概率发生转化。接下来,假设你能开展 \(10000\) 次实验。在每次实验中,你选取 \(100\) 名客户作为样本,给他们发送短邮件,然后观察平均转化率,这样你就能得到 \(10000\) 个转化率数据。这些实验得出的 \(10000\) 个转化率会围绕真实均值 \(0.08\) 分布(见\(图2-2\))。有些实验的转化率会低于真实值,有些则会高于真实值,但这 \(10000\) 个转化率的均值会非常接近真实均值。
列表推导式
每当我想要对序列中的每个元素应用一个函数时,我往往会大量使用列表推导式,而非 \(for\) 循环。列表推导式不过是映射型 \(for\) 循环的一种语法糖(即让代码更简洁、易读的语法形式 ) 。
import matplotlib.pyplot as plt
n = 100
conv_rate = 0.08
def run_experiment():
return np.random.binomial(1, conv_rate, size=n)
np.random.seed(42)
experiments = [run_experiment().mean() for _ in range(10000)]
plt.figure(figsize=(10,4))
freq, bins, img = plt.hist(experiments, bins=20, label="Experiment Means", color="0.6")
plt.vlines(conv_rate, ymin=0, ymax=freq.max(), linestyles="dashed", label="True Mean", color="0.5")
plt.legend()
这就是说,你永远无法确定实验的均值与真实的、柏拉图式的理想均值相匹配。不过,借助标准误,你能够构建一个区间,在你所开展的实验中,有 \(95\%\) 的实验里,这个区间会包含真实均值。
在现实生活中,你无法奢侈地用多个数据集模拟同一实验,通常你仅有一个数据集。但你可以借鉴模拟多次实验的思路来构建置信区间。置信区间带有与之关联的概率,最常见的是 \(95\%\)。该概率表明,如果你进行多次实验,并在每次实验中构建 \(95\%\) 置信区间,那么真实均值有 \(95\%\) 的概率会落在该区间内。
要计算置信区间,你会用到统计学中或许最令人惊叹的成果:中心极限定理。仔细看看你刚刚绘制的转化率分布。现在,记住转化率要么是 \(0\) 要么是 \(1\),因此遵循伯努利分布。要是你用直方图绘制这个伯努利分布,由于成功率仅为 \(8\%\),在 \(0\) 处会有一根很高的柱,在 \(1\) 处会有一根很矮的柱。这看起来一点也不像正态分布,对吧?
np.random.seed(42)
plt.figure(figsize=(10,4))
plt.hist(np.random.binomial(1, 0.08, 100), bins=20)
这就是那个令人惊叹的成果发挥作用之处了。尽管数据的分布并非正态分布(就像转化率的情况,其遵循伯努利分布 ),但数据的均值始终呈正态分布。要是你多次收集转化率数据,并且每次都计算平均转化率,这些均值会遵循正态分布。这太棒了,因为正态分布是众所周知的,你可以用它做各种各样有意思的事。比如,为了计算置信区间,你可以利用统计理论中的知识:正态分布 \(95\%\) 的概率质量落在均值上下 \(2\) 个标准差范围内(严格来说是 \(1.96\),但 \(2\) 是更容易记住的近似值 ),具体可参见\(图2-3\)。
exp_se = short_email.sem()
exp_mu = short_email.mean()
ci = (exp_mu - 2 * exp_se, exp_mu + 2 * exp_se)
print("95% CI for Short Email: ", ci)
## 95% CI for Short Email: (0.06436609374091676, 0.18563390625908324)
from scipy import stats
x = np.linspace(-4, 4, 100)
y = stats.norm.pdf(x, 0, 1)
plt.figure(figsize=(10,4))
plt.plot(x, y, linestyle="solid")
plt.fill_between(x.clip(-3, +3), 0, y, alpha=0.5, label="~99.7% mass", color="C2")
plt.fill_between(x.clip(-2, +2), 0, y, alpha=0.5, label="~95% mass", color="C1")
plt.fill_between(x.clip(-1, +1), 0, y, alpha=0.5, label="~68% mass", color="C0")
plt.ylabel("Density")
plt.legend()
回到你的交叉销售实验,现在你知道了,要是能开展多个类似实验,转化率会服从正态分布。对于那个(未知的)分布的均值,你所能得到的最佳估计值,就是你在小规模实验中得出的均值。此外,标准误可作为你对该样本均值未知分布的标准差的估计。所以,若你将标准误乘以 \(2\),再在实验均值的基础上加上和减去这个结果,就能构建出真实均值的\(95\%\)置信区间,步骤如下:
x = np.linspace(exp_mu - 4*exp_se, exp_mu + 4*exp_se, 100)
y = stats.norm.pdf(x, exp_mu, exp_se)
plt.figure(figsize=(10,4))
plt.plot(x, y, lw=3)
plt.vlines(ci[1], ymin=0, ymax=4, ls="dotted")
plt.vlines(ci[0], ymin=0, ymax=4, ls="dotted", label="95% CI")
plt.xlabel("Conversion")
plt.legend()
当然,你没必要局限于 \(95\%\) 置信区间。要是你想更谨慎些,也可以生成 \(99\%\) 置信区间。你只需用能涵盖正态分布 \(99\%\) 概率质量的系数去乘以标准差就行。
要找到这个系数,你可以使用 \(scipy\) 库中的 \(ppf\) 函数。该函数能给出标准正态分布累积分布函数\((CDF)\)的反函数。比如,\(ppf(0.5)\)会返回\(0.0\),这表明标准正态分布 \(50\%\) 的概率质量处于 \(0.0\) 以下。所以,对于任意显著性水平\(\alpha\),为得到\(1 - \alpha\)置信区间,你需要用以\(ppf((1 - \alpha)/2)\)绝对值为系数去乘以标准误\((SE)\) ,公式如下:
from scipy import stats
stats.norm.ppf((1-.99)/2)
## -2.5758293035489004
z = np.abs(stats.norm.ppf((1-.99)/2))
print(z)
## 2.5758293035489004
ci = (exp_mu - z * exp_se, exp_mu + z * exp_se)
ci
## (0.04690870373460816, 0.20309129626539185)
x = np.linspace(exp_mu - 4*exp_se, exp_mu + 4*exp_se, 100)
y = stats.norm.pdf(x, exp_mu, exp_se)
plt.figure(figsize=(10,4))
plt.plot(x, y, lw=3)
plt.vlines(ci[1], ymin=0, ymax=4, ls="dotted")
plt.vlines(ci[0], ymin=0, ymax=4, ls="dotted", label="99% CI")
ci_95 = (exp_mu - 1.96 * exp_se, exp_mu + 1.96 * exp_se)
plt.vlines(ci_95[1], ymin=0, ymax=4, ls="dashed")
plt.vlines(ci_95[0], ymin=0, ymax=4, ls="dashed", label="95% CI")
plt.xlabel("Conversion")
plt.legend()
这说的是短邮件的情况。你也可以给出与其他处理组相关的转化率的 \(95\%\) 置信区间,如下:
def ci(y: pd.Series):
return y.mean() - 2 * y.sem(), y.mean() + 2 * y.sem()
print("95% CI for Short Email:", ci(short_email))
## 95% CI for Short Email: (0.06436609374091676, 0.18563390625908324)
print("95% CI for Long Email:", ci(long_email))
## 95% CI for Long Email: (0.011153822341262012, 0.09893792077800405)
print("95% CI for No Email:", ci(no_email))
## 95% CI for No Email: (0.0006919679286838329, 0.08441441505003958)
import matplotlib
from cycler import cycler
default_cycler = (cycler(color=['0.1', '0.5', '1.0']))
color=['0.3', '0.5', '0.7', '0.9']
linestyle=['-', '--', ':', '-.']
marker=['o', 'v', 'd', 'p']
plt.rc('axes', prop_cycle=default_cycler)
matplotlib.rcParams.update({'font.size': 18})
plt.figure(figsize=(10,4))
x = np.linspace(-0.05, .25, 100)
short_dist = stats.norm.pdf(x, short_email.mean(), short_email.sem())
plt.plot(x, short_dist, lw=2, label="Short", linestyle=linestyle[0])
plt.fill_between(x.clip(ci(short_email)[0], ci(short_email)[1]), 0, short_dist, alpha=0.2, color="0.0")
long_dist = stats.norm.pdf(x, long_email.mean(), long_email.sem())
plt.plot(x, long_dist, lw=2, label="Long", linestyle=linestyle[1])
plt.fill_between(x.clip(ci(long_email)[0], ci(long_email)[1]), 0, long_dist, alpha=0.2, color="0.4")
no_email_dist = stats.norm.pdf(x, no_email.mean(), no_email.sem())
plt.plot(x, no_email_dist, lw=2, label="No email", linestyle=linestyle[2])
plt.fill_between(x.clip(ci(no_email)[0], ci(no_email)[1]), 0, no_email_dist, alpha=0.2, color="0.8")
plt.xlabel("Conversion")
plt.legend()
在此,你可以看到这\(3\)组的 \(95\%\) 置信区间相互重叠。要是它们不重叠,你就能得出结论:组间转化率的差异不单纯是偶然因素导致的。换句话说,你就能宣称,发送交叉销售邮件会使转化率产生统计上显著的差异。但由于这些区间确实存在重叠,你就无法这样断言,至少目前还不行。重要的是,置信区间重叠并不足以说明组间差异在统计上不显著;然而,要是它们不重叠,那就意味着这些差异在统计上是显著的。换句话说,不重叠的置信区间是统计显著性的保守性证据 。
概括而言,置信区间是一种为你的估计值赋予不确定性范围的方式。样本量越小,标准误就越大,相应地,置信区间也就越宽。鉴于置信区间计算极为简便,若缺少置信区间,要么意味着存在不良意图,要么单纯体现出知识欠缺,这两种情况都同样值得关注。最后,对于任何未给出不确定性度量的测量结果,你都应始终保持怀疑态度。
实际案例
新冠疫苗的有效性
随机对照试验\((Randomized \space control \space trials)\)对制药行业而言极为重要。鉴于新冠疫苗对地球上几乎所有人产生了巨大影响,或许最广为人知的例子就是为确定新冠疫苗有效性而开展的试验。以下是 \(2020\) 年发表的《Efficacy and Safety of the mRNA-1273 SARS-CoV-2 Vaccine 》研究中的结果部分:
该试验招募了 \(30420\) 名志愿者,以 \(1:1\) 的比例将他们随机分配,分别接受疫苗接种或注射安慰剂(每组 \(15210\) 名参与者 )。超过 \(96\%\) 的参与者完成了两剂接种,\(2.2\%\) 的参与者在基线时就有新冠病毒感染的证据(血清学、病毒学或两者兼有 )。安慰剂组中有 \(185\) 名参与者确诊出现有症状的新冠疾病(每千人年 \(56.5\) 例;\(95\%\) 置信区间为 \(48.7-65.3\) ),而 mRNA-1273 疫苗组中有 \(11\) 名参与者确诊(每千人年 \(3.3\) 例;\(95\%\) 置信区间为 \(1.7-6.0\) );疫苗效力为 \(94.1\%\)(\(95\%\) 置信区间为 \(89.3-96.8\) ;\(p\) 值 \(\leq0.001\) )。
结合你一直在学习的概念,以下是我对这些结果的解读。需注意,我并非健康领域专家,我的评论纯粹聚焦于统计与因果推断概念。
首先,他们明确了处理组(接种疫苗组 )和对照组(安慰剂组 ),说明处理(接种疫苗 )是随机分配的,这确保了处理与潜在结果相互独立。这样一来,他们就能从统计层面识别出疫苗的因果效应,具体就是识别出 \(E[Y|T = 0]\) 和 \(E[Y|T = 1]\) 这两个量(其中 \(T = 0\) 可表示对照组,\(T = 1\) 可表示处理组 )。接下来,他们把每千人年中出现有症状新冠疾病的情况定义为观测结果。最后,他们报告了对 \(E[Y|T = 0]\) 和 \(E[Y|T = 1]\) 估计值的 \(95\%\) 置信区间,分别为 \(48.7-65.3\) 以及 \(1.7-6.0\) 。这表明,与接种安慰剂的人群相比,接种疫苗人群中检测出有症状新冠疾病的情况要少得多。他们还报告了疫苗效力(即 \(E[Y|T = 0]/E[Y|T = 1]\) ),以及围绕该效力的 \(95\%\) 置信区间,为 \(89.3\%-96.8\%\)。
最后在此提醒一句。置信区间的解读,远比乍看之下要复杂。比如,我不该说某一个 \(95\%\) 置信区间有 \(95\%\) 的概率包含真实均值。在频率统计学中,总体均值被视为一个真实的总体常数。这个常数要么在某一特定置信区间内,要么在其外。换句话说,一个特定的置信区间要么包含真实均值,要么不包含。要是包含,包含的概率就是 \(100\%\),而非 \(95\%\);要是不包含,概率就是 \(0\%\)。相反,在置信区间的概念里,\(95\%\) 指的是:在众多研究中计算得出的这类置信区间,包含真实均值的频率。\(95\%\) 体现的是我们对用于计算 \(95\%\) 置信区间的算法的信心,而非对某个特定区间本身的信心 。
话虽如此,作为一名经济学家(统计学家们,现在请回避一下 ),我认为这种纯粹主义没什么太大用处。在实际情况中,你会看到人们说某个特定的置信区间有 \(95\%\) 的概率包含真实均值。虽然这种说法不对,但也没什么太大危害,因为它至少在你的估计结果中直观体现了不确定性程度。我的意思是,我宁愿你在估计值周围给出置信区间,哪怕解读有误,也不愿你因害怕误解而回避置信区间。我不在乎你是否说置信区间有 \(95\%\) 的概率包含真实均值。只是,求你了,永远不要忘记在你的估计值周围给出置信区间;否则,你会显得很外行。
可信区间
若你真心想给参数估计值落在某一区间内附上概率陈述,就应当了解贝叶斯可信区间 \((Bayesian \space credible \space intervals)\)。不过,以我的经验来看,在大多数情况下(尤其是样本量相对较大时 ),可信区间与频率主义置信区间的结果往往较为相似 。这也是我对置信区间被误读的情况往往更宽容的原因所在 。
纳入不确定性的另一种方式是进行假设检验:两组均值之差在统计上是否与零(或任何其他数值)存在差异?要回答这类问题,你需要回想一下:两个独立正态分布的和或差同样服从正态分布。所得分布的均值会是两个分布均值的和或差,而方差始终是两个分布方差的和,公式如下:
\(N(\mu_1, \sigma_1^2) - N(\mu_2, \sigma_2^2) = N(\mu_1 - \mu_2, \sigma_1^2 + \sigma_2^2)\)
\(N(\mu_1, \sigma_1^2) + N(\mu_2, \sigma_2^2) = N(\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2)\)
要是你不记得了,也没关系。你完全可以通过代码和模拟数据自行验证:
#%%
import seaborn as sns
from matplotlib import pyplot as plt
np.random.seed(123)
n1 = np.random.normal(4, 3, 30000)
n2 = np.random.normal(1, 4, 30000)
n_diff = n2 - n1
plt.figure(figsize=(10,4))
sns.distplot(n1, hist=False, label="$N(4,3^2)$", color="0.0", kde_kws={"linestyle":linestyle[0]})
sns.distplot(n2, hist=False, label="$N(1,4^2)$", color="0.4", kde_kws={"linestyle":linestyle[1]})
sns.distplot(n_diff, hist=False,label=f"$N(-3, 5^2) = N(1,4^2) - (4,3^2)$", color="0.8", kde_kws={"linestyle":linestyle[1]})
plt.legend();
如果你取两个组,每个组都对应一个分布,并用其中一个减去另一个,你最终会得到第三个分布。这个最终分布的均值将是两个均值的差,而标准差则是方差之和的平方根。
由于这里讨论的是实验平均值的分布,你可以把这些的标准差看作均值的标准误:
\[ \mu_{diff} = \mu_1 - \mu_2 \]
\[ SE_{diff} = \sqrt{SE_1^2 + SE_2^2} \]
你可以在比较交叉销售邮件实验转化率的问题中使用这个思路。如果你取两个组的估计分布——比如,简短邮件组和无邮件组——并用一个减去另一个,你就得到了差值的分布。利用这个分布,你可以很容易地构建均值差异的\(95\%\)置信区间:
exp_se = short_email.sem()
exp_mu = short_email.mean()
ci = (exp_mu - 2 * exp_se, exp_mu + 2 * exp_se)
print("95% CI for Short Email: ", ci)
## 95% CI for Short Email: (0.06436609374091676, 0.18563390625908324)
diff_mu = short_email.mean() - no_email.mean()
diff_se = np.sqrt(no_email.sem()**2 + short_email.sem()**2)
ci = (diff_mu - 1.96*diff_se, diff_mu + 1.96*diff_se)
print(f"95% CI for the differece (short email - no email):\n{ci}")
## 95% CI for the differece (short email - no email):
## (0.010239808474398426, 0.15465380854687816)
x = np.linspace(diff_mu - 4*diff_se, diff_mu + 4*diff_se, 100)
y = stats.norm.pdf(x, diff_mu, diff_se)
plt.figure(figsize=(10,3))
plt.plot(x, y, lw=3)
plt.vlines(ci[1], ymin=0, ymax=4, ls="dotted")
plt.vlines(ci[0], ymin=0, ymax=4, ls="dotted", label="95% CI")
plt.xlabel("Diff. in Conversion (Short - No Email)\n")
plt.legend()
plt.subplots_adjust(bottom=0.15)
有了这个区间,你就可以回答关于所谓原假设的问题。例如,你可以提出这样的假设:简短邮件和不发送邮件之间的转化率没有差异。你通常会用 \(H_0\) 来表示原假设:
\[ H_0 : \text{Conversion}_{\text{no email}} = \text{Conversion}_{\text{short email}} \]
一旦你有了这个假设,就该问自己:“如果原假设为真,我观察到如此大的差异的可能性有多大?”你会查看数据,判断其是否符合你的原假设。如果不符合,你就会说,如果原假设为真,出现这样的数据太奇怪了,因此你应该拒绝它。其中一种方法就是用你刚刚构建的置信区间来实现这一点。
注意前面的 \(95\%\) 置信区间不包含零。另外,回想一下,这是转化率差值的置信区间。由于原假设\((null \space hypothesis)\)表明这个差值为\(0\),但你可以看到置信区间完全处于\(0\)之外,那么可以说,如果原假设为真,观测到这样结果的概率会极低。因此,你可以以 \(95\%\) 的置信水平拒绝原假设。
显著性水平
显著性水平\(\alpha\),指的是在原假设为真时拒绝原假设的概率 —— 也就是犯第一类错误 \((Type \space 1 \space error )\)的概率 。显著性水平是在收集或分析数据之前就确定好的。要达到某一特定显著性水平(比如 \(5\%\)),在分析过程中,你需要围绕估计值构建一个置信水平为\(1 - \alpha\)(比如 \(95\%\))的置信区间 。
当然,除了 “完全没有差异” 这种原假设,你也可以设定其他原假设。举个例子,发送电子邮件是有成本的,这在现实中很常见。就算没有明显的货币成本,要是给客户发太多邮件,最终他们会把你标记为垃圾邮件发送者,这会切断与他们的沟通渠道,进而导致未来销量下降。在这种情况下,也许营销团队只有在转化率提升超过 \(1\%\) 时,才愿意推出交叉销售邮件。这时,你可以这样陈述原假设:“转化率差异为 \(1\%\)” 。要检验这个假设,你只需对置信区间进行调整,用均值差减去 \(1\%\) 即可 :
# shifting the CI
diff_mu_shifted = short_email.mean() - no_email.mean() - 0.01
diff_se = np.sqrt(no_email.sem()**2 + short_email.sem()**2)
ci = (diff_mu_shifted - 1.96*diff_se, diff_mu_shifted + 1.96*diff_se)
print(f"95% CI 1% difference between (short email - no email):\n{ci}")
## 95% CI 1% difference between (short email - no email):
## (0.00023980847439843134, 0.14465380854687815)
由于这个 \(95\%\) 置信区间也高于零,你同样可以拒绝这另一个原假设(\(\text{null hypothesis}\))。不过,现在这个 \(95\%\) 置信区间非常接近 \(0\),这意味着,对于 “效果等于 \(2\%\) 左右” 这样的原假设,你无法拒绝它,至少无法以 \(95\%\) 的置信水平拒绝 。
在本书中,大多数原假设会表述为等式(通常是等于 \(0\))。这类原假设的出发点是,仅当发现效应与 \(0\) 存在显著差异时才采取行动。然而,在某些情况下,你只想在处理效应等于 \(0\) 时采取行动。比如,考虑你想要叫停一场营销活动的情况。你只会在其效应可忽略不计(或不足以抵消其成本)时才这么做。在这些情形下,你需要将原假设表述为参数与某个值存在差异 。
这是因为,无法拒绝像 \(H_0 = 0\) 这样的原假设,和接受该假设为真不是一回事。这就对应那句著名的格言:“缺乏证据并非不存在证据\((absence \space of \space evidence \space is \space not \space evidence \space of \space absence)\)”。\(H_0 = 0\) 可能仅仅因为样本量过小(导致置信区间很宽 )就被拒绝,但这并不意味着该假设为真 。
为解决这一问题,统计学家创立了非劣效性检验。它是一种用于检验某种处理是否与另一种处理等效(或处理效应为\(0\))的方法。其基本思路是,查看置信区间是否包含\(0\),同时确保置信区间足够窄 。
除了置信区间,有时从检验统计量\((test \space statistic)\)的角度思考拒绝原假设也很有用。这些统计量的构造通常遵循 ”数值越大,越倾向于拒绝原假设” 的逻辑 。最常用的检验统计量之一是 \(t\) 统计量 。它可以通过对推导置信区间的分布进行标准化来定义:
\(t_\Delta = \frac{\mu_\Delta - H_0}{SE_\Delta} = \frac{(\mu_1 - \mu_2) - H_0}{\sqrt{\sigma_1^2/n_1 + \sigma_2^2/n_2}}\)
其中,\(H_0\) 是由你的原假设定义的数值 。
注意分子其实就是观测到的平均差值与原假设的差值。如果原假设为真,这个分子的期望值会是零:\(E[\mu_\Delta - H_0] = 0\) 。分母就是标准误,它对统计量进行标准化,使其方差为 \(1\) 。这样一来,若原假设为真,\(t_\Delta\) 会服从标准正态分布 ——\(N(0, 1)\) 。因为在原假设下,\(t_\Delta\) 以 \(0\) 为中心,所以数值大于 \(1.96\) 或小于 \(-1.96\) 的情况会极其罕见(出现概率低于 \(5\%\) )。这意味着,要是你观测到如此极端的 \(t\) 统计量,也能拒绝原假设 。在我们正在讨论的示例里,与 “无效应” 的 \(H_0\) 相关的统计量大于 \(2\),这表明你可以在 \(95\%\) 的置信水平下拒绝原假设 。
diff_mu = short_email.mean() - no_email.mean()
diff_se = np.sqrt(no_email.sem()**2 + short_email.sem()**2)
ci = (diff_mu - 1.96*diff_se, diff_mu + 1.96*diff_se)
print(f"95% CI for the differece (short email - no email):\n{ci}")
## 95% CI for the differece (short email - no email):
## (0.010239808474398426, 0.15465380854687816)
t_stat = (diff_mu - 0) / diff_se
t_stat
## 2.237951231871536
此外,由于在原假设下 \(t\) 统计量服从正态分布,你可以用它轻松计算 \(p\) 值 。
\(t\) 分布与正态分布
严格来讲,此处使用正态分布并不精确。实际上,你应该使用自由度等于样本量减去所估计参数数量(因为你在比较两个均值,所以是 \(2\) )的 \(t\) 分布 。不过,当样本量超过 \(100\) 时,这两种分布的差异在实际应用中就没什么重要意义了 。
之前提到过,要是未收到邮件和收到短邮件的客户转化率相同,观测到如此极端差异的概率小于 \(5\%\) 。但你能精确估算出这个概率到底是多少吗?观测到这样一个极端值的可能性有多大呢?这就该 \(p\) 值登场了!
和置信区间(实际上,和大多数频率统计学指标一样 )类似,\(p\) 值的准确定义可能非常令人困惑。所以,为了稳妥起见,我直接引用维基百科的定义:“\(p\) 值是在假定原假设正确的情况下,得到至少与实际观测到的测试结果一样极端的测试结果的概率”。
更简洁地说,\(p\) 值是在原假设为真时,观测到此类数据的概率(见图\(\boldsymbol{2 - 4}\) )。它衡量的是,在原假设为真的情况下,你所看到的测量结果有多不可能出现。当然,这常常会与 “原假设为真的概率” 相混淆。注意二者的区别:\(p\) 值不是\(P(H_0|\text{data})\) ,而是\(P(\text{data}|H_0)\) 。
x = np.linspace(-4, 4, 100)
y = stats.norm.pdf(x, 0, 1)
plt.figure(figsize=(10,4))
plt.plot(x, y, lw=2)
plt.vlines(t_stat, ymin=0, ymax=0.1, ls="dotted", label="T-Stat", lw=2)
plt.fill_between(x.clip(t_stat), 0, y, alpha=0.4, label="P-value")
plt.legend()
要计算 \(p\) 值,对于单侧原假设(“差值大于\(x\)” 或 “差值小于\(x\)” ),你只需计算标准正态分布中测试统计量左侧(或右侧,依情况而定 )的面积;对于双侧原假设(“差值为\(x\)” ),则将单侧计算结果乘以 \(2\),具体如下:
print("P-value:", (1 - stats.norm.cdf(t_stat))*2)
## P-value: 0.025224235562152142
\(p\) 值很有用,因为它让你无需指定像 \(95\%\) 或 \(99\%\) 这样的置信水平。不过,要是你想报告一个置信水平,通过 \(p\) 值,你能精准知晓你的检验在何种置信水平下会通过或不通过。比如,\(p\) 值为 \(0.025\) 时,在高达 \(2.5\%\) 的显著性水平下你都能得到显著结果。所以,虽然差值的 \(95\%\) 置信区间不包含零,但 \(99\%\) 置信区间会包含零。这个 \(p\) 值还意味着,要是差值真的为 \(0\),观测到这种极端检验统计量的概率只有 \(2.5\%\) 。
实际案例
面对面学习与在线学习
除病毒的直接影响外,\(2020\) 年疫情还带来了其他重要问题。其中最主要的是,孩子们无法上学,因此学习活动转入在线环境,且这种情况持续了长达两年之久。要估算这对一代人会产生何种影响颇具难度,因为从在线环境转回面对面学习的决策并非随机进行的。以巴西为例,公立学校恢复面授课程的时间就比私立学校要晚。
不过,尽管设计这样的实验绝非易事,但人们仍可设计实验来测试在线学习与面对面学习的影响,就像 Figlio、Rush 和 Yin 在《是实时授课还是网络授课?在线教学对学生学习效果的实验评估》(\(2013\) 年)中所做的那样。以下是该论文的摘要:
在一所顶尖研究型大学的大型微观经济学入门课程中,学生被随机分配去参加现场授课,或是在网络环境中观看相同的授课内容,且所有其他因素(如教学、补充材料等)都保持一致。与美国教育部近期针对高等教育网络教学的非实验性分析所做的元分析得出的结论相反,我们发现有适度证据表明,纯现场授课优于网络授课。这些结果对西班牙裔学生、男性学生以及成绩较差的学生而言尤为显著。我们还为其他情境下未来的实验提供了建议 。
注意,这项研究是在美国的一所大学开展的;很难说这些结果能推广到基础教育以及其他国家。用专业术语来讲,我们说该研究具有内部有效性,因为通过随机分组,处理组和对照组具有可比性。但就将结果推广到其他情境而言,这项研究可能不具备外部有效性,因为参与研究的对象并非来自总体的随机样本,而只是一所美国大学的经济学专业学生 。
到目前为止,你一直是从数据分析师的视角来看这些统计概念的,也就是拿到现有试验的数据,把数据当成既定事实来分析。但要是你被要求设计一项试验,而非只是解读已设计好的试验,那会怎样呢?这种情况下,你得确定每个变体所需的样本。比如,要是你还没开展交叉销售邮件试验,却需要决定给多少客户发长邮件、多少发短邮件,以及多少不发邮件,该怎么办?从这个角度看,目标就是要有足够大的样本,这样当 “无效应” 的原假设确实不成立时,你能正确拒绝它。一项检验正确拒绝原假设的概率,被称为该检验的效能\((power)\) 。这一概念不仅在你想确定试验所需样本量时很有用,在发现执行不佳的试验存在的问题时也很有帮助。
检验效能与统计显著性密切相关。\(\alpha\) 是原假设为真时拒绝原假设的概率,而检验效能(\(1 - \beta\) )是原假设为假时拒绝原假设的概率。从某种意义上讲,检验效能也可基于 \(\alpha\) 来定义,因为要正确拒绝原假设,你得明确拒绝所需的证据程度。
回想一下,\(95\%\) 置信区间意味着 \(95\%\) 的试验会包含你试图估计的真实参数。这也意味着 \(5\%\) 的试验不会包含,而这会导致你有 \(5\%\) 的概率错误地拒绝原假设。当 \(\alpha = 0.05\) 时,为了得出统计上显著的结论,参数估计值与原假设的差值 \(\delta\) 需至少距离零有 \(1.96SE\) 。这是因为 \(\delta - 1.96SE\) 是 \(95\%\) 置信区间的下限 。
好的,所以你需要\(\delta - 1.96SE > 0\)才能宣称结果显著。但你看到这种显著差异的可能性有多大呢?这时候就需要考虑检验效能了。检验效能是正确拒绝原假设的概率,即\(1 - \beta\),其中\(\beta\)是原假设为假时未拒绝原假设的概率(假阴性概率)。行业标准中,检验效能通常是 \(80\%\),也就是说,当原假设确实为假时,你只有 \(20\%\)(\(\beta = 0.2\))的概率不拒绝原假设。要达到 \(80\%\) 的检验效能,当原假设为假时,你需要有 \(80\%\) 的次数拒绝原假设。因为拒绝原假设意味着\(\delta - 1.96SE > 0\),所以你需要有 \(80\%\) 的次数得到这样大的差异。换句话说,你需要有 \(80\%\) 的次数让 \(95\%\) 置信区间的下限高于零 。
值得注意(或说也不算稀奇)的是,\(95\%\) 置信区间的下限同样服从正态分布。就像样本均值的分布那样,\(95\%\) 置信区间下限的分布,其方差等于标准误,但此时均值为\(\delta - 1.96SE\) 。它就是样本均值的分布,只是偏移了\(1.96SE\) 。因此,为了让\(\delta - 1.96SE > 0\)的情况在 \(80\%\) 的次数中出现(即 \(80\%\) 的检验效能 ),你需要让差值距离\(0\)达到\((1.96 + 0.84)SE\) :\(1.96\) 用于构建 \(95\%\) 置信区间,\(0.84\) 则用于让该区间的下限在 \(80\%\) 的次数中落在\(0\)以上 。
stats.norm.cdf(0.84)
## 0.7995458067395503
换个角度看,要意识到:若原假设为假,\(\delta\)(原假设与观测估计值的差值 )必须是可检测到的。当\(\alpha = 5\%\)且\(1 - \beta = 80\%\)时,可检测效应由\(2.8SE = (1.96SE + 0.84SE)\) 给出 。所以,要是你想设计一个交叉销售邮件实验,且想检测出 \(1\%\) 的差异,你必须保证样本量能让你至少得到\(1\% = 2.8SE\) 。若展开差值的 \(SE\) 公式,有\(SE_\Delta = \sqrt{SE_1^2 + SE_2^2}\) 。但要记住,此时你是从一名尚未开展实验、实际在尝试设计实验的分析师视角出发。这种情况下,你没有处理组的 \(SE\),但可假设处理组和对照组的方差相同,因此\(SE_\Delta = \sqrt{2SE^2} = \sqrt{2\sigma^2 / n} = \sigma\sqrt{2 / n}\) 。将其代入可检测差异中,若你希望有 \(80\%\) 检验效能和 \(95\%\) 显著性水平,最终会得到一个用于确定测试中每个变体样本量的相当简单的公式:
\(\begin{align*} \delta &= 2.8\sigma\sqrt{2/n}\\ n &= 2\times2.8^2\sigma^2/\delta^2 \approx 16\sigma^2/\delta^2 \end{align*}\)
其中,\(\delta\) 是可检测的差异,为保守起见,我把 \(2\times2.8^2\) 进行了四舍五入。将这个公式应用到你的数据上,把对照组的方差作为对 \(\sigma^2\) 的最佳估计,最终你会得到所需的样本量:
# in the book it is np.ceil(16 * no_email.std()**2/0.01), but it is missing the **2 in the denominator.
np.ceil(16 * (no_email.std()/0.08)**2)
## 103.0
data.groupby("cross_sell_email").size()
## cross_sell_email
## long 109
## no_email 94
## short 120
## dtype: int64
这在实验设计方面无疑极具价值,对我们当前开展的交叉销售实验而言也是好消息。在该实验中,两个处理组的样本量均超过 \(100\),对照组的样本量为 \(94\),这表明该实验具备恰当的检验效能 。
另见
这种计算样本量的简单方法源自 \(Ron \space Kohavi\) 等人 \(2022\) 年的论文《A/B 测试直觉误区:在线控制实验中的常见误解》(A/B Testing Intuition Busters: Common Misunderstandings in Online Controlled Experiments )。这个样本量公式只是论文中诸多有趣且实用内容之一,所以我强烈建议你去研读该论文 。
本章的主旨是将因果识别与估计相联系(同时回顾一些重要的统计概念 )。要记住,因果推断的目标是从数据中了解因果量。这一过程的第一步是识别,在此步骤中,你要利用关键假设,把不可观测的因果量转化为可从数据中估计的可观测统计量 。
例如,平均处理效应是一个因果量,它由不可观测的潜在结果定义:\(ATE = E[Y_1 - Y_0]\) 。为识别\((ATE)\),你要用到独立性假设\(T \perp (Y_0, Y_1)\) ,借助这一假设,你可以用可观测的量 \(E[Y|T = 1]\) 和 \(E[Y|T = 0]\) 来表示 \(ATE\) 。也就是说,在独立性假设下:
\(E[Y_1 - Y_0] = E[Y|T = 1] - E[Y|T = 0]\)
你还了解到,如何通过随机对照试验让这一假设更具合理性。要是你对处理进行随机分配,实际上就是强行让处理与潜在结果 \(Y_t\) 相互独立 。
但识别只是因果推断的第一步。一旦你能够用统计量来表述因果量,你仍需对那些统计量进行估计。比如,即便你可以用 \(E[Y|T = 1]\) 和 \(E[Y|T = 0]\)来表述平均处理效应,你还是得对它们进行估计。
本章第二部分涵盖了估计过程中用到的统计概念。具体而言,你了解了标准误:
\(SE = \sigma / \sqrt{n}\)
以及如何用它为估计值 \(\mu\) 构建置信区间:
\(\hat{\mu} \pm z^*SE\)
其中,\(z\) 是正态分布中包含 \(\alpha\%\) 概率质量的区间对应的临界值 。
你还学习了如何构建两组均值差的置信区间,其核心是将两组的方差相加,再求均值差的标准误:
\(SE_{diff} = \sqrt{SE_1^2 + SE_2^2}\)
最后,你了解了检验效能,以及如何用它计算你想要开展的实验所需的样本量。具体来说,对于 \(95\%\) 置信水平和 \(80\%\) 检验效能,你可以把样本量公式简化为:
\(N = 16*\sigma^2/\delta\)
其中,\(\sigma^2\) 是结果变量的方差,\(\delta\) 是可检测的差异 。