使用dplyr和条件语句在R中计算均值与置信区间

问题说明 #

原代码的核心错误：

• 计算置信区间时使用分组总样本量n()，而非对应变量剔除NA后的有效样本量
• 部分mean()和sd()计算未添加na.rm=TRUE，导致结果出现NA值

修正代码 #

library(dplyr)
library(tidyr)

# 第一步：分组计算各指标的均值、有效样本量及置信区间
summary_df <- df %>%
  group_by(name, region) %>%
  summarise(
    # 处理x变量（患病率）
    mean_x = mean(x == 1, na.rm = TRUE) * 100,
    n_x = sum(!is.na(x)),
    se_x = if(n_x > 1) sd(x == 1, na.rm = TRUE) / sqrt(n_x) else 0,
    lower_x = mean_x - qt(1 - 0.05/2, df = n_x - 1) * se_x,
    upper_x = mean_x + qt(1 - 0.05/2, df = n_x - 1) * se_x,
    
    # 处理y变量（患病率）
    mean_y = mean(y == 1, na.rm = TRUE) * 100,
    n_y = sum(!is.na(y)),
    se_y = if(n_y > 1) sd(y == 1, na.rm = TRUE) / sqrt(n_y) else 0,
    lower_y = mean_y - qt(1 - 0.05/2, df = n_y - 1) * se_y,
    upper_y = mean_y + qt(1 - 0.05/2, df = n_y - 1) * se_y,
    
    # 处理z变量（均值）
    mean_z = mean(z, na.rm = TRUE),
    n_z = sum(!is.na(z)),
    se_z = if(n_z > 1) sd(z, na.rm = TRUE) / sqrt(n_z) else 0,
    lower_z = mean_z - qt(1 - 0.05/2, df = n_z - 1) * se_z,
    upper_z = mean_z + qt(1 - 0.05/2, df = n_z - 1) * se_z,
    .groups = "drop"
  )

# 第二步：转换为目标整洁表格格式
df1 <- summary_df %>%
  pivot_longer(
    cols = -c(name, region),
    names_to = c(".value", "Indicator"),
    names_pattern = "(mean|lower|upper|n)_(.)"
  ) %>%
  select(name, region, Indicator, `mean/prevalence` = mean, `Lower interval` = lower, `Upper interval` = upper) %>%
  arrange(name, region, Indicator)

输出结果示例 #

运行上述代码后，df1的部分结果如下：

# A tibble: 18 × 5
   name  region Indicator `mean/prevalence` `Lower interval` `Upper interval`
   <chr> <chr>  <chr>                  <dbl>            <dbl>            <dbl>
 1 A     reg1   x                       33.3            -110.             176.
 2 A     reg1   y                      100.              100.             100.
 3 A     reg1   z                        2.33             -6.20            10.9
 4 A     reg2   x                       66.7             -38.8            172.
 5 A     reg2   y                      100.              100.             100.
 6 A     reg2   z                        2.33             -1.67             6.33

关键修正点 #

• 用sum(!is.na(x))计算各变量的有效样本量，替代原代码的n()
• 所有mean()和sd()添加na.rm=TRUE确保忽略NA值
• 添加if(n_x >1)判断避免样本量为1时自由度为0的报错
• 通过pivot_longer的names_pattern参数批量转换列名，简化格式整理

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