基于因果图的协变量数据缺失因果效应估计算法

Causal Effect Estimation Algorithm Based on Causal Graph with Missing Covariate Data

摘要

针对随机对照实验数据难获取的情况，基于观测数据的因果推断成为替代性研究范式。而观测数据中协变量的缺失往往引入偏差，导致因果效应估计失真。为此，提出一种在协变量缺失场景下三阶段加权因果效应估计算法。首先，基于因果图将协变量划分为混杂变量、暴露预测因子及结果预测因子三类，并针对混杂变量缺失项实施多重插补以保持因果结构完整性；其次，引入缺失模式，通过改进的协变量平衡策略构建倾向得分模型；最后，整合各缺失模式并采用逆概率加权估计因果效应。仿真实验结果表明，相较于梯度提升机，该算法在绝大多数场景下的估计均方根误差更低，能在数据缺失场景下有效估计因果效应。

Abstract

Given the challenges of obtaining randomized controlled trial data, causal inference based on observational data has emerged as an alternative research paradigm. However, missing covariates in observational data often introduce bias, leading to distorted estimates of causal effects. To address this problem, this paper proposes a three-stage weighted causal effect estimation algorithm for handling missing covariates. First, covariates are categorized into confounders, exposure predictors, and outcome predictors based on the causal graph, and a multiple imputation strategy is implemented for missing confounders to preserve the integrity of the causal structure. Second, a missing pattern is introduced, and propensity score models are constructed through an improved covariate balancing strategy. Finally, causal effects are estimated by integrating inverse probability weighting across missing patterns. Simulation results demonstrate that, compared with existing methods such as gradient boosting machines, the proposed method reduces the root mean square error (RMSE) in most scenarios, validating its effectiveness in estimating causal effects under missing data conditions.

Graphical abstract

关键词

因果推断 / 协变量缺失数据 / 处理效应 / 倾向得分

Key words

causal inference / missing data for covariates / treatment effect / propensity score

0　 引 言

使用随机对照实验数据进行因果效应估计被认为是因果推断的金标准^[1]。但是在社会经济领域研究中，由于现实条件的限制，通常不能开展随机对照实验^[2]，因此学者们开始研究如何使用观测性研究数据(以下简称“观测数据”)代替随机对照实验数据进行因果效应估计。这样做通常存在两个问题：一是在估计处理变量(Treatment)对结果变量(Outcome)的平均处理效应(Average Treatment Effect，ATE)时会受到混杂因素(Confounder)的影响，导致更大的偏差；二是观测数据中往往存在协变量(Covariate)数据缺失的情况。

为了控制因果效应估计中混杂因素带来的偏差，学者们基于混杂变量的可观测性提出了许多解决方法。当所有混杂变量都能被观测到时，常用方法有基于倾向得分的估计方法^[3-4]、基于回归的估计方法^[5]和将二者结合的双鲁棒方法^[6]。前两种方法对函数形式敏感，模型误设会导致偏差；双鲁棒方法要求至少一个模型的设定正确，估计效果对模型的设定仍有一定程度的依赖。为进一步克服估计对模型设定的依赖，学者们又提出了非参数方法^[7]和交错双重差分法^[8]。非参数方法计算复杂度高且受异常值的影响较大，交错双重差分法要求数据通过平行性检验。当部分混杂变量不能被观测到时，估计困难增大，常用方法有工具变量的估计方法^[9]、回归不连续的估计方法^[10]、阴性对照变量方法^[11]和基于数据融合的方法^[12-14]。工具变量法要求严格的外生性假设，现实工具变量往往难以满足；回归不连续设计仅适用于局部效应，外推性差；阴性对照变量方法需未观测混杂满足特定结构，适用性有限；数据融合方法依赖于多源数据集的互补以提升估计精度，而数据集的可获取性已成为制约该方法的瓶颈。综上，现有方法均依赖于不同程度的假设条件，其适用性受限于数据质量和实际问题的复杂性。尤其是模型存在不确定性或数据缺失的情况下，因果推断的稳健性受到影响。

当协变量数据存在缺失时，混杂因素的有效控制面临挑战，常用解决方法主要有三类：完整个体分析法、缺失指标法和插补法。完整个体分析法删除数据中存在缺失值的个体，仅使用协变量数据完整的个体进行估计，当缺失类型为完全随机缺失时，使用完整个体分析法能得到无偏的估计值。由于删除了部分观测个体，当不完整的个体数量较多时，该方法不适用^[15]。缺失指标法引入缺失指标，当某个个体的协变量存在缺失值时，其缺失指标记为1；无缺失值时，其缺失指标记为0^[16]。缺失指标法会在回归模型中引入偏差，导致因果效应估计偏差增大^[17]。插补方法相较于前两种方法应用更加普遍，通常分为单一插补(Single Imputation，SI)和多重插补(Multiple Imputations，MI)。单一插补方法相对简单，易于操作，适合缺失比例低的数据。对于缺失比例高的数据，多重插补方法更加适用，该方法由Rubin^[18]提出，给每个缺失值都构造多个替代值，产生多个完整数据集，综合所有数据集的处理结果，最终得到目标变量的估计。多重插补方法相较于单一插补方法稳健性更好，偏差更小^[19]。但插补法可能会导致数据分布的失真^[20]，因此有学者开始考虑无需删除或填补的缺失值处理方法。McCaffrey等^[21]采用梯度提升机(Gradient Boosting Machine，GBM)方法在数据缺失的情况下对倾向得分进行估计。相较于逻辑回归等传统方法，GBM方法由于需迭代构建多个弱学习器并优化损失函数，其算法复杂度显著提高，时间复杂度达O(T\cdot n\cdot d)，逻辑回归的时间复杂度为O(n\cdot d)，其中T为迭代次数，n为样本量，d为特征维度，在相同硬件条件下，计算耗时明显增加。Wang等^[22]在混合离散和非可忽略缺失协变量的情况下，基于贝叶斯方法进行分位数回归，但需要指定缺失协变量的分布。综上，现有处理协变量缺失的方法中，删除个体的方法会造成原始数据信息丢失；对全部协变量进行填补则会造成原始数据分布失真；无需删除和填补的GBM方法时间成本更高；贝叶斯方法则需要依赖于一定的假设条件。

