rm(list = ls())  # clear the memory
seat = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/seatbelt_table_6_2.dat", header = F, quote= "")
strength = seat[,-c(9:11)]
flash = seat[,-c(6:8)]

################## strength data, factors A,B,C only
y = c(as.matrix(strength[, 6:8]))
#y = c(as.matrix(flash[, 6:8]))
A = rep(strength[,2],3)
B = rep(strength[,3],3)
C = rep(strength[,4],3)

mod3 = function(x) {x - 3*floor(x/3)}
AB = factor(mod3(A+B))
AB2 = factor(mod3(A+2*B))
AC = factor(mod3(A+C))
AC2 = factor(mod3(A+2*C))
BC = factor(mod3(B+C))
BC2 = factor(mod3(B+2*C))
ABC = factor(mod3(A+B+C))
ABC2 = factor(mod3(A+B+2*C))
AB2C = factor(mod3(A+2*B+C))
AB2C2 = factor(mod3(A+2*B+2*C))

A = factor(A)
B = factor(B)
C = factor(C)

fit0 = aov(y ~ (A + B + C)^3)
summary(fit0)

fit = aov(y ~ A + B + C + AB + AB2 + AC + AC2 + BC + BC2 + ABC + ABC2 + AB2C + AB2C2)
anova(fit)

windows()
par(mfrow = c(1,2))
resid = fit$resid
qqnorm(resid)
qqline(resid)
text(0, .9*max(resid), paste("p = ", round(shapiro.test(resid)[[2]],3)))
plot(fit$fitted,fit$resid, xlab = "fits", ylab = "resids")

windows()
par(mfrow = c(3,3))

for(i in 0:2) {
resid = fit$resid[A==i]
qqnorm(resid)
qqline(resid)
text(0, .9*max(resid), paste("p = ", round(shapiro.test(resid)[[2]],3)))
title(sub = paste("A = ", i))
}


for(i in 0:2) {
resid = fit$resid[B==i]
qqnorm(resid)
qqline(resid)
text(0, .9*max(resid), paste("p = ", round(shapiro.test(resid)[[2]],3)))
title(sub = paste("B = ", i))
}



for(i in 0:2) {
resid = fit$resid[C==i]
qqnorm(resid)
qqline(resid)
text(0, .9*max(resid), paste("p = ", round(shapiro.test(resid)[[2]],3)))
title(sub = paste("C = ", i))
}

windows()
par(mfrow=c(1,3))
plot(A, y, xlab = "A")
plot(B, y, xlab = "B")
plot(C, y, xlab = "C")


windows()
par(mfrow=c(1,3))
stripchart(y ~ A, vertical = T, xlab = "A")
stripchart(y ~ B, vertical = T, xlab = "B")
stripchart(y ~ C, vertical = T, xlab = "C")




######################################################################
######################################################################

rm(list = ls())  # clear the memory
seat = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/seatbelt_table_6_2.dat", header = F, quote= "")
strength = seat[,-c(9:11)]
flash = seat[,-c(6:8)]


################## strength data, all factors, 3^(4-1) design
y = c(as.matrix(strength[, 6:8]))
A = rep(strength[,2],3)
B = rep(strength[,3],3)
C = rep(strength[,4],3)
D = rep(strength[,5],3) 

mod3 = function(x) {x - 3*floor(x/3)}
AB = factor(mod3(A+B))
AB2 = factor(mod3(A+2*B))
AC = factor(mod3(A+C))
AC2 = factor(mod3(A+2*C))
BC = factor(mod3(B+C))
BC2 = factor(mod3(B+2*C))
AD = factor(mod3(A+D))
BD = factor(mod3(B+D))
CD = factor(mod3(C+D))
#CD2 = factor(mod3(C+2*D))

A = factor(A)
B = factor(B)
C = factor(C)
D = factor(D)

smat = as.matrix(strength[, 6:8])
sbar = apply(smat, 1, mean)
sdisp = apply(smat, 1, var)
slogdisp = log(sdisp)

a = A[1:27]
b = B[1:27]
c = C[1:27]
d = D[1:27]

############## Location analysis

data = data.frame(sbar, a, b, c, d)
plot.design(data)
title(sub = "main effects; mean(strength)")

windows()
par(mfrow=c(3,2))
interaction.plot(b,a,sbar)
interaction.plot(c,a,sbar)
interaction.plot(d,a,sbar)
interaction.plot(c,b,sbar)
interaction.plot(d,b,sbar)
interaction.plot(d,c,sbar)
title(main = "interaction effects; mean(strength)", outer = T, line = -1)

fit0 = aov(y ~ (A + B + C + D)^3)
summary(fit0)

fit = lm(y ~ A + B + AB + AB2 + C + AC + AC2 + BC + BC2 + D + AD + BD + CD)
anova(fit)

windows()
par(mfrow=c(2,1))
# qq plot
eff = anova(fit)$"Sum Sq"[1:13]
S2 = sum(fit$resid^2)/fit$df.resid
eff = eff/S2
names(eff) = c("A","B","AB","AB2","C","AC","AC2","BC","BC2","D","AD","BD","CD")
plot(qchisq(rank(eff)/14, 2), eff, pch="", xlab = "SS vs. Exponential quantiles")
text(qchisq(rank(eff)/14, 2), eff, labels = names(eff))

# qq plot
eff = -2*log(anova(fit)$"Pr(>F)"[1:13])
names(eff) = c("A","B","AB","AB2","C","AC","AC2","BC","BC2","D","AD","BD","CD")
plot(qchisq(rank(eff)/14, 2), eff, pch="", xlab = "-log(p-value) vs. Exponential quantiles")
text(qchisq(rank(eff)/14, 2), eff, labels = names(eff))

title(main = "qqplots; mean(strength)", outer = T, line = -1)

############## Dispersion analysis

windows()
data = data.frame(slogdisp, a, b, c, d)
par(mfrow=c(1,1))
plot.design(data)
title(sub = "main effects; disp(strength)")


windows()
par(mfrow=c(3,2))
interaction.plot(b,a,slogdisp)
interaction.plot(c,a,slogdisp)
interaction.plot(d,a,slogdisp)
interaction.plot(c,b,slogdisp)
interaction.plot(d,b,slogdisp)
interaction.plot(d,c,slogdisp)
title(main = "interaction effects; disp(strength)", outer = T, line = -1)


# Revert to a single replicate with 27 runs:

A = rep(strength[,2],1)
B = rep(strength[,3],1)
C = rep(strength[,4],1)
D = rep(strength[,5],1) 

AB = factor(mod3(A+B))
AB2 = factor(mod3(A+2*B))
AC = factor(mod3(A+C))
AC2 = factor(mod3(A+2*C))
BC = factor(mod3(B+C))
BC2 = factor(mod3(B+2*C))
AD = factor(mod3(A+D))
BD = factor(mod3(B+D))
CD = factor(mod3(C+D))


A = factor(A)
B = factor(B)
C = factor(C)
D = factor(D)

dispfit = aov(slogdisp ~ A + B + AB + AB2 + C + AC + AC2 + BC + BC2 + D + AD + BD + CD)
summary(dispfit)

eff = anova(dispfit)$"Sum Sq"[1:13]
names(eff) = c("A","B","AB","AB2","C","AC","AC2","BC","BC2","D","AD","BD","CD")
windows()
plot(qchisq(rank(eff)/14, 2), eff, pch="", xlab = "Exponential quantiles")
text(qchisq(rank(eff)/14, 2), eff, labels = names(eff))
title(main = "qqplots; disp(strength)", outer = T, line = -1)