因果性
from datascience import *
%matplotlib inline
path_data = '../../../assets/data/'
import matplotlib.pyplot as plots
plots.style.use('fivethirtyeight')
import numpy as np我们比较两个样本的方法在分析随机对照实验中有着强大的用途。由于在这种实验中，治疗组和对照组是随机分配的，因此可以将它们结果的差异与如果治疗完全无效时纯粹由随机性造成的情况进行比较。如果观测到的差异比我们预期纯粹由随机性造成的差异更为显著，我们就有了“因果关系”的证据。由于个体被无偏地分配到治疗组和对照组，两组结果的差异可以归因于治疗。
分析随机对照实验的关键在于准确理解随机性是如何引入的。这有助于我们建立清晰的原假设和备择假设。一旦完成了这一步，我们就可以简单地使用前面章节的方法来完成分析。
让我们通过一个例子来看看如何操作。
治疗慢性背痛：一项随机对照试验
成人腰背痛可能非常顽固且难以治疗。常见方法从皮质类固醇到针灸，种类繁多。一项随机对照试验（RCT）研究了使用A型肉毒毒素（BTA）作为治疗手段的效果。肉毒毒素是一种引起肉毒杆菌病的神经毒性蛋白；维基百科称肉毒毒素是“已知最急性的致死毒素”。肉毒毒素有七种类型。A型肉毒毒素是能在人类中致病的类型之一，但它也被用于医学来治疗各种涉及肌肉的疾病。Foster、Clapp和Jabbari在2001年分析的这项随机对照试验将其作为腰背痛的治疗方法进行了研究。
31名腰背痛患者被随机分配到治疗组和对照组，其中治疗组15人，对照组16人。对照组接受生理盐水治疗，并且试验采用双盲进行，医生和患者都不知道自己在哪一组。
研究开始八周后，治疗组15人中有9人、对照组16人中有2人疼痛得到缓解（根据研究人员使用的精确定义）。这些数据在表格 bta 中，似乎表明治疗有明显的益处。
bta = Table.read_table(path_data + 'bta.csv')
bta.show()<IPython.core.display.HTML object>让我们看看每组有多少患者康复。请记住，计数等同于将0和1相加。对照组中1的总和就是对照组中缓解疼痛的患者人数。
bta.group('Group', sum)Group | Result sum
Control | 2
Treatment | 9由于计数等同于0和1的“总和”，因此缓解疼痛的患者“比例”就是0和1的“平均值”。它是总和除以每组患者总数。
bta.group('Group', np.average)Group | Result average
Control | 0.125
Treatment | 0.6在治疗组中，60%的患者疼痛缓解，而对照组仅为12.5%。没有患者出现任何副作用。
因此，迹象表明A型肉毒毒素比生理盐水效果更好。但这个结论还不是板上钉钉。患者是被随机分配到两组的，所以这种差异可能只是由随机性造成的？
要理解这意味着什么，我们必须考虑以下可能性：在研究的31名患者中，有些人即使没有治疗的帮助，也比其他人更容易康复。如果这些患者中异常大的比例碰巧被分配到了治疗组呢？那么，即使治疗的效果不比对照组的生理盐水更好，治疗组的结果看起来也可能比对照组好。
为了解释这种可能性，让我们先仔细设定随机模型。
潜在结果
在患者被随机分配到两组之前，我们本能地想象每个患者有两种可能的结果：患者被分配到治疗组时的结果，以及同一患者被分配到对照组时的结果。这被称为患者的两个“潜在结果”（potential outcomes）。
因此，有31个潜在治疗结果和31个潜在对照结果。问题在于这两组各31个结果的分布是否相同。
我们还不能回答这个问题，因为我们无法看到每组中的所有31个值。我们只能看到随机选择的16个潜在对照结果，以及其余15名患者的治疗结果。
这是一个很好可视化该场景的方法。每个患者都有一张双面票：
随机化之后，我们能看到随机选择的一组票的右半部分，以及剩余组的左半部分。
表格 observed_outcomes 收集了每个患者潜在结果的信息，将每张“票”未被观察到的一半留为空白。（这只是另一种看待 bta 表的方式，携带相同的信息。）
observed_outcomes = Table.read_table(path_data + "observed_outcomes.csv")
observed_outcomes.show()<IPython.core.display.HTML object>提出假设
问题是治疗是否有效。用表格 observed_outcomes 来说，问题是第1列中31个“治疗”值（包括未知值）的分布是否与第2列中31个“对照”值（同样包括未知值）的分布不同。
原假设： 所有31个潜在“治疗”结果的分布与所有31个潜在“对照”结果的分布相同。A型肉毒毒素与生理盐水没有区别；两个样本之间的差异仅由随机性造成。
备择假设： 31个潜在“治疗”结果的分布与31个对照结果的分布不同。治疗的效果与对照不同。
请注意，备择假设并未具体说明治疗是否有帮助——只是说它与对照不同。这在医学研究中是标准做法，因为它不预先判断结果可能的方向。但你也可以进行检验来判断治疗是否优于对照，只需相应调整你的检验统计量即可。
两组共有31个观测结果。如果原假设成立，那么这31个结果中哪些被标记为“治疗”、哪些被标记为“对照”就无关紧要了。31个值中任意16个的随机子集都可以被称为“对照”，剩余的15个称为“治疗”。
我们可以模拟这一点。我们可以随机置换这31个值，将它们分成16和15的两组，然后观察两个组的分布有多大差异。由于数据是0和1，我们只需看两个比例有多大差异即可。
这正是我们在上一节A/B测试中所做的。现在样本A是对照组，样本B是治疗组。我们将在下面进行检验，展示所有步骤的细节。你应该确认这些步骤与之前进行的步骤相同。...（内容已截断）
bta.group('Group', np.average)Group | Result average
Control | 0.125
Treatment | 0.6observed_proportions = bta.group('Group', np.average).column(1)
observed_distance = abs(observed_proportions.item(0) - observed_proportions.item(1))
observed_distance0.475和之前一样，我们将定义一个函数，接受以下两个参数：
- 数据表的名称
- 组标签的列标签
并返回两组比例之间的距离。
def distance(table, group_label):
reduced = table.select('Result', group_label)
proportions = reduced.group(group_label, np.average).column(1)
return abs(proportions.item(1) - proportions.item(0))distance(bta, 'Group')0.475在原假设下预测统计量
我们可以在原假设下模拟结果，以观察如果原假设成立，我们的检验统计量会是什么样子。
生成统计量的一个值
模拟过程与我们上一节使用的完全相同。我们首先随机置换所有组标签，然后将打乱后的标签附加到0/1结果上。
shuffled_labels = bta.sample(with_replacement=False).column(0)bta_with_shuffled_labels = bta.with_column('Shuffled Label', shuffled_labels)
bta_with_shuffled_labels.show()<IPython.core.display.HTML object>现在我们可以找出组标签被打乱后两个比例之间的距离。
distance(bta_with_shuffled_labels, 'Shuffled Label')0.08750000000000002这与两个原始比例之间的距离有很大不同。
distance(bta_with_shuffled_labels, 'Group')0.475置换检验
如果我们再次打乱标签，新的距离会有多大差异？为了回答这个问题，我们将定义一个函数，模拟在来自相同潜在分布的随机抽取假设下的一个距离模拟值。然后我们将在数组中收集20,000个这样的模拟值。
可以看到，我们正在做的就是我们在之前的置换检验示例中所做的。
def one_simulated_distance():
shuffled_labels = bta.sample(with_replacement = False
).column('Group')
shuffled_table = bta.select('Result').with_column(
'Shuffled Label', shuffled_labels)
return distance(shuffled_table, 'Shuffled Label') distances = make_array()
repetitions = 20000
for i in np.arange(repetitions):
new_distance = one_simulated_distance()
distances = np.append(distances, new_distance)检验的结论
数组 distances 包含在原假设下模拟的20,000个检验统计量值。以下是它们的经验直方图以及统计量的观测值。要找到检验的p值，请记住距离的大值支持备择假设。
Table().with_column('Distance', distances).hist(
bins = np.arange(0, 0.7, 0.1), left_end = observed_distance)
# Plotting parameters; you can ignore the code below
plots.ylim(-0.1, 5.5)
plots.scatter(observed_distance, 0, color='red', s=40, zorder=3)
plots.title('Prediction Under the Null Hypothesis')
print('Observed Distance', observed_distance)Observed Distance 0.475
Histogram titled 'Prediction Under the Null Hypothesis' with 'Distance' on the x-axis and 'Percent per unit' on the y-axis. The data ranges from 0 to 0.6, with the tallest bar between 0 and 0.1, the next two bars from 0.1 to 0.2 and 0.2 to 0.3 are a little under half the height of the first bar. The remaining bars are quite small. There is a red dot between 0.4 and 0.5 and the data on the right of the red dot is shaded gold. The gold shaded region is quite small.为了数值化地找到经验p值，我们必须找出等于或大于观测统计量的模拟统计量所占的比例。
empirical_p = np.count_nonzero(distances >= observed_distance) / repetitions
empirical_p0.00875这是一个很小的p值。观测统计量位于原假设下生成的检验统计量经验直方图的尾部。
结果是统计显著的。检验支持备择假设而非原假设。证据支持治疗有效的假设。
研究报告的p值为0.009（即0.9%），与我们的经验值相差不远。
因果性
由于试验是随机化的，该检验提供了治疗“导致”差异的证据。患者被随机分配到两组确保了没有混杂变量影响因果性的结论。
如果治疗不是随机分配的，我们的检验仍然会指向治疗与31名患者的背痛结果之间的“关联”。但请注意：没有随机化，这种关联并不意味着治疗导致了背痛结果的变化。例如，如果患者自己选择是否接受治疗，那么疼痛更剧烈的患者可能更倾向于选择治疗，并且即使没有药物也更容易经历一定程度的疼痛缓解。已有的疼痛就会成为分析中的一个“混杂因素”（confounding factor）。
荟萃分析
虽然这项随机对照试验确实提供了A型肉毒毒素治疗对患者有帮助的证据，但31名患者的研究不足以确立一种医疗方法的有效性。这不仅仅是因为样本量小。我们在本节中的结果对研究中的31名患者是有效的，但我们真正感兴趣的是所有可能患者的总体。
2011年，一组研究人员对这些治疗研究进行了一项荟萃分析（meta-analysis）。也就是说，他们识别了所有关于此腰背痛治疗的现有研究，并总结了整理后的结果。
虽然有若干项研究，但能够以科学严谨的方式纳入的并不多：“我们因非随机化、数据不完整或未发表而排除了19项研究的证据。”仅剩下三项随机对照试验，其中之一就是我们在本节中研究的这项。荟萃分析在所有研究中给予了它最高的评价（LBP指腰背痛）：“我们确定有三项研究考察了BoNT对LBP的疗效，但只有一项偏倚风险较低且评估了非特异性LBP患者（N = 31）。”
综合所有结果，荟萃分析得出结论：“有低质量证据表明BoNT注射在改善疼痛、功能或两者方面优于盐水注射，而有极低质量证据表明它们优于针灸或类固醇注射。……进一步的研究很可能对效果的估计及其可信度产生重要影响。未来的试验应规范患者群体、治疗方案和对照组，招募更多参与者，并纳入长期结果、成本效益分析和发现的临床相关性。”
确立一种医疗方法具有有益效果需要大量细致的工作。知道如何分析随机对照试验...（内容已截断）