Lab 5：数据分析进阶——模型、检验与因果推断实战

课程：数据库技术与应用：AI 时代的基石数据技能
授课教师：王健楠（计算机科学与技术系）
学期：2026 年（春季）

实验概述与学习目标

本实验涵盖从自动化特征工程到因果推断的完整方法链，共四个任务。

任务	主题	分值
任务 1	自动化特征工程	20 分
任务 2	可解释机器学习	30 分
任务 3	假设检验	25 分
任务 4	因果推断	25 分
合计		100 分

完成本实验后，你应能够：

1. 掌握自动化特征工程的基本方法： 能够使用 featuretools 对多张关联表进行深度特征合成（Deep Feature Synthesis, DFS），并通过特征选择提升模型表现。
2. 理解并实现机器学习可解释性技术： 能够从参数解读、置换重要性、部分依赖图、LIME 等多个维度解释模型行为，区分透明模型与事后解释方法的适用场景。
3. 掌握假设检验的核心思想与实现： 能够从零实现置换检验和卡方检验，理解 p 值的含义，并识别 p-hacking 的风险。
4. 理解因果推断的基本框架与估计方法： 能够绘制因果图，实现精确匹配、近邻匹配、倾向得分匹配和线性回归四种平均处理效应（ATE）估计方法。

第 0 步：环境准备

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

# 机器学习
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import MinMaxScaler
from sklearn.linear_model import LogisticRegression, LinearRegression
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import log_loss
from sklearn.inspection import PartialDependenceDisplay

# 绘图全局参数
plt.rcParams['figure.figsize'] = [14, 6]
plt.rcParams['font.size'] = 11

print("✅ 环境准备完成！")

---

任务 1：自动化特征工程（20 分）

任务背景

实际数据往往以多张关联表的形式存在，手工构造聚合特征耗时且容易遗漏有价值的组合。featuretools 提供了深度特征合成 方法：自动遍历表间关系，系统性地生成 COUNT、SUM、MEAN、MODE 等聚合特征及时间衍生特征，将繁琐的特征工程过程自动化。

本任务使用电商订单数据集（三张关联表），预测用户在未来 4 周内是否购买香蕉（二分类）。你需要构建 EntitySet、运行 DFS，并从候选特征中筛选 10 个，使随机森林 AUC 超过基线 0.61。

如尚未安装 featuretools，请运行：pip install featuretools

数据说明与加载

数据包含三张关联表：

• users：用户信息表
  • user_id：用户 ID（主键）
  • label：标签——1 表示该用户未来 4 周内购买了香蕉，0 表示没有
• orders：订单表
  • order_id：订单 ID（主键）
  • user_id：外键，关联到 users
  • order_time：下单时间
• order_products：订单商品明细表
  • order_product_id：主键
  • order_id：外键，关联到 orders
  • product_name：商品名称
  • aisle_id：货架编号
  • department：商品所属部门
  • reordered：是否为复购商品（1/0）
  • order_time：商品下单时间

import featuretools as ft

orders         = pd.read_csv("data/orders.csv")
order_products = pd.read_csv("data/order_products.csv")
users          = pd.read_csv("data/users.csv")

print(f"用户数：{len(users)}，正样本（会买香蕉）：{users['label'].sum()} 人")
print(f"订单数：{len(orders)}")
print(f"订单商品明细数：{len(order_products)}")
print()
print("标签分布（1=会买香蕉，0=不会）：")
print(users["label"].value_counts())

1.1 定义数据实体集（8 分）

featuretools 的第一步是将多张关联表组织成一个 EntitySet（实体集），明确声明每张表的主键、外键关系，以及时间列的数据类型。实体集是 DFS 的输入，只有声明正确，featuretools 才能正确推断跨表聚合的方向。

请实现函数 load_entityset(orders, order_products, users)，构建并返回一个正确的 EntitySet。

要求：

1. EntitySet 中包含三个实体：users、orders、order_products；
2. 正确声明各表的主键（user_id、order_id、order_product_id）；
3. 声明两条关系：orders.user_id → users.user_id，order_products.order_id → orders.order_id；
4. 将 order_time 列声明为时间类型（datetime）。

参考文档： featuretools 文档

