1. 背景和动机

因果推断引起了互联网头部企业的数据团队重视，并积极招聘相关领域的统计学家，数据科学家，机器学习科学家。无偏因果估计在实际场景中比较困难，但确是保证结论正确的必须任务。关于在互联网行业环境中使用 A/B 测试，与“理想化随机试验”相对应的 A/B 测试对比，注意到随机试验设计的观点通常过于简化。对于大规模进行的 A/B 测试，尤其是长时间进行的测试，很少能满足理想化随机试验的假设。

例如，某部门评估新功能设计对用户体验的影响。分析师任务是分析干预措施 Z（新功能设计）对结果 Y（点击率或其他成功指标）的平均因果效应的无偏估计。我们可以通过 IP 地址将用户随机分配到一个版本的 Z 来设计 A/B 测试（体验新功能设计与不体验新功能设计），进行实验以收集结果 Y 数据。

根据观察到的数据，我们将得到 Z 对 Y 的平均边际效应差异的经验估计。但这里可能直接假设 A/B 测试的“真正因果结构”对应于以下因果有向无环图（DAG）：

假设边际估计是 Z 对 Y 的因果效应差异均值的无偏估计：

在实际的 A/B 测试中，很少有上述的简单因果关系。大多数 A/B 测试的真正因果结构通常看起来像这样：

Z 为随机分组的干预变量，A 为干预的可观察变量，Y 为结果变量，L 为 A 到 Y 关系链路的中介变量，C 为结果变量 Y 删失信息后的情况，U1 和 U2 为不可观察的混淆变量。在以上 DAG 中，从 Z -> A -> [L] <- U2 -> Y 的开放路径，随机干预分配变量 Z 和结果变量 Y 之间显然没有 d-separation。因此，从未删失的数据信息中的条件均值直接估计 Z 对 Y 的因果效应是错误的。

在本文中，我们将给出示例，对应上述复杂的 DAG 进行大规模的长期 A/B 测试，基于逆概率加权的边际结构模型如何解决信息删失的调整。我们还将进行计算仿真以研究这些方法的效果。

2.样本不依从，干预交叉和治疗意向估计（ITT）

同样，从理想化的随机试验开始，我们具有以下因果 DAG：

鉴于 Z 是随机干预变量，要认识到 Z 实际上只是“随机干预分配”。可能有一小部分实验用户不遵守其随机分配，要么不参与实验，要么体验另一组实验的环境。所以在这一变化下我们观察到的干预是 A，其中 Z 是其直接原因。并给予一个未随机且影响 A 和 Y 的不可测量的因素 U1。我们的因果 DAG 如下：

我们对干预分配 Z 对结果 Y 的因果关系并不感兴趣（ITT估计）。而是想关注真正接受了干预 A 样本的结果变量 Y 的表现，即ATE估计^[1]

给定 A 和 Y 在上面的因果 DAG 中由于未测量的共同影响因子 U1 不能被 d-separation，我们该怎么办？我们合理地有两个选择：请注意，随机干预分配 Z 实际上是 A 对 Y 因果关系的工具变量（IV）。假设我们的用例满足一定的同质性条件，我们可以指定一个 IV 变量，以实现 A 对 Y 的平均因果效应的无偏估计。识别到我们的问题是优效试验，随机干预分配 Z 在 Y 上的平均因果效应是意向治疗（ITT）估计。在采用二元干预的优效试验的背景下，并且在没有信息删失的假设下，ITT 估计的期望值为“null-preserving”。从数学上讲，这意味着 Z 对 Y 的平均因果效应的期望值估计将始终比 A 对 Y 的真实平均因果效应的期望值更接近零值。因此，通过估计 Z 对 Y 的平均因果效应，我们还间接计算了感兴趣的 A 对 Y 的平均因果效应的下界。重要的是，由于 ITT 估计是“null-preserving”，因此如果A 对 Y 的真实平均因果效应期望值为零值，则 Z 与 Y 之间平均因果效应期望值的估计也是零值。