针对协变量数据完全随机缺失机制下的因果效应估计问题，本文提出三阶段加权的因果效应估计方法(Three-Stage Weighted causal inference，TSW)。该方法旨在减少插补数据带来的偏差，提高因果效应估计的准确性。TSW方法首先基于因果图的因果关系把协变量分为3组：混杂(Confounder)变量组、暴露预测(Exposure Predictor)因子组和结果预测(Outcome Predictor)因子组；接着对混杂变量组的缺失值使用多重插补方法进行插补；然后基于协变量分组构建一系列的子模型，通过加权得到倾向得分；最后使用加权倾向得分进行ATE估计。

1　 因果效应与TSW方法

1.1　基于潜在结果框架的因果效应定义

Rubin ^[23]提出的潜在结果框架在因果效应估计领域被广泛使用。本文基于潜在结果框架，将待研究的个体记为i,i=\mathrm{1,2},\cdots ,n，个体i的协变量记为{\mathit{x}}_i={\left(x_i^{\mathrm{1}},x_i^{\mathrm{2}},\dots ,x_i^q\right)}^{\mathrm{T}}，处理变量记为T_i，结果变量记为Y_i。本文仅考虑二元处理，即T_i=\left\{t\right\},t=\mathrm{0,1}。t=\mathrm{1}表示接受处理，个体i的两个潜在结果记为Y\left({\mathit{x}}_i|T_i=\mathrm{1}\right)和Y\left({\mathit{x}}_i|T_i=\mathrm{0}\right)，则个体处理效应（Individual Treatment Effect，ITE）为：

\mathrm{I}\mathrm{T}\mathrm{E}=Y\left({\mathit{x}}_i|T_i=\mathrm{1}\right)-Y\left({\mathit{x}}_i|T_i=\mathrm{0}\right)（1）

由于一个个体通常无法同时观测到两个潜在结果，所以个体处理效应无法通过观测数据直接得到。因此已有研究大多通过对个体处理效应求期望得到总体的ATE：

\mathrm{A}\mathrm{T}\mathrm{E}=E\left[\mathrm{I}\mathrm{T}\mathrm{E}\right]=

\begin{array}{c}E\left[Y\left({\mathit{x}}_i|T_i=\mathrm{1}\right)-Y\left({\mathit{x}}_i|T_i=\mathrm{0}\right)\right]=\end{array}

E\left[Y\left({\mathit{x}}_i|T_i=\mathrm{1}\right)\right]-E\left[Y\left({\mathit{x}}_i|T_i=\mathrm{0}\right)\right]（2）

本文中所提到的因果效应均为ATE。

1.2　基于倾向得分逆概率加权的因果效应估计

对ATE进行估计时，需要满足以下假设^[23]：

1) 个体稳定假设；

2) 没有未观测到的混杂变量(无混杂假设)；

3) 每个个体接受处理的概率非零(正性假设)。本文研究也遵循上述假设。

在估计ATE时，基于倾向得分的方法应用最为广泛。倾向得分记为p，第i个个体的倾向得分p_i定义如下：

p_i=P\left(T_i=\mathrm{1}|{\mathit{x}}_i\right)（3）

其含义为在给定协变量x_i的条件下，处理T_i为1的概率。

逻辑回归(Logistic Regression，LR)作为估计倾向得分的基准模型，在医学、社会科学、经济学等众多领域应用普遍，其原理清晰、结果易于解释。更复杂的加权策略，如模型平均、机器学习等方法可能会增加计算负担和结果的不确定性。逻辑回归为TSW方法提供了一个相对简洁、计算高效且易于与分组步骤结合的加权框架，便于验证该方法的有效性。在这种方法中，将处理变量T作为因变量，协变量 x 作为自变量，将LR输出结果作为倾向得分估计值。二元处理变量的倾向得分估计模型可以表示为：

\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(p)=\mathit{\alpha }+{\mathit{x}}^{\mathrm{T}}\mathit{\beta }（4）

其中，\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}为对数函数，\mathit{x}为协变量，\mathit{\alpha }为截距项，\mathit{\beta }={\left({\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},\dots ,{\beta }_q\right)}^{\mathrm{T}}为回归系数，q为协变量数。使用回归系数\mathit{\beta }计算个体i的倾向得分：

p_i=\frac{{\mathrm{e}}^{{\mathit{X}}_i^{\mathrm{T}}\mathit{\beta }}}{\mathrm{1}+{\mathrm{e}}^{{\mathit{X}}_i^{\mathrm{T}}\mathit{\beta }}}（5）

利用倾向得分控制混杂因素的影响，获得ATE的无偏估计^[15]。其中基于倾向得分的逆概率加权方法根据倾向得分估计值生成权重，平衡观测数据的分布。第i个个体的权重w_i计算公式如下：

w_i=\frac{T_i}{p_i}+\frac{\mathrm{1}-T_i}{\mathrm{1}-p_i}（6）

进而，可以利用下式估计ATE：