def load_entityset(orders, order_products, users):
    """构建并返回包含三个实体的 EntitySet。"""
    # YOUR CODE HERE
    raise NotImplementedError()

# ===== 任务 1.1 测试 =====
es = load_entityset(orders, order_products, users)
assert es is not None, "load_entityset 返回了 None"
es_str = str(es)
for name in ["users", "orders", "order_products"]:
    assert name in es_str, f"EntitySet 中缺少实体：{name}"
print("✅ 任务 1.1 测试通过！")
print(es)

1.2 深度特征合成与特征选择（12 分）

第一步：运行 DFS，生成候选特征

请调用 ft.dfs() 以 users 为目标实体（target_dataframe_name）执行深度特征合成，将结果存入 feature_matrix 和 feature_defs。

es = load_entityset(orders, order_products, users)

# YOUR CODE HERE
# feature_matrix, feature_defs = ft.dfs(...)
raise NotImplementedError()

# 查看生成的特征矩阵（前 5 行）
print(f"共生成候选特征数：{feature_matrix.shape[1]}")
feature_matrix.head()

第二步：特征选择，击败基线

从生成的所有特征中，筛选出恰好 10 个你认为最有预测价值的特征。

要求：

1. 将选出的 10 个特征存入 X（shape 应为 (767, 10)），标签存入 y；
2. 处理特征矩阵中可能存在的缺失值（如用 0 填充）；
3. 使用下方给出的随机森林和评估代码，确保 AUC 超过 0.61。

# YOUR CODE HERE
# 选择 10 个特征，构建 X 和 y
# X = feature_matrix[selected_columns].fillna(0)
# y = users.set_index("user_id")["label"]
raise NotImplementedError()

# ===== 任务 1.2 测试 =====
assert X.shape == (767, 10), f"X 的形状应为 (767, 10)，当前为 {X.shape}"
assert len(y) == 767, "y 的长度应为 767"
print("✅ 任务 1.2 形状检验通过！")

# 模型评估（超参数已固定，不要修改）
clf = RandomForestClassifier(n_estimators=400, n_jobs=-1, random_state=42)
scores = cross_val_score(clf, X, y, cv=3, scoring="roc_auc")
auc = scores.mean()
print(f"\n3 折交叉验证 AUC = {auc:.4f}")
assert auc > 0.61, f"AUC 为 {auc:.4f}，未超过基线 0.61，请重新选择特征"
print(f"✅ 任务 1.2 测试通过！AUC = {auc:.4f}，成功超越基线 0.61")
print()
print("你选出的 10 个特征：")
for i, col in enumerate(X.columns, 1):
    print(f"  {i:2d}. {col}")

思考题（必填）

1. 请列举 featuretools 的 三个优点 和 三个缺点。
2. 针对你列出的缺点，有没有可行的改进思路？

你的回答：

（请在此处作答）

---

任务 2：可解释机器学习（30 分）

任务背景

高精度模型往往是"黑盒"——梯度提升树和神经网络的内部结构难以直接解读。可解释机器学习使模型的预测结果可供人类理解与审查，是信贷审批、医疗诊断、公共政策等高风险决策场景的必要条件。

本任务使用 UCI 森林覆盖类型（Forest Covertype） 数据集，依次实践六种解释方法：

子任务	方法	解释粒度	模型依赖	分值
2.1	透明模型（逻辑回归系数）	全局	透明模型	8
2.2	置换重要性	全局	与模型无关	8
2.3	部分依赖图（PDP）	全局	与模型无关	6
2.4	全局代理模型	全局	与模型无关	4
2.5	LIME	局部	与模型无关	4

数据集说明： 54 个特征——前 10 列为数值特征（海拔、坡度、距水源距离等），其余为独热编码的分类特征（土壤类型、荒野类型），目标变量为覆盖类型（类别 1 vs 其他）。

数据加载与模型训练

下面的代码完成数据加载、划分和三个模型的训练。请直接运行，无需修改。

# 加载数据
data = pd.read_csv("data/bforest_sample.csv")
print(f"数据集：{data.shape[0]} 条样本，{data.shape[1]} 个特征列（含标签），正样本比例：{data.iloc[:, -1].mean():.3f}")

