看想看的做想做的电影网站东莞社保官方网站

看想看的做想做的电影网站,东莞社保官方网站,网站qq代码生成,小型企业门户网站源码本项目旨在详细解释一系列用于因果推断的Python代码实现。内容面向有一定机器学习基础#xff0c;但对因果推断领域尚不熟悉的读者。我们将深入探讨代码的每个阶段#xff0c;解释其背后的原理、目的和具体实现细节。通过本项目#xff0c;将能够理解这些复杂的因果模型是如…本项目旨在详细解释一系列用于因果推断的Python代码实现。内容面向有一定机器学习基础但对因果推断领域尚不熟悉的读者。我们将深入探讨代码的每个阶段解释其背后的原理、目的和具体实现细节。通过本项目将能够理解这些复杂的因果模型是如何工作的并具备将其应用于自己数据集的能力。核心概念什么是因果推断在机器学习中我们常常关注预测Prediction给定特征X预测结果Y。例如根据用户的浏览历史预测他是否会点击广告。然而在许多现实场景中我们更关心干预Intervention的效果如果我们改变某个变量TTreatment处理变量结果Y会如何变化这就是因果推断的核心问题。例如给病人用新药T1相比于用安慰剂T0病人的康复率Y会提高多少实施一项新的教学方法T1相比于传统方法T0学生的考试成绩Y会提升多少调整产品价格T会对销售量Y产生多大的影响挑战混淆Confounding直接比较接受处理T1和未接受处理T0的群体的结果Y通常是错误的。为什么因为决定个体是否接受处理的因素通常称为混淆变量 W, Confounders可能同时也会影响结果Y。例如在药物试验中如果医生倾向于给病情更严重的病人使用新药那么即使新药有效直接比较用药组和未用药组的康复率也可能会低估药效因为用药组的病人本身就“底子差”。这里的“病情严重程度”就是一个混淆变量W。因果推断的目标就是在使用观测数据而非严格的随机对照试验数据时通过各种统计和机器学习方法剥离混淆变量W的影响从而估计出处理T对结果Y的真实因果效应。异质性处理效应 (Heterogeneous Treatment Effects, HTE)同一个处理如新药、新政策对不同个体的影响可能不同。例如新药可能对年轻人更有效对老年人效果较差。这种因果效应随个体特征协变量 X, Covariates/Features变化的现象称为异质性处理效应。估计HTE对于个性化决策如精准医疗、个性化营销至关重要。本教程将详细介绍几种先进的模型来处理混淆问题并估计异质性因果效应。第一阶段环境配置与库导入目的:准备运行代码所需的环境导入所有必需的Python库。讲解:这一步是任何Python项目的标准起始点。我们导入了warnings: 用于控制警告信息的显示simplefilter(ignore)可以隐藏一些库版本更新或不推荐用法产生的警告让输出更整洁。numpy(np): Python中科学计算的基础库用于高效处理数组和矩阵。pandas(pd): 提供强大的数据结构如DataFrame方便进行数据处理和分析。matplotlib.pyplot(plt) 和seaborn(sns): 流行的Python数据可视化库用于绘制各种图表。sklearn: Scikit-learn库包含大量的机器学习算法如线性模型Lasso、LogisticRegression聚类KMeans降维PCA等和工具。keras: 一个高级神经网络API这里用于构建DeepIV模型中的神经网络部分。econml: 微软开发的因果推断库包含了本教程核心的多种先进因果模型如DeepIV,DROrthoForest,CausalForestDML等。这些模型专门设计用来处理混淆问题并估计因果效应。其他辅助库如datetime,time,itertools。# 导入基础库并忽略警告信息import warningswarnings.simplefilter(ignore) # 忽略运行中可能出现的警告信息保持输出干净import datetimeimport time# 导入数据处理和可视化相关库import numpy as np # 用于数值计算特别是数组操作from itertools import product # 用于生成笛卡尔积方便创建测试数据集from sklearn.linear_model import Lasso, LogisticRegression # 导入线性模型Lasso回归和逻辑回归import matplotlib.pyplot as plt # 用于绘制图表import pandas as pd # 用于数据处理特别是DataFrame结构# 导入因果推断相关库 (EconML)import econml # 导入EconML主库# 从EconML的orf(Orthogonal Random Forest)模块导入两种基于森林的因果模型from econml.orf import DMLOrthoForest, DROrthoForest# 从EconML的dml(Double Machine Learning)模块导入因果森林模型from econml.dml import CausalForestDML# 从EconML的sklearn扩展中导入带权重的Lasso模型用于某些因果模型的内部估计from econml.sklearn_extensions.linear_model import WeightedLassoCVWrapper, WeightedLasso# 从EconML的dr(Double Robust)模块导入线性双重稳健学习器 (虽然这里没直接用但属于常用工具)from econml.dr import LinearDRLearner# 从EconML的iv(Instrumental Variable)的nnet(Neural Network)模块导入DeepIV模型from econml.iv.nnet import DeepIVimport keras # 导入Keras深度学习库用于构建DeepIV的神经网络部分import matplotlib.pyplot as plt # 再次导入确保可用# 导入其他辅助库import seaborn as sns # 另一个强大的可视化库基于matplotlibfrom sklearn.linear_model import LassoCV # 带交叉验证的Lasso模型from sklearn.ensemble import GradientBoostingRegressor # 梯度提升回归模型可用作内部估计器from econml.inference import BootstrapInference # 用于进行自助法推断如计算置信区间# 再次忽略警告 (确保后续代码块也不会有警告干扰)import warningswarnings.simplefilter(ignore)第二阶段数据生成目的:创建一个模拟数据集其中我们人为设定了真实的因果关系和混淆机制。这使得我们能够评估后续因果模型的效果因为我们可以将模型估计出的因果效应与我们设定的“真实”效应进行比较。讲解:生成模拟数据的步骤如下定义真实处理效应函数 (exp_te): 我们设定真实的个体处理效应Treatment Effect, TE仅依赖于协变量X的第一个维度x[0]具体形式为exp(2*x[0])。这意味着处理效应是异质性的随着x[0]的增大而指数级增大。设置参数: 定义样本量n、控制变量W的维度n_w、协变量X的维度n_x等。support_size定义了W中实际影响结果Y和处理T的变量数量引入了稀疏性。生成变量TE * T: 处理的真实因果效应只有当T1时才产生影响。np.dot(W[:, support_Y], coefs_Y): 控制变量W对结果Y的直接影响混淆路径的一部分。epsilon_sample(n): 无法观测的随机噪声。 这个生成公式Y TE * T f(W) error是因果推断中典型的结果生成模型。W(控制变量/混淆变量): 从标准正态分布生成。这些变量会同时影响处理决策T和结果Y从而产生混淆。X(协变量/异质性来源): 从均匀分布生成。这些变量决定了个体处理效应的大小。instrument(Z, 工具变量): 从均匀分布生成。这个变量将在DeepIV模型中使用。工具变量需要满足特定条件与T相关与Y仅通过T相关这里我们直接生成一个假设它满足条件。注意在真实应用中找到有效的工具变量是极具挑战性的。TE(真实处理效应): 根据上面定义的exp_te函数和生成的X计算每个个体的真实处理效应。T(处理变量): 生成一个二元0或1处理变量。其生成过程依赖于控制变量Wnp.dot(W[:, support_T], coefs_T)) 和一些随机噪声eta_sample。通过sigmoid函数将线性组合转换为概率再用二项分布抽样得到0或1。这种依赖性引入了混淆W高的个体可能更倾向于接受处理或不接受取决于系数。Y(结果变量): 结果Y由三部分构成生成测试数据 (X_test): 创建一个覆盖X取值范围这里是0到1的等间隔数据点网格。我们将用这些点来评估模型估计出的处理效应函数CATE(x)的准确性。# 定义处理效应函数 (Ground Truth)def exp_te(x):定义真实的个体处理效应 (Individual Treatment Effect)这个效应只依赖于协变量X的第一个维度 x[0]。这是一个指数函数表示处理效果随 x[0] 增大而增强。return np.exp(2 * x[0])# 设置随机种子和数据集大小参数np.random.seed(1234) # 固定随机种子使得每次运行结果都一样方便复现n 1000 # 生成的样本数量n_w 30 # 控制变量 W (Confounders / Controls) 的维度support_size 5 # W 中实际影响 Y 和 T 的变量个数 (引入稀疏性)n_x 1 # 协变量 X (Heterogeneity source) 的维度# 设置结果变量 Y 的生成参数# 从 W 的 n_w 个维度中随机选出 support_size 个维度作为真正影响 Y 的变量support_Y np.random.choice(range(n_w), sizesupport_size, replaceFalse)# 为这些影响 Y 的 W 变量生成随机系数 (0到1之间)coefs_Y np.random.uniform(0, 1, sizesupport_size)def epsilon_sample(n):生成结果变量 Y 中的随机噪声项return np.random.uniform(-1, 1, sizen)# 设置处理变量 T 的生成参数support_T support_Y # 在这个模拟中假设影响 Y 和影响 T 的 W 变量是同一组 (制造混淆)# 为这些影响 T 的 W 变量生成随机系数 (0到1之间)coefs_T np.random.uniform(0, 1, sizesupport_size)def eta_sample(n):生成处理变量 T 决策过程中的随机噪声项return np.random.uniform(-1, 1, sizen)# --- 开始生成核心数据 ---# 生成控制变量 W: n个样本n_w个特征服从标准正态分布W np.random.normal(0, 1, size(n, n_w))# 生成协变量 X: n个样本n_x个特征服从[0, 1]均匀分布X np.random.uniform(0, 1, size(n, n_x))# 生成工具变量 Z (Instrumental Variable): n个样本n_x个特征服从[0, 1]均匀分布# 假设这个Z满足工具变量的条件 (与T相关与Y只通过T相关)instrument np.random.uniform(0, 1, size(n, n_x))# 生成真实的异质性处理效应 (Heterogeneous Treatment Effect, HTE)# 对每个样本的 X应用 exp_te 函数计算其真实的 TETE np.array([exp_te(x_i) for x_i in X])# 生成处理变量 T (二元变量, 0 或 1)# 1. 计算 W 对处理决策的线性影响 噪声log_odds np.dot(W[:, support_T], coefs_T) eta_sample(n)# 2. 使用 sigmoid 函数将 log-odds 转换为概率 P(T1 | W)T_sigmoid 1 / (1 np.exp(-log_odds))# 3. 根据计算出的概率为每个样本随机生成 T1 或 T0T np.array([np.random.binomial(1, p) for p in T_sigmoid])# 注意T 的生成依赖于 W这是引入混淆的关键步骤# 生成结果变量 Y# Y (真实的TE * 是否接受处理T) W对Y的直接影响 随机噪声Y TE * T np.dot(W[:, support_Y], coefs_Y) epsilon_sample(n)# 注意Y 同时受到 T 和 W 的影响# 设置后续模型 (如正交森林) 的参数和测试数据subsample_ratio 0.4 # 子采样率森林模型中常用参数# 生成一个测试点集 X_test覆盖 X 的取值范围 (0到1步长0.01)# 用于评估模型估计出的处理效应函数 CATE(x)X_test np.array(list(product(np.arange(0, 1, 0.01), repeatn_x)))第三阶段构建和训练DeepIV模型目的:使用DeepIV深度工具变量模型来估计处理效应。DeepIV特别适用于存在内生性Endogeneity问题的情况即处理变量T与结果模型Y的误差项相关通常由未观测到的混淆变量引起并且我们能找到一个有效的工具变量Z。讲解:DeepIV是一种基于神经网络的两阶段方法第一阶段 (Treatment Model):目的利用工具变量Z和控制变量W来预测处理变量T的分布而不是仅仅预测T的值。工具变量Z被认为是外生的不受未观测混淆的影响它可以帮助我们识别T中那些不受内生性污染的变动。实现这里使用一个 Keras 神经网络treatment_model来完成。输入是Z和W的拼接。模型的输出不是直接预测T而是被内部用于估计T的条件分布通过混合密度网络或类似技术DeepIV内部实现。第二阶段 (Outcome Model):目的基于第一阶段对T分布的估计以及控制变量W来预测结果Y。关键在于它使用的不是观测到的可能有内生性的T而是利用第一阶段模型从Z和W生成采样的、理论上与误差项无关的T的预测值/分布信息。实现使用另一个 Keras 神经网络outcome_model。输入是采样得到的T和W的拼接。DeepIV整合与训练:DeepIV类将这两个神经网络模型 (m对应第一阶段h对应第二阶段) 结合起来。n_components: 用于近似T的条件分布的混合高斯模型GMM的组件数。optimizer: 指定训练神经网络的优化器这里用Adagrad。keras_fit_options: 分别为两个阶段的神经网络训练指定参数如训练轮数(epochs)、验证集比例(validation_split)、早停(EarlyStopping防止过拟合)。fit(Y, T, XW, Zinstrument): 使用数据训练整个DeepIV模型。注意这里参数X传入的是我们的控制变量W而工具变量通过Z参数传入。预测与可视化:第一张图将预测值 (cate) 和真实值 (true y) 按预测值大小排序并计算滚动平均 (rolling(window15).mean())。这有助于观察预测值和真实值之间的整体趋势关系是否一致平滑掉个体噪声。第二张图直接绘制排序后的原始预测值和真实值展示更细节的对比。比较预测值均值和真实值均值。predict(T, W): 使用训练好的模型进行预测。这里预测的是给定T和W下的期望结果E[Y|T,W]。由于模型的目标是估计因果效应这里的预测值可以解释为条件平均处理效应 (Conditional Average Treatment Effect, CATE)的一种估计。注意DeepIV的predict方法通常用于估计E[Y|do(T), W]或相关的量具体解释可能需要参考EconML文档但在这里我们将其视为对Y的预测并与真实的Y比较。可视化代码绘制了两张图。DeepIV的适用性:当你怀疑存在未观测的混淆因素并且能找到一个强而有效的工具变量时DeepIV是一个强大的选择。但在实践中找到好的工具变量非常困难。# --- 构建神经网络模型 ---# 定义处理模型 (对应 DeepIV 的第一阶段)treatment_model keras.Sequential([# 输入层: 维度是 W 的维度 Z 的维度 (这里 Z 和 X 维度相同都是 n_x1)# 使用 ReLU 激活函数添加 L2 正则化防止过拟合keras.layers.Dense(128, activationrelu, input_shape(n_w n_x,), # W.shape[1]是n_w, Z的维度是n_xkernel_regularizerkeras.regularizers.l2(0.01)),keras.layers.Dropout(0.17), # Dropout 层随机失活一部分神经元增强泛化能力keras.layers.Dense(64, activationrelu,kernel_regularizerkeras.regularizers.l2(0.01)), # 隐藏层 1keras.layers.Dropout(0.17),keras.layers.Dense(32, activationrelu,kernel_regularizerkeras.regularizers.l2(0.01)), # 隐藏层 2keras.layers.Dropout(0.17),keras.layers.Dense(1) # 输出层 (DeepIV内部会处理这个输出来估计T的分布)], nametreatment_model) # 给模型命名# 定义结果模型 (对应 DeepIV 的第二阶段)outcome_model keras.Sequential([# 输入层: 维度是 T 的维度 (1) W 的维度 (n_w)keras.layers.Dense(128, activationrelu, input_shape(1 n_w,), # T.shape[1]是1 (或没有默认1), W.shape[1]是n_wkernel_regularizerkeras.regularizers.l2(0.01)),keras.layers.Dropout(0.17),keras.layers.Dense(64, activationrelu,kernel_regularizerkeras.regularizers.l2(0.01)), # 隐藏层 1keras.layers.Dropout(0.17),keras.layers.Dense(32, activationrelu,kernel_regularizerkeras.regularizers.l2(0.01)), # 隐藏层 2keras.layers.Dropout(0.17),keras.layers.Dense(1) # 输出层预测结果 Y], nameoutcome_model) # 给模型命名# --- 设置训练参数 ---# 第一阶段 (处理模型) 训练参数keras_fit_options_1 {epochs: 50, # 训练轮数validation_split: 0.20, # 从训练数据中分出 20% 作为验证集# 使用回调函数EarlyStopping - 如果验证集损失在 2 轮内没有改善则停止训练并恢复最佳权重callbacks: [keras.callbacks.EarlyStopping(patience2, restore_best_weightsTrue)]}# 第二阶段 (结果模型) 训练参数 (这里用了和第一阶段不同的options变量名但内容相同可以复用)keras_fit_options_2 {epochs: 50,validation_split: 0.20,callbacks: [keras.callbacks.EarlyStopping(patience2, restore_best_weightsTrue)]}# --- 创建并训练 DeepIV 模型 ---# 创建 DeepIV 估计器实例deepIvEst DeepIV(n_components15, # 混合密度网络 (MDN) 或 GMM 的组件数用于估计 T 的条件分布# m: 第一阶段模型函数接收 Z 和 X (这里是W)输出 T 的分布参数# keras.layers.concatenate 用于拼接 Z 和 W 作为输入mlambda z, x: treatment_model(keras.layers.concatenate([z, x])),# h: 第二阶段模型函数接收 T (采样得到) 和 X (这里是W)预测 Yhlambda t, x: outcome_model(keras.layers.concatenate([t, x])),n_samples2, # 从 T 的条件分布中采样的数量 (较小以加速)use_upper_bound_lossTrue, # 是否使用特定的损失函数上界 (有时能稳定训练)n_gradient_samples0, # 用于梯度计算的样本数 (0表示使用解析梯度如果可用)optimizerkeras.optimizers.Adagrad(learning_rate0.0009), # 指定优化器和学习率first_stage_optionskeras_fit_options_1, # 传递第一阶段训练参数second_stage_optionskeras_fit_options_2 # 传递第二阶段训练参数)# 训练 DeepIV 模型# Y: 结果变量# T: 处理变量# X: 这里传入的是控制变量 W (在DeepIV的语境下X常指 conditioning set)# Z: 工具变量 instrumentprint(Starting DeepIV training...)deepIvEst.fit(YY, TT, XW, Zinstrument)print(DeepIV training finished.)# --- 预测并可视化结果 ---# 使用训练好的模型进行预测# 注意这里的 T 和 W 是原始观测数据y_pred_per deepIvEst.predict(T, W)print(DeepIV prediction finished.)# 将预测结果和真实Y放入DataFrame方便处理和绘图te_df_per pd.DataFrame({cate: y_pred_per.flatten(), true y: Y.flatten()}) # 确保是一维数组# 按预测的 CATE 值排序te_df_per.sort_values(cate, inplaceTrue, ascendingTrue)te_df_per.reset_index(inplaceTrue, dropTrue) # 重置索引# 计算滚动平均 (窗口大小15)用于平滑曲线观察趋势z_per te_df_per.rolling(window15, centerTrue).mean()# 计算预测值和真实值的均值cate_mean np.mean(te_df_per[cate])y_mean np.mean(te_df_per[true y])# --- 绘制第一个图滚动平均后的 CATE vs True Y ---plt.figure(figsize(10, 6))plt.plot(z_per[cate],marker., linestyle-,linewidth0.8, markersize3,labelfPredicted Y (smoothed, mean: {cate_mean:.2f})) # 使用f-string格式化标签plt.plot(z_per[true y],marker., linestyle-,linewidth0.8, markersize3,labelfTrue Y (smoothed, mean: {y_mean:.2f}))plt.ylabel(Outcome (Y) / Predicted Outcome)plt.xlabel(Observations sorted by Predicted Outcome)plt.title(DeepIV: Smoothed Predicted vs True Outcome (Sorted))plt.legend()plt.grid(True, linestyle--, alpha0.6)plt.show()# --- 绘制第二个图原始排序后的 CATE vs True Y ---plt.figure(figsize(10, 6))# 先画 True Yplt.plot(te_df_per[true y],marker., linestyle, # 只画点不连线markersize2, alpha0.7,labelfTrue Y (mean: {y_mean:.2f}))# 再画 CATEplt.plot(te_df_per[cate],marker., linestyle, # 只画点markersize2, alpha0.7, colorblack,labelfPredicted Y (mean: {cate_mean:.2f}))plt.ylabel(Outcome (Y) / Predicted Outcome)plt.xlabel(Observations sorted by Predicted Outcome)plt.title(DeepIV: Raw Predicted vs True Outcome (Sorted))plt.legend()plt.grid(True, linestyle--, alpha0.6)plt.show()第四阶段DROrthoForest和CausalForestDML进行因果推断目的:使用基于双重机器学习 (Double Machine Learning, DML)和随机森林思想的先进模型来估计异质性处理效应 (HTE)。这些方法通常对混淆有较好的鲁棒性并且不需要工具变量尽管它们可以结合使用。讲解:双重机器学习 (DML) 的核心思想:DML 旨在减少估计因果效应时的偏差这种偏差可能源于使用复杂机器学习模型来拟合所谓的“讨厌的”参数nuisance parameters即结果模型E[Y|X, W]和倾向得分模型P(T1|X, W)。DML 通过以下步骤实现交叉拟合 (Cross-fitting):将数据分成几折例如K折。拟合讨厌模型:对每一折 i使用其他 K-1 折的数据训练结果模型和倾向得分模型。计算残差:使用在其他数据上训练好的模型在第 i 折上计算残差。结果残差大致为Y - E[Y|X, W]处理残差大致为T - P(T1|X, W)。最终阶段估计:在所有折的残差上使用一个相对简单的模型如岭回归、Lasso或另一个森林来估计处理残差对结果残差的影响。这个影响就是我们想要的 CATE。 这种交叉拟合的方式确保了用于计算残差的模型没有在当前数据点上“偷看”结果从而减少了过拟合偏差并使得理论性质更好如半参数有效性。本节使用的模型:DROrthoForest(Double Robust Orthogonal Random Forest):n_trees,min_leaf_size,max_depth,subsample_ratio: 标准的随机森林参数。propensity_model: 用于估计倾向得分P(T1|X, W)的模型这里用带L1惩罚的逻辑回归。model_Y: 用于估计结果模型E[Y|X, W]的模型这里用Lasso。propensity_model_final,model_Y_final: DML最后阶段用于估计CATE的模型这里用带权重的Lasso。fit(Y, T, XX, WW): 训练模型。注意这里X是异质性来源W是控制变量。effect(X_test): 估计在测试点X_test上的 CATE。effect_interval(X_test): 计算 CATE 的置信区间通常使用自助法 Bootstrapping of Little Bags, BLB。结合了 DML、正交随机森林Orthogonal Random Forest和双重稳健性Double Robustness的思想。正交化旨在使处理效应的估计对讨厌模型的微小误差不敏感。双重稳健性意味着只要结果模型或倾向得分模型中至少有一个被正确设定或估计得足够好那么处理效应的估计就是一致的渐近无偏。这是一个非常好的性质。CausalForestDML(Causal Forest Double Machine Learning):model_y,model_t: 分别指定内部用于结果模型和处理倾向得分模型的估计器这里用Lasso和逻辑回归。n_estimators,min_samples_leaf,max_depth,max_samples: 森林的参数。discrete_treatmentTrue: 指明处理变量是离散的0/1。fit(Y, T, XX, WW): 训练。effect(X_test),effect_interval(X_test): 功能同上。summary(): 提供模型的摘要信息可能包括平均处理效应(ATE)的估计等。也是基于 DML 和森林的方法实现上可能略有不同但目标相似估计 HTE。可视化比较:代码绘制了一个包含两个子图的图形。左图展示DROrthoForest估计的 CATE(x) 曲线、真实的TE(x)曲线以及 95% 置信区间。右图展示CausalForestDML的类似结果。通过比较模型曲线与真实曲线的接近程度以及置信区间的宽度可以评估模型的性能。第二个CausalForestDML实例与推断:创建了另一个CausalForestDML实例这次设置了不同的参数如cv交叉验证折数discrete_treatmentFalse虽然我们的T是二元的但有时设为False也能工作可能处理方式不同。const_marginal_effect_inference(X_test): 估计平均边际效应并进行推断点估计置信区间。这相当于假设效应是恒定的估计一个平均值。effect_inference(X_test, T00, T11): 专门估计从 T0 变化到 T1 的处理效应并提供推断点估计置信区间。这对于二元处理变量来说就是估计 CATE。结果被绘制出来。summary_frame(): 将推断结果点估计、标准误、置信区间、p值以 DataFrame 格式返回。构建特征矩阵 (total_frame):CATE 的 99% 置信区间下界 (lb)CATE 的 99% 置信区间上界 (ub)CATE 的点估计 (point)协变量X控制变量W为了进行后续的聚类分析代码将以下信息合并到一个 NumPy 数组total_frame中这个矩阵包含了每个样本的效应估计信息以及它们的特征是聚类分析的基础。# --- DROrthoForest ---# 设置子采样率和正则化参数 lambda_regsubsample_ratio .3# 计算 Lasso 和 Logistic Regression 的 L1 正则化系数 lambda# 这个公式基于理论目的是平衡模型复杂度和拟合优度避免过拟合# 它考虑了特征数量 (W.shape[1]) 和样本量 (W.shape[0])lambda_reg np.sqrt(np.log(W.shape[1]) / (10 * subsample_ratio * W.shape[0]))# 创建 DROrthoForest 估计器实例est DROrthoForest(n_trees100, # 森林中树的数量 (旧版参数名可能是 n_estimators)min_leaf_size10, # 每片叶子节点最少的样本数max_depth20, # 每棵树的最大深度subsample_ratiosubsample_ratio, # 构建每棵树时使用的样本比例# 倾向得分模型: 使用带 L1 惩罚的逻辑回归 (solversaga 支持 L1)# C 是正则化强度的倒数这里根据 lambda_reg 计算propensity_modelLogisticRegression(C1 / (X.shape[0] * lambda_reg), penaltyl1, solversaga),# 结果模型: 使用 Lasso 回归 (L1 惩罚)model_YLasso(alphalambda_reg),# DML 最后阶段的倾向得分模型 (可以与初始模型相同或不同)propensity_model_finalLogisticRegression(C1 / (X.shape[0] * lambda_reg), penaltyl1, solversaga),# DML 最后阶段的结果模型 (使用带权重的 Lasso)model_Y_finalWeightedLasso(alphalambda_reg),n_jobs1 # 使用的 CPU 核心数 (1 表示不并行避免某些环境下的问题))# 训练 DROrthoForest 模型print(Starting DROrthoForest training...)# Y: 结果变量, T: 处理变量, X: 异质性来源, W: 控制变量est.fit(Y, T, XX, WW)print(DROrthoForest training finished.)# 计算在测试点 X_test 上的处理效应 (CATE)treatment_effects est.effect(X_test)# 计算处理效应的 95% 置信区间 (使用 BLB 方法)te_lower, te_upper est.effect_interval(X_test, alpha0.05) # alpha0.05 表示 95% 置信度print(DROrthoForest effect estimation finished.)# --- CausalForestDML ---# 创建 CausalForestDML 估计器实例est2 CausalForestDML(model_yLasso(alphalambda_reg), # 结果模型 (Y ~ W X)model_tLogisticRegression(C1 / (X.shape[0] * lambda_reg), penaltyl1, solversaga), # 处理模型 (T ~ W X), 使用 penaltyl1 和 saga solvern_estimators200, # 森林中树的数量min_samples_leaf5, # 叶子节点最小样本数max_depth50, # 树的最大深度max_samplessubsample_ratio / 2, # 构建每棵树时抽取的样本比例 (这里设得比 DROrthoForest 小)discrete_treatmentTrue, # 指明 T 是离散变量 (0/1)random_state123, # 随机种子保证结果可复现n_jobs1 # 使用 CPU 核心数)# 训练 CausalForestDML 模型print(Starting CausalForestDML training...)# cache_valuesTrue 可以缓存中间计算结果加速后续 effect 等方法的调用est2.fit(Y, T, XX, WW, cache_valuesTrue)print(CausalForestDML training finished.)# 计算在测试点 X_test 上的处理效应 (CATE)treatment_effects2 est2.effect(X_test)# 计算 CATE 的 95% 置信区间te_lower2, te_upper2 est2.effect_interval(X_test, alpha0.05)print(CausalForestDML effect estimation finished.)# 打印第二个模型 (CausalForestDML) 的摘要信息# 可能包含估计的平均处理效应 (ATE) 及其标准误等print(\nCausalForestDML Summary:)try:print(est2.summary())except AttributeError:print(Summary method not available for this version or configuration.)# --- 可视化比较两个森林模型的结果 ---plt.figure(figsize(15, 6)) # 创建一个宽一点的图形容纳两个子图# 计算真实的 TE(x) 用于比较expected_te np.array([exp_te(x_i) for x_i in X_test])# 左侧子图DROrthoForest 结果plt.subplot(1, 2, 1) # 1行2列的第1个plt.title(DROrthoForest CATE Estimation)# 绘制 DROrthoForest 估计的 CATE 曲线plt.plot(X_test[:, 0], treatment_effects, labelDROrthoForest CATE Estimate, colortab:orange)# 绘制真实的 TE 曲线plt.plot(X_test[:, 0], expected_te, b--, labelTrue Effect (Ground Truth))# 填充置信区间区域plt.fill_between(X_test[:, 0], te_lower, te_upper, colortab:orange, alpha0.3, label95% Confidence Interval (BLB))plt.ylabel(Treatment Effect (CATE))plt.xlabel(Heterogeneity Feature (X))plt.legend(locupper left, prop{size: 10}) # 调整图例位置和大小plt.grid(True, linestyle--, alpha0.6)# 右侧子图CausalForestDML 结果plt.subplot(1, 2, 2) # 1行2列的第2个plt.title(CausalForestDML CATE Estimation)# 绘制 CausalForestDML 估计的 CATE 曲线plt.plot(X_test[:, 0], treatment_effects2, labelCausalForestDML CATE Estimate, colortab:green)# 绘制真实的 TE 曲线plt.plot(X_test[:, 0], expected_te, b--, labelTrue Effect (Ground Truth))# 填充置信区间区域plt.fill_between(X_test[:, 0], te_lower2, te_upper2, colortab:green, alpha0.3, label95% Confidence Interval)plt.ylabel(Treatment Effect (CATE))plt.xlabel(Heterogeneity Feature (X))plt.legend(locupper left, prop{size: 10})plt.grid(True, linestyle--, alpha0.6)plt.tight_layout() # 调整子图布局防止重叠plt.show()# 计算真实的平均处理效应 (Average Treatment Effect, ATE) 作为参考# 在我们的模拟中可以通过对真实 TE(x) 在 X 的分布上积分得到# 这里用测试点上的 TE 均值近似 ATEATE_for_cf np.mean(expected_te)print(f\nApproximate True Average Treatment Effect (ATE): {ATE_for_cf:.4f})# --- 创建并使用另一个 CausalForestDML 实例进行不同的推断 ---# 创建一个新的 CausalForestDML 实例使用不同的参数配置est_inf CausalForestDML(cv2, # 交叉拟合的折数criterionmse, # 森林分裂的标准 (均方误差)n_estimators400, # 更多的树min_var_fraction_leaf0.1, # 叶子节点样本占总样本最小比例 (控制复杂度)min_var_leaf_on_valTrue, # 在验证集上检查叶节点方差条件verbose0, # 不打印训练过程信息discrete_treatmentFalse, # !! 注意这里设为 False。对于二元T两种设置有时都可用但内部处理可能不同。# 设为 False 可能启用处理连续 T 的逻辑但对二元 T 仍能估计效应。n_jobs1, # CPU核心数random_state123 # 随机种子)print(\nStarting second CausalForestDML (for inference) training...)est_inf.fit(Y, T, XX, WW)print(Second CausalForestDML training finished.)# --- 边际效应推断 (Marginal Effect Inference) ---# 估计平均边际效应即假设效应对于所有X是恒定的print(Performing constant marginal effect inference...)res_me est_inf.const_marginal_effect_inference(X_test)point_me res_me.point_estimate # 平均效应的点估计lb_me, ub_me res_me.conf_int(alpha.01) # 99% 置信区间# 可视化平均边际效应 (常数线)plt.figure(figsize(10, 6))plt.plot(X_test[:, 0], point_me * np.ones_like(X_test[:, 0]), labelfConst. Marginal Effect Estimate ({point_me[0]:.2f}), colortab:red)plt.fill_between(X_test[:, 0], lb_me * np.ones_like(X_test[:, 0]), ub_me * np.ones_like(X_test[:, 0]), colortab:red, alpha0.3, label99% CI for Const. Marginal Effect)# 同时绘制真实的异质性效应作为对比plt.plot(X_test[:, 0], expected_te, b--, labelTrue Heterogeneous Effect)plt.xlabel(Heterogeneity Feature (X))plt.ylabel(Treatment Effect)plt.title(Constant Marginal Effect Inference vs True Effect)plt.legend(locupper left, prop{size: 10})plt.grid(True, linestyle--, alpha0.6)plt.show()# 打印边际效应推断的摘要信息print(\nConstant Marginal Effect Inference Summary:)print(res_me.summary_frame())# --- 处理效应推断 (Effect Inference for T0 vs T1) ---# 专门估计从 T0 变化到 T1 的处理效应 (即 CATE)print(\nPerforming effect inference (CATE)...)# T0 和 T1 定义了我们想比较的两个处理水平res_cate est_inf.effect_inference(X_test, T00, T11) # 对于二元TT00, T11 是标准做法point_cate res_cate.point_estimate # CATE(x) 的点估计lb_cate, ub_cate res_cate.conf_int(alpha0.01) # CATE(x) 的 99% 置信区间# 可视化 CATE 推断结果plt.figure(figsize(10, 6))plt.plot(X_test[:, 0], point_cate, labelCATE Estimate (Effect Inference), colortab:purple)plt.fill_between(X_test[:, 0], lb_cate, ub_cate, alpha.4, colortab:purple, label99% CI for CATE)plt.plot(X_test[:, 0], expected_te, b--, labelTrue Effect (Ground Truth))plt.xlabel(Heterogeneity Feature (X))plt.ylabel(Treatment Effect (CATE))plt.title(CATE Estimation and Inference (T0 vs T1))plt.legend(locupper left, prop{size: 10})plt.grid(True, linestyle--, alpha0.6)plt.show()第五阶段聚类分析目的:基于估计出的异质性处理效应 (CATE) 和/或其他特征将样本分成不同的组簇探索是否存在具有相似处理效应模式的子群体。讲解:这一阶段使用经典的K-Means 聚类算法来发现数据中的群体结构。寻找最佳聚类数量 (K):CH 指数衡量簇间离散度与簇内离散度的比率。分数越高通常表示聚类效果越好簇本身紧密簇之间分离度高。对于每个 K运行 K-Means 并计算 CH 分数。绘制 CH Score vs K 的曲线选择使得分最高的 K 值。注意:代码中计算 CH 分数时kmeans KMeans(n_clustersi, random_state1).fit(total_frame)使用了整个total_frame包含 CATE 置信区间、点估计、X 和 W作为聚类特征。这与计算 WCSS 时使用的特征不同需要注意。通常建议在同一组特征上应用这两种方法来选择 K。尝试不同的 K 值这里是从 2 到 24。对于每个 K运行 K-Means 并计算WCSS (Within-Cluster Sum of Squares)即每个簇内样本点到其质心距离的平方和。绘制 WCSS vs K 的曲线。理想情况下曲线会像一个手臂WCSS 随着 K 增大而减小在某个 K 值“肘部”之后下降速率显著减缓。这个“肘部”通常被认为是较好的 K 值。注意:代码中计算 WCSS 时kmeans_pca.fit(point.reshape(-1, 1))只使用了 CATE 点估计point_orig(来自训练数据的CATE) 作为聚类特征。K-Means 需要预先指定簇的数量 K。通常我们不知道最佳的 K 是多少。肘部法则 (Elbow Method):Calinski-Harabasz (CH) 指数:执行聚类:根据肘部法则或 CH 指数选择一个“最优”的 K。代码隐式地选择了使 CH 分数最高的 K并存储了对应的聚类标签hc_labels[np.argmax(hc_metrics)]。可视化聚类结果:代码绘制了一个散点图。X 轴是异质性特征Xtotal_frame[:, 3]因为前三列是效应估计。Y 轴是 CATE 的点估计total_frame[:, 2]。图中不同的颜色代表 K-Means 找到的不同簇。通过观察图我们可以看到具有不同X值和不同估计效应CATE的样本是如何被分到不同组的从而理解处理效应的异质性模式。聚类的意义:发现这些子群体有助于理解哪些人群对处理更敏感效应高或不敏感效应低。为不同群体制定差异化的干预策略个性化决策。进一步探索为什么这些群体的效应不同可能与X或W中的某些特定变量有关。# 导入聚类和降维相关库from sklearn.cluster import KMeans # K-Means 聚类算法from sklearn.decomposition import PCA # 主成分分析 (这里没用到但常与聚类结合)from sklearn import metrics # 包含多种聚类效果评估指标from sklearn.metrics import pairwise_distances # 用于计算距离 (这里没直接用到)import numpy as np# --- 评估不同聚类数量 K 的效果 ---wcss [] # 存储 WCSS (Within-Cluster Sum of Squares) 用于肘部法则hc_metrics [] # 存储 Calinski-Harabasz 分数hc_labels [] # 存储每个 K 值对应的聚类标签# 尝试 K 从 2 到 24print(\nEvaluating different numbers of clusters (K)...)k_range range(2, 25)for i in k_range:# 1. 使用 CATE 点估计进行 K-Means计算 WCSS (用于肘部法则)# 注意这里只用了 CATE 点估计 point_orig (来自训练集) 作为特征kmeans_pca KMeans(n_clustersi, initk-means, random_state42, n_init10) # n_init10 多次初始化选最优kmeans_pca.fit(point_orig.reshape(-1, 1)) # 需要 reshape 成 N x 1 的二维数组wcss.append(kmeans_pca.inertia_) # inertia_ 属性就是 WCSS# 2. 使用整个 total_frame 进行 K-Means计算 CH 分数和存储标签# 注意这里使用了所有特征 (CATE CI, CATE point, X, W)kmeans KMeans(n_clustersi, random_state1, n_init10) # 随机种子不同于上面kmeans.fit(total_frame)labels kmeans.labels_ # 获取每个样本所属的簇标签hc_labels.append(labels) # 存储这组标签# 计算 CH 分数需要原始数据和标签# 如果样本量很大这个计算可能较慢try:score metrics.calinski_harabasz_score(total_frame, labels)hc_metrics.append(score)except ValueError:# 如果只有一个簇或者数据有问题可能会报错hc_metrics.append(np.nan) # 添加 NaN 值print(fWarning: Could not compute CH score for K{i})score hc_metrics[-1]ch_score_str f{score:.2f} if not np.isnan(score) else NaNprint(f K{i}: WCSS{wcss[-1]:.2f}, CH Score{ch_score_str})# --- 绘制评估曲线 ---# 绘制肘部法则图 (WCSS vs K)plt.figure(figsize(10, 6))plt.plot(k_range, wcss, markero, linestyle-., colorred)plt.xlabel(Number of Clusters (K))plt.ylabel(WCSS (Within-Cluster Sum of Squares))plt.title(Elbow Method for Optimal K (based on CATE point only))plt.xticks(k_range) # 显示所有K值刻度plt.grid(True, linestyle--, alpha0.6)plt.show()# 绘制 Calinski-Harabasz 分数图 (CH Score vs K)# 过滤掉 NaN 值以防绘图错误valid_k [k for k, score in zip(k_range, hc_metrics) if not np.isnan(score)]valid_scores [score for score in hc_metrics if not np.isnan(score)]if valid_scores: # 确保有合法的分数可供绘制plt.figure(figsize(10, 6))plt.plot(valid_k, valid_scores, markero, linestyle-., colorblue) # 使用不同的颜色plt.xlabel(Number of Clusters (K))plt.ylabel(Calinski-Harabasz Score)plt.title(Calinski-Harabasz Score for Optimal K (based on all features))# 标记最高分对应的 Kbest_k_ch valid_k[np.argmax(valid_scores)]max_score np.max(valid_scores)plt.scatter(best_k_ch, max_score, colorred, s100, zorder5, labelfBest K {best_k_ch} (Score{max_score:.2f}))plt.xticks(valid_k)plt.legend()plt.grid(True, linestyle--, alpha0.6)plt.show()else:print(No valid Calinski-Harabasz scores were computed.)best_k_ch -1 # 标记为无效# --- 可视化最优聚类结果 (基于 CH 分数最高的 K) ---if best_k_ch ! -1:print(f\nVisualizing clusters for the best K {best_k_ch} (based on CH score)...)# 获取最佳 K 对应的标签best_labels hc_labels[valid_k.index(best_k_ch)] # 注意索引对应关系plt.figure(figsize(10, 8))# 使用 seaborn 的调色板获取不同颜色palette sns.color_palette(husl, best_k_ch)# 遍历每个簇 (0 到 best_k_ch-1)for i in np.unique(best_labels):# 找到属于当前簇 i 的所有样本的索引spots np.where(best_labels i)# 绘制这些样本点# X轴: total_frame[:, 3] (即第一个协变量 X[:, 0])# Y轴: total_frame[:, 2] (即 CATE 的点估计 point_orig)plt.scatter(total_frame[spots, 3], total_frame[spots, 2],labelfCluster {i}, colorpalette[i], alpha0.7, s30) # s控制点大小plt.xlabel(Heterogeneity Feature (X))plt.ylabel(Estimated CATE)plt.title(fK-Means Clustering Result (K {best_k_ch}))plt.legend(titleClusters)plt.grid(True, linestyle--, alpha0.6)plt.show()else:print(Cannot visualize clusters as no best K was determined.)应用到其他数据集的关键步骤要将本教程介绍的因果推断流程和模型应用到您自己的数据集您需要遵循以下关键步骤明确因果问题与数据准备:结果变量 Y (Outcome):你最终关心的那个结果指标是什么例如销售额、用户留存率、康复时间处理变量 T (Treatment):你想评估其效果的干预、行动或特征是什么它是二元的如用药/未用药还是连续的如广告投入金额协变量 X (Covariates / Features for Heterogeneity):你认为哪些个体特征可能会调节处理变量的效果例如用户年龄、用户地理位置、产品类别控制变量 W (Controls / Confounders):哪些变量同时影响了个体接受处理的可能性 T和结果 Y必须尽可能全面地识别并包含这些变量以控制混淆。例如用户的历史购买行为、病情严重程度、学生入学前的基础工具变量 Z (Instrumental Variable):仅当使用 IV 方法如 DeepIV 时需要是否存在一个变量 Z它能影响个体接受处理 T 的决策但它本身不会直接影响结果 Y除非通过 T 产生影响并且 Z 与 Y 的未观测影响因素无关寻找有效的 Z 非常困难需要很强的领域知识和严格的检验。定义变量:清晰地界定你的目标数据收集与清洗:收集包含上述所有变量的数据。进行必要的数据清洗、处理缺失值、转换数据类型等。数据合并:确保所有变量都在一个统一的数据结构中如 Pandas DataFrame。选择合适的因果模型:是否存在内生性/未观测混淆如果你强烈怀疑有重要的未观测混淆因素并且能找到有效的工具变量 Z那么DeepIV是一个可以考虑的强大选项。是否主要关心异质性处理效应 (HTE)如果你想了解处理效果如何随个体特征 X 变化并且假设你已经包含了足够多的控制变量 W 来处理主要的混淆那么基于 DML 的森林方法如DROrthoForest或CausalForestDML是非常好的选择。它们对 W 的建模有较好的鲁棒性。处理变量 T 是离散还是连续的确保所选模型的参数设置与你的 T 类型匹配例如CausalForestDML中的discrete_treatment参数。数据量大小:深度学习模型如 DeepIV通常需要较大数据量。森林模型在中小规模数据上表现也很好。模型配置与训练:正则化:对于 Lasso/Ridge 等线性模型正则化强度如alpha或C很重要可能需要通过交叉验证选择或者使用类似代码中基于理论的启发式计算。森林参数:树的数量 (n_estimators/n_trees)、树深 (max_depth)、叶节点大小 (min_samples_leaf/min_leaf_size)、子采样比例 (max_samples/subsample_ratio) 等都需要根据数据情况调整以平衡偏差和方差。神经网络参数:对于 DeepIV网络结构层数、节点数、激活函数、优化器、学习率、正则化、Dropout 比例等都需要仔细设计和调整。选择内部估计器:对于 DML 方法你需要为结果模型 (model_y) 和处理/倾向得分模型 (model_t) 选择合适的机器学习模型如 Lasso, Ridge, GradientBoostingRegressor/Classifier, 甚至其他森林。模型的选择会影响最终效果。调整超参数:执行训练:使用fit(Y, T, XX, WW, ZZ)根据模型需要传入相应变量来训练模型。估计与解释因果效应:估计 CATE:使用effect(X_new)方法估计在新数据点X_new或者原始数据X上的条件平均处理效应。估计 ATE:有些模型提供直接估计平均处理效应ATE的方法如est.ate()或在summary()中查看或者可以通过对 CATE(X) 在 X 的分布上求平均来近似得到 ATE。计算置信区间:使用effect_interval()或effect_inference()等方法获取效应估计的置信区间这对于理解估计的不确定性至关重要。一个很宽的置信区间意味着估计结果不太可靠。可视化:绘制 CATE(x) 随某个重要特征 x 变化的曲线图类似教程中的做法并带上置信区间。这有助于直观理解异质性。解读结果:结合业务或领域知识来解释估计出的因果效应的大小、方向和异质性模式是否有意义。效应是正向还是负向统计上是否显著置信区间是否包含0哪些人群受益最多/最少敏感性分析与模型检验 (重要但本项目未深入):检验假设:对于 IV 方法需要检验工具变量的有效性。对于 DML 方法虽然对模型设定有鲁棒性但仍需检查内部模型拟合情况。安慰剂检验 (Placebo Tests):例如将一个已知没有效应的变量作为“处理”变量看模型是否错误地估计出显著效应。对未观测混淆的敏感性分析:尝试量化结果对可能存在的、未包含在 W 中的混淆因素的敏感程度。EconML 也提供了一些相关工具。发现子群体 (可选):如果对发现具有相似效应模式的子群体感兴趣可以像教程中那样使用估计出的 CATE 或其他相关特征进行聚类分析如 K-Means。分析每个簇的特征X 和 W 的分布理解造成效应差异的原因。通过遵循这些步骤并结合对 EconML 库文档的进一步学习应该能够将这些强大的因果推断技术应用于您自己的实际问题中。记住因果推断往往比单纯的预测更具挑战性需要更仔细的思考、假设检验和结果解读。