论文笔记：Generalized Random Forests

本文暂且忽略渐近线证明部分，关注于GRF的prediction和split方法
ref：

• 知乎博客：https://zhuanlan.zhihu.com/p/448524822
• S. Athey, J. Tibshirani, and S. Wager, “Generalized random forests,” Ann. Statist., vol. 47, no. 2, Apr. 2019, doi: 10.1214/18-AOS1709.

Motivation

本文旨在找到一种general的forest-based的估计方法，是对random forest的泛化扩展。这也是该工作的最大贡献。具体而言，该工作所提出的General Object是：
$$\mathbb{E}[{\Psi }_{\theta (x),\nu (x)}(O_i)|X_i=x]=0\text{\hspace{1em}(1)}$$
其中，$\Psi (\cdot )$是scoring function，可以理解为loss function或者是优化目标，$\theta (x)$是我们期望估计的量，$\nu (x)$是可选的nuisance parameter，$O_i$是和$\theta (x)$有关的量。目的就是希望能够造一个森林使得Eq(1)成立。

为了实现上述目标（Eq(1)），我们需要解决如下优化问题：
$$(\hat{\theta }(x),\hat{\nu }(x))=\arg \underset{\theta ,\nu }{\min }\Vert \sum _{i=1}^n{\alpha }_i(x)\cdot {\Psi }_{\theta ,\nu }(O_i){\Vert }_2\text{\hspace{1em}(2)}$$
而该优化问题的最优解$\hat{\theta }(x)$就是我们的估计结果；

至于${\alpha }_i(x)$，它表示训练样本$i$与测试样本（由$x$表示）的相似度，起到加权的作用，其具体计算方法如下：
$${\alpha }_i(x)=\frac{1}{B}\cdot \sum _{b=1}^B{\alpha }_{b_i}(x)\text{\hspace{1em}(3)}$$
α b i ( x ) = 1 ( { X i ∈ L b ( x ) } ) ∣ L b ( x ) ∣ (4) \alpha_{b_i}(x)=\frac{1(\{X_i\in L_b(x)\})}{|L_b(x)|} \tag{4} αbi​​(x)=∣Lb​(x)∣1({Xi​∈Lb​(x)})​(4)
其中，$B$表示树的个数，$b$是代表第$b$棵树，$L_b(x)$表示与测试样本$x$落在第$b$棵树同一叶子上的训练样本集合。所以，${\alpha }_{b_i}(x)$就表示在第$b$棵树中第$i$个训练样本与测试样本$x$落在同一个叶子节点的频率【该频率反映相似度】。注意，${\sum }_{i=1}^n{\alpha }_i(x)=1$！

总结一下，Eq(1)和Eq(2)其实是等价的，在森林构建好之后可以根据Eq(1)或(2)进行预测评估，这俩式子就是GRF的核心，许多统计问题（如，最小二乘、最大似然、分位数回归等）都可以看作是Eq(1)的特例。