# 划分训练集和测试集（8:2）
train, test = train_test_split(data, test_size=0.2, random_state=733)
X_train, y_train = train.iloc[:, :-1], train.iloc[:, -1]
X_test,  y_test  = test.iloc[:, :-1],  test.iloc[:, -1]
feature_names = list(X_train.columns)

# 特征缩放到 [0, 1]
scaler = MinMaxScaler()
X_train = pd.DataFrame(scaler.fit_transform(X_train), columns=feature_names)
X_test  = pd.DataFrame(scaler.transform(X_test),      columns=feature_names)

# 训练三个模型
print("正在训练模型，请稍候…")
lr  = LogisticRegression(solver='liblinear').fit(X_train, y_train)
gb  = GradientBoostingClassifier(random_state=42).fit(X_train, y_train)
mlp = MLPClassifier(learning_rate_init=0.1, random_state=42).fit(X_train, y_train)

print(f"  逻辑回归   测试集准确率：{lr.score(X_test,  y_test):.4f}")
print(f"  梯度提升树 测试集准确率：{gb.score(X_test,  y_test):.4f}")
print(f"  多层感知机 测试集准确率：{mlp.score(X_test, y_test):.4f}")

2.1 透明模型解释（9 分）

逻辑回归是一个"透明模型"——它的参数直接揭示了每个特征的贡献方向和大小，不需要任何额外的解释工具。

子任务 A——全局解释： 实现 explain_logistic_regression(lr, feature_names)，以水平柱状图展示每个特征的回归系数，正负贡献用不同颜色区分，帮助理解模型整体依赖哪些特征、方向如何。

子任务 B——局部解释： 实现 explain_logistic_regression_prediction(lr, feature_names, sample)，对某条具体样本的预测进行解释：将每个特征的系数 × 该样本的特征值，得到该特征对这次预测的实际贡献量，以水平柱状图展示。

plt.rcParams['figure.figsize'] = [16, 8]

def explain_logistic_regression(lr, feature_names):
    """以水平柱状图展示逻辑回归各特征的系数（全局解释）。"""
    # YOUR CODE HERE
    raise NotImplementedError()

def explain_logistic_regression_prediction(lr, feature_names, sample):
    """以水平柱状图展示该样本中各特征对预测的贡献（系数 × 特征值）（局部解释）。"""
    # YOUR CODE HERE
    raise NotImplementedError()

# 全局解释
explain_logistic_regression(lr, feature_names)

# 局部解释：对测试集第 0 条样本
explain_logistic_regression_prediction(lr, feature_names, (X_test.iloc[0, :], y_test.iloc[0]))

请写下你从上图中观察到的 两个 最有趣的发现。

发现

1. [请填写]
2. [请填写]

2.2 置换重要性（9 分）

对于梯度提升树和神经网络，我们无法直接读取参数来理解特征重要性。置换重要性（Permutation Importance） 提供了一种与模型结构无关的解决方案：

对测试集中某个特征的值随机打乱（置换），如果模型误差显著上升，说明该特征非常重要；如果误差变化不大，说明这个特征对模型几乎没有贡献。

重要性分数定义为：$I_j = E'_j - E_0$（打乱第 $j$ 个特征后的误差 − 原始误差）。

请实现函数 permutation_importance(model, feature_names, X, y)：

1. 计算模型在原始测试集上的对数损失 $E_0$ 作为基准；
2. 对每个特征逐一进行置换，计算置换后的误差 $E'_j$；
3. 以 $I_j = E'_j - E_0$ 作为重要性分数；
4. 绘制前 5 个最重要特征的柱状图。

注意： 不得调用 sklearn 的 permutation_importance 函数，需自行实现。

def permutation_importance(model, feature_names, X, y):
    """计算置换重要性并绘制前 5 个特征的柱状图。"""
    # YOUR CODE HERE
    raise NotImplementedError()

print("=== 梯度提升树的置换重要性 ===")
permutation_importance(gb, feature_names, X_test.to_numpy(), y_test.to_numpy())

print("\n=== 多层感知机的置换重要性 ===")
permutation_importance(mlp, feature_names, X_test.to_numpy(), y_test.to_numpy())

