如何编写R函数批量生成子组与整体收入核密度对比图并保存PDF

现有如下数据框：

set.seed(1)
d <- data.frame(year= c(2001:2005,2001:2005,2001:2005),
                income = sample(2000:10000,15,replace = T),
                gender = sample(1:2,15,replace = T),
                education = sample(1:3,15,replace = T)
)

需求：编写函数，针对gender、education等所有指定分组变量的每个子组，绘制该子组与整体的收入核密度对比图，并将每个子组的图单独保存为PDF文件。

手动实现示例（以gender==1为例）：

male <- d %>% filter(gender == 1)

density_all <- density(d$income)
density_male <- density(male$income)

d_densisty <- data.frame(x = density_all$x, 
                      density_all = density_all$y, 
                      density_male = density_male$y)

plot <- ggplot(d_densisty, aes(x)) + 
  geom_line(aes(y = density_all), color = "red") +
  geom_line(aes(y = density_male), color = "blue")

ggsave("subgroup_name.pdf", plot, width = 300, height = 250, units = "mm")  

遇到的困惑：尝试过长格式转宽，但子组长度不一致；考虑过嵌套循环，但不确定最优方案及具体实现。

推荐使用嵌套循环+密度值对齐的方案，既直观又能解决子组密度x轴不一致的问题。核心思路：遍历每个分组变量→遍历该变量的每个子组→对齐整体与子组的密度x轴→绘图保存。

完整实现代码 #

library(dplyr)
library(ggplot2)

# 定义生成密度对比图的函数
plot_density_comparison <- function(data, group_vars, income_col = "income") {
  # 遍历每个分组变量
  for (var in group_vars) {
    # 获取当前变量的所有唯一取值
    group_values <- unique(data[[var]])
    # 遍历每个子组取值
    for (val in group_values) {
      # 提取子组数据
      subgroup_data <- data %>% filter(.data[[var]] == val)
      # 计算整体收入的密度
      density_all <- density(data[[income_col]])
      # 计算子组收入的密度，并对齐到整体的x轴上
      density_subgroup <- approx(density(subgroup_data[[income_col]]), xout = density_all$x)
      
      # 合并成绘图用的数据框
      plot_data <- data.frame(
        x = density_all$x,
        整体 = density_all$y,
        子组 = density_subgroup$y
      )
      
      # 绘制密度对比图
      p <- ggplot(plot_data, aes(x = x)) +
        geom_line(aes(y = 整体), color = "red", linewidth = 1) +
        geom_line(aes(y = 子组), color = "blue", linewidth = 1) +
        labs(title = paste0("收入密度对比：", var, " = ", val),
             x = "收入", y = "密度") +
        theme_minimal()
      
      # 保存PDF，文件名包含变量名和子组值
      ggsave(paste0(var, "_", val, ".pdf"), p, width = 300, height = 250, units = "mm")
    }
  }
}

# 调用函数，指定要处理的分组变量
plot_density_comparison(d, group_vars = c("gender", "education"))

关键细节说明 #

1. 密度x轴对齐：用approx()函数将子组的密度值插值到整体密度的x轴上，解决不同子组密度x序列长度不一致的问题，确保能合并到同一数据框绘图。
2. 动态变量引用：用.data[[var]]实现对分组变量的动态引用，避免硬编码，让函数支持任意数量的分组变量。
3. 命名规范：PDF文件名采用变量名_子组值.pdf格式，方便后续识别和管理。

扩展建议 #

• 若需调整线条样式、颜色或标题，直接修改ggplot()代码块即可。
• 若子组样本量过小导致密度估计不稳定，可在density()中指定adjust参数调整平滑程度，比如density(..., adjust = 2)。

内容的提问来源于stack exchange，提问作者MicL