data(PlantGrowth)
d <- PlantGrowth                                 # weight ~ group, 3 treatments x 10 reps, balanced

treatment <- factor(d$group)                     # single treatment factor
y         <- d$weight                            # response variable

# One-way linear model and ANOVA table (mirrors the module; one model per response)
model  <- lm(y ~ treatment)
result <- anova(model)
result                                           # F = 4.85, p = 0.0159

# Error (residual) mean square and df, reused by the post-hoc tests
MSE   <- result[2, 3]                            # 0.3886
dferr <- result[2, 1]                            # 27

# Cohen's f effect size, derived from eta-squared
eta2 <- lsr::etaSquared(model)[1, 1]             # 0.2641
sqrt(eta2 / (1 - eta2))                          # 0.5991

alpha <- 0.05                                    # default significance level

# Least Significant Difference (default): agricolae uses the pooled MSE and
# residual df from the one-way model; grouping letters ordered with mixedsort.
out <- agricolae::LSD.test(y, treatment, DFerror = dferr, MSerror = MSE,
                           alpha = alpha, p.adj = "none", group = TRUE)

out$groups[gtools::mixedsort(rownames(out$groups)), ]
out$statistics                                   # MSerror, Df, Mean, CV, t.value, LSD (= CD)

# Grouping letters are printed only when the treatment effect is significant:
if (result[1, 5] > alpha) out$groups[, 2] <- ""

# 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, treatment, 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, treatment, DFerror = dferr, MSerror = MSE, alpha = alpha)
#   DMRT      -> agricolae::duncan.test(y, treatment, DFerror = dferr, MSerror = MSE, alpha = alpha)

# Standard errors from the pooled MSE (r = replications per treatment)
r   <- out$means[1, 3]                           # 10
SEM <- sqrt(MSE / r)                             # 0.1971
SED <- sqrt(2 * MSE / r)                         # 0.2788

# The app's MANOVA mode requires >= 2 response columns (one per measured trait)
# and more treatments than response variables. Illustrated with `iris`, a
# completely randomized layout of 3 treatments x 50 units and 4 measured traits.
data(iris)
trt <- factor(iris$Species)
Y   <- as.matrix(iris[, 1:4])

res.man <- manova(Y ~ trt)
broom::tidy(res.man)                             # Pillai's trace for the treatment term
effectsize::eta_squared(res.man)                 # partial eta-squared