请写下你从上图中观察到的两个最有趣的发现。（两个模型的重要特征有何异同？）

发现

1. [请填写]
2. [请填写]

2.3 部分依赖图（6 分）

置换重要性告诉我们"哪个特征重要"，但没有告诉我们"这个特征怎么影响预测"。部分依赖图（Partial Dependence Plot, PDP） 填补了这个空白：在控制其他所有特征不变的前提下，单独展示某一特征取值变化对预测概率的边际影响。

请分别为梯度提升树和多层感知机绘制前 10 个数值特征（feature_names[:10]）的 PDP 图。

使用 sklearn.inspection.PartialDependenceDisplay.from_estimator() 实现。

plt.rcParams['figure.figsize'] = [18, 10]

# YOUR CODE HERE
# 绘制梯度提升树对前 10 个数值特征的 PDP
raise NotImplementedError()

# YOUR CODE HERE
# 绘制多层感知机对前 10 个数值特征的 PDP
raise NotImplementedError()

请写下你从上图中观察到的 两个 最有趣的发现。（对比两个模型的 PDP 形态，哪些特征呈线性关系，哪些呈非线性？）

发现

1. [请填写]
2. [请填写]

2.4 全局代理模型（4 分）

全局代理模型（Global Surrogate） 是一种间接的解释思路：

1. 用黑盒模型（梯度提升树 / 神经网络）对训练集进行预测，得到"软标签"（预测概率）；
2. 以软标签为目标，训练一个透明的逻辑回归模型；
3. 用 2.1 中实现的 explain_logistic_regression() 解读代理模型的系数。

代理模型的系数近似反映了黑盒模型的决策逻辑。

请分别对梯度提升树和多层感知机各训练一个逻辑回归代理模型，并可视化系数。

plt.rcParams['figure.figsize'] = [16, 8]

# YOUR CODE HERE
# 对 gb 和 mlp 分别训练全局代理逻辑回归，调用 explain_logistic_regression() 展示
raise NotImplementedError()

请写下你从上图中观察到的 两个 最有趣的发现。（与 2.1 中真实逻辑回归的系数相比，代理模型的系数有何异同？说明了什么？）

发现

1. [请填写]
2. [请填写]

2.5 LIME 局部解释（2 分）

LIME（Local Interpretable Model-agnostic Explanations）是另一种局部解释方法。LIME 不依赖模型内部结构，适用于任何黑盒模型。其核心思想是：

在待解释样本的附近生成大量扰动样本，用这些扰动样本的黑盒预测结果训练一个简单的线性模型，用该线性模型近似黑盒在该样本周围的行为。

请安装 lime 库，使用 lime.lime_tabular.LimeTabularExplainer 解释多层感知机对测试集第 0 条样本的预测。

import lime
import lime.lime_tabular

plt.rcParams['figure.figsize'] = [12, 5]

# YOUR CODE HERE
# 使用 LimeTabularExplainer 解释 mlp 对测试集第 0 条样本的预测
raise NotImplementedError()

思考题： 多次重新运行上方代码单元（每次 Kernel 重新执行该单元），你会发现 LIME 的解释结果每次略有不同。

1. 你如何评价 LIME 的稳定性？
2. 结合 LIME 的算法原理，请解释为什么会出现这种现象。

你的回答：

1. [请填写]
2. [请填写]

---

任务 3：假设检验（25 分）

任务背景

清华网络学堂对课程搜索框进行了 A/B 测试：A 组用户看到原有的空白搜索框，B 组看到带有提示文字的搜索框（如下图）。实验结束后，B 组用户平均搜索次数比 A 组高 0.135 次。

这一差异是界面改版的真实效果，还是随机分组带来的统计波动？假设检验提供了严格的决策框架：

1. 零假设 $H_0$：两组差异由随机分组引起，界面设计无效果。
2. 备择假设 $H_1$：差异由界面改版引起，具有统计显著性。
3. p 值：在 $H_0$ 成立的前提下，出现不低于观测差异的概率（本实验采用单尾检验，因为 $H_1$ 是有方向的——B 组搜索量高于 A 组）。
4. *决策准则：若 p < 0.05，拒绝 $H_0$，认为差异在统计上显著。

