Custom emmeans contrasts

38 Views Asked by At

I have a longitudinal study in which there are two treatments on day -3 but then individuals in each of these two treatments are further split into four treatments on day 0 and onward into day 2. In other words, the day -3 treatments no longer exist on days 0–2 and the day 0–2 treatments did not exist on day -3.

I would like to specify custom pairwise contrasts such that it removes all irrelevant contrasts that don't truly exist.

However, I can't quite figure out how to do this the long-hand way specifying each custom contrast based on the mean matrix nor using 'by' variables.

lm.hf <- lm(y ~ day*xtrt + sex, data=df)
emm <- emmeans(lm.hf, ~day*xtrt)
emm

 day xtrt       emmean   SE  df lower.CL upper.CL
 -3  MS only     18.54 1.06 132    16.44    20.64
 0   MS only    nonEst   NA  NA       NA       NA
 2   MS only    nonEst   NA  NA       NA       NA
 -3  HF only     15.92 1.04 132    13.87    17.97
 0   HF only    nonEst   NA  NA       NA       NA
 2   HF only    nonEst   NA  NA       NA       NA
 -3  MS:Control nonEst   NA  NA       NA       NA
 0   MS:Control  13.63 1.20 132    11.26    16.00
 2   MS:Control  10.65 1.20 132     8.29    13.02
 -3  HF:Control nonEst   NA  NA       NA       NA
 0   HF:Control  14.81 1.34 132    12.16    17.46
 2   HF:Control  13.45 1.47 132    10.55    16.35
 -3  MS:WT      nonEst   NA  NA       NA       NA
 0   MS:WT        5.66 1.16 132     3.37     7.95
 2   MS:WT        9.66 1.34 132     7.00    12.32
 -3  HF:WT      nonEst   NA  NA       NA       NA
 0   HF:WT        9.49 1.34 132     6.84    12.14
 2   HF:WT       13.90 1.34 132    11.25    16.54

Results are averaged over the levels of: sex 
Results are given on the sqrt (not the response) scale. 
Confidence level used: 0.95 

I've added a figure to provide additional clarity to the treatment design, if needed.

1

There are 1 best solutions below

0
Stéphane Laurent On

Instead of specifying custom contrasts, you can simply subset.

For example, to subset the non-estimable main effects, you can do:

library(emmeans)
fit <- lm(breaks ~ wool*tension, data = warpbreaks)
fit2 <- update(fit, subset = -(37:48))
emm <- emmeans(fit2, ~ wool*tension)
rg <- ref_grid(fit2)
emm[which(rg@grid$.wgt. != 0)]

To subset the pairwise comparisons table, you can do:

ps <- pairs(emmeans(fit2, ~ wool*tension))
sps <- summary(ps)
ps[which(!is.na(sps$estimate))]