这里得出结论，可以通过估计 Z 对 Y 的平均因果效应来确定 A 对 Y 的真实平均因果效应的下界，解决了不可观测的常见因素比如 U2 的影响。但是，我们还没有解决信息删失的问题。回想一下我之前所说的，ITT 估计仅在没有信息删失，或者针对信息删失进行了调整的前提下有效。这是我们接下来要解决的问题。

3. 边际结构建模如何调整信息删失场景下的 A/B 测试结论分析

假定有下面的删失 DAG 模型：

• 变量 L 是观察结果 A 的直接决定因素

• 有不可观测的常见因素 U2 影响 L和结果 Y

• 变量 L 是变量 C 的直接唯一原因，在分析中被“删失”。所谓“删失”，是指失去了后续的数据记录。对于删失组（C = 1），我们没有观察并测量它们的结果 Y。我们仅观察C = 0 的受试者的结果信息，仅限于进行该人群的分析。

这里给出一个该场景下的 A/B 测试案例，假设在部署我们的 A/B 测试之前，我们没有彻底检查新功能是否已正确配置或与所有流行的浏览器兼容。假设 Internet Explorer 的默认设置无法使用该新功能。Explorer 用户首先需要删除默认的防火墙块（由变量 L 表示）。Explorer 的主要用户年龄较大（未测量的变量 U2），这也是结果 Y（点击率）的直接决定因素。因此，即使我们白名单分配了一部分用户可使用新功能（由变量 Z 表示），但许多 Explorer 用户根本没有经历过干预 A，我们也没有测量他们的结果 Y，并且他们的数据信息从分析中被“删失”（C = 1）。我们称这种结构为信息删失。

现在，让我们讨论如何针对信息删失问题进行调整。如本文前面所述，由于从 Z -> A -> [L] <- U2 -> Y 的路径，在以上 DAG 中，随机干预分配Z和结果Y显然没有 d-separation。因此，我们有什么选择？

在第2节中，ITT 估计解决了 U1 这类不能观察到的混淆变量问题，但未解决信息删失问题，IV 估计同理。

现在，让我们讨论如何在给定条件下针对信息删失进行无偏因果效应估计的分析调整。如果我们已经在 U1 和 U2 上记录了数据，则可以使用 standardizing 方法对这些变量进行调整。但是，由于 U1 和 U2 未测量，因此standardizing 无法估计。另外，在现实场景中，我们只是知道 U1、U2 存在，不知道它们代表什么。

使用 Standardization 方法或者使用基于逆概率加权的边际结构模型进行调整，但是这里因为信息删失无法观察到样本删失的结果值 Y，因此这里 Standardization 方法不适用，但是，只要我们观察到所有受试者干预的数据，就可以将 IPW 与边际结构建模一起利用。我们可以通过指定以下因果 SWIG（Single World Intervention Graph）来恢复每个观察样本的 IPW：

推导过程：

先建立预测模型，通过模型计算每个观察样本的倾向得分 ：

表示已删失样本的倾向得分值。

由于我们只想在所有未删失的观测值（c = 0）的 pseudo-population 中估计真实总体的平均因果效应，因此我们为未删失的观测值指定了 IPW，对于所有经过删失的观测值指定了 IPW 等于 0（C = 1）：

为观察数据集中的每个样本计算适当的 IPW，则可以通过将每个观察的加权值与对应的 IPW 加权来生成 pseudo-population。该 pseudo-population 将对应于上图中的因果 DAG 关系。然后，我们通过加权最小二乘回归（WLS）将以下边际结构模型拟合到 pseudo-population 上，来分析 Z 对 Y 影响效应的无偏估计：

对于参数 ，即是在给定的未删失观察数据下，对于处理组 和对照组 之间因果效应差异均值的无偏估计：

根据第2节中得出的关系：

因此，A 对 Y 总体的因果效应估计，由此得到了因果效应差异均值无偏估计的下界值：

接下来进行计算仿真看看模型效果。

4. 边际结构模型的实例计算仿真

我们将在 Python 中进行计算仿真，以分析信息删失的 A/B 测试分析结果是否会受到影响，并查看通过边际结构模型进行调整后的结果差异。

