Quick summary
Summarize this blog with AI
Introduction
One-way ANOVA tests whether all group population means are equal. A small p-value indicates evidence against that joint hypothesis; it does not identify every differing pair or measure practical importance. This complete example uses three fictional independent groups and keeps one fitted model throughout.
Tested with R 4.5.1. Run the code blocks in order. Install the named packages before starting; inputs are supplied below.
Supply the data and define the analysis rows
Each row is one independent unit, with a numeric response and a categorical group. Repeated measurements from one person are not independent units. Missing rows are counted before removal; the decision to analyze complete cases must be justified in a real study.
raw <- data.frame(group = factor(rep(c("A", "B", "C"), each = 5)),
score = c(8, 9, 10, 11, 12, 10, 11, 12, 13, 14,
14, 15, 16, 17, 18))
keep <- complete.cases(raw[c("score", "group")])
cat("Excluded rows:", sum(!keep), "
")
dat <- droplevels(raw[keep, ])
stopifnot(is.numeric(dat$score), all(is.finite(dat$score)),
nlevels(dat$group) >= 2, all(table(dat$group) >= 2))
print(aggregate(score ~ group, dat, mean))
No rows are excluded; the group means are A = 10, B = 12 and C = 16. The response’s units should accompany these estimates in a report.
Fit the model and verify the F statistic
In aov(score ~ group), group is a factor, so the model estimates a separate group mean. The classical F statistic divides the between-group mean square by the within-group mean square. It uses 2 and 12 degrees of freedom here.
fit <- aov(score ~ group, data = dat, na.action = na.fail)
tab <- summary(fit)[[1]]
print(tab)
group_mean <- ave(dat$score, dat$group, FUN = mean)
ss_between <- sum((group_mean - mean(dat$score))^2)
ss_within <- sum((dat$score - group_mean)^2)
f_manual <- (ss_between / 2) / (ss_within / 12)
stopifnot(abs(f_manual - tab[1, "F value"]) < 1e-12,
abs(tab[1, "F value"] - 18.6666666667) < 1e-8)
cat("F:", round(f_manual, 3), "p:", signif(tab[1, "Pr(>F)"], 4), "
")
F = 18.667 and p is approximately 0.000207. Under the model this is evidence that at least one population mean differs. It does not mean the null hypothesis has that probability.
Inspect assumptions without turning checks into certificates
The classical model assumes independent errors with a common variance and normally distributed errors within groups for exact small-sample inference. Normality does not protect against outliers. Independence comes from design; a nonsignificant diagnostic test does not prove the assumptions. Inspect group distributions and residual patterns, including data collection order when relevant.
draw_anova_checks <- function() {
old <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
on.exit(par(old))
boxplot(score ~ group, data = dat, xlab = "Group", ylab = "Score",
main = "Within-group spread", col = "#b6d9e6")
qqnorm(residuals(fit), main = "Residual Q-Q plot", pch = 19, cex = 0.7)
qqline(residuals(fit), col = "#b33d24")
}
draw_anova_checks()
print(tapply(dat$score, dat$group, var))
All three sample variances are 2.5 in this deliberately simple fixture. The plot illustrates the checks; fifteen constructed observations do not demonstrate that a real study meets its assumptions.
Use Tukey comparisons for the classical model
TukeyHSD() reports simultaneous intervals and adjusted p-values for pairwise differences. Its familywise interpretation belongs to the model assumptions; it is not a guarantee about arbitrary searches across analyses.
pairs <- TukeyHSD(fit, "group")$group
print(round(pairs, 4))
stopifnot(isTRUE(all.equal(unname(pairs[, "diff"]), c(2, 6, 4))),
pairs["B-A", "lwr"] < 0, pairs["B-A", "upr"] > 0,
pairs["C-A", "lwr"] > 0, pairs["C-B", "lwr"] > 0)
The observed B − A, C − A and C − B differences are 2, 6 and 4 points. The B − A interval includes zero; the other two do not. Report the intervals, not just a significant/nonsignificant label.
Unequal variances and invalid groups
Welch’s one-way test in oneway.test() allows unequal variances. It does not solve dependence or every outlier problem. Do not append classical Tukey intervals to a Welch analysis as if they share its variance model. Kruskal–Wallis is a rank-based alternative with its own assumptions and is not generally a test of equal means.
welch <- oneway.test(score ~ group, data = dat, var.equal = FALSE)
cat("Welch p:", signif(welch$p.value, 4), "
")
check_groups <- function(d) {
if (!is.numeric(d$score) || any(!is.finite(d$score)) ||
anyNA(d$group) || length(unique(d$group)) < 2 ||
any(table(droplevels(factor(d$group))) < 2))
stop("Need finite scores and at least two observations in each of two groups")
invisible(TRUE)
}
fails <- function(expr) inherits(tryCatch(expr, error = identity), "error")
stopifnot(check_groups(dat), fails(check_groups(dat[FALSE, ])),
fails(check_groups(dat[1:5, ])),
fails(check_groups(transform(dat, score = NA_real_))))
Plan sample size from a meaningful difference
For a balanced future study, power.anova.test() uses the variance of anticipated group means and the within-group variance. Those inputs should come from substantive planning or external evidence. This illustrative plan targets means 10, 12 and 14 with within-group SD 3, 80% power and a 5% test level.
planned <- power.anova.test(groups = 3, between.var = var(c(10, 12, 14)),
within.var = 3^2, sig.level = 0.05, power = 0.8)
n_each <- ceiling(planned$n)
cat("Units per group:", n_each, "total:", 3 * n_each, "
")
achieved <- power.anova.test(groups = 3, n = n_each,
between.var = 4, within.var = 9, sig.level = 0.05)$power
stopifnot(achieved >= 0.8)
The result rounds up to 12 units per group, or 36 total. Allow separately for expected attrition. This calculation assumes balanced independent groups and the classical variance model; observed post-hoc power is not extra evidence for a fitted result.