数据说明与加载

data/searchlog.json 包含清华网络学堂搜索功能 A/B 测试期间每位用户的行为记录：

字段	说明
uid	用户 ID
is_instructor	是否为教师用户（True）或学生用户（False）
search_ui	该用户被分配到的界面版本（'A' 或 'B'）
search_count	该用户在实验期间的总搜索次数

search_df = pd.read_json("data/searchlog.json", orient='records', lines=True)
print(f"实验参与用户数：{len(search_df)} 人")
print(f"A 组：{(search_df['search_ui']=='A').sum()} 人，B 组：{(search_df['search_ui']=='B').sum()} 人")
print(f"教师用户：{search_df['is_instructor'].sum()} 人，学生用户：{(~search_df['is_instructor']).sum()} 人")
search_df.head()

3.1 计算检验统计量（5 分）

在进行假设检验之前，先选定一个检验统计量来量化我们观察到的效应。

这里我们使用 B 组与 A 组搜索次数均值之差：

$$\delta = \bar{X}_B - \bar{X}_A$$

请计算 $\delta$，将结果存入变量 delta。

# 分别提取两组的搜索次数
A_search_count = search_df[search_df['search_ui'] == 'A']['search_count']
B_search_count = search_df[search_df['search_ui'] == 'B']['search_count']

# YOUR CODE HERE
# delta = ...
raise NotImplementedError()

print(f"A 组平均搜索次数：{A_search_count.mean():.4f}")
print(f"B 组平均搜索次数：{B_search_count.mean():.4f}")
print(f"均值差 δ = {delta:.6f}")

# ===== 任务 3.1 测试 =====
assert isinstance(delta, (int, float, np.floating)), "delta 应为数值类型"
assert abs(delta - 0.135) < 0.01, f"delta 计算有误，期望约 0.135，实际为 {delta:.6f}"
print("✅ 任务 3.1 测试通过！")

3.2 置换检验（10 分）

$\delta = 0.135$ 看起来不小，但真的有统计意义吗？

置换检验（Permutation Test） 的核心逻辑：如果零假设成立（两组没有真实差异），那么 A/B 的分组标签对用户来说只是一个随机标记，把标签打乱重新分配后，得到的 $\delta^*$ 应该和真实观测值差不多。

如果在成千上万次随机置换中，模拟出的 $\delta^*$ 几乎从来不会达到观测值 $\delta$ 那么大，说明 $\delta$ 在零假设下极不可能发生——p 值很小，可以拒绝零假设。

算法步骤：

1. 将 A 组和 B 组的搜索次数合并；
2. 重复 10000 次：随机打乱合并数据，按原始比例重新分成两组，计算新的均值差 $\delta^*$；
3. p 值 = $\delta^* \geq \delta_{\text{observed}}$ 的比例（单尾）。

单尾 vs 双尾： 本实验的 $H_1$ 是"B 组搜索量高于 A 组"，方向明确， 故使用单尾检验，只统计 $\delta^* \geq \delta$ 的比例。 若 $H_1$ 改为"两组存在任何差异"（无方向），则应使用双尾检验： p 值 = $|\delta^*| \geq |\delta|$ 的比例，通常比单尾 p 值更大（更保守）。

请实现函数 permutation_test(A_data, B_data, num_samples=10000)，返回 p 值。

注意： 不得调用任何现有的假设检验库函数（如 scipy.stats），需自行实现上述算法。

def permutation_test(A_data, B_data, num_samples=10000):
    """
    置换检验：计算观察到的均值差 delta 对应的 p 值。

    参数：
        A_data: A 组数据（array-like）
        B_data: B 组数据（array-like）
        num_samples: 置换次数
    返回：
        p_value (float)
    """
    # YOUR CODE HERE
    raise NotImplementedError()

np.random.seed(42)
p_value = permutation_test(A_search_count.values, B_search_count.values)
print(f"p 值 = {p_value:.4f}")
if p_value < 0.05:
    print("结论：p < 0.05，拒绝零假设——界面改版带来的搜索量提升在统计上显著。")
else:
    print("结论：p ≥ 0.05，不能拒绝零假设——观察到的差异可能只是随机波动。")

