翻譯：陳超

校對：趙茹萱

圖片來自作者

比較同一變量在不同組別之間的經(jīng)驗分布是數(shù)據(jù)科學(xué)當(dāng)中的常見問題，尤其在因果推斷中，我們經(jīng)常在需要評估隨機化質(zhì)量時遇到上述問題。
我們想評估某一政策的效果（或者用戶體驗功能，廣告宣傳，****物，……），因果推斷當(dāng)中的金標(biāo)準(zhǔn)就是隨機對照試驗，也叫作A/B測試。在實際情況下，我們會選擇一個樣本進(jìn)行研究，隨機分為對照組和實驗組，并且比較兩組之間結(jié)果差異。隨機化能夠確保兩組間唯一的差異是是否接受治療，平均而言，以便于我們可以將結(jié)果差異歸因于治療效應(yīng)。
問題是，盡管進(jìn)行了隨機化，兩組也不會完全相同。有時，他們甚至不是“相似的”。例如，我們可能會在一組中有更多男性或年齡更大的人，等等（我們通常把這些叫做特質(zhì)協(xié)變量或控制變量）。這種情況發(fā)生時，我們再也無法確定結(jié)果的差異僅僅是由治療的效果導(dǎo)致，也不能將其完全歸因于不平衡的協(xié)變量。因此，隨機化之后非常重要的一步就是檢查是否所有觀測變量都是組間平衡的，是否不存在系統(tǒng)性差異。另外一個選擇是分層抽樣，額可以事先確保特定協(xié)變量是平衡的。
在本文中，我們將通過不同方式比較兩組（或多組）分布并評估他們之間差異的量級和顯著性水平。可視化和統(tǒng)計角度這兩種方法通常是嚴(yán)謹(jǐn)性和直覺的權(quán)衡：從圖上，我們可以迅速評估和探究差異，但是很難區(qū)分這些差異是否是系統(tǒng)性的還是僅僅由于噪聲導(dǎo)致。
例子
假設(shè)我們需要將一組人隨機分到處理組和對照組。我們需要讓兩組盡可能地相似，以便于將組間差異歸因于治療效應(yīng)。我們也需要將處理組分成幾個亞組來測試不同治療的影響（例如，同一種****物的細(xì)微變化）。
對于這個例子來說，我們已經(jīng)模擬了1000個被試數(shù)據(jù)集，我從src.dgp導(dǎo)入了數(shù)據(jù)生成過程dgp_rnd_assignment()，并從src.utils導(dǎo)入了一些繪圖函數(shù)和庫，從而觀測到一系列特征。

from src.utils import *from src.dgp import dgp_rnd_assignment
df = dgp_rnd_assignment().generate_data()df.head()

數(shù)據(jù)快照，圖片來自作者

我們有1000個被試的信息，從中可以觀測到性別、年齡和周收入。每個被試被分配到處理組或?qū)φ战M，被分到處理組的被試又被分到四種不同的治療亞組當(dāng)中去。
兩組-圖
讓我們從最簡單的情況開始：比較處理組和對照組的收入分布。首先用可視化方法來進(jìn)行探究，然后再使用統(tǒng)計方法?？梢暬椒ǖ膬?yōu)勢在于直觀，而統(tǒng)計方法方法的優(yōu)勢則在于嚴(yán)謹(jǐn)。
對大多數(shù)可視化來說，我會使用python當(dāng)中的searborn庫。

箱線圖
第一種可視化方法是箱線圖。箱線圖是統(tǒng)計概要和數(shù)據(jù)可視化之間的很好的兌易。箱體的中心表征中位數(shù)，上下邊界則表征第1和第3百分位數(shù)。須體延長到超過箱體四分位數(shù)（Q3-Q1）1.5倍的第一個數(shù)據(jù)點。落在須體之外的點則分別繪制，且通常被視作異常值。
因此，箱線圖提供了統(tǒng)計概要（箱體和須體）以及直觀的數(shù)據(jù)可視化（異常值）。

