数据分析师必看：如何用Python实战案例拆解辛普森悖论（附完整代码）

e1f2g

466人浏览 · 2026-03-05 03:31:55

e1f2g · 2026-03-05 03:31:55 发布

数据分析师必看：如何用Python实战案例拆解辛普森悖论（附完整代码）

你是否曾在分析一份看似完美的数据报告时，发现整体结论与分组结论完全相反？那种感觉，就像你精心策划的营销活动，总点击率上升了，但每个渠道的点击率却都在下降，让人百思不得其解。这不是数据在“说谎”，而是你可能遇到了统计学中一个经典的陷阱——辛普森悖论。对于数据分析师和Python开发者而言，仅仅理解其概念是远远不够的。真正的挑战在于，如何在日常的海量数据清洗、聚合与分析流程中，敏锐地识别它，并用可靠的代码工具去验证和规避它。本文将带你跳出纯理论的窠臼，通过几个亲手可运行的Python案例，深入数据肌理，拆解这一悖论的形成机制，并为你装备一套即拿即用的分析“防坑”工具箱。

1. 从经典案例到数据生成：亲手“制造”一个悖论

理解悖论最好的方式，不是阅读它，而是创造它。当我们能用自己的代码模拟出辛普森悖论时，对其内在驱动力的把握才会真正深刻。

1.1 重温经典：药物试验数据的Python再现

让我们回到1951年那个著名的药物试验案例。与其盯着静态表格，不如用pandas和numpy将其动态构建出来，观察数据是如何“层叠”成悖论的。

首先，我们创建这个数据集。注意，这里的关键不仅在于痊愈率，更在于各组样本量的严重不均衡。

import pandas as pd
import numpy as np

# 创建数据字典，精确对应经典案例
data = {
    '组别': ['服药'] * 350 + ['未服药'] * 350,
    '性别': ['男性'] * 87 + ['女性'] * 263 + ['男性'] * 270 + ['女性'] * 80,
    '是否痊愈': [1]*81 + [0]*6 + [1]*192 + [0]*71 + [1]*234 + [0]*36 + [1]*55 + [0]*25
}

df = pd.DataFrame(data)
print("数据前10行预览：")
print(df.head(10))

运行这段代码，我们就得到了一个包含700条记录的数据框。接下来，让我们进行分层计算。这里，聚合计算的顺序和维度是揭示悖论的核心。

# 1. 计算总体痊愈率（忽略性别）
total_rate = df.groupby('组别')['是否痊愈'].mean()
print("\n总体痊愈率：")
print(total_rate.apply(lambda x: f"{x:.2%}"))

# 2. 计算按性别分层的痊愈率
stratified_rate = df.groupby(['组别', '性别'])['是否痊愈'].mean().unstack()
print("\n按性别分层的痊愈率：")
print(stratified_rate.applymap(lambda x: f"{x:.2%}"))

# 3. 计算各性别组内的样本量
group_counts = df.groupby(['组别', '性别']).size().unstack()
print("\n各分组样本数量：")
print(group_counts)

执行后，你会看到类似下面的输出：

总体痊愈率：
组别
服药      78.00%
未服药    83.00%

按性别分层的痊愈率：
性别       男性       女性
组别                    
服药    93.10%   73.00%
未服药  86.67%   68.75%

各分组样本数量：
性别     男性  女性
组别              
服药      87   263
未服药   270    80

悖论出现了：从总体看，未服药组的痊愈率（83%）高于服药组（78%）。但分别看男性和女性，服药组的痊愈率均高于未服药组。这个矛盾结果的“罪魁祸首”，就藏在最后那个样本数量表里：服药组中痊愈率较低的女性样本占了大头（263 vs 87），而未服药组中痊愈率较高的男性样本占了大头（270 vs 80）。这种样本结构的混杂，导致了聚合结果的逆转。

提示：在数据分析中，永远不要满足于一个“总计”数字。你的第一个动作应该是问：“这个总计是由哪些部分构成的？它们的权重如何？”

1.2 悖论生成器：理解混杂变量的力量

为了举一反三，我们设计一个更通用的“辛普森悖论生成器”。它可以帮助我们理解，需要满足哪些条件，数据才会“背叛”我们的直觉。

