配对数据均值/几何均值计算异常：mutate调用结果错误求助

需求背景 #

有一组按剂量随机分组的患者数据，两天内测量了50余项参数，部分参数仅在部分患者中测量。需生成箱线图：

• 按剂量比较天数、按天数比较剂量
• 图中添加均值/几何均值点：剂量间比较用总体均值/几何均值，天数间比较用配对均值（仅计算两天均有测量值的患者）

自行编写了means_long函数，单独调用时结果正确，但在dplyr的mutate分组调用时，配对均值始终被计算为总体均值，存在错误。

模拟数据 #

subject <- c(paste0("sub", 1:40),
             paste0("sub", 1:40))

day <- c(rep("day 1", times = 40),
         rep("day 2", times = 40))

ps <- rep(1, times = 80)

dose <- c(rep("low", times = 10),
          rep("high", times = 10),
          rep("low", times = 10),
          rep("high", times = 10),
          rep("low", times = 10),
          rep("high", times = 10),
          rep("low", times = 10),
          rep("high", times = 10))

param1 <- c(NA, NA, 3:20, rep(NA, times = 20),
            1.8*(1:20), rep(NA, times = 20))

param2 <- c(0.75*(1:40), 1.05*(1:40))

s_data <- data.frame(subject, day, ps, dose, param1, param2)

参数设置（问题场景） #

param <- "param1"  # 要绘图的参数
by <- "dose"       # 子图拆分依据
what <- "day"      # 比较的维度
by2 <- "low"       # 目标剂量组
geometric <- FALSE # 是否计算几何均值
paired <- TRUE     # 是否用配对均值（day比较时设为TRUE，dose比较时设为FALSE）

问题函数：means_long #

means_long <- function(values, group, id_col, grouping_var, is_geometric, is_paired)  {
    
  dummy_df <- tibble(id_col, values, grouping_var) %>% 
  pivot_wider(id_cols = id_col,
              names_from = grouping_var,
              values_from = values)
  
  if (is_paired) {
    dummy_df <- na.omit(dummy_df)
  }
   
  if (any(na.omit(dummy_df[[group]]) <= 0) & is_geometric) {
    return(NA)
  } else {
    est <- if_else(is_geometric,
                   exp(mean(log(dummy_df[[group]]), na.rm = TRUE)),
                   mean(dummy_df[[group]], na.rm = TRUE)
                   )
    return(est)
  }
}

调用差异与错误表现 #

1. mutate分组调用（结果错误） #

plot <- s_data %>% 
  filter(.data[[by]] == by2) %>%
  group_by(!!sym(what)) %>% 
  mutate(no = paste("n = ", sum(!is.na(!!sym(param)))),
         mean_v = means_long(!!sym(param),
                             cur_group()[[what]],
                             subject,
                             !!sym(what),
                             geometric,
                             paired)
  ) %>%
  ggplot() + 
  # 绘图图层省略

得到错误结果：day1均值6.5，day2均值9.9（应为配对均值11.7）

2. 单独调用（结果正确） #

s_data_dose_l <- s_data %>%
  filter(dose == "low")

mean_p_day1 <- means_long(s_data_dose_l[[param]], 
                        "day 1",
                        s_data_dose_l$subject, 
                        s_data_dose_l[[what]],
                        geometric,
                        paired)

mean_p_day2 <- means_long(s_data_dose_l[[param]], 
                        "day 2",
                        s_data_dose_l$subject, 
                        s_data_dose_l[[what]],
                        geometric,
                        paired)

得到正确配对均值：day1为6.5，day2为11.7

诊断发现 #

当is_paired=TRUE时，两种调用方式下na.omit处理后的数据集行数不同：

• mutate分组调用时，day2对应的dummy_df有10行
• 单独调用时，day2对应的dummy_df有8行

尝试过na.omit、complete.cases、drop_na等多种去NA方法，问题依旧。

问题原因 #

在group_by(day)后使用mutate调用函数时，每个分组（day1、day2）仅传入当前组的子集数据。而配对均值需要整个剂量分组下的两天数据来筛选出两天均有值的患者，再计算对应组的均值。但分组后mutate只传递当前组数据，导致pivot_wider只能得到当前组的患者ID和对应值，无法跨组筛选配对患者，na.omit无法剔除仅单天有值的患者，最终计算的是当前组的总体均值而非配对均值。

解决办法 #

提前在完整的剂量分组数据中计算好配对均值，再与原数据合并后绘图，避免分组后的数据限制。

修改后的代码示例 #

# 1. 提前计算配对均值（筛选两天都有值的患者）
mean_data <- s_data %>%
  filter(dose == by2) %>%
  # 宽转长，保留患者ID和两天的参数值
  pivot_wider(id_cols = subject, names_from = day, values_from = all_of(param)) %>%
  drop_na() %>% # 仅保留两天都有值的患者
  # 转回长格式，方便分组计算均值
  pivot_longer(cols = -subject, names_to = what, values_to = param) %>%
  group_by(!!sym(what)) %>%
  summarise(
    mean_v = if_else(geometric,
                     exp(mean(log(.data[[param]]), na.rm = TRUE)),
                     mean(.data[[param]], na.rm = TRUE)
    ),
    no = paste("n = ", n()) # 配对样本量
  )

# 2. 合并原数据与均值数据
plot_data <- s_data %>%
  filter(dose == by2) %>%
  left_join(mean_data, by = what)

# 3. 绘图
ggplot(plot_data, aes(x = !!sym(what), y = !!sym(param))) +
  geom_boxplot() +
  geom_point(aes(y = mean_v), color = "red", size = 3) + # 添加均值点
  geom_text(aes(label = no), y = Inf, vjust = 1.5) # 添加样本量文本

补充说明 #

如果需要支持剂量间的总体均值计算，只需调整mean_data的计算逻辑：去掉drop_na()步骤，直接按剂量分组计算总体均值即可。

内容的提问来源于stack exchange，提问作者Radek Jaźwiec