Rcall:Stata 与 R 的无缝对接

发布时间:2020-07-28 阅读 308

Stata 连享会   主页 || 视频 || 推文

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

课程详情 https://gitee.com/arlionn/Course   |   lianxh.cn

课程主页 https://gitee.com/arlionn/Course

谢雁翔 (南开大学,xyxmask1995@163.com)
郭晓慧 (山东大学,xiaohuiguo161@sina.com)

本文主要参考如下论文,特此致谢。

Haghish E F. Seamless interactive language interfacing between R and Stata[J]. The Stata Journal, 2019, 19(1): 61-82. -Link-


目录


1. 简介

R 是一门免费、开源的语言,具有体积小 (R4.0.0 为 83.5M)、功能多 (外部包多达 16052 个) 等优势,其被广泛应用于数据挖掘、机器学习、数据可视化、计量经济学和空间统计等领域。那么如何在 Stata 中运行 R 语言? --- 快使用 Rcall

Rcall 安装包可将 R 集成在 Stata 中,并允许 R 与 Stata 之间的数据衔接。它能够将 Stata 中的宏、标量、矩阵和数据集自动传输到 R,并将 R 中不同类别的对象 (数据框 data.frame、列表 list、矩阵 matrix、向量 vector、逻辑 logical 等) 自动导入到 Stata。

Stata 与 R 的结合可使 Stata 用户在 R 环境中交互执行的其他编程语言:例如使用 Rcpp 的 C ++,或使用 V8 软件包的 JavaScript。此外,Rcall 允许在 Stata 的 ado 程序中嵌入 R 语言,从而不仅可以调用 R 函数和程序包,还可以通过 Stata 和 R 之间的转换来对 Stata 程序包进行编程。

2. R 和 Rcall 命令的安装

2.1 R 的下载

Rcall 命令安装之前首先要安装 R。由于 R 是免费开源的软件,下载十分方便,可直接在其「官网下载」。R 的安装相对简单,这里不再赘述。具体操作如下图所示:

进入官网下载
进入官网下载
选择对应镜像
选择对应镜像
选择对应系统
选择对应系统
选择初次安装
选择初次安装
下载最新版本
下载最新版本

2.2 Rcall 命令的安装

github package 是安装 rcall 的唯一推荐方式。由于 github 网络不稳定,该安装过程需要花费大量时间。

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

一旦 github 被安装,输入:

github install haghish/rcall, stable 

3. Rcall 命令及 R 语言初识

3.1 Rcall 命令基本语法

Rcall 语法如下 (具体参见help rcall):

rcall [mode] [:] [R-command]

mode 模式具体有:vanillasyncinteractiveconsole 四种模式,默认交互模式。

此外,rcall 还包括一些子命令 subcommand,以方便将 R 集成到 Stata 中,语法如下:

rcall [subcommand]

以下函数可用于将数据从 Stata 传输到 R :

函数 英文描述 中文描述
st.scalar(name) passes a scalar to R 传输标量
st.matrix(name) passes a matrix to R 传输矩阵
st.var(varname) passes a numeric or string variable to R 传输数值或字符变量
st.data(filename) passes Stata data to R. without filename, the currently loaded data is used 将Stata数据传输到R。不带文件名,使用当前加载数据
st.load(dataframe) loads data from R dataframe to Stata 将数据从 R 的数据框加载到 Stata

3.2 R 语言基本的数据类型

数据类型 示例
数值型 (numeric) x = 2
字符型 (character) x = "2"
复数型 (complex) x= 1 + 3i
逻辑型 (logical) x=1; y=x>2
因子型 (factor) genders <- factor(c("男", "女", "女", "男", "男"))

3.3 R 语言基本的数据结构

向量 vector (一维):数据类型都可取,不允许出现不同数据类型

b = c(1,2,3,4)

矩阵 matrix (二维):数据类型都可取,不允许出现不同数据类型

matrix(c(1,2,3,4,5,6),nrow=3,byrow=T/F)  

数组 array (三维):数据类型都可取,不允许出现不同数据类型

aa=array(1:24,dim=c(3,4,2))  

数据框 data.frame (二维):数据类型都可取,不同列的数据类型可不同

M = data.frame(a1=c(1,2,3),b=c(“a”,”b”,”c”))

列表 list (“万能胶”):数据类型都可取,任何元素数据类型均可不同

x = 1:3 
y = c("A", "B", "c") 
z = c(TRUE, FALSE) 
a = c(x, y, z) 
L = list(x, y, z, a)

Note: 以上内容部分摘自「游万海:R 语言初识」

4. Stata实例

4.1 Stata 与 R 语言的数据转换

