【发布时间】：2020-03-05 16:33:11
【问题描述】：
我正在尝试使用以下数据执行混合模型拟合。
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")
然后我计算了数据中所有四种治疗（CC、CCS、CS 和 SCS）的成对 p 值，得到 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)
)
我希望我的最终结果类似于带有字母符号的最小二乘均值表，以表示如下所示的重要治疗组（注意：我在下面的预期表格中没有收益率数据） .我该怎么做才能得到这种格式的最小二乘均值？提前感谢您的帮助！
【问题讨论】：
-
类似
stargazer::stargazer(LS_MEAN, type="text", summary=FALSE)? -
@jay.sf 我并不是要格式化表格，而是要根据 p_table 使用字母来分隔重要的 LS_MEAN 值。
-
啊哈，你能详细说明一下 a-b-logic 吗？
-
@jay.sf 任何相互显着的治疗组都有不同的字母（字母）。不显着的治疗组将具有相同的字母。比如Protein_U列，在预期结果表中，CC(a)-CCS(ab)不显着，但CC(a)-CS(b)显着。
-
我想鼓励您找到一种不同的方式来显示这些比较。将字母分组的问题在于，它们倾向于强调不“显着”的差异——因此强调非发现而不是发现。你真的不能说两种方法是一样的。
标签： r statistics lme4 mixed-models p-value