CausalML
是一个 Python 包,它使用基于最新研究的机器学习算法提供一套提升建模和因果推理方法。它提供了一个标准界面,允许用户从实验或观察数据中估计条件平均处理效应 (CATE),也称为个体治疗效应 (ITE)。从本质上讲,它估计了干预 W 对具有观察到特征 X 的用户的结果 Y 的因果影响,而无需对模型形式进行强有力的假设。
在本节中,我们将更深入地探究CausalML中所实现的算法。
我们采用内曼 - 鲁宾潜在结果框架,并假定Y代表结果,W代表处理分配, XiX_iXi 代表观测到的协变量。
支持的算法
CausalML
目前支持以下方法:
- 基于树的算法
- 基于KL散度、欧几里得距离和卡方的提升随机森林
- 基于上下文处理选择的提升随机森林
- 基于delta-delta-p(ΔΔP)准则的提升随机森林(仅适用于二叉树和二分类问题)
- 基于IDDP的提升随机森林(仅适用于二叉树和二分类问题)
- 交互树(仅适用于二叉树和二分类问题)
- 因果推断树(仅适用于二叉树和二分类问题)
- 元学习器算法
- S - 学习器
- T - 学习器
- X - 学习器
- R - 学习器
- 双稳健(DR)学习器
- 工具变量算法
- 两阶段最小二乘法(2SLS)
- 双稳健工具变量(DRIV)学习器
- 基于神经网络的算法
- CEVAE
- DragonNet
- 处理优化算法
- 反事实单元选择
- 反事实价值估计器
元学习器算法
元算法(或元学习器)是一种框架,它利用任意机器学习估计器(称为基学习器)来估计条件平均处理效应(CATE)。
元算法要么使用单个基学习器,并将处理指示变量作为一个特征(例如S - 学习器),要么为处理组和对照组分别使用多个基学习器(例如T - 学习器、X - 学习器和R - 学习器)。
S - 学习器
S - 学习器使用单个机器学习模型按如下方式估计处理效应:
- 阶段1:使用机器学习模型,借助协变量X和处理指示变量W来估计平均结果μ(x)\mu(x)μ(x):
μ(x,w)=E[Y∣X=x,W=w]\mu(x, w)=E[Y|X = x, W = w]μ(x,w)=E[Y∣X=x,W=w] - 阶段2:将CATE估计定义为:
τ^(x)=μ^(x,W=1)−μ^(x,W=0)\hat{\tau}(x)=\hat{\mu}(x, W = 1)-\hat{\mu}(x, W = 0)τ^(x)=μ^(x,W=1)−μ^(x,W=0)
在模型中纳入倾向得分可减少由正则化引发的混杂偏差。
当对照组和处理组在协变量方面差异很大时,单个线性模型不足以对对照组和处理组特征的不同相关维度及平滑度进行编码。
T - 学习器
T - 学习器由以下两个阶段组成:
- 阶段1:使用机器学习模型估计平均结果μ0(x)\mu_0(x)μ0(x)和μ1(x)\mu_1(x)μ1(x):
μ0(x)=E[Y(0)∣X=x]\mu_0(x)=E[Y(0)|X = x]μ0(x)=E[Y(0)∣X=x]
μ1(x)=E[Y(1)∣X=x]\mu_1(x)=E[Y(1)|X = x]μ1(x)=E[Y(1)∣X=x] - 阶段2:将CATE估计定义为:
τ^(x)=μ^1(x)−μ^0(x)\hat{\tau}(x)=\hat{\mu}_1(x)-\hat{\mu}_0(x)τ^(x)=μ^1(x)−μ^0(x)
X - 学习器
X - 学习器是T - 学习器的扩展,由以下三个阶段组成:
- 阶段1:使用机器学习模型估计平均结果μ0(x)\mu_0(x)μ0(x)和μ1(x)\mu_1(x)μ1(x):
μ0(x)=E[Y(0)∣X=x]\mu_0(x)=E[Y(0)|X = x]μ0(x)=E[Y(0)∣X=x]
μ1(x)=E[Y(1)∣X=x]\mu_1(x)=E[Y(1)|X = x]μ1(x)=E[Y(1)∣X=x] - 阶段2:基于μ0(x)\mu_0(x)μ0(x)为处理组中的用户i估算用户级处理效应Di1D_{i1}Di1,基于μ1(x)\mu_1(x)μ1(x)为对照组中的用户j估算用户级处理效应Dj0D_{j0}Dj0:
Di1=Yi1−μ^0(Xi1)D_{i1}=Y_{i1}-\hat{\mu}_0(X_{i1})Di1=Yi1−μ^0(Xi1)
Di0=μ^1(Xi0)−Yi0D_{i0}=\hat{\mu}_1(X_{i0})-Y_{i0}Di0=μ^1(Xi0)−Yi0
然后使用机器学习模型估计τ1(x)=E[D1∣X=x]\tau_1(x)=E[D_1|X = x]τ1(x)=E[D1∣X=x]和τ0(x)=E[D0∣X=x]\tau_0(x)=E[D_0|X = x]τ0(x)=E[D0∣X=x]。
- 阶段3:通过τ1(x)\tau_1(x)τ1(x)和τ0(x)\tau_0(x)τ0(x)的加权平均来定义CATE估计:
τ(x)=g(x)τ0(x)+(1−g(x))τ1(x)\tau(x)=g(x)\tau_0(x)+(1 - g(x))\tau_1(x)τ(x)=g(x)τ0(x)+(1−g(x))τ1(x)
其中g∈[0,1]g\in[0, 1]g∈[0,1]。我们可以将倾向得分用于g(x)g(x)g(x)。
R - 学习器
R - 学习器使用结果的交叉验证折外估计m^(−i)(xi)\hat{m}^{(-i)}(x_i)m^(−i)(xi)和倾向得分e^(−i)(xi)\hat{e}^{(-i)}(x_i)e^(−i)(xi)。它由以下两个阶段组成:
- 阶段1:使用交叉验证通过机器学习模型拟合m^(x)\hat{m}(x)m^(x)和e^(x)\hat{e}(x)e^(x)。
- 阶段2:通过最小化R - 损失L^n(τ(x))\hat{L}^n(\tau(x))L^n(τ(x))来估计处理效应:
L^n(τ(x))=1n∑i=1n((Yi−m^(−i)(Xi))−(Wi−e^(−i)(Xi))τ(Xi))2\hat{L}^n(\tau(x))=\frac{1}{n}\sum_{i = 1}^{n}((Y_i-\hat{m}^{(-i)}(X_i))-(W_i-\hat{e}^{(-i)}(X_i))\tau(X_i))^2L^n(τ(x))=n1i=1∑n((Yi−m^(−i)(Xi))−(Wi−e^(−i)(Xi))τ(Xi))2
其中e^(−i)(Xi)\hat{e}^{(-i)}(X_i)e^(−i)(Xi)等表示未使用第i个训练样本做出的折外留出预测。
双稳健(DR)学习器
DR - 学习器通过在两个阶段交叉拟合双稳健得分函数来估计CATE,步骤如下。我们首先将数据{Y,X,W}\{Y, X, W\}{Y,X,W}随机划分为3个分区{Yi,Xi,Wi}\{Y_i, X_i, W_i\}{Yi,Xi,Wi},i={1,2,3}i = \{1, 2, 3\}i={1,2,3}。
- 阶段1:使用{X1,W1}\{X_1, W_1\}{X1,W1}通过机器学习拟合倾向得分模型e^(x)\hat{e}(x)e^(x),并使用{Y2,X2,W2}\{Y_2, X_2, W_2\}{Y2,X2,W2}通过机器学习为已处理和未处理用户拟合结果回归模型m^0(x)\hat{m}_0(x)m^0(x)和m^1(x)\hat{m}_1(x)m^1(x)。
- 阶段2:使用{Y3,X3,W3}\{Y_3, X_3, W_3\}{Y3,X3,W3}通过机器学习从以下伪结果拟合CATE模型τ^(X)\hat{\tau}(X)τ^(X):
ϕ=W−e^(X)e^(X)(1−e^(X))(Y−m^W(X))+m^1(X)−m^0(X)\phi=\frac{W-\hat{e}(X)}{\hat{e}(X)(1 - \hat{e}(X))}(Y-\hat{m}_W(X))+\hat{m}_1(X)-\hat{m}_0(X)ϕ=e^(X)(1−e^(X))W−e^(X)(Y−m^W(X))+m^1(X)−m^0(X) - 阶段3:再次重复阶段1和阶段2两次。首先使用{Y2,X2,W2}\{Y_2, X_2, W_2\}{Y2,X2,W2}、{Y3,X3,W3}\{Y_3, X_3, W_3\}{Y3,X3,W3}和{Y1,X1,W1}\{Y_1, X_1, W_1\}{Y1,X1,W1}分别用于倾向得分模型、结果模型和CATE模型。然后使用{Y3,X3,W3}\{Y_3, X_3, W_3\}{Y3,X3,W3}、{Y2,X2,W2}\{Y_2, X_2, W_2\}{Y2,X2,W2}和{Y1,X1,W1}\{Y_1, X_1, W_1\}{Y1,X1,W1}分别用于倾向得分模型、结果模型和CATE模型。最终的CATE模型是这3个CATE模型的平均值。
双稳健工具变量(DRIV)学习器
我们将DR - 学习器的思路与局部平均处理效应(LATE)的双稳健得分函数相结合,以估计条件LATE。为此,我们首先将数据{Y,X,W,Z}\{Y, X, W, Z\}{Y,X,W,Z}随机划分为3个分区{Yi,Xi,Wi,Zi}\{Y_i, X_i, W_i, Z_i\}{Yi,Xi,Wi,Zi},i={1,2,3}i = \{1, 2, 3\}i={1,2,3}。
- 阶段1:使用{X1,W1,Z1}\{X_1, W_1, Z_1\}{X1,W1,Z1}为已分配和未分配用户拟合倾向得分模型e^0(x)\hat{e}_0(x)e^0(x)和e^1(x)\hat{e}_1(x)e^1(x),并使用{Y2,X2,Z2}\{Y_2, X_2, Z_2\}{Y2,X2,Z2}通过机器学习为已分配和未分配用户拟合结果回归模型m^0(x)\hat{m}_0(x)m^0(x)和m^1(x)\hat{m}_1(x)m^1(x)。分配概率pZp_ZpZ既可以由用户提供,也可以来自一个简单模型,因为在大多数用例中,分配从设计上就是随机的。
- 阶段2:使用{Y3,X3,W3}\{Y_3, X_3, W_3\}{Y3,X3,W3}通过最小化以下损失函数,通过机器学习拟合条件LATE模型τ^(X)\hat{\tau}(X)τ^(X):
L(τ^(X))=E^[(m^1(X)−m^0(X)+Z(Y−m^1(X))pZ−(1−Z)(Y−m^0(X))1−pZ−(e^1(X)−e^0(X)+Z(W−e^1(X))pZ−(1−Z)(W−e^0(X))1−pZ)τ^(X))2]L(\hat{\tau}(X))=\hat{E}[(\hat{m}_1(X)-\hat{m}_0(X)+\frac{Z(Y-\hat{m}_1(X))}{p_Z}-\frac{(1 - Z)(Y-\hat{m}_0(X))}{1 - p_Z}-(\hat{e}_1(X)-\hat{e}_0(X)+\frac{Z(W-\hat{e}_1(X))}{p_Z}-\frac{(1 - Z)(W-\hat{e}_0(X))}{1 - p_Z})\hat{\tau}(X))^2]L(τ^(X))=E^[(m^1(X)−m^0(X)+pZZ(Y−m^1(X))−1−pZ(1−Z)(Y−m^0(X))−(e^1(X)−e^0(X)+pZZ(W−e^1(X))−1−pZ(1−Z)(W−e^0(X)))τ^(X))2] - 阶段3:与DR - 学习器类似,使用不同的分区排列再次重复阶段1和阶段2两次进行估计。最终的条件LATE模型是这3个条件LATE模型的平均值。
基于树的算法
提升树
提升树方法由一组使用基于树的算法的方法构成,其分裂准则基于提升差异。提出了三种不同方式来量化分裂导致的散度增益:
Dgain=Daftersplit(PT,PC)−Dbeforesplit(PT,PC)D_{gain}=D_{aftersplit}(P_T, P_C)-D_{beforesplit}(P_T, P_C)Dgain=Daftersplit(PT,PC)−Dbeforesplit(PT,PC)
其中DDD衡量散度,PTP_TPT和PCP_CPC分别指处理组和对照组中感兴趣结果的概率分布。该软件包中实现了三种量化散度的不同方式:KL散度、欧几里得距离和卡方。
KL散度
Kullback - Leibler(KL)散度由下式给出:
KL(P:Q)=∑k=left,rightpklogpkqkKL(P:Q)=\sum_{k = left, right}p_k\log\frac{p_k}{q_k}KL(P:Q)=k=left,right∑pklogqkpk
其中ppp是处理组中的样本均值,qqq是对照组中的样本均值,kkk表示计算ppp和qqq所在的叶子节点。
欧几里得距离
欧几里得距离由下式给出:
ED(P:Q)=∑k=left,right(pk−qk)2ED(P:Q)=\sum_{k = left, right}(p_k - q_k)^2ED(P:Q)=k=left,right∑(pk−qk)2
其中符号与上述相同。
卡方
最后,χ2\chi^2χ2散度由下式给出:
χ2(P:Q)=∑k=left,right(pk−qk)2qk\chi^2(P:Q)=\sum_{k = left, right}\frac{(p_k - q_k)^2}{q_k}χ2(P:Q)=∑k=left,rightqk(pk−qk)2
其中符号与上述相同。
DDP
delta-delta-p(ΔΔP)方法,其样本分裂准则定义如下:
ΔΔP=∣(PT(y∣a0)−PC(y∣a0)−(PT(y∣a1)−PC(y∣a1)))∣\Delta\Delta P = |(P_T(y|a_0)-P_C(y|a_0)-(P_T(y|a_1)-P_C(y|a_1)))|ΔΔP=∣(PT(y∣a0)−PC(y∣a0)−(PT(y∣a1)−PC(y∣a1)))∣
其中a0a_0a0和a1a_1a1是分裂A的结果,yyy是选定的类别,PTP_TPT和PCP_CPC分别是处理组和对照组的响应率。换句话说,我们首先计算每个分支中的响应率差异(ΔPleft\Delta P_{left}ΔPleft和ΔPright\Delta P_{right}ΔPright),随后计算它们之间的差异(ΔΔP=∣ΔPleft−ΔPright∣\Delta\Delta P = |\Delta P_{left}-\Delta P_{right}|ΔΔP=∣ΔPleft−ΔPright∣)。
IDDP
基于ΔΔP方法,实现了IDDP方法,其样本分裂准则定义如下:
IDDP=ΔΔP×I(ϕ,ϕl,ϕr)IDDP = \Delta\Delta P \times I(\phi, \phi_l, \phi_r)IDDP=ΔΔP×I(ϕ,ϕl,ϕr)
其中,ΔΔP∗\Delta\Delta P^*ΔΔP∗定义为ΔΔP−∣E[Y(1)−Y(0)]∣X∈ϕ∣\Delta\Delta P - |E[Y(1) - Y(0)]|_{X \in \phi}|ΔΔP−∣E[Y(1)−Y(0)]∣X∈ϕ∣,I(ϕ,ϕl,ϕr)I(\phi, \phi_l, \phi_r)I(ϕ,ϕl,ϕr)定义为:
I(ϕ,ϕl,ϕr)=H(nt(ϕ)n(ϕ),nc(ϕ)n(ϕ))×21+ΔΔP∗×3+nt(ϕ)n(ϕ)H(nt(ϕl)n(ϕ),nt(ϕr)n(ϕ))+nc(ϕ)n(ϕ)×H(nc(ϕl)n(ϕ),nc(ϕr)n(ϕ))+12I(\phi, \phi_l, \phi_r) = H(\frac{n_t(\phi)}{n(\phi)}, \frac{n_c(\phi)}{n(\phi)}) \times \frac{2}{1 + \Delta\Delta P^*} \times \frac{3 + \frac{n_t(\phi)}{n(\phi)}H(\frac{n_t(\phi_l)}{n(\phi)}, \frac{n_t(\phi_r)}{n(\phi)}) + \frac{n_c(\phi)}{n(\phi)} \times H(\frac{n_c(\phi_l)}{n(\phi)}, \frac{n_c(\phi_r)}{n(\phi)}) + 1}{2}I(ϕ,ϕl,ϕr)=H(n(ϕ)nt(ϕ),n(ϕ)nc(ϕ))×1+ΔΔP∗2×23+n(ϕ)nt(ϕ)H(n(ϕ)nt(ϕl),n(ϕ)nt(ϕr))+n(ϕ)nc(ϕ)×H(n(ϕ)nc(ϕl),n(ϕ)nc(ϕr))+1
其中,熵HHH定义为H(p,q)=(−p×log2(p))+(−q×log2(q))H(p, q) = (-p \times \log_2(p)) + (-q \times \log_2(q))H(p,q)=(−p×log2(p))+(−q×log2(q)),ϕ\phiϕ是与当前决策节点相关联的特征空间的一个子集,ϕl\phi_lϕl和ϕr\phi_rϕr分别是左子节点和右子节点。nt(ϕ)n_t(\phi)nt(ϕ)是处理样本的数量,nc(ϕ)n_c(\phi)nc(ϕ)是对照样本的数量,n(ϕ)n(\phi)n(ϕ)是当前(父)节点中所有样本的数量。
IT(Interaction Tree)
交互树(IT),其样本分裂准则是在所有可行的分裂中最大化GGG统计量:
G(s∗)=maxG(s)G(s^*) = \max G(s)G(s∗)=maxG(s)
其中G(s)=t2(s)G(s) = t^2(s)G(s)=t2(s),t(s)t(s)t(s)定义为:
t(s)=(y1L−y0L)−(y1R−y0R)σ×1n1+1n2+1n3+1n4t(s) = \frac{(y_{1L} - y_{0L}) - (y_{1R} - y_{0R})}{\sigma \times \sqrt{\frac{1}{n_1} + \frac{1}{n_2} + \frac{1}{n_3} + \frac{1}{n_4}}}t(s)=σ×n11+n21+n31+n41(y1L−y0L)−(y1R−y0R)
其中σ=∑i=14wisi2\sigma = \sum_{i = 1}^{4} w_i s_i^2σ=∑i=14wisi2是常数方差的合并估计量,wi=(ni−1)/∑j=14(nj−1)w_i = (n_i - 1) / \sum_{j = 1}^{4}(n_j - 1)wi=(ni−1)/∑j=14(nj−1)。此外,y1Ly_{1L}y1L、s12s_{1}^{2}s12和n1n_1n1分别是左子节点中处理组的样本均值、样本方差和样本大小。其他量也采用类似的表示方法。
请注意,此实现与原始实现的不同之处在于:(1) 剪枝技术;(2) 确定最佳树大小的验证方法。
CIT(Causal Inference Tree)
因果推断树(CIT),其样本分裂准则计算似然比检验统计量:
LRT(s)=−nτL2×ln(nτLSSEτL)−nτR2×ln(nτRSSEτR)+nτL1lnnτL1+nτL0lnnτL0+nτR1lnnτR1+nτR0lnnτR0LRT(s) = -\frac{n_{\tau L}}{2} \times \ln(\frac{n_{\tau L}}{SSE_{\tau L}}) - \frac{n_{\tau R}}{2} \times \ln(\frac{n_{\tau R}}{SSE_{\tau R}}) + n_{\tau L1} \ln n_{\tau L1} + n_{\tau L0} \ln n_{\tau L0} + n_{\tau R1} \ln n_{\tau R1} + n_{\tau R0} \ln n_{\tau R0}LRT(s)=−2nτL×ln(SSEτLnτL)−2nτR×ln(SSEτRnτR)+nτL1lnnτL1+nτL0lnnτL0+nτR1lnnτR1+nτR0lnnτR0
其中nτn_{\tau}nτ、nτ0n_{\tau 0}nτ0和nτ1n_{\tau 1}nτ1分别是节点τ\tauτ中的观测总数、分配到对照组的观测数以及分配到处理组的观测数。SSEτSSE_{\tau}SSEτ定义为:
SSEτ=∑i∈τ:ti=1(yi−y^t1)2+∑i∈τ:ti=0(yi−y^t0)2SSE_{\tau} = \sum_{i \in \tau : t_i = 1}(y_i - \hat{y}_{t1})^2 + \sum_{i \in \tau : t_i = 0}(y_i - \hat{y}_{t0})^2SSEτ=i∈τ:ti=1∑(yi−y^t1)2+i∈τ:ti=0∑(yi−y^t0)2
y^t0\hat{y}_{t0}y^t0和y^t1\hat{y}_{t1}y^t1分别是节点τ\tauτ中对照组和处理组的样本平均响应。
请注意,此实现与原始实现的不同之处在于:(1) 剪枝技术;(2) 确定最佳树大小的验证方法。
CTS(Contextual Treatment Selection)
上下文处理选择(CTS)方法,其样本分裂准则定义如下:
Δ^μ(s)=p^(ϕl∣ϕ)×maxt=0,…,Ky^t(ϕl)+p^(ϕr∣ϕ)×maxt=0,…,Ky^t(ϕr)−maxt=0,…,Ky^t(ϕ)\hat{\Delta}\mu(s) = \hat{p}(\phi_l | \phi) \times \max_{t = 0, \ldots, K} \hat{y}^t(\phi_l) + \hat{p}(\phi_r | \phi) \times \max_{t = 0, \ldots, K} \hat{y}^t(\phi_r) - \max_{t = 0, \ldots, K} \hat{y}^t(\phi)Δ^μ(s)=p^(ϕl∣ϕ)×t=0,…,Kmaxy^t(ϕl)+p^(ϕr∣ϕ)×t=0,…,Kmaxy^t(ϕr)−t=0,…,Kmaxy^t(ϕ)
其中ϕl\phi_lϕl和ϕr\phi_rϕr分别指左叶和右叶中的特征子空间,p^(ϕj∣ϕ)\hat{p}(\phi_j | \phi)p^(ϕj∣ϕ)表示在给定ϕ\phiϕ的情况下,主体处于ϕj\phi_jϕj的估计条件概率,y^t(ϕj)\hat{y}^t(\phi_j)y^t(ϕj)是在处理ttt下的条件期望响应。
价值优化方法
当处理成本较高时,该软件包支持用于分配处理组的方法。为了理解这个问题,将总体分为以下四类会有所帮助:
- 依从者:这类人只有在接受处理时才会有有利的结果。
- 总是接受者:这类人无论是否接受处理都会有有利的结果。
- 从不接受者:这类人无论是否接受处理都不会有有利的结果。
- 违抗者:这类人只有在不接受处理时才会有有利的结果。
反事实单元选择
提出了一种使用反事实逻辑选择处理单元的方法。假设选择上述不同类别单元有以下收益:
- 依从者:β\betaβ
- 总是接受者:γ\gammaγ
- 从不接受者:θ\thetaθ
- 违抗者:δ\deltaδ
如果XXX表示个体的特征集,单元选择问题可以表述如下:
argmaxXβP(complier∣X)+γP(always - taker∣X)+θP(never - taker∣X)+δP(defier∣X)\underset{X}{\text{argmax}} \beta P(\text{complier} | X) + \gamma P(\text{always - taker} | X) + \theta P(\text{never - taker} | X) + \delta P(\text{defier} | X)XargmaxβP(complier∣X)+γP(always - taker∣X)+θP(never - taker∣X)+δP(defier∣X)
这个问题可以用反事实逻辑重新表述。假设W=wW = wW=w表示个体接受处理,W=w′W = w'W=w′表示个体未接受处理。类似地,令F=fF = fF=f表示个体有有利结果,F=f′F = f'F=f′表示个体有不利结果。那么优化问题变为:
argmaxXβP(fw,fw′′∣X)+γP(fw,fw′∣X)+θP(fw′,fw′′∣X)+δP(fw′,fw′∣X)\underset{X}{\text{argmax}} \beta P(f_w, f_{w'}' | X) + \gamma P(f_w, f_{w'} | X) + \theta P(f_{w'}, f_{w'}' | X) + \delta P(f_{w'}, f_{w'} | X)XargmaxβP(fw,fw′′∣X)+γP(fw,fw′∣X)+θP(fw′,fw′′∣X)+δP(fw′,fw′∣X)
请注意,上述内容仅仅是根据相关用户细分的定义得出的。随后使用反事实逻辑在某些条件下解决上述优化问题。
注意:该软件包中当前的实现具有很强的实验性。
反事实价值估计器
该软件包中实现的反事实价值估计方法使用标准机器学习模型预测一个单元在不同处理条件下的结果。将一个单元分配到特定处理的预期值由下式给出:
E[(v−ccw)Yw−icw]E[(v - cc_w)Y_w - ic_w]E[(v−ccw)Yw−icw]
其中YwY_wYw是在给定处理www下有利事件(如转化)发生的概率,vvv是有利事件的价值,ccwcc_wccw是有利事件发生时触发的处理成本,icwic_wicw是无论结果是否有利与处理相关的成本。
因果概率
如果一个结果在没有某个原因的情况下不会发生,那么这个原因对于该结果来说被称为是“必要的”。如果一个结果在有某个原因的情况下会发生,那么这个原因对于该结果来说被称为是“充分的”。如果上述两个条件都成立,那么这个原因对于该结果来说被称为是“必要且充分的”。Jin Tian and Judea Pearl表明,我们可以计算出一个原因属于上述三种类型中每一种的概率界限。
为了理解如何计算因果概率的界限,我们需要特殊的符号来表示反事实量。令yty_tyt表示命题“如果处理组设置为‘处理’,yyy会发生”,yc′y_{c}'yc′表示命题“如果处理组设置为‘对照’,yyy不会发生”,对于(假设为)二元结果和处理变量的其余两种组合也类似表示。
那么,处理对于yyy的发生是“充分的”概率可以定义为:
PS=P(yt∣c,y′)P_S = P(y_t | c, y')PS=P(yt∣c,y′)
这是指在实际上处理设置为对照且结果未发生的情况下,如果处理设置为ttt,yyy会发生的概率。
处理对于yyy的发生是“必要的”概率可以定义为:
PN=P(yc′∣t,y)P_N = P(y_{c}' | t, y)PN=P(yc′∣t,y)
这是指在实际上yyy发生且处理发生的情况下,如果处理设置为对照,yyy不会发生的概率。
最后,处理既是必要又是充分的概率定义为:
PNS=P(yt,yc′)P_{NS} = P(y_t, y_{c}')PNS=P(yt,yc′)
它表示如果处理发生,yyy会发生;如果处理不发生,yyy不会发生。PNSP_{NS}PNS与PNP_NPN和PSP_SPS的关系如下:
PNS=P(t,y)PN+P(c,y′)PSP_{NS} = P(t, y)P_N + P(c, y')P_SPNS=P(t,y)PN+P(c,y′)PS
在界定上述三个量时,除了实验数据,我们还利用观测数据。观测数据通过联合概率来表征:
PTY=P(t,y),P(c,y),P(t,y′),P(c,y′)P_{TY} = P(t, y), P(c, y), P(t, y'), P(c, y')PTY=P(t,y),P(c,y),P(t,y′),P(c,y′)
这个程序的主要思想是将界定任务转化为一个线性规划问题(关于他们方法的现代实现
使用线性规划方法,并结合某些约束条件和观测数据,Jin Tian and Judea Pearl发现PNSP_{NS}PNS的精确下限为:
max{0,P(yt)−P(yc),P(y)−P(yc),P(yt)−P(y)}\max\{0, P(y_t) - P(y_c), P(y) - P(y_c), P(y_t) - P(y)\}max{0,P(yt)−P(yc),P(y)−P(yc),P(yt)−P(y)}
精确上限为:
min{P(yt),P(yc′),P(t,y)+P(c,y′),P(yt)−P(yc)+P(t,y′)+P(c,y)}\min\{P(y_t), P(y_{c}'), P(t, y) + P(c, y'), P(y_t) - P(y_c) + P(t, y') + P(c, y)\}min{P(yt),P(yc′),P(t,y)+P(c,y′),P(yt)−P(yc)+P(t,y′)+P(c,y)}
他们使用类似的程序来找到PSP_SPS和PNP_NPN的界限。get_pns_bounds()
函数使用Jin Tian and Judea Pearl中的结果计算这三个因果概率中每一个的界限。
选定的传统方法
该软件包支持选定的传统因果推断方法。这些方法通常用于使用观测(非实验)数据进行因果推断。在这类研究中,处理组和对照组之间观察到的差异通常不等于“潜在结果”E[Y(1)−Y(0)]E[Y(1) - Y(0)]E[Y(1)−Y(0)]之间的差异。因此,以下方法试图以不同方式处理这个问题。
匹配法
匹配法的总体思路是找到在相关特征方面尽可能相似的处理单元和未处理单元。因此,匹配方法可以被视为试图模仿随机对照试验的因果推断方法家族的一部分。
虽然有许多不同的方法来匹配处理单元和未处理单元,但最常见的方法是使用倾向得分:
ei(Xi)=P(Wi=1∣Xi)e_i(X_i) = P(W_i = 1 | X_i)ei(Xi)=P(Wi=1∣Xi)
然后根据e(X)e(X)e(X)使用某种距离准则(如k:1k:1k:1最近邻)来匹配处理单元和未处理单元。由于匹配通常在处理群体和对照群体之间进行,这种方法估计处理组的平均处理效应(ATT):
E[Y(1)∣W=1]−E[Y(0)∣W=1]E[Y(1) | W = 1] - E[Y(0) | W = 1]E[Y(1)∣W=1]−E[Y(0)∣W=1]
处理逆概率加权法
处理逆概率加权(IPTW)方法使用倾向得分eee,通过实际处理WWW的概率的倒数来对处理群体和未处理群体进行加权。对于二元处理W∈{1,0}W \in \{1, 0\}W∈{1,0}:
We+1−W1−e\frac{W}{e} + \frac{1 - W}{1 - e}eW+1−e1−W
通过这种方式,IPTW方法可以被视为创建了一个人工总体,其中处理单元和未处理单元在其观测特征XXX方面是相似的。
与匹配法相比,IPTW的一个可能优点是,由于处理单元和未处理单元之间缺乏重叠而丢弃的数据可能更少。该方法的一个已知问题是,极端的倾向得分可能会产生高度可变的估计量。已经提出了不同的方法来修剪和标准化IP权重。
两阶段最小二乘法(2SLS)
识别WWW对YYY的处理效应的基本要求之一是,在给定协变量XXX的条件下,WWW与YYY的潜在结果正交。如果WWW和YYY都受到一个未观测到的变量(从YYY中去除WWW的真实效应后的误差项)的影响,而这个变量不在XXX中,那么这个条件可能会被违反。在这种情况下,工具变量方法试图借助与WWW相关但与误差项不相关的第三个变量ZZZ来估计WWW对YYY的效应。换句话说,工具变量ZZZ仅通过经过WWW的有向路径与YYY相关。如果这些条件得到满足,在没有协变量的情况下,WWW对YYY的效应可以使用以下样本类似公式进行估计:
Cov(Yi,Zi)Cov(Wi,Zi)\frac{Cov(Y_i, Z_i)}{Cov(W_i, Z_i)}Cov(Wi,Zi)Cov(Yi,Zi)
工具变量估计最常用的方法是两阶段最小二乘法(2SLS)。在这种方法中,原因变量WWW首先对工具变量ZZZ进行回归。然后,在第二阶段,将感兴趣的结果YYY对第一阶段模型的预测值进行回归。直观地说,通过仅使用由ZZZ的变化引起的WWW的变化比例来估计WWW对YYY的效应。具体来说,假设我们有线性模型:
Y=Wα+Xβ+u=Ξγ+uY = W\alpha + X\beta + u = \Xi\gamma + uY=Wα+Xβ+u=Ξγ+u
为方便起见,这里我们令Ξ=[W,X]\Xi = [W, X]Ξ=[W,X]且γ=[α′,β′]′\gamma = [\alpha', \beta']'γ=[α′,β′]′,假设我们有工具变量 ZZZ,其列数至少与 WWW 的列数相同,令 Ω=[Z,X]\Omega = [Z, X]Ω=[Z,X],2SLS 估计量如下:
γ^2SLS=[Ξ′Ω(Ω′Ω)−1Ω′Ξ]−1[Ξ′Ω′(Ω′Ω)−1Ω′Y]\hat{\gamma}_{2SLS} = [\Xi'\Omega(\Omega'\Omega)^{-1}\Omega'\Xi]^{-1}[\Xi'\Omega'(\Omega'\Omega)^{-1}\Omega'Y]γ^2SLS=[Ξ′Ω(Ω′Ω)−1Ω′Ξ]−1[Ξ′Ω′(Ω′Ω)−1Ω′Y]
局部平均处理效应(LATE)
在许多情况下,处理措施 WWW 可能取决于受试者自己的选择,并且无法在实验环境中直接实施。然而,可以将用户随机分配到处理组/对照组,以便促使处理组中的用户接受处理。这就是不依从的情况,即用户可能不遵守他们被分配的状态 ZZZ,决定是否接受处理。与价值优化方法部分类似,通常在这种情况下有 3 种类型的用户:
- 依从者:这些人只有在被分配到处理组时才会接受处理。
- 总是接受者:这些人无论被分配到哪个组都会接受处理。
- 从不接受者:这些人无论被分配到哪个组都不会接受处理。
然而,为了便于识别,我们假设不存在违抗者,即那些只有在被分配到对照组时才会接受处理的人。
在这种情况下,可以测量依从者的处理效应:
τ^Complier=E[Y∣Z=1]−E[Y∣Z=0]E[W∣Z=1]−E[W∣Z=0]\hat{\tau}_{Complier} = \frac{E[Y|Z = 1] - E[Y|Z = 0]}{E[W|Z = 1] - E[W|Z = 0]}τ^Complier=E[W∣Z=1]−E[W∣Z=0]E[Y∣Z=1]−E[Y∣Z=0]
这就是局部平均处理效应(LATE)。如果我们将分配状态 ZZZ 作为工具变量,该估计量也等同于 2SLS。
针对平均处理效应(ATE)的靶向最大似然估计(TMLE)
靶向最大似然估计(TMLE)提供了一种双稳健半参数方法,借助机器学习算法直接“针对”平均处理效应。与其他方法(包括结果回归和处理逆概率加权)相比,TMLE 通常表现更好,特别是在处理有偏处理和异常值时。
给定二元处理 WWW、协变量 XXX 和结果 YYY,针对 ATE 的 TMLE 按以下步骤进行:
-
步骤 1
使用交叉拟合通过机器学习估计倾向得分 e^(x)\hat{e}(x)e^(x)、处理组的预测结果 m^1(x)\hat{m}_1(x)m^1(x) 以及对照组的预测结果 m^0(x)\hat{m}_0(x)m^0(x)。 -
步骤 2
将 YYY 缩放为 Y~=Y−minYmaxY−minY\widetilde{Y} = \frac{Y - \min Y}{\max Y - \min Y}Y=maxY−minYY−minY,使得 Y~∈[0,1]\widetilde{Y} \in [0, 1]Y∈[0,1]。使用相同的缩放函数将 m^i(x)\hat{m}_i(x)m^i(x) 转换为 m~i(x)\widetilde{m}_i(x)mi(x),i=0,1i = 0, 1i=0,1。裁剪缩放后的函数,使其值保持在单位区间内。 -
步骤 3
令 Q=log(m~W(X)1−m~W(X))Q = \log(\frac{\widetilde{m}_W(X)}{1 - \widetilde{m}_W(X)})Q=log(1−mW(X)mW(X))。最大化以下伪对数似然函数:
maxh0,h1−1N∑i[Y~ilog(1+exp(−Qi−h0(1−W)1−e^(Xi)−h1We^(Xi)))+(1−Y~i)log(1+exp(Qi+h0(1−W)1−e^(Xi)+h1We^(Xi)))]\max_{h_0, h_1} -\frac{1}{N} \sum_{i} [\widetilde{Y}_i \log(1 + \exp(-Q_i - \frac{h_0(1 - W)}{1 - \hat{e}(X_i)} - \frac{h_1W}{\hat{e}(X_i)})) + (1 - \widetilde{Y}_i) \log(1 + \exp(Q_i + \frac{h_0(1 - W)}{1 - \hat{e}(X_i)} + \frac{h_1W}{\hat{e}(X_i)}))]h0,h1max−N1i∑[Yilog(1+exp(−Qi−1−e^(Xi)h0(1−W)−e^(Xi)h1W))+(1−Yi)log(1+exp(Qi+1−e^(Xi)h0(1−W)+e^(Xi)h1W))]
- 步骤 4
令:
Q~0=11+exp(−Q−h01−e^(X))\widetilde{Q}_0 = \frac{1}{1 + \exp(-Q - \frac{h_0}{1 - \hat{e}(X)})}Q0=1+exp(−Q−1−e^(X)h0)1
Q~1=11+exp(−Q−h1e^(X))\widetilde{Q}_1 = \frac{1}{1 + \exp(-Q - \frac{h_1}{\hat{e}(X)})}Q1=1+exp(−Q−e^(X)h1)1
将 Q~1\widetilde{Q}_1Q1 和 Q~0\widetilde{Q}_0Q0 的差异重新缩放回原始范围后,取样本平均值作为 ATE 的估计值。