R语言与Stata等价命令-statar
Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
New！
lianxh命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc,ihelp,rdbalance,gitee,installpkg
⛳ Stata 系列推文：
全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令 计量专题 | 论文写作 | 数据分享 | 专题课程 结果输出 | Stata绘图 | 数据处理 | [Stata程序](https://www.lianxh.cn/b：郭盼亭 (厦门大学 ) logs/26.html) 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法 Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA 文本分析-爬虫 | Python-R-Matlab | 机器学习 Markdown | 工具软件 | 其它
作者：郭盼亭 (厦门大学)
邮箱：gpting2020@163.com
目录
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命令 专题：倍分法DID 专题：Python-R-Matlab 专题：最新课程
相关课程
免费公开课
直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版 Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版 Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页 部分直播课 课程资料下载 (PPT，dofiles等)
最新课程-直播课
| 专题 | 嘉宾 | 直播/回看视频 |
|---|---|---|
| ⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
| 研究设计 | 连玉君 | 我的特斯拉-实证研究设计，-幻灯片- |
| 面板模型 | 连玉君 | 动态面板模型，-幻灯片- |
| 面板模型 | 连玉君 | 直击面板数据模型 [免费公开课，2小时] |
Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。
⛳ 课程主页
⛳ 课程主页
关于我们
Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。 连享会-主页 和 知乎专栏，700+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词： 课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法等
连享会小程序：扫一扫，看推文，看视频……
扫码加入连享会微信群，提问交流更方便
✏ 连享会-常见问题解答：
✨ https://gitee.com/lianxh/Course/wikis
New！
lianxh和songbl命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh
热门资讯
