data(warpbreaks)
env <- warpbreaks$wool                              # 2 environments (A, B)

# one CRD per environment; Bartlett across their error variances
per_env <- lapply(split(warpbreaks, env),
                  function(d) aov(breaks ~ tension, data = d))
bt <- bartlett.test(per_env)                        # homogeneity across environments
bt

# per-environment MSE (used by Aitken when Bartlett is significant)
mse <- sapply(per_env, function(m) sum(resid(m)^2) / df.residual(m))
mse
# if bt$p.value <= 0.05: divide each environment's values by sqrt(its MSE) before pooling

library(car)
fit  <- lm(breaks ~ wool + tension + wool:tension,  # environment + treatment + interaction
           data = warpbreaks)
atab <- Anova(fit, type = "II")                     # Type II sums of squares
atab

MS <- atab[, "Sum Sq"] / atab[, "Df"]               # mean square = SS / df
MS
# F-ratios use the error term implied by the chosen pooled model; p-values via:
# pf(Fvalue, df1, df2, lower.tail = FALSE)

library(agricolae)
dfe      <- df.residual(fit)
mse_pool <- sum(resid(fit)^2) / dfe
alpha    <- 0.05                                    # default (0.01 optional)

# LSD is the default; TUKEY and DMRT are the alternatives
LSD.test(warpbreaks$breaks, warpbreaks$tension,
         DFerror = dfe, MSerror = mse_pool, alpha = alpha, p.adj = "none")$groups
HSD.test(warpbreaks$breaks, warpbreaks$tension,
         DFerror = dfe, MSerror = mse_pool, alpha = alpha)$groups
duncan.test(warpbreaks$breaks, warpbreaks$tension,
            DFerror = dfe, MSerror = mse_pool, alpha = alpha)$groups

x <- warpbreaks$breaks

log_x  <- if (min(x) <= 0) log10(x - min(x) + 1) else log10(x)   # logarithmic
sqrt_x <- if (any(x == 0)) sqrt(x + 0.5)         else sqrt(x)    # square root
p      <- (x - min(x)) / (max(x) - min(x))
asin_x <- asin(sqrt(p))                                          # arcsine (proportions)

# Aitken (automatic, applied per environment only when Bartlett is significant):
#   x_within_environment / sqrt(MSE_of_that_environment)

library(scales)

# illustrative multi-trait treatment means (rows = treatments, cols = traits)
means <- aggregate(cbind(Sepal.Length, Sepal.Width, Petal.Length) ~ Species,
                   data = iris, FUN = mean)
X   <- scale(means[, -1])
pca <- prcomp(X, center = TRUE, scale. = TRUE)      # requires >= 2 traits

index <- rescale(pca$x[, 1], to = c(0, 1))          # PCA-based selection index in [0, 1]
data.frame(Treatment = means$Species, Index = round(index, 3))

cor(iris[, 1:4])                                     # trait correlation matrix