因果中介分析：双机器学习与非参数估计框架解析

1. 因果中介分析：从理论困境到非参数估计的破局之路

在社会科学、流行病学、经济学乃至机器学习公平性研究中，我们常常需要回答一个核心问题：一个干预（比如一项新政策、一种药物、一个算法调整）是如何起作用的？它的效果是直接作用于最终结果，还是通过改变某个中间变量（比如收入水平、生物标志物、模型特征）来间接实现的？这就是因果中介分析要解决的难题。简单来说，它试图量化“通过中介路径的效应”和“不通过中介路径的效应”各自有多大。

过去几十年，以自然直接效应和自然间接效应为核心的理论框架占据了主导地位。它们的概念直观：想象给同一个人同时施加治疗和对照两种状态，一种状态下让其中介变量“冻结”在未接受治疗时的自然水平，另一种状态下则让其“自然发展”，两者的差异就分别代表了直接效应和间接效应。然而，这个看似完美的构想建立在“跨世界反事实”这一哲学上无法验证、实践中极易被违反的假设之上。一旦存在既受治疗影响、又同时影响中介和结果的变量——即“中间混杂因素”，整个分析的基础就崩塌了。更棘手的是，在基因组学、微生物组研究或高维特征分析中，中介变量本身可能就是成百上千个连续或高维的指标，传统方法几乎束手无策。

因此，领域内催生了一系列试图绕过“跨世界假设”的新参数，如随机干预效应、可分离效应等。但理论上的突破并未完全解决实践中的估计难题，尤其是面对复杂、高维的数据时。本文将深入探讨如何运用前沿的双机器学习与非参数估计技术，构建一套能够稳健处理连续/高维中介变量和中间混杂因素的通用估计框架。这不仅是一次方法学的梳理，更是面向实际数据分析挑战的一次“工具箱”升级。

2. 核心概念与参数演化：超越自然效应

在深入技术细节前，我们必须厘清几个核心的因果中介参数，理解它们各自的优劣与适用场景。这就像医生开药前，得先搞清楚各种药物的机理和副作用。

2.1 自然直接效应与间接效应：理想与现实的鸿沟

自然直接效应与间接效应是因果中介分析的基石。其定义依赖于一个思想实验：比较个体在两种反事实情形下的结果。

• 自然直接效应：将治疗A设置为1（如用药），但同时将中介变量M“固定”在该个体如果未接受治疗（A=0）时会自然呈现的水平M(0)。这衡量了治疗“绕过”中介变量直接对结果Y产生的影响。
• 自然间接效应：将治疗A设置为1，但比较中介变量取其在治疗下的自然水平M(1)与在对照下的自然水平M(0)时，结果的差异。这衡量了治疗通过改变中介变量而对结果产生的效应。

注意：这里的关键在于“自然水平”M(0)和M(1)是同一个体在不同治疗下的潜在结果。NDE的计算需要同时考虑A=1和M=M(0)，这构成了一个“跨世界”的陈述，因为现实中我们无法同时观测到一个个体既接受治疗又拥有其未接受治疗时的中介变量值。这引出了识别NDE/NIE所需的“跨世界独立性”假设，该假设在存在中间混杂因素Z时必然不成立。

2.2 主流替代方案：为现实世界设计的参数

为了克服NDE/NIE的缺陷，研究者提出了多种不依赖跨世界假设的替代参数。它们可以看作是在不同约束条件下，对“直接”和“间接”机制的不同近似和量化方式。

路径特异性锐零准则是评估一个中介参数是否真正捕捉到“机制”的重要标准。它要求：如果没有任何一个个体经历通过某条特定路径的效应，那么该路径的效应估计值应该恰好为零。RIDE/RIIE可能违反这一准则，而可分离效应和反证双胞胎效应则满足它。

2.3 统一识别框架：两大核心统计泛函

尽管上述参数定义各异，但一个令人振奋的发现是，对于二元处理（A=0或1）的情况，这六大类参数的识别公式都可以归结为两个核心的统计泛函：ψ^N 和 ψ^R。这为构建统一的估计器奠定了理论基础。

这两个泛函本质上是经过特定条件期望嵌套处理的平均结果函数：

