Как получить разделение букв для обозначения значимых наименьших квадратов средних на основе таблицы значений p

Question

Я пытаюсь выполнить смешанную подборку моделей с моими данными ниже.

df.urbana <- structure(list(Location = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Monmouth", 
"Urbana"), class = "factor"), treatment = structure(c(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), .Label = c("1", 
"2", "3", "4"), class = "factor"), B.glucosidase = c(0.845077, 
1.011463, 0.857032, 0.989803, 0.859022, 0.919467, 1.01717, 0.861689, 
0.972332, 0.952922, 0.804431, 0.742634, 1.195837, 1.267285, 1.08571, 
1.20097), Protein = c(7933.333333, 7000, 6352.982456, 8153.684211, 
6077.894737, 4939.649123, 5002.807018, 6489.122807, 4694.035088, 
5901.052632, 4303.859649, 6768.421053, 6159.298246, 6090.526316, 
4939.649123, 5262.45614), POX.C = c(683.3528, 595.9173, 635.4315, 
672.4234, 847.2944, 745.5665, 778.3548, 735.8141, 395.2647, 570.4148, 
458.0383, 535.3851, 678.0293, 670.7419, 335.2923, 562.5674), 
    yield = c(5583L, 5442L, 5693L, 5739L, 5045L, 4902L, 5006L, 
    5086L, 4639L, 4781L, 4934L, 4857L, 4537L, 4890L, 4842L, 4608L
    )), row.names = 17:32, class = "data.frame")

Затем я вычислил парное значение p для всех четырех обработок (CC, CCS, CS и SCS) в моих данных, чтобы получить таблицу значений p (p-table) следующим образом:

mod.yield.U <- lmerTest::lmer(yield ~ treatment + (1|block),data=df.urbana)
summary(mod.yield.U)
p.value.yield.U <- emmeans::emmeans(mod.yield.U, pairwise ~ treatment)
p.value.yield.U <- (summary(p.value.yield.U, adjust = "none") %>%   # default adjust is tukey
                      purrr::pluck("contrasts") %>% 
                      as.data.frame())[c(1,6)]
colnames(p.value.yield.U) <- c("contrasts.contrast", "Yield_Urbana")

mod.B.glucosidase.U <- lmerTest::lmer(B.glucosidase ~ treatment + (1|block),data=df.urbana)
summary(mod.B.glucosidase.U)
p.value.B.glucosidase.U <- emmeans::emmeans(mod.B.glucosidase.U, pairwise ~ treatment)
p.value.B.glucosidase.U <- (summary(p.value.B.glucosidase.U, adjust = "none") %>%   # default adjust is tukey
                              purrr::pluck("contrasts") %>% 
                              as.data.frame())[c(1,6)]
colnames(p.value.B.glucosidase.U) <- c("contrasts.contrast", "B.glucosidase_Urbana")

mod.Protein.U <- lmerTest::lmer(Protein ~ treatment + (1|block), data=df.urbana)
summary(mod.Protein.U)
p.value.Protein.U <- emmeans::emmeans(mod.Protein.U, pairwise ~ treatment)
p.value.Protein.U <- (summary(p.value.Protein.U, adjust = "none") %>%   # default adjust is tukey
                        purrr::pluck("contrasts") %>% 
                        as.data.frame())[c(1,6)]
colnames(p.value.Protein.U) <- c("contrasts.contrast", "Protein_Urbana")
mod.POX.C.U <- lmerTest::lmer(POX.C ~ treatment + (1|block),data=df.urbana)
summary(mod.POX.C.U)
p.value.POX.C.U <- emmeans::emmeans(mod.POX.C.U, pairwise ~ treatment)
p.value.POX.C.U <- (summary(p.value.POX.C.U, adjust = "none") %>%   # default adjust is tukey
                      purrr::pluck("contrasts") %>% 
                      as.data.frame())[c(1,6)]
colnames(p.value.POX.C.U) <- c("contrasts.contrast", "POX.C_Urbana")

