论文复现：使用因果森林估计处理效应

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
  . ssc install lianxh
详情参见帮助文件 (有惊喜)：
  . help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg
 

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

• 全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
• 计量专题 | 论文写作 | 数据分享 | 专题课程
• 结果输出 | Stata绘图 | 数据处理 | Stata程序
• 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
• 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
• Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
• 文本分析-爬虫 | Python-R-Matlab | 机器学习
• Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：王颖 (四川大学)
邮箱：wangyingchn@outlook.com

编者按：本文主要整理自下文，特此致谢！
Source：Athey, S., & Wager, S. (2019). Estimating treatment effects with causal forests: An application. Observational Studies, 5(2), 37-51. -PDF- -Replication File-

目录

1. 研究概要

本文将因果森林应用于 National Study of Learning Mindsets 数据集。面对可能存在选择偏差和聚类偏误的数据，因果森林能够通过估计倾向得分和使用稳健的聚类方法较好的解决这两个挑战。

2. 数据背景

数据背景：国家学习心态研究 (National Study of Learning Mindsets) 是一项在美国高中进行的研究。这项研究对学生进行心态干预，向学生灌输成长型学习心态 (相对与固定型学习心态)，并评估这种干预对学生成绩的影响。

数据来源：来自 76 所美国高中 10,391 位学生的调查研究。

研究问题：

• 这种学习心态干预是否能有效提高学生成绩？
• 干预效果是否受到学校平均成绩水平 (X2) 或学生先前学习心态 (X1) 的调节？
• 干预效果是否受到其他协变量的调节？

变量说明：

需要注意的是，

• 为保护学生隐私，作者提供的复现数据不是原始数据，而是作者提供的具有类似结构的合成数据。
• 论文原文中处理变量使用名称 W，但在论文提供的数据中，处理变量名称为 Z。
• 论文对有两个以上水平的分类变量 (C1 和 XC) 进行了哑变量处理，因此最终有 28 个协变量。
• 论文中使用的参数在新版本的包中已经取消，修改为其他参数，因此复现结果可能与论文中的结果不同，此处直接展示论文原文中的结果。

3. 面临的两个挑战及解决方法

3.1 存在选择偏差：估计倾向得分

在因果推断中，要估计因果效应，需要满足条件独立性假设。即在给定的协变量 Xi 下，个体是否接受干预分配 (Wi) 与潜在的结果 Y 不相关。然而，通过数据发现，本文中的学习心态研究并不满足该假设，因为成功期望更高的学生 (S3，以往成绩的代理变量) 更可能接受干预，存在选择偏差的问题。当该假设不满足时，需要采用方法来消除这种潜在的干扰。

在因果森林中，通过估计倾向得分的方式来处理这个问题。倾向得分是某个个体接受干预的可能性得分。具体而言，在分析过程中考虑两个调整。一是用结果变量 Y 的真实值 (即观测值) 减去 Y 的期望值 (即无论个体是否接受干预时可能的 Y 值)；二是用处理变量 W 的真实值减去倾向得分。

总结而言，本文的数据集中可能存在选择偏差问题，处理方法是在因果森林中估计倾向得分，从而得到更加稳健的结果。

3.2 存在聚类偏误：聚类稳健分析

本文的数据中各个样本不是独立的，样本之间可能存在关联。由于学生来自 76 个学校，且不同的学校可能存在较大的异质性，可能存在有一些不可观测的学校层面变量影响干预效果，因此本文假设同一所学校内各个学生的结果变量 Yi 之间是相关的，并采用稳健的聚类分析 (cluster-robust analysis)。

当涉及来自不同群体的样本时，需要考虑以下两个问题：

• 是应该拟合一个能够准确估计数据中 76 所学校异质性的模型？还是拟合一个应用到其他学校的模型？
• 对于那些样本更多的学校，是否应该在模型中赋予更多的权重？

在本研究中，我们希望结果能推广到其他学校，且赋予每个学校一样的权重，即希望模型能够准确预测一个来自新学校的新学生的结果。

因果森林是将随机森林应用到异质处理效应估计上的一种算法。首先看看如何在随机森林中实现聚类稳健的分析。传统的随机森林，先随机抽取一组样本，然后构造树模型，重复多次构建多个树模型形成森林，最后进行预测。这里的预测使用袋外预测 (out-of-bag prediction)，即用于预测的样本是训练集之外的样本。