sns.boxplot(data=df, x='Group', y='Income');plt.title("Boxplot");

處理組合對照組的收入分布，圖片來自作者

看起來處理組的收入分布更加分散：橘色箱體更大，須體覆蓋范圍更廣。然而，箱線圖的問題在于它隱藏了數(shù)據(jù)的形態(tài)，僅僅告訴我們統(tǒng)計概要而未向我們展示真實的數(shù)據(jù)分布情況。
直方圖
直方圖是展示分布最直觀的方式，它將數(shù)據(jù)分成同等寬度的組，將每組觀測值數(shù)量畫出來。

sns.histplot(data=df, x='Income', hue='Group', bins=50);plt.title("Histogram");

處理組和對照組的收入分布情況，圖片來自作者

該圖也存在很多問題：
因為兩組觀測值數(shù)量不同，兩個直方圖不具備可比性。
分組的數(shù)量是武斷的。
我們可以通過stat選項來解決第一個方法，繪制density而非計數(shù)，將common_norm選項設(shè)置為False來分別對每個直方圖進(jìn)行歸一化。

sns.histplot(data=df, x='Income', hue='Group', bins=50, stat='density', common_norm=False);plt.title("Density Histogram");

處理組和對照組的分布，圖片來自作者
現(xiàn)在兩組直方圖就可比較了！
然而，一個重要的問題仍然存在：分組的大小是武斷的。在極端情況下，如果我們把更少的數(shù)據(jù)捆綁在一起，最后會得到每組至多一條觀測數(shù)據(jù)，如果我們把更多的數(shù)據(jù)捆綁在一起，我們最終可能會得到一個組。在兩種情況下，如果我們夸大，圖就會損失信息量。這就是經(jīng)典的偏差-變異兌易。
核密度圖
一種可能的解決方法是使用核密度函數(shù)，使用核密度估計(KDE)用連續(xù)函數(shù)近似直方圖。

sns.kdeplot(x='Income', data=df, hue='Group', common_norm=False);plt.title("Kernel Density Function");

處理組和對照組收入分布，圖片來自作者
從圖上可以看出，似乎處理組的收入的估計核密度有“更胖的尾巴”（更高的方差），但組間均值更為相似。
核密度估計的問題自安于它是一個黑箱，可能會掩蓋數(shù)據(jù)的相關(guān)特征。
累積分布圖
一種更為透明的表征兩個分布的方法是累積分布函數(shù)。在x軸的每個點（收入）我們繪制出數(shù)值相等或更低的數(shù)據(jù)點的百分比。累積分布函數(shù)的優(yōu)勢在于：

• 我們不需要做出任何武斷的決策（例如，分組數(shù)量）
• 我們不需要做任何近似（例如：KDE），但是我們可以表征所有的數(shù)據(jù)點

sns.histplot(x='Income', data=df, hue='Group', bins=len(df), stat="density",element="step", fill=False, cumulative=True, common_norm=False);plt.title("Cumulative distribution function");

處理組和對照組的累積分布圖，圖片來自作者

我們應(yīng)該如何解釋這幅圖？

• 兩條線在0.5（y軸）附近交叉，意味著他們的中位數(shù)相似
• 在左側(cè)橘色線在藍(lán)色線上，而右側(cè)則相反，意味著處理組分布的尾部更胖（極端值更多）

Q-Q圖
一個相關(guān)的方法是Q-Q圖，其中Q代表分位數(shù)。Q-Q圖將兩個分布的分位數(shù)相互繪制出來。如果分布相同，就會得到45度的直線。
Python中沒有本地的Q-Q圖函數(shù)，雖然statmodels包提供了一個qqplot函數(shù)，但它相當(dāng)麻煩。因此，我們需要手動完成。
首先，我們需要使用percentile函數(shù)計算兩組的四分位數(shù)。