Case of Regression
以回归问题为例，证明random forest是GRF的一个特例：
对于回归问题而言，我们关心的估计量$\mu (x)=\mathbb{E}[Y_i|X_i=x]$（这里$\mu (x)$就是$\theta (x)$)，对应的scoring function就是${\Psi }_{\mu (x)}(O_i)=Y_i-\mu (x)$。
同时，我们知道random forest在森林构建好之后，给定测试样本$x$，其预测值是$x$所在叶子节点的训练样本集合的Y均值，形式化表示如下：
μ ^ ( x ) = 1 B ⋅ ∑ b = 1 B μ ^ b ( x ) ,   μ ^ b ( x ) = ∑ { i : X i ∈ L b ( x ) } Y i ∣ L b ( x ) ∣ (5) \hat\mu(x)=\frac{1}{B}\cdot\sum_{b=1}^B\hat\mu_b(x), \ \hat\mu_b(x)=\frac{\sum_{\{i:X_i\in L_b(x)\}}Y_i}{|L_b(x)|} \tag{5} μ^​(x)=B1​⋅b=1∑B​μ^​b​(x), μ^​b​(x)=∣Lb​(x)∣∑{i:Xi​∈Lb​(x)}​Yi​​(5)
现在，我们只需要证明当scoring function为${\Psi }_{\mu (x)}(O_i)=Y_i-\mu (x)$时，Eq(5)成立是Eq(1)成立的充要条件。 证明如下：
Eq(1)成立等价于Eq(6)成立：
$$\sum _{i=1}^n{\alpha }_i(x)\cdot (Y_i-\hat{\mu }(x))=0\text{\hspace{1em}(6)}$$
又由于${\sum }_{i=1}^n{\alpha }_i(x)=1$成立，所以Eq(6)可以转化成Eq(7)：
μ ^ ( x ) = ∑ i = 1 n α i ( x ) ⋅ Y i = ∑ i = 1 n 1 B ⋅ ∑ b = 1 B α b i ( x ) ⋅ Y i = 1 B ⋅ ∑ b = 1 B ⋅ ∑ i = 1 n 1 ( { X i ∈ L b ( x ) } ) ∣ L b ( x ) ∣ ⋅ Y i = 1 B ⋅ ∑ b = 1 B μ ^ b ( x ) (7) \begin{aligned} \hat\mu(x) &=\sum_{i=1}^n\alpha_i(x)\cdot Y_i \\ &=\sum_{i=1}^n \frac{1}{B}\cdot\sum_{b=1}^B\alpha_{b_i}(x) \cdot Y_i \\ &=\frac{1}{B}\cdot\sum_{b=1}^B\cdot\sum_{i=1}^n\frac{1(\{X_i\in L_b(x)\})}{|L_b(x)|} \cdot Y_i \\ &=\frac{1}{B}\cdot\sum_{b=1}^B \hat\mu_b(x) \end{aligned} \tag{7} μ^​(x)​=i=1∑n​αi​(x)⋅Yi​=i=1∑n​B1​⋅b=1∑B​αbi​​(x)⋅Yi​=B1​⋅b=1∑B​⋅i=1∑n​∣Lb​(x)∣1({Xi​∈Lb​(x)})​⋅Yi​=B1​⋅b=1∑B​μ^​b​(x)​(7)

由此可见，当Eq(1)成立时，能推出Eq(6)成立，因此，random forest是GRF的一个特例。

split criterion

最原始的思想是，最小化子节点评估值与真实值之间的误差，也就是最小化$err(C_1,C_2)$：
$$err(C_1,C_2)=\sum _{j=1}^2\mathbb{P}[X\in C_j|X\in P]\cdot \mathbb{E}[({\hat{\theta }}_{C_j}(\mathcal{J})-\theta (x){)}^2|X\in C_j]\text{\hspace{1em}(8)}$$
但是，由于真实值$\theta (x)$并不可知，因此，进过一番推导，我们将最小化$err(C_1,C_2)$转化成最大化$Delta(C_1,C_2)$：
$$\Delta (C_1,C_2)=\frac{n_{c_1}\cdot n_{c_2}}{n_p^2}\cdot ({\hat{\theta }}_{C_1}(\mathcal{J})-{\hat{\theta }}_{C_2}(\mathcal{J}){)}^2\text{\hspace{1em}(9)}$$
经过转化，可以发现，最大化Eq(7)的含义就是最大化子节点之间的异质性。

至此，我们知道了节点的分裂标准，但在实际操作中，由于${\hat{\theta }}_{C_j}(\mathcal{J})$的计算开销较大，因此，作者提出了基于gradient的近似求解方法：

gradient tree algorithm

