data(CO2, package = "datasets")

# CO2 is a real split-plot-style trial: 12 plants (2 Types x 2 Treatments x 3
# replicates), each measured at 7 CO2 concentrations. Keeping two
# concentrations gives a perfectly balanced pooled split plot (1,1):
#   Location = Type       (2: Quebec, Mississippi)
#   Block    = replicate  (3, nested within Location)
#   FactorA  = Treatment  (2: nonchilled, chilled)   - main plot
#   FactorB  = conc       (2: 175, 675 mL/L)         - sub plot
d <- subset(CO2, conc %in% c(175, 675))
d <- data.frame(
  Location   = factor(d$Type, labels = c("Quebec", "Mississippi")),
  FactorA    = factor(d$Treatment, labels = c("A_nonchilled", "A_chilled")),
  FactorB    = factor(ifelse(d$conc == 175, "B_175", "B_675")),
  Block      = factor(substr(as.character(d$Plant), 3, 3)),
  Uptake     = d$uptake,
  Efficiency = d$uptake / d$conc * 100          # second trait, needed for PCA
)
nrow(d)                                                         # 24
all(table(d$Location, d$Block, d$FactorA, d$FactorB) == 1)      # TRUE

chars <- c("Uptake", "Efficiency")
d_tr  <- d                                       # Aitken-transformed copy

for (ch in chars) {
  # One aov per Location. Block:FactorA is the whole-plot error, so the
  # residuals are the sub-plot error of that Location.
  fits <- lapply(split(d, d$Location), function(sub) {
    sub <- droplevels(sub)
    aov(as.formula(paste(ch, "~ Block + FactorA + Block:FactorA + FactorB + FactorA:FactorB")),
        data = sub)
  })
  mse <- sapply(fits, function(m) sum(resid(m)^2) / df.residual(m))

  # Residual-based Bartlett test: all residuals in one vector, labelled by Location
  res  <- lapply(fits, residuals)
  test <- bartlett.test(unlist(res), rep(names(fits), sapply(res, length)))
  print(c(chisq = unname(test$statistic), p = test$p.value, mse))

  # Aitken fires on a HARDCODED 0.05, not the user's alpha
  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])   # divide by sqrt(MSE)
    }
  }
}
#> Uptake:     chisq = 0.56  p = 0.4538  MSE Quebec = 6.97, Mississippi = 4.38
#> Efficiency: chisq = 0.83  p = 0.3616  MSE Quebec = 1.43, Mississippi = 0.81
#> -> both homogeneous, no Aitken transformation

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 + B + A:B + L:A + L:B + L:A:B + L:Block:A)
at  <- as.data.frame(car::Anova(fit, type = "II"))
rownames(at)
#> "L" "A" "B" "L:Block" "A:B" "L:A" "L:B" "L:A:B" "L:Block:A" "Residuals"
MS  <- at[, "Sum Sq"] / at[, "Df"]

MS_Ea <- MS[9];  df_Ea <- at[9, "Df"]            # L:Block:A  -> Pooled Error (a)
MS_Eb <- MS[10]; df_Eb <- at[10, "Df"]           # Residuals  -> Pooled Error (b)
c(df_Ea = df_Ea, MS_Ea = MS_Ea, df_Eb = df_Eb, MS_Eb = MS_Eb)
#> df_Ea = 4  MS_Ea = 14.78   df_Eb = 8  MS_Eb = 5.67
# matching 2*(3-1)*(2-1) = 4 and 2*2*(3-1)*(2-1) = 8

# Each row is tested against the stratum the module assigns to it
Fp <- function(k, MS_err, df_err) {
  Fv <- MS[k] / MS_err
  c(F = Fv, p = pf(Fv, at[k, "Df"], df_err, lower.tail = FALSE))
}
Fp(2, MS_Ea, df_Ea)      # Main plot A      -> Error (a)   F = 19.33   p = 0.0117
Fp(1, MS_Eb, df_Eb)      # Location main    -> Error (b)   F = 161.29  p < 0.0001
Fp(3, MS_Eb, df_Eb)      # Sub plot B       -> Error (b)   F = 98.82   p < 0.0001
Fp(7, MS_Eb, df_Eb)      # Location x Sub   -> Error (b)   F = 8.00    p = 0.0222
Fp(9, MS_Eb, df_Eb)      # Error (a) itself -> Error (b)   F = 2.60    p = 0.1161