而要考虑聚类效应时，并不是随机抽取样本。以本文的数据为例，先从 76 所学校中随机抽取 x 个学校，然后随机从选定的 x 个学校中每个学校抽取 k 个样本形成 x*k 个样本构建训练集，基于该训练集构造树模型，并重复多次构建森林，最后做出预测。

同样，考虑到学校内部可能存在的关联，这里的预测使用类外预测，即用于预测的样本是来自用于训练的 x 个学校之外的样本。通过以上方法，可以在随机森林中考虑聚类效应，进行聚类稳健的分析。

总结而言，本文的数据集中可能存在聚类偏误问题。处理方法是，在构建随机森林的抽样、预测两个阶段都进行调整，进行稳健的聚类分析。同时，为了使结果能够外推到其他类别，对每个类别都赋予一样的权重。

4. 分析过程及结果 (R 实现)

在了解数据背景及要解决的问题后，使用 R 语言进行数据分析，通过 grf 包构建因果森林模型。

4.1 数据准备和模型构建

首先是数据准备，包括数据读取、构造变量等。

# 基本设置
set.seed(1)  # 设定随机数
rm(list = ls())  # 清空工作空间

# 安装相关包
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  grf  # 用于构造因果森林模型
  )

# 数据预处理
data.all <- read.csv('synthetic_data.csv')
data.all$schoolid = factor(data.all$schoolid)

DF <- data.all[,-1]
school.id = as.numeric(data.all$schoolid)  # 提取 school.id 作为后续聚类的标签
school.mat = model.matrix(~ schoolid + 0, data = data.all) 
school.size = colSums(school.mat) 

# 准备用于构建因果森林模型的数据集，其中 W 为处理变量，Y 结果变量，X 包含所有协变量
W <- DF$Z
Y <- DF$Y
X.raw <- DF[, -(1:2)]

C1.exp <- model.matrix(~ factor(X.raw$C1) + 0)  # 转化为哑变量
XC.exp <- model.matrix(~ factor(X.raw$XC) + 0) 

X <- cbind(X.raw[, -which(names(X.raw) %in% c("C1", "XC"))], C1.exp, XC.exp)

其次，构建因果森林模型。

# 构建结果变量 Y 和处理变量 W 的预测模型。基于 school.id 聚类，且每个类别权重一样
Y.forest <- regression_forest(X, Y, clusters = school.id, equalize.cluster.weights = TRUE)  
Y.hat <- predict(Y.forest)$predictions
W.forest <- regression_forest(X, W, clusters = school.id, equalize.cluster.weights = TRUE)
W.hat <- predict(W.forest)$predictions

# 初步构建因果森林模型
cf.raw <- causal_forest(X, Y, W,
                       Y.hat = Y.hat, W.hat = W.hat,
                       clusters = school.id,
                       equalize.cluster.weights = TRUE)  
# 原文中 samples_per_cluster 的参数在新版本 grf 包中已经取消，通过 equalize.cluster.weights 设置。

# 选择重要特征，再次构建因果森林模型
varimp <- variable_importance(cf.raw)
selected.idx <- which(varimp > mean(varimp))
cf <- causal_forest(X[,selected.idx], Y, W,
                   Y.hat = Y.hat, W.hat = W.hat,
                   clusters = school.id,
                   equalize.cluster.weights = TRUE, 
                   tune.parameters = "all")

# 进行预测，得到估计值
tau.hat <- predict(cf)$predictions 

4.2 估计平均处理效应

基于因果森林模型，估计平均处理效应 (average treatment effect，ATE)。结果显示，心态干预的平均处理效应是积极的 (ATE 为 0.247)。

ATE <- average_treatment_effect(cf)
paste("95% CI for the ATE:", round(ATE[1], 3),
      "+/-", round(qnorm(0.975) * ATE[2], 3))

4.3 估计整体的异质性

接下来，估计平均处理效应的整体异质性，即心态干预的处理效应是否有显著的异质性，通过两种方法检验。第一种方法，基于估计值的中位数分为高低两组，比较高组和低组的 ATE 是否有显著差异。结果显示，两组之间的差异很小 (差值为 0.053)。

high_effect <- tau.hat > median(tau.hat)
ate.high <- average_treatment_effect(cf, subset = high_effect)
ate.low <- average_treatment_effect(cf, subset = !high_effect)
paste("95% CI for difference in ATE:",
      round(ate.high[1] - ate.low[1], 3), "+/-",
      round(qnorm(0.975) * sqrt(ate.high[2]^2 + ate.low[2]^2), 3))

