MLRtime：如何在 Stata 调用 R 的机器学习包？

👇 连享会 · 推文导航 | www.lianxh.cn

🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令
📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享
💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归
⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析
🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
🔨 工具：工具软件 | Markdown | Python-R-Stata
🎧 课程：公开课-直播 | 计量专题 | 关于连享会

连享会课程 · 文本分析专题

作者：吕卓阳 (厦门大学)
邮箱：lvzy20@163.com

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

1. 配置环境
2. MLRtime 详解

2.1 causal_forest 与 ranger
2.2 MLRtime主要函数介绍

3. 具体示例
4. 结论
5. 参考资料
6. 相关推文

R 是统计和数据分析应用中的一件利器，它具有免费开源、语言易懂、大量成熟的库与工具等优点。在「Stata 与 R 的无缝对接」一文中，我们已经介绍了如何通过 Rcall 在 Stata 中运行 R 语言。在本文中，我们将继续为大家介绍新的在 Stata 中调用 R 的包 MLRtime。该包允许在 Stata 中使用 R 的 grf 包与 ranger 包，从而实现随机森林的因果推断与预测功能。

1. 配置环境

我们需要配置环境以方便在 Stata 中调用 R 以实现机器学习。首先，需要通过以下命令安装 MLRtime 包：

*方法一
cnssc install lxhget, replace
lxhget mlrtime.pkg, install replace

*方法二
net install MLRtime, from("https://raw.githubusercontent.com/NickCH-K/MLRtime/master/")

其次，我们需要在 Stata 中配置 rcall 环境，主要步骤如下：

在电脑中安装 R 或 Rstudio 软件，可以在官网中直接免费下载；
使用 github package 下载 rcall:

net install github, from("https://haghish.github.io/github/")
github install haghish/rcall, stable

最后，需要在 Stata 命令框中输入 MLRtimesetup, go，以确保 Stata 所需的包与 rcall (Haghish, 2019) 都已正确下载安装。完成以上设置后，就可以在 Stata 中调用 R 了！我们每次只需要在 Stata 的命令框中输入 rcall，就可以开始在 Stata 中运行 R 代码。由于 rcall 的基本功能已在「Stata 与 R 的无缝对接」中有详细介绍，本文将着重介绍使用 MLRtime 来进行机器学习建模。

rcall clear
sysuse auto.dta, clear
causal_forest price foreign mpg rep78 headroom trunk weight length turn, ///
  clearR replace seed(1002) pred(effect) opts(num.trees=50)              ///
  varreturn(effectSE = sqrt(predict(CF, X, estimate.variance = TRUE)@@variance.estimates)) ///
  return(ate = average_treatment_effect(CF)[1];ateSE = average_treatment_effect(CF)[2])

2. MLRtime 详解

MLRtime 的两个主要函数是 causal_forest 与 ranger，本质是调用 R 中的 causal_forest 函数进行随机森林因果推断与 ranger 函数进行随机森林预测。首先，我们先介绍一下R中的这两个函数。

R 的 causal_forest 是 grf 包中用于随机森林因果推断的一个函数，本质是使用随机森林的算法来进行回归，估计得到平均处理效应。grf 包是基于随机森林进行统计估计与推断的一个集成包，它为最小二乘回归、分位数回归、生存分析与处理效应估计提供非参数方法。R 的 ranger 是 ranger 包中用于随机森林因果预测的一个函数。ranger 包是 Wright 和 Ziegler (2017) 开发的高维数据随机森林的一个 R 包，可以高效的实现随机森林预测，支持分类、回归、生存森林的实现。

2.1 causal_forest 与 ranger

我们可以在 R document 的「官网」搜索 causal_forest，点击进入函数介绍页面，causal_forest 因果森林可以用于估计平均处理效应。下面将简单介绍 R 的 grf 包中 causal_forest 函数的主要参数，其余参数详见「R document」。

causal_forest(
  X,
  Y,
  W,
  Y.hat = NULL,
  W.hat = NULL,
  num.trees = 2000,
  sample.weights = NULL,
  clusters = NULL
)

X：用于估计因果森林的协变量；
Y：结果变量，必须是无 NA 的数值型变量；
W：处理变量，必须是无 NA 的二元或数值型变量；
Y.hat ：处理效应的边际影响，默认值为 NULL；
W.hat：处理效应的倾向匹配，默认值为 NULL；
num.trees：随机森林数量，默认值为 2000；
sample.weights：给予每个样本的权重，默认值为 NULL；
clusters：标准误集聚，默认值为 NULL。

