ANOVA in R: aov(), oneway.test() and Tukey

To run a one-way ANOVA in R, fit aov(outcome ~ group, data = ...) and print summary() of the fit. The summary gives the one-way ANOVA F value and its p-value for one question: are the group means all equal? A complete analysis adds three more steps: a check of the equal-variance assumption, an effect size with its interval, and a post hoc test such as TukeyHSD() that says which groups differ.

The worked example below runs all of it on one dataset that ships with R and reads each line of output. The numbers come from a deterministic Celeus engine run of the same analysis, and the R script that run generated is below the article. For the wider picture of when to choose which ANOVA, start from the guide to ANOVA tests and the overview of t-tests and their ANOVA extension.

The data and the question

The example uses PlantGrowth, a data frame in R’s datasets package from Dobson’s 1983 textbook. It records the dried weight of 30 plants: 10 under a control condition (ctrl) and 10 under each of two treatment conditions (trt1 and trt2). Each plant belongs to one condition and is weighed once, so the three groups are independent. That is a one-way between-groups design, the case aov() is built for.

The question is: does mean dried weight differ between the three conditions, and if it does, which conditions differ from which?

The design decides the function. Plants measured repeatedly under several conditions would need a repeated-measures model, and two grouping factors would need a two-way ANOVA. A single outcome measured once per unit in independent groups is the situation everything below assumes.

Look at the data first

Plot the individual values before any test. A boxplot with the raw points on top shows the centre and spread of each condition, any outlying plant, and whether one group is far more variable than the others, which is the first thing that would make the classic F a poor choice.

boxplot(weight ~ group, data = PlantGrowth, ylab = "Dried weight")
stripchart(weight ~ group, data = PlantGrowth, vertical = TRUE,
           method = "jitter", add = TRUE, pch = 16)
aggregate(weight ~ group, data = PlantGrowth,
          FUN = function(x) c(M = mean(x), SD = sd(x)))

The summary table from aggregate() belongs in the report as well, since the ANOVA output itself never prints the group means. A plot also catches data-entry mistakes that no test would: a weight typed with a misplaced decimal point shows up as an isolated point long before it shows up in an F value.

The whole script

The analysis is a handful of calls. aov(), oneway.test() and TukeyHSD() are in base R’s stats package, so nothing needs installing for them; the effect size with an interval needs the effectsize package from CRAN.

library(effectsize)

data <- PlantGrowth
data$group <- factor(data$group)   # already a factor here; needed if read from a CSV

# 1. Fit the model and print the ANOVA table
fit <- aov(weight ~ group, data = data)
summary(fit)

# 2. The same F test through oneway.test(), classic form
oneway.test(weight ~ group, data = data, var.equal = TRUE)

# 3. Check the equal-variance and normality assumptions
car::leveneTest(weight ~ group, data = data, center = median)   # needs the car package
shapiro.test(residuals(fit))

# 4. Effect size: eta squared with a two-sided 95% interval
eta_squared(fit, partial = FALSE, ci = 0.95, alternative = "two.sided")

# 5. Post hoc: Tukey's honestly significant difference
TukeyHSD(fit)

The engine’s own script, under the article, does the same work on a copy of the data read from a CSV file. It cleans the cell text, drops incomplete rows explicitly and gets the F from oneway.test() with var.equal = TRUE. The engine’s assumption checks are not in that script; they come from the engine run itself, per group instead of on the residuals, which is why they differ from the shapiro.test() call above. On the same package versions the headline numbers match.

Three ways to run the same F test

R offers three routes to the one-way F test, and two of them print the same ANOVA table. In a check under R 4.6.1 on this data, aov() and the anova() table of an lm() fit returned the same F value and p-value, because aov() is a wrapper that fits a linear model and arranges its sums of squares into an ANOVA table.

summary(aov(weight ~ group, data = PlantGrowth))        # classic F, ANOVA table
anova(lm(weight ~ group, data = PlantGrowth))           # the same table through lm()
oneway.test(weight ~ group, data = PlantGrowth, var.equal = TRUE)  # the same F, no table

The third call, oneway.test(), is a different function with a different default. Without var.equal = TRUE it does not assume equal variances: it runs Welch’s ANOVA (Welch, 1951), whose denominator degrees of freedom are fractional. In a check on this data, leaving the argument out gave a larger F and a denominator df that was not a whole number. The guide to ANOVA tests works through that Welch result on the same plants.

