R语言与Stata等价命令-statar
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会课程 · 2023 暑期班
作者：郭盼亭 (厦门大学)
邮箱：gpting2020@163.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 前言
2. 前期准备工作
2.1 安装 statar 包
2.2 调用需要的 R 语言包
2.3 导入数据
3. R 函数与 Stata 命令对应关系
3.1 描述性统计数据
3.2 合并数据框
3.3 计算变量的分位数和对数据分组
3.4 对数据进行缩尾
3.5 计时器功能
3.6 判断是否面板数据
3.7 填补缺失值
3.8 画图功能
4. 小结
5. 相关推文
1. 前言
古语有言：工欲善其事，必先利其器。作为一名科研工作人员，熟练掌握并应用一门或者几门统计和数据分析的软件是必不可少的。有些科研工作者可能擅长应用 Stata，有些擅长应用 R 语言，也有些擅长应用 Python 等等。
而当一门统计和分析软件不能够满足科研工作需求的时候，我们就要再花时间和精力去学习和掌握其他的软件。为了助力那些从 Stata 转向 R 语言的使用者们，本文将简单介绍一下 Stata 和 R 语言的联系，即一组 Stata 常用命令的 R 函数。具体来说，这是一组对应于 Stata 常用命令的 R 语言包——statar 包。
关于 R 和 RStudio 的安装教程，详见连享会推文：R和RStudio安装教程。
2. 前期准备工作
2.1 安装 statar 包
类似 Stata 在使用某个命令之前需要 ssc install 命令，我们在 R 语言里使用这组 Stata 常用命令的 R 函数之前，也需要先安装 statar 包，具体操作如下：
install.packages("statar") # 安装 statar 包
install.packages("haven") # 安装应用数据范例分析时可能会应用到的包
install.packages("tidyverse")
2.2 调用需要的 R 语言包
不同于 Stata 在安装了某个命令之后，接下来就可以直接使用，R 语言里还多了一步——调用需要使用的 R 语言包。具体操作如下，调用名称为 statar、haven、tidyverse、magrittr 的四个 R 语言包。
library(statar) # 调用 statar 包
library(haven)
library(tidyverse)
library(magrittr) # 为了使用该版本 R 识别的管道符
2.3 导入数据
范例数据是 Stata 官方自带的 nlswork.dta 数据集。具体操作是使用管道符 (%>%) 把存储路径中的数据赋值给nlswork。(若无法顺利通过以下代码下载数据，也可以通过把网址复制到浏览器，手动下载数据并保存至电脑的某个路径，然后把 url 地址替换成本机电脑路径。)
# 使用管道运算符 (%>%) 从指定的 url 中读取数据框，并将其赋值给 nlswork 变量
nlswork <- "http://www.stata-press.com/data/r17/nlswork.dta" %>%
haven::read_stata()
为了后续分析，这里先分别筛选 nlswork 数据集里边的不同变量，生成 data0、data1、data2 三个数据框。
# 筛选部分变量用于后续分析，并生成三个数据框
data0 <-nlswork %>%
dplyr::select(idcode, year, hours)
data1 <-nlswork %>%
dplyr::select(idcode, year, ln_wage)
data2 <- nlswork %>%
dplyr::select(idcode, year, ln_wage, hours, ind_code)
3. R 函数与 Stata 命令对应关系
3.1 描述性统计数据
一般来说，一篇科研论文的第一个图表是关于文章使用的数据的一个描述性统计表格，其中包含着对各个变量的名称、数据量、均值、标准差等数据特征的统计。在 Stata 中，我们一般通过 sum 命令来实现，对应到 R 语言来，就是 R 语言的 statar 包里的 sum_up。如下示例是对 data2 这个数据框的描述性统计。
首先，对 data2 这个数据框进行描述性统计：
> statar::sum_up(data2)
Variable │ Obs Missing Mean StdDev Min Max
─────────┼───────────────────────────────────────────────────────
hours │ 28467 67 36.5596 9.86962 1 168
idcode │ 28534 0 2601.28 1487.36 1 5159
ind_code │ 28193 341 7.69297 2.99403 1 12
ln_wage │ 28534 0 1.67491 0.47809 0 5.26392
year │ 28534 0 77.9586 6.38388 68 88
另外，还可以对数据框中某一个变量 hours 单独进行描述性统计：
> statar::sum_up(data2, hours)
Variable │ Obs Missing Mean StdDev Min Max
─────────┼───────────────────────────────────────────────────────
hours │ 28467 67 36.5596 9.86962 1 168
除了上述对数据总体情况的一个描述性统计，我们还可以求单个变量的频数、百分比、累积占比。这在 Stata 里常用到的命令是 tabulate，在 R 语言里，可以对应到 statar 包里的 tab() 函数。具体示例如下：
> statar::tab(data2, year) # 对数据框data2中的year变量的描述 (频数、百分比、累积占比)
year │ Freq. Percent Cum.
─────────┼────────────────────────────
68 │ 1375 4.82 4.82
69 │ 1232 4.32 9.14
70 │ 1686 5.91 15.05
71 │ 1851 6.49 21.53
72 │ 1693 5.93 27.47
73 │ 1981 6.94 34.41
75 │ 2141 7.50 41.91
77 │ 2171 7.61 49.52
78 │ 1964 6.88 56.40
80 │ 1847 6.47 62.88
82 │ 2085 7.31 70.18
83 │ 1987 6.96 77.15
85 │ 2085 7.31 84.45
87 │ 2164 7.58 92.04
88 │ 2272 7.96 100.00
3.2 合并数据框
在前文中，本文分别生成了 data0、data1、data2 三个数据框。其中，data0 和 data1 数据框分别含有 idcode、year、hours 和 idcode、year、ln_wage 各三个向量。在 Stata 中，我们可以通过 merge 命令将这两个数据集横向连接，而在 R 语言里，我们可以通过 statar 包里的 join() 函数实现。具体示例如下：
# 合并数据框data0和data1，注意：不同于stata的merge命令，join函数默认是不生成_merge变量的
> statar::join(data1, data0, kind = "full", check = m~1)
# Joining with `by = join_by(idcode, year)` #根据共有的变量idcode,year来匹配
# A tibble: 28,534 x 4
idcode year ln_wage hours
<dbl> <dbl> <dbl> <dbl>
1 1 70 1.45 20
2 1 71 1.03 44
3 1 72 1.59 40
4 1 73 1.78 40
5 1 75 1.78 10
6 1 77 1.78 32
7 1 78 2.49 52
8 1 80 2.55 45
9 1 83 2.42 49
10 1 85 2.61 42
# i 28,524 more rows
# i Use `print(n = ...)` to see more rows
# 由于数据量太大，这里仅展示前面几行的数据。
若要生成 _merge 变量，需要在函数里标注，具体示例如下：
> # merge m:1 v1, gen(_merge)
> statar::join(data0, data1, kind = "full", gen = "_merge")
# Joining with `by = join_by(idcode, year)` #根据共有的变量idcode,year来匹配
# A tibble: 28,534 x 5
idcode year hours ln_wage `_merge`
<dbl> <dbl> <dbl> <dbl> <int>
1 1 70 20 1.45 3
2 1 71 44 1.03 3
3 1 72 40 1.59 3
4 1 73 40 1.78 3
5 1 75 10 1.78 3
6 1 77 32 1.78 3
7 1 78 52 2.49 3
8 1 80 45 2.55 3
9 1 83 49 2.42 3
10 1 85 42 2.61 3
# i 28,524 more rows
# i Use `print(n = ...)` to see more rows
# 由于数据量太大，这里仅展示前面几行的数据。
3.3 计算变量的分位数和对数据分组
statar 包里有一组函数，可以用来计算变量的分位数，或者根据不同的规则对数据进行分组。为了方便展示该函数的应用效果，我们这里不使用上述的数据集，而采取一个包含缺失值 NA 和 1 到 10 的整数的向量来进行举例。
3.3.1 pctile() 函数-计算数据的分位数和加权分位数
首先，给向量 v 赋值为“缺失值 NA 和 1 到 10 的整数”。
> v <- c(NA, 1:10) # v是包含缺失值NA和1到10的整数的向量
其次，在 pctile() 函数中，用 probs = c(0.3, 0.7) 表示指定要计算第 30 和第 70 百分位数，na.rm = TRUE 代表删除缺失值。返回值如下所示，分别为 3.5 和 7.5。
> statar::pctile(v, probs = c(0.3, 0.7), na.rm = TRUE)
30% 70%
3.5 7.5
3.3.2 xtile() 函数-用来对数据分组
xtile() 函数的功能类似于 Stata 的 xtile 命令，可以用来对数据分组。如下是根据百分位点对向量 v 平均分成 3 组，给出的结果是每个数据所属的组别：
> statar::xtile(v, n = 3) # 3 groups based on terciles
[1] NA 1 1 1 1 2 2 2 3 3 3
根据指定的分位数，例子中 probs = c(0.3, 0.7)，将数据分成三组：
> statar::xtile(v, probs = c(0.3, 0.7)) # 3 groups based on two quantiles
[1] NA 1 1 1 2 2 2 2 3 3 3
根据切点 cutpoints = c(2, 3)，将数据分成三组：
> statar::xtile(v, cutpoints = c(2, 3)) # 3 groups based on two cutpoints
[1] NA 1 1 2 3 3 3 3 3 3 3
3.4 对数据进行缩尾
类似 Stata 里的 winsor 命令，R 语言的 statar 包里的 winsorize() 函数也可以用来对数据进行缩尾。例如我们常用的对数据进行上下 1% 的缩尾处理，用 winsorize() 函数可以表示为如下的例子：
# 为展示winsorize()的效果，这里我们将v变量赋值为从1到100的整数。
> v <- c(1:100)
# 对v向量进行上下1%分位数的缩尾处理，返回结果显示如下：
> statar::winsorize(v, probs = c(0.01, 0.99))
0.00 % observations replaced at the bottom
1.00 % observations replaced at the top
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 1
8 19 20 21 22 23 24 25 26 27 28 29 30 31 32
[33] 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
[65] 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
[97] 97 98 99 99
3.5 计时器功能
除了上述的几个功能，statar 包在应对面板数据的时候，还有着计时器、判断是否为面板数据、填补缺失值等功能，这里首先介绍计时器功能。如下示例，首先我们将 "01/04/1992"，"03/15/1992"，"04/03/1992" 这三个时间数据赋值给 date 变量，然后使用 statar 包里的 as.monthly() 函数，把它们定义成一个月份时间变量。
同时，我们将 (4.1、4.5、3.3 ) 赋值给 value 变量。然后，使用 statar 包里的 tlag() 函数，选择滞后一期的 value 变量。
# 由于是需要用到时间序列的数据，所以这里先安装和调用一下timeSeries包和lubridate包,
install.packages("timeSeries")
library(timeSeries)
library(lubridate) # 非常强大,能够识别各种类型的日期 (字符型和时间型数据 )
# 定义一个时间序列变量
> date <- mdy(c("01/04/1992", "03/15/1992", "04/03/1992"))
> datem <- statar::as.monthly(date)
> datem
[1] "1992m1" "1992m3" "1992m4"
# 定义一个value变量
> value <- c(4.1, 4.5, 3.3)
# 应用tlag()函数，以date这个时间变量为尺度， (默认 ) 选择滞后一期的value变量
> statar::tlag(value, 1, time = datem)
[1] NA NA 4.5
3.6 判断是否面板数据
我们首先定义了一个数据框，然后用 is.panel() 函数判断该数据框是否为面板数据。假如根据 id1 来识别，该数据框不是一个面板数据，返回结果为 FALSE，原因是表示时间的 year 变量有缺失值，并且变量 (id1，year) 在第 4，5 两行是重复的。
紧接着，我们删除年份缺失的行的数据，并赋值给 df1，依旧以 id1 为分组依据，判断 df1 是否为面板数据，返回值显示 FALSE，不是面板数据。然后，我们以 (id1，id2) 为分组依据，判断 df1 是否属于面板数据,返回值为TRUE。
由此可见，is.panel() 函数可以帮助我们在开启研究的早期阶段，检验数据质量。
# 定义df的取值
> df <- tibble(
+ id1 = c(1, 1, 1, 2, 2),
+ id2 = 1:5,
+ year = c(1991, 1993, NA, 1992, 1992),
+ value = c(4.1, 4.5, 3.3, 3.2, 5.2)
+ )
# 用statar::is.panel()函数判断df是否为面板数据，返回值为FALSE
> df %>% group_by(id1) %>% statar::is.panel(year)
Variable year has missing values in 1 row(s): 3
Variables (id1 , year) have duplicates for rows (4,5)
[1] FALSE
# 删除年份缺失的行的数据,并赋值给df1
> df1 <- df %>% dplyr::filter(!is.na(year))
# 判断df1是否为面板数据,返回值显示FALSE,不是面板数据
> df1 %>% group_by(id1) %>% statar::is.panel(year)
Variables (id1 , year) have duplicates for rows (3,4)
[1] FALSE
# 以(id1,id2)为分组依据,判断df1是否属于面板数据,返回值为TRUE
> df1 %>% group_by(id1, id2) %>% statar::is.panel(year)
[1] TRUE
3.7 填补缺失值
还有一个在科研的过程中经常遇到的问题，那就是缺失值的问题。我们的解决方案一般有填补缺失值或者删除缺失值，fill_gap() 函数可以在填补面板数据的缺失值的时候助我们一臂之力。具体示例如下：
> df <- tibble(
id = c(1, 1, 1, 2),
datem = statar::as.monthly(mdy(c("04/03/1992", "01/04/1992", "03/15/1992", "05/11/1992"))),
value = c(4.1, 4.5, 3.3, 3.2)
) # 为展示fill_gap()函数的功能，这里先定义一个4 X 3的数据框以便后续举例
> df # 可以看出，这里的datem变量缺少1992m2的取值
# A tibble: 4 x 3
id datem value
<dbl> <monthly> <dbl>
1 1 1992m4 4.1
2 1 1992m1 4.5
3 1 1992m3 3.3
4 2 1992m5 3.2
# 使用fill_gap()函数来填补缺失值
# fill_gap()函数里的第一个参数是要填充缺失值的列 (这里是 datem 列 )
# fill_gap()函数里的第二个参数 roll = "nearest" 表示使用最近的非缺失值来填充缺失值。
> df %>% group_by(id) %>% statar::fill_gap(datem, roll = "nearest")
[1] "datem"
[1] "ok3"
id datem
1: 1 1992m1
2: 1 1992m2
3: 1 1992m3
4: 1 1992m4
5: 2 1992m5
# A tibble: 5 x 3
# Groups: id [2] # 填补之后的结果显示,填补了1992m2的id为1,value(数据值)为4.5
id datem value
<dbl> <monthly> <dbl>
1 1 1992m1 4.5
2 1 1992m2 4.5
3 1 1992m3 3.3
4 1 1992m4 4.1
5 2 1992m5 3.2
3.8 画图功能
stat_binmean() 函数的功能类似 Stata 里用来画分仓散点图的命令——binscatter。该函数将 x 分组后，对于每个组，计算该组中 x 的平均值和 y 的平均值，然后绘制 y 的平均值随 x 的平均值的变化曲线。
接下来介绍如何在 R 语言里应用 statar 包的 stat_binmean() 函数画分仓散点图。具体示例中使用的数据是 R 语言中自带的 iris 数据集。iris 数据集是一个经典的多分类数据集，包含了 3 种不同的鸢尾花 (Setosa、Versicolor、Virginica) 各 50 个样本，每个样本包含 4 个特征：花萼长度 (Sepal Length)、花萼宽度 (Sepal Width)、花瓣长度 (Petal Length)、花瓣宽度 (Petal Width)。
首先，在画图之前，我们还需要调用 ggplot2 包：
> library(ggplot2)
紧接着，我们使用 iris 数据集，应用 ggplot() 函数和 stat_binmean() 函数画了三张图片，每张图片的参数设置见如下具体介绍。
> ggplot(iris, aes(x = Sepal.Width , y = Sepal.Length))+ statar::stat_binmean()
# 这段代码使用R中的ggplot2包创建了一个散点图，
# 显示了鸢尾花数据集中花萼宽度 (Sepal Width) 和花萼长度 (Sepal Length) 之间的关系。
# statar::stat_binmean()默认分组是n=20。
> (ggplot(iris, aes(x = Sepal.Width , y = Sepal.Length, color = Species))
+ statar::stat_binmean(n=10) )
# 类似上一段代码，使用ggplot2包在R中创建散点图，
# 显示了鸢尾花数据集中花萼宽度 (Sepal Width) 和花萼长度 (Sepal Length) 之间的关系。
# 不同的是，这个代码使用了color=Species来根据不同物种 (Species) 对数据点进行着色，
# 将不同物种的数据点区分开来。同时设置stat_binmean()分组为n=10。
> (ggplot(iris, aes(x = Sepal.Width , y = Sepal.Length, color = Species))
+ statar::stat_binmean(n=10) + stat_smooth(method = "lm", se = FALSE))
# 这段代码也是使用ggplot2包在R中创建散点图，
# 显示了鸢尾花数据集中花萼宽度 (Sepal Width) 和花萼长度 (Sepal Length) 之间的关系。
# 与前面的代码相似，这个代码使用了color=Species来根据不同物种 (Species) 对数据点进行着色，
# 将不同物种的数据点区分开来。同时设置stat_binmean()分组为n=10。
# 不同的是，这个代码还使用了stat_smooth(method = "lm", se = FALSE)函数来添加一条线性回归线到图中，
# 用于显示花萼宽度和花萼长度之间的趋势关系。
# 其中method = "lm"表示使用线性回归方法，se = FALSE表示不显示置信区间带。
4. 小结
根据上述介绍，我们可以发现有许多经常在 Stata 中使用的命令都可以在 R 语言的 statar 包里找到对应的函数，上文只是大概介绍了各个函数的应用方法。类似地，Stata 里的每个命令都有许多选项 (option) 设置，在这里每个函数也会有不同的参数设置。对同一个函数设置不同的参数，返回的结果可能就不尽相同，这还需要大家继续探索！
最后，让我们以表格的形式，概括一下 Statar 包各个函数与对应 Stata 常用命令的一一对应关系吧。
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh r语言, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Stata与R语言等价命令 专题：倍分法DID 如何在R语言中实现多期DID 专题：Python-R-Matlab R语言随机指数图分析-statnet R语言绘制社会网络图 专题：最新课程 ⏫ 连享会：R语言初级课程来啦！
课程推荐：2023 暑期班
主讲老师：连玉君，王群勇
🍓 课程主页：https://www.lianxh.cn/news/fdc69c3695aec.html
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。