data(CO2, package = "datasets")
d <- subset(CO2, conc %in% c(175, 350, 675))
d <- data.frame(
  Location   = factor(d$Type, labels = c("Quebec", "Mississippi")),
  FactorA    = factor(d$Treatment, labels = c("A_nonchilled", "A_chilled")),   # row strips
  FactorB    = factor(paste0("B_", d$conc)),                                     # column strips
  Block      = factor(substr(as.character(d$Plant), 3, 3)),
  Uptake     = d$uptake,
  Efficiency = d$uptake / d$conc * 100
)
nrow(d)                                                         # 36

d_tr <- d
for (ch in c("Uptake", "Efficiency")) {
  fits <- lapply(split(d, d$Location), function(sub) {
    sub <- droplevels(sub)
    aov(as.formula(paste(ch, "~ Block + FactorA + Block:FactorA + FactorB + Block:FactorB + FactorA:FactorB")),
        data = sub)
  })
  mse  <- sapply(fits, function(m) sum(resid(m)^2) / df.residual(m))
  res  <- lapply(fits, residuals)
  test <- bartlett.test(unlist(res), rep(names(fits), sapply(res, length)))
  if (test$p.value <= 0.05) {
    for (loc in names(mse)) {
      rows <- d$Location == loc
      d_tr[rows, ch] <- d[rows, ch] / sqrt(mse[loc])
    }
  }
}
#> Uptake:     chisq = 33.09  p < 0.0001  MSE Quebec = 4.38, Mississippi = 0.16  -> Aitken applied
#> Efficiency: chisq = 4.06   p = 0.0439  MSE Quebec = 1.50, Mississippi = 0.55  -> Aitken applied

library(car)
L <- d_tr$Location; A <- d_tr$FactorA; B <- d_tr$FactorB; Block <- d_tr$Block
y <- d_tr$Uptake

fit <- lm(y ~ L + L:Block + A + L:A + L:Block:A + B + L:B + L:Block:B + A:B + L:A:B)
at  <- as.data.frame(car::Anova(fit, type = "II"))
rownames(at)
#> "L" "A" "B" "L:Block" "L:A" "L:B" "A:B" "L:Block:A" "L:Block:B" "L:A:B" "Residuals"

MS <- at[, "Sum Sq"] / at[, "Df"]
Ea <- 8; Eb <- 9; Ec <- 11                       # Pooled Error (a), (b), (c)
Fp <- function(k, e) {
  Fv <- MS[k] / MS[e]
  c(F = Fv, p = pf(Fv, at[k, "Df"], at[e, "Df"], lower.tail = FALSE))
}
MS[Ec]       # 1 - after the Aitken transformation the pooled residual is exactly 1
Fp(2,  Ea)   # Row treatments (A)    F = 17.42     p = 0.0140
Fp(5,  Ea)   # Location x Row        F = 12.20     p = 0.0250
Fp(3,  Eb)   # Column treatments (B) F = 32.24     p = 0.0001
Fp(6,  Eb)   # Location x Column     F = 7.53      p = 0.0145
Fp(1,  Ec)   # Location              F = 12916.56  p < 0.0001
Fp(7,  Ec)   # Row x Column          F = 79.85     p < 0.0001
Fp(10, Ec)   # Location x Row x Col  F = 93.60     p < 0.0001

library(agricolae)
alpha <- 0.05

# Row treatments -> Pooled Error (a)
outA <- LSD.test(y, A, DFerror = at[Ea, "Df"], MSerror = MS[Ea], alpha = alpha)
# (values are on the Aitken-transformed scale)
outA$groups                          # A_nonchilled 42.69 a, A_chilled 28.55 b
outA$statistics[1, 6]                # CD = 9.40
r <- outA$means[1, 3]
sqrt(outA$statistics[1, 1] / r)      # SE(m) = 2.39
sqrt(2 * outA$statistics[1, 1] / r)  # SE(d) = 3.39
outA$statistics[1, 4]                # CV (%) = 28.52

# Column treatments -> Pooled Error (b)
outB <- LSD.test(y, B, DFerror = at[Eb, "Df"], MSerror = MS[Eb], alpha = alpha)
outB$groups                          # B_675 40.10 a, B_350 38.32 a, B_175 28.44 b
outB$statistics[1, 6]                # CD = 3.61
HSD.test(y, B, DFerror = at[Eb, "Df"], MSerror = MS[Eb], alpha = alpha)$statistics[1, 5]   # 4.47
duncan.test(y, B, DFerror = at[Eb, "Df"], MSerror = MS[Eb], alpha = alpha)$duncan          # critical ranges 3.61, 3.76

# Location -> Pooled Error (c)
LSD.test(y, L, DFerror = at[Ec, "Df"], MSerror = MS[Ec], alpha = alpha)$statistics[1, 6]   # CD = 0.77

library(gtools); library(factoextra); library(scales)
d$AxBxL <- interaction(d$FactorA, d$FactorB, d$Location, sep = "x")
means <- aggregate(cbind(Uptake, Efficiency) ~ AxBxL, data = d, FUN = mean)
rownames(means) <- means$AxBxL
means <- means[mixedorder(rownames(means)), -1]          # 12 combinations

res.pca <- prcomp(means, center = TRUE, scale = TRUE)
res.pca$x        <- -res.pca$x                           # sign convention used by the app
res.pca$rotation <- -res.pca$rotation
get_eigenvalue(res.pca)                                  # PC1 = 1.122 (56.08%), PC2 = 0.878 (43.92%)
rescale(res.pca$x[, 1], to = c(0, 1))                    # scaled Index 1