我有一个随机区组设计的实验，自变量包括处理、地点、区组和天数，因变量是一个计数变量。我进行了方差分析（ANOVA），结果显示处理因素显著。希望在每个天数对应的取值上得到每个处理相对于对照的p 值。我可以先建立一个线性模型，然后对其运行LSD.test。但输出并未提供对比的p 值。

library(agricolae) # for LSD.test

df = data.frame(days = c(rep(14, 9), rep(21, 9), 
                         rep(14, 9), rep(21, 9), rep(42, 18)), 
                 site = c(rep(2, 18), rep(1, 18), rep(2, 9), rep(1, 9)), 
                 block = rep(c(rep('1', 3), rep('2', 3), rep('3', 3)),6), 
                 trt = rep(c('pink', 'yellow', 'control'), 18), 
                 measurement = c(1, 1, 1, 0, 0, 1, 0, 0, 2, 0, 5, 5, 0, 3, 6, 
                                 0, 0, 7, 1, 0, 0, 0, 5, 5, 0, 2, 4, 5, 0, 19, 
                                 3, 0, 9, 2, 0, 3, 0, 0, 0, 0, 4, 4, 0, 4, 0, 
                                 1, 0, 4, 5, 0, 5, 0, 0, 0))

df$blockUID <- paste0(df$site, "_", df$block)

summary(aov(measurement ~ trt + blockUID + days, data = df))

linmod <- lm(measurement ~ trt + blockUID + days, data = df)

LSD.test(linmod, trt = 'trt')

解决方案

如果你想要每个对比的p 值，在运行检验时需要设置 group=FALSE，然后再提取比较表。

result <- LSD.test(linmod, trt = 'trt', group=FALSE)
result$comparison
#                  difference pvalue signif.        LCL      UCL
# control - pink    3.1666667 0.0031      **  1.1286894 5.204644
# control - yellow  2.8333333 0.0075      **  0.7953561 4.871311
# pink - yellow    -0.3333333 0.7434         -2.3713106 1.704644