# 十六、比較兩個樣本 > 原文：[Comparing Two Samples](https://github.com/data-8/textbook/tree/gh-pages/chapters/16) > 譯者：[飛龍](https://github.com/wizardforcel) > 協議：[CC BY-NC-SA 4.0](http://creativecommons.org/licenses/by-nc-sa/4.0/) > 自豪地采用[谷歌翻譯](https://translate.google.cn/) 最近鄰分類方法的動機是這樣的，個體可能像最近的鄰居。 從另一個角度來看，我們可以說一個類別的個體不像另一個類別中的個體。 機器學習為我們提供了一種有力的方法來發現這種相似性的缺乏，并將其用于分類。 它揭示了一種模式，通過一次檢查一兩個屬性，我們不一定能發現它。 但是，我們可以從屬性中學到很多東西。 為了了解它，我們將比較兩個類中的屬性分布。 讓我們來看看 Brittany Wenger 的乳腺癌數據，看看是否只用一個屬性，就有希望生成一個合理的分類器。 和以前一樣，我們將在隨機選擇的訓練集上進行探索，然后在剩余的保留集上測試我們的分類器。 ```py patients = Table.read_table('breast-cancer.csv').drop('ID') shuffled_patients = patients.sample(with_replacement=False) training_set = shuffled_patients.take(np.arange(341)) test_set = shuffled_patients.take(np.arange(341, 683)) training_set ``` | Clump Thickness | Uniformity of Cell Size | Uniformity of Cell Shape | Marginal Adhesion | Single Epithelial Cell Size | Bare Nuclei | Bland Chromatin | Normal Nucleoli | Mitoses | Class | | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | | 5 | 1 | 1 | 1 | 2 | 1 | 2 | 1 | 1 | 0 | | 5 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | | 4 | 1 | 1 | 1 | 2 | 1 | 1 | 1 | 1 | 0 | | 5 | 1 | 2 | 1 | 2 | 1 | 3 | 1 | 1 | 0 | | 4 | 10 | 8 | 5 | 4 | 1 | 10 | 1 | 1 | 1 | | 7 | 2 | 4 | 1 | 3 | 4 | 3 | 3 | 1 | 1 | | 9 | 4 | 5 | 10 | 6 | 10 | 4 | 8 | 1 | 1 | | 3 | 1 | 1 | 1 | 2 | 2 | 3 | 1 | 1 | 0 | | 3 | 2 | 1 | 1 | 2 | 1 | 2 | 2 | 1 | 0 | | 6 | 3 | 3 | 5 | 3 | 10 | 3 | 5 | 3 | 0 | （省略了 331 行） 讓我們看看第二個屬性`Uniformity of Cell Size`，能告訴我們患者分類的什么事情。 ```py training_cellsize = training_set.select('Class', 'Uniformity of Cell Size').relabel(1, 'Uniformity') training_cellsize ``` | Class | Uniformity | | --- | --- | | 0 | 1 | | 0 | 1 | | 0 | 1 | | 0 | 1 | | 1 | 10 | | 1 | 2 | | 1 | 4 | | 0 | 1 | | 0 | 2 | | 0 | 3 | （省略了 331 行） `Class`和`Uniformity`列顯示為數字，但他們真的都是類別值。 這些類別是“癌癥”（1）和“非癌癥”（0）。 `Uniformity`為 1-10，但是這些標簽是由人確定的，他們也可能有十個標簽，如“非常一致”，“不一致”等等。 （一致性的 2 不一定是 1 的兩倍。）所以我們比較兩個類別分布，每個分類一個。 對于每個類別和每個一致評分，我們都需要訓練集的患者數量。`pivot`方法將為我們計數。 ```py training_counts = training_cellsize.pivot('Class', 'Uniformity') training_counts ``` | Uniformity | 0 | 1 | | --- | --- | --- | | 1 | 181 | 3 | | 2 | 21 | 2 | | 3 | 16 | 15 | | 4 | 4 | 18 | | 5 | 0 | 17 | | 6 | 0 | 8 | | 7 | 0 | 8 | | 8 | 1 | 13 | | 9 | 1 | 4 | | 10 | 0 | 29 | 我們現在有了一些東西，類似于每個類別的一致評分的分布。 而這兩者看起來相當不同。 但是，我們要小心 - 這兩個類別的患者總數是 341（訓練集的大小），超過一半的人在類別 0 里面。 ```py np.sum(training_counts.column('0')) 224 ``` 所以為了比較兩個分布，我們應該把計數轉換成比例然后可視化。 ```py def proportions(array): return array/np.sum(array) training_dists = training_counts.select(0).with_columns( '0', proportions(training_counts.column('0')), '1', proportions(training_counts.column('1')) ) training_dists.barh('Uniformity') ```  這兩個分布看起來不一樣！ 事實上，它們看起來相當不同，我們應該能夠基于對這種差異的直截了當的觀察來構建一個非常合理的分類器。 一個簡單的分類規則是：“如果一致性大于 3，類別就是 1，也就是說這個單元格就有癌癥的，否則類別就是 0。 這么粗糙的東西有什么好處嗎？ 讓我們試試看。 對于測試集中的任何個體，我們所要做的就是，查看一致評分是否大于 3。例如，對于前 4 名患者，我們將得到一組四個布爾值： ```py test_set.take(np.arange(4)).column('Uniformity of Cell Size') > 3 array([ True, False, False, False], dtype=bool) ``` 請記住，`True`等于`1`，如果一致性大于 3，那么這是我們要劃分的分類。因此，為了測量粗分類器的準確性，我們所要做的就是，求得測試集患者的比例， 其中分類與患者已知的分類相同。 我們將使用上一節中寫的`count_equal`函數。 ```py classification = test_set.column('Uniformity of Cell Size') > 3 count_equal(classification, test_set.column('Class'))/test_set.num_rows 0.935672514619883 ``` 這相當準確，即使我們只使用單個屬性單行代碼的分類器！ 這是否意味著上一章中最近鄰的方法是不必要的？ 不，因為那些更準確，并且對于癌癥診斷，任何患者都想要盡可能精確的方法。 但是看到簡單的方法并不壞，這是令人欣慰的。 ## 兩個類別分布 為了查看兩個數值變量如何相關，可以使用相關系數來衡量線性關聯。 但是，我們應該如何確定兩個分類變量是否相關？ 例如，我們如何決定一個屬性是否與個體的類別有關？ 這是一個很重要的問題，因為如果不相關的話，你可以把它從你的分類器中刪除。 在乳腺癌數據中，我們來看看有絲分裂活動是否與這個類別有關。 我們已經標記了“癌癥”和“非癌癥”的類別，以便以后參考。 ```py classes = Table().with_columns( 'Class', make_array(0, 1), 'Class Label', make_array('Not Cancer', 'Cancer') ) patients = Table.read_table('breast-cancer.csv').drop('ID').join('Class', classes) patients = patients.drop('Class').relabel('Class Label', 'Class') mitoses = patients.select('Class', 'Mitoses') mitoses ``` | Class | Mitoses | | --- | --- | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 5 | | Not Cancer | 1 | | Not Cancer | 1 | （省略了 673 行） 我們可以使用`pivot`和`proportions`（在前面的章節中定義）來顯示兩類中`Mitoses`的分布。 ```py counts = mitoses.pivot('Class', 'Mitoses') counts ``` | Mitoses | Cancer | Not Cancer | | --- | --- | --- | | 1 | 132 | 431 | | 2 | 27 | 8 | | 3 | 31 | 2 | | 4 | 12 | 0 | | 5 | 5 | 1 | | 6 | 3 | 0 | | 7 | 8 | 1 | | 8 | 7 | 1 | | 10 | 14 | 0 | ```py dists = counts.select(0).with_columns( 'Cancer', proportions(counts.column(1)), 'Not Cancer', proportions(counts.column(2)) ) dists.barh(0) ```  與“非癌癥”類別的分布相比，“癌癥”類別的`Mitoses`都集中于最低評分。 所以看起來類別和有絲分裂活動是相關的。 但是，這可能只是由于偶然嘛？ 為了了解偶然來自哪里，請記住，數據就像是來自更大總體的隨機樣本 - 總體包含我們可能要分類的新個體。 可能在總體中，類別和有絲分裂是相互獨立的，只是由于偶然與樣本相關。 ### 假設 我們試著通過對以下假設進行測試來回答這個問題。 原假設。 在總體中，類別和有絲分裂評分是相互獨立的；換句話說，這兩個類別的有絲分裂的分布是一樣的。 由于偶然性，樣本分布是不同的。 備選假說。 在總體中，類別和有絲分裂評分是相關的。 為了了解如何測試它，我們再看一下數據。 ```py mitoses ``` | Class | Mitoses | | --- | --- | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 1 | | Not Cancer | 5 | | Not Cancer | 1 | | Not Cancer | 1 | （省略了 673 行） ### 隨機排列 如果類別和有絲分裂評分是不相關的，那么`Mitoses`值出現的順序并不重要，因為它們與類別的值無關，所有的重新排列應該是等可能的。 這與我們在分析足球`Deflategate`數據時采用的方法相同。 所以讓我們將所有的`Mitoses`值整理到一個名為`shuffled_mitoses`的數組中。 你可以看到下面的第一項，但它包含 683 個項目，因為它是整個`Mitoses`列的排列（即重新排列）。 ```py shuffled_mitoses = mitoses.select('Mitoses').sample(with_replacement=False).column(0) shuffled_mitoses.item(0) 1 ``` 讓我們擴展`mitoses`表，添加一列亂序的值。 ```py mitoses = mitoses.with_column('Shuffled Mitoses', shuffled_mitoses) mitoses ``` | Class | Mitoses | Shuffled Mitoses | | --- | --- | --- | | Not Cancer | 1 | 1 | | Not Cancer | 1 | 1 | | Not Cancer | 1 | 1 | | Not Cancer | 1 | 1 | | Not Cancer | 1 | 7 | | Not Cancer | 1 | 1 | | Not Cancer | 1 | 1 | | Not Cancer | 5 | 3 | | Not Cancer | 1 | 1 | | Not Cancer | 1 | 2 | （省略了 673 行） 讓我們看看亂序數據的有絲分裂的分布，使用與原始數據相同的過程。 ```py shuffled = mitoses.select('Class', 'Shuffled Mitoses') shuffled_counts = shuffled.pivot('Class', 'Shuffled Mitoses') shuffled_counts ``` | Shuffled Mitoses | Cancer | Not Cancer | | --- | --- | --- | | 1 | 199 | 364 | | 2 | 12 | 23 | | 3 | 12 | 21 | | 4 | 5 | 7 | | 5 | 2 | 4 | | 6 | 0 | 3 | | 7 | 3 | 6 | | 8 | 3 | 5 | | 10 | 3 | 11 | 這兩個類中的亂序數據的分布可以展示為條形圖，就像原始數據一樣。 ```py shuffled_dists = shuffled_counts.select(0).with_columns( 'Cancer', proportions(shuffled_counts.column(1)), 'Not Cancer', proportions(shuffled_counts.column(2)) ) shuffled_dists.barh(0) ```  這與原始條形圖看起來有點不同，為方便起見，再次展示如下。 ```py dists.barh(0) ```  ### 檢驗統計量：總變異距離 我們需要一個測試統計量來衡量藍色和金色分布之間的差異。 回想一下，總變異距離可以用來量化兩個類別分布的差異。 ```py def tvd(dist1, dist2): return 0.5*(np.sum(np.abs(dist1 - dist2))) ``` 在原始樣本中，兩個類別的有絲分裂的分布的 TVD 約為 0.4： ```py observed_tvd = tvd(dists.column(1), dists.column(2)) observed_tvd 0.41841946549059517 ``` 但是在亂序的樣本中，它比較小： ```py tvd(shuffled_dists.column(1), shuffled_dists.column(2)) 0.022173847487655045 ``` 隨機排列的有絲分裂評分和原始評分似乎表現不一樣。 但是如果我們再次運行，隨機打亂可能會有所不同。 讓我們重新打亂并重新計算總變異距離。 ```py shuffled_mitoses = mitoses.select('Mitoses').sample(with_replacement=False).column(0) shuffled = mitoses.select('Class').with_column('Shuffled Mitoses', shuffled_mitoses) shuffled_counts = shuffled.pivot('Class', 'Shuffled Mitoses') tvd(proportions(shuffled_counts.column(1)), proportions(shuffled_counts.column(2))) 0.039937426966715643 ``` 總變異距離仍然比我們從原始數據得到的 0.42 小很多。 為了看看它變化了多少，我們不得不重復多次隨機打亂過程，在它現在已經變得很熟悉了。 ### 原假設下 TVD 的經驗分布 如果原假設是真的，則有絲分裂評分的所有排列都是等可能的。 有很多可能的排列；讓我們做 5000 次，看看我們的檢驗統計量的變化。 代碼與上面的代碼完全一樣，只是現在我們將收集所有 5000 個距離并繪制經驗直方圖。 ```py repetitions = 5000 tvds = make_array() for i in np.arange(repetitions): shuffled_mitoses = mitoses.select('Mitoses').sample(with_replacement=False).column(0) shuffled = mitoses.select('Class').with_column('Shuffled Mitoses', shuffled_mitoses) shuffled_counts = shuffled.pivot('Class', 'Shuffled Mitoses') new_tvd = tvd(proportions(shuffled_counts.column(1)), proportions(shuffled_counts.column(2))) tvds = np.append(tvds, new_tvd) Table().with_column('TVD', tvds).hist(bins=20) plots.title('Empirical Distribution Under the Null') print('Observed TVD:', observed_tvd) Observed TVD: 0.418419465491 ```  觀察到的總變異距離 0.42 根本不接近于假設零假設為真所產生的分布。 數據支持備選假設：有絲分裂評分與類別有關。 ### 兩個類別分布的相等性的排列檢驗 我們上面所做的檢驗被稱為原假設的排列檢驗，即兩個樣本是從相同的底層分布中抽取的。 為了定義一個執行檢驗的函數，我們可以復制前一個單元格的代碼，并更改表和列的名稱。函數`permutation_test_tvd`接受數據表的名稱，包含類別變量的列標簽，它的分布要檢驗，包含二元類別變量的列標簽，以及要運行的隨機排列的數量。 在我們上面的例子中，我們沒有計算 P 值，因為觀測值遠離原假設下統計量的分布。但是，一般來說，我們應該計算 P 值，因為在其他例子中統計量可能不是那么極端。 P 值是“假設原假設為真，所得距離大于等于觀測距離”的幾率，因為備選假設比原假設預測了更大的距離。 ```py def permutation_test_tvd(table, variable, classes, repetitions): """Test whether a categorical variable is independent of classes: table: name of table containing the sample variable: label of column containing categorical variable whose distribution is of interest classes: label of column containing binary class data repetitions: number of random permutations""" # Find the tvd between the distributions of variable in the two classes counts = table.select(classes, variable).pivot(classes, variable) observed_tvd = tvd(proportions(counts.column(1)), proportions(counts.column(2))) # Assuming the null is true, randomly permute the variable and collect all the new tvd's tvds = make_array() for i in np.arange(repetitions): shuffled_var = table.select(variable).sample(with_replacement=False).column(0) shuffled = table.select(classes).with_column('Shuffled Variable', shuffled_var) shuffled_counts = shuffled.pivot(classes, 'Shuffled Variable') new_tvd =tvd(proportions(shuffled_counts.column(1)), proportions(shuffled_counts.column(2))) tvds = np.append(tvds, new_tvd) # Find the empirical P-value: emp_p = np.count_nonzero(tvds >= observed_tvd)/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('TVD', tvds).hist(bins=20) plots.title('Empirical Distribution Under the Null') print('Observed TVD:', observed_tvd) print('Empirical P-value:', emp_p) permutation_test_tvd(patients, 'Clump Thickness', 'Class', 5000) Observed TVD: 0.638310905047 Empirical P-value: 0.0 ```  同樣，觀測距離 0.64 離原假設預測的分布很遠。 經驗 P 值為 0，所以準確的 P 值將接近于零。 因此，如果類別和有絲分裂評分是不相關的，那么觀測的數據是極不可能的。 所以得出的結論是，有絲分裂評分與類別有關，不僅在樣本中，而且在總體中。 我們使用排列檢驗來幫助我們確定，類別屬性的分布是否與類別相關。 一般來說，排列檢驗可以這樣使用來確定，兩個類別分布是否從相同的基本分布隨機抽樣。 ## A/B 測試 我們使用隨機排列來查看，兩個樣本是否從相同的基本分類分布抽取。 如果樣本是數值的，則可以使用相同的方法；檢驗統計量的選擇通常比較簡單。 在我們使用`Deflategate`數據的例子中，我們使用了不同的方法來測試愛國者隊和小馬隊用球是否來自相同的基本分布。 在現代數據分析中，決定兩個數值樣本是否來自相同的基本分布稱為 A/B 測試。 名稱是指兩個樣本 A 和 B 的標簽。 ### 吸煙者和不吸煙者 我們對隨機抽樣的母親及其新生兒進行了許多不同的分析，但是我們還沒有查看母親是否吸煙的數據。 研究的目的之一，是看母親吸煙是否與出生體重有關。 ```py baby = Table.read_table('baby.csv') baby ``` | Birth Weight | Gestational Days | Maternal Age | Maternal Height | Maternal Pregnancy Weight | Maternal Smoker | | --- | --- | --- | --- | --- | --- | | 120 | 284 | 27 | 62 | 100 | False | | 113 | 282 | 33 | 64 | 135 | False | | 128 | 279 | 28 | 64 | 115 | True | | 108 | 282 | 23 | 67 | 125 | True | | 136 | 286 | 25 | 62 | 93 | False | | 138 | 244 | 33 | 62 | 178 | False | | 132 | 245 | 23 | 65 | 140 | False | | 120 | 289 | 25 | 62 | 125 | False | | 143 | 299 | 30 | 66 | 136 | True | | 140 | 351 | 27 | 68 | 120 | False | （省略了 1164 行） 我們首先選擇`Birth Weight`和`Maternal Smoker`。 樣本中有 715 名非吸煙者，459 名吸煙者。 ```py weight_smoke = baby.select('Birth Weight', 'Maternal Smoker') weight_smoke.group('Maternal Smoker') ``` | Maternal Smoker | count | | --- | --- | | False | 715 | | True | 459 | 下面的第一個直方圖顯示了樣本中非吸煙者的嬰兒出生體重的分布。 第二個顯示了吸煙者的嬰兒出生體重。 ```py nonsmokers = baby.where('Maternal Smoker', are.equal_to(False)) nonsmokers.hist('Birth Weight', bins=np.arange(40, 181, 5), unit='ounce') ```  ```py smokers = baby.where('Maternal Smoker', are.equal_to(True)) smokers.hist('Birth Weight', bins=np.arange(40, 181, 5), unit='ounce') ```  兩種分布都大致是鐘形，中心在 120 盎司附近。 當然，這些分布并不相同，這就產生了這樣一個問題，即差異是否僅僅反映了機會變異，還是反映了總體分布的差異。 這個問題可以通過假設檢驗來回答。 原假設：在總體中，不吸煙的母親的嬰兒出生體重的分布和吸煙的母親相同。 樣本中的差異是偶然的。 備選假設：兩種分布在總體中是不同的。 檢驗統計量：出生體重是一個定量變量，所以用均值的絕對差作為檢驗統計量是合理的。 檢驗統計量的觀測值約為 9.27 盎司。 ```py means_table = weight_smoke.group('Maternal Smoker', np.mean) means_table ``` | Maternal Smoker | Birth Weight mean | | --- | --- | | False | 123.085 | | True | 113.819 | ```py nonsmokers_mean = means_table.column(1).item(0) smokers_mean = means_table.column(1).item(1) nonsmokers_mean - smokers_mean 9.266142572024918 ``` ## 排列檢驗 為了看看原假設下是否有可能出現這種差異，我們將使用排列檢驗，就像我們在前一節中所做的那樣。 我們必須為檢驗統計量改變代碼。 為此，我們將像上面那樣計算平均值的差，然后取絕對值。 請記住，在原假設下，出生體重的所有排列與`Maternal Smoker `列等可能出現。 所以，就像以前一樣，每次重復都是打亂正在比較的變量。 ```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(bins=20) plots.title('Empirical Distribution Under the Null') print('Observed statistic:', obs_stat) print('Empirical P-value:', emp_p) permutation_test_means(baby, 'Birth Weight', 'Maternal Smoker', 5000) Observed statistic: 9.266142572024918 Empirical P-value: 0.0 ```  原始樣本中的觀測差異約為 9.27 盎司，與此分布不一致：經驗 P 值為 0，這意味著確切的 P 值確實非常小。 因此，測試的結論是，在總體中，不吸煙者和吸煙者的嬰兒出生體重的分布是不同的。 ## 差值的自舉置信區間 我們的 A/B 測試得出結論，這兩個分布是不同的，但有點不盡人意。他們有多么不同？哪一個均值更大？這些自然是測試無法回答的問題。 回想一下，我們之前已經討論過這個問題了：不僅僅是問“兩個分布是否不同”的是與否的問題，我們可以通過不作任何假設，并簡單地估計均值之間的差異，來學到更多。 觀測差異（不吸煙者減去吸煙者）約為 9.27 盎司；這個正面跡象表明，不吸煙的母親通常有更大的嬰兒。但由于隨機性，樣本可能會有所不同。為了了解有多么不同，我們必須生成更多的樣本；為了生成更多的樣本，我們將使用`bootstrap`，就像我們以前做過的那樣。自舉過程不會假設這兩個分布是否相同。它只是復制原始隨機樣本并計算統計量的新值。 函數`bootstrap_ci_means`返回總體中兩組均值之間差異的自舉置信區間。在我們的例子中，置信區間將估計總體中吸煙和不吸煙的母親的嬰兒的平均出生體重之間的差異。 + 表名稱，它包含原始樣本中的數據 + 列標簽，它包含數值變量 + 列標簽，它包含兩個樣本的名稱 + 自舉的重復次數 該函數使用自舉百分比方法，返回兩個均值之間的差異的約 95% 置信區間。 ```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(baby, 'Birth Weight', 'Maternal Smoker', 5000) Observed difference between means: 9.266142572024918 Approximate 95% CI for the difference between means: 7.23940878698 to 11.3907887554 ```  不吸煙的母親的嬰兒比吸煙的母親的嬰兒平均重 7.2 盎司到 11.4 盎司。 這比“兩個分布不同”更有用。 由于置信區間不包含 0，它也告訴我們這兩個分布是不同的。 所以置信區間估計了我們的均值之間的差異，也讓我們決定兩個基本分布是否相同。 不吸煙的母親比吸煙的母親平均年齡稍大。 ```py bootstrap_ci_means(baby, 'Maternal Age', 'Maternal Smoker', 5000) Observed difference between means: 0.8076725017901509 Approximate 95% CI for the difference between means: 0.154278698588 to 1.4701157656 ```  但毫不奇怪，證據并沒有指出，他們的平均身高與不吸煙的母親不同。 零在均值之間差異的置信區間中。 ```py bootstrap_ci_means(baby, 'Maternal Height', 'Maternal Smoker', 5000) Observed difference between means: 0.09058914941267915 Approximate 95% CI for the difference between means: -0.390841928035 to 0.204388297872 ```  總之： 如果你想知道兩個基本分布是否相同，則可以使用帶有適當檢驗統計量的排列檢驗。 當分布是類別時，我們使用總變異距離，而分布是數值時，我們使用均值之間的絕對差。 為了比較兩個數值分布，將假設檢驗替換為估計，通常更富有信息。 只需估計一個差異，比如兩組均值之間的差異。 這可以通過構建自舉置信區間來完成。 如果零不在這個區間內，你可以得出這樣的結論：這兩個分布是不同的，你也可以估計均值有多么不同。 ## 因果 我們用于比較兩個樣本的方法在隨機對照實驗的分析中具有強大的用途。由于在這些實驗中，實驗組和對照組被隨機分配，因此如果實驗完全沒有效果，結果中的任何差異，可以與僅僅由于分配中的隨機性而發生的情況進行比較。如果觀察到的差異比我們預測的，純粹由于偶然的差異更為顯著，我們就會有因果關系的證據。由于個體無偏分配到實驗組和對照組，兩組結果中的差異可歸因于實驗。 隨機對照實驗分析的關鍵是，了解偶然因素如何出現。這有助于我們設定明確的原假設和備選假設。一旦完成，我們可以簡單地使用前一節的方法來完成分析。 讓我們看看如何在一個例子中實現。 ### 治療慢性背痛：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 ```  基本分布的均值之間的差異的約 95% 置信區間，范圍是約 -80% 到 -20%。換句話說，實驗組好轉了 20% 到 80% 左右。 注意這個變化很大的估計。那是因為每個組的樣本量只有 15 個左右。雖然這些作用于這些數值而沒有進一步的假設，但結果并不十分精確。 ### 元分析 雖然 RCT 確實真名了肉毒桿菌毒素 A 實驗幫助了患者，但對 31 名患者進行的研究不足以確定治療的有效性。這不僅僅是因為樣本量小。我們在這一部分的結果對于研究中的 31 位患者是有效的，但我們對所有可能患者的總體真正感興趣。如果 31 名患者是來自較大總體的隨機樣本，那么我們的置信區間對該總體是有效的。但他們不是隨機樣本。 2011 年，一組研究人員對實驗的研究進行了元分析。也就是說，他們確定了所有被痛治療的可用研究，并總結了整理后的結果。 有幾項研究，但沒有多少可以納入科學合理的方式：“由于非隨機性，不完整或未發表的數據，我們排除了 19 項研究的證據。只剩下三個隨機對照試驗，其中之一是我們在本節研究的。元分析給予它所有研究的最高評價（LBP 代表背痛）：“我們確定了三項研究，它們調查了 BoNT 治療 LBP 的優點，但只有一項的偏差風險低，并且使用非特異性 LBP（N = 31）來評價患者”。 元分析得出的結論是：“有一些低質量的證據表明，BoNT 注射劑能改善疼痛，功能，或者兩者都比注射生理鹽水更好，而且質量很低的證據表明，它比針灸或類固醇注射更好。進一步的研究很可能會對效果評估和我們的信心產生重要影響，未來的試驗應該對患者總體，實驗方案和比較組進行標準化，爭取更多的參與者，并包括長期結果，成本效益分析和臨床相關性的發現”。 為了確定醫療有好處，需要很多精心的工作。了解如何分析隨機對照試驗是這項工作的重要組成部分。現在你們知道了如何實現，你們有條件幫助醫療和其他行業建立因果關系。