R与Stata在计量经济学中的比较

发布时间:2022-04-29 阅读 2515

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 系列推文:

PDF下载 - 推文合集

作者:陈卓然(中山大学)
邮箱chenzhr25@mail2.sysu.edu.cn

编者按:本文主要参考自 Paul Goldsmith-Pinkham 教授的博客「Comparing tidyverse R to Stata」,特此致谢!


目录


基础的 R 语言在计量分析当中往往十分笨拙,但是在一些外部包中,这一问题已经部分得到了解决,例如 tidyverse 包。本推文将对 R 和 Stata 进行比较,探究二者在计量分析中各自存在的优劣势。

1. 数据集简介

本推文使用的数据集为 nlswork.dta,它是对美国在 1968 年时 14 到 24 岁年轻女性的 21 年调查数据。该数据集共有 4711 位女性,在这段时间内她们都处于工作状态,并且小时工资在 1-700 美元之间。我们研究的是在所有已婚人士中,相比于其他人种,白种人是否会更加努力勤奋工作。

2. 数据清洗与估计

首先用 Stata 进行估计,具体命令如下:

. * load Data + create variables
. lxhuse nlswork, clear
. keep if msp == 1 // 只保留已婚人士
. qui sum wks_work
. gen HardWorker = wks_work > r(mean) if ~missing(wks_work) 
. //wks_work指去年总共工作的周数,生成勤奋工作者的虚拟变量
. gen white = race == 1 if ~missing(race) // 生成白种人
. keep HardWorker white idcode year grade

. * Basic summary stats by groups
. mean white
. tab year, sum(white)             // 按照年份进行分组统计
. tab idcode in 1/100, sum(white)  // 按照个体进行分组统计

. * Simple Panel Regression
. areg HardWorker    white, absorb(year) cluster(idcode)
. * 或者采用下述命令
. reghdfe HardWorker white, absorb(year) cluster(idcode)  

HDFE Linear regression                            Number of obs   =     16,763
Absorbing 1 HDFE group                            F(   1,   3612) =       2.69
Statistics robust to heteroskedasticity           Prob > F        =     0.1012
                                                  R-squared       =     0.5468
                                                  Adj R-squared   =     0.5464
                                                  Within R-sq.    =     0.0002
Number of clusters (idcode)  =      3,613         Root MSE        =     0.3272
                             (Std. err. adjusted for 3,613 clusters in idcode)
------------------------------------------------------------------------------
             |               Robust
  HardWorker | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       white |     -0.011      0.007    -1.64   0.101       -0.024       0.002
       _cons |      0.390      0.006    66.59   0.000        0.379       0.402
------------------------------------------------------------------------------

接下来我们使用 R 进行同样的操作。在展示最终的结果之前,有必要了解一下在 R 中如何进行高维固定效应的估计。我们使用的命令是 felm, 这一命令是内嵌在 lfe 包中的,因此在使用之前需要先加载 lfe,即 library(lfe)。这一命令的语法格式是:

felm(
  formula,
  data,
  exactDOF = FALSE,
  subset,
  na.action,
  contrasts = NULL,
  weights = NULL,
  ...
)

这里主要介绍 formula 部分,其余部分使用不多,此处暂且略过。formula 这一部分是由四个子部分构成的,第一部分是普通的协变量,也就是控制变量;第二部分是需要控制的固定效应;第三部分是 IV 的设定;第四部分是标准误的聚类设定。也就是类似如下形式:

y ~ x1 + x2 | f1 + f2 | (Q|W ~ x3 + x4) | clu1 + clu2

其中,y 是被解释变量,x1x2 是普通的协变量,f1f2 是需要控制的固定效应,QW 是内生变量,它们的工具变量是 x3x4clu1clu2 是用来计算聚类标准误的变量。上述四个部分中,如果某一部分不存在,就用 0 代替,除非是最后一部分 (0 可以省略)。下面是使用 R 进行上述相同操作的代码:

library("readxl")
library("dplyr")
library("tidyr")
library("broom")
library("lfe")
nlswork <- read_excel(path = "D:\\Stata17\\personal\\连享会推文\\推文blog\\推文6-reghdfe\\nlswork.xlsx")
## 保留已婚人士
## 生成勤奋工作者
## 生成白种人

reghdfeData <- nlswork %>% filter(msp == 1) %>%
  mutate(Hardworkers = wks_work > mean(wks_work, na.rm = TRUE),
         White = race == "White") %>%
  select(idcode, year, Hardworkers, White)
summary(reghdfeData)

## Basic Summary Statistics (仅能给出均值,但是在Stata中可以给出
## 标准误以及置信区间)
reghdfeData %>% summarise(mean(White)) 
## 分组回归的结果 (Stata 中不能将mean 与 bys连用,因此只能使用sum)
### 按年份分组
reghdfeData %>% group_by(year) %>% summarise(mean(White))
### 按个体分组
reghdfeData %>% group_by(idcode) %>% summarise(mean(White))
## 回归:探究是否为白人对是否努力工作的影响 (高维固定效应,类似于
## Stata中的reghdfe)
est <- felm(Hardworkers ~ White | year | 0 | idcode,data = reghdfeData)
tidy(est)
summary(est)

