# MinutePlot validation -- Two-way ANOVA + Tukey letters (npk, N x P)
# Install once (skip on later runs):
#   install.packages(c("emmeans", "multcomp", "multcompView"))   # exactly what THIS script uses

# Load packages (every session):
library(emmeans); library(multcomp); library(multcompView)

# Versions (compare against those stated on this page):
cat(R.version.string, "\n")
for (p in c("emmeans", "multcomp", "multcompView")) cat(p, as.character(packageVersion(p)), "\n")

# Load data:
d <- npk                         # built-in; or: d <- read.csv("npk_r.csv"); d$N <- factor(d$N); d$P <- factor(d$P)

# Analysis (one print per table shown on the page, in page order):
fit <- aov(yield ~ N * P, data = d)
print(summary(fit))
print(multcomp::cld(emmeans(fit, ~ N), adjust = "tukey", Letters = letters, alpha = 0.05))
print(multcomp::cld(emmeans(fit, ~ P), adjust = "tukey", Letters = letters, alpha = 0.05))
print(multcomp::cld(emmeans(fit, ~ N * P), adjust = "tukey", Letters = letters, alpha = 0.05))