\widehat{\mathrm{A}\mathrm{T}\mathrm{E}}=\frac{\mathrm{1}}{n}\sum _i{w}_i\bullet Y_i=\frac{\mathrm{1}}{n}\sum _i\frac{T_i\bullet Y_i}{p_i}-\frac{\mathrm{1}}{n}\sum _i\frac{\left(\mathrm{1}-T_i\right)\bullet Y_i}{\mathrm{1}-p_i}（7）

1.3　三阶段加权因果效应估计算法

本研究针对缺失数据集S开展研究，该数据集包含表示个体特征的协变量W、处理变量T和结果变量Y，其中协变量W中存在缺失值。协变量又可以根据是否与处理变量、结果变量存在直接影响关系，分为混杂变量、暴露预测因子和结果预测因子。TSW算法使用子模型加权方法估计ATE，如算法1。

第一阶段 填补混杂变量缺失值

在构建因果分析模型时，首先需要明确处理变量与结果变量之间的因果关系网络。因此根据与处理变量、结果变量的因果关系，构建如图1所示的因果图，将q维原始协变量W分为3组：混杂变量组C_{\mathrm{1}}=\left\{w_{\mathrm{11}},w_{\mathrm{12}},\dots ,w_{\mathrm{1}{q}_{\mathrm{1}}}\right\}，暴露预测因子组C_{\mathrm{2}}=\left\{w_{\mathrm{21}},w_{\mathrm{22}},\dots ,w_{\mathrm{2}{q}_{\mathrm{2}}}\right\}和结果预测因子组C_{\mathrm{3}}=\left\{w_{\mathrm{31}},w_{\mathrm{32}},\dots ,w_{\mathrm{3}{q}_{\mathrm{3}}}\right\}。其中，混杂变量既对处理变量T有直接影响，也对结果变量Y有直接影响；暴露预测因子只对处理变量T有直接影响，对结果变量Y有间接影响；结果预测因子只对结果变量Y有直接影响。

对不同类型的协变量采用不同的缺失值处理办法，对混杂变量使用插补方法补全缺失的协变量值，处理方式见图2；暴露预测因子和结果预测因子则使用缺失模式下倾向得分加权的方法进行处理，处理方式见表1。当原始协变量数据个数较多时，直接根据每个协变量缺失情况划分模式会导致缺失模式数量过多，进而出现某些缺失模式下样本量较少的情况，影响建模分析。因此对原始协变量进行分类后再划分缺失模式，这样一方面控制了缺失模式数量和每个缺失模式下个体的数量；另一方面能给三类协变量赋予不同的影响权重。

在估计ATE时，混杂变量是否完整对于因果效应估计十分重要。为了得到完整的混杂变量，在第一阶段先使用MI多重插补方法填补混杂变量组C_{\mathrm{1}}的缺失值。具体做法如下：

将原始数据集S^{\mathrm{0}}中变量\left\{C_{\mathrm{1}},T,Y\right\}对应的数据子集记为S_{\mathrm{1}}，其中处理变量T和结果变量Y作为辅助变量^[24]；对S_{\mathrm{1}}使用基于链式方程的多重插补方法进行插补，该方法通过迭代预测模型插补数据集中的缺失值，在每次迭代中使用数据集中的其他变量来估算待处理缺失变量，直至满足收敛条件。TSW方法中迭代预测模型采用预测均值匹配和线性回归预测，分别对应基于多重均值插补的三阶段加权因果效应估计算法(TSW-M)和基于多重回归插补的三阶段加权因果效应估计算法(TSW-R)。将混杂变量组C_{\mathrm{1}}的缺失值填补完整后，得到新数据集S。

第二阶段 计算加权倾向得分

1） 划分缺失子集

为了减少插补方法带来的干扰，对C_{\mathrm{2}}和C_{\mathrm{3}}不再使用插补的方法处理缺失值，而是通过划分缺失子集处理缺失。划分缺失子集是指将变量缺失模式相同的个体划分到同一模式中，构成具有相同缺失特征的缺失子集。

首先根据数据集S的缺失情况将其划分为4种缺失模式，每种缺失模式分别对应缺失子集S_r,r=\mathrm{1,2},\mathrm{3,4}，每个缺失子集中包含具有相同缺失模式的个体，见表1。参考朱建平等^[25]对缺失模式进行定义，将缺失函数定义为R\left(x_i,C_m\right)，表示个体i的第m组变量的缺失情况，其中i=\mathrm{1,2},\dots ,n，m=\mathrm{1,2},\mathrm{3}。若存在缺失，则R\left(x_i,C_m\right)=\mathrm{1}，不存在缺失时，R\left(x_i,C_m\right)=\mathrm{0}。对于任意两个个体i和j，如果满足：

R\left(i,C_m\right)=R\left(j,C_m\right),i\ne j（8）

则个体i和j具有相同缺失模式。

缺失子集S_r中包含的个体i_{s_r}都具有相同的缺失模式，缺失子集S_{\mathrm{1}}中所有个体i_{s_{\mathrm{1}}}的协变量均无缺失，因此均满足：

\left\{i_{s_{\mathrm{1}}}:R\left(i_{s_{\mathrm{1}}},C_m\right)=\mathrm{0},m=\mathrm{1,2},\mathrm{3}\right\}（9）

缺失子集S_{\mathrm{2}}中个体i_{s_{\mathrm{2}}}，若C_{\mathrm{2}}组某一个或多个协变量存在缺失，C_{\mathrm{1}}组和C_{\mathrm{3}}组协变量无缺失，则：

\left\{i_{s_{\mathrm{2}}}:R\left(i_{s_{\mathrm{2}}},C_m\right)=\mathrm{0},m=\mathrm{1,3};R\left(i_{s_{\mathrm{2}}},C_m\right)=\mathrm{1},m=\mathrm{2}\right\}（10）

