data(warpbreaks)
d <- warpbreaks                                  # breaks ~ wool (A) x tension (B), 2 x 3, balanced

A  <- factor(d$wool)                             # Factor A
B  <- factor(d$tension)                          # Factor B
AB <- factor(paste(d$wool, d$tension, sep = ":"))# A x B interaction cell (AxB)
y  <- d$breaks                                   # response variable

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

# Error (residual) mean square, reused by the post-hoc tests
MSE <- result[4, 1] / result[4, 2]

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

alpha <- 0.05                                    # default significance level
dferr <- result[4, 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)
outAB <- agricolae::LSD.test(y, AB, DFerror = dferr, MSerror = MSE, alpha = alpha)

outA$groups[gtools::mixedsort(rownames(outA$groups)), ]

# 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)

# Standard errors from the pooled MSE (r = replications per factor level)
r   <- outA$means[1, 3]
SEM <- sqrt(MSE / r)
SED <- sqrt(2 * MSE / r)

# 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(breaks     = warpbreaks$breaks,
           log_breaks = log1p(warpbreaks$breaks))

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