第二种方法，是通过最佳线性预测 (best linear) 的方法，检验处理效应的异质性是否显著。结果显示，异质性并不显著 ($p$ = 0.294)。

test_calibration(cf)

4.4 分析 X1 与 X2 的影响

虽然整体异质性分析的结果不显著，但不代表不存在异质性。因此，进一步分析 X1 和 X2 对平均处理的影响。首先，基于 X1 的中位数分为高低两组比较。结果显示，平均处理效应在 X1 的两个水平上有显著性差异 ($p$ = 0.00349)。

dr.score = tau.hat + W / cf$W.hat *
  (Y - cf$Y.hat - (1 - cf$W.hat) * tau.hat) -
  (1 - W) / (1 - cf$W.hat) * (Y - cf$Y.hat + cf$W.hat * tau.hat)
school.score = t(school.mat) %*% dr.score / school.size

# 基于 X1 中位数分为高低两组比较
school.X1 = t(school.mat) %*% X$X1 / school.size
high.X1 = school.X1 > median(school.X1)
t.test(school.score[high.X1], school.score[!high.X1])

然后，分析 X2 的影响，将 X2 划分为两组或三组比较。结果都显示，平均处理效应在 X2 上没有显著的异质性。

# 基于 X2 中位数分为高低两组比较
school.X2 = (t(school.mat) %*% X$X2) / school.size
high.X2 = school.X2 > median(school.X2)
t.test(school.score[high.X2], school.score[!high.X2])

# 分为三组比较
school.X2.levels = cut(school.X2,
  breaks = c(-Inf, quantile(school.X2, c(1/3, 2/3)), Inf))
summary(aov(school.score ~ school.X2.levels))

由于 X1 和 X2 都是学校层面的变量，如果只关注学校层面的变量是否可能更准确呢？文章进一步探索，只使用学校层面的协变量构建模型，结果与上面的结果一致。同时，文章使用学校层面的变量进行线性回归，发现只有 X1 在 0.10 的水平上显著。因此，如果考虑更多学校层面的变量，X1 的影响可能也是不存在的。

5. 补充分析

根据前文的说明，由于数据可能存在潜在的选择偏差和聚类偏误问题，前面的因果森林分析考虑了这两个问题，使用了正交化 (即通过估计倾向得分，使得数据满足独立性假设)，以及稳健的聚类方法。在补充分析中，作者比较了使用和不使用这两个处理方法的模型，以突出这两种方法的价值。

5.1 考虑聚类的价值

前文的分析基于 school.id 进行聚类，且每个学校权重一致，这里我们不考虑聚类，重新训练因果森林模型。

# 不考虑聚类，重新训练因果森林模型
cf.noclust = causal_forest(X[,selected.idx], Y, W,
                           Y.hat = Y.hat, W.hat = W.hat,
                           tune.parameters = "all")

ATE.noclust = average_treatment_effect(cf.noclust)
paste("95% CI for the ATE:", round(ATE.noclust[1], 3),
      "+/-", round(qnorm(0.975) * ATE.noclust[2], 3))

test_calibration(cf.noclust)

结果显示，处理效应具有显著的异质性。然而，这种显著的异质性主要是由于学校产生的，没有考虑到不同的学校内学生的相关性，不同学校间的差异很大，而且这种学校效应是特殊的，无法通过协变量来解释。因此，如果要将模型应用到新学校，就需要按学校分组，而此时学校之间的异质性就会像噪音一样干扰结果。

5.2 考虑正交的价值

前面的因果森林中，基于协变量进行了倾向得分估计，这里我们不估计倾向得分，重新训练因果森林模型，并比较结果的差异。

# 不估计倾向得分 (即 W.hat) 
cf.noprop = causal_forest(X[,selected.idx], Y, W,
                          Y.hat = Y.hat, W.hat = mean(W),
                          tune.parameters = "all",
                          equalize.cluster.weights = TRUE,
                          clusters = school.id)
tau.hat.noprop = predict(cf.noprop)$predictions

ATE.noprop = average_treatment_effect(cf.noprop)
paste("95% CI for the ATE:", round(ATE.noprop[1], 3),
      "+/-", round(qnorm(0.975) * ATE.noprop[2], 3))

结果显示，估计了倾向得分和未估计倾向得分的模型结果差异不大。原因可能是影响倾向的因素 (S3) 同时也是预测结果变量的重要因素，也就是说，在这个数据中，没有只对干预分配 (倾向) 有显著影响的协变量。

