d <- MASS::immer                                 # 5 varieties x 6 locations (blocks), 30 plots

treatment <- factor(d$Var)                       # treatment factor
block     <- factor(d$Loc)                       # block factor
y         <- d$Y1                                # response variable

# Blocked linear model and ANOVA table (mirrors the module; one model per response)
model  <- lm(y ~ treatment + block)
result <- anova(model)
result                                           # treatment F = 4.23, p = 0.0121
                                                 # block     F = 21.89, p = 1.75e-07

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

# Partial Cohen's f effect size; the module reports the treatment row
effectsize::cohens_f(model)[1, 2]                # 0.9199

alpha <- 0.05                                    # default significance level

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

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

# Letters and the CD value 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 = number of blocks)
r   <- out$means[1, 3]                           # 6
SEM <- sqrt(MSE / r)                             # 5.2104
SED <- sqrt(2 * MSE / r)                         # 7.3686

# The app's MANOVA mode requires >= 2 response columns (one per measured trait)
# and more treatments than response variables. Both factors are carried over
# from the univariate model, so the block effect is estimated as well.
Y <- as.matrix(d[, c("Y1", "Y2")])

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