# dataset1.csv : Location, FactorA, FactorB, Yield, Char1 ... Char6
data <- read.csv("dataset1.csv", header = TRUE, na.strings = c(""))

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

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 ANOVA per environment -> list of fitted models
models_by_loc <- lapply(split(data, data$Location), function(sub) {
  sub <- droplevels(sub)
  if (nlevels(factor(sub$FactorA)) < 2 || nlevels(factor(sub$FactorB)) < 2) return(NULL)
  aov(as.formula(paste(trait, "~ 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)
Y <- data_t[[trait]]

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

# mean squares = SS / df
MS <- atab[, 1] / atab[, 2]
names(MS) <- c("L", "A", "B", "AxB", "LxA", "LxB", "LxAxB", "Error")

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

dfABL <- atab[7, 2]     # df of L x A x B
dfP   <- atab[8, 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)