r-通过sjPlots tab_model和emmeans对比度函数进行调整后，同一LMM的不同调整p值

Question

我想我错过了一些基本知识，可能忽略了一些重要的东西。。。

背景：我有一个数据集，其中来自4个不同组(1个对照组和3个治疗组(的动物接受了握力测试。每个试验包括约5次测量。为了了解治疗组与对照组的差异，使用lme4的lmer函数进行了线性混合模型，如下所示：

model1<-lmer(griPASTA~group+(1|ID), data=all)

griPASTA是表示握力的变量，组表示治疗，ID是每只动物的ID。

sjPlot::tab_model(model1, p.adjust="none")

创建此：[OHDA 6 vs CTR]p=0.001在此处输入图像描述

包括holm调整

sjPlot::tab_model(model1, p.adjust="holm")

提供了以下内容：[OHDA 6 vs CTR]p=0.004在此处输入图像描述

然而，当我尝试自己指定对比度并使用holm的相同校正时，我得到了非常不同的p值。

例如

m1_emm <- emmeans(model1, specs = c("group"))
CTR <- c(1,0,0,0)
OHDA_12 <- c(0,1,0,0)
OHDA_6 <- c(0,0,1,0)
OHDA_6X <- c(0,0,0,1)
m1_simple<-contrast(m1_emm, method = list("CTR - OHDA 12" = CTR - OHDA_12,
"CTR - OHDA 6"= CTR - OHDA_6,
"CTR - OHDA 6 w/o REB"= CTR - OHDA_6X),
adjust = "holm")%>% summary (infer=T)

[OHDA 6 vs CTR]p=0.0100在此处输入图像描述

为了进行比较，这是我删除更正时得到的：

m1_simple<-contrast(m1_emm, method = list("CTR - OHDA 12" = CTR - OHDA_12,
"CTR - OHDA 6"= CTR - OHDA_6,
"CTR - OHDA 6 w/o REB"= CTR - OHDA_6X),
adjust = "none")%>% summary (infer=T)

[OHDA 6 vs CTR]p=0.0033在此处输入图像描述

因此，使用带有p.adjust的tab_model，CTR和OHDA 6之间比较的holm校正p值为0.004，如果我在emmeans对比度函数中引入校正，则为0.01。

你知道这种差异的来源是什么以及如何处理它吗？我想我错过了什么，而在第一种和第二种方法中，调整实际上考虑了不同的因素？

无论如何，我想我错过了一些基本的统计/R知识，我现在太愚蠢了，无法理解这一点，所以任何帮助都将不胜感激。

Answer 1

我尝试了一个移动部件较少的类似示例：

data(PlantGrowth)
mod <- lm(weight ~ group, data = PlantGrowth)
## sjPlot::tab_model(mod, show.se = TRUE, show.stat = TRUE, p.adjust = "none")
## Yields the same results as...
summary(mod)$coefficients
#>             Estimate Std. Error   t value     Pr(>|t|)
#> (Intercept)    5.032  0.1971284 25.526514 1.936575e-20
#> grouptrt1     -0.371  0.2787816 -1.330791 1.943879e-01
#> grouptrt2      0.494  0.2787816  1.771996 8.768168e-02

如注释中所述，上面的sjPlot::tab_model()调用仅生成模型的summary()中系数的HTML版本。

所以我要得到p值，然后尝试对它们进行Holm调整：

p.unadj <- summary(mod)$coefficients[, 4]
# Unadjusted P values...
round(p.unadj, 4)
#> (Intercept)   grouptrt1   grouptrt2 
#>      0.0000      0.1944      0.0877
# Holm-adjusted P values...
round(p.adjust(p.unadj, "holm"), 4)
#> (Intercept)   grouptrt1   grouptrt2 
#>      0.0000      0.1944      0.1754
# Holm adjustment, excluding the intercept...
round(p.adjust(p.unadj[-1], "holm"), 4)
#> grouptrt1 grouptrt2 
#>    0.1944    0.1754
# Replacing the intercept's P value with a larger one...
round(p.adjust(c(0.15, p.unadj[-1]), "holm"), 4)
#>           grouptrt1 grouptrt2 
#>     0.300     0.300     0.263 

我验证了sjPlot::tab_model(mod, p.adjust = "holm")产生了与上述相同的调整p值(0.0000、0.1944、0.1754(。在这个特定的例子中，当我们排除截距时，调整后的P值是相同的，但对于不同的P值模式，这不一定是真的，如最终结果所示

现在来看emmeans结果：

library(emmeans)
EMM <- emmeans(mod, "group")
CON <- contrast(EMM, "trt.vs.ctrl1")
# Unadjusted P values
test(CON, adjust = "none")
#>  contrast    estimate    SE df t.ratio p.value
#>  trt1 - ctrl   -0.371 0.279 27 -1.331  0.1944 
#>  trt2 - ctrl    0.494 0.279 27  1.772  0.0877
# Holm-adjusted P values
test(CON, adjust = "holm")
#>  contrast    estimate    SE df t.ratio p.value
#>  trt1 - ctrl   -0.371 0.279 27 -1.331  0.1944 
#>  trt2 - ctrl    0.494 0.279 27  1.772  0.1754 
#> 
#> P value adjustment: holm method for 2 tests

^{创建于2021-05-10由reprex包(v1.0.0(}

这些结果与上面使用来自summary()的未调整的p值并且排除Holm调整之前的截距获得的结果相匹配。

我确信OP中提到的差异是由于一些对象是"；"陈腐"；，即，在进行一些改变之前获得，而在进行一些更改之后获得其他结果。我建议重新开始：重新拟合模型，并重新计算所有正在比较的结果。我有时也会遇到同样的事情，因为人们很容易忘记在尝试一大堆其他事情时迈出的一步。

在混合模型的情况下，自由度问题可能会导致一些差异，但我预计这些差异会很小，除非d.f.非常小。