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

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

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

连享会 · 2022 面板数据因果推断专题 

👉 10月30日前均可报名！

作者：陈卓然(中山大学)
邮箱：[e-mail]

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

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

---

目录

• 1. 数据集简介

• 2. 数据清洗与估计

• 3. 加入工具变量

• 4. 其他方面的比较

• 5. 参考文献

• 6. 相关推文

---

基础的 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 是被解释变量，x1 和 x2 是普通的协变量，f1 和 f2 是需要控制的固定效应，Q 和 W 是内生变量，它们的工具变量是 x3和 x4，clu1 和 clu2 是用来计算聚类标准误的变量。上述四个部分中，如果某一部分不存在，就用 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

• 专题：计量专题
• 主成分分析-交互固定效应基础：协方差矩阵的几何意义
• 专题：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：二者配合的基本设定
• 专题：Stata命令
• Stata-Mata系列 (二)：Mata与Stata的交互
• 专题：面板数据
• regife：面板交互固定效应模型-Interactive Fixed Effect
• 专题：交乘项-调节
• med4way：中介效应和交互效应分析
• interactplot：图示交乘项-交互项-调节效应
• Stata：图示交互效应-调节效应
• 专题：Probit-Logit
• Logit-Probit：非线性模型中交互项的边际效应解读
• 专题：Python-R-Matlab
• Stata交互：Python-与-Stata-对比

课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp

 

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

🍏 关于我们

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