R 的 ranger 包中的ranger函数可以用于随机森林的快速实现，尤其适用于高维数据，主要参数介绍如下，其余参数可以详见「R document」。

ranger(
  formula = NULL,
  data = NULL,
  num.trees = 500,
  mtry = NULL,
  importance = "none"
)

formula：随机森林拟合模型，仅适用于数值型变量；
data：训练数据集，一般为 dataframe、matrix、dgCMatrix、gwaa.data；
num.trees：随机森林数量；
mtry：每个节点中可能拆分的变量数，默认值为平方根；
importance：变量重要性，可选 none、impurity、impurity_corrected、permutation。

2.2 MLRtime主要函数介绍

随机森林是基于树模型的一种分类方法，Susan 等 (2019) 基于 Breiman (2001) 的随机森林方法提出了一种广义随机森林方法。该方法可以用于估计局部矩条件识别的任意，并基于广义随机森林方法开发了非参数分位数回归、条件平均处理效应估计和使用工具变量的异质平均处理效应。Susan 等 (2019) 开发的 R 的 grf 包中causal_forest 就是用于估计随机森林的平均处理效应的函数。我们假设是二元无混淆的分配变量，定义平均处理效应如下：

当连续时，我们可以有效的估计平均偏效应作为无混淆的平均处理效应：

我们可以在 Stata 中调用 MLRtime 包的 causal_forest 函数，来调用 R 的 grf 包中的 causal_forest 以实现随机森林因果推断。函数及参数含义如下：

causal_forest depvar indvar varlist [if] [in] [weight], [options]

clearR：运行前启动一个新的 R 环境；
replace：替换内存中同名变量；
seed(integer)：在 R 中设置随机数种子；
pred(string)：存储预测效果的变量名，需要满足 R 和 Stata 的命名条件；
predopts(string)：一系列预测命令的选项，用逗号分割，主要参数如下，可以使用 rcall:help(predict.causal_forest,package='grf') 命令来进一步查看：



    
predict(object,newdata = NULL,linear.correction.variables = NULL,
ll.lambda = NULL,ll.weight.penalty = FALSE,
num.threads = NULL,estimate.variance = FALSE)

opts(string)：一系列参数，用逗号连接，使用 @@ 代替 $，否则 Stata 将返回全局变量；
varreturn(string)：返回新变量，以分隔，可以将随机森林对象 CF、协变量、处理和结果矩阵存储到 Stata 中，使用 @@ 代替 $，否则 Stata 将返回全局变量；
return(string)：返回变量，以分隔，格式为，为在 Stata 创建的新变量名称，为随机森林函数返回的 CF、协变量、处理和结果矩阵，譬如如果返回平均处理效应及其标准误，则输入以下命令：

(ate = average_treatment_effect(CF)[1]; ateSE = average_treatment_effect(CF)[2])

varsend(varlist)：发送到R中到变量列表
precedes：其他参数，详见rcall:help(grf::causal_forest)
follows：希望在R中运行并返回Stata前执行到操作，详见 rcall:help(grf::causal_forest)

Wright 和 Ziegler (2017) 提出了高维数据随机森林，提高了随机森林的运算效率，并开发了 R 的 ranger 包。ranger 函数可以高效率的实现随机森林算法，并支持分类和回归森林 (Breiman 2001)、生存森林 (Ishwaran 等 2008)。除了在 R 中直接调用 ranger 函数，我们还可以在 Stata 中通过 MLRtime 来调用这一函数，Stata 中 ranger 函数命令介绍如下：




    
ranger depvar varlist [if] [in] [weight], [options]

clearR：运行前启动一个新的 R 环境；
replace：替换内存中同名变量；
seed(integer)：在 R 中设置随机数种子；
pred(string)：存储预测效果的变量名，需要满足 R 和 Stata 的命名条件；
predopts(string)：一系列预测命令的选项，用逗号分割，主要参数如下，可以使用 rcall:help(predict.ranger,package='ranger') 命令来进一步查看：

predict(object,data = NULL,predict.all = FALSE,num.trees = object$num.trees,
type = "response",se.method = "infjack",
quantiles = c(0.1, 0.5, 0.9),what = NULL,seed = NULL,
num.threads = NULL,verbose = TRUE)

opts(string) 传递给 R 一系列参数，用逗号连接，使用 @@ 代替 $，否则 Stata 将返回全局变量；
varreturn(string)：返回新变量，以分隔，可以将随机森林对象 CF、协变量、处理和结果矩阵存储到 Stata 中，使用 @@ 代替 $，否则 Stata 将返回全局变量；
return(string)：返回变量，以分隔，格式为，为在 Stata 创建的新变量名称，为随机森林函数返回的 CF、协变量、处理和结果矩阵，譬如如果返回平均处理效应及其标准误，则输入：

