#problem 11
dat <- data.frame(
setup = 1:8,
flasher = c("A","B","A","B","A","B","A","B"),
inertia = c("Low","Low","High","High","Low","Low","High","High"),
task = c("Y","Y","Y","Y","Z","Z","Z","Z"),
y = c(11,12,10,11,16,14,15,19),
order = c(1,4,5,3,2,6,7,8)
)
# factors as -1 and +1
dat$F <- ifelse(dat$flasher == "B", 1, -1)
dat$I <- ifelse(dat$inertia == "High", 1, -1)
dat$T <- ifelse(dat$task == "Z", 1, -1)
# Interaction columns
dat$FI <- dat$F * dat$I
dat$FT <- dat$F * dat$T
dat$IT <- dat$I * dat$T
dat$FIT <- dat$F * dat$I * dat$T
#coded table
dat[, c("setup", "F", "I", "T", "FI", "FT", "IT", "FIT", "y", "order")]
## setup F I T FI FT IT FIT y order
## 1 1 -1 -1 -1 1 1 1 -1 11 1
## 2 2 1 -1 -1 -1 -1 1 1 12 4
## 3 3 -1 1 -1 -1 1 -1 1 10 5
## 4 4 1 1 -1 1 -1 -1 -1 11 3
## 5 5 -1 -1 1 1 -1 -1 1 16 2
## 6 6 1 -1 1 -1 1 -1 -1 14 6
## 7 7 -1 1 1 -1 -1 1 -1 15 7
## 8 8 1 1 1 1 1 1 1 19 8
# factorial effects
effects <- sapply(c("F","I","T","FI","FT","IT","FIT"), function(v) {
sum(dat[[v]] * dat$y) / 4
})
effects
## F I T FI FT IT FIT
## 1.0 0.5 5.0 1.5 0.0 1.5 1.5
abs_eff <- sort(abs(effects))
labels <- names(sort(abs(effects)))
n_eff <- length(abs_eff)
half_norm_quant <- qnorm((1:n_eff + n_eff - 0.5) / (2 * n_eff))
plot(half_norm_quant, abs_eff,
pch = 19,
xlab = "Half-normal quantiles",
ylab = "Absolute factorial effects",
main = "Half-normal Plot of Effects")
text(half_norm_quant, abs_eff, labels = labels, pos = 4, cex = 0.8)
abs.effects <- abs(effects)
s0 <- 1.5 * median(abs.effects)
PSE <- 1.5 * median(abs.effects[abs.effects < 2.5 * s0])
# Appendix H critical value for m = 7 and alpha = 0.05
c_ier <- 2.30
ME <- c_ier * PSE
s0
## [1] 2.25
PSE
## [1] 2.25
ME
## [1] 5.175
abs.effects > ME
## F I T FI FT IT FIT
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE
abs.effects > ME
## F I T FI FT IT FIT
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#Day 1: randomized orders 1 to 4
#Day 2: randomized orders 5 to 8
dat$day <- ifelse(dat$order <= 4, 1, 2)
dat$sigma2 <- ifelse(dat$day == 1, 1^2, 4^2)
#Variance of task effect
var_task <- sum((dat$T / 4)^2 * dat$sigma2)
var_task
## [1] 4.25
# problem 16
dat <- data.frame(
run = 1:8,
A = c(-1,-1,-1,-1, 1, 1, 1, 1),
B = c(-1,-1, 1, 1,-1,-1, 1, 1),
C = c(-1, 1,-1, 1,-1, 1,-1, 1),
y1 = c(54.6, 86.2, 41.4, 62.8, 59.6, 82.0, 43.4, 65.6),
y2 = c(73.0, 66.2, 51.2, 64.8, 52.8, 72.8, 49.0, 65.0),
y3 = c(139.2,79.2, 42.6, 74.6, 55.2, 76.6, 48.6, 64.2),
y4 = c(55.4, 86.0, 58.6, 74.6, 61.0, 73.4, 49.6, 60.8),
y5 = c(52.6, 82.6, 58.4, 64.6, 61.0, 75.0, 55.2, 77.4)
)
# interactions
dat$AB <- dat$A * dat$B
dat$AC <- dat$A * dat$C
dat$BC <- dat$B * dat$C
dat$ABC <- dat$A * dat$B * dat$C
# run means, variances, standard deviations, and log variances
dat$mean <- rowMeans(dat[, c("y1","y2","y3","y4","y5")])
dat$var <- apply(dat[, c("y1","y2","y3","y4","y5")], 1, var)
dat$sd <- sqrt(dat$var)
dat$logvar <- log(dat$var)
#summary table
dat[, c("run","A","B","C","AB","AC","BC","ABC","mean","var","sd","logvar")]
## run A B C AB AC BC ABC mean var sd logvar
## 1 1 -1 -1 -1 1 1 1 -1 74.96 1356.928 36.836504 7.212979
## 2 2 -1 -1 1 1 -1 -1 1 80.04 68.068 8.250333 4.220507
## 3 3 -1 1 -1 -1 1 -1 1 50.44 68.428 8.272122 4.225782
## 4 4 -1 1 1 -1 -1 1 -1 68.28 33.892 5.821684 3.523179
## 5 5 1 -1 -1 -1 -1 1 1 57.92 13.852 3.721828 2.628430
## 6 6 1 -1 1 -1 1 -1 -1 75.96 13.588 3.686190 2.609187
## 7 7 1 1 -1 1 -1 -1 -1 49.16 17.548 4.189033 2.864940
## 8 8 1 1 1 1 1 1 1 66.60 39.900 6.316645 3.686376
effects_mean <- sapply(c("A","B","C","AB","AC","BC","ABC"), function(v) {
sum(dat[[v]] * dat$mean) / 4
})
effects_mean
## A B C AB AC BC ABC
## -6.02 -13.60 14.60 4.54 3.14 3.04 -3.34
abs_mean <- sort(abs(effects_mean))
labels_mean <- names(sort(abs(effects_mean)))
n_eff <- length(abs_mean)
hnq <- qnorm((1:n_eff + n_eff - 0.5) / (2*n_eff))
plot(hnq, abs_mean,
pch = 19,
xlab = "Half-normal quantiles",
ylab = "Absolute location effects",
main = "Half-normal Plot: Location Effects")
text(hnq, abs_mean, labels = labels_mean, pos = 4, cex = 0.8)
effects_logvar <- sapply(c("A","B","C","AB","AC","BC","ABC"), function(v) {
sum(dat[[v]] * dat$logvar) / 4
})
effects_logvar
## A B C AB AC BC ABC
## -1.8483785 -0.5927063 -0.7232202 1.2495561 1.1243171 0.7826368 -0.3622973
abs_disp <- sort(abs(effects_logvar))
labels_disp <- names(sort(abs(effects_logvar)))
plot(hnq, abs_disp,
pch = 19,
xlab = "Half-normal quantiles",
ylab = "Absolute dispersion effects",
main = "Half-normal Plot: Dispersion Effects")
text(hnq, abs_disp, labels = labels_disp, pos = 4, cex = 0.8)
#Problem 20
# multiply factorial effects
multiply_effects <- function(words) {
chars <- unlist(strsplit(words, ""))
tab <- table(chars)
keep <- names(tab[tab %% 2 == 1])
paste(sort(keep), collapse = "")
}
#generating the full block defining contrast subgroup
block_subgroup <- function(gens) {
out <- c()
q <- length(gens)
for (r in 1:q) {
combs <- combn(gens, r, simplify = FALSE)
for (cmb in combs) {
out <- c(out, multiply_effects(cmb))
}
}
out <- unique(out)
out[order(nchar(out), out)]
}
#computing gi(b)
gi_values <- function(effects) {
table(factor(nchar(effects), levels = 1:6))
}
# First scheme
gens1 <- c("126", "136", "346", "456")
scheme1 <- block_subgroup(gens1)
scheme1
## [1] "14" "23" "25" "35" "126" "136" "156" "246" "346"
## [10] "456" "1234" "1245" "1345" "12356" "23456"
gi_values(scheme1)
##
## 1 2 3 4 5 6
## 0 4 6 3 2 0
# Second scheme from Table 4.1, k = 6, q = 4
gens2 <- c("12", "34", "56", "135")
scheme2 <- block_subgroup(gens2)
scheme2
## [1] "12" "34" "56" "135" "136" "145" "146" "235"
## [9] "236" "245" "246" "1234" "1256" "3456" "123456"
gi_values(scheme2)
##
## 1 2 3 4 5 6
## 0 3 8 3 0 1