# STAT_516_Lec_06_inclass

y <- c(68,63,65,71,66,66,126,128,133,107,110,116,
       93,101,98,63,60,59,56,59,57,40,41,44)
agg <- as.factor(rep(c(rep("B",3),rep("S",3)),4))
comp <- as.factor(c(rep("st",6),rep("r",6),
                    rep("l",6),rep("vl",6)))


lm_out <- lm(y ~ agg + comp + agg:comp)
summary(lm_out)

a <- 2
b <- 4
n <- 3
sum(lm_out$residuals^2)/( a*b*(n-1))


plot(lm_out,which = 2)



boxplot(y ~ agg:comp)

aggregate(y, by = list(agg, comp), FUN = mean)

aggregate(y, by = list(agg), FUN = mean)

aggregate(y, by = list(comp), FUN = mean)

mean(y)

anova(lm_out)

interaction.plot(agg,comp,y)

interaction.plot(comp,agg,y)

# Tukey's for comparing all pairs of means in two-way factorial 

lm_tensile <- lm(y ~ agg + comp + agg:comp)
TukeyHSD(aov(lm_tensile))

# create one factor whose levels are all the combinations of the two factors
agg_comp <- as.factor(paste(agg,comp,sep="_"))
data.frame(agg,comp,agg_comp)

library(DescTools)
DunnettTest(y ~ agg_comp, control = "B_st")

# compare mean of Basalt and low versus Basalt and static with 
# Tukey adjustment for comparing all pairs of means WITHIN the 
# Basalt level of the aggregate factor
yBl <- mean(y[agg == "B" & comp == "l"])
yBst <- mean(y[agg == "B" & comp == "st"])
alpha <- 0.05
me <- qtukey(1-alpha/2,4,16) * 3.082 * sqrt(1/n)
lo <- yBl - yBst - me
up <- yBl - yBst + me
c(lo,up)


# Dunnett's for comparing all compaction methods to a 
# baseline compaction method for just the aggregate type Basalt

yBl <- mean(y[agg == "B" & comp == "l"])
yBst <- mean(y[agg == "B" & comp == "st"])
me <- 2.59 * 3.082 *sqrt( 2 / 3)
lo <- yBl - yBst - me
up <- yBl - yBst + me
c(lo,up)

# Glucose example

y <- c(42.5,43.3,42.9,138.4,144.4,142.7,180.9,180.5,183.0,
       39.8,40.3,41.2,132.4,132.4,130.3,176.8,173.6,174.9)
method <- as.factor(c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2))
glc <- as.factor(c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3))

boxplot(y ~ method + glc)
logy <- log(y)

lm_glucose <- lm(logy ~ method + glc + method:glc)
plot(lm_glucose,which=1)
plot(lm_glucose,which=2)

anova(lm_glucose)

interaction.plot(method,glc,logy)

interaction.plot(glc,method,logy)

# compare means for method 1 versus method 2
yM1 <- mean(logy[method == 1])
yM2 <- mean(logy[method == 2])

summary(lm_glucose)
me <- qtukey(.95,2,12) * 0.0135 /sqrt(9)
lo <- yM1 - yM2 - me
up <- yM1 - yM2 + me
c(lo,up)

summary(lm_glucose)