最近阅读了文章：Using Machine Learning to Target Treatment The Case of Household Energy Use，写下此文，记录文章的主要内容和R语言grf包causal forest方法的简单过程。

一、文章的基本内容和方法

引言

（1）背景

Machine-learning algorithms have traditionally been built for prediction of from y, rather than parameter estimation of treatment effect β.（机器学习方法主要用于预测y值，而不是β值。）即普通的机器学习方法不能用于因果推断。因此，本文希望可以提出一种可以用于因果推断的机器学习方法，即因果森林算法（causal forest algorithm），依托于R语言generalized random forest 包（grf包，广义随机森林模型）。

（2）目的

文章基于美国家庭的90万能源使用数据，识别“节能宣传行为-在家庭能源报告中提供节能信息”是否可以推动家庭节约能源使用，量化政策效果，探究家庭特征对结果的影响。

方法

（1）数据

美国新英格兰最大的电力公司Eversource提供家庭能源报告（Home Energy Report, HER），在报告中，公司会不定期的提供节能的相关信息。研究地区和用能结构如下方左图。结合了三类数据：来自Eversource的家庭月用电量；Eversource政策刺激的时间安排（treatment）；以及消费者人口和社会经济特征。

本次研究的数据利用了15波政策刺激数据，覆盖了902581个Eversource住宅客户。观察了2013-2018年的每月家庭用电量以及与家庭及其住户相关的横截面特征。其中，部分月份，该公司会提供如下的信息（下方右图），告诉消费者有效能耗、你的能耗和你邻居的能耗的差别，并提供一些节能方案（被视为政策刺激）。

 15波政策刺激的情况如下：

（2）DID

使用传统的双重差分法（DID）回归，将住户分配到实验组和对照组，来估计HER政策刺激对电力消耗的影响：

kWh-用电量

T-treatment vector，实验组为1，对照组为0

X-家庭的特征变量

θ-wave的固定效应

w-年-月固定效应

还可以估计HER计划的累计影响： 

j-the beginning of treatment in the relevant wave，j∈[-12,37]

D-treatment matrix，实验组为1，对照组为0

（3）Causal forest

为了测试HER计划的异质性影响，并调查家庭特征在预测中的作用，使用了因果森林算法。因果森林算法（Athey等人，2019年）是对随机森林（Breiman，2001年）的一种调整，用于测量因果效应。传统的回归树算法使用相同的数据集来扩展树结构和估计每个节点的average treatment effect（ATE）。然而，Athey和Imbens（2016）表明，这种做法往往夸大了拟合优度，因此他们引入了“honest estimate”，即将全部随机子样本一分为二，一个子集用于生长树结构，另一个子集用于估计ATE。因果随机森林的相关介绍见：使用随机森林进行因果推断 - 知乎

结果和讨论

（1）DID的结果

HER政策Wave-specific ATE的估值范围为每月-1.6至-17.7 kWh，Pooled ATE为每月-8.85 kWh，或-1%。HER政策伴随时间的影响如下图，政策开始前的12个月，没有估计值显著区别于0，政策开始的1-2个月没有显著影响，之后愈发显著，一直到第三年也有增效。

 （2）Causal forest的结果

因为采用了“森林”的计算方法，可以根据各个树的结果得到ATE估值的概率分布函数。其中第一年的效果集中在-10kWh和0，说明一部分奏效，一部分没有影响。第二、三年的政策效果的范围更加广泛，但是出现了18%的用户增加了能源消耗（回力镖效应，因为看到自己比较节能，反而在下一个时段放松了警惕）。不同用户特征对于估值结果的影响，也得到了很好的量化（略）。

下图显示了预计减少消费的家庭与预计增加消费的家庭的平均差异。森林中使用的每一个特征都存在显著差异：“reducers”更有可能拥有自己的住所；他们的房子往往更大、更新、更有价值；而且他们往往有更多的孩子、更高的收入、更年轻、受教育程度更高的户主。最后，那些减少用电量的人的平均基线用电量显著高于其他人，而且这种差异在数量上占主导地位。这一发现验证了Opower针对高消耗用户的寄送HER的策略，并且与基线消耗量较大时有更多节能“空间”的想法一致。

 

二、grf包初探

参考资料：Introduction to grf • grf

如下为一个简单的案例，其中X代表自变量（训练集），Y代表因变量（训练集），W为treatment变量（一般为0/1），X.test为测试集。

library(grf)
n <- 2000
p <- 10
X <- matrix(rnorm(n * p), n, p)
X.test <- matrix(0, 101, p)
X.test[, 1] <- seq(-2, 2, length.out = 101)

W <- rbinom(n, 1, 0.4 + 0.2 * (X[, 1] > 0))
Y <- pmax(X[, 1], 0) * W + X[, 2] + pmin(X[, 3], 0) + rnorm(n)
tau.forest <- causal_forest(X, Y, W)
tau.forest
#> GRF forest object of type causal_forest 
#> Number of trees: 2000 
#> Number of training samples: 2000 
#> Variable importance: 
#>     1     2     3     4     5     6     7     8     9    10 
#> 0.691 0.033 0.041 0.046 0.029 0.029 0.034 0.027 0.033 0.037

使用OOB预测估计训练集的treatment效果：

tau.hat.oob <- predict(tau.forest)
hist(tau.hat.oob$predictions)

  评估测试集的treatment效果：

tau.hat <- predict(tau.forest, X.test)
plot(X.test[, 1], tau.hat$predictions, ylim = range(tau.hat$predictions, 0, 2), xlab = "x", ylab = "tau", type = "l")
lines(X.test[, 1], pmax(0, X.test[, 1]), col = 2, lty = 2)

 估计全样本的conditional average treatment effect (CATE):

average_treatment_effect(tau.forest, target.sample = "all")
#>   estimate    std.err 
#> 0.37345648 0.04978157

估计treatment样本的conditional average treatment effect (CATT):

average_treatment_effect(tau.forest, target.sample = "treated")
#>   estimate    std.err 
#> 0.47136673 0.05136118

增加置信区间:

tau.forest <- causal_forest(X, Y, W, num.trees = 4000)
tau.hat <- predict(tau.forest, X.test, estimate.variance = TRUE)
sigma.hat <- sqrt(tau.hat$variance.estimates)
plot(X.test[, 1], tau.hat$predictions, ylim = range(tau.hat$predictions + 1.96 * sigma.hat, tau.hat$predictions - 1.96 * sigma.hat, 0, 2), xlab = "x", ylab = "tau", type = "l")
lines(X.test[, 1], tau.hat$predictions + 1.96 * sigma.hat, col = 1, lty = 2)
lines(X.test[, 1], tau.hat$predictions - 1.96 * sigma.hat, col = 1, lty = 2)
lines(X.test[, 1], pmax(0, X.test[, 1]), col = 2, lty = 1)