随机森林环境因子重要性排序分析 随机森林 因果推断

Paper Review

• 1. Estimation and Inference of Heterogeneous Treatment Effects using Random Forest

• 1.1 Asymptotic analysis
• 1.2 Double-Sample Trees

• 2. Generalized Random Forests

• 2.1 Algorithm

• 1. Forest-based local estimation
• 2. Splitting to maximize heterogeneity
• 3. The gradient tree algorithm

• 2.2 Asymptotic analysis
• 2.3 Experiments

• 1. CAPE
• 2. Quantile Regression Forest

• 3. Orthogonal Random Forest for Causal Inference

• 3.1 Introduction
• 3.2 Algorithm

• 1. first stage
• 2. second stage

• 3.3 Experiments

• 4. Decision trees for uplift modeling with single and multiple treatments

• 4.1 Single Treatment
• 4.2 Multiple treatment

因果森林总结：基于树模型的异质因果效应估计

Uplift model with multiple treatments

1. Estimation and Inference of Heterogeneous Treatment Effects using Random Forest

二元干预情形下估计$\tau(x)=E[Y^1-Y^0|X=x]$

1.1 Asymptotic analysis

• Under some condition,
  $(\hat{\tau}(x)-\tau(x)) / \sqrt{\operatorname{Var}[\hat{\tau}(x)]} \Rightarrow \mathcal{N}(0,1)$
• $\operatorname{Var}[\hat{\tau}(x)]$可以用infinitesimal jackknife估计$\widehat{V}_{I J}(x) / \operatorname{Var}[\hat{\tau}(x)] \rightarrow 1$
  $\widehat{V}_{I J}(x)=\frac{n-1}{n}\left(\frac{n}{n-s}\right)^2 \sum_{i=1}^n \operatorname{Cov}_*\left[\hat{\tau}_b^*(x), N_{i b}^*\right]^2$
  其中，系数项$(n-1)n/(n-s)^2$只能对无放回的子抽样做修正

证明过程分为两步：

• 先证明偏差$E[\hat \mu_n(x)-\mu(x)]$的bound

• 再证明$\hat \mu_n(x)-E[\hat \mu_n(x)]$近似正态

利用Hajek projection和k-PNN先证明T is ν-incremental
$\stackrel{\circ}{T}=\mathbb{E}[T]+\sum_{i=1}^n\left(\mathbb{E}\left[T \mid Z_i\right]-\mathbb{E}[T]\right)$

1.2 Double-Sample Trees

回归树T分裂准则为最小化MSE， $\hat{\mu}(x)=\frac{1}{\left|\left\{i: X_i \in L(x)\right\}\right|} \sum_{\left\{i: X_i \in L(x)\right\}} Y_i=\bar Y_L$

$\sum_{i \in \mathcal{J}}\left(\hat{\mu}\left(X_i\right)-Y_i\right)^2=\sum_{i \in \mathcal{J}} Y_i^2-\sum_{i \in \mathcal{J}} \hat{\mu}\left(X_i\right)^2$

考虑到$\sum_{i \in \mathcal{J}} \hat{\mu}\left(X_i\right)=\sum_{i \in \mathcal{J}} Y_i$，上式等价于最大化$\hat{\mu}(X_i)$的方差

2. Generalized Random Forests

2.1 Algorithm

1. Forest-based local estimation

目的：给定$(O_i, X_i)$， 估计$\theta(\cdot)$，如估计HTE时，$O_i=(Y_i, W_i)$。
方法：求解方程$\mathbb{E}\left[\psi_{\theta(x), v(x)}\left(O_i\right) \mid X_i=x\right]=0$，其中，$\theta(x), v(x)$分别是感兴趣的参数和无关参数

• 权重估计阶段：$\alpha_i(x)$衡量$x_i$和$x$的相似程度，将同一叶子结点中的"“共现频率”"作为其权重
  $\alpha_{b i}(x)=\frac{\mathbf{1}\left(\left\{X_i \in L_b(x)\right\}\right)}{\left|L_b(x)\right|}, \quad \alpha_i(x)=\frac{1}{B} \sum_{b=1}^B \alpha_{b i}(x)$
  其中$L_b(x)$为第b棵树$x$所在叶子结点的所有数据
