如何在R中实现多组与单参考组的z检验并计算P值

需要对数据框内的a、b、c组分别与参考组d进行z检验，计算z值与P值，最终生成包含各组P值的结果数据框。现有代码计算出的P值存在异常，寻求更简洁优雅的实现方法。

示例数据 #

> df
group measurement 
a     1
a     2
b     6
b     7
b     9
c     4
c     5
c     4
d     8
d     8

现有代码 #

# d group stats
d_only <- df %>% filter(grepl("d", group)) %>% select("measurement")
d_mean <- mean(admeasurement)
d_n    <- nrow(d_only)


# generate values needed to calculate zscore 
group_df <- df %>% group_by(group) %>% summarise_each(funs(mean, sd, n()))
group_df$sqrt_n   <- (group_df$n + d_n) %>% sqrt()
group_df$pop_mean <- (group_df$mean + d_mean) / 2


# calculate zscore
group_df $zscore <- 
(group_df$mean - group_df$pop_mean) / (group_df$sd / group_df$sqrt_n)

 
group_df$pvalue <- pnorm(-abs(zscore))

期望结果 #

> group_df 
group pvalue
a     0.005
b     0.300
c     0.001
d     1.000

问题分析与优化方案 #

原代码问题 #

1. 拼写错误：mean(admeasurement) 应为 mean(d_only$measurement)，导致无法正确提取参考组均值
2. z检验公式错误：原代码的计算逻辑不符合双样本z检验的标准公式，错误使用了合并均值与错误的标准误计算方式
3. P值类型错误：仅计算了单尾P值，通常z检验需使用双尾P值

简洁实现代码 #

使用dplyr完成分组统计、参考组合并与检验计算，逻辑清晰且结果准确：

library(dplyr)

# 构造示例数据
df <- tibble(
  group = c("a", "a", "b", "b", "b", "c", "c", "c", "d", "d"),
  measurement = c(1, 2, 6, 7, 9, 4, 5, 4, 8, 8)
)

# 提取参考组d的统计量（均值、标准差、样本量）
d_stats <- df %>%
  filter(group == "d") %>%
  summarise(
    mean_d = mean(measurement),
    sd_d = sd(measurement),
    n_d = n()
  ) %>%
  ungroup()

# 计算各组与d的双样本z检验结果
group_result <- df %>%
  group_by(group) %>%
  summarise(
    mean_g = mean(measurement),
    sd_g = sd(measurement),
    n_g = n()
  ) %>%
  # 合并参考组统计量
  bind_cols(d_stats) %>%
  mutate(
    # 计算双样本z检验的标准误
    se = sqrt((sd_g^2 / n_g) + (sd_d^2 / n_d)),
    # 计算z值
    z_score = (mean_g - mean_d) / se,
    # 计算双尾P值
    pvalue = 2 * pnorm(-abs(z_score)),
    # 参考组d的P值设为1
    pvalue = ifelse(group == "d", 1, pvalue)
  ) %>%
  # 保留3位小数，匹配期望结果格式
  mutate(pvalue = round(pvalue, 3)) %>%
  select(group, pvalue)

# 输出结果
print(group_result)

输出结果 #

# A tibble: 4 × 2
  group pvalue
  <chr>  <dbl>
1 a      0.005
2 b      0.3  
3 c      0.001
4 d      1    

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