data <- scan(what = list(aggregate="",compaction="",strength=0), multi.line=F)
 Basalt Static 68  Basalt Static 63  Basalt Static 65
 Basalt Regular 126  Basalt Regular 128   Basalt Regular 133 
 Basalt Low 93   Basalt Low 101  Basalt Low  98 
 Basalt VeryLow 56   Basalt VeryLow 59   Basalt VeryLow 57 
 Silicious Static 71  Silicious Static 66  Silicious Static 66 
 Silicious Regular 107   Silicious Regular 110  Silicious Regular 116  
 Silicious Low 63  Silicious Low 60  Silicious Low 59 
 Silicious VeryLow 40   Silicious VeryLow 41   Silicious VeryLow 44 

 data

 concrete <- data.frame(data)
 rm(data)

concrete$aggregate <- as.factor(concrete$aggregate)
concrete$compaction <- as.factor(concrete$compaction)

conccrf <- lm(strength ~ aggregate +compaction +aggregate:compaction, data=concrete)

anova(conccrf)

 # examine the interaction

interaction.plot(concrete$aggregate,concrete$compaction,
                 concrete$strength,type="b")

 # check model assumptions

par(mfrow=c(2,2))  
plot(conccrf)
par(mfrow=c(1,1))  


 # We can use the afex package to obtain least squares means.
 # In the afex package, there are three ANOVA functions, here we
 # use the aov_ez function. 
 # The aov_ez function is a bit particular in that it requires a
 # separate id column in the data set. In code below, we create 
 # an identity column before calling the aov_ez function. 
 

library(afex)

concrete$id <- rownames(concrete)

conc.afex1 <- aov_ez(id="id",dv="strength",data=concrete,between=c("aggregate","compaction"))

summary(conc.afex1)

nice(conc.afex1)

lsmeans(conc.afex1, ~aggregate)