income = df['Income'].valuesincome_t = df.loc[df.Group=='treatment', 'Income'].valuesincome_c = df.loc[df.Group=='control', 'Income'].valuesdf_pct = pd.DataFrame()df_pct['q_treatment'] = np.percentile(income_t, range(100))df_pct['q_control'] = np.percentile(income_c, range(100))

現(xiàn)在，我們可以將兩個分位數(shù)分布相互對照，加上45度線，表示基準(zhǔn)的完美擬合。

plt.figure(figsize=(8, 8))plt.scatter(x='q_control', y='q_treatment', data=df_pct, label='Actual fit');sns.lineplot(x='q_control', y='q_control', data=df_pct, color='r', label='Line of perfect fit');plt.xlabel('Quantile of income, control group')plt.ylabel('Quantile of income, treatment group')plt.legend()plt.title("QQ plot");

Q-Q 圖, 圖片來自作者
Q-Q圖提供了與累積分布圖非常相似的見解:處理組的收入有相同的中位數(shù)(在中心交叉的線)，但更寬的尾部(點在左邊的線以下，右邊的線以上)。
兩組——檢驗
到目前為止，我們已經(jīng)看到了可視化分布之間差異的不同方法。可視化的主要優(yōu)點是直觀:我們可以通過肉眼觀察差異并直觀地評估它們。
然而，我們可能想要更嚴(yán)格地評估分布之間的差異的統(tǒng)計意義，即回答這個問題“觀察到的差異是系統(tǒng)的還是由于采樣噪聲?”
我們現(xiàn)在將分析不同的測試來辨別兩個分布。

T檢驗
第一個也是最常見的檢驗是學(xué)生t檢驗。t檢驗通常用于比較平均值。在這種情況下，我們希望測試兩組的收入分配均值是否相同。兩均值比較檢驗的檢驗統(tǒng)計量為：

T檢驗統(tǒng)計，圖片來自作者
式中為樣本均值，s為樣本標(biāo)準(zhǔn)差。在較溫和的條件下，檢驗統(tǒng)計量是漸近分布的Student t分布。
我們使用scipy中的ttest_ind函數(shù)來執(zhí)行t檢驗。該函數(shù)返回測試統(tǒng)計數(shù)據(jù)和隱含的p值。

from scipy.stats import ttest_indstat, p_value = ttest_ind(income_c, income_t)print(f"t-test: statistic={stat:.4f}, p-value={p_value:.4f}")t-test: statistic=-1.5549, p-value=0.1203

from causalml.match import create_table_onedf['treatment'] = df['Group']=='treatment'create_table_one(df, 'treatment', ['Gender', 'Age', 'Income'])

在前兩列中，我們可以看到處理組和對照組不同變量的平均值，括號中是標(biāo)準(zhǔn)誤差。在最后一列，SMD的值表明所有變量的標(biāo)準(zhǔn)化差異大于0.1，表明兩組可能是不同的。

Mann–Whitney U 檢驗
另一種可選的檢驗是Mann–Whitney U 檢驗。零假設(shè)是兩組有相同的粉不，而備擇假設(shè)是一組的值比另一組更大（或更?。?。
不同于我們之前看過的檢驗，Mann–Whitney U 檢驗不關(guān)注異常值，而把注意力放在分布的中心上。
檢驗流程如下。
1.將所有數(shù)據(jù)點合并排序（升序或降序）2.計算U? = R? ? n?(n? + 1)/2, R?是第一組的秩和，n?是第一組數(shù)據(jù)的數(shù)量。3.用相似的方法計算第二組的U?4.統(tǒng)計檢驗量是stat = min(U?, U?)
在兩個分布之間沒有系統(tǒng)秩差(即中位數(shù)相同)的零假設(shè)下，檢驗統(tǒng)計量在均值和方差已知的情況下，是漸近正態(tài)分布的。
計算R和U的直觀方法是:如果第一個樣品的值都大于第二個樣品的值，那么R?= n?(n?+ 1)/2，因此，U?將為零(可得到的最小值)。否則，如果兩個樣本相似，U?和U?就會非常接近n?n?/ 2(可得到的最大值)。
我們使用來自scipy的mannwhitneyu函數(shù)執(zhí)行測試。

