使用how()进行重复测量时Adonis2阻塞错误的排查与代码验证

我开展了一项试验设计：在5个区块（B1-B5）中种植3种不同植物（S1-S3），设置2种处理（T1、T2），并在2个时间点（Y1-Y2）采集土壤微生物样本。目标是探究植物×处理×年份对微生物群落组成的影响（需分析所有三向及双向交互作用）。根据试验设计，区块（block）和样方（plot，即区块×植物×处理组合，作为重复测量的受试对象）均需设为随机效应。

调研后采用permute包构建模型，但Adonis2运行出现问题，以下是代码及三种尝试：

library(vegan)
library(tidyr)
library(dplyr)
library(tibble) #i think these would do.. if not, I apologize!

meta.raw<- as.factor(c(1:60))

meta.raw %>%
  as.data.frame() %>%
  `colnames<-`("sample") %>%
  mutate(block = rep(paste0("B", sample(1:5)), 2, each = 6),
         crop = rep(c("S1","S2","S3"), 10, each = 2),
         trt = rep (c("T1", "T2"), 30),
         year = rep(c("Y1", "Y2"), 2, each = 15),
         plot = paste0(block,crop,trt)) -> meta
  
matrix(round(runif(n=6000, min=0, max=2000), 0), nrow=60) %>%
  as.data.frame() %>%
  mutate(sample = as.factor(c(1:60))) %>%
  column_to_rownames("sample")-> df

# how1: 参考permute文档语法的写法

h1 <- how(within = Within(type = "series"),
         plots = Plots(strata = meta$plot),
         blocks = Blocks(strata = meta$block),
         nperm = 499)

adonis2(df ~ crop*trt*year,
        data = meta,
        dist = "bray",
        perm = h1,
        by = "margin")
# 报错：
# Error in check(sn, control = control, quietly = quietly) : 
#   Number of observations and length of Block 'strata' do not match.

# how 2: 在how()中设置within和plots，在adonis2()的strata参数中传入区块变量

h2 <- how(within = Within(type = "series"),
         plots = Plots(strata = meta$plot),
         nperm = 499) # 这里不设置blocks

adonis2(df ~ crop*trt*year,
        data = meta,
        dist = "bray",
        perm = h2,
        strata = meta$block, # 这里传入block
        by = "margin")
# 运行成功，但不确定是否正确：
# Permutation test for adonis under reduced model
# Marginal effects of terms
# Blocks:  strata 
# Plots: meta$plot, plot permutation: none
# Permutation: series
# Number of permutations: 499
# 
# adonis2(formula = df ~ crop * trt * year, data = meta, permutations = h2, by = "margin", strata = meta$block, dist = "bray")
#               Df SumOfSqs      R2      F Pr(>F)
# crop:trt:year  2   0.1027 0.03085 0.9117      1
# Residual      48   2.7043 0.81219              
# Total         59   3.3296 1.00000  

#how 3: 在how()中传入区块变量，但写法和permute文档不同

h3 <- how(within = Within(type = "series"),
         plots = Plots(strata = meta$plot),
         blocks = meta$block,
         nperm = 499)

adonis2(df ~ crop*trt*year,
        data = meta,
        dist = "bray",
        perm = h3,
        by = "margin")
# 也运行成功，但不确定是否正确：
# Permutation test for adonis under reduced model
# Marginal effects of terms
# Blocks:  meta$block 
# Plots: meta$plot, plot permutation: none
# Permutation: series
# Number of permutations: 499
# 
# adonis2(formula = df ~ crop * trt * year, data = meta, permutations = h3, by = "margin", dist = "bray")
#               Df SumOfSqs      R2      F Pr(>F)
# crop:trt:year  2   0.1007 0.02936 0.8709      1
# Residual      48   2.7760 0.80897              
# Total         59   3.4316 1.00000

尝试1运行报错，尝试2和3可运行，但不确定统计方法是否正确。数据是平衡设计，却出现类似非平衡数据的错误，想知道哪里操作有误，以及正确的实现方式。

1. 尝试1报错原因 #

报错Number of observations and length of Block 'strata' do not match是因为Blocks(strata = meta$block)的用法错误：

• Blocks()的strata参数需要的是分组后的层级标识（而非样本级向量），直接传入meta$block会让permute包误将每个样本视为独立的block层级，和实际观测数不匹配。你的试验中每个区块对应12个样本，正确的做法是直接用样本级的分组向量赋值给blocks参数（无需嵌套Blocks()函数）。

2. 尝试2和3的正确性判断 #

两种方法本质等价且均正确，核心逻辑都是符合你的试验设计的：

• 你的试验是重复测量设计：每个plot（区块×植物×处理）在两个年份被重复采样，Within(type = "series")确保同一plot内的样本不被跨时间点置换。
• 区块是更高层级的随机效应，两种方法都限制了置换仅在区块内部进行，避免跨区块的样本混淆。

3. 推荐的正确实现方式 #

使用更清晰的permute语法构建置换方案，同时确保模型覆盖所有交互项：

# 构建符合试验设计的置换方案
h_correct <- how(
  blocks = meta$block,  # 最高层级分组：区块，置换仅在区块内进行
  plots = Plots(strata = meta$plot),  # 重复测量单位：样方，保证同一样方的样本不被拆分
  within = Within(type = "series"),  # 重复测量置换规则：固定同一样方内的时间点顺序
  nperm = 499
)

# 运行包含所有交互项的adonis2模型（边际效应对应Type III平方和，适合交互项分析）
result <- adonis2(
  df ~ crop * trt * year,
  data = meta,
  dist = "bray",
  permutations = h_correct,
  by = "margin"
)

print(result)

关键说明 #

• blocks = meta$block：明确区块为随机效应，置换被限制在区块内部，符合试验的随机化逻辑。
• Plots(strata = meta$plot)：指定样方为重复测量的基本单位，确保每个样方的两个年份样本不会被置换到其他样方。
• Within(type = "series")：针对重复测量数据的标准置换策略，避免破坏同一受试对象的时间序列关联。
• by = "margin"：计算边际效应（Type III平方和），能准确评估每个交互项的独立贡献，满足你分析所有三向、双向交互作用的需求。

4. 额外注意事项 #

• 模拟数据是随机生成的，因此结果中P值为1是正常现象，实际数据会呈现真实效应。
• 需确保meta数据框中的plot变量准确对应每个重复测量单位，即每个plot恰好包含2个样本（Y1和Y2）。

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