如何基于汇总数据（均值、标准差）计算Scheffé检验的P值？

问题背景 #

需要基于汇总数据（均值、标准差）计算Scheffé检验的P值，实现类似R语言DescTools包中sheffeTest的功能，完成所有可能的成对比较，但未找到相关文档。

现有汇总数据Data_sum：

Location  n       mean       sd
<fctr>    <dbl>   <dbl>      <dbl>
TM    5   0.08016    0.0145622
NP    5   0.07550    0.0091643
PB    5   0.10434    0.0195800

所需成对对比表（与emmeans(model, pairwise ~ Location, adjust="scheffe")结果一致）：

Contrast.Name     TM   NP    PB
 TMvsNP           -1    1     0
 TMvsPB           -1    0     1
 NPvsPB            0   -1     1

已使用HH包的anovaMean()完成方差分析，能手动计算Scheffé检验，但无法计算P值，需解决思路。

问题解决 #

解决步骤 #

1. 基于汇总统计量（均值、标准差、样本量）计算方差分析
2. 构建所有成对对比项
3. 计算Scheffé检验的P值及相关统计量

注意： 由于代码中使用pooled_n = mean(n_precalcs)，非平衡设计下结果可能不准确

实现代码 #

library(HH)
library(magrittr)
library(rstatix)

scheffe.test <- function(iv,n,ybar,sd, alpha = 0.05){
  if (sd(n) != 0) warning("非平衡设计，结果可能存在偏差!\n")
  iv = factor(iv,levels = unique(iv)) # 将自变量转为因子并设置水平
  anova_sum = anovaMean(iv,n,ybar,sd) # 基于均值和标准差执行方差分析
  
  # 计算单个对比项的核心函数
  calculate_contrasts <- function(lvl1,lvl2,iv,obs,ybar){
    m <- function(lvl){ybar[iv == lvl] %>% mean()}
    n <- function(lvl){obs[iv == lvl] %>% sum()}
    
    return(list(contrast   = paste(lvl1, "-", lvl2), # 生成对比项名称
                mean_diffs = m(lvl1) - m(lvl2),      # 两组均值差
                n_precalc  = 1/n(lvl1) + 1/n(lvl2)))}# t统计量的前置计算项
    
  # 循环生成所有成对对比项
  contr      = character()
  mean_diffs = numeric()
  n_precalcs = numeric()
    
  for (i in 1:(nlevels(iv)-1)) {
    for (j in (i+1):nlevels(iv)) {
      lvl1 = levels(iv)[i]
      lvl2 = levels(iv)[j]
      res  = calculate_contrasts(lvl1,lvl2,iv,n,ybar)
      contr      = append(contr, res$contrast)
      mean_diffs = append(mean_diffs, res$mean_diffs)
      n_precalcs = append(n_precalcs, res$n_precalc)}}
  
  print(contr)

  # 计算Scheffé检验的各项统计量
  n_total    = sum(n)      # 总观测数（各组样本量之和）
  k_groups   = nlevels(iv) # 分组数/自变量水平数
  rank       = k_groups-1  # Scheffé检验的秩
  mse        = anova_sum$`Mean Sq`[2] # 从方差分析结果提取均方误差
  pooled_n   = mean(n_precalcs)       # 平衡设计下所有值一致
  se         = sqrt(mse*pooled_n)     # 均值的标准误
  t_ratio    = mean_diffs/se          # t比值
  p_values   = pf(t_ratio**2/rank, rank, n_total - k_groups, lower.tail = F)
  t_crit     = sqrt(rank*qf(1-alpha, rank, n_total - k_groups)) # 临界t值
  lwr_ci     = mean_diffs - t_crit * se # 置信区间下限
  upr_ci     = mean_diffs + t_crit * se # 置信区间上限
  
  return(data.frame(contrast  = contr,
                    estimates = mean_diffs,
                    SE        = se,
                    df        = n_total - k_groups,
                    t.ratio   = t_ratio,
                    p.value   = p_values,
                    lwr.ci    = lwr_ci,
                    upr.ci    = upr_ci))}

使用示例 #

scheffe.test(iv    = Data_sum$Location, # 分组变量
             n     = Data_sum$n,        # 各组样本量
             ybar  = Data_sum$mean,     # 各组均值
             sd    = Data_sum$sd,       # 各组标准差
             alpha = 0.05) %>%          # 可自定义alpha水平
  add_significance("p.value") # 给P值添加显著性标记

结果说明 #

该函数输入参数与HH包的anovaMean()一致，默认alpha水平为0.05，运行后返回类似DescTools::sheffeTest()格式的DataFrame，包含置信区间，完全基于汇总统计量计算。

注意： 若要生成与DescTools::sheffeTest()完全一致的对比项，需对原始数据的自变量执行factor(levels = unique())处理。

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