So the choice is not aov() against oneway.test() but classic F against Welch. Pick it before you look at the output, and write the choice into the call. The guide to one-way ANOVA, Welch’s ANOVA and Kruskal-Wallis explains how to decide, and Delacre and colleagues (2019) argue for Welch’s version by default. This post uses the classic F because the next sections read aov() output, and TukeyHSD() needs an aov object.

Reading the summary(aov) output line by line

The first step printed the following when the script ran in base R:

            Df Sum Sq Mean Sq F value Pr(>F)  
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Df counts the independent pieces of information in each row. The group row has 2, the number of groups minus one, and the Residuals row has 27, the 30 plants minus the 3 group means that were estimated. If the group row shows 1 instead of 2, the group column is numeric: R fitted a straight line through the codes 1, 2 and 3, not a comparison of three means. Wrap the column in factor().

Sum Sq splits the total variation in weight into the part explained by the group means and the part left within the groups. Mean Sq divides each by its degrees of freedom, which turns the sums into variance estimates. The group row estimates how far the three condition means spread out, and the Residuals row estimates the plant-to-plant variance inside conditions.

F value is the ratio of those two, here 4.846. If the three conditions had the same population mean, the group mean square would be about as large as the residual mean square and F would sit near 1. The engine reports the same F, 4.846, on 2 and 27 degrees of freedom.

Pr(>F) is the p-value: the probability of an F at least this large if all three population means were equal and the model’s assumptions held. It is 0.0159, so the data are poorly compatible with equal means. The p-value is not the probability that the means are equal, and it says nothing about how large the differences are.

The stars and the Signif. codes line are decoration. Report the p-value itself, not the asterisks.

What the table does not show is which groups differ, the direction of any difference, the group means, or the sample size. A significant F says at least one mean differs from the others. The post hoc step and the effect size below supply the rest, and the group means come from aggregate(weight ~ group, data = PlantGrowth, FUN = mean). Count the rows R used: aov() drops rows with a missing outcome or group without a message in a default session, so compare nrow(na.omit(data[, c("weight", "group")])) with the plants you expect. The engine’s record shows n = 30 with no rows dropped.

Reading the oneway.test() output

The classic form of the same test prints a shorter block:

	One-way analysis of means

data:  weight and group
F = 4.8461, num df = 2, denom df = 27, p-value = 0.01591

The title reads “One-way analysis of means”. With var.equal = FALSE, the default, it reads “One-way analysis of means (not assuming equal variances)”, and that line is the quickest way to see which test ran. The data line names your outcome and group, so check that they are the right way round. The numbers line carries the same F, the numerator df of 2 and the denominator df of 27 as the table above, with the p-value to more digits.

The output has no sums of squares, so it cannot give you an effect size on its own. That is a practical reason to prefer aov() for the classic F: the pieces the effect size needs are in the fitted object.

Checking the assumptions the F test relies on

The classic F assumes independent observations, residuals that are roughly normal, and equal variances across groups. Independence comes from the design, not from a test: each plant is a separate unit weighed once.

The engine ran the other two checks on the data. A Shapiro-Wilk test in each condition found nothing to flag (p = 0.747 for ctrl, 0.452 for trt1 and 0.564 for trt2), and a median-centred Levene test on the variances also found nothing (statistic 1.119, p = 0.341). In base R, shapiro.test(residuals(fit)) pools the residuals and runs one normality test, while car::leveneTest() is the usual Levene call and bartlett.test() is a base R alternative that is more sensitive to non-normality.

Read these as diagnostics, not as a verdict. With 10 plants per group both tests have little power, so a clean result is not evidence that the weights are normal or the variances equal. It only means the tests found nothing to flag. Rochon and colleagues (2012) show that choosing a test from a preliminary test distorts the error rate of what follows, in the two-sample case they studied. Choose the test from the design and say in the plan which one is primary. The post on homogeneity of variance covers what unequal spreads do to the classic F, and the parametric test assumptions post walks through the checks and their fallbacks.

Effect size: eta squared with an interval

Eta squared is the between-groups sum of squares divided by the total sum of squares: the share of the variation in weight that the grouping accounts for. With one factor it equals the R-squared of the equivalent lm() fit, and it is the same as partial eta squared, a distinction that matters only once a second factor appears. The post on eta squared explains how the versions diverge.

