### R code for sections 7.3 & 7.4
## R code for Example 7.3-2 (Table 7.3-5 on page 411)
x1 = rep(c(-1,1),8)
x2 = rep(c(-1,-1,1,1),4)
x3 = rep(c(-1,-1,-1,-1,1,1,1,1),2)
x4 = c(rep(-1,8),rep(1,8))
Y = c(42,31,45,29,39,28,46,32,40,30,50,25,40,25,50,23)
effect = rep(0, 16) # estimated effects
names(effect) = c("mu","1","2","3","4","12","13","14","23","24","34","123","124","134","234","1234")
effect["mu"] = sum(Y)/16
effect["1"] = sum(Y*x1)/16
effect["2"] = sum(Y*x2)/16
effect["3"] = sum(Y*x3)/16
effect["4"] = sum(Y*x4)/16
effect["12"] = sum(Y*x1*x2)/16
effect["13"] = sum(Y*x1*x3)/16
effect["14"] = sum(Y*x1*x4)/16
effect["23"] = sum(Y*x2*x3)/16
effect["24"] = sum(Y*x2*x4)/16
effect["34"] = sum(Y*x3*x4)/16
effect["123"] = sum(Y*x1*x2*x3)/16
effect["124"] = sum(Y*x1*x2*x4)/16
effect["134"] = sum(Y*x1*x3*x4)/16
effect["234"] = sum(Y*x2*x3*x4)/16
effect["1234"] = sum(Y*x1*x2*x3*x4)/16
round(effect,2)
# Normal probability plot (QQ-plot)
# see Table 7.3-6 on page 413
x = sort(effect[2:16]) # ordered estimated effects
P = (seq(1:15)-0.5)/15
z = qnorm(P)
round(z, 2)
# similar to Figure 7.3-1 on page 414
# note that it is not exactly same as Figure 7.3-1
# for example, (134),(4),(3) are three separated points with same x value
# while in Figure 7.3-1, they are ploted as the same position
plot(x, z, xlab="Estimated effects", ylab="Normal scores")
text(x+0.25, z, labels=names(x)[order(x)],cex=0.5) # updated on 04/25/2010
# QQ-plot
qqnorm(effect[2:16])