创建具有真实因果 DAG 的模拟数据集，如下图所示，创建两个数据集，一是 A 对 Y 的真实因果效应差为 0，二是 A 对 Y 的真实因果效应差为 0.11。

导入所需的库，并模拟数据集，这里生成函数作为指定模拟数据集用的 A 对 Y 的实际因果效应差异，先设置这个场景 A 对 Y 的真实因果效应差为零，即实验组和对照组之间没有显著差异。

######################
## import libraries ##
######################
import numpy as np
np.random.seed(1081565)
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smfdef simulate_dataset(n=2000000, A_on_Y=0):## Z: “randomized” binary intervention assignment (0, 1)    ## A: observed binary intervention (0, 1)## L: direct binary determinant of intervention A (0, 1)## U1: binary unmeasured common cause of A and Y## U2: binary unmeasured common cause of L and Y## C: censoring variable## Y: countinuous outcome## specify dataframedf = pd.DataFrame()## specify variables Z, U1, and U2Z_split = 0.4U1_split = 0.52U2_split = 0.23df['Z'] = np.random.choice([0, 1], size=n, replace=True, p=[Z_split, (1-Z_split)])df['U1'] = np.random.choice([0, 1], size=n, replace=True, p=[U1_split, (1-U1_split)])df['U2'] = np.random.choice([0, 1], size=n, replace=True, p=[U2_split, (1-U2_split)])## specify variable Alambda_0 = -5.2lambda_1 = 9.3lambda_2 = 3.45   df['A'] = 0df.loc[(df['Z']==0) & (df['U1']==0), 'A'] = np.random.binomial(1, (1/(1+np.exp(-lambda_0))), size=df.loc[(df['Z']==0) & (df['U1']==0), 'A'].shape[0])df.loc[(df['Z']==1) & (df['U1']==0), 'A'] = np.random.binomial(1, (1/(1+np.exp(-(lambda_0+lambda_1)))), size=df.loc[(df['Z']==1) & (df['U1']==0), 'A'].shape[0])df.loc[(df['Z']==0) & (df['U1']==1), 'A'] = np.random.binomial(1, (1/(1+np.exp(-(lambda_0+lambda_2)))), size=df.loc[(df['Z']==0) & (df['U1']==1), 'A'].shape[0])df.loc[(df['Z']==1) & (df['U1']==1), 'A'] = np.random.binomial(1, (1/(1+np.exp(-(lambda_0+lambda_1+lambda_2)))), size=df.loc[(df['Z']==1) & (df['U1']==1), 'A'].shape[0])#model = smf.glm('A ~ Z + U1', data=df, family=sm.families.Binomial()).fit()#print(model.summary())#pd.crosstab(df['Z'], df['A'])#pd.crosstab(df['Z'], df['A'], normalize=True)## specify variable LBeta_0 = -0.52Beta_1 = 2.32Beta_2 = 1.98df['L'] = 0df.loc[(df['A']==0) & (df['U2']==0), 'L'] = np.random.binomial(1, (1/(1+np.exp(-Beta_0))), size=df.loc[(df['A']==0) & (df['U2']==0), 'L'].shape[0])df.loc[(df['A']==1) & (df['U2']==0), 'L'] = np.random.binomial(1, (1/(1+np.exp(-(Beta_0+Beta_1)))), size=df.loc[(df['A']==1) & (df['U2']==0), 'L'].shape[0])df.loc[(df['A']==0) & (df['U2']==1), 'L'] = np.random.binomial(1, (1/(1+np.exp(-(Beta_0+Beta_2)))), size=df.loc[(df['A']==0) & (df['U2']==1), 'L'].shape[0])df.loc[(df['A']==1) & (df['U2']==1), 'L'] = np.random.binomial(1, (1/(1+np.exp(-(Beta_0+Beta_1+Beta_2)))), size=df.loc[(df['A']==1) & (df['U2']==1), 'L'].shape[0])#model = smf.glm('L ~ A + U2', data=df, family=sm.families.Binomial()).fit()#print(model.summary())## specify variable Cpsi_0 = 0.15psi_1 = -5.8df['C'] = 0df.loc[df['L']==0, 'C'] = np.random.binomial(1, (1/(1+np.exp(-psi_0))), size=df.loc[df['L']==0, 'C'].shape[0])df.loc[df['L']==1, 'C'] = np.random.binomial(1, (1/(1+np.exp(-(psi_0+psi_1)))), size=df.loc[df['L']==1, 'C'].shape[0])#model = smf.glm('C ~ L', data=df, family=sm.families.Binomial()).fit()#print(model.summary())## specify variable Ytheta_0 = -0.5theta_1 = 5.78theta_2 = A_on_Ytheta_3 = 1.42df['Y'] = theta_0 + (theta_1*df['U2']) + (theta_2*df['A']) + (theta_3*df['U1']) + np.random.normal(0, 0.5, df.shape[0])df_observed = df[['Z', 'A', 'L', 'C', 'Y']].copy()df_observed.loc[df['C']==1, 'Y'] = Nonereturn(df, df_observed)#################################################################
## simulate dataset with no causal effect difference of A on Y ##
#################################################################
df, df_observed = simulate_dataset(n=2000000, A_on_Y=0)