For this example the engine reported 0.2641 with a 95% interval from 0.01099 to 0.4902. You can check the point estimate by hand: divide the group Sum Sq in the table above by the total of the two Sum Sq values, and the result matches to the precision of the printed sums.

# Effect Size for ANOVA (Type I)

Parameter | Eta2 |       95% CI
-------------------------------
group     | 0.26 | [0.01, 0.49]

“Type I” in the header means sequential sums of squares, the kind aov() reports. With a single factor there is nothing to be sequential about, so Type I, II and III sums of squares give the same table; the choice only matters once a model has two or more factors and unequal group sizes.

The condition accounts for about a quarter of the variation in this sample, which is a large effect by common benchmarks. The interval is the more honest summary. It runs from compatible-with-almost-nothing to compatible-with-about-half, which is what 30 plants can say about size.

Two cautions apply. The interval comes from the noncentral F distribution and assumes the model’s assumptions hold, so treat it as approximate when the variances differ. And eta squared computed from a sample tends to overstate the population value; Lakens (2013) explains the less biased omega squared. The engine asked for a two-sided interval on purpose, and the script above does the same with alternative = "two.sided", so set it explicitly rather than rely on a package default. The post on effect sizes covers how the mean-based measures relate to the rank-based ones.

Which groups differ: TukeyHSD()

A significant omnibus F calls for a pairwise follow-up. Tukey’s honestly significant difference (Tukey, 1949) compares every pair of group means using the pooled within-group variance from the ANOVA, and adjusts for all the comparisons together, so the family of intervals has the stated coverage. TukeyHSD() takes the aov object, which is why oneway.test() is no use here. Below its header lines, the base R printout shows:

$group
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

The row label reads “first minus second” in the order written. The row trt2-trt1 is the mean of trt2 minus the mean of trt1, so a positive difference means the first-named group is heavier. R orders the levels alphabetically, so the order can surprise you; set it with relevel() or factor(levels = ...) before fitting if you want the control as the baseline.

diff, lwr and upr are the difference in means and its simultaneous 95% interval. p adj is the adjusted p-value. They agree with each other: the interval excludes zero exactly when the adjusted p-value is below 0.05.

Only one pair separates after adjustment: trt2 against trt1. Neither treatment differs clearly from the control. The engine summarised the same result as compact letters, ctrl (ab), trt1 (a) and trt2 (b); groups that share a letter are not significantly different at the adjusted threshold, which is not the same as being shown equal.

The Hedges’ g column needs care. The g interval for trt2 against ctrl, 0.01252 to 1.793, just excludes zero, while its adjusted p-value is 0.198. They do not conflict. Each g interval is computed for its pair alone with no adjustment for the other two comparisons, so base the claim about which pairs differ on the adjusted results and label the g intervals as unadjusted. The post on post hoc tests after ANOVA compares Tukey with Games-Howell and Dunn’s test, and the guide to multiple comparisons covers when other corrections fit.

What this study could detect

The engine also reports a sensitivity figure, which describes the design and not the result. With three groups of 10 at the 0.05 level, the study had 80% power only for effects of Cohen’s f = 0.5997 or larger, a very large effect, and the engine flags it as underpowered. Detecting a conventional medium effect of f = 0.25 would need 53 plants per group.

The sensitivity figure is not “observed power”, which only restates the p-value. In practice it means the non-significant comparisons against the control are inconclusive, not evidence of no difference, and the significant one is more likely to overstate the true gap than understate it.

When another test fits better

A one-way ANOVA compares means. When the outcome is ordinal, badly skewed or a count, a rank-based comparison may suit the question better, and the post on the Kruskal-Wallis test in R runs that route with the same structure: omnibus test, effect size, pairwise follow-up. When a measured covariate, such as a baseline value, explains part of the outcome, the guide to ANCOVA shows how adjusting for it changes the comparison. Neither is a repair for a failed assumption check; each answers a different question, so choose it from the question and the design.

Pitfalls and silent behaviours

A numeric group column. Codes such as 1, 2 and 3 make aov() fit a slope with one degree of freedom for the group. The F value and p-value then answer a different question. Check that the Df of the group row is the number of groups minus one.

Missing values. Rows with a missing outcome or group are dropped without a message in a default session (R Core Team, 2026). Count the complete rows yourself and report that n.

pairwise.t.test() is a different follow-up. Base R’s pairwise.t.test() with its default pooled standard deviation runs a t test per pair with a p-value adjustment you choose, Holm by default. It is a valid procedure, but it is not Tukey’s, and its adjusted p-values differ. Name the one you ran.

