# Self-contained copy of dataset1 (Location x FactorA x FactorB x Block, balanced).
# In the app these columns come from the user's uploaded CSV.
data <- data.frame(
  Location = rep(c("A", "B"), each = 16),
  FactorA  = rep(rep(c("S1", "S2"), each = 8), times = 2),
  FactorB  = rep(rep(c("C1", "C2"), each = 4), times = 4),
  Block    = rep(1:4, times = 8),
  Yield    = c(1.35, 1.27, 0.93, 1.04, 1.39, 0.63, 0.78, 1.23,
               1.13, 1.64, 0.85, 0.68, 1.87, 0.88, 0.87, 0.83,
               1.67, 1.13, 1.31, 1.03, 1.43, 1.22, 1.02, 0.90,
               1.06, 0.77, 1.11, 1.29, 1.30, 0.83, 0.87, 0.91),
  stringsAsFactors = FALSE
)

data$Location <- as.factor(data$Location)   # A, B      (2 environments)
data$FactorA  <- as.factor(data$FactorA)    # S1, S2
data$FactorB  <- as.factor(data$FactorB)    # C1, C2
data$Block    <- as.factor(data$Block)      # 1..4      (blocks within each location)

trait <- "Yield"                            # analysis is repeated for every trait

# combination labels used by the interaction sources and mean comparisons
data$AxB   <- interaction(data$FactorA, data$FactorB, sep = "x")
data$AxL   <- interaction(data$FactorA, data$Location, sep = "x")
data$BxL   <- interaction(data$FactorB, data$Location, sep = "x")
data$AxBxL <- interaction(data$FactorA, data$FactorB, data$Location, sep = "x")

# one within-location RBD ANOVA per environment -> list of fitted models
models_by_loc <- lapply(split(data, data$Location), function(sub) {
  sub <- droplevels(sub)
  if (nlevels(factor(sub$Block))   < 2 ||
      nlevels(factor(sub$FactorA)) < 2 ||
      nlevels(factor(sub$FactorB)) < 2) return(NULL)
  aov(as.formula(paste(trait, "~ Block + FactorA + FactorB")), data = sub)
})
models_by_loc <- Filter(Negate(is.null), models_by_loc)

# Bartlett test across the per-location residuals
bt <- bartlett.test(models_by_loc)
bt$statistic   # chi-square
bt$p.value     # significance of variance heterogeneity

# per-location error mean square
mse_by_loc <- sapply(models_by_loc, function(m) sum(resid(m)^2) / df.residual(m))

# Aitken transform ONLY when Bartlett is significant (p <= 0.05)
data_t <- data
if (bt$p.value <= 0.05) {
  for (loc in names(mse_by_loc)) {
    rows <- data$Location == loc
    data_t[rows, trait] <- data[rows, trait] / sqrt(mse_by_loc[loc])
  }
}
# data_t holds the (possibly transformed) response used in the pooled ANOVA

L     <- factor(data_t$Location)
A     <- factor(data_t$FactorA)
B     <- factor(data_t$FactorB)
Block <- factor(data_t$Block)
Y     <- data_t[[trait]]

fit   <- lm(Y ~ L + L:Block + A + B + A:B + L:A + L:B + L:A:B)
atab  <- car::Anova(fit, type = "II")   # rows 1..8 = sources, row 9 = Residuals

# mean squares = SS / df, in the module's row order
MS <- atab[, 1] / atab[, 2]
names(MS) <- c("L", "A", "B", "LxBlock", "AxB", "LxA", "LxB", "LxAxB", "Error")

MSABL <- MS["LxAxB"]    # three-way interaction mean square
MSP   <- MS["Error"]    # pooled error mean square

dfABL <- atab[8, 2]     # df of L x A x B
dfP   <- atab[9, 2]     # df of Pooled Error

model <- "model_4"      # "model_1", "model_2", "model_3", "model_4" (default)

# choose the denominator MS and df for L, A and B under each pooling model
denom <- switch(model,
  # A, B, Location all random: test all three against L x A x B
  model_1 = list(L = c(MSABL, dfABL), A = c(MSABL, dfABL), B = c(MSABL, dfABL)),
  # factors fixed, location random: L vs Error; A, B vs L x A x B
  model_2 = list(L = c(MSP,  dfP),   A = c(MSABL, dfABL), B = c(MSABL, dfABL)),
  # location fixed, factors random: L vs L x A x B; A, B vs Error
  model_3 = list(L = c(MSABL, dfABL), A = c(MSP,  dfP),   B = c(MSP,  dfP)),
  # all fixed (default): everything vs Pooled Error
  model_4 = list(L = c(MSP,  dfP),   A = c(MSP,  dfP),   B = c(MSP,  dfP))
)

F_and_p <- function(ms, df_num, ms_err, df_err) {
  Fv <- ms / ms_err
  c(F = Fv, p = pf(Fv, df_num, df_err, lower.tail = FALSE))
}

resL <- F_and_p(MS["L"], atab[1, 2], denom$L[1], denom$L[2])
resA <- F_and_p(MS["A"], atab[2, 2], denom$A[1], denom$A[2])
resB <- F_and_p(MS["B"], atab[3, 2], denom$B[1], denom$B[2])

# interactions are always tested against the Pooled Error (car::Anova default F, p)
# star codes: ** p <= 0.01, * 0.01 < p <= 0.05, NS otherwise
star <- function(p) ifelse(p <= 0.01, "**", ifelse(p <= 0.05, "*", "NS"))

alpha <- 0.05
padj  <- "none"     # "none" (default), "bonferroni", "holm" or "BH"; LSD only

# error term for a main effect follows the pooling model chosen above;
# interaction sources use the Pooled Error (MSP, dfP)
outA <- agricolae::LSD.test(Y, A, DFerror = denom$A[2], MSerror = denom$A[1],
                            alpha = alpha, p.adj = padj)

statA <- outA$statistics        # includes CV and the CD (LSD) column
rA    <- outA$means[1, 3]       # replications behind each Factor A mean

CV_A  <- statA[1, 4]                       # coefficient of variation (%)
CD_A  <- statA[1, 6]                       # critical difference at alpha (p.adj = "none")
SEm_A <- sqrt(statA[1, 1] / rA)            # standard error of a mean
SEd_A <- sqrt(2 * statA[1, 1] / rA)        # standard error of a difference

# compact-letter grouping (ordered naturally); shown only when significant
groupsA <- outA$groups[gtools::mixedsort(rownames(outA$groups)), ]

# --- interaction example: A x B against the Pooled Error ---
AB   <- factor(data_t$AxB)
outAB <- agricolae::LSD.test(Y, AB, DFerror = dfP, MSerror = MSP,
                             alpha = alpha, p.adj = padj)