set.seed(123)
design <- expand.grid(A = factor(c("a1", "a2")),
                       B = factor(c("b1", "b2")),
                       C = factor(c("c1", "c2", "c3")),
                       rep = 1:4)                      # balanced 2 x 2 x 3 factorial, 4 reps/cell
design$y <- with(design, 10 + 2 * (A == "a2") + 1.5 * (B == "b2") +
                    1 * (C == "c2") + 2 * (C == "c3") +
                    rnorm(nrow(design), sd = 2))

A <- design$A; B <- design$B; C <- design$C; y <- design$y  # Factors A, B, C and response
AB  <- factor(paste(A, B, sep = ":"))                  # A x B interaction cell
AC  <- factor(paste(A, C, sep = ":"))                  # A x C interaction cell
BC  <- factor(paste(B, C, sep = ":"))                  # B x C interaction cell
ABC <- factor(paste(A, B, C, sep = ":"))                # A x B x C interaction cell

# Full factorial linear model and Type II ANOVA (mirrors the module)
model  <- lm(y ~ A + B + C + A:B + A:C + B:C + A:B:C)
result <- as.data.frame(car::Anova(model, type = "II"))
result

# Error (residual) mean square, reused by the post-hoc tests
MSE <- result[8, 1] / result[8, 2]                     # residual is the 8th row (7 terms + residual)

# Partial Cohen's f effect size for A, B, C and all interaction terms
effectsize::cohens_f(model)

alpha <- 0.05                                    # default significance level
dferr <- result[8, 2]                            # residual degrees of freedom

# Least Significant Difference (default): agricolae uses the pooled MSE and
# residual df from the factorial model; grouping letters ordered with mixedsort.
outA   <- agricolae::LSD.test(y, A,   DFerror = dferr, MSerror = MSE, alpha = alpha)
outB   <- agricolae::LSD.test(y, B,   DFerror = dferr, MSerror = MSE, alpha = alpha)
outC   <- agricolae::LSD.test(y, C,   DFerror = dferr, MSerror = MSE, alpha = alpha)
outAB  <- agricolae::LSD.test(y, AB,  DFerror = dferr, MSerror = MSE, alpha = alpha)
outAC  <- agricolae::LSD.test(y, AC,  DFerror = dferr, MSerror = MSE, alpha = alpha)
outBC  <- agricolae::LSD.test(y, BC,  DFerror = dferr, MSerror = MSE, alpha = alpha)
outABC <- agricolae::LSD.test(y, ABC, DFerror = dferr, MSerror = MSE, alpha = alpha)

outA$groups[gtools::mixedsort(rownames(outA$groups)), ]
outABC
# Optional p-value adjustment for LSD (default "none"; selectable: bonferroni, holm, BH).
# agricolae applies the correction internally and recomputes the grouping letters.
agricolae::LSD.test(y, A, DFerror = dferr, MSerror = MSE, alpha = alpha,
                    p.adj = "bonferroni")$groups

# Alternatives selectable in the app (same DFerror / MSerror / alpha arguments):
#   Tukey HSD -> agricolae::HSD.test(y, A, DFerror = dferr, MSerror = MSE, alpha = alpha)
#   DMRT      -> agricolae::duncan.test(y, A, DFerror = dferr, MSerror = MSE, alpha = alpha)

# The app's MANOVA mode requires >= 2 response columns (one per measured trait).
# Illustrated here with two response columns built from the same design.
Y <- cbind(y1 = y, y2 = y + rnorm(length(y), sd = 1))

res.man <- manova(as.matrix(Y) ~ A + B + C + A:B + A:C + B:C + A:B:C)
broom::tidy(res.man)                             # Pillai's trace per term
effectsize::eta_squared(res.man)                 # partial eta-squared per term