def generate_simpson_paradox(n_samples=1000, effect_strength=0.2):
    """
    生成一个模拟辛普森悖论的数据集。
    n_samples: 总样本量
    effect_strength: 处理组（如服药）在子群中的真实正向效应强度
    """
    np.random.seed(42)  # 确保可重复性

    # 模拟一个混杂变量，例如“病情严重程度”（0=轻，1=重）
    # 假设处理组（服药）更倾向于接收病情更重的病人（常见于观察性研究）
    severity = np.random.binomial(1, 0.5, n_samples)
    treatment = np.where(severity == 1,
                         np.random.binomial(1, 0.7, n_samples),  # 病情重的人，70%概率被分到处理组
                         np.random.binomial(1, 0.3, n_samples))  # 病情轻的人，30%概率被分到处理组

    # 模拟结果（例如痊愈）。基线痊愈率，病情重的人更低。
    base_rate = np.where(severity == 0, 0.8, 0.4)
    # 处理效应：处理组在各自子群内增加 effect_strength 的痊愈概率
    treatment_effect = treatment * effect_strength
    # 最终痊愈概率，并加入一些随机噪声
    prob = np.clip(base_rate + treatment_effect + np.random.normal(0, 0.05, n_samples), 0, 1)
    outcome = np.random.binomial(1, prob)

    sim_df = pd.DataFrame({
        '治疗': treatment,
        '病情严重': severity,
        '痊愈': outcome
    })
    sim_df['治疗'] = sim_df['治疗'].map({0: '对照组', 1: '处理组'})
    sim_df['病情严重'] = sim_df['病情严重'].map({0: '轻', 1: '重'})

    return sim_df

# 生成数据并验证悖论
paradox_df = generate_simpson_paradox()
print("模拟数据概览：")
print(paradox_df.groupby(['治疗', '病情严重']).agg(痊愈率=('痊愈', 'mean'), 样本数=('痊愈', 'count')))

这个生成器的核心逻辑在于，它人为地创建了一个混杂变量（病情严重程度），这个变量同时影响了患者被分配到“处理组”的概率，也直接影响了最终的“痊愈”结果。正是这种双重影响，使得在简单比较“处理组”和“对照组”时，结论被扭曲。

2. 诊断与探测：在你的数据中寻找悖论信号

当面对一份新的数据集时，我们如何系统性地筛查辛普森悖论的风险？以下是一套可操作的诊断流程。

2.1 可视化诊断：从聚合到分布的透视

可视化是发现数据异样的第一道防线。使用seaborn或matplotlib，我们可以快速绘制出数据的多层视图。

import matplotlib.pyplot as plt
import seaborn as sns

# 使用之前生成的经典案例数据 df
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 图1：总体痊愈率对比
sns.barplot(data=df, x='组别', y='是否痊愈', ax=axes[0], errorbar=None)
axes[0].set_title('总体痊愈率对比')
axes[0].set_ylabel('痊愈率')

# 图2：按性别分层的痊愈率对比
sns.barplot(data=df, x='组别', y='是否痊愈', hue='性别', ax=axes[1], errorbar=None)
axes[1].set_title('按性别分层的痊愈率')
axes[1].set_ylabel('痊愈率')

# 图3：样本量构成对比（堆叠条形图）
count_data = df.groupby(['组别', '性别']).size().unstack()
count_data.plot(kind='bar', stacked=True, ax=axes[2], color=['steelblue', 'salmon'])
axes[2].set_title('各组样本的性别构成')
axes[2].set_ylabel('样本数量')

plt.tight_layout()
plt.show()

这三张图放在一起，构成了一个强大的诊断面板：

左图告诉你“总体故事”。
中图将这个故事拆解到关键子群（这里是性别）。
右图揭示了子群样本量的分布，这是理解总体与分组结论为何矛盾的关键。

如果左图与中图的趋势相反，而右图又显示出明显的样本不均衡，那么辛普森悖论的警报就该拉响了。

2.2 量化探测：计算潜在混杂影响

除了看图，我们还需要定量的指标来评估风险。一个简单有效的方法是计算加权平均与简单平均的差异。

def check_simpson_risk(df, group_col, stratify_col, outcome_col):
    """
    检查数据集中是否存在由指定分层变量导致的辛普森悖论风险。
    返回分层调整后的率与总体率的对比。
    """
    # 计算总体平均
    overall_rate = df.groupby(group_col)[outcome_col].mean()

    # 计算分层后的加权平均（以分层组样本量为权）
    # 即先计算每个 (group_col, stratify_col) 组合内的平均，再按 stratify_col 的总体分布加权
    stratum_rate = df.groupby([group_col, stratify_col])[outcome_col].mean().unstack(level=0)
    stratum_weight = df[stratify_col].value_counts(normalize=True)

    # 确保权重与分层率的维度对齐
    adjusted_rate = (stratum_rate.T * stratum_weight).T.sum()

    result = pd.DataFrame({
        '总体率': overall_rate,
        '分层调整率': adjusted_rate,
        '差异': adjusted_rate - overall_rate
    })
    return result

# 对经典案例数据进行检测
risk_report = check_simpson_risk(df, '组别', '性别', '是否痊愈')
print("辛普森悖论风险检测报告：")
print(risk_report.applymap(lambda x: f"{x:.4f}" if isinstance(x, float) else x))

