R中pwr包执行HLM功效分析报错：非数值参数传入数学函数

错误原因 #

你遇到的Error in ceiling(result$g): non-numeric argument to a mathematical function，核心原因是**pwr.anova.test的返回对象中不存在g这个属性**，调用result$g会得到NULL，自然无法用ceiling()处理。此外代码还有两个关键问题：

• 错误设置k=4：你的研究是干预/对照2组比较，应设k=2
• 未考虑HLM的ICC：普通ANOVA功效分析不适合分层数据，必须纳入集群内相关系数(ICC)的方差膨胀效应

修正后的代码（适配2-level HLM场景） #

library(pwr)

# 定义核心参数
n_per_clin_group <- 9  # 每个临床组的护理学生数
icc <- 0.1             # 集群内相关系数
target_power <- 0.8    # 目标检验效能
sig_level <- 0.05      # 显著性水平
effect_sizes <- c(0.1, 0.2, 0.3)

# 计算设计效应：集群抽样导致的方差膨胀系数
design_effect <- 1 + (n_per_clin_group - 1) * icc

# 定义HLM功效分析函数
calculate_hlm_power <- function(effect_size) {
  # 2组比较下，将Cohen's f转换为Cohen's d（适配pwr.t.test的参数要求）
  cohen_d <- 2 * effect_size / sqrt(1 - effect_size^2)
  
  # 计算独立样本t检验所需的总个体数（无集群效应的理想情况）
  t_test_result <- pwr.t.test(d = cohen_d, power = target_power, sig.level = sig_level, type = "two.sample")
  total_independent_n <- ceiling(t_test_result$n * 2)
  
  # 调整为集群设计下的总样本量（考虑方差膨胀）
  total_cluster_n <- ceiling(total_independent_n * design_effect)
  
  # 计算所需的最少临床组数
  num_clin_groups <- ceiling(total_cluster_n / n_per_clin_group)
  
  # 输出结果
  cat(sprintf("效应量f=%.1f:\n", effect_size))
  cat(sprintf("  达到80%功效所需最少临床组数: %d\n", num_clin_groups))
  cat(sprintf("  对应总样本量: %d\n", num_clin_groups * n_per_clin_group))
  cat(sprintf("  设计效应（方差膨胀系数）: %.2f\n\n", design_effect))
}

# 遍历所有效应量计算
lapply(effect_sizes, calculate_hlm_power)

代码说明 #

1. 设计效应：用来量化集群抽样导致的方差膨胀，ICC越大，需要的样本量越多，公式为DE = 1 + (每组个体数-1)*ICC
2. 效应量转换：因为pwr.t.test使用Cohen's d，而你提供的是Cohen's f，通过2组比较下的转换公式完成适配
3. 样本量调整：先计算无集群效应的理想样本量，再乘以设计效应得到分层模型下的实际需求，最后换算成临床组数

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