其他
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日前均可报名！
作者：陈卓然(中山大学)
邮箱：chenzhr25@mail2.sysu.edu.cn
编者按：本文主要参考自 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和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。