• ψ^N(a1, a2, a3)：其计算涉及对中介变量Z和M的分布进行积分，但积分是在“自然”条件下进行的，即Z和M的分布与某个特定的处理水平a2, a3相关联。这对应了自然效应、决策理论效应和有机效应。
• ψ^R(a1, a2, a3, a4)：其关键创新在于引入了一个置换变量Z^π。这个变量满足：在给定(A,W)的条件下，Z^π与Z同分布，但Z^π与M条件独立。这相当于“打破”了Z到M的因果箭头，从而在计算积分时，可以将对Z分布的复杂积分，转化为一个关于Z^π的条件期望问题，极大简化了高维情况下的估计。这对应了随机干预效应和基于置换思想的路径效应（如可分离效应、反证双胞胎效应的部分路径）。

通过为(a1,a2,a3,a4)选择不同的取值组合，我们可以从ψ^N和ψ^R中恢复出表1中的所有效应。例如：

• NDE = ψ^N(1,0,0) - ψ^N(0,0,0)
• RIIE = ψ^R(1,1,1,1) - ψ^R(1,1,0,0)
• 可分离效应中的路径P2效应：SE2 = ψ^R(0,1,1,1) - ψ^R(0,0,1,1)

这种统一性意味着，我们只需要集中精力为ψ^N和ψ^R这两个“母函数”开发出强大、稳健的估计方法，就能一劳永逸地估计所有类型的中介效应。

3. 非参数估计的挑战与双机器学习框架

面对高维中介变量Z和M，传统的参数化模型（如线性回归、逻辑回归）会因模型误设而产生严重偏差。而非参数方法虽然灵活，但直接估计ψ^N和ψ^R中涉及的多重嵌套积分和条件密度，在计算上是灾难性的。这就是双机器学习框架大显身手的地方。

3.1 插值估计器的问题与一阶偏误修正

最直观的估计方法是“插值法”：用机器学习模型（如梯度提升树、神经网络）依次拟合数据中的条件期望函数（如E[Y|A,Z,M,W]），然后将估计出的函数代入ψ^N或ψ^R的公式中计算。然而，即便使用了最先进的机器学习算法，这种插值估计器通常也存在一阶偏误。这是因为机器学习模型追求的是总体预测精度，而非我们关心的特定因果参数估计的无偏性。其偏差的阶数可能达不到O_p(n^{-1/2})，导致基于中心极限定理的置信区间覆盖概率失效。

双机器学习的核心思想，就是刻画并校正这个一阶偏误。它依赖于半参数效率理论中的一个关键工具：有效影响函数。

3.2 有效影响函数：估计的“罗盘”

有效影响函数是参数ψ在分布P处的某种“导数”，它描述了当数据分布发生微小扰动时，参数估计值应该如何变化。对于我们的泛函ψ^N和ψ^R，其有效影响函数φ^N(X; η^N)和φ^R(X; η^R)具有类似以下的结构（以φ^N为例）：

φ^N = α_1 * (伪结局 - θ_1) + α_2 * (伪结局 - θ_2) + α_3 * (Y - θ_3) + (θ_1的预测值) - ψ^N

其中：

• θ_1, θ_2, θ_3是我们需要估计的回归函数（即一系列条件期望）。
• α_1, α_2, α_3被称为Riesz表示子，它们本质上是处理机制（倾向得分）和中介变量条件密度比的复杂组合，其作用是给每一步回归的残差进行加权，以纠偏。

有效影响函数有两个至关重要的性质：

1. 无偏性：无论我们用什么样的估计值η_hat（包含对θ和α的估计）去计算φ(X; η_hat)，只要对于每一步k，我们要么正确估计了回归函数θ_k，要么正确估计了Riesz表示子α_k，那么φ的期望就等于真值ψ。这就是双重稳健性的根源。
2. 方差界：基于有效影响函数构造的估计器，其渐近方差可以达到非参数模型下的效率下界，即理论上最“准”的估计所能达到的精度。

3.3 Riesz学习：绕过密度估计的巧思

估计Riesz表示子α是最大的技术难点。以α^N_2为例，其理论形式包含棘手的密度比P(z|a3,w)/P(z|a,w)。当Z是高维连续变量时，直接估计条件密度P(z|a,w)几乎是不可能的任务。

Riesz学习提供了一种巧妙的回避方案。它不直接估计密度，而是将α视为一个函数逼近问题。根据Riesz表示定理，α是某个线性泛函的唯一表示子。我们可以通过最小化一个特定的均方误差损失函数来直接学习α。

例如，对于α^N_1，其损失函数为：L(α) = E[α(A,W)^2 - 2 * b^N_1(W) * α(A,W)]其中b^N_1是一个已知的、由前一步回归函数构造的伪结局。我们可以使用任何灵活的回归器（如神经网络）来最小化这个损失函数的样本版本，从而得到α^N_1的估计。α^N_2, α^N_3以及ψ^R对应的α^R_k都可以通过类似定义的序列损失函数进行估计。