因此，为了进一步突出因果森林中估计倾向得分的价值，作者利用了新的合成数据来比较估计倾向得分与不估计倾向得分的模型。结果显示，两个模型估计的结果有很大差异，估计倾向得分的模型计算的 ATE 值为 0.125，而不估计倾向得分的模型 ATE 为 0.220。也就是说，如果数据中存在影响干预倾向的协变量，那么不考虑倾向得分可能会导致结果有很大的偏差。

# 生成新的数据
n.synth = 1000
p.synth = 10
X.synth = matrix(rnorm(n.synth * p.synth), n.synth, p.synth)
W.synth = rbinom(n.synth, 1, 1 / (1 + exp(-X.synth[,1])))
Y.synth = 2 * rowMeans(X.synth[,1:6]) + rnorm(n.synth)

# 结果变量与处理变量的预测模型
Y.forest.synth = regression_forest(X.synth, Y.synth)
Y.hat.synth = predict(Y.forest.synth)$predictions
W.forest.synth = regression_forest(X.synth, W.synth)
W.hat.synth = predict(W.forest.synth)$predictions

# 考虑倾向得分的因果森林模型
cf.synth = causal_forest(X.synth, Y.synth, W.synth,
                         Y.hat = Y.hat.synth, W.hat = W.hat.synth)
ATE.synth = average_treatment_effect(cf.synth)
paste("95% CI for the ATE:", round(ATE.synth[1], 3),
      "+/-", round(qnorm(0.975) * ATE.synth[2], 3))

# 不考虑倾向得分的因果森林模型
cf.synth.noprop = causal_forest(X.synth, Y.synth, W.synth,
                                Y.hat = Y.hat.synth, W.hat = mean(W.synth))
ATE.synth.noprop = average_treatment_effect(cf.synth.noprop)
paste("95% CI for the ATE:", round(ATE.synth.noprop[1], 3),
      "+/-", round(qnorm(0.975) * ATE.synth.noprop[2], 3))

6. 讨论

本文将因果森林应用到国家学习心态研究的数据集上，以分析异质处理效应。在这个数据集中，面临两个分析的挑战，一是样本受到干预分配的倾向不一样，二是可能受到聚类结果的影响。因果森林能够解决这两个挑战，一是通过估计倾向得分来调整，二是通过稳健的聚类分析来调整。

在本文的数据中，稳健的聚类分析有很明显的作用，如果不考虑聚类效应可能会产生有误的异质性分析结果。因此，在考虑异质处理效应时，需要更小心、更深入的讨论聚类问题。更好的理解进行聚类数据的异质处理效应估计的不同方法是很有意义的。

7. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 因果推断 随机森林, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：论文写作

  • 论文中因果推断的经典图形
• 专题：计量专题

  • 因果推断：哪本教材适合我？
  • 因果推断新书在线读：Causal Inference-The Mixtape
• 专题：Stata程序

  • Stata：随机森林算法简介与实现
• 专题：内生性-因果推断

  • 因果推断：双重机器学习-ddml
  • 一组动图读懂因果推断
  • 因果推断：混杂因素敏感性分析实操(下)-tesensitivity
  • 因果推断：混杂因素敏感性分析理论(上)
  • Stata因果推断：hettreatreg-用OLS估计异质性处理效应
  • Stata：因果推断方法综述和Stata操作
  • fect：基于面板数据的因果推断（上）-T218a
  • fect：基于面板数据的因果推断（下）-T218b
  • 因果推断：未测量混杂因素的敏感性分析-T249
  • 用FE-固定效应模型能做因果推断吗？
  • 经典文献回顾：政策评价-因果推断的计量方法
  • 因果推断好书：Causal-Inference-Measuring-the-Effect-of-X-on-y
  • Stata因果推断新书：The-SAGE-Handbook-of-Regression-Analysis-and-Causal-Inference
  • Stata新命令：konfound - 因果推断的稳健性检验
• 专题：其它

  • 50问-T2：面板数据因果推断常见问题-对话徐轶青老师
• 专题：最新课程

  • ✅ 江艇专场：深度因果推断
  • ⏩ 视频课：邱嘉平-因果推断实用计量方法
  • ⏫NEW-面板数据因果推断：从入门到精通
• 专题：机器学习

  • 机器学习：随机森林算法的Stata实现

相关课程

免费公开课

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版
• Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

最新课程-直播课

• Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

⛳ 课程主页

⛳ 课程主页

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 连享会-主页 和 知乎专栏，700+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会-常见问题解答：
✨ https://gitee.com/lianxh/Course/wikis

New！ lianxh 和 songbl 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh