含缺失值的多分组计数/比例显著性检验方法问询

给定如下示例数据集：

data.frame(
  Treatment = c("A", "A", "A", "A", "A", "A",
                "A", "A", "A", "A", "A", "A",
                "B", "B", "B", "B", "B", "B",
                "B", "B", "B", "B", "B", "B"),
  Patient = c(1, 1, 1, 1, 1, 1,
              2, 2, 2, 2, 2, 2,
              3, 3, 3, 3, 3, 3,
              4, 4, 4, 4, 4, 4),
  Timepoint = c("PRE", "PRE", "PRE", "POST", "POST", "POST",
                "PRE", "PRE", "PRE", "POST", "POST", "POST",
                "PRE", "PRE", "PRE", "POST", "POST", "POST",
                "PRE", "PRE", "PRE", "POST", "POST", "POST"),
  Phenotype = c("NK", "T Cell", "Macrophage", "NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage", "NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage", "NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage", "NK", "T Cell", "Macrophage"),
  Count = c(523,235,2352,352,646,234,
            3463,525,646,234,725,264,
            1636,3153,455,134,646,253,
            464,252,464,276,364,353)
)

需要完成两层显著性对比：

1. 患者层面：针对每个患者的每种细胞表型，对比PRE与POST时间点的计数/比例，输出格式如下：

data.frame(
  Patient = c(1, 1, 1,
              2, 2, 2,
              3, 3, 3,
              4, 4, 4),
  Phenotype = c("NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage"),
  Pvalue = c(0, 0, 0,
             0, 0, 0,
             0, 0, 0,
             0, 0, 0)
)

2. 治疗组层面：基于Treatment分组进行更高层面的对比，输出格式如下：

data.frame(
  Treatment = c("A", "A", "A",
              "B", "B", "B"),
  Phenotype = c("NK", "T Cell", "Macrophage",
                "NK", "T Cell", "Macrophage"),
  Pvalue = c(0, 0, 0,
             0, 0, 0)
)

实际数据中存在部分患者的某表型仅在一个时间点有计数（另一个为0或NA）的情况，导致常规卡方检验代码执行失败，同时不确定比例检验或卡方检验的适配性，需要针对含缺失值的大规模数据集，给出批量完成两层检验的高效方案。

一、检验方法选择 #

• 患者层面：优先用McNemar检验，这是专门针对配对二分类数据的检验，适配同一患者前后两次观测的比例差异；
• 治疗组层面：可选择卡方检验/Fisher精确检验（聚合组内数据后），或用混合效应模型整合患者个体差异，结果更稳健。

核心前提：先补全缺失的时间点记录（将0填充进无数据的单元格），再处理检验的边界情况。

二、患者层面显著性检验实现 #

结合dplyr和tidyr补全数据，自定义函数处理0计数/缺失的边界情况：

library(dplyr)
library(tidyr)

# 自定义McNemar检验函数，处理0计数和缺失值
run_mcnemar <- function(pre_count, post_count, total_pre, total_post) {
  # 补全0值
  pre <- ifelse(is.na(pre_count), 0, pre_count)
  post <- ifelse(is.na(post_count), 0, post_count)
  
  # 构建配对列联表：该表型计数 vs 其他表型计数
  table <- matrix(
    c(pre, total_pre - pre, post, total_post - post),
    nrow = 2,
    dimnames = list(PRE = c("Target", "Others"), POST = c("Target", "Others"))
  )
  
  # 边界情况处理：所有观测一致时p值为1；行列全0时直接返回1
  if (all(table[,1] == 0) || all(table[,2] == 0) || all(table[1,] == 0) || all(table[2,] == 0)) {
    return(1)
  }
  
  # 期望频数<5时用精确检验，否则用常规McNemar检验
  if (any(table < 5)) {
    test <- mcnemar.test(table, correct = FALSE)
  } else {
    test <- mcnemar.test(table)
  }
  return(test$p.value)
}

# 患者层面检验流程
patient_level_pvals <- df %>%
  # 补全每个患者-表型的PRE/POST记录，缺失Count填0
  complete(Patient, Phenotype, Timepoint, fill = list(Count = 0)) %>%
  # 计算每个患者-时间点的总细胞数
  group_by(Patient, Timepoint) %>%
  mutate(total_time_count = sum(Count)) %>%
  ungroup() %>%
  # 按患者和表型分组
  group_by(Patient, Phenotype) %>%
  summarise(
    pre_count = Count[Timepoint == "PRE"],
    post_count = Count[Timepoint == "POST"],
    total_pre = total_time_count[Timepoint == "PRE"],
    total_post = total_time_count[Timepoint == "POST"],
    .groups = "drop"
  ) %>%
  # 计算p值
  mutate(Pvalue = mapply(run_mcnemar, pre_count, post_count, total_pre, total_post)) %>%
  select(Patient, Phenotype, Pvalue)

# 查看结果
patient_level_pvals

三、治疗组层面显著性检验实现 #

方法1：聚合后卡方/Fisher检验 #

适合快速组间对比，处理边界情况：

# 自定义组间检验函数
run_group_test <- function(treatment, phenotype) {
  # 提取目标组-表型的PRE/POST聚合数据
  agg_data <- df %>%
    filter(Treatment == treatment, Phenotype == phenotype) %>%
    group_by(Timepoint) %>%
    summarise(target_total = sum(Count), .groups = "drop") %>%
    complete(Timepoint, fill = list(target_total = 0))
  
  # 计算该治疗组对应时间点的总细胞数
  total_pre <- df %>% filter(Treatment == treatment, Timepoint == "PRE") %>% summarise(sum(Count)) %>% pull()
  total_post <- df %>% filter(Treatment == treatment, Timepoint == "POST") %>% summarise(sum(Count)) %>% pull()
  
  # 构建列联表
  table <- matrix(
    c(agg_data$target_total[1], total_pre - agg_data$target_total[1],
      agg_data$target_total[2], total_post - agg_data$target_total[2]),
    nrow = 2,
    dimnames = list(Timepoint = c("PRE", "POST"), Type = c("Target", "Others"))
  )
  
  # 边界情况处理
  if (all(table[,1] == 0) || all(table[,2] == 0) || all(table[1,] == 0) || all(table[2,] == 0)) {
    return(1)
  }
  
  # 选择检验方法
  if (any(table < 5)) {
    test <- fisher.test(table)
  } else {
    test <- chisq.test(table, correct = FALSE)
  }
  return(test$p.value)
}

# 治疗组层面检验流程
treatment_level_pvals <- df %>%
  distinct(Treatment, Phenotype) %>%
  mutate(Pvalue = mapply(run_group_test, Treatment, Phenotype)) %>%
  select(Treatment, Phenotype, Pvalue)

# 查看结果
treatment_level_pvals

方法2：混合效应模型（更稳健） #

适合大规模数据集，整合患者个体差异：

library(lme4)
library(purrr)
library(tibble)

# 预处理数据：补全缺失值，计算比例
model_data <- df %>%
  complete(Treatment, Patient, Phenotype, Timepoint, fill = list(Count = 0)) %>%
  group_by(Patient, Timepoint) %>%
  mutate(total_count = sum(Count), Prop = Count / total_count) %>%
  ungroup()

# 针对每个表型拟合混合效应模型，提取交互项p值
treatment_level_pvals_model <- map_dfr(unique(model_data$Phenotype), function(p) {
  sub_data <- filter(model_data, Phenotype == p)
  # 拟合二项混合效应模型，加入患者随机效应
  model <- glmer(Prop ~ Treatment * Timepoint + (1|Patient), data = sub_data, family = binomial)
  # 检验Treatment与Timepoint的交互效应
  anova_result <- anova(model, test = "Chisq")
  p_val <- anova_result$`Pr(>Chisq)`[3]
  
  tibble(
    Treatment = unique(sub_data$Treatment),
    Phenotype = p,
    Pvalue = p_val
  )
})

# 查看结果
treatment_level_pvals_model

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