Unequal variances after the fact. If the spreads differ a lot, especially with unequal group sizes, switch the plan to Welch’s ANOVA and Games-Howell, and report the classic F as a sensitivity analysis rather than quietly replacing one by the other.

Several outcomes. Each outcome is its own ANOVA, and running many of them is a multiplicity problem of its own.

Reading F as an effect size. A large F in a large sample can come from a tiny difference. Always report the effect size and its interval.

Reporting the result

Rounded to the usual conventions, the finding reads:

Mean dried weight differed across the three conditions, F(2, 27) = 4.85, p = .016, eta squared = .26, 95% CI [.01, .49]. Tukey’s HSD showed that plants under the second treatment were heavier than plants under the first, a difference of 0.87, 95% CI [0.17, 1.56], adjusted p = .012. Neither treatment differed clearly from the control.

Put the group means and standard deviations in front of that sentence or in a table. The post on reporting ANOVA results in APA style has fill-in sentences for the other ANOVA designs, and setting seeds and pinning package versions covers the reproducibility side. The example is research statistics, not clinical guidance.

Rerunning the example

The R script and the engine’s output sit under the article, along with the seed. The script reads its data from data/anova-in-r.csv next to it; dir.create("data") followed by write.csv(PlantGrowth, "data/anova-in-r.csv", row.names = FALSE) creates that file from the copy that ships with R. It also needs the readr, effectsize and pwr packages. On the same package versions, the classic F, the interval for eta squared, the Tukey table and the sensitivity figure reproduce to the digits shown.

# Results

- Test: One-way ANOVA
- n: 30
- statistic: 4.846
- p-value: 0.0159
- effect size: eta-squared = 0.2641 (95% CI 0.01099 to 0.4902)
  (two-sided CI at the run's confidence level)
- seed: 20261111

## Assumption checks
- normality by group (Shapiro-Wilk): passed=TRUE
- variance homogeneity (Levene): passed=TRUE (no evidence of unequal variances)

## Post-hoc (Tukey HSD; correction: single-step (Tukey))
- trt1 vs ctrl: adj p = 0.3909 | diff -0.371 [-1.062, 0.3202] | Hedges' g -0.5102 [-1.359, 0.3521]
- trt2 vs ctrl: adj p = 0.198 | diff 0.494 [-0.1972, 1.185] | Hedges' g 0.9139 [0.01252, 1.793]
- trt2 vs trt1: adj p = 0.01201 | diff 0.865 [0.1738, 1.556] | Hedges' g 1.289 [0.3371, 2.213]

groups: ctrl (ab), trt1 (a), trt2 (b) - groups sharing a letter are not significantly different at the adjusted threshold.

## Power (sensitivity, not observed power)
- minimum detectable effect at n=30: Cohen's f = 0.5997 (severely underpowered)
- to detect a conventional medium effect: 53 per group to detect f=0.25
- alpha = 0.05, target power = 0.8

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: One-way ANOVA  |  engine 0.53.2  |  catalog c40-dc08
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261111)
.data_path <- "data/anova-in-r.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("weight", "group")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- oneway.test(`weight` ~ `group`, data = data, var.equal = TRUE)
print(result)

# Effect size (Eta2) and its confidence interval.
effect <- { .es_d <- data; .es_d[["group"]] <- as.factor(.es_d[["group"]]); .es_d[["weight"]] <- .es_d[["weight"]] - stats::median(.es_d[["weight"]], na.rm = TRUE); effectsize::eta_squared(stats::aov(`weight` ~ `group`, data = .es_d), partial = FALSE, ci = 0.95, alternative = "two.sided") }
print(effect)

# Post-hoc pairwise comparisons (Tukey HSD; correction: single-step (Tukey)).
data[["group"]] <- as.factor(data[["group"]])
posthoc <- { .ph_d <- data.frame(.y = data[["weight"]], .g = as.factor(data[["group"]])); .ph_tk <- stats::TukeyHSD(stats::aov(.y ~ .g, data = .ph_d)); names(.ph_tk) <- "group"; .ph_tk }
print(posthoc)

# Sensitivity: minimum detectable effect at n=30 (power=0.8, alpha=0.05). NOT observed power.
power <- pwr::pwr.anova.test(k = 3, n = 10, sig.level = 0.05, power = 0.8)
print(power)