rm(list = ls())  # clear the memory
height = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/LeafSpring_table_5_2.dat",
	 header = F, quote= "")
B = as.factor(height[1:16,1])
C = as.factor(height[1:16,2])
D = as.factor(height[1:16,3])
E = as.factor(height[1:16,4])
Q = as.factor(height[1:16,5])

ymat = as.matrix(height[,6:8])
y = c(ymat)

ybar = apply(ymat, 1, mean)
s2 = apply(ymat, 1, var)
logs2 = log(s2)
data = data.frame(y, B, C, D, E, Q)

fit = aov(y ~ (B + C + D + E + Q)^5, data)
summary(fit)

xB =  2*as.numeric(data$B=="+")-1
xC =  2*as.numeric(data$C=="+")-1
xD =  2*as.numeric(data$D=="+")-1
xE =  2*as.numeric(data$E=="+")-1
xQ =  2*as.numeric(data$Q=="+")-1

# Response = y
fit2 = lm(y ~(xB + xC + xD + xE + xQ)^5)
anova(fit2) #Q, CQ, C, B most significant, E somewhat so, maybe BQ
fit2a = lm(y ~ xB + xE + xC + xQ + xB:xQ + xC:xQ)
summary(fit2a)

# Response is log(S^2)

xB = xB[1:16]
xC = xC[1:16]
xD = xD[1:16]
xE = xE[1:16]
xQ = xQ[1:16]
fit3 = lm(logs2 ~(xB + xC + xD + xE + xQ)^5)

# load the BsMD package
library(BsMD)
windows()
par(mfrow=c(1,2))
DanielPlot(fit2, half = T, sub = "reponse = Y")
DanielPlot(fit3, half = T, sub = "reponse = log S^2")


fit3a = lm(logs2 ~ xB + xD:xQ + xB:xC:xQ + xQ) # hierarchy
summary(fit3a)


library(car)
leveneTest(y, group = as.factor(rep(1:16,3)))


# Pooled estimate of variance:
S=sqrt(mean(s2))
S # compare with 'fit'

# Nominal the best analysis
tab = model.tables(aov(fit3a), type = "effects")
tab
# Should have B = -1, DQ = -1, BCQ = +1, Q = -1; hence C = +1, D = +1: (B,C,D,Q) = (-,+,+,-)
# fit2a$coef = (Int, B, E, C, Q, BQ, CQ)
sum(fit2a$coef*c(1,-1,0,1,-1,1,-1)); fit2a$coef[3]
# At these levels y =  7.868333 + 0.051875 xE

XE = (8 - 7.868333)/0.051875  # Need xE = 2.538159; 
2.5 + .5*XE		      #  so E = 3.76908; out of range