每当执行 rcall 时,Stata 都会自动将 R 对象作为 rclass 接收。如果 R 以交互方式运行 (即不使用 vanilla 子命令),则先前的对象仍可访问 Stata,除非将它们从 R 中更改或删除。此外,从 Stata 中加载的 R 软件包将保持加载状态,直到终止。在 Stata 中访问 R 对象是同时进行的,这使得使用 rcall 变得方便。例如,可以在 Stata 中访问 R 中定义的数字或字符串向量,就像使用 rclass 调用该对象的名称一样简单,即 r(objectname)

数值型案例:

rcall: a <- 100 
display r(a)      
100

如果没有 vanilla 子命令,则定义的对象将保留在 R 的内存中,因此,只要调用 R 便将其返回给 Stata。

rcall: a 
[1] 100 

字符型案例:

rcall: str <- "Hello World"  // <-为赋值符号
display r(str)  
Hello World   
rcall: str <- c("Hello", "World") 
display r(str) 
"Hello"  "World" 

向量案例:

rcall: v <- c(1,2,3,4,5)
display r(v)      
1 2 3 4 5

矩阵案例:

rcall: A = matrix(1:6, nrow=2, byrow = TRUE) 
mat list r(A) 
r(A)[2,3]

  c1  c2  c3     
r1  1  2  3
r2  4  5  6 

列表案例:

rcall: mylist <- list(a=c(1:10))
display r(mylist_a) 
1 2 3 4 5 6 7 8 9 10

逻辑性数据案例:

rcall: l <- T 
display r(l) 
TRUE     

空值型数据案例:

rcall: n <- NULL
display r(n)
NULL

为了实现 Stata 和 R 之间的理想交互,rcall 可以将 Stata 中的变量传递给 R。

可以在 R 代码中传递局部宏和全局宏,如下例所示:

global a 99 
rcall: (a <- $a)    
[1] 99 

为了将标量从 Stata 传递到 R ,可以使用 st.scalar() 函数,如下所示:

scalar a = 50 
rcall: (a <- st.scalar(a))   
[1] 50     

同样,可以使用 st.matrix() 函数将 Stata 矩阵传输到 R ,如下所示:

matrix A = (1,2\3,4) 
matrix B = (96,96\96,96)        
rcall: C <- st.matrix(A) + st.matrix(B)
rcall: C 
     [,1] [,2]
[1,]  97  98
[2,]  99  100   

当然,也可以调用 Stata 中的 R 矩阵:

mat list r(C)
r(C)[2,2]

    c1  c2
r1  97  98
r2  99  100

使用 st.var(varname) 函数将变量从 Stata 传递到 R 很方便。因此,只需将分析所需的变量从 Stata 传递到 R,即可在 R 中执行任何分析:

sysuse auto, clear 
(1978 Automobile Data)

rcall: dep <- st.var(price)      
rcall: pre <- st.var(mpg)   
rcall: lm(dep~pre)
Call:
lm(formula = dep ~ pre)

Coefficients:
(Intercept)          pre  
    11253.1       -238.9  

rcall 包还允许在 st.data(filename) 函数中将 Stata 数据传递到 R。此功能依赖于 R 中 readstata13 软件包来加载 Stata 数据集,而无需将其转换为 csv 或类似格式。readstata13 的 R 软件包可以按以下方式安装在 Stata 中:

rcall: install.packages("readstata13", repos="http://cran.uk.r-project.org")

指定数据集的相对或绝对路径,以将数据从 Stata 传输到 R。例如:

rcall: data <- st.data(/Applications/Stata/ado/base/a/auto.dta) 
rcall: dim(data)

若未指定文件名,则该函数将当前加载的数据传递给 R。

sysuse auto, clear 
(1978 Automobile Data)
rcall: data <- st.data() 
rcall: dim(data) 
[1] 74 12

最后,可以使用 st.load(dataframe) 函数将数据从 R 自动导入到 Stata 中。此功能将自动从 R 中保存一个 Stata 数据集,并通过清除当前数据集 (如果有的话) 将其加载到 Stata 中。当然,如果在 R 中编写了适当的代码以导出 Stata 数据集,则可以更好地控制转换变量类型。该功能在大多数情况下应该可以正常工作:

clear 
rcall: st.load(cars) 
list in 1/2
     +--------------+
     | speed   dist |
     |--------------|
  1. |     4      2 |
  2. |     4     10 |
     +--------------+

4.2 在 Stata 中运行 R

本部分以 Kleiber 和 Zeileis (2008) 的 Grunfeld.dta 数据集为例,分别使用 Stata 和 R 软件进行固定效应分析。对于面板数据的理解,请参考「Stata: 面板数据模型一文读懂」