因此个体i_{s_{\mathrm{2}}}保留\left\{C_{\mathrm{1}},C_{\mathrm{3}}\right\}协变量参与后续模型的建立。

缺失子集S_{\mathrm{3}}中个体i_{s_{\mathrm{3}}}，若C_{\mathrm{3}}组某一个或多个协变量存在缺失，C_{\mathrm{1}}组和C_{\mathrm{2}}组协变量无缺失，则：

\left\{i_{s_{\mathrm{3}}}:R\left(x_{s_{\mathrm{3}}},C_m\right)=\mathrm{0},m=\mathrm{1,2};R\left(i_{s_{\mathrm{3}}},C_m\right)=\mathrm{1},m=\mathrm{3}\right\}（11）

因此个体i_{s_{\mathrm{3}}}保留\left\{C_{\mathrm{1}},C_{\mathrm{2}}\right\}协变量参与后续模型的建立。

缺失子集S_{\mathrm{4}}中个体i_{s_{\mathrm{4}}}，若C_{\mathrm{2}}组和C_{\mathrm{3}}组某一个或多个协变量存在缺失，仅有C_{\mathrm{1}}组协变量无缺失，则：

\left\{i_{s_{\mathrm{4}}}:R\left(i_{s_{\mathrm{4}}},C_m\right)=\mathrm{0},m=\mathrm{1};R\left(i_{s_{\mathrm{4}}},C_m\right)=\mathrm{1},m=\mathrm{2,3}\right\}（12）

因此个体i_{s_{\mathrm{4}}}仅保留无缺失的\left\{C_{\mathrm{1}}\right\}协变量参与后续模型的建立。设n_{s_r}为S_r中个体i_{s_r}的数量，则有：

n={\displaystyle \sum _{r=\mathrm{1}}^{\mathrm{4}}}{n}_{s_r}（13）

根据上述定义，将数据集S中的个体按照缺失模式划分到缺失子集S_r(r=\mathrm{1,2},\mathrm{3,4})中。

2） 构建子模型组

在每种缺失模型下构建子模型组，使用子模型加权法给协变量组赋予不同权重。首先将子模型加权方法分别应用于不同缺失模式，针对缺失子集S_r，根据无缺失协变量构建不同变量组合U_{rj}，进而根据不同的变量组合建立子模型组M_r=\left\{m_{r\mathrm{1}},m_{r\mathrm{2}},\dots ,m_{r{k}_r}\right\}，缺失模式下的子模型组如表1所示。

以缺失模式1为例，因为该模式保留了所有协变量组，所以在C_{\mathrm{1}}保持不变的情况下，无缺失的协变量的组合形式共有4种：\left\{C_{\mathrm{1}},C_{\mathrm{2}},C_{\mathrm{3}}\right\}，\left\{C_{\mathrm{1}},C_{\mathrm{3}}\right\}，\left\{C_{\mathrm{1}},C_{\mathrm{2}}\right\}，\left\{C_{\mathrm{1}}\right\}。每个协变量组的组合形式对应一个子模型，因此该模式下缺失子集S_{\mathrm{1}}对应子模型组包含4个子模型，记为\{m_{\mathrm{11}},m_{\mathrm{12}},m_{\mathrm{13}},m_{\mathrm{14}}\}。

3） 计算加权倾向得分

对S_r的子模型组M_r=\left\{m_{r\mathrm{1}},m_{r\mathrm{2}},\dots ,m_{r{k}_r}\right\}，对每个子模型m_{rj}利用LR估计倾向得分{\mathit{P}}_{rj}=(p_{rj\mathrm{1}},p_{rj\mathrm{2}},\dots ,p_{rj{l}_r}{)}^{\mathrm{T}}，则缺失子集S_r的倾向得分矩阵{\mathit{P}}_r：

{\mathit{P}}_r=[{\mathit{P}}_{r\mathrm{1}},{\mathit{P}}_{r\mathrm{2}},\dots ,{\mathit{P}}_{r{k}_r}]=\left[\begin{array}{cccc}{p}_{r\mathrm{11}}& p_{r\mathrm{21}}& \cdots & p_{r{k}_r\mathrm{1}}\\ p_{r\mathrm{12}}& p_{r\mathrm{22}}& \cdots & p_{r{k}_r\mathrm{2}}\\ ⋮& ⋮& \ddots & ⋮\\ p_{r\mathrm{1}{l}_r}& p_{r\mathrm{2}{l}_r}& \cdots & p_{r{k}_r{l}_r}\end{array}\right]

以缺失模式1为例，设该模式下缺失子集S_{\mathrm{1}}有l_{\mathrm{1}}=n_{s_{\mathrm{1}}}个个体，S_{\mathrm{1}}对应子模型组包含子模型\left\{m_{\mathrm{11}},m_{\mathrm{12}},m_{\mathrm{13}},m_{\mathrm{14}}\right\}。子模型m_{\mathrm{11}}使用\left\{C_{\mathrm{1}},C_{\mathrm{2}},C_{\mathrm{3}}\right\}变量建立LR模型计算倾向得分，得到{\mathit{P}}_{\mathrm{11}}=(p_{\mathrm{111}},p_{\mathrm{112}},\dots ,p_{\mathrm{11}{l}_{\mathrm{1}}}{)}^{\mathrm{T}}，子模型m_{\mathrm{12}}使用\left\{C_{\mathrm{1}},C_{\mathrm{3}}\right\}变量建立LR模型计算倾向得分，得到{\mathit{P}}_{\mathrm{12}}=(p_{\mathrm{121}},p_{\mathrm{122}},\dots ,p_{\mathrm{12}{l}_{\mathrm{1}}}{)}^{\mathrm{T}}，以此类推，得到子模型m_{\mathrm{13}}、m_{\mathrm{14}}倾向得分{\mathit{P}}_{\mathrm{13}}和{\mathit{P}}_{\mathrm{14}}，最终得到S_{\mathrm{1}}的倾向得分矩阵{\mathit{P}}_{\mathrm{1}}：

