From ce0c096efdccee085bab8aacf58190b6bb28b754 Mon Sep 17 00:00:00 2001 From: wizardforcel <562826179@qq.com> Date: Wed, 24 Jan 2018 22:23:27 +0800 Subject: [PATCH] ch16. --- 16.md | 226 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 226 insertions(+) diff --git a/16.md b/16.md index d34f03c..e6677f8 100644 --- a/16.md +++ b/16.md @@ -612,3 +612,229 @@ Approximate 95% CI for the difference between means: 为了比较两个数值分布，将假设检验替换为估计，通常更富有信息。 只需估计一个差异，比如两组均值之间的差异。 这可以通过构建自举置信区间来完成。 如果零不在这个区间内，你可以得出这样的结论：这两个分布是不同的，你也可以估计均值有多么不同。 +## 因果 + +我们用于比较两个样本的方法在随机对照实验的分析中具有强大的用途。由于在这些实验中，实验组和对照组被随机分配，因此如果实验完全没有效果，结果中的任何差异，可以与仅仅由于分配中的随机性而发生的情况进行比较。如果观察到的差异比我们预测的，纯粹由于偶然的差异更为显着，我们就会有因果关系的证据。由于个体无偏分配到实验组和对照组，两组结果中的差异可归因于实验。 + +随机对照实验分析的关键是，了解偶然因素如何出现。这有助于我们设定明确的原假设和备选假设。一旦完成，我们可以简单地使用前一节的方法来完成分析。 + +让我们看看如何在一个例子中实现。 + +### 治疗慢性背痛：RCT + +成年人的背痛可能非常顽固，难以治疗。常见的方法从皮质类固醇到针灸。随机对照试验（RCT）检验了使用肉毒毒素 A 来治疗的效果。肉毒杆菌毒素是一种神经毒蛋白，会导致肉毒中毒的疾病；维基百科说肉毒杆菌是“已知最致命的毒素”。有七种类型的肉毒杆菌毒素。肉毒杆菌毒素 A 是可导致人类疾病的类型之一，但也被用于治疗涉及肌肉的各种疾病。福斯特（Foster），克拉普（Clapp）和贾巴里（Jabbari）在 2001 年分析的随机对照试验（RCT）将其用作一种治疗背痛的方法。 + +将 31 名背痛患者随机分为实验组和对照组，实验组 15 例，对照组 16 例。对照组给予生理盐水，试验是双盲的，医生和病人都不知道他们在哪个组。 + +研究开始 8 周后，实验组 15 名中的 9 名和对照组 16 名中的 2 名缓解了疼痛（由研究人员精确定义）。这些数据在`bta`表中，似乎表明实验有明显的益处。 + +```py +bta = Table.read_table('bta.csv') +bta +``` + +| Group | Result | +| --- | --- | +| Control | 1 | +| Control | 1 | +| Control | 0 | +| Control | 0 | +| Control | 0 | +| Control | 0 | +| Control | 0 | +| Control | 0 | +| Control | 0 | +| Control | 0 | + +（省略了 21 行） + +```py +bta.group('Group', np.mean) +``` + +| Group | Result mean | +| --- | --- | +| Control | 0.125 | +| Treatment | 0.6 | + +在实验组中，60% 的患者缓解了疼痛，而对照组只有 12.5%。没有一个患者有任何副作用。 + +因此，这表示是 A 型肉毒毒素比盐水更好。但结论还没确定。病人随机分配到两组，所以也许差异可能由于偶然？ + +为了理解这意味着什么，我们必须考虑这样的可能性，即在研究中的 31 个人中，有些人可能比其他人恢复得更好，即使没有任何治疗的帮助。如果他们中的大部分不正常地分配到实验组，只是偶然呢？即使实验组仅仅给予对照组的生理盐水，实验组的结果可能会好于对照组。 + +为了解释这种可能性，我们首先仔细建立机会模型。 + +### 潜在的结果 + +在患者随机分为两组之前，我们的大脑本能地想象出每个患者的两种可能的结果：患者分配到实验组的结果，以及分配给对照组的结果。这被称为患者的两个潜在的结果。 + +因此有 31 个潜在的实验结果和 31 个潜在的对照结果。问题关于 31 个结果的这两组的分布。他们是一样的，还是不一样？ + +我们还不能回答这个问题，因为我们没有看到每个组中的所有 31 个值。我们只能看到随机选择的 16 个潜在的对照结果，以及其余 15 个患者的实验结果。 + +这是一个展示设定的好方法。每个病人都有一张双面票： + +随机化之后，我们可以看到随机选择的一组票的右半部分，以及剩余分组的左半部分。 + +`observed_outcomes`表收集每个患者潜在结果的信息，每张“票”的未观察的一半留空。 （这只是考虑`bta`表的另一种方式，它载有的信息相同。） + +```py +observed_outcomes = Table.read_table("observed_outcomes.csv") +observed_outcomes.show() +``` + +| Group | Outcome if assigned treatment | Outcome if assigned control | +| --- | --- | --- | +| Control | Unknown | 1 | +| Control | Unknown | 1 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Control | Unknown | 0 | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 1 | Unknown | +| Treatment | 0 | Unknown | +| Treatment | 0 | Unknown | +| Treatment | 0 | Unknown | +| Treatment | 0 | Unknown | +| Treatment | 0 | Unknown | +| Treatment | 0 | Unknown | + +### 假设检验 + +问题是实验是否有用什么。根据观察得出的结果，问题在于第 2 列（包括未知数）的 31 个“实验”值的分布是否与第 3 列 31 个“对照”值的分布不同（同样包括未知数）。 + +原假设：所有 31 个潜在“实验”结果的分布与所有 31 个潜在“对照”结果的分布相同。实验与对照没有任何不同。两个样本的差异只是偶然而已。 + +备选假设：31 个潜在“实验”结果的分布与 31 个对照结果的分布不同。治疗做了一些不同于控制。 + +为了检验这些假设，请注意，如果原假设是真实的，那么 31 个观察结果的所有分布，对于标记为“对照”的 16 个结果和另外标记为“实验”的 15 个结果将具有相等的可能性。所以我们可以简单地对这些值进行排列，看看这两个组的分布是多么不同。更简单地说，由于数据是数值的，我们可以看到两个均值有多么不同。 + +这正是我们在上一节中为 A/B 测试所做的。样本 A 现在是对照组，样本 B 是实验组。我们的检验统计量是两组平均值的绝对差。 + +让我们为均值之间的差异运行我们的排列检验。只有 31 个观测值，所以我们可以运行大量的排列，而不必等待太久的结果。 + +```py +def permutation_test_means(table, variable, classes, repetitions): + + """Test whether two numerical samples + come from the same underlying distribution, + using the absolute difference between the means. + table: name of table containing the sample + variable: label of column containing the numerical variable + classes: label of column containing names of the two samples + repetitions: number of random permutations""" + + t = table.select(variable, classes) + + # Find the observed test statistic + means_table = t.group(classes, np.mean) + obs_stat = abs(means_table.column(1).item(0) - means_table.column(1).item(1)) + + # Assuming the null is true, randomly permute the variable + # and collect all the generated test statistics + stats = make_array() + for i in np.arange(repetitions): + shuffled_var = t.select(variable).sample(with_replacement=False).column(0) + shuffled = t.select(classes).with_column('Shuffled Variable', shuffled_var) + m_tbl = shuffled.group(classes, np.mean) + new_stat = abs(m_tbl.column(1).item(0) - m_tbl.column(1).item(1)) + stats = np.append(stats, new_stat) + + # Find the empirical P-value: + emp_p = np.count_nonzero(stats >= obs_stat)/repetitions + + # Draw the empirical histogram of the tvd's generated under the null, + # and compare with the value observed in the original sample + Table().with_column('Test Statistic', stats).hist() + plots.title('Empirical Distribution Under the Null') + print('Observed statistic:', obs_stat) + print('Empirical P-value:', emp_p) +permutation_test_means(bta, 'Result', 'Group', 20000) +Observed statistic: 0.475 +Empirical P-value: 0.00965 +``` + +经验 P 值非常小（研究报告 P 值为 0.009，这与我们的计算一致），因此证据倾向于备选假设：潜在实验和控制分布是不同的。 + +这是实验导致差异的证据，因为随机化确保了没有影响结论的混淆变量。 + +如果实验没有被随机分配，我们的测试仍然会指出我们 31 位患者的实验和背痛结果之间的关联。但要小心：没有随机化，这种关联并不意味着，实验会导致背痛结果的改变。例如，如果患者自己选择是否进行实验，则可能疼痛更严重的患者更可能选择实验，并且甚至在没有药物治疗的情况下，更可能减轻疼痛。预先存在的疼痛将成为分析中的混淆因素。 + +### 效果的置信区间 + +正如我们在上一节中指出的那样，只是简单地得出结论，说治疗是有用的，还不够。我们还想知道它做了什么。 + +因此，不要对两个基本分布的假设进行是与否的测试，而是仅仅估计它们之间的差异。具体来说，我们查看所有 31 个对照结果的平均值减去所有 31 个实验结果的平均值。这是未知的参数，因为我们只有 16 个对照值和 15 个实验值。 + +在我们的样本中，平均值的差异是 -47.5%。对照组平均为 12.5%，而治疗组平均为 60%。差异的负面信号表明实验组效果更好。 + +```py +group_means = bta.group('Group', np.mean) +group_means +``` + +| Group | Result mean | +| --- | --- | +| Control | 0.125 | +| Treatment | 0.6 | + +```py +group_means.column(1).item(0) - group_means.column(1).item(1) +-0.475 +``` + +但这只是一个样本的结果；样本可能会有所不同。 因此，我们将使用`bootstrap`复制样本，并重新计算差异。 这正是我们在前一节所做的。 + +```py +def bootstrap_ci_means(table, variable, classes, repetitions): + + """Bootstrap approximate 95% confidence interval + for the difference between the means of the two classes + in the population""" + + t = table.select(variable, classes) + + mean_diffs = make_array() + for i in np.arange(repetitions): + bootstrap_sample = t.sample() + m_tbl = bootstrap_sample.group(classes, np.mean) + new_stat = m_tbl.column(1).item(0) - m_tbl.column(1).item(1) + mean_diffs = np.append(mean_diffs, new_stat) + + left = percentile(2.5, mean_diffs) + right = percentile(97.5, mean_diffs) + + # Find the observed test statistic + means_table = t.group(classes, np.mean) + obs_stat = means_table.column(1).item(0) - means_table.column(1).item(1) + + Table().with_column('Difference Between Means', mean_diffs).hist(bins=20) + plots.plot(make_array(left, right), make_array(0, 0), color='yellow', lw=8) + print('Observed difference between means:', obs_stat) + print('Approximate 95% CI for the difference between means:') + print(left, 'to', right) +bootstrap_ci_means(bta, 'Result', 'Group', 20000) +Observed difference between means: -0.475 +Approximate 95% CI for the difference between means: +-0.759090909091 to -0.162393162393 +``` + -- GitLab