data <- matrix(scan(),ncol=3,byrow=TRUE)
 1 1 106  
 1 1 108  
 1 2  93 
 1 2  101 
 1 2  98 
 1 3 56  
 2 1 107 
 2 1 110
 2 1 116
 2 2  63
 2 2  60
 2 3 40 
 2 3 41
 2 3 44 
 

data
  
UnbalConcrete <- data.frame(data)
colnames(UnbalConcrete) <- c("aggregate","compaction","strength")
UnbalConcrete$aggregate  <- as.factor(UnbalConcrete$aggregate)
UnbalConcrete$compaction <- as.factor(UnbalConcrete$compaction)
rm(data)

levels(UnbalConcrete$aggregate) <- c("Basalt","Silicious") 
levels(UnbalConcrete$compaction) <- c("Regular","Low","VeryLow") 

UnbalConcrete

 # the anova function applied to an lm object calculates Type I SS

anova(lm(strength ~ aggregate + compaction + aggregate:compaction, data = UnbalConcrete))



 # the car package contains a function 'Anova' that can calculate Type II or Type III SS
 # it calculates Type II SS by default

library(car)

 # this produces the results using Type II SS 
Anova(lm(strength ~ aggregate + compaction + aggregate:compaction, data = UnbalConcrete),type="II")

 # the Type III SS produced by using the type="III" argument do not match SAS results


library(phia)

# interactionMeans gives least squares means

interactionMeans(unbal.lm,factors="aggregate")
interactionMeans(unbal.lm,factors="compaction")
interactionMeans(unbal.lm)


 # 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)

UnbalConcrete$id <- rownames(UnbalConcrete)

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

summary(Unbalconc.afex1)

nice(Unbalconc.afex1)

lsmeans(Unbalconc.afex1, ~aggregate)

lsmeans(Unbalconc.afex1, ~compaction)