I'm trying to exclude contrasts from simple pairs.
Made up data:
library(emmeans)
library(lme4)
set.seed(1234)
dat <- data.frame(
dv = c(
# Add in t1 ctrl
rlnorm(meanlog=0.8, sdlog=0.1, n=1000), #Long Healthy
rlnorm(meanlog=0.7, sdlog=0.1, n=1000), #Lat Healthy
rlnorm(meanlog=0.6, sdlog=0.1, n=500), #Long Damaged
rlnorm(meanlog=0.5, sdlog=0.1, n=500), #Lat Damaged
# Add in t2 ctrl
rlnorm(meanlog=0.7, sdlog=0.1, n=1000), #Long Healthy
rlnorm(meanlog=0.6, sdlog=0.1, n=1000), #Lat Healthy
rlnorm(meanlog=0.5, sdlog=0.1, n=500), #Long Damaged
rlnorm(meanlog=0.4, sdlog=0.1, n=500), #Lat Damaged
# Add in t1 trt
rlnorm(meanlog=0.8, sdlog=0.1, n=1000), #Long Healthy
rlnorm(meanlog=0.7, sdlog=0.1, n=1000), #Lat Healthy
rlnorm(meanlog=0.6, sdlog=0.15, n=500), #Long Damaged
rlnorm(meanlog=0.5, sdlog=0.15, n=500), #Lat Damaged
# Add in t2 trt
rlnorm(meanlog=0.7, sdlog=0.1, n=1000), #Long Healthy
rlnorm(meanlog=0.6, sdlog=0.1, n=1000), #Lat Healthy
rlnorm(meanlog=0.65, sdlog=0.15, n=500),#Long Damaged
rlnorm(meanlog=0.55, sdlog=0.15, n=500) #Lat Damaged
),
id=c(rep(c("subj_1", "subj_2"), times=c(6000, 6000))),
intervention=c(rep(c("ctrl", "trt"), times=c(6000, 6000))),
timepoint=c(rep(rep(c("t1", "t2"), times=c(3000, 3000)),2)),
direction=c(rep(rep(c("long", "lat", "long", "lat"), times=c(1000, 1000, 500, 500)),4)),
region=c(rep(rep(c("healthy", "damaged"), times=c(2000, 1000)),4))
)
Model:
lmm_1 <- lmer(log(dv) ~ intervention*timepoint*region + direction + (1|id), data=dat)
emmeans reference grid:
lmm_1_emm <- emmeans(lmm_1,
pairwise ~ intervention*region*timepoint,
type = "response",
bias.adj=TRUE,
sigma=sigma(lmm_1))
Simple pairs:
pairs(regrid(lmm_1_emm,
bias.adjust = TRUE,
sigma = sigma(lmm_1)),
simple = "each",
combine = TRUE)
I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
See '? emm_list' for more information
region timepoint intervention contrast estimate SE df z.ratio p.value
damaged t1 . ctrl - trt -0.00263 0.06162 Inf -0.043 1.0000
healthy t1 . ctrl - trt -0.00554 0.07468 Inf -0.074 1.0000
damaged t2 . ctrl - trt -0.26043 0.06057 Inf -4.300 0.0002
healthy t2 . ctrl - trt -0.00388 0.06743 Inf -0.057 1.0000
. t1 ctrl damaged - healthy -0.37968 0.01228 Inf -30.928 <.0001
. t1 trt damaged - healthy -0.38259 0.01234 Inf -31.001 <.0001
. t2 ctrl damaged - healthy -0.33741 0.01099 Inf -30.691 <.0001
. t2 trt damaged - healthy -0.08085 0.00814 Inf -9.932 <.0001
damaged . ctrl t1 - t2 0.16378 0.00907 Inf 18.057 <.0001
damaged . trt t1 - t2 -0.09403 0.00906 Inf -10.382 <.0001
healthy . ctrl t1 - t2 0.20605 0.00863 Inf 23.868 <.0001
healthy . trt t1 - t2 0.20771 0.00867 Inf 23.957 <.0001
Results are averaged over some or all of the levels of: direction
Degrees-of-freedom method: inherited from asymptotic when re-gridding
P value adjustment: bonferroni method for 12 tests
How would you exclude the 4 "damaged - healthy" contrasts? I could compute the simple effects separately and leave out "region", but I like that the 'simple = "each" & combine = TRUE' options roll in the p-value correction for all the tests.
I referred to the documentation and attempted to just pull out one contrast with the following code but it failed:
contrast(lmm_1_emm, method = list("ctrl damaged t2" - "trt damaged t2") )
Error in "ctrl damaged t2" - "trt damaged t2" :
non-numeric argument to binary operator
p.s: What does "I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong. See '? emm_list' for more information" mean?
First, for custom contrasts, it is always best to not have a left-hand side in the formula in emmeans(). Just get the means you want, then do the contrasts separately, e.g.,
emm <- emmeans(lmm_1,
~ intervention * region * timepoint,
type = "response",
bias.adj = TRUE,
sigma = sigma(lmm_1))
Second, I don't see how contrast(lmm_1_emm, method = list("ctrl damaged t2" - "trt damaged t2") ) could have been suggested by anything in the documentation. The method argument needs to be a named list or the name of a contrast family.
We can see the means, and also the order in which they are presented, just by listing the object, e.g.
emm
In this example there will be 8 means listed. Then you can pick the ones you want to compare by putting 1 and -1 in the positions in the order of listing and 0 elsewhere, e.g.,
contrast(emm, list(
`t d 1 - c h 2` = c(0, 1, -1, 0, 0, 0, 0, 0),
`t h 1 - t d 1` = c(0, -1, 0, 1, 0, 0, 0, 0)) )
This can be rather tedious, so when possible it is better to use by variables and named contrasts families where possible. For example, to compare successive intervention x timepoint combinations separately for each region, do
contrast(emm, "consec", by = "region")