如何从LSD.test中获取各处理与对照组的对比p值？

1. 先修正模型：加入处理与天数的交互项 #

你的原模型只纳入了trt、blockUID和days的主效应，但要分析不同天数下处理组与对照组的差异，必须加入trt:days交互项——否则只能得到处理因素的整体效应，无法拆分到不同天数层面讨论。

# 带交互项的线性模型
linmod_interact <- lm(measurement ~ trt + blockUID + days + trt:days, data = df)
summary(aov(linmod_interact))

另外，因变量是计数变量，线性模型可能违背正态性假设，更推荐用泊松/负二项回归（如果数据存在过度离散）：

library(MASS)
# 负二项回归（适配计数数据的过度离散情况）
pois_mod <- glm.nb(measurement ~ trt + blockUID + days + trt:days, data = df)
summary(pois_mod)

2. 分天数提取处理组与对照组的对比p值 #

这里推荐用emmeans包，它能高效对交互效应做简单检验，直接按days分组，输出每个处理组与对照组的对比p值。

操作步骤： #

# 首次使用需安装包
# install.packages("emmeans")
library(emmeans)

# 基于模型生成按days分组的边际均值
emm <- emmeans(linmod_interact, ~ trt | days)

# 指定以control为参照，输出每个处理组与对照组的对比结果（含p值）
contrast(emm, method = "trt.vs.ctrl", ref = "control")

如果用泊松/负二项模型，只需把linmod_interact替换成对应的模型对象即可，emmeans会自动适配广义线性模型。

3. 用agricolae包分天数做LSD检验 #

如果你更习惯用agricolae的工具，可以按days分组循环处理：

library(agricolae)

# 获取所有天数水平
days_levels <- unique(df$days)

# 循环分天数做LSD检验，提取与对照组的对比p值
for (d in days_levels) {
  cat("=== 天数", d, "的检验结果 ===\n")
  # 筛选对应天数的子集
  df_sub <- subset(df, days == d)
  # 拟合子集模型（区组+处理）
  mod_sub <- lm(measurement ~ trt + blockUID, data = df_sub)
  # 执行LSD检验
  lsd_result <- LSD.test(mod_sub, trt = "trt", p.adj = "none")
  # 提取仅与control对比的结果
  ctrl_comp <- lsd_result$comparison[grepl("control", lsd_result$comparison$trt1), ]
  print(ctrl_comp)
  cat("\n")
}

关键提示 #

• 必须加入trt:days交互项，否则无法得到不同天数下处理组的差异结果。
• 计数变量优先选择广义线性模型（泊松/负二项），线性模型可能因违背假设导致结果偏差。
• emmeans的trt.vs.ctrl方法比分组循环更高效，且能直接输出标准化的对比结果。

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