data(CO2, package = "datasets")

# CO2 is a real split-plot-style trial: 12 plants (2 Types x 2 Treatments x 3
# replicates), each measured at 7 CO2 concentrations. Restricting to four
# concentrations and reading them as a 2 x 2 factorial gives a perfectly
# balanced pooled split-plot (2,1):
#   Location = Type       (2: Quebec, Mississippi)
#   Block    = replicate  (3, nested within Location)
#   FactorA  = Treatment  (2: nonchilled, chilled)   - main plot 1
#   FactorB  = conc range (2: low {95,175}, high {350,675}) - main plot 2
#   FactorC  = conc within range (2: {95,350}, {175,675})   - sub plot
d <- subset(CO2, conc %in% c(95, 175, 350, 675))
d <- transform(
  d,
  Location = factor(Type, labels = c("Quebec", "Mississippi")),
  Block    = factor(substr(as.character(Plant), 3, 3)),
  FactorA  = factor(Treatment, labels = c("A_nonchilled", "A_chilled")),
  FactorB  = factor(ifelse(conc <= 175, "B_low", "B_high"),
                    levels = c("B_low", "B_high")),
  FactorC  = factor(ifelse(conc %in% c(95, 350), "C1", "C2")),
  Y        = uptake
)
d <- d[, c("Location", "Block", "FactorA", "FactorB", "FactorC", "Y")]

nrow(d)                                                    # 48
all(table(d$Location, d$Block, d$FactorA, d$FactorB, d$FactorC) == 1)  # TRUE

# The app renames the user's columns to exactly these five structural names,
# then builds the interaction factors it needs:
L  <- d$Location; Block <- d$Block
A  <- d$FactorA;  B <- d$FactorB; C <- d$FactorC
AB <- interaction(A, B, sep = "x", lex.order = TRUE)

data_temp <- data.frame(Y = d$Y, L = L, Block = Block, A = A, B = B, C = C, AB = AB)

# Whole-plot fit: the LAST term, L:Block:A, is Error(a).
# Block:A is deliberately absent from the model, so this term absorbs both
# Block:A and L:Block:A  ->  df_Ea = nLoc * (nBlock - 1) * (nA - 1)
model_wp <- lm(Y ~ L + Block + L:Block + A + B + A:B + L:A + L:B + L:A:B + L:Block:A,
               data = data_temp)
aov_wp   <- anova(model_wp)
rownames(aov_wp)[10]                       # "L:Block:A"
df_Ea <- aov_wp[10, "Df"];  MS_Ea <- aov_wp[10, "Mean Sq"]

# Block-first fit: its Residuals row (the L:Block:AB pooling) is Error(b).
#   df_Eb = nLoc * nA * nB * (nBlock - 1) * (nC - 1)
model_block <- lm(Y ~ Block + L + L:Block + AB + L:AB + L:Block:AB + C + AB:C + L:AB:C,
                  data = data_temp)
aov_block <- anova(model_block)
df_Eb <- aov_block["Residuals", "Df"];  MS_Eb <- aov_block["Residuals", "Mean Sq"]

c(df_Ea = df_Ea, MS_Ea = MS_Ea, df_Eb = df_Eb, MS_Eb = MS_Eb)
#> df_Ea = 4  MS_Ea = 19.368   df_Eb = 16  MS_Eb = 3.959
# matching 2*(3-1)*(2-1) = 4 and 2*2*2*(3-1)*(2-1) = 16

# Main-plot partition: A, B, A:B
aov_main <- anova(lm(Y ~ A + B + A:B, data = data_temp))

# Sub-plot / higher-order partition: C, A:C, B:C, A:B:C, C:L, A:C:L, B:C:L, A:B:C:L
aov_sub  <- anova(lm(Y ~ A * B * C * L, data = data_temp))

# Each effect is tested against ITS OWN stratum - never a pooled residual.
f_and_p <- function(SS, df_eff, MS_error, df_error) {
  MS <- SS / df_eff
  Fv <- MS / MS_error
  c(MS = MS, F = Fv, p = pf(Fv, df_eff, df_error, lower.tail = FALSE))
}

# Main plot A -> Error(a)
f_and_p(aov_main["A", "Sum Sq"], aov_main["A", "Df"], MS_Ea, df_Ea)
#> MS = 460.660  F = 23.785  p = 0.0082

# Sub plot C -> Error(b)
f_and_p(aov_sub["C", "Sum Sq"], aov_sub["C", "Df"], MS_Eb, df_Eb)
#> MS = 383.635  F = 96.893  p < 0.0001

# Location MAIN effect -> Error(b), NOT Error(a). Location's INTERACTIONS with
# the main-plot factors (L:A, L:B, L:A:B) are tested against Error(a) instead.
f_and_p(aov_wp[1, "Sum Sq"], aov_wp[1, "Df"], MS_Eb, df_Eb)
#> MS = 1396.442  F = 352.693  p < 0.0001

# One lm per Location, in BOTH strata. The main-plot model is fitted to cell
# means averaged over FactorC; the sub-plot model uses the raw observations.
per_loc <- lapply(split(d, d$Location), function(sub) {
  sub <- droplevels(sub)
  m_b <- lm(Y ~ Block + FactorA * FactorB * FactorC, data = sub)
  cm  <- aggregate(Y ~ Block + FactorA + FactorB, data = sub, FUN = mean)
  m_a <- lm(Y ~ Block + FactorA * FactorB, data = cm)
  if (df.residual(m_a) <= 0 || df.residual(m_b) <= 0) return(NULL)
  list(m_a = m_a, m_b = m_b)
})
valid <- Filter(Negate(is.null), per_loc)          # needs >= 2 Locations