from scipy.stats import mannwhitneyustat, p_value = mannwhitneyu(income_t, income_c)print(f" Mann–Whitney U Test: statistic={stat:.4f}, p-value={p_value:.4f}")Mann–Whitney U Test: statistic=106371.5000, p-value=0.6012

sample_stat = np.mean(income_t) - np.mean(income_c)stats = np.zeros(1000)for k in range(1000):    labels = np.random.permutation((df['Group'] == 'treatment').values)    stats[k] = np.mean(income[labels]) - np.mean(income[labels==False])p_value = np.mean(stats > sample_stat)print(f"Permutation test: p-value={p_value:.4f}")Permutation test: p-value=0.0530

plt.hist(stats, label='Permutation Statistics', bins=30);plt.axvline(x=sample_stat, c='r', ls='--', label='Sample Statistic');plt.legend();plt.xlabel('Income difference between treatment and control group')plt.title('Permutation Test');

置換間的均值差異分布，圖片來自作者
正如我們所看到的，相對于排列后的樣本值，樣本統(tǒng)計值是相當(dāng)極端的，但也合理。
卡方檢驗
卡方檢驗是一個效力很強的檢驗，常用于檢驗頻率差異。
卡方檢驗最不為人知的應(yīng)用之一是檢驗兩個分布之間的相似性。把兩組觀測值分組。如果這兩個分布是相同的，我們將期望在每個組中有相同的觀測頻率。重要的是，我們需要每個組內(nèi)有足夠多的觀測值，以保證測試的有效性。
我生成對應(yīng)于對照組收入分布十分位數(shù)的組，然后計算處理組中每個組別的預(yù)期觀察值頻數(shù)，來確定兩種分布是否相同。

# Init dataframedf_bins = pd.DataFrame()# Generate bins from control group_, bins = pd.qcut(income_c, q=10, retbins=True)df_bins['bin'] = pd.cut(income_c, bins=bins).value_counts().index# Apply bins to both groupsdf_bins['income_c_observed'] = pd.cut(income_c, bins=bins).value_counts().valuesdf_bins['income_t_observed'] = pd.cut(income_t, bins=bins).value_counts().values# Compute expected frequency in the treatment groupdf_bins['income_t_expected'] = df_bins['income_c_observed'] / np.sum(df_bins['income_c_observed']) * np.sum(df_bins['income_t_observed'])df_bins

from scipy.stats import chisquarestat, p_value = chisquare(df_bins['income_t_observed'], df_bins['income_t_expected'])print(f"Chi-squared Test: statistic={stat:.4f}, p-value={p_value:.4f}")Chi-squared Test: statistic=32.1432, p-value=0.0002

原因在于兩個分布有一個相似的中心，但尾部不同。而卡方檢驗檢驗的是整個分布的相似性，而不是像之前檢驗?zāi)菢又辉谥行摹?/span>
這個結(jié)果告訴我們:在從p值得出盲目結(jié)論之前，了解您實際測試的是什么是非常重要的！

Kolmogorov-Smirnov檢驗
Kolmogorov-Smirnov檢驗可能是比較分布最流行的非參數(shù)檢驗。Kolmogorov-Smirnov檢驗的思想是比較兩組的累積分布。特別是，Kolmogorov-Smirnov檢驗統(tǒng)計量是兩個累積分布之間的最大絕對差值。
Kolmogorov-Smirnov檢驗統(tǒng)計量,圖片來自作者
其中F?和F?為兩個累積分布函數(shù)，x為基礎(chǔ)變量的值。Kolmogorov- smirnov檢驗統(tǒng)計量的漸近分布是Kolmogorov分布。
為了更好地理解檢驗，讓我們畫出累積分布函數(shù)和檢驗統(tǒng)計量。首先，我們計算累積分布函數(shù)。

