箱线图
第一种视觉方法是箱线图。箱线图是汇总统计和数据可视化之间的良好折衷。框的中心代表中位数,而边框分别代表第1(Q1)和第3四分位数(Q3)。扩展线延伸到框外超过四分位距 (Q3 - Q1) 1.5 倍的第一个数据点。落在扩展线之外的点是单独绘制的,通常会被认为是异常值。
因此,箱线图提供了汇总统计数据(方框和扩展线)和直接数据可视化(异常值)。
sns.boxplot(data=df, x='Group', y='Income');
plt.title("Boxplot");
实验组和对照组的收入分配情况
实验组的收入分配更加分散:橙色盒子更大,它的扩展线覆盖范围更广。但是箱线图的问题是它隐藏了数据的形状,它告诉我们一些汇总的统计数据,但没有显示实际的数据分布。
直方图
绘制分布图最直观的方法是直方图。直方图将数据分组到同等宽的容器(bin)中,并绘制出每个容器中的观察数据的数量。
sns.histplot(data=df, x='Income', hue='Group', bins=50);
plt.title("Histogram");
直方图也存在一些问题
-
由于两组的观察次数不同,因此两个直方图不具有可比性
-
bin的数量是任意的
我们可以使用 stat 选项来绘制密度而不是计数来解决第一个问题,并将 common_norm 设置为 False 分别对每个直方图进行归一化。
sns.histplot(data=df, x='Income', hue='Group', bins=50, stat='density', common_norm=False);
plt.title("Density Histogram");
这样这两个直方图就具有可比性!
但是一个重要的问题仍然存在:bin的大小是任意的。在极端特殊的情况下,如果我们将数据更少分组,最终会得到最多只有一个观察值的 bin,如果我们将数据分组更多,我们最终会只得到一个 bin。在这两种情况下,我们都无法判断。这是一个经典的偏差-方差权衡的问题。
核密度
一种可能的解决方案是使用核密度函数,该函数尝试使用核密度估计 (KDE) 用连续函数逼近直方图。
sns.kdeplot(x='Income', data=df, hue='Group', common_norm=False);
plt.title("Kernel Density Function");
从图中可以看到,收入核密度似乎在实验组中具有更高的方差,但是各组的平均值却是相似的。核密度估计的问题在于它有点像一个黑匣子,可能会掩盖数据的相关特征。
累积分布
两种分布更透明的表示是它们的累积分布函数(Cumulative Distribution Function)。在 x 轴(收入)的每个点,我们绘制具有相等或更低值的数据点的百分比。累积分布函数的主要优点是
-
不需要做出任何的选择(例如bin的数量)
-
不需要执行任何近似(例如使用 KDE),图形代表所有数据点
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");
我们应该如何解释这张图?
-
由于这两条线在0.5 (y轴)处或多或少交叉,这意味着它们的中值是相似的
-
因为橙色线在左边的蓝线之上,在右边的蓝线之下,这意味着实验组的分布是fatter tails(肥尾)
QQ图
一种相关的方法是 QQ 图,其中 q 代表分位数。QQ 图绘制了两个分布的分位数。如果分布相同应该得到一条 45 度线。
Python 中没有原生的 QQ 图功能,而 statsmodels 包提供了 qqplot 功能,但相当麻烦。因此,我们将手动完成。
首先,我们需要使用 percentile 函数计算两组的四分位数。
income = df['Income'].values
income_t = df.loc[df.Group=='treatment', 'Income'].values
income_c = df.loc[df.Group=='control', 'Income'].values
df_pct = pd.DataFrame()
df_pct['q_treatment'] = np.percentile(income_t, range(100))
df_pct['q_control'] = np.percentile(income_c, range(100))
现在我们可以绘制两个分位数分布,加上 45 度线,代表基准完美拟合。
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");
QQ 图在累积分布图方面提供了非常相似的信息:实验组的收入具有相同的中位数(线在中心交叉)但尾部更宽(点在左边的线以下,右边的线以上)。
T检验
第一个也是最常见的是学生 t 检验。T 检验通常用于比较均值。我们要检验两组的收入分配均值是否相同。两均值比较检验的检验统计量由下式给出:
其中 x̅ 是样本均值,s 是样本标准差。在较温和的条件下,检验统计量作为学生 t 分布渐近分布
我们使用 scipy 中的 ttest_ind 函数来执行 t 检验。该函数返回检验统计量和隐含的 p 值。
from scipy.stats import ttest_ind
stat, 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
检验的p值为0.12,因此我们不拒绝实验组和对照组平均值无差异的零假设。
标准化平均差 (SMD)
一般来说,当我们进行随机对照试验或 A/B 测试时,最好对实验组和对照组中所有变量的均值差异进行检验。
然而,由于 t 检验统计量的分母取决于样本量,因此 t 检验的 p 值难以跨研究进行比较。所以我们可能在一个差异非常小但样本量很大的实验中获得显着的结果,而在差异很大但样本量小的实验中我们可能会获得不显着的结果。
解决这个问题的一种解决方案是标准化平均差 (SMD)。顾名思义,这不是一个适当的统计量,而只是一个标准化的差异,可以计算为:
通常,低于0.1的值被认为是一个“小”的差异。
最将实验组和对照组的所有变量的平均值以及两者之间的距离度量(t 检验或 SMD)收集到一个称为平衡表的表中。可以使用causalml库中的create_table_one函数来生成它。正如该函数的名称所显示的那样,在执行A/B测试时,平衡表应该是你希望看到的的第一个表。
from causalml.match import create_table_one
df['treatment'] = df['Group']=='treatment'
create_table_one(df, 'treatment', ['Gender', 'Age', 'Income'])
平衡表
在前两列中,我们可以看到实验组和对照组的不同变量的平均值,括号中是标准误差。在最后一列中,SMD 的值表示所有变量的标准化差异均大于 0.1,这表明两组可能不同。
Mann–Whitney U检验
另一种检验是 Mann-Whitney U 检验,它比较两个分布的中位数。该检验的原假设是两组具有相同的分布,而备择假设是一组比另一组具有更大(或更小)的值。
与上面我们看到的其他检验不同,Mann-Whitney U 检验对异常值不可知的。
检验过程如下。
-
合并所有数据点并对它们进行排名(按升序或降序排列)
-
计算 U₁ = R₁ - n₁(n₁ + 1)/2,其中 R₁ 是第一组数据点的秩和,n₁ 是第一组数据点的数量。
-
类似地计算第二组的 U₂。
-
检验统计量由 stat = min(U₁, U₂) 给出。
在两个分布之间没有系统等级差异的原假设下(即相同的中位数),检验统计量是渐近正态分布的,具有已知的均值和方差。
计算 R 和 U 背后的理论如下:如果第一个样本中的值都大于第二个样本中的值,则 R₁ = n₁(n₁ + 1)/2 并且作为结果,U 1 将为零(可达到的最小值)。否则如果两个样本相似,U1 和 U2 将非常接近 n1 n2 / 2(可达到的最大值)。
我们使用 scipy 的 mannwhitneyu 函数。
from scipy.stats import mannwhitneyu
stat, 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
我们得到的p值为0.6,这意味着我们不拒绝实验组和对照组的中位数没有差异的零假设。
置换检验
一种非参数替代方法是置换检验。在原假设下,两个分布应该是相同的,因此打乱组标签不应该显着改变任何统计数据。
可以选择任何统计数据并检查其在原始样本中的值如何与其在组标签排列中的分布进行比较。例如使用实验组和对照组之间样本均值的差异作为检验统计。
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
置换检验得到的 p 值为 0.053,这意味着在 5% 的水平上对原假设的弱不拒绝。那么应该如何解释 p 值?这意味着数据中均值的差异大于置换样本中均值差异的 1–0.0560 = 94.4%。
我们可以通过绘制检验统计在排列中的分布与其样本值的分布来可视化。
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');
排列间的平均差分布
正如我们所看到的,样本统计量相对于置换样本中的值是相当极端的,但并不过分。
卡方检验
卡方检验是一种非常强大的检验,主要用于检验频率差异。
卡方检验最不为人知的应用之一是检验两个分布之间的相似性。这个想法是对两组的观察结果进行分类。如果两个分布相同,我们会期望每个 bin 中的观察频率相同。这里重要的一点是需要在每个 bin 中进行足够的观察,以使检验有效。
生成与对照组中收入分布的十分位数相对应的bin,然后如果两个分布相同,我计算实验组中每个bin中的预期观察数。
# Init dataframe
df_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 groups
df_bins['income_c_observed'] = pd.cut(income_c, bins=bins).value_counts().values
df_bins['income_t_observed'] = pd.cut(income_t, bins=bins).value_counts().values
# Compute expected frequency in the treatment group
df_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
bin和频率
为了计算检验统计量和检验的 p 值,我们使用 scipy 的卡方函数。
from scipy.stats import chisquare
stat, 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
与上面介绍的所有其他检验不同,卡方检验强烈拒绝两个分布相同的原假设。这是为什么?
原因在于这两个分布具有相似的中心但尾部不同,并且卡方检验测试了整个分布的相似性,而不仅仅是中心,就像我们在之前的检验中所做的那样。
这个结果讲述了一个警示:在从 p 值得出盲目结论之前,了解实际检验的内容非常重要!
Kolmogorov-Smirnov 检验
Kolmogorov-Smirnov检验的思想是比较两组的累积分布。特别是,Kolmogorov-Smirnov 检验统计量是两个累积分布之间的最大绝对差。
其中 F₁ 和 F₂ 是两个累积分布函数,x 是基础变量的值。Kolmogorov-Smirnov 检验统计量的渐近分布是 Kolmogorov 分布。
为了更好地理解,让我们绘制累积分布函数和检验统计量。首先计算累积分布函数。
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])/2
plt.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检验统计量
从图中我们可以看出,检验统计量的值对应于收入~650 时的两个累积分布之间的距离。对于该收入值在两组之间存在最大的不平衡。
我们可以使用 scipy 中的 kstest 函数执行实检验。
from scipy.stats import kstest
stat, 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% 的置信度拒绝两个分布相同的原假设。
箱线图
当我们有多组时,箱线图可以很好地扩展,因为我们可以并排放置不同的框。
sns.boxplot(x='Arm', y='Income',
data=df.sort_values('Arm'));
plt.title("Boxplot, multiple groups");
各部门的收入分配箱线图
从图中可以看出,不同实验组的收入分配不同,编号越高的组平均收入越高。
提琴图
结合汇总统计和核密度估计的箱线图的一个非常好的扩展是小提琴图。小提琴图沿 y 轴显示不同的密度,因此它们不会重叠。默认情况下,它还在里面添加了一个微型箱线图
sns.violinplot(x='Arm', y='Income',
data=df.sort_values('Arm'));
plt.title("Violin Plot, multiple groups");
各部门的收入分配提琴图
小提琴图表明不同实验组的收入不同。
山脊图
山脊图沿 x 轴绘制了多个核密度分布,它比小提琴图更直观。在 matplotlib 和 seaborn 中都没有默认的山脊线图。素以需要joypy包。
from joypy import joyplot
joyplot(df, by='Arm', column='Income',
colormap=sns.color_palette("crest", as_cmap=True));
plt.xlabel('Income');
plt.title("Ridgeline Plot, multiple groups");
各部门的收入分配山脊图
山脊图表明,编号越高的实验组收入越高。从这个图中也更容易理解分布的不同形状。
F检验
对于多个组最流行的检验方法是 F 检验。F 检验比较不同组间变量的方差。这种分析也称为方差分析。
F 检验统计量由下式给出
其中 G 是组数,N 是观察数,x̅ 是总体平均值,x̅g 是组 g 内的平均值。在组独立性的原假设下,f 统计量是 F 分布的。
from scipy.stats import f_oneway
income_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 值基本上为零,这意味着强烈拒绝实验组之间收入分配没有差异的原假设。
作者:Matteo Courthoud 来源:https://medium.com/towards-data-science 数据STUDIO 侵删
如果本文对您有帮助,希望您可以点赞+收藏+评论,您的支持是我更新的动力~
文章评论