{\mathit{P}}_{\mathrm{1}}=[{\mathit{P}}_{\mathrm{11}},{\mathit{P}}_{\mathrm{12}},\dots ,{\mathit{P}}_{\mathrm{14}}]=\left[\begin{array}{cccc}{p}_{\mathrm{111}}& p_{\mathrm{121}}& \cdots & p_{\mathrm{141}}\\ p_{\mathrm{112}}& p_{\mathrm{122}}& \cdots & p_{\mathrm{142}}\\ ⋮& ⋮& \ddots & ⋮\\ p_{\mathrm{11}{l}_{\mathrm{1}}}& p_{\mathrm{12}{l}_{\mathrm{1}}}& \cdots & p_{\mathrm{14}{l}_{\mathrm{1}}}\end{array}\right]

参考广义线性模型平均方法对子模型组进行加权，设子模型组M_r的权重W_r=\left\{w_{r\mathrm{1}},w_{r\mathrm{2}},\dots ,w_{r{k}_r}\right\}，使用非线性逻辑回归模型对矩阵{\mathit{P}}_r进行加权：

{\mathit{P}}_r^{*}=P\left(T|{\mathit{P}}_{rj}\right)=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left\{w_{rj}{\mathit{P}}_{rj}+\mathit{b}\right\}}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}\left\{w_{rj}{\mathit{P}}_{rj}+\mathit{b}\right\}}（14）

其中，b是截距项。使用逻辑回归模型求解上述方程，用极大似然估计方法进行估计：

l\left(w_{rj},\mathit{b}\right)={\displaystyle \sum _{j=\mathrm{1}}^{k_r}}\mathrm{l}\mathrm{n}P\left(T|{\mathit{P}}_{rj};w_{rj},\mathit{b}\right)（15）

得到长度为l_r的加权倾向得分{\mathit{P}}_r^{*}。对k_r个子模型加权赋予不同缺失模式下的子数据集不同的变量组权重，可使结果更加稳定。

合并全部缺失子集S_r(r=\mathrm{1,2},\mathrm{3,4})的加权倾向得分{\mathit{P}}_r^{*}(r=\mathrm{1,2},\mathrm{3,4})，得到全部样本点总加权倾向得分{\mathit{P}}^{\mathrm{*}}=({\mathit{P}}_{\mathrm{1}}^{\mathrm{*}},{\mathit{P}}_{\mathrm{2}}^{\mathrm{*}},{\mathit{P}}_{\mathrm{3}}^{\mathrm{*}},{\mathit{P}}_{\mathrm{4}}^{\mathrm{*}}{)}^{\mathrm{T}}。

第三阶段 计算ATE

基于加权倾向得分{\mathit{P}}^{*}，使用逆概率加权的方法计算ATE估计。

\mathrm{A}\mathrm{T}\mathrm{E}=\frac{\mathrm{1}}{n}\sum _i\frac{T_i\bullet Y_i}{{\mathit{P}}^{*}\left(T_i=\mathrm{1}|{\mathit{x}}_i\right)}-\frac{\mathrm{1}}{n}\sum _i\frac{\left(\mathrm{1}-T_i\right)\bullet Y_i}{\mathrm{1}-{\mathit{P}}^{*}\left(T_i=\mathrm{1}|{\mathit{x}}_i\right)}（16）

使用加权倾向得分估计ATE值的示意图如图3所示。

2　 蒙特卡罗模拟实验

本节使用蒙特卡罗模拟数据验证TSW算法的估计效果。TSW算法结合均值插补和回归插补，分别构成基于多重均值插补的三阶段加权因果效应估计算法（TSW-M）和基于多重回归插补的三阶段加权因果效应估计算法（TSW-R）。其中，多重回归插补的插补数据集数量设置为5，迭代次数设置为10次。在模拟数据集上将TSW算法与常用处理协变量缺失的插补方法和GBM方法进行比较。

2.1　评价指标

模拟实验使用以下指标进行评价：偏差BIAS，偏差比RBIAS，标准误差SE，均方根误差RMSE。

偏差BIAS计算公式为：

\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}=\left|\overline{\widehat{\gamma }}-\gamma \right|（17）

其中，\overline{\widehat{\gamma }}为所有ATE估计的均值，\gamma为真实ATE值，在本实验中设置为-\mathrm{0.4}。

RBIAS是在BIAS的基础上得到的，计算公式为：

\mathrm{R}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}=\frac{\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}}{\left|\gamma \right|}\times \mathrm{100}\mathrm{\%}（18）

标准误差\mathrm{S}\mathrm{E}表示预测值与其均值之间的差异。\mathrm{S}\mathrm{E}越小，样本的可靠性就越高。标准误差SE由标准差s得来，需要先计算s：

s=\sqrt[]{\frac{\mathrm{1}}{n-\mathrm{1}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\left({\widehat{\gamma }}_i-\overline{\gamma }\right)}^{\mathrm{2}}}

得到s后再计算\mathrm{S}\mathrm{E}，计算公式为：

\mathrm{S}\mathrm{E}=\frac{\mathrm{1}}{\sqrt[]{H}}\bullet s（19）

其中，H为实验重复次数。

均方根误差RMSE用于衡量预测值与真实值之间的偏差，计算每个估计值与真实值之间差异的均方根，计算公式为：

