all 2-way interactions for 3-way interaction

require(lsmeans)
require(afex)
data(obk.long, package = "afex")
#lsm.options(disable.pbkrtest = TRUE)  pbkrtest needed for F-test

# let's fit this example data with 2 between. and two within-subjects factors via the ANOVA equivalent using lme4:
m1 <- mixed(value ~ treatment*gender*phase*hour + (1|id) + (1|id:phase) + (1|id:hour), obk.long)
##                         Effect     F ndf   ddf F.scaling p.value
## 1                    treatment  3.94   2 10.00      1.00     .05
## 2                       gender  3.66   1 10.00      1.00     .08
## 3                        phase 16.13   2 20.00      1.00 <.0001
## 4                         hour 16.69   4 40.00      1.00 <.0001
## 5             treatment:gender  2.86   2 10.00      1.00     .10
## 6              treatment:phase  4.85   4 20.00      1.00    .007
## 7                 gender:phase  0.28   2 20.00      1.00     .76
## 8               treatment:hour  0.09   8 40.00      1.00 >.99
## 9                  gender:hour  0.45   4 40.00      1.00     .77
## 10                  phase:hour  1.18   8 80.00      1.00     .32
## 11      treatment:gender:phase  0.64   4 20.00      1.00     .64
## 12       treatment:gender:hour  0.62   8 40.00      1.00     .76
## 13        treatment:phase:hour  0.35  16 80.00      1.00     .99
## 14           gender:phase:hour  0.93   8 80.00      1.00     .50
## 15 treatment:gender:phase:hour  0.74  16 80.00      1.00     .75

# I am interested in the phase:hour interaction per levels of treatment. It has 8 numerator df.

# set up the interaction per treatment:
(ls1 <- lsmeans(m1$full.model, ~phase:hour |treatment))
# alternatively:
# ls1 <- lsmeans(m1, ~phase:hour|treatment, data = obk.long)

# (following code by Russel Lenth)

#  create the contrasts for phase – 2 d.f. for each combination of the other factors:
(con1 <- contrast(ls1, "trt.vs.ctrl1", by = c("hour","treatment")))

# Now contrast these contrasts for each contrast (label in above results) and treatment combination, and summarize them by treatment
con2 <- contrast(con1, "trt.vs.ctrl1", by = c("contrast","treatment"))
summary(con2, by = "treatment")

# Finally, use the new capability of test to obtain the joint tests
test(con2, joint = TRUE, by = "treatment")
##              F df1 df2 p.value
## control 0.4090   8  80  0.9123
## A       0.9879   8  80  0.4518
## B       0.3838   8  80  0.9263

# Finally, use the new capability of test to obtain the joint tests
test(con2, joint = TRUE, rows = 1:8)

test(con2, joint = TRUE, rows = 9:16)

test(con2, joint = TRUE, rows = 17:24)


## Alternatively in next versions of lsmeans:
# test(con2, joint = TRUE, by = "treat")