如何用R的purrr与dplyr包为多列分组计算置信区间

需求说明 #

在starwars数据集中，按skin_color分组，为mass、birth_year等多列变量（实际场景有50+变量）计算置信区间的上下限。

已完成步骤 #

1. 按skin_color分组，计算各变量的非NA观测数

pacman::p_load("tidyr","purrr", "dplyr")

data <- starwars
data_obs <- 
  data %>%
  dplyr::select(mass,birth_year,skin_color) %>%
  dplyr::group_by(skin_color)%>%
  dplyr::summarise_all(funs(sum(!is.na(.))))%>%
  dplyr::ungroup()

2. 按skin_color分组，计算各变量的均值和标准差

data_stats <- 
  data %>%
  dplyr::select(mass,birth_year,skin_color)%>%
  dplyr::group_by(skin_color) %>%
  dplyr::summarise_all(., list(mean,sd), na.rm=T)%>%
  dplyr::ungroup()

3. 合并数据集，得到各变量的观测数、均值、标准差

data_complete <-
  dplyr::inner_join(data_obs,data_stats, by="skin_color")

4. 手动计算单个变量（mass）的置信区间

data_complete <-
  dplyr::mutate(data_complete,
                mass_se = mass_sd/sqrt(mass), 
                mass_ci_upper = mass_mean + qt(1 - (0.05 / 2), mass - 1)*mass_se,
                mass_ci_lower = mass_mean - qt(1 - (0.05 / 2), mass - 1)*mass_se)

遇到的问题 #

实际数据集有50+变量，手动逐个计算效率极低。尝试用purrr批量处理但代码失败：

list_vectors <- list(data$mass,data$birth_year)
list_ready <- map(list_vectors, 
                  ~ data %>%
                    group_by(skin_color)%>%
                    dplyr::summarise_all(funs(sum(!is.na(.))))%>%
                    dplyr::summarise_all(., list(mean,sd), na.rm=T) %>%
                    dplyr::ungroup()%>%
                    dplyr::mutate(var_se=var_sd/sqrt(var_n))) 

vector_1 <- list_ready[[1]]

解决方案 #

推荐两种高效的批量处理方式，适配不同的输出格式需求：

方法1：嵌套数据框+map（长格式输出，适合分析/可视化） #

这种方式结构清晰，便于扩展到任意数量的变量：

# 加载所需包
pacman::p_load(tidyverse)

# 1. 整理数据：保留目标变量和分组变量，按skin_color嵌套
nested_data <- starwars %>%
  select(skin_color, mass, birth_year) %>% # 替换为你的变量列表，比如matches("your_var_pattern")
  group_by(skin_color) %>%
  nest()

# 2. 定义置信区间计算函数
calc_ci <- function(df, vars) {
  map_dfr(vars, function(var) {
    # 提取变量数据并过滤NA
    var_data <- df[[var]] %>% na.omit()
    n <- length(var_data)
    
    # 处理样本量不足的情况
    if(n < 2) {
      return(tibble(
        variable = var,
        n = n,
        mean = NA_real_,
        sd = NA_real_,
        se = NA_real_,
        ci_lower = NA_real_,
        ci_upper = NA_real_
      ))
    }
    
    # 计算统计量和置信区间
    mean_val <- mean(var_data)
    sd_val <- sd(var_data)
    se_val <- sd_val / sqrt(n)
    t_crit <- qt(0.975, df = n - 1)
    ci_lower <- mean_val - t_crit * se_val
    ci_upper <- mean_val + t_crit * se_val
    
    tibble(
      variable = var,
      n = n,
      mean = mean_val,
      sd = sd_val,
      se = se_val,
      ci_lower = ci_lower,
      ci_upper = ci_upper
    )
  })
}

# 3. 批量应用函数并展开结果
result <- nested_data %>%
  mutate(ci_results = map(data, ~ calc_ci(.x, vars = c("mass", "birth_year")))) %>%
  unnest(ci_results) %>%
  select(-data)

方法2：purrr循环变量名（宽格式输出，匹配手动计算的结构） #

如果需要保持每个变量的统计量为单独列，用此方法：

# 定义单个变量的统计量计算函数
get_var_stats <- function(var_name) {
  starwars %>%
    select(skin_color, all_of(var_name)) %>%
    group_by(skin_color) %>%
    summarise(
      !!paste0(var_name, "_n") := sum(!is.na(.data[[var_name]])),
      !!paste0(var_name, "_mean") := mean(.data[[var_name]], na.rm = TRUE),
      !!paste0(var_name, "_sd") := sd(.data[[var_name]], na.rm = TRUE),
      !!paste0(var_name, "_se") := sd(.data[[var_name]], na.rm = TRUE)/sqrt(sum(!is.na(.data[[var_name]]))),
      !!paste0(var_name, "_ci_lower") := mean(.data[[var_name]], na.rm = TRUE) - qt(0.975, sum(!is.na(.data[[var_name]]))-1)*sd(.data[[var_name]], na.rm = TRUE)/sqrt(sum(!is.na(.data[[var_name]]))),
      !!paste0(var_name, "_ci_upper") := mean(.data[[var_name]], na.rm = TRUE) + qt(0.975, sum(!is.na(.data[[var_name]]))-1)*sd(.data[[var_name]], na.rm = TRUE)/sqrt(sum(!is.na(.data[[var_name]]))),
      .groups = "drop"
    )
}

# 批量处理变量并合并结果
vars_to_calc <- c("mass", "birth_year")
result_wide <- map(vars_to_calc, get_var_stats) %>%
  reduce(inner_join, by = "skin_color")

说明 #

• 两种方法都支持扩展到50+变量，只需修改vars_to_calc或select中的变量列表即可。
• 函数中加入了样本量不足的判断，避免因n<1导致的计算错误。

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