# ===== 任务 3.2 测试 =====
assert isinstance(p_value, (int, float, np.floating)), "p_value 应为数值类型"
assert 0 <= p_value <= 1, f"p_value 应在 [0, 1] 之间，实际为 {p_value}"
print("✅ 任务 3.2 测试通过！")

3.3 p-hacking 风险讨论（5 分）

若基于同一份数据集，对教师用户子群体（is_instructor == True）单独重复上述假设检验，是否存在 p-hacking 的风险？若存在，应如何避免？

你的回答： （请在此处作答，至少 100 字）

3.4 卡方检验（5 分）

卡方检验（Chi-squared Test） 用于检验两个类别变量是否统计独立。这里我们想验证：实验中 is_instructor（教师/学生）和 search_ui（A/B 分组）是否相关？

如果两者相关，说明 A/B 随机分组可能存在问题（例如教师用户被过多分入 B 组），实验结果就不可信了。

卡方统计量的计算方法：

1. 构建列联表：行为 is_instructor 的取值，列为 search_ui 的取值，单元格为频数；
2. 计算期望频数：$E_{ij} = \dfrac{\text{行和}_i \times \text{列和}_j}{\text{总计}}$；
3. 计算统计量：$\chi^2 = \sum_{i,j} \dfrac{(O_{ij} - E_{ij})^2}{E_{ij}}$；
4. 自由度 $df = (\text{行数} - 1) \times (\text{列数} - 1)$，计算 p 值。

请实现函数 chi_squared_test(df, col1, col2)，返回卡方统计量和 p 值。

注意： 不得调用 scipy.stats.chi2_contingency 等现成函数，需自行实现。可以使用 scipy.stats.chi2.sf(chi2_stat, df) 来计算 p 值。

from scipy.stats import chi2

def chi_squared_test(df, col1, col2):
    """
    卡方检验：检验 col1 和 col2 两个类别变量是否统计独立。

    参数：
        df: DataFrame
        col1, col2: 列名
    返回：
        (chi2_stat, p_value)
    """
    # YOUR CODE HERE
    raise NotImplementedError()

chi2_stat, chi2_pval = chi_squared_test(search_df, 'is_instructor', 'search_ui')
print(f"卡方统计量 χ² = {chi2_stat:.4f}")
print(f"p 值          = {chi2_pval:.4f}")
if chi2_pval < 0.05:
    print("结论：is_instructor 和 search_ui 存在显著关联——分组可能不均衡，需要关注。")
else:
    print("结论：is_instructor 和 search_ui 无显著关联——随机分组质量良好。")

# ===== 任务 3.4 测试 =====
assert isinstance(chi2_stat, (int, float, np.floating)), "chi2_stat 应为数值类型"
assert chi2_stat > 0, "卡方统计量应为正数"
assert 0 <= chi2_pval <= 1, "p 值应在 [0, 1] 之间"
print("✅ 任务 3.4 测试通过！")

请简要解释：上述卡方检验的结果对我们的 A/B 测试结论意味着什么？

你的回答： （请在此处作答）

---

任务 4：因果推断（25 分）

任务背景

相关性 ≠ 因果性是数据分析中最重要的原则之一。简单比较两组样本的结果均值，往往受混淆变量的干扰——那些同时影响干预变量与结果变量的因素。

本任务使用 Lalonde 职业培训数据集（因果推断领域的标准基准数据集，收集于 1970 年代美国），研究职业培训项目（treat）对 1978 年收入（re78）的因果效应。

核心概念——反事实（Counterfactual）： 对每个个体，我们希望同时知道其接受和不接受培训时的潜在收入，但每个人只能经历其中一条路径。平均处理效应（Average Treatment Effect, ATE） 定义为：

$$\text{ATE} = \mathbb{E}[Y(1) - Y(0)]$$

其中 $Y(1)$、$Y(0)$ 分别为接受和不接受培训时的潜在收入。本任务将实现四种控制混淆变量的 ATE 估计方法。

import pandas as pd
import numpy as np
pd.options.mode.chained_assignment = None

lalonde = pd.read_csv("data/lalonde.csv", index_col=0)

