emmeans 中的计划对比

Question

我想使用 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

药物 A，时间 3 对比药物 B，时间 3

药物 A，时间 1 与时间 2

药物 A，时间 2 与时间 3

药物 B，时间 1 与时间 2

药物 B，时间 2 与时间 3

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

Answer 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")

Answer 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

emmeans 中的计划对比

Planned contrasts in emmeans

r

lme4

emmeans