# MinutePlot validation -- Two-way ANOVA -- unbalanced (ToothGrowth 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("toothgrowth_unbalanced_r.csv")   # same folder as script: ToothGrowth minus rows 1-4 and 55-56 of the hosted toothgrowth.csv
d$dose <- factor(d$dose); d$supp <- factor(d$supp)

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