如何从拟合模型中获取成对p值表

Question

我下面有一个数据dt，并且我已经为此数据拟合了模型。我想获取模型中所有可能的对（所有级别）的成对p值表。

数据：

dt <- structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L), .Label = c("CC", "CCS", 
"CS", "SCS"), class = "factor"), block = structure(c(1L, 2L, 
3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 
3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("1", 
"2", "3", "4"), class = "factor"), yield = c(5156L, 5157L, 5551L, 
5156L, 4804L, 4720L, 4757L, 5021L, 4826L, 4807L, 4475L, 4596L, 
4669L, 4588L, 4542L, 4592L, 5583L, 5442L, 5693L, 5739L, 5045L, 
4902L, 5006L, 5086L, 4639L, 4781L, 4934L, 4857L, 4537L, 4890L, 
4842L, 4608L)), row.names = c(NA, -32L), class = "data.frame")

fit <- lmerTest::lmer(yield ~ treatment + (1|block), data = dt)
summary(fit)$coeff

预期结果将是（全部填充）：

               treatmentCC  treatmentCCS    treatmentCS treatmentSCS
treatmentCC                 2.44E-06         1.25E-08    1.34E-09
treatmentCCS    2.44E-06            
treatmentCS     1.25E-08            
treatmentSCS    1.34E-09

我可以重新调整因子并多次拟合模型，并手动填写p值表，但是原始数据中有很多层次。有简单的方法吗？

Answer 1

为什么您认为多重比较的p值与简单比较相同？ （如果您需要更多信息，也许交叉验证是一个更好的网站）。我同意斯特凡·洛朗的评论。

fit <- lmerTest::lmer(yield ~ treatment + (1|block), data = dt)
emmeans_mod <- emmeans::emmeans(fit, pairwise ~ treatment)

## the results of lmerTest

summary(fit) %>% 
  pluck("coefficients") %>% 
  as.data.frame()
#              Estimate Std. Error df   t value     Pr(>|t|)
# treatmentCCS -517.000   87.73629 28 -5.892658 2.443788e-06


## the result of summary method of emmeans

summary(emmeans_mod, adjust = "none") %>%   # default adjust is tukey
  purrr::pluck("contrasts") %>% 
  as.data.frame()
#    contrast estimate       SE df   t.ratio      p.value
# 1  CC - CCS  517.000 87.73629 25 5.8926583 3.783457e-06


pt(-5.892658, df = 28, lower.tail = TRUE) * 2  # 2.44379e-06
pt(-5.892658, df = 28 - 3, lower.tail = TRUE) * 2  # 3.78346e-06