Как получить попарно таблицу значений p из подобранной модели

Question

У меня есть данные dt ниже, и я подобрал модель для этих данных. Я хотел бы получить попарно таблицу значений p для всех возможных пар (для всех уровней) в модели.

data:

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-значения множественного сравнения такие же, как и простое сравнение? (Если вам нужна дополнительная информация, возможно, лучше обратиться к Cross Validated). Я согласен с комментарием Стефана Лорана.

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