首先，PROPOSITION1指出，${\hat{\theta }}_C$有如下近似解${\tilde{\theta }}_C$：
θ ~ C = θ ^ p − 1 ∣ { i : X i ∈ C } ∣ ⋅ ∑ { i : X i ∈ C } ξ T ⋅ A p − 1 Ψ θ ^ p , ν ^ p ( O i ) (10) \tilde\theta_{C}=\hat \theta_p-\frac{1}{|\{i:X_i\in C\}|}\cdot\sum_{\{i:X_i\in C\}}\xi^T\cdot A_p^{-1}\Psi_{\hat\theta_p,\hat\nu_p}(O_i) \tag{10} θ~C​=θ^p​−∣{i:Xi​∈C}∣1​⋅{i:Xi​∈C}∑​ξT⋅Ap−1​Ψθ^p​,ν^p​​(Oi​)(10)
其中，${\hat{\theta }}_p$表示父节点$P$上$\theta$的估计值，可以由Eq(1)orEq(2)求得；至于$\xi$，论文中说它是从$(\theta ,\nu )$向量中筛选出$\theta$-coordinate的向量，但我在其他论文中看到大家都省略了这个玩意儿；而$A_p$的含义是${\Psi }_{{\hat{\theta }}_p,{\hat{\nu }}_p}(O_i)$的期望的梯度，计算公式如下：
A p = ∇ E [ Ψ θ ^ p , ν ^ p ( O i ) ∣ X i ∈ P ] = 1 ∣ { i : X i ∈ C } ∣ ⋅ ∑ { i : X i ∈ P } ∇ Ψ θ ^ p , ν ^ p ( O i ) (11) A_p=\nabla\mathbb{E}[\Psi_{\hat\theta_p,\hat\nu_p}(O_i)|X_i\in P]=\frac{1}{|\{i:X_i\in C\}|}\cdot\sum_{\{i:X_i\in P\}}\nabla\Psi_{\hat\theta_p,\hat\nu_p}(O_i) \tag{11} Ap​=∇E[Ψθ^p​,ν^p​​(Oi​)∣Xi​∈P]=∣{i:Xi​∈C}∣1​⋅{i:Xi​∈P}∑​∇Ψθ^p​,ν^p​​(Oi​)(11)
但我不太理解这里的导数是对谁求的。

当${\hat{\theta }}_C$有近似解${\tilde{\theta }}_C$之后，可以推出$\Delta (C_1,C_2)$也有相应的近似解$\tilde{\Delta }(C_1,C_2)$：【这一步的推导暂时省略】
Δ ~ ( C 1 , C 2 ) = ∑ j = 1 2 1 ∣ { i : X i ∈ C j } ∣ ⋅ ( ∑ { i : X i ∈ C j } ρ i ) 2 (12) \tilde\Delta(C_1,C_2)=\sum_{j=1}^2\frac{1}{|\{i:X_i\in C_j\}|}\cdot(\sum_{\{i:X_i\in C_j\}}\rho_i)^2 \tag{12} Δ~(C1​,C2​)=j=1∑2​∣{i:Xi​∈Cj​}∣1​⋅({i:Xi​∈Cj​}∑​ρi​)2(12)
其中，${\rho }_i=-{\xi }^T\cdot A_p^{-1}\cdot {\Psi }_{{\hat{\theta }}_p,{\hat{\nu }}_p}(O_i)$，表示第i个样本在计算${\hat{\theta }}_p$时的影响。

至此，我们就可以将节点分裂总结为以下两个步骤：
1. labeling step
这一步，首先需要计算${\hat{\theta }}_p$和$A_p$，进而计算${\rho }_i$；注意，每次分裂时，只需要计算一个${\rho }_i$【因为父节点已经确定了】
2. regreession step
寻找子节点，使得$\tilde{\Delta }(C_1,C_2)$最大。这一步可通过标准CART回归分裂实现。

GRF for CATE