treat_group   = lalonde[lalonde['treat'] == 1]
control_group = lalonde[lalonde['treat'] == 0]

print(f"总样本量：{len(lalonde)} 人")
print(f"处理组（参加培训，NSW）：{len(treat_group)} 人，平均 re78：${treat_group['re78'].mean():.2f}")
print(f"对照组（未参加，PSID）：{len(control_group)} 人，平均 re78：${control_group['re78'].mean():.2f}")
print(f"\n朴素均值差（有偏估计）：${treat_group['re78'].mean() - control_group['re78'].mean():.2f}")
lalonde.head()

4.1 绘制因果图（5 分）

在进行 ATE 估计之前，我们首先需要通过因果图（Causal Graph） 明确变量之间的因果关系。

因果图用有向无环图（DAG）表示：节点为变量，有向箭头代表"直接影响"的方向。

在 Lalonde 数据集中：

• treat（参加培训）是干预变量；
• re78（1978 年收入）是结果变量；
• age（年龄）和 married（婚姻状况）是混淆变量——它们既影响一个人是否参加培训，也直接影响收入。

请用 graphviz 绘制这四个变量之间正确的因果图。

from graphviz import Digraph

dot = Digraph(comment='Lalonde 因果图')
dot.attr(rankdir='LR')

# 声明节点
dot.node('treat',   'treat（参加培训）')
dot.node('re78',    're78（1978年收入）')
dot.node('age',     'age（年龄）')
dot.node('married', 'married（婚姻状况）')

# YOUR CODE HERE
# 添加正确的有向边，反映变量之间的因果关系
# 提示：混淆变量应各自有两条出边
raise NotImplementedError()

dot

4.2 精确匹配（5 分）

精确匹配（Perfect Matching） 是最直接的反事实估计：对每个对照组（未参加培训）成员，在处理组中寻找所有协变量完全相同的人，以该人的收入作为"如果他也参加了培训会怎样"的估计。

协变量：age、educ、black、hispan、married、nodegree、re74、re75（共 8 个）。

请计算：（1）能进行精确匹配的对照组成员数量和比例；（2）这些成员的 ATE。若处理组中有多个匹配，取第一个。

lalonde_0 = lalonde[lalonde['treat'] == 0].reset_index()  # 对照组（PSID）
lalonde_1 = lalonde[lalonde['treat'] == 1].reset_index()  # 处理组（NSW）
covariates = ['age', 'educ', 'black', 'hispan', 'married', 'nodegree', 're74', 're75']

# YOUR CODE HERE
# 对每个对照组成员，在处理组中寻找协变量完全相同的人
# 统计 count（匹配数量）和 percentage（占所有对照组成员的比例）
raise NotImplementedError()

print(f"能进行精确匹配的对照组成员数：{count}")
print(f"占所有对照组成员的比例：{percentage:.4f}")

# YOUR CODE HERE
# 计算精确匹配的 ATE
# ATE = 平均（处理组匹配成员的 re78 - 对照组成员的 re78）
raise NotImplementedError()

# ===== 任务 4.2 测试 =====
print(f"精确匹配 ATE = ${ATE:.2f}")
assert isinstance(ATE, (int, float, np.floating)), "ATE 应为数值"
assert 500 < ATE < 3000, f"精确匹配 ATE 超出合理范围（500~3000），实际为 {ATE:.2f}"
print("✅ 任务 4.2 测试通过！")

4.3 近邻匹配（5 分）

精确匹配的覆盖率往往很低（因为现实中很难找到各方面都完全一样的人）。近邻匹配（Nearest Neighbor Matching） 放宽限制，用欧氏距离衡量相似性，为每个对照组成员找到最接近的处理组成员。

要求：

• 使用 8 个协变量，用 np.linalg.norm(x1 - x2) 计算距离；
• 距离阈值设为 1000（超过阈值视为无法匹配）；
• 统计能匹配的数量和比例，并计算 ATE。

# YOUR CODE HERE
# 实现近邻匹配（阈值 = 1000）
# 输出 count、percentage
raise NotImplementedError()

print(f"能进行近邻匹配的对照组成员数：{count}，比例：{percentage:.4f}")

# YOUR CODE HERE
# 计算近邻匹配的 ATE
raise NotImplementedError()