实操心得：在实践中，我们通常使用深度学习框架（如PyTorch, TensorFlow）来实现Riesz学习，因为它们允许自定义损失函数。对于α的模型，一个包含若干隐藏层的全连接神经网络通常是个不错的起点。需要确保对α的输出不做任何范围限制（如sigmoid），因为Riesz表示子理论上可以取任意实数值。

4. 完整估计流程与实现细节

下面，我们以估计ψ^N为例，拆解整个双机器学习估计器的构建步骤。估计ψ^R的流程类似，但步骤更多一层。

4.1 步骤一：数据准备与样本分割

1. 数据：观测到 i.i.d. 样本(W_i, A_i, Z_i, M_i, Y_i), i=1,...,n。其中W是基线协变量，A是二元处理，Z是中间混杂因素，M是中介变量，Y是结局。
2. 交叉拟合：为了理论上的严谨性（避免Donsker类条件），我们采用K折交叉拟合。将数据随机分为K份（通常K=5或10）。对于每一折k，用其余所有折的数据作为训练集T_k，来估计模型参数，然后用这些参数对第k折的数据进行预测。这保证了用于构造估计量的“模型拟合”和“影响函数计算”两部分数据是独立的。

4.2 步骤二：序列回归估计θ

对于每一折的训练集T_k，我们进行以下三步回归（对应ψ^N中的三层期望）：

1. 估计 θ^N_3：用T_k中的数据，以(A, Z, M, W)为特征，Y为结局，训练一个机器学习模型g_3（如梯度提升树、随机森林、神经网络）。这个模型用于估计E[Y | A, Z, M, W]。
2. 构造伪结局并估计 θ^N_2：
   • 对于T_k中的每个观测i，计算b^N_3i = g_3(a1, Z_i, M_i, W_i)。这里a1是我们在ψ^N中设定的目标处理值（例如，对于NDE中的ψ^N(1,0,0)，a1=1）。
   • 然后，以(A, Z, W)为特征，以刚计算出的b^N_3为新的结局变量，在T_k上训练模型g_2，用于估计E[b^N_3 | A, Z, W]。
3. 构造伪结局并估计 θ^N_1：
   • 对于T_k中的每个观测i，计算b^N_2i = g_2(a2, Z_i, W_i)（例如，对于NDE，a2=0）。
   • 以(A, W)为特征，以b^N_2为结局，训练模型g_1，用于估计E[b^N_2 | A, W]。

至此，我们得到了回归函数序列θ^N_3, θ^N_2, θ^N_1的估计。

4.3 步骤三：Riesz学习估计α

同样在每一折的训练集T_k上，我们序列化地学习Riesz表示子：

1. 估计 α^N_1：
   • 计算目标值：对于T_k中每个i，计算b^N_1i = g_1(a3, W_i)（对于NDE，a3=0）。
   • 定义损失函数：L_1(α) = (1/|T_k|) Σ_{i in T_k} [α(A_i, W_i)^2 - 2 * b^N_1i * α(A_i, W_i)]。
   • 用一个模型h_1（如神经网络）来最小化L_1，得到的函数就是α^N_1的估计。
2. 估计 α^N_2：
   • 利用上一步估计的α^N_1，计算权重：weight_i = α^N_1(A_i, W_i)。
   • 定义损失函数：L_2(α) = (1/|T_k|) Σ_{i in T_k} [α(A_i, Z_i, W_i)^2 - 2 * weight_i * b^N_2i * α(A_i, Z_i, W_i)]。
   • 训练模型h_2最小化L_2，得到α^N_2的估计。
3. 估计 α^N_3：
   • 利用α^N_2的估计值作为权重。
   • 定义损失函数：L_3(α) = (1/|T_k|) Σ_{i in T_k} [α(A_i, Z_i, M_i, W_i)^2 - 2 * α^N_2(A_i, Z_i, W_i) * b^N_3i * α(A_i, Z_i, M_i, W_i)]。
   • 训练模型h_3最小化L_3，得到α^N_3的估计。

4.4 步骤四：构造一步估计量及其方差

对于测试折V_k中的每一个观测i：

1. 使用在训练集T_k上拟合的所有模型(g_1, g_2, g_3, h_1, h_2, h_3)，计算该观测的φ^N(X_i; η_hat)值。这需要代入该观测的特征值，并利用上述模型计算相应的预测值和权重。
2. 计算b^N_1(W_i)的预测值，即g_1(a3, W_i)。

