在R中实现各时段下多会话的比例两两比较求助

针对每个时段的三组会话比例两两比较需求，我们采用Fisher精确检验（适配小样本计数场景）来计算p值，步骤如下：

1. 加载依赖包并整理数据 #

先加载dplyr处理数据，清理原始数据中的冗余列，同时定义各组总样本量：

library(dplyr)

# 原始数据
data <- tibble::tribble(
  ~Horaire, ~SessionM1, ~NbM1, ~SessionM2, ~NbM2, ~SessionM3, ~NbM3,
"18h-19h","M1", 6, "M2", 7, "M3", 6,
"19h-20h","M1", 6, "M2", 10, "M3", 17,
"20h-21h","M1", 14, "M2", 10, "M3", 16,
"21h-22h","M1", 23, "M2", 5, "M3", 9,
"22h-23h","M1", 15, "M2", 9, "M3", 4,
"23h-00h","M1", 19, "M2", 4, "M3", 6,
"00h-01h","M1", 19, "M2", 5, "M3", 10,
"01h-02h","M1", 21, "M2", 3, "M3", 7,
"02h-03h","M1", 8, "M2", 5, "M3", 7,
"03h-04h","M1", 15, "M2", 4, "M3", 10,
"04h-05h","M1", 15, "M2", 4, "M3", 9,
"05h-06h","M1", 20, "M2", 1, "M3", 5)

# 定义各组总样本量，清理冗余列
total_M1 <- 181
total_M2 <- 67
total_M3 <- 106

clean_data <- data %>% select(Horaire, NbM1, NbM2, NbM3)

2. 定义比例比较函数 #

编写函数实现两组比例的Fisher检验，返回p值：

compare_props <- function(n_success, n_total, n_success2, n_total2) {
  # 构建2x2列联表：[目标时段计数, 其他时段计数]
  contingency_table <- matrix(
    c(n_success, n_total - n_success, n_success2, n_total2 - n_success2),
    nrow = 2
  )
  # 执行Fisher精确检验
  test_result <- fisher.test(contingency_table)
  return(test_result$p.value)
}

3. 批量计算各时段两两对比p值 #

为每个时段添加三组对比的原始p值，同时可选添加Bonferroni校正p值（控制多重检验的Ⅰ类错误）：

result_data <- clean_data %>%
  mutate(
    # 原始p值
    p_M1M2 = compare_props(NbM1, total_M1, NbM2, total_M2),
    p_M2M3 = compare_props(NbM2, total_M2, NbM3, total_M3),
    p_M1M3 = compare_props(NbM1, total_M1, NbM3, total_M3),
    # Bonferroni校正p值（每个时段3次检验，校正因子为3）
    p_M1M2_corrected = p.adjust(p_M1M2, method = "bonferroni", n = 3),
    p_M2M3_corrected = p.adjust(p_M2M3, method = "bonferroni", n = 3),
    p_M1M3_corrected = p.adjust(p_M1M3, method = "bonferroni", n = 3)
  )

# 查看结果（保留4位小数）
print(result_data, digits = 4)

关键说明 #

• Fisher精确检验比卡方检验更适合小计数场景，避免近似偏差；
• 若需严格控制整体显著性水平，建议使用校正后的p值；探索性分析可直接查看原始p值；
• p值小于0.05通常认为两组比例差异具有统计学显著性。

内容的提问来源于stack exchange，提问作者H. Ravao