emmeans 中的计划对比答案

【问题标题】：Planned contrasts in emmeansemmeans 中的计划对比
【发布时间】：2020-05-05 20:31:58
【问题描述】：

我想使用 emmeans 计算计划对比度的特定子集，但在编码时遇到了麻烦。

在我的示例数据集中，我有两个条件，“drugA”和“drugB”。共6只动物A-F，每只动物在每种药物作用下称重3次。

id <- rep(c("A","B","C","D","E","F"),6)
drug <- c(rep(c("drugA"), 18), rep(c("drugB"), 18))
time <- rep(rep(1:3, each = 6),2)
value <- c(rnorm(6, 1, 0.4), rnorm(6, 3, 0.5), rnorm(6, 6, 0.8), rnorm(6, 1.1, 0.4), rnorm(6, 0.8, 0.2), rnorm(6, 1, 0.6))
df <- data.frame(id,drug, time, value)

df$id <- as.factor(df$id) 
df$drug <- as.factor(df$drug)
df$time <- as.factor(df$time)
stats <- lmer(value ~ drug*time + drug + time + (1|id), data = df)
summary(stats)

emm <- emmeans(stats, list(pairwise ~ drug + time), adjust = "tukey") 
emm

但是，我只想计算以下对比：

药物 A，时间 1 与药物 B，时间 1

药物 A，时间 2 与药物 B，时间 2

DrugA，time3 与 DrugB，time3

DrugA，时间 1 与时间 2

DrugA，时间 2 与时间 3

DrugB，时间 1 与时间 2

DrugB，时间 2 与时间 3

我必须如何对这些对比进行编码？非常感谢您的建议。

【问题讨论】：

您是否尝试过写出代表每种药物：时间组合的组均值的 0 和 1 向量？或者那是你被困的地方？您将这些向量基于 emmeans 的输出。我会在没有“pairwise”的情况下制作 emm 并从那里开始构建我的向量。
感谢您的建议。是的，我认为挑战是从输出（网格有 3 列和 16 行，第三列是什么？）到对比...
看来你取得了不错的进展！是的，我看到了困难。我认为计算加/减 1 或对另一个因素求平均值可能很困难，这就是为什么我教学生制作一个表示每个组合平均值的向量，然后用这些向量进行数学运算以表示他们想要的比较。 :)
嗯...我仍然无法解决以下问题：估计值、SE 和 p 值有时相同似乎令人费解。有什么问题吗？
由于您不允许模型中的交互，因此对我来说，时间 1 的 A 与 B 的差异与时间 2 的差异是有意义的。您的模型表示无论时间长短，A 和 B 之间存在一个总体差异。

标签： r lme4 emmeans

【解决方案1】：

这里的线索是看emmeans的结果网格。前两列是形成对比的基础，每一行代表一个因素的组合。

emm <- emmeans(stats, list(~ drug + time)) # not used afterwards, but to check result Grid

con <- list(
  DrugA1_DrugB1 = c(1,-1,0,0,0,0),
  DrugA2_DrugB2 = c(0,0,1,-1,0,0),
  DrugA3_DrugB3 = c(0,0,0,0,-1,1),

  DrugA1_DrugA2 = c(1,0,-1,0,0,0),
  DrugA2_DrugA3 = c(0,0,1,0,-1,0),

  DrugB1_DrugB2 = c(0,1,0,-1,0,0),
  DrugB2_DrugB3 = c(0,0,0,1,0,-1)
)

接下来的内容正是为您提供了这些对比。

 emm <- emmeans(stats, list(~ drug + time), contr = con, adjust = "mvt")

【讨论】：

另一个问题是我想使用 multcomp 包进行更强大的 p-val 调整。 emm 对比错误。 emmGrid(X[[i]], ...) : 对比度系数的不一致数
在这个答案中，``adjust = “tukey”` 将被忽略，因为除了成对比较之外，该调整是不合适的。使用adjust = “mvt” 获得与multcomp 使用的单步调整完全相同的调整。您也可以使用as.glht() 让 multcomp 直接执行此操作。请参阅文档。

【解决方案2】：

您可以为此使用contrast()：

library(lme4)
library(emmeans)

id <- rep(c("A","B","C","D","E","F"),6)
drug <- c(rep(c("drugA"), 18), rep(c("drugB"), 18))
time <- rep(rep(1:3, each = 6),2)
value <- c(rnorm(6, 1, 0.4), rnorm(6, 3, 0.5), rnorm(6, 6, 0.8), rnorm(6, 1.1, 0.4), rnorm(6, 0.8, 0.2), rnorm(6, 1, 0.6))
df <- data.frame(id,drug, time, value)

df$id <- as.factor(df$id) 
df$drug <- as.factor(df$drug)
df$time <- as.factor(df$time)
stats <- lmer(value ~ drug*time + drug + time + (1|id), data = df)

emm <- emmeans(stats, ~ drug + time)
emm
#>  drug  time emmean    SE   df lower.CL upper.CL
#>  drugA 1      1.16 0.187 28.8    0.778     1.55
#>  drugB 1      1.05 0.187 28.8    0.666     1.43
#>  drugA 2      3.30 0.187 28.8    2.917     3.68
#>  drugB 2      0.84 0.187 28.8    0.457     1.22
#>  drugA 3      5.99 0.187 28.8    5.602     6.37
#>  drugB 3      1.30 0.187 28.8    0.920     1.69
#> 
#> Degrees-of-freedom method: kenward-roger 
#> Confidence level used: 0.95

contrast(emm, method = "pairwise", by = "time")
#> time = 1:
#>  contrast      estimate    SE df t.ratio p.value
#>  drugA - drugB    0.113 0.253 25  0.446  0.6593 
#> 
#> time = 2:
#>  contrast      estimate    SE df t.ratio p.value
#>  drugA - drugB    2.461 0.253 25  9.745  <.0001 
#> 
#> time = 3:
#>  contrast      estimate    SE df t.ratio p.value
#>  drugA - drugB    4.683 0.253 25 18.545  <.0001 
#> 
#> Degrees-of-freedom method: kenward-roger

contrast(emm, method = "consec", by = "drug")
#> drug = drugA:
#>  contrast estimate    SE df t.ratio p.value
#>  2 - 1       2.139 0.253 25  8.471  <.0001 
#>  3 - 2       2.685 0.253 25 10.634  <.0001 
#> 
#> drug = drugB:
#>  contrast estimate    SE df t.ratio p.value
#>  2 - 1      -0.209 0.253 25 -0.828  0.6244 
#>  3 - 2       0.463 0.253 25  1.834  0.1383 
#> 
#> Degrees-of-freedom method: kenward-roger 
#> P value adjustment: mvt method for 2 tests

【讨论】：