\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt[]{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{\left({\widehat{\gamma }}_i-\gamma \right)}^{\mathrm{2}}}（20）

其中，{\widehat{\gamma }}_i为ATE估计值，\gamma为真实ATE值。

2.2　模拟实验

2.2.1 模拟数据生成

1） 数据模型设置

实验使用的模拟数据和模型场景设置参考Setoguchi等^[26]的设置，变量间的相关关系如图4所示。

图4中\left\{w_{\mathrm{1}},w_{\mathrm{3}},w_{\mathrm{5}},w_{\mathrm{6}},w_{\mathrm{8}},w_{\mathrm{9}}\right\}为二元变量，\left\{w_{\mathrm{2}},w_{\mathrm{4}},w_{\mathrm{7}},w_{\mathrm{10}}\right\}为连续变量。w_{\mathrm{2}}和w_{\mathrm{6}}、w_{\mathrm{3}}和w_{\mathrm{8}}之间存在弱相关关系，相关系数{\rho }_{\mathrm{2,6}}=\mathrm{0.2}，{\rho }_{\mathrm{3,8}}=\mathrm{0.2}，w_{\mathrm{1}}和w_{\mathrm{5}}、w_{\mathrm{4}}和w_{\mathrm{9}}之间存在强相关关系，{\rho }_{\mathrm{1,5}}=\mathrm{0.9}，{\rho }_{\mathrm{4,9}}=\mathrm{0.9}。

处理变量T的设置方式如表2所示，在场景A~G下，使用LR构建真实倾向评分模型得到p\left[T|w_i\right]。然后从均匀分布中生成一个\left(\mathrm{0,1}\right)之间的随机数，如果随机数小于估计的真实倾向评分p\left[T|w_i\right]，T设为1，否则设为0。

结果变量Y对协变量\left\{w_{\mathrm{1}},\dots ,w_{\mathrm{4}},w_{\mathrm{8}},\dots ,w_{\mathrm{10}}\right\}和处理变量T使用LR构造结果模型，如（21）式所示：

\mathrm{P}\mathrm{r}\left[Y=\mathrm{1}|w_i,T_i\right]={\left(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\left(-\mathrm{3.85}+\mathrm{0.3}{w}_{\mathrm{1}}-\mathrm{0.36}{w}_{\mathrm{2}}-\mathrm{0.73}{w}_{\mathrm{3}}-\mathrm{0.2}{w}_{\mathrm{4}}+\mathrm{0.71}{w}_{\mathrm{8}}-\mathrm{0.19}{w}_{\mathrm{9}}+\mathrm{0.26}{w}_{\mathrm{10}}-\mathrm{0.4}T\right)\right\}\right)}^{-\mathrm{1}}（21）

其中，T为二元处理变量；Y为连续结果变量；协变量\left\{w_{\mathrm{1}},w_{\mathrm{2}},\dots ,w_{\mathrm{10}}\right\}可以分为：混杂变量组C_{\mathrm{1}}=\left\{w_{\mathrm{1}},w_{\mathrm{2}},w_{\mathrm{3}},w_{\mathrm{4}}\right\}，暴露预测因子组C_{\mathrm{2}}=\left\{w_{\mathrm{5}},w_{\mathrm{6}},w_{\mathrm{7}}\right\}，结果预测因子组C_{\mathrm{3}}=\left\{w_{\mathrm{8}},w_{\mathrm{9}},w_{\mathrm{10}}\right\}。

再从均匀分布中生成一个\left(\mathrm{0,1}\right)之间的随机数，如果随机数小于p\left[Y|w_i,T_i\right]，Y设为1，否则设为0。

2） 数据缺失机制

得到完整数据后，需要生成具有缺失的数据集。模拟实验使用完全随机缺失机制对完整数据中的\left\{w_{\mathrm{1}},w_{\mathrm{2}},\dots ,w_{\mathrm{10}}\right\}进行处理，完全随机缺失指某变量的缺失数据与其他任何观测或未观测变量都不相关。当样本量为n时，10个协变量共有\mathrm{10}\times n个协变量值，设协变量的缺失率为\theta，则缺失的协变量值个数为\mathrm{10}\times n\times \theta。在A~G共7个场景下，生成样本量n分别为2 000、5 000，缺失率\theta分别为5%、10%、20%的数据集，进行6组模拟实验，每组实验重复H=\mathrm{100}次。

2.2.2 实验结果

在协变量数据缺失的情况下，将TSW-M和TSW-R与单一插值、多重插值、广义提升等方法的结果进行对比。其中，单一均值插值方法记为SI-M，多重均值插值和多重回归插值方法分别记为MI-M和MI-R。

首先将协变量数据无缺失情况下的ATE估计结果作为基线值。表3展示了样本量为2 000时，对无缺失数据进行因果效应估计的模拟结果。在使用LR计算倾向得分、使用逆概率加权计算ATE估计值的情况下，不同场景的偏差值存在很大的区别，A场景下RBIAS最小，仅有1.25%；F场景下RBIAS最大，为32.21%。

表4为当样本量为2 000、缺失率为20%时，不同处理缺失数据方法的ATE估计及其各项指标结果。其中，Complete为基线值，表示对无缺失数据使用LR方法估计得到的结果。TSW方法在A、B、C、E、F、G共6个场景下BIAS和RMSE最小，在A、B、C、F、G共5个场景下SE最小。TSW-M和TSW-R方法相比，指标BIAS稍大，但是指标SE和RMSE表现较好。实验结果显示，TSW类方法明显优于其他方法。