这个函数的核心思想是：如果忽略分层变量（性别）直接计算总体率，与考虑分层变量后、按照分层变量分布重新加权计算出的调整率存在显著差异，那么就说明该分层变量是一个潜在的混杂因子，数据存在辛普森悖论的风险。差异的绝对值越大，风险越高。

3. 破解之道：因果推断与分层分析技术

识别出风险只是第一步，更重要的是如何得到更可靠的结论。这里介绍两种在业务分析中实用性极强的技术。

3.1 分层分析（Stratification）与标准化（Standardization）

这是最直接的方法。既然总体结果被扭曲，那我们就在每个同质的子群内进行比较，然后再进行合理的汇总。常用的汇总方法有标准化。

def stratified_analysis(df, treatment_col, outcome_col, stratify_col):
    """
    执行分层分析，并计算标准化率差（风险差异）。
    """
    # 计算各层内的效应（例如，服药 vs 未服药的痊愈率差）
    strata_effects = {}
    strata_weights = {}
    unique_strata = df[stratify_col].unique()

    for stratum in unique_strata:
        sub_df = df[df[stratify_col] == stratum]
        # 计算该层内处理组和对照组的结果均值
        effect_by_treatment = sub_df.groupby(treatment_col)[outcome_col].mean()
        if len(effect_by_treatment) == 2: # 确保该层内两组都存在
            # 效应量：处理组率 - 对照组率
            effect = effect_by_treatment.iloc[1] - effect_by_treatment.iloc[0] # 假设第二项是处理组
            strata_effects[stratum] = effect
            # 以该层样本量作为权重
            strata_weights[stratum] = len(sub_df)

    # 计算加权平均的效应（标准化）
    total_weight = sum(strata_weights.values())
    standardized_effect = sum(effect * strata_weights[s] / total_weight for s, effect in strata_effects.items())

    # 输出结果
    result_df = pd.DataFrame({
        '层': list(strata_effects.keys()),
        '层内效应': list(strata_effects.values()),
        '层样本量': list(strata_weights.values()),
        '层权重': [w/total_weight for w in strata_weights.values()]
    })
    result_df.loc['标准化汇总'] = ['-', standardized_effect, total_weight, 1.0]

    return result_df

# 对df进行分析，假设‘服药’是处理组
# 为了方便，我们先创建一个数值型的治疗变量
df['治疗数值'] = df['组别'].map({'服药': 1, '未服药': 0})
strat_result = stratified_analysis(df, '治疗数值', '是否痊愈', '性别')
print("分层分析结果：")
print(strat_result.to_string(float_format=lambda x: f'{x:.4f}'))

通过这个分析，我们得到的“标准化汇总”效应，是在控制了性别变量后，服药相对于未服药的平均效应。它通常比简单的总体比较更能反映真实情况。

3.2 倾向得分匹配（PSM）实战简介

在观察性研究中，处理组和对照组的分配并非随机（如病情重的病人更可能服药），直接比较会有偏差。倾向得分匹配是一种模拟随机试验、控制混杂变量的高级方法。其核心思想是为每个处理组的个体，在对照组中找到一个“背景”相似的个体进行配对。

下面是一个使用statsmodels和sklearn进行简单PSM分析的示例框架：

# 注意：这是一个简化示例，真实PSM分析需要更严谨的步骤和评估。
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors

def propensity_score_matching(df, treatment_col, outcome_col, covariate_cols):
    """
    简单的倾向得分匹配实现。
    """
    # 1. 估计倾向得分（给定协变量，个体进入处理组的概率）
    X = df[covariate_cols]
    y = df[treatment_col]
    ps_model = LogisticRegression(max_iter=1000)
    ps_model.fit(X, y)
    df['倾向得分'] = ps_model.predict_proba(X)[:, 1]

    # 2. 对处理组每个样本，在对照组中寻找倾向得分最接近的样本（1:1最近邻匹配）
    treated = df[df[treatment_col] == 1]
    control = df[df[treatment_col] == 0]

    nbrs = NearestNeighbors(n_neighbors=1, metric='euclidean').fit(control[['倾向得分']])
    distances, indices = nbrs.kneighbors(treated[['倾向得分']])

    # 3. 构建匹配后的数据集
    matched_control = control.iloc[indices.flatten()].copy()
    matched_df = pd.concat([treated, matched_control])

    # 4. 比较匹配后的处理效应
    matched_effect = matched_df[matched_df[treatment_col]==1][outcome_col].mean() - \
                     matched_df[matched_df[treatment_col]==0][outcome_col].mean()

    return matched_effect, matched_df