打印模拟数据集的前几个观察结果：

#####################
## the "true" data ##
#####################
print(df.head(n=10))

但在真实的场景中我们无法测量变量 U1 和 U2 的值，所以这里模拟真实场景的数据集如下：

##############################################
## data we get to "observe" in our analysis ##
##############################################
print(df_observed.head(n=10))

如上所示，在我们的“现实世界”观察数据集中，我们无法得到变量 U1 和 U2 的数据，也无法观察删失后的结果数据Y（C = 1）。先来看看删失数据的百分比：

#####################################################
## print percentage of censored observations (c=1) ##
#####################################################
print(df_observed['C'].value_counts(normalize=True))

如上所示，大约有8％的数据集被删失（C=1），即这部分样本没有观察到结果Y。接下来，让我们用本文中描述的逆概率权重（IPW）边际结构模型来拟合数据：

def get_IPWs(df):PS_model = smf.glm('C ~ L', data=df, family=sm.families.Binomial()).fit()df['PS'] = PS_model.predict(df)df['IPW'] = 0df.loc[df['C']==1, 'IPW'] = 0df.loc[df['C']==0, 'IPW'] = 1 / (1 - df.loc[df['C']==0, 'PS'])return(df)######################################################################
## get Inverse Probability Weights for Marginal Structural Modeling ##
######################################################################
df_observed = get_IPWs(df_observed)

现在先用我们常规中用到的普通回归分析方法，直接在未删失数据上估计，这里强调下我们确定自己模拟的此数据集，A 对 Y 的真实因果效应差异为零。因此，让我们看看常规的分析返回什么：

###################################
## ITT estimator of Z on Y (not adjusting for informative censoring) ##
###################################
model = smf.ols(formula='Y ~ Z', data=df_observed).fit()
print(model.summary())

如下所示，我们的普通回归分析得出 Z 对 Y 的显著因果差异为 -0.2035，置信区间为（-0.211，-0.196）。置信区间不包含 0，说明有显著差异。当然不同于我们模拟设置的真实差异（等于0）。说明我们仅仅忽略8％的删失观测值出现了错误的结论！现在，让我们拟合边际结构模型，返回的结果如下：

###################################
## ITT estimator of Z on Y, adjusting for informative censoring with MSM ##
###################################
MSM_model = smf.wls('Y ~ Z', data=df_observed, weights=np.array(df_observed['IPW'], dtype=np.float64)).fit()
print(MSM_model.summary())

现在，我们找到 Z 对 Y 因果效应的无偏估计，点估计为 0.001，置信区间为（-0.007，0.009），置信区间包含 0，说明无显著差异。这实际上是正确的结果！

其中在信息删失的情况下，当 A 对 Y 的真正因果效应为 0 时， 边际结构模型正确返回了 0 。当 A 对 Y 的真实因果差异不为 0 时该怎么办？让我们再模拟一个数据集，其中 A 对 Y 的真实因果差异为 0.11：