GBM方法仅次于TSW类方法，在场景D上表现最优，SI-M的BIAS总体表现差于MI-M方法，但是优于MI-R方法，这是由于MI-R在不同场景下的BIAS波动较大，稳定性差；并且SI-M方法的SE和RMSE均小于MI类方法，表现更好。

表5为当样本量为5 000、缺失率为20%时，不同处理缺失数据方法的ATE估计及各项指标结果。当样本量增大后，TSW方法的优势依旧明显，在A、B、C、E、F、G共6个场景下获得了最小的BIAS，在B、C、E、F、G共5个场景下获得了最小的SE和RMSE。TSW-M和TSW-R方法相比，BIAS和RMSE稍差，SE表现较好。GBM在D场景下获得最小BIAS、SE和RMSE。其余3种方法与样本量为2 000时表现差别不大，MI类的在BIAS上的表现更好，但MI-R稳定性差；SI的SE和BISA的值更小，表现更好。可以看出，当样本量增加到5 000后，BIAS的变化不大，但SE和RMSE的值更小，ATE估计的效果更好。

鉴于篇幅限制，不再展示样本量分别为2 000和5 000、缺失率分别为5%和10%情形下的全部指标运算结果，仅展示RMSE值，结果如图5所示。在4种不同组合的情况下，TSW-M和TSW-R都能获得较低且更稳定的RMSE值。

总体来看，在不同的模拟数据场景下，TSW类方法的表现总体优于插补类和GBM方法，ATE估计效果更好。由于TSW方法结合倾向得分加权与回归调整，其稳健性随模型误设非线性类型和程度而异，并受数据生成模型中可加性与非线性交互作用的影响。在轻度非线性和中度非线性模型中，TSW因倾向得分与结果模型双重误设而表现不稳定，估计效果下降；但在中度非可加性模型下，其稳定性优势凸显，能抑制高维交互产生的影响。值得注意的是，在中度非可加性与非线性模型中，TSW的精度反而优于轻度模型及中度非可加性模型。这可能源于其机制对误设的非线性响应，轻度复杂度下双重误设敏感导致算法脆弱，高复杂度下权重平均化效应和回归调整潜力使算法相对稳定，而中度非可加性与非线性模型中大量交互项使倾向得分分布更均匀，降低了极端权重的风险，使其稳定性增益部分抵消了回归调整的不足。

3　 结 语

本文主要研究了基于协变量数据缺失下的因果效应估计方法，提出了三阶段加权因果效应估计算法，提升了估计ATE的效果。三阶段加权因果效应估计算法是一种新的协变量缺失的因果效应估计方法，在使用观测数据进行因果效应估计时，减少了混杂因素和数据缺失带来的影响。该方法在保留原始样本信息的情况下，仅对混杂变量使用插补方法补全缺失，减少插补数据带来的偏差，得到更好的ATE估计值。模拟实验验证了三阶段加权因果效应估计算法，相较于其他插值和GBM处理缺失值的方法，估计效果更好，总体表现也更好。

本文提出的算法扩充了对缺失数据进行因果效应估计的研究，此算法是基于数据完全随机缺失模式对因果效应展开估计的。数据非随机缺失的情况更为复杂，未来的研究会针对非随机缺失情况进一步探索。

参考文献

[1] 范菊逸， 詹铭峰， 蔡宗武， 等.带有变量选择的协变量平衡倾向得分的估计： 基于GMM-LASSO方法［J］.系统工程理论与实践 2021， 41（10）：2631-2639. DOI：10.12011/SETP2020-0037 .

[2] FAN J Y， ZHAN M F， CAI Z W， et al. Covariate balancing in propensity score estimation with variable selection： Based on GMM-LASSO approach［J］. Systems Engineering-Theory & Practice， 2021， 41（10）：2631-2639. DOI：10.12011/SETP2020-0037（Ch ）.

[3] BOTTOU L， PETERS J， QUIÑONERO-CANDELA J， et al. Counterfactual reasoning and learning systems： The example of computational advertising［J］ Journal of Machine Learning Research，2013， 14（1）： 3207-3260. DOI：10.5555/2567709.2567766 .

[4] ROSENBAUM P R， RUBIN D B. The central role of the propensity score in observational studies for causal effects［J］. Biometrika， 1983， 70（1）： 41-55. DOI: 10.1093/biomet/70.1.41 .

[5] 蒋青嬗， 马佳羽， 黄灿， 等. 因果推断中基于能源距离的协变量分布平衡［J］. 统计研究， 2023， 40（5）： 144-151. DOI：10.19343/j.cnki.11-1302/c.2023.05.011 .

[6] JIANG Q S， MA J Y， HUANG C， et al. Covariate distribution balance via energy distance for causal inference［J］. Statistical Research， 2023， 40（5）： 144-151. DOI：10.19343/j.cnki.11-1302/c.2023.05.011（Ch ）.

[7] ROBINS J M， ROTNITZKY A， ZHAO L P. Estimation of regression coefficients when some regressors are not always observed［J］. Journal of the American Statistical Association， 1994， 89（427）： 846-866. DOI：10.1080/01621459.1994.10476818 .

[8] DUDÍK M， LANGFORD J， LI L H. Doubly robust policy evaluation and learning［EB/OL］.［2025-06-07］. DOI: 10.1214/14-sts500 .

[9] RODRIGUEZ DUQUE D， STEPHENS D A， MOODIE E E M， et al. Semiparametric Bayesian inference for optimal dynamic treatment regimes via dynamic marginal structural models［J］. Biostatistics， 2023， 24（3）： 708-727. DOI：10.1093/biostatistics/kxac007 .