• 加权求解
  $(\hat{\theta}(x), \hat{v}(x)) \in \underset{\theta, v}{\operatorname{argmin}}\left\{\left\|\sum_{i=1}^n \alpha_i(x) \psi_{\theta, v}\left(O_i\right)\right\|_2\right\}$
  例子：求解$\mu(x)=\mathbb{E}\left[Y_i \mid X_i=x\right]=0$,取 $\psi_{u(x)}\left(Y_i\right)=Y_i-\mu(x)$，则有$\sum_{i=1}^n \frac{1}{B} \sum \alpha_{b i}(x)\left(Y_i-\hat{\mu}(x)\right)=0$，方程的解为$\hat{\mu}(x)=\frac{1}{B} \sum \hat{\mu}_b(x)$

2. Splitting to maximize heterogeneity

针对某一节点P和数据J，参数的估计方法为
$\left(\hat{\theta}_P, \hat{v}_P\right)(J) \in \arg \min _{\theta, v}\left\|\sum_{\left\{i \in J: X_i \in P\right\}} \psi_{\theta, v}\left(O_i\right)\right\|_2$
将结点P分裂为两个子节点$C_1, C_2$，目标为最小化$err(C_1, C_2)$
$\operatorname{err}\left(C_1, C_2\right)=\sum_{j=1,2} P\left(X \in C_j \mid X \in P\right) E\left[\left(\hat{\theta}_{C_j}(J)-\theta(X)\right)^2 \mid X \in C_j\right]$
在某些条件下， $\operatorname{err}\left(C_1, C_2\right)=K(P)-\mathbb{E}\left[\Delta\left(C_1, C_2\right)\right]+o\left(r^2\right)$，所以分裂等价于最大化节点间的异质性，即
$\triangle\left(C_1, C_2\right)=\frac{n_{C_1} n_{C_2}}{n_{C_P}^2}\left(\hat{\theta}_{C_1}(J)-\hat{\theta}_{C_2}(J)\right)^2$

3. The gradient tree algorithm

为减少计算量，采用梯度近似
$\tilde{\theta}_C=\hat{\theta}_P-\frac{1}{\left|\left\{i: X_i \in C\right\}\right|} \sum_{\left\{i: X_i \in C\right\}} \xi^{\top} A_P^{-1} \psi_{\hat{\theta}_P, \hat{v}_P}\left(O_i\right)$
其中，$\xi$取出$\theta$的值，消去无关参数，$A_P$近似$\nabla \mathbb{E}\left[\psi_{\hat{\theta}_P, \hat{v}_P}\left(O_i\right) \mid X_i \in P\right]$
$A_P=\frac{1}{\left|\left\{i: X_i \in P\right\}\right|} \sum_{\left\{i: X_i \in P\right\}} \nabla \psi_{\hat{\theta}_P, \hat{v}_P}\left(O_i\right)$
注：当$\psi$不可导时，可以采用分位数回归

故分裂阶段可以分为以下2步

• labeling step：计算父节点的 $\hat{\theta}_P, \hat{v}_P, A_P^{-1}$ ，以及每个样本的伪值
  $\rho_i=-\xi^T A_P^{-1} \psi_{\hat{\theta}_P, \hat{v}_P}\left(O_i\right)$
• regression step：最大化近似分裂准则
  $\tilde{\triangle}\left(C_1, C_2\right)=\sum_{j=1}^2 \frac{1}{\left|\left\{i: X_i \in C_j\right\}\right|}\left(\sum_{\left\{i: X_i \in C_j\right\}} \rho_i\right)^2$

回归过程中，分裂准则的近似误差在一定范围内

2.2 Asymptotic analysis

定义expected score function
$M_{\theta, \nu}(x):=\mathbb{E}\left[\psi_{\theta, v}(O) \mid X=x\right]$

• consistency
• approximate normality

2.3 Experiments

1. CAPE

