data <- scan(what = list(E="",F="",y=0), multi.line=F)
                lowlow      lowlow      26.1
                lowlow      lowlow      27.5
                lowlow      lowhigh     23.5
                lowlow      lowhigh     21.1
                lowlow      highlow     22.8
                lowlow      highlow     23.8
                lowlow      highhigh    30.6
                lowlow      highhigh    32.5
                lowhigh     lowlow      22.0
                lowhigh     lowlow      20.2
                lowhigh     lowhigh     28.1
                lowhigh     lowhigh     29.9
                lowhigh     highlow     30.0
                lowhigh     highlow     29.3
                lowhigh     highhigh    38.3
                lowhigh     highhigh    38.5
                highlow     lowlow      11.4
                highlow     lowlow      11.0
                highlow     lowhigh     20.4
                highlow     lowhigh     22.0
                highlow     highlow     22.3
                highlow     highlow     20.2
                highlow     highhigh    28.7
                highlow     highhigh    28.8
                highhigh    lowlow      18.9
                highhigh    lowlow      16.4
                highhigh    lowhigh     26.6
                highhigh    lowhigh     26.5
                highhigh    highlow     29.6
                highhigh    highlow     29.8
                highhigh    highhigh    34.5
                highhigh    highhigh    34.9


data

oex9_2 <- data.frame(data)

oex9_2$E <- as.factor(oex9_2$E)
oex9_2$F <- as.factor(oex9_2$F)

oex9_2.lm1  <- lm(y ~ E +F +E:F, data=oex9_2)

anova(oex9_2.lm1)


 # examine the interaction

interaction.plot(oex9_2$E,oex9_2$F,
                 oex9_2$y,type="b")

 # check model assumptions

x11()

par(mfrow=c(2,2))  
plot(oex9_2.lm1)
par(mfrow=c(1,1))  


 # ----------------------------------------------------------------------------
 # excluding the E = "lowlow" group

oex9_2b   <- subset(oex9_2, E != "lowlow")


 # augmenting the visual approach from the profile plot

oex9_2b.lm1  <- lm(y ~ E +F +E:F, data=oex9_2b)

anova(oex9_2b.lm1)



 # ----------------------------------------------------------------------------
 # use the phia package to test simple main effects

library(phia)

testInteractions(oex9_2.lm1, fixed="E", across="F", adjustment="none")

testInteractions(oex9_2.lm1, fixed="F", across="E", adjustment="none")

 # here are the combination means

interactionMeans(oex9_2.lm1)