clear
webuse grunfeld,clear //利用webuse从网络读取数据
list in 1/10          // 显示该数据集的前10行
xtset company year, yearly //设置面板数据格式
xtreg invest mvalue kstock, fe  //fe表示固定效应  
     +--------------------------------------------------+
     | company   year   invest   mvalue   kstock   time |
     |--------------------------------------------------|
  1. |       1   1935    317.6   3078.5      2.8      1 |
  2. |       1   1936    391.8   4661.7     52.6      2 |
  3. |       1   1937    410.6   5387.1    156.9      3 |
  4. |       1   1938    257.7   2792.2    209.2      4 |
  5. |       1   1939    330.8   4313.2    203.4      5 |
     |--------------------------------------------------|
  6. |       1   1940    461.2   4643.9    207.2      6 |
  7. |       1   1941      512   4551.2    255.2      7 |
  8. |       1   1942      448   3244.1    303.7      8 |
  9. |       1   1943    499.6   4053.7    264.1      9 |
 10. |       1   1944    547.5   4379.3    201.6     10 |
     +--------------------------------------------------+

       panel variable:  company (strongly balanced)
        time variable:  year, 1935 to 1954
                delta:  1 year

Fixed-effects (within) regression               Number of obs     =        200
Group variable: company                         Number of groups  =         10

R-sq:                                           Obs per group:
     within  = 0.7668                                         min =         20
     between = 0.8194                                         avg =       20.0
     overall = 0.8060                                         max =         20

                                                F(2,188)          =     309.01
corr(u_i, Xb)  = -0.1517                        Prob > F          =     0.0000

------------------------------------------------------------------------------
      invest |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      mvalue |   .1101238   .0118567     9.29   0.000     .0867345    .1335131
      kstock |   .3100653   .0173545    17.87   0.000     .2758308    .3442999
       _cons |  -58.74393   12.45369    -4.72   0.000    -83.31086     -34.177
-------------+----------------------------------------------------------------
     sigma_u |  85.732501
     sigma_e |  52.767964
         rho |  .72525012   (fraction of variance due to u_i)
------------------------------------------------------------------------------
F test that all u_i=0: F(9, 188) = 49.18                     Prob > F = 0.0000
rcall: install.packages("plm",repos ="http://cran.us.r-project.org")   //安装面板数据回归包

rcall:library(plm)   //加载包
rcall:data("Grunfeld", package="plm")  //调用自带数据集
rcall:head(Grunfeld,10)  //列出前十行

rcall:zz <- plm(inv ~ value + capital ,data = Grunfeld, index = c("firm")) 
rcall:summary(zz)
   firm year   inv  value capital
1     1 1935 317.6 3078.5     2.8
2     1 1936 391.8 4661.7    52.6
3     1 1937 410.6 5387.1   156.9
4     1 1938 257.7 2792.2   209.2
5     1 1939 330.8 4313.2   203.4
6     1 1940 461.2 4643.9   207.2
7     1 1941 512.0 4551.2   255.2
8     1 1942 448.0 3244.1   303.7
9     1 1943 499.6 4053.7   264.1
10    1 1944 547.5 4379.3   201.6

Oneway (individual) effect Within Model

Call:
plm(formula = inv ~ value + capital, data = Grunfeld, index = c("firm"))

Balanced Panel: n = 10, T = 20, N = 200

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-184.00857  -17.64316    0.56337   19.19222  250.70974 

Coefficients:
        Estimate Std. Error t-value  Pr(>|t|)    
value   0.110124   0.011857  9.2879 < 2.2e-16 ***
capital 0.310065   0.017355 17.8666 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2244400
Residual Sum of Squares: 523480
R-Squared:      0.76676
Adj. R-Squared: 0.75311
F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16

也可以利用 R 通过利用 write.csv() 函数直接将数据转换为 .csv 格式的数据储存,再导入 Stata 进行面板数据回归进一步与 R 的回归结果进行对比。使用 R 软件 (panel linear model, plm) 包自带数据集 Produc (Munnell, 1990) 。

*rcall:install.packages("plm",repos ="http://cran.us.r-project.org")
rcall:library(plm) 
rcall:data(Produc,package="plm") 
rcall:aa <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,data = Produc, index = c("state","year")) 
rcall:summary(aa)
Oneway (individual) effect Within Model

Call:
plm(formula = log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, 
    data = Produc, index = c("state", "year"))

Balanced Panel: n = 48, T = 17, N = 816

Residuals:
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-0.120456 -0.023741 -0.002041  0.018144  0.174718 