test_a <- bartlett.test(lapply(valid, `[[`, "m_a"))  # chisq_Ea / pval_Ea
test_b <- bartlett.test(lapply(valid, `[[`, "m_b"))  # chisq_Eb / pval_Eb
c(pval_Ea = test_a$p.value, pval_Eb = test_b$p.value)
#> pval_Ea = 0.4148   pval_Eb = 0.9446   -> homogeneous, no transformation

# Aitken fires on a HARDCODED 0.05, not the user's alpha, and on EITHER stratum.
if (test_a$p.value <= 0.05 || test_b$p.value <= 0.05) {
  mse_b <- sapply(valid, function(v) sum(resid(v$m_b)^2) / df.residual(v$m_b))
  for (loc in names(mse_b)) {
    rows <- d$Location == loc
    d$Y[rows] <- d$Y[rows] / sqrt(mse_b[loc])      # divide by sqrt(MSE_b)
  }
}

library(agricolae); library(gtools)
alpha <- 0.05                                   # effective shipped default

# Main-plot factor A -> Error(a)
outA <- LSD.test(data_temp$Y, A, DFerror = df_Ea, MSerror = MS_Ea, alpha = alpha)
outA$groups[mixedsort(rownames(outA$groups)), ]
outA$statistics[1, 6]                           # CD = 3.5273  (column 6 for LSD)

# Sub-plot factor C -> Error(b)
outC <- LSD.test(data_temp$Y, C, DFerror = df_Eb, MSerror = MS_Eb, alpha = alpha)
outC$statistics[1, 6]                           # CD = 1.2177

# Tukey: the CD equivalent is in column 5, not 6
HSD.test(data_temp$Y, A, DFerror = df_Ea, MSerror = MS_Ea, alpha = alpha)$statistics[1, 5]

# DMRT: tabulated values and critical range values (only produced under DMRT)
dun <- duncan.test(data_temp$Y, A, DFerror = df_Ea, MSerror = MS_Ea, alpha = alpha)
dun$duncan[, 1]                                 # tabulated (studentized range) = 3.9265
dun$duncan[, 2]                                 # critical range               = 3.5273

# Reported dispersion statistics, all from the same agricolae object
r <- outA$means[1, 3]                           # replications per mean
sqrt(outA$statistics[1, 1] / r)                 # SE(m)  = 0.8983
sqrt(2 * outA$statistics[1, 1] / r)             # SE(d)  = 1.2704
outA$statistics[1, 4]                           # CV(%)  = 18.1184

# Partial eta-squared uses ONLY the effect's own error stratum in the denominator.
SS_A  <- aov_main["A", "Sum Sq"]
SS_Ea <- aov_wp[10, "Sum Sq"]

eta_sq_A <- SS_A / (SS_A + SS_Ea)
cohens_f <- sqrt(eta_sq_A / (1 - eta_sq_A))
round(cohens_f, 2)                              # 2.44, rounded to split_digit (default 2)

library(factoextra); library(scales)

# PCA needs at least two response variables. Uptake and a derived efficiency
# measure (uptake per unit concentration) are used here.
p <- subset(CO2, conc %in% c(95, 175, 350, 675))
resp <- data.frame(
  LxAxBxC    = interaction(d$Location, d$FactorA, d$FactorB, d$FactorC,
                           sep = "x", lex.order = TRUE),
  Uptake     = p$uptake,
  Efficiency = p$uptake / p$conc * 100
)

# Means per L x A x B x C treatment combination, row-ordered by mixedorder().
means <- aggregate(cbind(Uptake, Efficiency) ~ LxAxBxC, data = resp, FUN = mean)
rownames(means) <- means$LxAxBxC
means <- means[mixedorder(rownames(means)), -1, drop = FALSE]   # 16 combinations

res.pca <- prcomp(means, center = TRUE, scale = TRUE)   # correlation PCA via SVD

# The app reverses both sign conventions - reproduce this to match its output.
res.pca$x        <- -res.pca$x
res.pca$rotation <- -res.pca$rotation

get_eigenvalue(res.pca)          # PC1 = 1.185 (59.23%), PC2 = 0.815 (40.77%)
res.pca$rotation                 # loadings
fviz_eig(res.pca, addlabels = TRUE)
fviz_pca_biplot(res.pca, axes = c(1, 2), repel = TRUE)

index1 <- res.pca$x[, 1]; index2 <- res.pca$x[, 2]
cbind(`Index Score` = index1, `Scaled Index` = rescale(index1, to = c(0, 1)))
cbind(`Index Score` = index2, `Scaled Index` = rescale(index2, to = c(0, 1)))

round(cor(means, use = "pairwise.complete.obs"), 2)
PerformanceAnalytics::chart.Correlation(means, histogram = TRUE, pch = 19)

z <- na.omit(means$Uptake); n <- length(z)
if (n >= 3 && n <= 5000 && var(z) > 0) shapiro.test(z)
if (n >= 7 && var(z) > 0)              nortest::ad.test(z)
if (n >= 2 && var(z) > 0)              ks.test(z, "pnorm", mean(z), sd(z))