$Y_i=W_i \cdot b_i+\varepsilon_i, \beta(x)=\mathbb{E}\left[b_i \mid X_i=x\right]$目标是估计$\theta(x)=\xi \cdot \beta(x)，\xi \in \mathbb{R}^p$，score function
$\psi_{\beta(x), c(x)}\left(Y_i, W_i\right)=\left(Y_i-\beta(x) \cdot W_i-c(x)\right)\left(1 \ W_i^{\top}\right)^{\top}$此时 $\arg \min _\theta\left\|E\left[\psi_{\theta(x), v(x)}\left(O_i\right)\right]\right\|_2$ 相当于:
$\theta(x)=\xi^{\top} \operatorname{Var}\left[W_i \mid X_i=x\right]^{-1} \operatorname{Cov}\left[W_i, Y_i \mid X_i=x\right\rceil$

• Forest
  $\hat{\theta}(x)=\xi^{\top}\left(\sum_{i=1}^n \alpha_i(x)\left(W_i-\bar{W}_\alpha\right)^{\otimes 2}\right)^{-1} \sum_{i=1}^n \alpha_i(x)\left(W_i-\bar{W}_\alpha\right)\left(Y_i-\bar{Y}_\alpha\right)$其中，$\bar{W}_\alpha=\sum \alpha_i(x) W_i，\bar{Y}_\alpha=\sum \alpha_i(x) Y_i， v^{\otimes 2}=v v^{\top}$
  GRF算法实施时，权重可自动求解，但需要计算对应的伪结果，注意$\rho_i和A_P$只关注$\theta(x)$
  $\psi_i=W_i \left(Y_i-\bar{Y}_P-\left(W_i-\bar{W}_P\right) \hat{\beta}_P\right)$$\nabla \psi_i = W_i^{\otimes 2}$
  对比实际代入的表达式，实际对W做centering
  $\rho_i =\xi^{\top} A_P^{-1}\left(W_i-\bar{W}_P\right)\left(Y_i-\bar{Y}_P-\left(W_i-\bar{W}_P\right) \hat{\beta}_P\right)$ $A_P =\frac{1}{\left|\left\{i: X_i \in P\right\}\right|} \sum_{\left\{i: X_i \in P\right\}}\left(W-\bar{W}_P\right)^{\otimes 2}$
• Local Centering
  提前对Y和W做中心化处理，类似残差，使得估计效果更好
  $\begin{aligned} \theta(x)= & \xi^{\top} \operatorname{Var}\left[\left(W_i-\mathbb{E}\left[W_i \mid X_i\right]\right) \mid X_i \in \mathcal{S}\right]^{-1} \\ & \times \operatorname{Cov}\left[\left(W_i-\mathbb{E}\left[W_i \mid X_i\right]\right),\left(Y_i-\mathbb{E}\left[Y_i \mid X_i\right]\right) \mid X_i \in \mathcal{S}\right] \end{aligned}$

2. Quantile Regression Forest

$\psi_\theta\left(Y_i\right)=q \mathbf{1}\left(\left\{Y_i>\theta\right\}\right)-(1-q) \mathbf{1}\left(\left\{Y_i \leq \theta\right\}\right)$

3. Orthogonal Random Forest for Causal Inference

3.1 Introduction

DML的优势：即使第一阶段的估计有误差，第二阶段的估计仍可以近似正态；劣势：HTE预设参数形式。CF的优势：非参数估计；劣势：很大程度上要求低维度W。ORF在GRF的基础上，参考DML新增对无关参数的正交估计（First stage），减少误差。

At a high level, ORF can be viewed as an orthogonalized version of GRF that is more robust to the nuisance estimation error. The key modification to GRF’s tree learner is our incorporation of orthogonal nuisance estimation in the splitting criterion.

3.2 Algorithm

建树时每次分裂的过程two-stage

1. first stage

$\hat \nu(x)=\arg \min \sum a_i l (Z_i;\nu)+\lambda ||\nu||_1$

2. second stage

• split

具体执行分裂的算法类似GRF的gradient tree algorithm，但考虑到honesty，集合略有改动，其中$S^1$是用于分裂的数据，$\hat h_P$是first stage估计的无关参数
$\tilde{\theta}_C=\hat{\theta}_P-\frac{1}{\left|C \cap S^1\right|} \sum_{i \in C_j \cap S^1} A_P^{-1} \psi\left(Z_i ; \hat{\theta}_P, \hat{h}_P\left(X_i, W_i\right)\right)$
where $A_P=\frac{1}{\left|P \cap S^1\right|} \sum_{i \in P \cap S_b^1} \nabla_\theta \psi\left(Z_i ; \hat{\theta}_P, \hat{h}_P\left(X_i, W_i\right)\right)$

• labeling step：计算父节点的 $\hat{\theta}_P, \hat{h}_P, A_P^{-1}$ ，以及每个样本的伪值
  $\rho_{t, i}=A_P^{-1} \psi_t\left(Z_i ; \hat{\theta}_P, \hat{h}_P\left(X_i, W_i\right)\right)$
• regression step：maximize proxy heterogeneity score
  $\tilde{\Delta}_t\left(C_1, C_2\right)=\sum_{j=1}^2 \frac{1}{\left|C_j \cap S^1\right|}\left(\sum_{i \in C_j \cap S^1} \rho_{t, i}\right)$
• Predict

$a_{ib}$同样仅限于$S^2$的估计样本
$a_{i b}=\frac{\mathbf{1}\left[\left(X_i \in L_b(x)\right) \wedge\left(Z_i \in S_b^2\right)\right]}{\left|L_b(x) \cap S_b^2\right|}, \quad a_i=\frac{1}{B} \sum_{b=1}^B a_{i b}$

以下定理保证了$a_{ib}$在x邻域内非零

3.3 Experiments

• DML Partially Linear Regression(PLR, Robinson, 1988)
  $\begin{array}{cl} Y=D \theta_0+g_0(X)+U, & \mathrm{E}[U \mid X, D]=0 \\ D=m_0(X)+V, & \mathrm{E}[V \mid X]=0 \end{array}$则score function为$\psi(W ; \theta, \eta)=(Y-D \alpha-g(X))(D-m(X))$
• ORF
  数据$D = \{Z_i = (T_i, Y_i, W_i, X_i)\}_{i=1}^{2n}$，其中T是连续或离散的Treatment，Y是outcome，$W \in [-1,1]^{d_\nu}$是potential confounders/controls，$X \in [0,1]^d$是特征
  $\begin{array}{cl} Y=\left\langle\mu_0(X, W), T\right\rangle+f_0(X, W)+\varepsilon, & \mathbb{E}[\varepsilon \mid W, X, T]=0 \\ T=g_0(X, W)+\eta, & \mathbb{E}[\eta \mid X, W, \varepsilon]=0 \end{array}$confounders分别通过$f_0$和$g_0$影响outcome和treatment
  $\mu_0: \mathbb{R}^d \times \mathbb{R}^{d_\nu} \rightarrow[-1,1]^p$为treatment effect function，目标是估计CATE
  $\theta_0(x)=\mathbb{E}\left[\mu_0(X, W) \mid X=x\right]$
  基于DML思想，残差化
  $Y-\mathbb{E}[Y \mid X, W]=\left\langle\mu_0(X, W), T-\mathbb{E}[T \mid X, W]\right\rangle+\varepsilon$定义$q_0(X, W)=\mathbb{E}[Y \mid X, W]$，$\tilde{Y}=Y-q(X, W)$，$\tilde{T}=T-$ $g_0(X, W)=\eta$， 则有
  $\mathbb{E}[\tilde{Y} \mid X, \tilde{T}]=\mathbb{E}\left[\mu_0(X, W) \mid X\right] \cdot \tilde{T}=\theta(X) \cdot \tilde{T}$则score function$\psi(Z ; \theta, h(X, W)) = \{Y-q(X, W)-\theta(T-g(X, W)\rangle)\}(T-g(X, W))$
  其中$q, g$是$q_0, g_0$的估计

4. Decision trees for uplift modeling with single and multiple treatments

4.1 Single Treatment

• Split rule：maximize the differences between class distributions
  $D_{\text {gain }}(A)=D\left(P^T(Y): P^C(Y) \mid A\right)-D\left(P^T(Y): P^C(Y)\right)$
• Normalising：C4.5对gain除以info避免bias，而本文的norm主要惩罚两边子节点中treatment和control组比例不平衡的，这和随机试验的假设相悖
  下式第一项系数考虑比例不平衡，后两项考虑相对样本大小
  (1) D=KL:
  $\begin{aligned} I(A)= & H\left(\frac{N^T}{N}, \frac{N^C}{N}\right) K L\left(P^T(A): P^C(A)\right) \\ & +\frac{N^T}{N} H\left(P^T(A)\right)+\frac{N^C}{N} H\left(P^C(A)\right)+\frac{1}{2}\end{aligned}$
  (2) D=欧式/卡方
  $\begin{aligned} & J(A)=\operatorname{Gini}\left(\frac{N^T}{N}, \frac{N^C}{N}\right) D\left(P^T(A): P^C(A)\right) \\ &+ \frac{N^T}{N} \operatorname{Gini}\left(P^T(A)\right)+\frac{N^C}{N} \operatorname{Gini}\left(P^C(A)\right)+\frac{1}{2}\end{aligned}$

4.2 Multiple treatment

• Split rule

• Normalizing
  $\begin{aligned} I(A)=\alpha H & \left(\frac{N^T}{N}, \frac{N^C}{N}\right) K L\left(P^T(A): P^C(A)\right) \\ & +(1-\alpha) \sum_{i=1}^k H\left(\frac{N^{T_i}}{N^{T_i}+N^C}, \frac{N^C}{N^{T_i}+N^C}\right) K L\left(P^{T_i}(A): P^C(A)\right) \\ & +\sum_{i=1}^k \frac{N^{T_i}}{N} H\left(P^{T_i}(A)\right)+\frac{N^C}{N} H\left(P^C(A)\right)+\frac{1}{2} \end{aligned}$