獨家 | 如何比較兩個(gè)或多個(gè)分布形態(tài)（附鏈接）

翻譯：陳超

校對：趙茹萱

圖片來(lái)自作者

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

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

數據快照，圖片來(lái)自作者

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

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

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

處理組合對照組的收入分布，圖片來(lái)自作者

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

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

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

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

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

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

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

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

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

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");

處理組和對照組的累積分布圖，圖片來(lái)自作者

我們應該如何解釋這幅圖？

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

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

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))

現在，我們可以將兩個(gè)分位數分布相互對照，加上45度線(xià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 圖, 圖片來(lái)自作者
Q-Q圖提供了與累積分布圖非常相似的見(jiàn)解:處理組的收入有相同的中位數(在中心交叉的線(xiàn))，但更寬的尾部(點(diǎn)在左邊的線(xiàn)以下，右邊的線(xiàn)以上)。
兩組——檢驗
到目前為止，我們已經(jīng)看到了可視化分布之間差異的不同方法?？梢暬闹饕獌?yōu)點(diǎn)是直觀(guān):我們可以通過(guò)肉眼觀(guān)察差異并直觀(guān)地評估它們。
然而，我們可能想要更嚴格地評估分布之間的差異的統計意義，即回答這個(gè)問(wèn)題“觀(guān)察到的差異是系統的還是由于采樣噪聲?”
我們現在將分析不同的測試來(lái)辨別兩個(gè)分布。

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

T檢驗統計，圖片來(lái)自作者
式中為樣本均值，s為樣本標準差。在較溫和的條件下，檢驗統計量是漸近分布的Student t分布。
我們使用scipy中的ttest_ind函數來(lái)執行t檢驗。該函數返回測試統計數據和隱含的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'])

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

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

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');

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

# 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

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

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()

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])

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");

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

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

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

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");

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

參考文獻

[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.