New！ 搜推文，找资料，用
lianxh命令：
安装：ssc install lianxh, replace
使用：lianxh 合成控制
lianxh DID + 多期, w
作者：刘晓飞 (四川大学)
邮箱：liuxiaofei@stu.scu.edu.cn
- Title: 基本的数据处理和统计分析实例
- Keywords: R，数据处理，正则表达式，管道操作
编者按：本文整理自以下内容，特此致谢！
Source：Davies Ben, 2021, Blog, Female representation and collaboration at the NBER. -Link- -R-
1. 背景介绍
在实证研究中，数据处理和统计分析是理解和解释数据的关键步骤。这些步骤通常包括数据清洗、数据转换和数据分析。在这一过程中，正则表达式和管道操作起着至关重要的作用。本推文旨在介绍一些基础的管道操作和正则表达式操作，并通过文章详细展示如何在 R 中进行基本的数据处理和统计分析。
2. R 实例
2.1 管道操作命令
在 R 语言中，管道操作符 (pipe operator) 是一种非常方便的语法，允许你将前一个函数的输出直接传递给下一个函数，这使得代码更加清晰易读。R 语言中最著名的管道操作符是 %>%，在 dplyr 等包中被广泛使用。
2.1.1 基本用法
管道操作符 %>% 的基本语法如下：
input %>% function1() %>% function2() %>% function3()这里的意思是，input 的输出将成为 function1() 的输入，function1() 的输出又传递给 function2()，依此类推。让我们通过几个具体的示例，来看看如何在 R 中使用管道操作处理 mtcars 数据集：
> library(dplyr)
> # 查看数据的前几行
> mtcars %>%
+ head()
mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1> mtcars %>%
+ group_by(cyl) %>% # 根据cyl列分类
+ summarise(
+ mean_hp = mean(hp), # 对hp求均值
+ mean_disp = mean(disp) # 对disp求均值
+ )
# A tibble: 3 × 3
cyl mean_hp mean_disp
<dbl> <dbl> <dbl>
1 4 82.6 105.
2 6 122. 183.
3 8 209. 353.2.1.2 进阶用法
在 R 的管道操作中，通常情况下，每个函数都会自动接收前一步的输出作为其第一个参数。但有时，我们需要在管道操作中进行一些不那么直接的操作，这时就可以用花括号 {} 和点符号 . 来手动控制输入。
result <- data %>%
filter(column1 > 10) %>%
{subset(., select = c(column1, column2))} # 使用{}和.来处理更复杂的操作- 花括号
{}：它允许你在管道中执行复杂的表达式或多步操作。在花括号内，你可以写几乎任何有效的 R 代码。 - 点符号
.：这是管道操作中的占位符，代表前一步的结果。在这个例子中，.代表了filter(column1 > 10)的输出。
让我们通过下边这个示例，来看看如何在 R 中使用 {} 和 . 来处理 mtcars 数据集：
> mtcars %>%
+ {
+ . <- mutate(., efficiency = ifelse(mpg > 25, "High", ifelse(mpg > 20, "Medium", "Low")))
+ . <- head(.) # 返回修改后的数据集的前几行
+ .
+ }
mpg cyl disp hp drat wt qsec vs am gear carb efficiency
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 Medium
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 Medium
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 Medium
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 Medium
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 Low
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1 Low在这个例子中，ifelse() 函数用来基于 mpg 值设置 efficiency 列的值。我们使用 . 来表示 mtcars 数据集，并将其传递给 mutate() 函数。
2.2 正则表达式介绍
R 语言中的正则表达式是一种强大的工具，用于文本处理、数据清洗和模式识别。正则表达式通过特定的模式来匹配字符串中的字符组合。在 R 中，你可以使用多个函数来实现正则表达式的匹配和操作，这些函数主要包括 grep()、grepl()、gsub()、sub()、regexpr() 和 gregexpr()。
2.2.1 基本匹配
grep()：返回输入向量中与模式匹配元素的索引。grepl()：返回一个逻辑向量，指示输入向量中的每个元素是否与模式匹配。
> # 找出包含"mpg"的行的索引
> grep("mpg", names(mtcars))
[1] 1
> # 检查向量中的元素是否包含数字
> grepl("\\d", c("abc", "def1", "ghi2"))
[1] FALSE TRUE TRUE2.2.2 替换修改
sub() 和 gsub()：sub() 替换字符串中的第一个匹配项，gsub() 替换所有匹配项。
> # 替换字符串中的第一个数字
> sub("\\d", "X", "abc1def2")
[1] "abcXdef2"
> # 替换字符串中的所有数字
> gsub("\\d", "X", "abc1def2")
[1] "abcXdefX"2.2.3 正则表达式的语法
在正则表达式中，有几个常用的特殊字符和模式：
.：匹配任意单个字符。^：匹配字符串的开始。$：匹配字符串的结束。[...]：匹配方括号内的任何一个字符。[^...]：匹配不在方括号内的任何一个字符。|：匹配两个子模式中的一个 (或)。\*：匹配前面的模式零次或多次。+：匹配前面的模式一次或多次。?：匹配前面的模式零次或一次。{n}：匹配前面的模式恰好 n 次。{n,}：匹配前面的模式至少 n 次。{n,m}：匹配前面的模式至少 n 次，但不超过 m 次。
2.2.4 使用案例
假设你有一组文本数据，你想提取其中的日期信息 (格式如 dd/mm/yyyy)：
> texts <- c("This event happened on 12/05/2021 and it was significant.")
> dates <- regmatches(texts, gregexpr("\\b\\d{2}/\\d{2}/\\d{4}\\b", texts))
> dates
[[1]]
[1] "12/05/2021"2.3 案例复现
接下来，我们将通过复现《Female representation and collaboration at the NBER》一文来探讨如何运用 R 语言进行基础的数据处理与统计分析。这篇文章分析了过去四十年来 NBER 工作论文的女性作者的代表性以及她们之间的合作。
2.3.1 逐年女性作者发表论文比例
首先，我们将着手处理和分析 papers 与 authors 这两个数据集。在这一过程中，我们的焦点放在对作者名字进行性别估计的技术上，并据此对数据进行深入的统计分析。以下是这一过程的详细步骤分解：
步骤1: 数据预处理
- 转换年份为十年段：使用
mutate()和floor()函数将year列中的每个年份转换为对应的十年段，例如 1987 年转换为 1980 年。 - 筛选和选择列：选择
paper和新计算的decade列，并过滤出那些其十年段存在于decades向量中的记录。 - 数据连接：使用
left_join()将papers数据与authors数据框进行左连接，这假定两个数据集中都有可以匹配的键。 - 名字提取：通过
mutate()和case_when()添加一个新列name，根据不同的正则表达式模式提取作者的名字。 - 名字存在性标记：添加
has_first_name列，标记是否成功从作者字段中提取到名字。
步骤2: 性别估计与数据整合
- 筛选有名字的记录：过滤出那些成功提取出名字的记录。
- 性别估计：使用
gender函数对名字进行性别估计，考虑 1940 至 1995 年的数据，并将估计结果作为新的左连接合并到数据中。 - 选择和重命名列：选择相关列，并将
proportion_female列重命名为p_female。 - 数据去重和去 NA：使用
distinct()去除重复记录，使用drop_na()去除含有 NA 的行。 - 重新连接原始数据：将处理过的数据与原始
tmp数据进行左连接。
步骤3: 数据分析和报告
- 数据分组和变量计算：按十年段分组，并计算每个十年段中的论文数和作者数。
- 去除不需要的列和去重：去除
paper列并再次去除重复的记录。 - 统计计算：分组计算每个十年段的论文数和作者数，计算女性作者的百分比。
- 性别可估计率计算：计算可估计性别的作者百分比。
- 数据展示：使用
kable()函数输出格式化的表格，用于报告或展示。
# 需要安装bldr和nberwp(version=0.1.0) package
# install.packages("remotes")
# remotes::install_github("bldavies/bldr")
library(bldr)
library(dplyr)
library(gender)
library(ggplot2)
library(igraph)
library(knitr)
library(nberwp)
library(purrr)
library(tidygraph)
library(tidyr)
decades <- c(1980, 1990, 2000, 2010) # 定义您想要筛选的十年段
tmp = papers %>%
# 将年份转换为十年段，例如1987年将转换为1980年
mutate(decade = floor(year / 10) * 10) %>%
# 选择 'paper' (论文标题) 和新计算的 'decade' 列
select(paper, decade) %>%
# 过滤出那些其十年段存在于 'decades' 向量中的记录
filter(decade %in% decades) %>%
# 与 'authors' 数据框进行左连接，假定两个数据框中都有匹配的键 (通常是作者名或ID)
left_join(authors) %>%
# 添加一个新的列 'name'，基于作者名的格式提取名字，并检查是否成功提取出名字
mutate(
name = case_when(
# 如果作者名以两个任意字母开头，提取第一个以空格分隔的单词作为名字
grepl('^[A-Za-z]{2}', author) ~ sub('^([A-Za-z-]+) .*', '\\1', author),
# 如果作者名以大写字母开头，后跟空格和其他字符，提取第一个单词作为名字
grepl('^[A-Z] .* .*', author) ~ sub('^[A-Z] ([A-Za-z]+) .*', '\\1', author),
# 如果上述条件都不符合，返回空字符串作为名字
TRUE ~ ''
),
# 添加 'has_first_name' 列，如果 'name' 列非空，则为TRUE，否则为FALSE
has_first_name = name != ''
)
df = tmp %>%
# 筛选出那些成功提取出名字的记录
filter(has_first_name) %>%
# 使用gender包的gender函数来估计名字对应的性别比例，考虑1940至1995年的数据
left_join(gender(.$name, years = c(1940, 1995))) %>%
# 选择相关的列并重命名proportion_female为p_female
select(paper, decade, author, p_female = proportion_female) %>%
# 删除重复记录
distinct() %>%
# 删除包含NA的行
drop_na() %>%
# 将原始的tmp数据与上述处理过的数据进行左连接
{left_join(tmp, .)}
df %>%
# 按十年段分组并添加's'后缀
group_by(Decade = paste0(decade, 's')) %>%
# 计算每个十年段中的论文和作者数量
mutate(Papers = n_distinct(paper),
Authors = n_distinct(author)) %>%
# 删除含有NA的行
drop_na() %>%
# 去除paper列
select(-paper) %>%
# 删除重复记录
distinct() %>%
# 再次按Decade, Papers, Authors分组
group_by(Decade, Papers, Authors) %>%
# 对当前分组计算数量并计算女性作者的百分比
summarise(tmp = n(),
`% authors female` = 100 * mean(round(p_female))) %>%
# 解除分组
ungroup() %>%
# 计算可估计性别的作者的百分比
mutate(`% authors with estimable sex` = 100 * tmp / Authors) %>%
# 选择显示的列
select(c(1:3, 5:6)) %>%
# 输出为一个表格，设置对齐、数字格式等
kable(align = 'c', digits = 1, format.args = list(big.mark = ','))最终统计结果如下表所示：
| Decade | Papers | Authors | % authors female | % authors with estimable sex |
|:------:|:------:|:-------:|:----------------:|:----------------------------:|
| 1980s | 2,820 | 972 | 14.1 | 93.9 |
| 1990s | 4,213 | 2,211 | 19.7 | 88.3 |
| 2000s | 8,188 | 5,118 | 24.0 | 85.5 |
| 2010s | 10,970 | 9,519 | 27.0 | 84.1 |上表列出了 1980 年代、1990 年代、2000 年代和 2010 年代 NBER 工作论文和作者的数量。它还报告了估计为女性作者的百分比，以及可以估计性别的作者百分比。作者数量大约每十年翻一番，女性作者比例在从 1980 年代到 2010 年代几乎翻了一番。
2.3.2 跨研究项目的代表性
步骤1: 项目合并统计
- 数据合并：
left_join(programs)将先前的df数据框与programs数据框合并，添加相关的项目类型信息。 - 去重：
distinct(program, paper)选择唯一的项目和论文对，以避免重复计算同一篇论文对应多个项目。 - 计数添加：
add_count(paper)为每篇论文添加一个计数，统计其出现的频率。 - 加权计数：
count(program, wt = 1 / n, sort = T)计算每个项目的加权计数，权重为论文频率的倒数，sort = T表示按照计数排序。 - 项目调整：
mutate(program_adj = ifelse(dense_rank(-n) <= 10, program, 'Other'))根据频率对项目进行排名，只保留前 10 个最频繁的项目，其余归类为 Other。
步骤2: 统计调整后的项目计数
- 计数：根据调整后的项目名
program_adj计数，权重为上一步计算的加权计数n。 - 排序：结果首先确保 All 项在最后显示，其余按照计数逆序排序。
步骤3: 数据分析和报告
- 数据合并：再次将
df与programs和program_counts合并。 - 添加全总计行：使用
bind_rows()添加一个表示全部项目的行，标记为 "All"。 - 去重并调整因子级别：去除重复行，并将
program_adj转换为因子，级别为先前统计的项目加 "All"。 - 分组和汇总：按十年段和项目类型分组，计算女性作者百分比平均值并格式化。
- 数据重排：使用
spread()将数据从长格式转为宽格式，不同的十年作为列展示。 - 条件修改：对于 "Asset Pricing" 在 1980s 的数据 (只有一篇论文)，显示为 "-"。
- 添加描述和选择列：结合项目类型描述，并选择显示的列。
- 表格输出：使用
kable()生成 Markdown 或 HTML 表格。
# 加载dplyr包用于数据处理和tidyr包用于数据整理
library(dplyr)
library(tidyr)
# 处理原始数据，计算每个项目的加权计数，并调整罕见项目为'Other'
program_counts = df %>%
left_join(programs) %>% # 左连接programs数据框，添加项目信息
distinct(program, paper) %>% # 移除重复的项目和论文对
add_count(paper) %>% # 对每篇论文进行计数
# 对项目进行加权计数，权重为1/论文数量
count(program, wt = 1 / n, sort = T) %>%
# 如果项目不在前10，则归类为'Other'
mutate(program_adj = ifelse(dense_rank(-n) <= 10, program, 'Other'))
# 对调整后的项目进行再计数，并按特定顺序排列
program_adj_counts = program_counts %>%
# 对调整后的项目名进行计数
count(program_adj, wt = n) %>%
# 确保'Other'分类在末尾显示，其余按计数逆序排列
arrange(program_adj == 'Other', -n)
# 加载kableExtra包用于美化表格输出
library(kableExtra)
# 合并数据并计算每个项目每个十年的女性作者百分比
df %>%
left_join(programs) %>% # 左连接programs数据框
left_join(program_counts) %>% # 左连接项目计数数据
{bind_rows(., mutate(., program_adj = 'All'))} %>% # 向数据框添加总计行'All'
distinct(program_adj, decade, author, p_female) %>% # 移除重复行
# 将program_adj转为因子
mutate(program_adj = factor(program_adj, c(program_adj_counts$program_adj, 'All'))) %>%
# 按十年和项目分组
group_by(decade = paste0(decade, 's'), program = program_adj) %>%
# 计算并格式化女性作者的平均百分比
summarise(p_female = sprintf('%.1f', 100 * mean(round(p_female), na.rm = T))) %>%
# 将数据重塑为宽格式，填充缺失值为'-'
spread(decade, p_female, fill = '-') %>%
# 对特定情况进行调整
mutate(`1980s` = ifelse(program == 'CF', '-', `1980s`)) %>%
# 左连接项目描述数据框
left_join(program_descriptions) %>%
# 根据描述修改项目名称
mutate(Program = ifelse(!is.na(program_desc), paste0(program_desc, ' (', program, ')'), program)) %>%
# 选择需要的列
select(Program, 2:5) %>%
# 使用kable生成表格并进行格式化
kable(align = 'lrrrr', format = "html", escape = FALSE) %>%
# 添加kableExtra样式选项
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))最终结果如下表所示：
|Program | 1980s| 1990s| 2000s| 2010s|
|:---------------------------------------------------|-----:|-----:|-----:|-----:|
|Labor Studies (LS) | 19.1| 26.2| 27.1| 29.9|
|Economic Fluctuations and Growth (EFG) | 5.9| 9.2| 17.4| 18.9|
|Public Economics (PE) | 8.5| 16.7| 21.9| 26.3|
|International Finance and Macroeconomics (IFM) | 14.5| 13.8| 16.5| 20.7|
|International Trade and Investment (ITI) | 14.7| 15.9| 23.6| 23.2|
|Monetary Economics (ME) | 5.3| 11.8| 13.9| 17.5|
|Asset Pricing (AP) | -| 10.7| 16.5| 18.1|
|Productivity, Innovation, and Entrepreneurship (PR) | 15.4| 23.4| 22.2| 24.1|
|Corporate Finance (CF) | -| 12.9| 22.4| 20.6|
|Health Economics (HE) | 20.4| 23.3| 33.9| 33.5|
|Other | 11.2| 22.7| 26.4| 28.4|
|All | 14.1| 19.7| 24.0| 27.0|分析女性代表性的另一种方法是比较各项目中女性撰写的工作论文的密度。我在下面的图表中展示了这一对比，重点是 2010 年代发表的论文。横轴表示每个项目中女性作者发表的工作论文的百分比。
plot_df = df %>%
# 过滤出数据集中“decades”列最大值对应的行，即最近的十年数据
filter(decade == max(decades)) %>%
# 创建新列“female”，标识女性作者的比例大于50%
mutate(female = p_female > 0.5) %>%
# 计算每篇论文的作者数量，结果存储在新列 n_authors 中
add_count(paper, name = 'n_authors') %>%
# 移除“female”列中有缺失值的行
filter(!is.na(female)) %>%
# 选择特定的列进行进一步分析
select(paper, n_authors, author, female) %>%
# 左连接“programs”数据集，以获取每个作者所属的研究项目信息
left_join(programs) %>%
# 在数据集中添加一个总体项目“All”，用于后续全体数据的比较
{bind_rows(., mutate(., program = 'All'))} %>%
distinct() %>% # 删除重复行
# 移除 program 列中有缺失值的行
filter(!is.na(program)) %>%
# 计算每篇论文和每位作者的研究项目数量
add_count(paper, author, name = 'n_programs') %>%
# 计算权重，为每位作者的论文数和所参与的项目数的倒数
mutate(wt = 1 / (n_authors * n_programs)) %>%
# 根据性别和项目对论文进行加权计数
count(female, program, wt = wt) %>%
group_by(program) %>% # 按照项目分组
# 计算每个项目中女性作者论文的百分比
summarise(p = 100 * sum(female * n) / sum(n)) %>%
ungroup() %>% # 取消分组
# 左连接 program_descriptions 数据集，以获取项目的描述信息
left_join(program_descriptions)
plot_df %>%
# 选取不包括“全部项目”的数据
filter(program != 'All') %>%
# 对女性作者的百分比进行排序，并生成新列pos
mutate(pos = dense_rank(p)) %>%
# 使用ggplot2包绘制图形，横坐标为排名，纵坐标为百分比
ggplot(aes(pos, p)) +
# 添加一条虚线，表示所有项目的平均值
geom_hline(yintercept = filter(plot_df, program == 'All')$p, lty = 'dashed') +
# 添加柱状图，颜色填充表示项目类别
geom_col(aes(fill = program_category), alpha = 0.33) +
# 在每个柱状图中心添加项目描述文本
geom_text(aes(y = p / 2, label = program_desc), size = 3) +
coord_flip(clip = 'off') + # 翻转坐标轴，使柱状图水平展示
labs(x = NULL, # 不显示x轴标签
y = '% of (fractionally counted) papers by female authors', # 设置y轴标签
# 设置图形标题
title = 'NBER working papers authored by females during the 2010s,\nby research program',
# 设置图形副标题
subtitle = 'Dashed line represents female authorships across all programs',
fill = 'Category') + # 设置图例标题
# 调整图例位置和对齐方式
guides(fill = guide_legend(label.position = 'left', title.hjust = 1)) +
# 设置柱状图颜色
scale_fill_discrete(breaks = sort(unique(program_descriptions$program_category))) +
scale_x_discrete(expand = c(0, 0)) + # 调整x轴展开设置
scale_y_continuous(expand = c(0, 0)) + # 调整y轴展开设置
theme(legend.justification = c(1, 0), # 调整图例的布局
legend.position = c(1, 0)) # 设置图例的位置
# 保存图形为PNG文件，设置大小和分辨率
ggsave('_media/female-authorships.png', width = 6, height = 6, dpi = 300)最终结果如下图。总体而言，在 2010 年代发表的工作论文中，女性约占 21%。这些论文相对集中在应用微观经济学的项目中，而不是宏观经济学或金融学。
2.3.3 共同创作模式
接下来，从每十年的工作论文合作网络中推断出 NBER 作者之间的合作模式。在每个网络中，节点对应于在该十年内至少发表一篇工作论文的作者，而边则连接在该十年内共同撰写了至少一篇工作论文的作者。
# 节点数据准备
nodes_df = df %>%
# 左连接 programs 数据框，根据相同列合并数据
left_join(programs) %>%
# 在每个数据组中加入一个总体类别 All
{bind_rows(., mutate(., program = 'All'))} %>%
# 过滤掉项目信息缺失的行
filter(!is.na(program)) %>%
# 去除 paper 列
select(-paper) %>%
distinct() %>% # 删除重复的行
# 按照 decade 和 program 进行分组
group_by(decade, program) %>%
# 将每个组的数据嵌套在一个新列 data 中
nest() %>%
ungroup() %>% # 取消分组
rename(nodes = data) # 将 data 列重命名为 nodes
# 边缘数据准备
edges_df = df %>%
# 移除 p_female 列
select(-p_female) %>%
# 数据自身连接，根据 paper 和 decade 列匹配
{left_join(., ., by = c('paper', 'decade'))} %>%
# 仅保留 author.x 小于 author.y 的行，避免重复边
filter(author.x < author.y) %>%
# 左连接 programs 数据框，添加项目信息
left_join(programs) %>%
# 加入一个总体项目类别 All
{bind_rows(., mutate(., program = 'All'))} %>%
# 过滤掉项目信息缺失的行
filter(!is.na(program)) %>%
# 删除重复的行
distinct(decade, program, author.x, author.y) %>%
# 按照 decade 和 program 进行分组
group_by(decade, program) %>%
# 将每个组的数据嵌套在一个新列 data 中
nest() %>%
# 取消分组
ungroup() %>%
# 将 data 列重命名为 edges
rename(edges = data)
# 网络构建函数
func = function(x, y) {
x %>%
# 使用数据框创建无向图
graph_from_data_frame(directed = F, vertices = y$author) %>%
# 转换为 tbl_graph 对象
as_tbl_graph() %>%
# 将顶点数据与原始数据左连接
left_join(y, by = c('name' = 'author'))
}
# 网络数据集创建与分析
nets_df = nodes_df %>%
# 内连接节点和边缘数据框
inner_join(edges_df) %>%
# 对每对边缘和节点数据应用网络构建函数
mutate(net = pmap(list(edges, nodes), func),
# 在 decade 列数据后添加's'
decade = paste0(decade, 's'))
get_network_properties = function(G) {
tibble(
# 计算图中的顶点数
Nodes = gorder(G),
# 计算图中的边数
Edges = gsize(G),
# 计算图的边密度
`Edge density (%)` = 100 * graph.density(G),
# 计算图的平均度数
`Mean degree` = mean(degree(G))
)
}
nets_df %>%
# 选取包含所有项目的数据行
filter(program == 'All') %>%
# 对每个网络计算网络属性
mutate(res = map(net, get_network_properties)) %>%
# 选择 decade 和结果列
select(Decade = decade, res) %>%
# 展开 res 列中的嵌套数据
unnest('res') %>%
# 使用kable格式化表格输出结果
kable(align = 'c', digits = 2, format.args = list(big.mark = ','))下表总结了这个合作网络的特点。
网络特征的变化：随着时间的推移，这些合作网络变得更大但密度变低。密度低意味着虽然网络中的作者数量增加，但每个作者的直接联系 (即直接合作关系) 相对减少。
平均度数的上升：平均度数，即平均合作者数量的上升，反映了在其他研究中记录的经济学家之间合作增加的趋势。
| Decade | Nodes | Edges | Edge density (%) | Mean degree |
|:------:|:-----:|:------:|:----------------:|:-----------:|
| 1980s | 972 | 1,197 | 0.25 | 2.46 |
| 1990s | 2,211 | 3,062 | 0.13 | 2.77 |
| 2000s | 5,118 | 8,890 | 0.07 | 3.47 |
| 2010s | 9,519 | 21,455 | 0.05 | 4.51 |下图比较了不同性别的共同作者网络度分布。女性的共同作者往往比男性少，但平均差异很小，并且随着时间的推移而下降。
# 数据处理与网络度数计算
plot_df = nets_df %>%
# 筛选出所有包含 "All" 项目的数据行，即汇总了所有项目的数据
filter(program == 'All') %>%
# 使用 map 函数对每个网络应用 gorder 函数，计算每个网络的节点数，并存储在新列 order 中
mutate(order = map(net, gorder),
# 对每个网络应用 degree 函数计算节点的度数，结果以数据框的形式存储在新列 res 中
res = map(net, ~tibble(deg = degree(.)))) %>%
# 选择 decade、nodes 和 res 列，为后续操作准备数据
select(decade, nodes, res) %>%
# 展开 nodes 和 res 列中的嵌套数据，使每行对应一个节点和其度数
unnest(c('nodes', 'res')) %>%
# 对每个十年的数据进行计数，计算该十年中每个网络的节点数量，并以 order 命名
add_count(decade, name = 'order') %>%
# 过滤掉 p_female 中有缺失值的行，确保性别数据的完整性
filter(!is.na(p_female)) %>%
# 创建新列 female，根据 p_female 判断节点代表的是男性还是女性，并标记
mutate(female = ifelse(p_female > 0.5, 'Females', 'Males'))
# 计算每个性别每个十年的平均度数
mean_degrees = plot_df %>%
# 按十年和性别分组
group_by(decade, female) %>%
# 计算每组的平均度数
summarise(mean = mean(deg)) %>%
ungroup() %>% # 解除分组
# 将不同性别的平均度数分布到不同的列
spread(female, mean) %>%
# 计算男性和女性平均度数的差异
mutate(diff = Males - Females) %>%
# 如果是数值类型，格式化为保留两位小数的字符串
mutate_if(is.numeric, ~sprintf('%.2f', .))
# 度数分布的可视化
plot_df %>%
# 将度数限制在10以内，超过10的视为10
mutate(deg = pmin(deg, 10)) %>%
# 对不同十年、顺序、度数、性别的组合计数
count(decade, order, deg, female) %>%
# 填充缺失的组合，使得每个组合至少出现一次
complete(nesting(decade, order), deg, female, fill = list(n = 0)) %>%
drop_na() %>% # 去除含有NA值的行
# 按十年和性别分组
group_by(decade, female) %>%
# 按度数排序
arrange(deg) %>%
# 计算每个度数的相对频率
mutate(p = n / sum(n)) %>%
ungroup() %>% # 解除分组
# 创建一个ggplot图，x轴为度数，y轴为频率
ggplot(aes(deg, p)) +
# 使用柱状图展示数据，不同性别用不同颜色表示，并使柱子并排显示
geom_col(aes(fill = female), position = 'dodge') +
facet_wrap(~decade) + # 按十年分面展示
labs(x = 'Co-authors during decade', # 设置x轴标签
y = 'Relative frequency', # 设置y轴标签
# 设置图表标题
title = 'Co-authorship degree distributions by decade',
# 设置图表副标题
subtitle = 'Females had slightly fewer co-authors than males',
# 不显示图例的标题
fill = NULL) +
# 防止绘图区域之外的内容被剪切
coord_cartesian(clip = 'off') +
# 设置图例的标签位置
guides(fill = guide_legend(label.position = 'left')) +
# 设置x轴的刻度和标签
scale_x_continuous(expand = c(0, 0), breaks = 2 * (0:5), labels = c(2 * (0:4), '10+')) +
# 设置y轴的扩展方式
scale_y_continuous(expand = c(0, 0)) +
# 设置图例的对齐方式
theme(legend.justification = c(1, 1),
# 设置图例的位置
legend.position = c(1, 1),
# 移除x轴的主要网格线
panel.grid.major.x = element_blank(),
# 设置面板之间的x轴间距
panel.spacing.x = unit(1, 'line'))
# 保存图表为PNG文件，指定尺寸和分辨率
ggsave('_media/Co-authorship_degree.png', width = 6, height = 6, dpi = 300)
接下来的三个表格将描述基于作者估计性别的每个十年的共作者网络的结构属性。这表明这些表格会展示不同十年间，根据作者性别估计得出的共作者网络的一些关键结构特征，如包括网络密度、聚类系数、中心性等，反映不同性别作者之间合作的模式和密度。这种分析有助于理解性别如何影响学术合作的结构和动态。
get_structural_properties = function(G, random = F) {
# 从图G中选取p_female非空的子图
G = induced_subgraph(G, !is.na(V(G)$p_female))
n = gorder(G) # 获取图的顶点数
if (random) {
# 如果是随机的，生成n个[0,1]均匀分布的随机数
t = runif(n)
} else {
# 否则为每个顶点分配0.5
t = rep(0.5, n)
}
# 生成女性作者的子图
Hf = induced_subgraph(G, t < V(G)$p_female)
# 生成非女性作者的子图
Hn = induced_subgraph(G, t >= V(G)$p_female)
tibble(
# 计算女性子图的聚类系数
clust_coeff_f = transitivity(Hf),
# 计算非女性子图的聚类系数
clust_coeff_nf = transitivity(Hn),
# 计算图G的同质性，基于性别
assort_coeff = assortativity(G, t < V(G)$p_female)
)
}
get_structural_property_confidence_intervals = function(G, n_reps = 1000, coverage = 0.95) {
# 重复n_reps次随机化过程
lapply(1:n_reps, function(x) get_structural_properties(G, random = T)) %>%
bind_rows() %>% # 绑定所有结果行
gather(key, x) %>% # 转换数据格式，准备进行统计
group_by(key) %>% # 按性质分组
summarise(lb = quantile(x, (1 - coverage) / 2, na.rm = T), # 计算下限
ub = quantile(x, (1 + coverage) / 2, na.rm = T)) %>% # 计算上限
ungroup() %>% # 解除分组
# 创建置信区间的字符串表示
mutate(ci = sprintf('(%.2f, %.2f)', lb, ub)) %>%
select(key, ci) %>% # 选择关键结果和置信区间
spread(key, ci) # 将数据转换为宽格式
}
### 设置随机种子保证可复制性
set.seed(0)
spci = nets_df %>%
# 选取特定项目和 All
filter(program %in% c(program_adj_counts$program_adj, 'All')) %>%
# 设置项目为因子变量
mutate(program = factor(program, c(program_adj_counts$program_adj, 'All')),
# 计算每个网络的聚类系数
clust_coeff = map_dbl(net, transitivity),
# 计算结构性质的置信区间
spci = map(net, get_structural_property_confidence_intervals)) %>%
# 移除不需要的列
select(-nodes, -edges, -net) %>%
unnest('spci') # 展开嵌套的数据框
### 计算女性作者比例的标准差是否大于0，按十年和论文分组
tmp = df %>%
filter(!is.na(p_female)) %>% # 移除p_female为空的行
mutate(female = p_female > 0.5) %>% # 创建性别标签
add_count(paper) %>% # 计算每篇论文的计数
filter(n > 1) %>% # 选取有多于一个作者的论文
group_by(decade, paper) %>% # 按十年和论文分组
summarise(x = sd(female) > 0) %>% # 计算性别的标准差是否大于0
group_by(decade) %>% # 按十年分组
summarise(x = round(100 * mean(x))) # 计算百分比并四舍五
### 使用kable展示聚类系数的表格
spci %>%
filter(program == 'All') %>%
mutate_if(is.numeric, ~sprintf('%.2f', .)) %>% # 格式化数字为两位小数
select(decade, starts_with('clust')) %>%
gather(key, value, -decade) %>% # 转换数据格式以便使用
spread(decade, value) %>% # 转换为宽格式
# 修改键的名称
mutate(key = c('Overall', 'Among females (95% CI)', 'Among males (95% CI)')) %>%
rename(`Clustering coefficient` = key) %>% # 重命名列
kable(align = 'lcccc') # 使用kable生成表格
### 展示配对系数
spci %>%
filter(program == 'All') %>%
select(Decade = decade, 4) %>%
# 生成配对系数表格
kable(align = 'c', col.names = c('Decade', 'Assort. coeff. (95% CI)'))
### 展示不同十年和项目的配对系数表格
spci %>%
select(decade, Program = program, assort_coeff) %>%
spread(decade, assort_coeff, fill = '-') %>% # 转换为宽格式
kable(align = 'c') # 使用kable生成表格聚类系数 (Clustering coefficient)：下表比较了每十年完整的合作作者网络的聚类系数与估计为女性和男性的作者集合所引起的子网络的聚类系数。女性子网络比整体网络和男性网络的聚类程度要高，这表明女性作者之间更倾向于与共享至少一个共同女性作者的其他女性合作，形成闭合的三元组 (triads)。此外，随着时间的推移，女性之间的聚类程度有所下降，这可能是由于跨性别合作的增加。
|Clustering coefficient | 1980s | 1990s | 2000s | 2010s |
|:----------------------|:------------:|:------------:|:------------:|:------------:|
|Overall | 0.17 | 0.18 | 0.21 | 0.24 |
|Among females (95% CI) | (0.39, 0.50) | (0.41, 0.50) | (0.30, 0.35) | (0.32, 0.35) |
|Among males (95% CI) | (0.16, 0.17) | (0.17, 0.17) | (0.20, 0.21) | (0.23, 0.23) |同质性系数 (Assortativity coefficient)：同质性系数用来衡量作者倾向于与同性别的其他作者合作的程度。系数值为 1 时表示完全同质 (没有跨性别合作)，为 -1 时表示完全异质 (没有同性别合作)，为 0 时则表示网络中的合作关系随机分布，无明显性别偏向。所有网络的同质性系数均为正值，这意味着同性别之间的合作比随机情况下更为常见。
| Decade | Assort. coeff. (95% CI) |
|:------:|:-----------------------:|
| 1980s | (0.05, 0.09) |
| 1990s | (0.08, 0.11) |
| 2000s | (0.07, 0.09) |
| 2010s | (0.08, 0.10) |跨项目同质性系数 (Assortativity coefficient)：进一步分析了 NBER 中十个最大研究项目的同质性系数，并提供了 95% 置信区间。分析发现，例如在劳动研究 (LS) 项目中，随着时间的推移，网络中的同质性趋向减少，而在健康经济学 (HE) 项目中，则趋向增加。这些变化与女性在这些领域的代表性增加有关，但这并不直接等同于促进女性合作的机制。
| Program | 1980s | 1990s | 2000s | 2010s |
|:-------:|:--------------:|:-------------:|:--------------:|:------------:|
| LS | (0.16, 0.27) | (0.13, 0.19) | (0.05, 0.09) | (0.09, 0.11) |
| EFG | (-0.07, 0.08) | (-0.07, 0.01) | (-0.02, 0.02) | (0.02, 0.06) |
| PE | (0.04, 0.14) | (-0.01, 0.05) | (0.03, 0.07) | (0.05, 0.07) |
| IFM | (-0.05, 0.04) | (-0.01, 0.08) | (-0.01, 0.05) | (0.03, 0.08) |
| ITI | (-0.06, 0.04) | (0.01, 0.09) | (0.00, 0.07) | (0.05, 0.10) |
| ME | (-0.07, 0.04) | (-0.03, 0.06) | (-0.10, -0.03) | (0.03, 0.09) |
| AP | - | (-0.06, 0.07) | (-0.01, 0.05) | (0.00, 0.05) |
| PR | (-0.15, -0.01) | (0.12, 0.22) | (0.02, 0.09) | (0.07, 0.11) |
| CF | - | (-0.04, 0.07) | (-0.03, 0.04) | (0.03, 0.09) |
| HE | (-0.14, -0.01) | (0.01, 0.09) | (0.01, 0.05) | (0.07, 0.10) |
| All | (0.05, 0.09) | (0.08, 0.11) | (0.07, 0.09) | (0.08, 0.10) |2.3.4 附录
下表 (按分数) 按计划和十年对工作文件进行了计数。按照四个十年中相关论文的降序排列项目。
|Program | 1980s| 1990s| 2000s| 2010s|
|:---------------------------------------------------|-----:|-----:|-----:|------:|
|Labor Studies (LS) | 454| 635| 868| 1,081|
|Economic Fluctuations and Growth (EFG) | 458| 471| 921| 1,083|
|Public Economics (PE) | 445| 557| 827| 993|
|International Finance and Macroeconomics (IFM) | 374| 466| 731| 662|
|International Trade and Investment (ITI) | 370| 517| 631| 525|
|Monetary Economics (ME) | 418| 327| 389| 514|
|Asset Pricing (AP) | 0| 221| 610| 627|
|Productivity, Innovation, and Entrepreneurship (PR) | 96| 231| 371| 563|
|Corporate Finance (CF) | 1| 131| 436| 560|
|Health Economics (HE) | 82| 115| 355| 497|
|Development of the American Economy (DAE) | 45| 78| 311| 379|
|Industrial Organization (IO) | 0| 82| 300| 432|
|Economics of Aging (AG) | 33| 126| 237| 340|
|Health Care (HC) | 0| 100| 248| 316|
|Environment and Energy Economics (EEE) | 1| 6| 138| 483|
|Economics of Education (ED) | 0| 1| 209| 411|
|Children (CH) | 2| 35| 246| 297|
|Political Economics (POL) | 0| 0| 141| 415|
|Law and Economics (LE) | 20| 57| 188| 231|
|Development Economics (DEV) | 0| 0| 0| 462|
|Technical Working Papers (TWP) | 0| 0| 25| 95|
|None | 24| 58| 7| 0|
|Total | 2,820| 4,213| 8,188| 10,970|
3. 结语
在本篇推文中，我们向大家介绍了R语言中的一些基本管道操作和正则表达式的使用方法，并通过复现一篇论文详细展示了如何在 R 中进行基础的数据处理和统计分析。希望这篇文章能帮助大家提升基础数据分析技能！
4. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh r语言, m
安装最新版lianxh命令：
ssc install lianxh, replace
- 专题：Stata命令
- 任建辉, 2022, Stata与R语言等价命令, 连享会 No.1083.
- 专题：交乘项-调节-中介
- 余坚, 2023, 中介效应：有序因果中介分析的半参数估计B-实操, 连享会 No.1179.
- 专题：面板数据
- 刘欣妍, 史柯, 2022, Stata：双向固定效应模型中是否要控制公司年龄？, 连享会 No.942.
- 专题：Python-R-Matlab
- 刘潍嘉, 2023, R语言：如何理解F检验, 连享会 No.1273.
- 吴小齐, 2023, R语言：L2 Boosting 在经济学中的应用, 连享会 No.1288.
- 吴小齐, 2023, R语言：将tikz图片转换为PNG, 连享会 No.1289.
- 专题：倍分法DID
- 张梓瑶, 2021, 如何在R语言中实现多期DID, 连享会 No.829.
- 专题：工具软件
- 李明来, 2021, R和RStudio安装教程, 连享会 No.773.
- 温世彬, 2022, R：在几十种计量软件之间转换数据格式-rio, 连享会 No.1014.
- 罗银燕, 2023, 如何在 R 中安装并使用 chatgpt 包？, 连享会 No.1171.
- 范思妤, 2024, R/RStudio下载安装指南, 连享会 No.1333.
- 专题：专题课程
- 连享会, 2023, Stata2R：针对 Stata 用户的 R 课程, 连享会 No.1265.
- 连玉君, 2023, 刺猬与狐狸：Stata与R语言常用命令等价对比, 连享会 No.1270.
- 邱一崎, 2024, R语言：因果推断中的混淆变量、交互项和线性回归, 连享会 No.1369.
- 郑晓雪, 2022, R语言绘制社会网络图, 连享会 No.953.
- 郑晓雪, 2022, R语言随机指数图分析-statnet, 连享会 No.963.
- 郭盼亭, 2023, R语言与Stata等价命令-statar, 连享会 No.1228.
- 闫钊鹏, 2023, R语言：gtExtras包-绘制漂亮的表格和图形, 连享会 No.1311.
- 陈卓然, 2022, Stata+R：合成DID原理及实现-sdid, 连享会 No.970.
- 陈点点, 2023, R：如何使用逆概率加权和贝叶斯倾向得分, 连享会 No.1279.
- 雷诺, 2023, R语言：拆解因果推断中的后门法则, 连享会 No.1323.
- 专题：生存分析
- 雷诺, 2023, 生存分析及R语言实现案例, 连享会 No.1324.
资源共享
- 连享会资料 ……
- 在线视频：lianxh-class.cn
- Stata 33 讲，100 万+ 播放，Stata 入门必备，公开课
- 直击面板数据模型，10 万+ 播放，白话面板模型，公开课
- … more …
- 论文复现和数据
- 主题分类
- 热门推文
尊敬的老师 / 亲爱的同学们：
连享会致力于不断优化和丰富课程内容，以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求，我们诚挚地邀请您参与到我们的课程规划中来。
请您在下面的问卷中，分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源，将直接影响到我们课程的改进和创新。
我们期待您的反馈，因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间，与我们共同塑造更加精彩的学习旅程！https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见！
关于我们
- Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。more……
- 扫码加入连享会微信群，提问交流更方便