# ===== 任务 4.3 测试 =====
print(f"近邻匹配 ATE = ${ATE:.2f}")
assert isinstance(ATE, (int, float, np.floating)), "ATE 应为数值"
assert -500 < ATE < 2000, f"近邻匹配 ATE 超出合理范围，实际为 {ATE:.2f}"
print("✅ 任务 4.3 测试通过！")

4.4 倾向得分匹配（5 分）

在高维协变量下，欧氏距离的意义减弱（"维数灾难"）。倾向得分匹配（Propensity Score Matching, PSM） 将多维协变量"压缩"成一个标量：接受处理的条件概率 $P(\text{treat}=1 | X)$，称为倾向得分。

再在倾向得分维度上做近邻匹配，可以有效缓解维数问题。

步骤：

1. 用逻辑回归拟合 $P(\text{treat}=1 | X)$，将结果添加为 lalonde 的新列 psm；
2. 对每个对照组成员，在处理组中找倾向得分差异最小的匹配，阈值 0.01；
3. 统计匹配数量、比例，计算 ATE。

# YOUR CODE HERE
# 第一步：用逻辑回归计算倾向得分，添加 "psm" 列
raise NotImplementedError()

lalonde.head()

# ===== 任务 4.4 倾向得分测试 =====
assert 'psm' in lalonde.columns, "lalonde 中缺少 'psm' 列"
assert ((lalonde['psm'] >= 0) & (lalonde['psm'] <= 1)).all(), "倾向得分应在 [0, 1] 之间"
print("✅ 倾向得分计算通过！")

# YOUR CODE HERE
# 第二步：基于 psm 进行近邻匹配（阈值 = 0.01），计算 count、percentage 和 ATE
raise NotImplementedError()

# ===== 任务 4.4 测试 =====
print(f"PSM 匹配数：{count}，比例：{percentage:.4f}")
print(f"倾向得分匹配 ATE = ${ATE:.2f}")
assert isinstance(ATE, (int, float, np.floating)), "ATE 应为数值"
assert -500 < ATE < 2000, f"PSM ATE 超出合理范围，实际为 {ATE:.2f}"
print("✅ 任务 4.4 测试通过！")

4.5 线性回归（5 分）

线性回归方法 不进行匹配，而是分别对处理组和对照组训练两个线性回归模型：

• $M_1$：在处理组上训练，学习 $P(Y | \text{treat}=1, X)$
• $M_0$：在对照组上训练，学习 $P(Y | \text{treat}=0, X)$

然后对全体样本估计反事实，计算 ATE：

$$\text{ATE} = \frac{1}{N} \sum_{i=1}^{N} \left[ M_1(X_i) - M_0(X_i) \right]$$

# YOUR CODE HERE
# 训练两个线性回归模型，对全体样本计算反事实，估计 ATE
raise NotImplementedError()

# ===== 任务 4.5 测试 =====
print(f"线性回归 ATE = ${ATE:.2f}")
assert isinstance(ATE, (int, float, np.floating)), "ATE 应为数值"
assert -500 < ATE < 3000, f"线性回归 ATE 超出合理范围，实际为 {ATE:.2f}"
print("✅ 任务 4.5 测试通过！")

汇总与分析

请简要比较四种方法的优缺点，并解释为什么它们给出的 ATE 估计值不同。

你的回答： （请在此处作答）

---

实验圆满结束！

任务	主题	核心技术
任务 1	自动化特征工程	featuretools 深度特征合成
任务 2	可解释机器学习	系数解读、Permutation、PDP、LIME
任务 3	假设检验	置换检验、卡方检验、p-hacking 识别
任务 4	因果推断	精确匹配、近邻匹配、PSM、线性回归

四个任务依次覆盖预测 → 理解 → 验证 → 决策的数据科学能力矩阵。

---

提交前最终检查：

• 建议点击顶部菜单栏的 Kernel → Restart & Run All，确保代码环境从头到尾运行顺畅无报错。
• 确认全部任务的测试单元均已通过。
• 确认每个"发现"部分均已填写。
• 确认无误后，将完成后的notebook文件与实验所需原始 data 文件一并压缩打包后上传至网络学堂。