# MinutePlot validation -- Three-way ANOVA -- unbalanced (npk with rows removed)
# Install once (skip on later runs):
#   install.packages(c("car", "emmeans", "multcomp", "multcompView"))   # exactly what THIS script uses

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

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

# Load data:
d <- read.csv("npk_unbalanced_r.csv")   # same folder as script: npk minus rows 1, 6, 15 and 24 of the hosted npk.csv
d$N <- factor(d$N); d$P <- factor(d$P); d$K <- factor(d$K)

# Analysis (one print per table shown on the page, in page order):
fit <- aov(yield ~ N * P * K, data = d)
print(car::Anova(fit, type = 2))                  # Type II sums of squares -- MinutePlot's type for factorial ANOVA
print(multcomp::cld(emmeans(fit, ~ N * P * K), adjust = "tukey", Letters = letters, alpha = 0.05))