# 示例：假设我们有包含更多协变量的数据集 `obs_df`
# covariate_cols 应包含所有可能影响治疗分配和结果的变量，如年龄、基础疾病等。
# 此处仅为流程演示
# avg_treatment_effect, matched_data = propensity_score_matching(obs_df, 'treatment', 'outcome', covariate_cols)

注意：倾向得分匹配是一个强大的工具，但实施起来需要谨慎。必须检查匹配后的数据是否在关键协变量上达到了平衡（即处理组和对照组在匹配后背景相似），并且要注意对匹配方法、匹配比例、有无放回等细节的选择。

4. 实战演练：商业分析中的悖论识别与应对

让我们将上述技术应用到一个更贴近业务的场景中：评估一个网站改版（如按钮颜色从蓝色改为绿色）对用户点击率的影响。

4.1 场景构建与问题定义

假设我们有A/B测试数据，但流量分配不是完全均匀的：新版页面更多地推给了新用户，而老版页面更多地留给了老用户。已知新用户的点击率天然低于老用户。原始数据可能如下表所示：

用户类型	版本	访问人数	点击人数	点击率
新用户	新版	10,000	400	4.00%
新用户	旧版	2,000	60	3.00%
老用户	新版	3,000	300	10.00%
老用户	旧版	15,000	1350	9.00%

总体计算：

新版总点击率 = (400+300) / (10000+3000) ≈ 5.38%
旧版总点击率 = (60+1350) / (2000+15000) ≈ 8.29%

结论：旧版点击率更高，改版失败？

分层计算：

在新用户中，新版点击率（4.00%）> 旧版（3.00%）。
在老用户中，新版点击率（10.00%）> 旧版（9.00%）。

结论：在每一个用户群体中，新版表现都更好！

4.2 Python代码分析与决策

我们用Python来复现这个分析，并做出正确决策。

# 构建业务数据
business_data = pd.DataFrame({
    '用户类型': ['新用户', '新用户', '老用户', '老用户'],
    '版本': ['新版', '旧版', '新版', '旧版'],
    '访问人数': [10000, 2000, 3000, 15000],
    '点击人数': [400, 60, 300, 1350]
})
business_data['点击率'] = business_data['点击人数'] / business_data['访问人数']

print("业务数据明细：")
print(business_data.to_string(index=False))

# 计算总体点击率（错误的方法）
total_clicks = business_data.groupby('版本').agg(总访问=('访问人数', 'sum'), 总点击=('点击人数', 'sum'))
total_clicks['总体点击率'] = total_clicks['总点击'] / total_clicks['总访问']
print("\n**忽略用户类型的总体比较（可能误导）：**")
print(total_clicks[['总体点击率']].applymap(lambda x: f"{x:.2%}"))

# 计算分层点击率
print("\n**按用户类型分层的点击率（正确比较）：**")
print(business_data.pivot(index='用户类型', columns='版本', values='点击率').applymap(lambda x: f"{x:.2%}"))

# 计算标准化点击率（以用户类型整体分布为权）
user_type_distribution = business_data.groupby('用户类型')['访问人数'].sum() / business_data['访问人数'].sum()
print(f"\n用户类型总体分布：\n{user_type_distribution.apply(lambda x: f'{x:.2%}')}")

# 计算新版和旧版的标准化点击率
def calculate_standardized_rate(version):
    # 对于指定版本，计算其在各用户类型的点击率，然后按总体用户分布加权
    version_data = business_data[business_data['版本'] == version].set_index('用户类型')['点击率']
    # 确保顺序与权重分布一致
    standardized_rate = (version_data * user_type_distribution).sum()
    return standardized_rate

std_rate_new = calculate_standardized_rate('新版')
std_rate_old = calculate_standardized_rate('旧版')
print(f"\n**标准化后的点击率比较（控制用户类型结构）：**")
print(f"新版标准化点击率：{std_rate_new:.2%}")
print(f"旧版标准化点击率：{std_rate_old:.2%}")
print(f"新版相对于旧版的提升：{(std_rate_new/std_rate_old - 1):.2%}")

运行这段代码，你会清晰地看到，在控制了“用户类型”这个混杂变量后（通过标准化），新版的点击率其实是高于旧版的。正确的商业决策应该是推进新版上线，而不是被总体数据误导而回滚。

这个案例给你的核心启发是：在评估任何策略效果时，尤其是当实验组和对照组存在系统性差异时（如本例中新老用户比例不同），必须进行分层分析或因果推断，剥离出策略本身的净效应。

在日常工作中，养成一个习惯：看到任何“总体提升”或“总体下降”，立即思考并检查是否存在重要的子群体，以及这些子群体的样本构成。你的分析工具箱里，pandas的groupby、pivot_table，以及上面演示的标准化函数，应该成为条件反射般的首选武器。数据不会主动欺骗你，但忽略其内在结构，它就会给你讲一个完全相反的故事。