那么，对于该折数据，ψ^N的估计值为该折所有观测的φ^N值的平均。最终，将K折的估计结果平均，就得到了基于全部数据的一步估计量ψ_hat^N。

方差估计与置信区间： 根据定理2，ψ_hat^N是渐近正态的，其渐近方差为Var[φ^N(X; η)]/n。因此，我们可以用样本方差来估计：σ_hat^2 = (1/n) Σ_{i=1}^n [φ^N(X_i; η_hat_{j(i)}) - ψ_hat^N]^2其中η_hat_{j(i)}是用于预测第i个观测的、在对应训练折上拟合的模型集合。那么，ψ_hat^N的标准误为SE = sqrt(σ_hat^2 / n)，95%置信区间为ψ_hat^N ± 1.96 * SE。

4.5 处理连续或多变量治疗：修正治疗策略

上述框架针对二元处理A。对于连续或多变量的治疗（如药物剂量、多维度政策包），我们需要新的因果参数定义。一个强大的工具是修正治疗策略。

其核心思想是，我们不比较“治疗 vs. 不治疗”，而是比较两种不同的治疗分配规则。例如，规则d1(a,w)可能是在当前剂量a基础上增加一个固定量，而规则d2(a,w)可能是根据个体特征w设定一个个性化剂量。我们关心的效应是：如果整个人群都遵循规则d1，与都遵循规则d2相比，结果Y的期望差异。

通过将MTD与可分离效应的思想结合，我们可以为连续治疗定义路径特异性效应。这需要假设治疗A可以概念上分解为只影响中介的A_M和只影响结果的A_Y。然后，我们定义两种MTD：一种只改变A_M成分，另一种只改变A_Y成分。比较这两种策略下的结果，就能分离出通过中介的效应和直接效应。

估计时，前述的双机器学习框架依然适用，但需要将处理变量A替换为根据MTD规则生成的“处理值”，并相应地调整倾向得分模型（即P(A|W)）为处理机制模型P(A|W)在MTD下的对应形式。Riesz学习部分也需要进行相应调整，但核心逻辑不变。

5. 实操要点、常见陷阱与解决方案

在实际应用这套方法时，以下几个环节最容易出错，需要格外留意。

5.1 模型选择与过拟合控制

• 回归模型（θ）的选择：对于θ_3（预测Y），如果Y是连续变量，可以使用回归树、梯度提升、神经网络或弹性网络。对于分类结局，使用相应的分类算法。关键在于使用交叉验证选择超参数，避免过拟合。
• Riesz表示子（α）模型的选择：神经网络因其灵活性成为首选。需要注意：
  • 网络结构：不宜过深，2-4个隐藏层通常足够。宽度应适中，确保有足够容量但不过度参数化。
  • 激活函数：输出层使用线性激活，因为α可取任意实数值。隐藏层使用ReLU、Tanh等。
  • 正则化：强烈推荐使用权重衰减或Dropout。因为损失函数α^2 - 2bα中，如果不对α的大小加以约束，在b很大时，模型可能倾向于学习一个数值上巨大且不稳定的α，导致估计方差爆炸。
  • 优化器：Adam或AdamW优化器通常表现良好，学习率需要仔细调优。

踩坑实录：在早期尝试中，我曾未对α网络施加任何正则化，结果在模拟数据中出现了极端的α估计值（某些点上的值超过1e6），导致最终效应估计的置信区间宽到毫无意义。加入L2正则化后，估计立刻稳定下来。

5.2 中间混杂因素Z的置换构造

对于ψ^R的估计，我们需要生成置换变量Z^π。其要求是：在给定(A, W)的条件下，Z^π的分布与Z相同，但Z^π与M条件独立。一个简单可靠的实践方法是分层置换：

1. 将数据按照处理水平A和（离散化的）协变量W的 strata 进行分组。
2. 在每个 strata 内部，对Z的值进行随机重排（permutation）。
3. 将重排后的Z值作为该观测的Z^π。

这样，在每个 strata 内，Z^π的边际分布与Z相同，但由于是随机重排，Z^π与原始的M在给定(A,W)下是独立的。这完美满足了Z^π的构造要求。

5.3 双重稳健性的实际含义与诊断

双重稳健性是我们的“安全网”。它意味着，对于每一步k（例如，估计E[b_3 | A,Z,W]的这一步），我们允许θ_k的模型或α_k的模型有误设，只要其中一个正确，最终估计就是一致的。

