嵌套PERMANOVA排列数计算疑惑：permute包结果与文献不符求解

核心问题 #

你当前的代码设置的置换逻辑，和检验因子A效应（统计量MSA/MSB(A)）所需的置换策略完全不符，这是导致排列数与Anderson & ter Braak提到的10种不一致的原因。

在Anderson & ter Braak的场景中：

• 因子A含2个水平，每个A水平下嵌套3个B水平（共6个独立B组）
• 检验A效应时，置换对象是B组所属的A水平：将6个B组随机分配到A的两个水平中，每组3个。由于交换A1和A2的分组会得到完全等价的统计量值，因此有效唯一排列数为 C(6,3)/2 = 20/2 = 10，对应最小p值1/10=0.1。

而你当前的代码是按B分层，对每个B组内的观测值做自由置换，这是检验B(A)效应（即B在A内的变异）的置换策略，和检验A效应的逻辑无关。同时原数据中obs.rep的重复值会导致allPerms自动去重，进一步减少了返回的排列数。

修正后的代码 #

要枚举检验A效应的所有有效排列，需将B组作为置换单元，在A水平间重新分配（每个A水平保持3个B组），代码如下：

library(permute)

# 构造数据：给B组添加A前缀明确嵌套关系，用随机观测值避免重复去重
tmp <- data.frame(
  A = rep(c('A.1','A.2'), each = 9),
  B = rep(paste0(c('A.1_','A.2_'), rep(1:3, each=3)), each=3),
  obs = rnorm(18)
)

# 提取所有唯一B组
unique_B <- unique(tmp$B)
# 生成所有从6个B组选3个分配给A.1的组合
combos <- combn(length(unique_B), 3, simplify = FALSE)

# 生成每种组合对应的排列（B组归属A水平的置换）
permutations <- lapply(combos, function(idx) {
  new_A <- ifelse(tmp$B %in% unique_B[idx], 'A.1', 'A.2')
  data.frame(tmp, new_A = new_A)
})

# 统计有效唯一排列数（去除A水平互换的重复组合）
length(permutations)/2 # 结果为10，与文献描述一致

代码说明 #

1. 用随机观测值替代重复的obs.rep，避免allPerms因数据重复自动去重。
2. 通过组合生成的方式直接枚举所有B组分配到A水平的可能，再除以2得到10种有效排列（每个组合与其补集对应等价的统计量，视为同一排列）。

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