# merge all
p_table <- Reduce(
  function(x, y, ...) merge(x, y, by = "contrasts.contrast", ...),
  list(p.value.yield.U, p.value.Protein.U, p.value.POX.C.U, p.value.B.glucosidase.U)
)

Я также рассчитал среднее значение по методу наименьших квадратов (LS_MEAN таблица) следующим образом:

mod.yield.U <- lmerTest::lmer(yield ~ treatment + (1|block),data=df.urbana)
summary(mod.yield.U)
LS_MEAN.yield.U <- emmeans::emmeans(mod.yield.U, pairwise ~ treatment)
LS_MEAN.yield.U <- as.data.frame(LS_MEAN.yield.U[1])
LS_MEAN.yield.U <- as.data.frame(LS_MEAN.yield.U[c("emmeans.treatment", "emmeans.emmean")])
colnames(LS_MEAN.yield.U) <- c("emmeans.treatment", "Yield_Urbana")

mod.B.glucosidase.U <- lmerTest::lmer(B.glucosidase ~ treatment + (1|block),data=df.urbana)
summary(mod.B.glucosidase.U)
LS_MEAN.B.glucosidase.U <- emmeans::emmeans(mod.B.glucosidase.U, pairwise ~ treatment)
LS_MEAN.B.glucosidase.U <- as.data.frame(LS_MEAN.B.glucosidase.U[1])
LS_MEAN.B.glucosidase.U <- as.data.frame(LS_MEAN.B.glucosidase.U[c("emmeans.treatment", "emmeans.emmean")])
colnames(LS_MEAN.B.glucosidase.U) <- c("emmeans.treatment", "B.glucosidase_Urbana")

mod.Protein.U <- lmerTest::lmer(Protein ~ treatment + (1|block), data=df.urbana)
summary(mod.Protein.U)
LS_MEAN.Protein.U <- emmeans::emmeans(mod.Protein.U, pairwise ~ treatment)
LS_MEAN.Protein.U <- as.data.frame(LS_MEAN.Protein.U[1])
LS_MEAN.Protein.U <- as.data.frame(LS_MEAN.Protein.U[c("emmeans.treatment", "emmeans.emmean")])
colnames(LS_MEAN.Protein.U) <- c("emmeans.treatment", "Protein_Urbana")

mod.POX.C.U <- lmerTest::lmer(POX.C ~ treatment + (1|block),data=df.urbana)
summary(mod.POX.C.U)
LS_MEAN.POX.C.U <- emmeans::emmeans(mod.POX.C.U, pairwise ~ treatment)
LS_MEAN.POX.C.U <- as.data.frame(LS_MEAN.POX.C.U[1])
LS_MEAN.POlX.C.U <- as.data.frame(LS_MEAN.POX.C.U[c("emmeans.treatment", "emmeans.emmean")])
colnames(LS_MEAN.POX.C.U) <- c("emmeans.treatment", "POX.C_Urbana")

# merge all
LS_MEAN <- Reduce(
  function(x, y, ...) merge(x, y, by = "emmeans.treatment", ...),
  list(LS_MEAN.yield.U, LS_MEAN.Protein.U, LS_MEAN.POX.C.U, LS_MEAN.B.glucosidase.U)
)

Что бы я хотел получить в качестве конечного результата, это что-то вроде таблицы средних по наименьших квадратовс буквенными обозначениями для обозначения значимых групп лечения, как показано ниже (Примечание: у меня нет данных по урожайности в ожидаемой таблице ниже). Что я могу сделать, чтобы получить значение наименьших квадратов в этом формате? Заранее спасибо за помощь!

Answer 1

Это, наверное, веселее, если вы автоматизируете процесс. Вы можете использовать lapply() для обхода различных зависимых переменных.

Поместить имена зависимых переменных в вектор Y и создать базовую формулу fo. В каждой итерации используйте update() для изменения зависимой переменной в формуле. lapply() будет циклически проходить четыре y с, и в результате будет получен список Mod.U.

