Code:

################################################################################
## TABLE 16.5 (p 715) #
## Randomization Samples and Test Statistics--Quality Control Example #
## FIGURE 16.8 #
## Randomization Distribution of F* and Corresponding #
## F Distribution--Quality Control Example #
## #
## Since there is no algorithm to compute this example we had to devise one. #
## It should come as rather straight-forward. The Xi's are as in the above #
## examples. The 'y' will hold the 1,680 cases of 9-sequences consisting of #
## the response variables. The 'ti' implies the treatment group. In this case #
## t1 is the first group (3-sequence) and t12 is the composite of t1 and t2. #
## The 'remainder' function is a wrapper for grabing a subset of 'set' based #
## on those values not in 'x'. The 'seq6' is the 6-sequence remainder after t1 #
## is defined. The whole process took less than 10 seconds on a 2.4 GHz #
## processor. As for the output, the columns are arbitrarily labeled 1-9. #
## Clearly they represent the three treatment groups based on groups of three. #
## The function 'f' uses the matrix algebra discussed in Ch. 5. It is possible #
## to get away with merely fitting an 'lm' object, and then extract the #
## f-statistic in a single call. However, this requires a lot of additional #
## work for each of the 1680 rows. It took somewhere between 30-60 seconds to #
## produce the same result. #
################################################################################
remainder <- function(x, set) set[!set %in% x]
f <- function(Y, X) {
Y <- matrix(Y) ## Turn row-vector into column
p <- ncol(X); n <- nrow(X)
J <- matrix(1, n, n) ## (5.18)
H <- X %*% solve(t(X) %*% X) %*% t(X) ## (5.73a)
SSE <- t(Y) %*% (diag(n) - H) %*% Y ## (5.89b)
SSR <- t(Y) %*% (H - (1/n)*J) %*% Y ## (5.89c)
fstar <- (SSR / (p - 1)) / (SSE / (n - p)) ## (6.39b)
}
base <- c(1.1, 0.5, -2.1, 4.2, 3.7, 0.8, 3.2, 2.8, 6.3)
t2 <- t12 <- t123 <- list()
y <- NULL
X <- cbind(
X1 = c(1, 1, 1, 0, 0, 0, 0, 0, 0),
X2 = c(0, 0, 0, 1, 1, 1, 0, 0, 0),
X3 = c(0, 0, 0, 0, 0, 0, 1, 1, 1)
);
t1 <- t(combn(base, 3))
seq6 <- t(combn(base, 3, remainder, set = base))
for (i in 1:84) t2[[i]] <- t(combn(seq6[i, ], 3))
for (i in 1:84) t12[[i]] <- cbind(t1[i, 1], t1[i, 2], t1[i, 3], t2[[i]])
for (i in 1:84)
t123[[i]] <- cbind(t12[[i]], t(apply(t12[[i]], 1, remainder, set = base)))
for (i in 1:84) y <- rbind(y, t123[[i]])
fstar <- apply(y, 1, function(Y) f(Y, X))
cbind(y, data.frame(f = fstar))
hist(fstar, freq = FALSE, ylim = c(0, 1), col = "gray90", main = "")
curve(df(x, 2, 6), add = TRUE, lwd = 2)
rm(base, fstar, i, remainder, seq6, t1, t2, t12, t123, f, X, y)