[10] 刘冲， 沙学康， 张妍. 交错双重差分： 处理效应异质性与估计方法选择［J］. 数量经济技术经济研究， 2022， 39（9）： 177-204. DOI：10.13653/j.cnki.jqte.20220805.001 .

[11] LIU C， SHA X K， ZHANG Y. Staggered difference-in-differences method： Heterogeneous treatment effects and choice of estimation［J］. Journal of Quantitative & Technological Economics， 2022， 39（9）： 177-204. DOI：10.13653/j.cnki.jqte.20220805.001（Ch ）.

[12] ANGRIST J D， IMBENS G W， RUBIN D B. Identification of causal effects using instrumental variables［J］. Journal of the American Statistical Association， 1996， 91（434）： 444-455. DOI：10.1080/01621459.1996.10476902 .

[13] ANDERSON M， MAGRUDER J. Learning from the crowd： Regression discontinuity estimates of the effects of an online review database［J］. The Economic Journal， 2012， 122（563）： 957-989. DOI：10.1111/j.1468-0297.2012.02512.x .

[14] MIAO W， TCHETGEN TCHETGEN E. Invited commentary： Bias attenuation and identification of causal effects with multiple negative controls［J］. American Journal of Epidemiology， 2017， 185（10）： 950-953. DOI：10.1093/aje/kwx012 .

[15] LI H K， MIAO W， CAI Z， et al. Causal data fusion methods using summary-level statistics for a continuous outcome［J］. Statistics in Medicine， 2020， 39（8）： 1054-1067. DOI：10.1002/sim.8461 .

[16] JOSEY K P， YANG F， GHOSH D， et al. A calibration approach to transportability and data-fusion with observational data［J］. Statistics in Medicine， 2022， 41（23）： 4511-4531. DOI：10.1002/sim.9523 .

[17] YANG S， DING P. Combining multiple observational data sources to estimate causal effects［J］. Journal of the American Statistical Association， 2020， 115（531）： 1540-1554. DOI：10.1080/01621459.2019.1609973 .

[18] WHITE I R， CARLIN J B. Bias and efficiency of multiple imputation compared with complete-case analysis for missing covariate values［J］. Statistics in Medicine， 2010， 29（28）： 2920-2931. DOI：10.1002/sim.3944 .

[19] ROSENBAUM P R， RUBIN D B. Reducing bias in observational studies using subclassification on the propensity score［J］. Journal of the American Statistical Association， 1984， 79（387）： 516-524. DOI：10.1080/01621459.1984.10478078 .

[20] D’AGOSTINO R， LANG W， WALKUP M， et al. Examining the impact of missing data on propensity score estimation in determining the effectiveness of self-monitoring of blood glucose （SMBG）［J］. Health Services and Outcomes Research Methodology， 2001， 2（3）： 291-315. DOI：10.1023/A： 1020375413191 .

[21] RUBIN D B. Multiple Imputation for Nonresponse in Surveys［M］. Hoboken： Wiley， 1987. DOI：10.1002/9780470316696 .

[22] 邓建新， 单路宝， 贺德强， 等. 缺失数据的处理方法及其发展趋势［J］. 统计与决策， 2019， 35（23）： 28-34. DOI：10.13546/j.cnki.tjyjc.2019.23.005 .

[23] DENG J X， SHAN L B， HE D Q， et al. Processing method of missing data and its developing tendency［J］. Statistics & Decision， 2019， 35（23）： 28-34. DOI：10.13546/j.cnki.tjyjc.2019.23.005（Ch ）.

[24] GRAHAM J W. Missing data analysis： Making it work in the real world［J］. Annual Review of Psychology， 2009， 60： 549-576. DOI：10.1146/annurev.psych.58.110405.085530 .

[25] MCCAFFREY D F， GRIFFIN B A， ALMIRALL D， et al. A tutorial on propensity score estimation for multiple treatments using generalized boosted models［J］. Statistics in Medicine， 2013， 32（19）： 3388-3414. DOI：10.1002/sim.5753 .

[26] WANG Z Q， TANG N S. Bayesian quantile regression with mixed discrete and nonignorable missing covariates［J］. Bayesian Analysis， 2020， 15（2）： 579-604. DOI：10.1214/19-ba1165 .

[27] RUBIN D B. Estimating causal effects of treatments in randomized and nonrandomized studies［J］. Journal of Educational Psychology， 1974， 66（5）： 688-701. DOI：10.1037/h0037350 .

[28] CROWE B J， LIPKOVICH I A， WANG O H. Comparison of several imputation methods for missing baseline data in propensity scores analysis of binary outcome［J］. Pharmaceutical Statistics， 2010， 9（4）： 269-279. DOI：10.1002/pst.389 .

[29] 朱建平， 郑陈璐， 方匡南. 缺失数据下的两阶段信用评分模型——基于互联网消费金融数据的研究［J］. 数理统计与管理， 2021， 40（4）： 613-624. DOI：10.13860/j.cnki.sltj.20201219-006 .

[30] ZHU J P， ZHENG C L， FANG K N. Two-stage credit scoring model with missing data： Based on the analysis of Internet consumer credit data［J］. Journal of Applied Statistics and Management， 2021， 40（4）： 613-624. DOI：10.13860/j.cnki.sltj.20201219-006（Ch ）.

[31] SETOGUCHI S， SCHNEEWEISS S， BROOKHART M A， et al. Evaluating uses of data mining techniques in propensity score estimation： A simulation study［J］. Pharmacoepidemiology and Drug Safety， 2008， 17（6）： 546-555. DOI：10.1002/pds.1555 .

基金资助

国家社会科学基金(23ATJ005)

湖北省教育厅科学研究计划项目(D20222202)

数字金融创新湖北省重点实验室开放基金(DFIK2024Y06)