Coefficients:
             Estimate  Std. Error t-value  Pr(>|t|)    
log(pcap) -0.02614965  0.02900158 -0.9017    0.3675    
log(pc)    0.29200693  0.02511967 11.6246 < 2.2e-16 ***
log(emp)   0.76815947  0.03009174 25.5273 < 2.2e-16 ***
unemp     -0.00529774  0.00098873 -5.3582 1.114e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    18.941
Residual Sum of Squares: 1.1112
R-Squared:      0.94134
Adj. R-Squared: 0.93742
F-statistic: 3064.81 on 4 and 764 DF, p-value: < 2.22e-16
rcall:library(plm) 
rcall:data(Produc,package="plm") 
rcall:write.csv(Produc,file="D:/produc.csv")  //数据储存至D盘
import delimited "D:\produc.csv", encoding(ISO-8859-9) clear 

encode state, gen (id)
gen lgsp = log(gsp)
gen lpcap = log(pcap)
gen lpc = log(pc)
gen lemp = log(emp)

xtset id year
xtreg lgsp lpcap lpc lemp unemp ,fe
       panel variable:  id (strongly balanced)
        time variable:  year, 1970 to 1986
                delta:  1 unit

Fixed-effects (within) regression               Number of obs     =        816
Group variable: id                              Number of groups  =         48

R-sq:                                           Obs per group:
     within  = 0.9413                                         min =         17
     between = 0.9921                                         avg =       17.0
     overall = 0.9910                                         max =         17

                                                F(4,764)          =    3064.81
corr(u_i, Xb)  = 0.0608                         Prob > F          =     0.0000

------------------------------------------------------------------------------
        lgsp |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       lpcap |  -.0261493   .0290016    -0.90   0.368    -.0830815    .0307829
         lpc |   .2920067   .0251197    11.62   0.000     .2426949    .3413185
        lemp |   .7681595   .0300917    25.53   0.000     .7090872    .8272318
       unemp |  -.0052977   .0009887    -5.36   0.000    -.0072387   -.0033568
       _cons |   2.352898   .1748131    13.46   0.000     2.009727    2.696069
-------------+----------------------------------------------------------------
     sigma_u |  .09057293
     sigma_e |  .03813705
         rho |   .8494045   (fraction of variance due to u_i)
------------------------------------------------------------------------------
F test that all u_i=0: F(47, 764) = 75.82                    Prob > F = 0.0000

4.3 拓展:在 R 中运行 Stata

install.packages("RStata") ##在 R 中安装包“RStata”,首次运行需要选择合适镜像地址
library(RStata)  ##加载包
chooseStataBin() ##找出stata具体位置确定打开并复制显示路径
options("RStata.StataPath" =  "\"D:\\Stata16\\StataMP-64\"") 
options("RStata.StataVersion" = 16)  ##设置版本
stata_src <- " 
+ sysuse auto, clear
+ reg mpg weight 
+ "
stata(stata_src)
[1] "\"D:\\Stata16\\StataMP-64\""
. sysuse auto, clear
(1978 Automobile Data)
. reg mpg weight 

      Source |       SS           df       MS      Number of obs   =        74
-------------+----------------------------------   F(1, 72)        =    134.62
       Model |   1591.9902         1   1591.9902   Prob > F        =    0.0000
    Residual |  851.469256        72  11.8259619   R-squared       =    0.6515
-------------+----------------------------------   Adj R-squared   =    0.6467
       Total |  2443.45946        73  33.4720474   Root MSE        =    3.4389

------------------------------------------------------------------------------
         mpg |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      weight |  -.0060087   .0005179   -11.60   0.000    -.0070411   -.0049763
       _cons |   39.44028   1.614003    24.44   0.000     36.22283    42.65774
------------------------------------------------------------------------------

相关课程

连享会-直播课 上线了!
http://lianxh.duanshu.com

免费公开课:


课程一览

支持回看,所有课程可以随时购买观看。

专题 嘉宾 直播/回看视频
最新专题 DSGE, 因果推断, 空间计量等
Stata数据清洗 游万海 直播, 2 小时,已上线
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]

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


关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。直播间 有很多视频课程,可以随时观看。
  • 连享会-主页知乎专栏,300+ 推文,实证分析不再抓狂。
  • 公众号推文分类: 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类,主流方法介绍一目了然:DID, RDD, IV, GMM, FE, Probit 等。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会主页  lianxh.cn
连享会主页 lianxh.cn

连享会小程序:扫一扫,看推文,看视频……

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

✏ 连享会学习群-常见问题解答汇总:
https://gitee.com/arlionn/WD