这是一个非常简单的小例子,我们并不需要做很多数据清理的工作,看上去 Stata 和 R 并没有显著区别。但是当我们需要做很多数据清洗工作时,尤其是涉及到合并数据集,我们最好还是选择 R。R 中可以创建不同的矩阵或者向量存储数据,因而可以很方便的纵向合并或者横向合并数据。在 Stata 当中,当涉及到合并数据的操作时,我们不得不反复加载和存储数据到内存中,显得很麻烦。

从结果来看,尽管二者完全一样,但是很明显,Stata 中的结果输出地更加友好,相比之下,R 输出的结果略显凌乱,因此在结果输出方面,Stata 要略胜一筹。

3. 加入工具变量

我们不妨假设教育背景与是否勤奋工作无关,而与是否为白种人相关。白种人往往可以获得比其他人种更多的学习机会,因而拥有更高的教育背景,而是否勤奋工作与教育背景不相关。注意此处仅作为演示,上述工具变量相关性和外生性的论述可能不太恰当。具体地,使用 Stata 估计的命令如下:

. * 加入工具变量
. ivreghdfe HardWorker (white = grade), absorb(year) cluster(idcode)

IV (2SLS) estimation
------------------------------------------------------------------------------
             |               Robust
  HardWorker | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       white |      0.348      0.077     4.50   0.000        0.196       0.499
------------------------------------------------------------------------------
Underidentification test (Kleibergen-Paap rk LM statistic):             46.819
                                                   Chi-sq(1) P-val =    0.0000
------------------------------------------------------------------------------
Weak identification test (Cragg-Donald Wald F statistic):              336.528
                         (Kleibergen-Paap rk Wald F statistic):         52.745
Stock-Yogo weak ID test critical values: 10% maximal IV size             16.38
                                         15% maximal IV size              8.96
                                         20% maximal IV size              6.66
                                         25% maximal IV size              5.53
Source: Stock-Yogo (2005).  Reproduced by permission.
NB: Critical values are for Cragg-Donald F statistic and i.i.d. errors.
------------------------------------------------------------------------------
Hansen J statistic (overidentification test of all instruments):         0.000
                                                 (equation exactly identified)
------------------------------------------------------------------------------
Instrumented:         white
Excluded instruments: grade
Partialled-out:       _cons
                      nb: total SS, model F and R2s are after partialling-out;
                          any small-sample adjustments include partialled-out
                          variables in regressor count K
------------------------------------------------------------------------------

使用 R 估计的命令如下:

## 加入工具变量之后进行回归
est2 <- felm(Hardworkers ~ 0 | year | (White ~ grade) | idcode, data = reghdfeData)
tidy(est2)
summary(est2)

4. 其他方面的比较

下面我将简要介绍一些常用的 R 和 Stata 处理相同任务时的命令,以方便小伙伴们在两种语言之间快速切换。

Stata R Description
clear all rm(list = ls()) Clears data, value labels, etc from memory
insheet using "foo.csv", comma names mydata <- read.csv("foo.csv", header = TRUE) Change working directories
pwd getwd() Display the working directory
reg y x1 x2 summary(lm(y ~ x1 + x2, data = mydata)) Ordinary least squares with constant
reg y x1 x2, nocon summary(lm(y~x1+x2-1, data = mydata)) Ordinary least squares without constant
reg y x if x>0 lm(y~x,data = subset(mydata, x>0)) Select a conditional subset of data
forvalues i = 1/100 {...} for (i in 1:100) {...} Loop through integer values of i from 1 to 100
foreach i of varlist "a b c" {...} for (i in c("a","b","c")) {...} Loop through a list of items
logit y x summary(glm(y~x,data = mydata, family = "binomial")) Perform logit maximum likelihood estimatio
probit y x summary(glm(y~x, data = mydata, family = binomial(link = "probit"))) Perform probit maximum likelihood estimation
sort x y mydata[order(mydata\$x, mydata\$y), ] Sort the data frame by variable x
summarize summary(mydata) Provide summary values for data
replace x = y1 + y2 mydata\$x<-with(mydata, y1+y2) Change the x value of data to be equal to y1+y2
drop x mydata\$x<-NULL Drop variable x from the data
append using "mydata2.dta" mydata<-rbind(mydata,mydata2) Append mydata2 to mydata
merge 1:1 index using "mydata2.dta" merge(mydata, mydata2, index) Merge two data sets together by index variable(s)
set obs 1000 /// gen x = rnormal() mydata\$x<-rnorm(1000) Generate 1000 random normal draws
count nrow(mydata) Count the number of observations in the data
rename oldvar newvar colnames(dataframe)[colnames(dataframe) == "oldvar"]<-"newvar" Rename Variables
egen id = group(x y) mydata\$ID <- with(mydata, ave(ID, list(x,y), FUN = seq_along)) Create an identifier ID from variables x and y
egen smallist = min(x y) mydata %>% rowwise() %>% mutate(smallest = min(c(x,y))) For each row, find the smallest value within a set of columns
egne newvar = fcn(x y) mydata %>% rowwise() %>% mutate(newvar = fcn(c(x,y))) For each row, apply a function to variables

5. 参考文献

  • Dictionary: Stata to R -Link-
  • Comparing tidyverse R to Stata -Link-
  • Muenchen R A, Hilbe J. R for Stata users[M]. New York, NY: Springer, 2010. -PDF-
  • Stata & R Command Dictionary. -PDF-
  • Getting Started in R-Stata Notes on Exploring Data. -PDF-

6. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 交互, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

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

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

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

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

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