library(agricolae); library(gtools)
alpha <- 0.05                                    # effective shipped default

# Means table first, in natural (mixedsort) order - the letters follow it
meanL <- aggregate(y ~ L, FUN = mean); rownames(meanL) <- meanL$L
meanL <- meanL[mixedsort(rownames(meanL)), , drop = FALSE]

# Location -> Error (b)
outL <- LSD.test(y, L, DFerror = df_Eb, MSerror = MS_Eb, alpha = alpha)
outL$groups[rownames(meanL), ]                   # Mississippi 20.94 b, Quebec 33.29 a
outL$statistics[1, 6]                            # CD = 2.2424  (column 6 for LSD)

# Main plot -> Error (a)
outA <- LSD.test(y, A, DFerror = df_Ea, MSerror = MS_Ea, alpha = alpha)
outA$statistics[1, 6]                            # CD = 4.3571
r <- outA$means[1, 3]                            # replications per mean
sqrt(outA$statistics[1, 1] / r)                  # SE(m) = 1.1097
sqrt(2 * outA$statistics[1, 1] / r)              # SE(d) = 1.5693
outA$statistics[1, 4]                            # CV(%) = 14.1757

# Tukey: the CD equivalent is in column 5, not 6
HSD.test(y, A, DFerror = df_Ea, MSerror = MS_Ea, alpha = alpha)$statistics[1, 5]

# DMRT: the module's duncan_test() (duncan_test.R) - same output as
# agricolae::duncan.test() whenever qtukey() is finite
dun <- duncan.test(y, B, DFerror = df_Eb, MSerror = MS_Eb, alpha = alpha)
dun$duncan[, 1]                                  # tabulated value  = 3.2612
dun$duncan[, 2]                                  # critical range   = 2.2424

library(dplyr)
# Computed on the uploaded (untransformed) values, grouped by one factor at a time
d %>% group_by(FactorA) %>%
  summarise(Mean = mean(Uptake), SD = sd(Uptake), SE = SD / sqrt(n()),
            Min = min(Uptake), Max = max(Uptake), CV = SD / Mean * 100,
            Skewness = moments::skewness(Uptake), Kurtosis = moments::kurtosis(Uptake))
#> A_nonchilled  Mean 30.6  SD 8.09  SE 2.34  Min 19.2  Max 43.9  CV 26.5  Skew 0.136  Kurt 2.05
#> A_chilled     Mean 23.7  SD 9.47  SE 2.73  Min 11.4  Max 39.6  CV 40.0  Skew 0.506  Kurt 1.97

library(factoextra); library(scales)

# Means per Location x Main x Sub combination, row-ordered by mixedorder().
# PCA uses the uploaded (user-transformed) values, not the Aitken-scaled copy.
d$AxBxL <- interaction(d$Location, d$FactorA, d$FactorB, sep = "x")
means <- aggregate(cbind(Uptake, Efficiency) ~ AxBxL, data = d, FUN = mean)
rownames(means) <- means$AxBxL
means <- means[mixedorder(rownames(means)), -1]           # 8 combinations

res.pca <- prcomp(means, center = TRUE, scale = TRUE)     # correlation PCA via SVD

# The app reverses both sign conventions - reproduce this to match its output
res.pca$x        <- -res.pca$x
res.pca$rotation <- -res.pca$rotation

get_eigenvalue(res.pca)          # PC1 = 1.114 (55.68%), PC2 = 0.886 (44.32%)
fviz_eig(res.pca, addlabels = TRUE)
fviz_pca_biplot(res.pca, axes = c(1, 2), repel = TRUE)

index1 <- res.pca$x[, 1]
cbind(`Index Score` = index1, `Scaled Index` = rescale(index1, to = c(0, 1)))

# The Q-Q panel uses this model, which adds L:Block:B to the ANOVA model
m <- lm(y ~ L + L:Block + A + L:A + L:Block:A + B + L:B + L:Block:B + A:B + L:A:B)
z <- residuals(m); n <- length(z)
if (n >= 3 && n <= 5000 && var(z) > 0) shapiro.test(z)
if (n >= 7 && var(z) > 0)              nortest::ad.test(z)
if (n >= 2 && var(z) > 0)              ks.test(z, "pnorm", mean(z), sd(z))