Y <- c("yield", "Protein", "POX.C", "B.glucosidase")
fo <- y ~ treatment + (1|block)
Mod.U <- lapply(Y, function(y) lmerTest::lmer(update(fo, paste(y, "~ .")), data=df.urbana))

(Обратите внимание, что приходит сообщение boundary (singular) fit: see ?isSingular , и вы можете иметь единственное числоfit куда-нибудь.)

Расчет p-значения и среднего значения может быть выполнен одинаково. lapply() теперь перебирает элементы списка Mod.U. setNames() создает имена ваших столбцов. (Имена немного скрыты, но str(Mod.U) показывает, что их можно найти в names(mod.U@frame)[1]).

P.value.U <- lapply(Mod.U, function(mod.U)
  (emmeans::emmeans(mod.U, pairwise ~ treatment) %>% 
     summary(adjust = "none") %>%
     purrr::pluck("contrasts") %>% 
     as.data.frame)[c(1, 6)] %>%
    setNames(c("contrasts.contrast", paste0(names(mod.U@frame)[1], "_urbana"))))

Теперь также легко набрать Reduce(), поскольку у нас уже есть список.

p_table <- Reduce(function(x, y, ...) merge(x, y, by = "contrasts.contrast", ...), P.value.U)
#   contrasts.contrast yield_urbana Protein_urbana POX.C_urbana B.glucosidase_urbana
# 1           CC - CCS 7.977935e-05    0.005002692  0.084083738         0.8581554114
# 2            CC - CS 7.676418e-06    0.002530699  0.043628964         0.3799550686
# 3           CC - SCS 3.416131e-06    0.004771013  0.235380612         0.0023712583
# 4           CCS - CS 4.505872e-02    0.664550854  0.002028628         0.4785341647
# 5          CCS - SCS 9.777084e-03    0.976231386  0.010602656         0.0018094630
# 6           CS - SCS 3.726498e-01    0.685994660  0.310848786         0.0006411248

Мы поступаем аналогичным образом со средствами.

LS_MEAN.U <- lapply(Mod.U, function(mod.U) 
  (emmeans::emmeans(mod.U, pairwise ~ treatment)[1] %>%
     as.data.frame)[c("emmeans.treatment", "emmeans.emmean")] %>%
    setNames(c("emmeans.treatment", paste0(names(mod.U@frame)[1], "_urbana"))))

LS_MEAN <- Reduce(function(x, y, ...) merge(x, y, by = "emmeans.treatment", ...), LS_MEAN.U)
#   emmeans.treatment yield_urbana Protein_urbana POX.C_urbana B.glucosidase_urbana
# 1                CC      5614.25       7360.000     646.7812            0.9258438
# 2               CCS      5009.75       5627.368     776.7575            0.9143370
# 3                CS      4802.75       5416.842     489.7757            0.8680797
# 4               SCS      4719.25       5612.982     561.6577            1.1874505

Теперь я все еще не уверен, как таблица ожидаемых результатов (назовем ее out) относится к p_table таблица, но вы можете использовать paste0() для создания суффиксов по желанию. formatC() обрезает до нужного количества цифр.

out <- LS_MEAN[-2]
out$Protein_urbana <- paste0(formatC(LS_MEAN$Protein_urbana, format="f", digits=1), 
                             c("a", "ab", "b", "ab"))
out$POX.C_urbana <- paste0(formatC(LS_MEAN$POX.C_urbana, format="f", digits=2), 
                           c("ab", "a", "b", "b"))
out$B.glucosidase_urbana <- paste0(formatC(LS_MEAN$B.glucosidase_urbana, format="f", digits=3), 
                                   c("b", "b", "b", "a"))

names(out) <- c("Treatment", "Protein_U", "POX.C_U", "B.glucosidase_U")
out
#   Treatment Protein_U  POX.C_U B.glucosidase_U
# 1        CC   7360.0a 646.78ab          0.926b
# 2       CCS  5627.4ab  776.76a          0.914b
# 3        CS   5416.8b  489.78b          0.868b
# 4       SCS  5613.0ab  561.66b          1.187a