如何利用这一性质？

1. 模型堆叠：对于关键的θ模型（尤其是第一步的θ_3），可以使用超级学习器来组合多个基学习器，以提升正确设定的概率。
2. 敏感性分析：尝试使用完全不同类型的模型来估计θ序列和α序列（例如，用树模型估计θ，用神经网络估计α）。如果两种组合得到的点估计和置信区间相近，则结果更可信。
3. 部分验证：虽然无法直接验证α_k是否正确，但可以验证θ_k的预测性能（通过交叉验证的R²或AUC）。如果θ_k的预测性能很差，那么我们对α_k正确设定的依赖就更大，此时需要更仔细地调优α模型或尝试不同的函数形式。

5.4 高维中介变量的特殊处理

当M的维度极高（如基因表达数据）时，直接将其作为回归特征会导致维数灾难。此时有几种策略：

1. 降维：先对M进行主成分分析、因子分析或使用自编码器提取低维表征，然后用这个低维表征替代原始的M进入模型。
2. 结构化模型：如果M的不同维度之间有已知的结构（如网络结构），可以在神经网络中引入相应的结构（如图卷积层）。
3. 简化参数：考虑使用更聚合的中介效应参数，例如只关心通过M的某个标量汇总统计量（如M的第一主成分）的效应，而不是通过每一个维度的效应。

6. 模拟与实例：从理论到数据

为了验证方法的有效性，我们通常需要进行模拟研究。模拟数据生成过程应包含复杂的非线性关系、交互效应以及中间混杂因素Z。例如，可以这样生成数据：

W ~ Normal(0, I) A ~ Bernoulli(expit(δ * W)) # 处理受W影响 Z = f_z(W, A) + ε_z # 中间混杂受W和A影响 M = f_m(W, A, Z) + ε_m # 中介受W, A, Z影响 Y = θ_y + β_a * A + β_z * Z + β_m * M + g(W, A, Z, M) + ε_y # 结果，包含直接效应和通过Z、M的间接效应

其中f_z,f_m,g可以是复杂的非线性函数。在这样的数据上，比较我们提出的双机器学习估计器与传统的参数化G-computation、基于错误模型的估计器的表现，评估其偏差、均方误差和置信区间覆盖率。

在一个真实的慢性疼痛对阿片类药物使用障碍影响的研究中，W可能包括人口学特征和共病，A是慢性疼痛诊断，Z可能是疼痛导致的社会功能下降或抑郁症状（它既受疼痛影响，又影响后续的止痛药使用和成瘾风险），M是阿片类药物的处方模式或使用剂量，Y是阿片类药物使用障碍的诊断。应用本文方法，可以量化慢性疼痛通过改变处方模式（M）的间接效应，与直接导致成瘾（或通过抑郁症状Z）的直接效应，即使在处方模式M是多维、连续的情况下。

7. 总结与展望

从依赖脆弱跨世界假设的自然效应，到基于随机干预、可分离成分等更务实定义的效应，因果中介分析的理论框架日益丰富和稳健。然而，只有配以同样强大的估计工具，这些理论才能在实践中焕发生机。本文阐述的基于双机器学习和Riesz学习的非参数估计框架，正是这样一座连接稳健理论与复杂数据的桥梁。

这套方法的核心优势在于其统一性与灵活性。一套估计流程，通过设定不同的参数(a1,a2,a3,a4)，就能适配六大类中介效应参数。通过Riesz学习规避了高维密度估计的诅咒，通过交叉拟合和双重稳健性保障了统计推断的可靠性。面对当今研究中日益常见的连续、多维度中介变量和复杂的混杂结构，它提供了一条可行的分析路径。

当然，工具越强大，对使用者的要求也越高。它要求分析师不仅理解因果图模型和识别假设，还要熟悉机器学习模型的训练与调优。计算成本也显著高于传统方法。未来的发展方向可能包括开发更高效稳定的Riesz表示子学习算法、将其扩展到纵向中介分析、以及开发用户友好的软件包以降低应用门槛。

在我自己的应用经验中，最关键的一步永远是因果图的绘制与审视。清晰界定W、A、Z、M、Y，认真论证是否存在未测量的混杂，是任何高级估计方法都无法替代的基础。在这个基础上，本文的方法才能发挥其最大价值，帮助我们从观测数据中，更清晰、更可靠地剥离出那些交织在一起的因果路径。