接着，我们看一下GRF是如何应用于CATE的评估的。
在这一应用中，作者仍以Partially Linear model为基础来构造$\Psi (\cdot )$，所谓的Partially Linear Model是指数据满足以下结构：
$$Y=\theta (x)\cdot T+g(x)+ϵ,T=f(x)+\eta \text{\hspace{1em}(13)}$$
所谓的”部分线性“主要体现在Y的结构上。
放在CATE评估问题中，$\theta (x)$就表示$x$条件下的处理效应，形式化表述为$\theta (x)=\mathbb{E}[Y(T=1)-Y(T=0)|X=x]$。
基于Partially Linear Model，作者构造的scoring function为${\Psi }_{\theta (x),\nu (x)}(O_i)=Y_i-\theta (x)\cdot T_i-\nu (x)$，可以理解为这个scoring function的目的是寻求一个$(\hat{\theta }(x),\hat{\nu }(x))$使得$Y_i$与$\theta (x)\cdot T_i+\nu (x)$尽可能接近【本质就是拟合问题】。

在这个设定下，各值求解如下：
$$\hat{\theta }(x)={\xi }^T\cdot \frac{Cov(T_i,Y_i|Xi=x)}{Var(T_i|X_i=x)}\text{\hspace{1em}(14)}$$
A p = 1 ∣ { i : X i ∈ P } ∣ ⋅ ∑ { i : X i ∈ C j } ( T i − T ˉ p ) ⨂ 2 (15) A_p=\frac{1}{|\{i:X_i\in P\}|}\cdot\sum_{\{i:X_i\in C_j\}}(T_i-\bar T_p)^{\bigotimes 2} \tag{15} Ap​=∣{i:Xi​∈P}∣1​⋅{i:Xi​∈Cj​}∑​(Ti​−Tˉp​)⨂2(15)
$${\rho }_i={\xi }^T\cdot A_p^{-1}\cdot (Y_i-{\bar{Y}}_p-(T_i-{\bar{T}}_p)\cdot {\hat{\theta }}_p)\text{\hspace{1em}(16)}$$
关于这些值的推导，目前只理解了Eq(14)的来源：
考虑$Y=\theta (x)\cdot T+g(x)$，最优$\theta (x)$的求解可以看作是求解一元一次方程$y=ax+b$的斜率，而这一斜率可以由方差及协方差表示【参考资料】
需要注意的是，这里的均值、方差、协方差都是加权计算的，而权重就是${\alpha }_i$。

CausalForestDML

顾名思义，CausalForestDML是融合了CausalForest和DML。DML在估计CATE时的核心思想是基于如下等式：
$$Y-\mathbb{E}[Y|X]=\theta (x)\cdot (T-\mathbb{E}[T|X])等价于\tilde{Y}=\theta (x)\cdot \tilde{T}\text{\hspace{1em}(17)}$$
也就是，将CATE的评估问题，转化成用T的残差去拟合Y的残差，而这个回归系数就是CATE。【计算残差的过程，其实就是正交化的过程】
基于DML的思想，CausalForestDML构造了如下的scoring function：${\Psi }_{\theta (x),\nu (x)}(O_i)=Y_i-\mathbb{E}[Y_i|X]-\theta (x)\cdot (T_i-\mathbb{E}[T_i|X])-\nu (x)$。
相应的最优$\theta (x)$就变成了$\hat{\theta }(x)={\xi }^T\cdot \frac{Cov(Y_i-\mathbb{E}[Y_i|X_i],T_i-\mathbb{E}[T_i|X_i]|Xi=x)}{Var(T_i-\mathbb{E}[Y_i|X_i]|X_i=x)}$。

注意：原文将这种方法称为Centered GRF！且表1结果表明GRF with centering效果比GRF without centering效果要好。

GRF vs Causal Forest

这两篇都是Susan的，GRF相比于Causal Forest最大的区别在于

1. Casual Forest使用exact loss criterion（即Eq(9)）进行分裂而非gradient-based loss criterion（即Eq(12))；
2. Causal Forest再计算treatment effect时没有采用加权平均；