• 
  

• 
  

• 
  

• 
  

df_ks = pd.DataFrame()df_ks['Income'] = np.sort(df['Income'].unique())df_ks['F_control'] = df_ks['Income'].apply(lambda x: np.mean(income_c<=x))df_ks['F_treatment'] = df_ks['Income'].apply(lambda x: np.mean(income_t<=x))df_ks.head()

累計分布數(shù)據(jù)集快照，圖片來自作者
我們現(xiàn)在需要找到累積分布函之間的絕對距離最大的點。

• 
  

• 
  

k = np.argmax( np.abs(df_ks['F_control'] - df_ks['F_treatment']))ks_stat = np.abs(df_ks['F_treatment'][k] - df_ks['F_control'][k])

我們可以通過繪制兩個累積分布函數(shù)和測試統(tǒng)計量的值來可視化測試統(tǒng)計量的值。

• 
  

• 
  

• 
  

• 
  

• 
  

y = (df_ks['F_treatment'][k] + df_ks['F_control'][k])/2plt.plot('Income', 'F_control', data=df_ks, label='Control')plt.plot('Income', 'F_treatment', data=df_ks, label='Treatment')plt.errorbar(x=df_ks['Income'][k], y=y, yerr=ks_stat/2, color='k',capsize=5, mew=3, label=f"Test statistic: {ks_stat:.4f}")plt.legend(loc='center right');plt.title("Kolmogorov-Smirnov Test");

Kolmogorov-Smirnov檢驗統(tǒng)計量,圖片來自作者
從圖中，我們可以看到檢驗統(tǒng)計量的值對應(yīng)于收入~650處兩個累積分布之間的距離。就收入價值而言，兩組之間的不平衡是最大的。
現(xiàn)在我們可以使用來自scipy的kstest函數(shù)執(zhí)行實際的測試。

• 
  

• 
  

• 
  

from scipy.stats import ksteststat, p_value = kstest(income_t, income_c)print(f" Kolmogorov-Smirnov Test: statistic={stat:.4f}, p-value={p_value:.4f}")Kolmogorov-Smirnov Test: statistic=0.0974, p-value=0.0355

p值低于5%:我們以95%的置信度拒絕兩個分布相同的零假設(shè)。
注1:KS檢驗過于保守，很少拒絕零假設(shè)。Lilliefors檢驗使用測試統(tǒng)計量的不同分布(Lilliefors分布)校正了這一偏差。
注2:KS測試使用的信息很少，因為它只比較在一點上的兩個累積分布:最大距離的一個。Anderson-Darling檢驗和Cramér-von Mises檢驗通過積分來比較整個域上的兩個分布(兩者之間的差異在于平方距離的加權(quán))。
多組-圖
到目前為止，我們只考慮了兩組的情況:處理組和對照組。但如果我們有多個組呢?我們看到的一些方法可以很好地擴展，而另一些則不行。
作為一個可行的例子，我們現(xiàn)在要檢查不同處理組的收入分布是否相同。
箱線圖
當(dāng)我們有許多個位數(shù)的組時，箱線圖可以很好地縮放，因為我們可以把不同的盒子并排放在一起。

• 
  

• 
  

sns.boxplot(x='Arm', y='Income', data=df.sort_values('Arm'));plt.title("Boxplot, multiple groups");

不同處理亞組收入分布，圖片來自作者
從圖中可以看出，不同組別的收入分布是不同的，編號越高的組別平均收入越高。
小提琴圖
結(jié)合了匯總統(tǒng)計和核密度估計的箱線圖的一個很好的擴展是小提琴圖。小提琴圖顯示了沿y軸的獨立密度，所以他們不會重疊。默認(rèn)情況下，它還在內(nèi)部添加一個微型箱線圖。

• 
  

• 
  

sns.violinplot(x='Arm', y='Income', data=df.sort_values('Arm'));plt.title("Violin Plot, multiple groups");

