使用 R 语言绘制分段回归线+拟合方程文本标签以及如何导出可以插入 word 文档中的矢量图

最近有个小伙伴想绘制一幅类似这样的图表：

这里我们以 1952～2020 年中国历年 GDP 的数据为例进行演示。

读取数据：

library(tidyverse)
haven::read_dta("中国历年GDP.dta") -> df
df %>%
  select(year = 年份, gdp = 地区生产总值_亿元) -> df
df

#> # A tibble: 69 × 2
#>     year   gdp
#>    <dbl> <dbl>
#>  1  1952  679.
#>  2  1953  824.
#>  3  1954  860.
#>  4  1955  912.
#>  5  1956 1031.
#>  6  1957 1071.
#>  7  1958 1312.
#>  8  1959 1448.
#>  9  1960 1470.
#> 10  1961 1232.
#> # ℹ 59 more rows

根据 1978 和 2001 两个时间点把数据分成 3 段：

df %>%
  mutate(gdp = gdp / 10000,
         group = cut(year, breaks = c(1951, 1978, 2001, 2021))) -> df1

df1 %>%
  count(group)

#> # A tibble: 3 × 2
#>   group                   n
#>   <fct>               <int>
#> 1 (1.95e+03,1.98e+03]    27
#> 2 (1.98e+03,2e+03]       23
#> 3 (2e+03,2.02e+03]       19

ggpmisc 包提供了直接添加拟合方程文本标签的图层，ggrepel 提供的 geom_text_repel 图层可以避免文本标签重叠问题。

library(ggpmisc)
library(ggrepel)

把三个时间组记录下，方便等下写代码，因为上面用了科学技术法，所以这里也用上：

group1 <- "(1.95e+03,1.98e+03]"
group2 <- "(1.98e+03,2e+03]"
group3 <- "(2e+03,2.02e+03]"

然后就可以使用 ggplot2 绘图了，首先分别使用 geom_point、geom_smooth、stat_poly_eq 三种图层分别绘制三个时间组的图：

ggplot() +
  geom_point(data = subset(df1, group == group1),
             aes(x = year, y = gdp, color = "#e64b35"),
             shape = 15, size = 1.2) +
  geom_smooth(data = subset(df1, group == group1),
              aes(x = year, y = gdp), color = "#e64b35",
              method = "lm", formula = y ~ x,
              linetype = "dashed",
              linewidth = 0.8) +
  stat_poly_eq(data = subset(df1, group == group1),
               aes(x = year, y = gdp + 5,
                   label =  paste(after_stat(eq.label),
                                  after_stat(rr.label),
                                  after_stat(p.value.label),
                                  sep = "*\", \"*")),
               formula = y ~ x, parse = TRUE, color = "#e64b35",
               geom = "text_repel", size = 3) +
  geom_point(data = subset(df1, group == group2),
             aes(x = year, y = gdp, color = "#00a087"),
             shape = 15, size = 1.2) +
  geom_smooth(data = subset(df1, group == group2),
              aes(x = year, y = gdp), color = "#00a087",
              method = "lm", formula = y ~ x,
              linetype = "dashed",
              linewidth = 0.8) +
  stat_poly_eq(data = subset(df1, group == group2),
               aes(x = year + 12, y = gdp + 3.5,
                   label = paste(after_stat(eq.label),
                                 after_stat(rr.label),
                                 after_stat(p.value.label),
                                 sep = "*\", \"*")),
               formula = y ~ x, parse = TRUE, color = "#00a087",
               geom = "text_repel", size = 3, angle = 4) +
  geom_point(data = subset(df1, group == group3),
             aes(x = year, y = gdp, color = "#3c5488"),
             shape = 15, size = 1.2) +
  geom_smooth(data = subset(df1, group == group3),
              aes(x = year, y = gdp), color = "#3c5488",
              method = "lm", formula = y ~ x,
              linetype = "dashed",
              linewidth = 0.8) +
  stat_poly_eq(data = subset(df1, group == group3),
               aes(x = year + 5.5, y = gdp - 35,
                   label = paste(after_stat(eq.label),
                                 after_stat(rr.label),
                                 after_stat(p.value.label),
                                 sep = "*\", \"*")),
               formula = y ~ x, parse = TRUE, color = "#3c5488",
               geom = "text_repel", size = 3, angle = 60) -> p1
p1

然后添加总曲线，使用二次多项式拟合：

# 这里的 cnfont 是我在 profile 里面设置的字体，大家可以参考系列课程 R 语言数据科学的第一次课进行设置。
p1 +
  geom_smooth(data = df1, formula = y ~ poly(x, 2),
              aes(x = year, y = gdp), linewidth = 1) +
  scale_x_continuous(breaks = seq(1952, 2022, by = 5)) +
  scale_color_manual(values = c("#e64b35" = "#e64b35",
                                "#00a087" = "#00a087",
                                "#3c5488" = "#3c5488"),
                     labels = c("1952～1978年，线性显著",
                                "1979～2001年，线性显著",
                                "2002～2020年，线性显著"),
                     name = "") +
  theme_classic(base_family = cnfont) +
  theme(legend.position = c(0.15, 0.8),
        plot.background = element_rect(fill = "white", color = "white")) +
  labs(x = "年份", y = "GDP(万亿元)") -> p

ggsave("pic.png", width = 8.3, height = 5, device = png)
ggsave("pic.pdf", width = 8.3, height = 5, device = pdf)

另外大家经常会使用 word 撰写论文，如何在 word 中插入矢量图呢？word 中支持的矢量图格式主要是 emf，在 R 语言中可以使用 devEMF 包提供的 emf 设备导出：

library(devEMF)
emf(file = "pic.emf", width = 20, height = 12, units = "cm")
source("pic.R")
p
dev.off()

#> quartz_off_screen
#>                 2

另外这些长宽参数可能需要反复测试才能找到适合自己的。pic.R 文件存储的是经过调整文本大小等内容的绘图代码。

也可以使用下面的格式保存：

source("pic.R")
ggsave("pic2.emf", device = devEMF::emf, width = 20, height = 12, units = "cm")

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制分段回归线+拟合方程文本标签以及如何导出可以插入 word 文档中的矢量图