如何提取 LSD.test 的 p 值并按天比较处理组与对照组？

2026-11-04数据小许

将 LSD.test() 的参数设为 group = FALSE，就可以从返回对象的 $comparison 中读取 pvalue。当前模型没有处理与天数的交互项，因此无法据此获得随天数变化的处理效应。如果要分别比较第 14、21、42 天的处理组与对照组，可以将 days 设为因子，拟合包含 trt * days 的模型，再用 emmeans 构造各天的对照比较。

提取 LSD.test 的 p 值

沿用问题中的 linmod：

library(agricolae)

lsd_result <- LSD.test(
  linmod,
  trt = "trt",
  group = FALSE,
  p.adj = "none"
)

lsd_result$comparison
lsd_result$comparison$pvalue

默认的 group = TRUE 返回显著性字母分组，此时 $comparison 为 NULL。设为 FALSE 后，返回的比较表包含 difference、pvalue 等列。如果希望调用时直接打印结果，还可以加上 console = TRUE。LSD.test 源码

要只保留涉及 control 的比较，可以这样筛选：

comparisons <- lsd_result$comparison

control_comparisons <- comparisons[
  grepl("control", rownames(comparisons), fixed = TRUE),
  ,
  drop = FALSE
]

control_comparisons

这种筛选方式适用于示例中的组名。p.adj = "none" 表示不进行多重比较校正。这段代码给出的是跨天汇总的处理比较，没有按天拆分。

为什么原模型不能给出不同天数的处理效应？

原公式是：

measurement ~ trt + blockUID + days

这个公式有两个限制：

days 是数值变量，模型假定其效应为线性趋势。
没有 trt:days 交互项，模型假定处理组之间的差异在各天相同。

因此，即使计算各天的预测均值，原模型中同一处理与对照的差值也不会随天数变化。如果要分别估计三个观测日的处理效应，可以将 days 设为因子，并加入交互项。

按天比较 pink、yellow 与 control

下面的代码接在问题中的 df 定义之后运行。示例沿用线性模型，并将 blockUID 作为固定区组效应。这些假设是否合适，需要结合实验设计和残差诊断判断。

library(emmeans)

dat <- df
dat$days <- factor(dat$days, levels = c(14, 21, 42))
dat$trt <- relevel(factor(dat$trt), ref = "control")
dat$blockUID <- factor(dat$blockUID)

fit <- lm(
  measurement ~ blockUID + trt * days,
  data = dat
)

emm <- emmeans(fit, ~ trt | days)

ctrl_contrasts <- contrast(
  emm,
  method = "trt.vs.ctrl",
  ref = 1
)

result <- as.data.frame(
  summary(
    ctrl_contrasts,
    infer = c(TRUE, TRUE),
    adjust = "holm"
  )
)

result[, c("days", "contrast", "estimate", "SE", "df", "p.value")]

relevel() 将 control 设为第一个水平，因此 ref = 1 指定它为对照。~ trt | days 表示在每个天数内分别比较处理组；"trt.vs.ctrl" 只构造处理组与对照组的比较。emmeans 使用指南

这里有三个天数和两个非对照处理，共需进行六项比较。estimate 表示处理组减去对照组的估计均值差，p.value 是对应检验的 p 值。正差值表示处理组的预测平均计数更高。

明确多重比较校正的范围

上面的代码使用 adjust = "holm"，并保留 days 分组，因此只在每一天的两个比较内分别校正，没有同时覆盖全部六项比较。

如果推断目标要求将所有天数的六项比较视为同一个检验族，可以取消分组后统一校正：

result_all_days <- as.data.frame(
  summary(
    ctrl_contrasts,
    by = NULL,
    infer = c(TRUE, TRUE),
    adjust = "holm"
  )
)

result_all_days

如果研究方案明确要求未校正的比较，可以改为 adjust = "none"，并在报告结果时说明。by 会影响多重比较校正的检验族范围。emmeans 结果汇总文档

计数数据和实验设计还需要检查什么？

measurement 是计数变量，示例中也包含不少零值。要判断是否适合使用 lm()，需要检查残差分布、方差齐性和观测独立性。若这些假设明显不合适，可以根据计数分布考虑 Poisson 或负二项模型。仅凭现有描述，无法确定最终应采用哪一种。

如果同一实验单元在多个天数被重复测量，模型还需要处理单元内相关性。仅加入固定的 blockUID 效应，并不等于已经处理了重复测量。确定模型时，应说明实验单元的标识、区组是否作为随机效应，以及地点与区组的设计关系。

ANOVA 中处理项显著，也不意味着每个天数下的每个处理与对照比较都会显著。应以合适模型下的具体 contrasts、效应估计和 p 值为依据。