1

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

在我的示例数据集中，我有两个条件，“drugA”和“drugB”。有 6 只动物 AF，每只动物的体重在每种药物的影响下测量了 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

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

DrugA, time1 与 DrugB, time1

DrugA, time2 与 DrugB, time2

DrugA, time3 与 DrugB, time3

DrugA，时间 1 与时间 2

DrugA，时间 2 与时间 3

DrugB，时间 1 与时间 2

DrugB，时间 2 与时间 3

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

4

2 回答 2

0

这里的线索是查看 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")
于 2020-01-22T08:23:48.630 回答
0

您可以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
于 2020-06-24T10:42:43.180 回答