set.seed(1)
data(Soils, package = "carData")

A     <- Soils$Contour                             # Factor A: slope position (3 levels)
B     <- Soils$Depth                               # Factor B: sampling depth (4 levels)
AB    <- factor(paste(Soils$Contour, Soils$Depth, sep = ":"))  # A x B cell (AxB)
Block <- Soils$Block 
# response variable
Ca <- Soils$Ca   
# Full RBD factorial linear model and Type II ANOVA (mirrors the module)
model  <- lm(Ca ~ Block + A + B + A:B)
result <- as.data.frame(car::Anova(model, type = "II"))
result

# Error (residual) mean square; rows are Block, A, B, A:B, Residuals -> row 5
MSE <- result[5, 1] / result[5, 2]

# Partial Cohen's f for A, B and A:B (row 1 of the model is Block)
effectsize::cohens_f(model)

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

# Least Significant Difference (default): agricolae uses the pooled MSE and
# residual df from the RBD factorial model; grouping letters ordered with mixedsort.
outA  <- agricolae::LSD.test(Ca, A,  DFerror = dferr, MSerror = MSE, alpha = alpha)
outB  <- agricolae::LSD.test(Ca, B,  DFerror = dferr, MSerror = MSE, alpha = alpha)
outAB <- agricolae::LSD.test(Ca, 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(Ca, 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)

# Four soil cations measured on every plot of the Contour x Depth RBD.
Y <- as.matrix(Soils[, c("Ca", "Mg", "K", "Na")])

res.man <- manova(Y ~ Block + A + B + A:B)
broom::tidy(res.man)                             # Pillai's trace per term
summary(res.man, test = "Pillai")                # full multivariate test table
effectsize::eta_squared(res.man)                 # partial multivariate eta-squared