(ate = average_treatment_effect(CF)[1]; ateSE = average_treatment_effect(CF)[2])

varsend(varlist)：发送到R中到变量列表；
precedes：其他参数，详见 rcall:help(grf::causal_forest)；
follows：希望在 R 中运行并返回 Stata 前执行到操作，详见 rcall:help(grf::causal_forest)。

可以看到 ranger 函数在 MLRtime 的形式与 causal_forest 完全相同，但两个函数的目的不同。前者是用于因果森林，后者是用于随机森林预测，尤其是对于高维数据。通过这两个函数，就可以实现在 Stata 中调用 R 进行随机森林算法的估计，其本质是调用 R 的 ranger 包和 grf 包。

3. 具体示例

我们使用 titanic 数据集进行因果森林与随机森林预测估计：

*使用 titanic 数据集样本内估计 ate
lxhget titanic.csv, replace
import delimited using titanic.csv, clear
sum
tabstat age pclass sibsp parch,by(survived)
gen female = 0
replace female = 1 if sex == "female"
gen name_length = strlen(name)
rcall clear
causal_forest survived age pclass sibsp female name_length,  ///
  clearR replace seed(1002) pred(effect) opts(num.trees=50) ///
  varreturn(effectSE = sqrt(predict(CF, X, estimate.variance = TRUE)@@variance.estimates)) ///
  return(ate = average_treatment_effect(CF)[1];ateSE = average_treatment_effect(CF)[2])
di r(ate)
di r(ateSE)


*使用 titanic 数据集样本内预测生存概率
lxhget titanic.csv, replace
import delimited using titanic.csv, clear

sum
tabstat age pclass sibsp parch,by(survived)
gen female = 0
replace female = 1 if sex == "female"
gen name_length = strlen(name)
rcall clear
ranger survived age pclass sibsp female name_length, clearR replace ///
  seed(1002) pred(pred_survived) predopts(data=rangerdata)          ///
  opts(num.trees=50, importance = 'permutation')                    ///
  return(importance = as.matrix(importance(RF)))
pwcorr survived pred_survived

4. 结论

MLRtime 是使用 Stata 调用 R 机器学习包的一个途径。尽管目前能实现的只有因果森林与随机森林预测，并且无法像 R 与 Python 一样可以方便的划分训练集与测试集，但是我们相信在未来的发展中，Stata 与 R 或其他软件的交互会越来越深入，可以实现的功能也越来越强大。

5. 参考资料

Athey S, Tibshirani J, Wager S. Generalized random forests[J]. The Annals of Statistics, 2019, 47(2): 1148-1178. -PDF-
Wright M N, Ziegler A. ranger: A fast implementation of random forests for high dimensional data in C++ and R[J]. arXiv preprint arXiv:1508.04409, 2015. -PDF-

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 随机森林机器学习交互, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题：Stata教程

Stata-Python交互-9：将python数据导入Stata
Stata-Python交互-8：将Stata数据导入Python
Stata-Python交互-7：在Stata中实现机器学习-支持向量机
Stata-Python交互-6：调用APIs和JSON数据
Stata-Python交互-5：边际效应三维立体图示
Stata-Python交互-4：如何调用Python宏包
Stata-Python交互-3：如何安装Python宏包
Stata-Python交互-2：在Stata中调用Python的三种方式
Stata-Python交互-1：二者配合的基本设定
Rcall：Stata 与 R 的无缝对接

专题：Stata程序

Stata：随机森林算法简介与实现

专题：Python-R-Matlab

Stata交互：Python-与-Stata-对比

专题：其它

知乎热议：纠结-计量经济、时间序列和机器学习

专题：机器学习

机器学习如何用？金融+能源经济学文献综述
机器学习：随机森林算法的Stata实现
Stata：机器学习分类器大全

课程推荐：计量和因果推断 · 强基班
主讲老师：司继春
课程时间：2023 年 11 月 5/12/19 (三个周日)
🍓 课程主页：https://www.lianxh.cn/news/b8936ea77cced.html

New！ Stata 搜索神器：lianxh 和 songbl GIF 动图介绍
搜：推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all

🍏 关于我们

连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
直通车： 👉【百度一下：连享会】即可直达连享会主页。亦可进一步添加「知乎」,「b 站」,「面板数据」,「公开课」等关键词细化搜索。