不同處理亞組收入分布，圖片來自作者
相較于箱線圖，小提琴圖能表明不同處理組的收入不同。
脊線圖
最后，脊線圖沿x軸繪制多個核密度分布，比小提琴圖更直觀，但部分重疊。不幸的是，在matplotlib和seaborn中都沒有默認(rèn)的脊線圖。我們需要從joypy導(dǎo)入它。

• 
  

• 
  

from joypy import joyplotjoyplot(df, by='Arm', column='Income', colormap=sns.color_palette("crest", as_cmap=True));plt.xlabel('Income');plt.title("Ridgeline Plot, multiple groups");

不同處理亞組收入分布，圖片來自作者
同樣，脊線圖表明，編號越高的治療組收入越高。從這個圖中，也更容易理解分布的不同形狀。
多組-檢驗
最后，讓我們考慮假設(shè)檢驗來比較多個組。為了簡單起見，我們將集中討論最常用的一個:F檢驗。
F-檢驗
對于多個組，最常用的測試是f測試。f檢驗比較一個變量在不同組之間的方差。這種分析也被稱為方差分析，或ANOVA。
在實際應(yīng)用中，F(xiàn)檢驗統(tǒng)計量由
F檢驗統(tǒng)計量，圖片來自作者
其中G為組數(shù)，N為觀察次數(shù)，為總體均值，g為g組內(nèi)均值。在組獨立性零假設(shè)下，f統(tǒng)計量為F分布。

• 
  

• 
  

• 
  

• 
  

from scipy.stats import f_onewayincome_groups = [df.loc[df['Arm']==arm, 'Income'].values for arm in df['Arm'].dropna().unique()]stat, p_value = f_oneway(*income_groups)print(f"F Test: statistic={stat:.4f}, p-value={p_value:.4f}")F Test: statistic=9.0911, p-value=0.0000

檢驗p值基本為零，這意味著強烈拒絕零假設(shè)，即各治療組之間的收入分配沒有差異。
結(jié)論
在這篇文章中，我們已經(jīng)看到了大量不同的方法來比較兩個或多個分布，無論是視覺上的還是統(tǒng)計上的。這是許多應(yīng)用的主要關(guān)注點，在因果推斷中尤其如此，我們使用隨機化方法使處理組和對照組盡可能具有可比性。
我們還看到不同的方法可能更適合不同的情況。可視化方法有助于建立直覺，但統(tǒng)計方法對于決策至關(guān)重要，因為我們需要能夠評估差異的量級和統(tǒng)計意義。

參考文獻(xiàn)

[1] Student, The Probable Error of a Mean (1908), Biometrika.

[2] F. Wilcoxon, Individual Comparisons by Ranking Methods (1945), Biometrics Bulletin.

[3] B. L. Welch, The generalization of “Student’s” problem when several different population variances are involved (1947), Biometrika.

[4] H. B. Mann, D. R. Whitney, On a Test of Whether one of Two Random Variables is Stochastically Larger than the Other (1947), The Annals of Mathematical Statistics.

[5] E. Brunner, U. Munzen, The Nonparametric Behrens-Fisher Problem: Asymptotic Theory and a Small-Sample Approximation (2000), Biometrical Journal.

[6] A. N. Kolmogorov, Sulla determinazione empirica di una legge di distribuzione (1933), Giorn. Ist. Ital. Attuar..

[7] H. Cramér, On the composition of elementary errors (1928), Scandinavian Actuarial Journal.

[8] R. von Mises, Wahrscheinlichkeit statistik und wahrheit (1936), Bulletin of the American Mathematical Society.

[9] T. W. Anderson, D. A. Darling, Asymptotic Theory of Certain “Goodness of Fit” Criteria Based on Stochastic Processes (1953), The Annals of Mathematical Statistics.

*博客內(nèi)容為網(wǎng)友個人發(fā)布，僅代表博主個人觀點，如有侵權(quán)請聯(lián)系工作人員刪除。