CODEWITHSUNDEEP
rlme4emmeans
Christine Blume
Christine Blume

Reputation: 53

Planned contrasts in emmeans

I would like to compute a specific subset of planned contrasts using emmeans, but have trouble coding these.

In my sample dataset, I have two conditions, "drugA" and "drugB". There are 6 animals A-F and the weight of each animal has been measured 3 times under the influence of each drug.

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

However, I would only like to calculate the following contrasts:

DrugA, time1 vs. DrugB, time1

DrugA, time2 vs. DrugB, time2

DrugA, time3 vs. DrugB, time3

DrugA, time1 vs. time2

DrugA, time2 vs. time3

DrugB, time1 vs. time2

DrugB, time2 vs. time3

How do I have to code these contrasts? Thank you very much for your suggestions.

Upvotes: 2

Views: 2232

Answers (2)

MSB
MSB

Reputation: 136

You can use contrast() for that:

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

Upvotes: 1

Christine Blume
Christine Blume

Reputation: 53

The clue here is to look at the result grid of emmeans. The two first columns are the basis for the formulation of contrasts, each line represents one combination of factors.

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

The following then gives you exactly those contrasts.

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

Upvotes: 1

Related Questions