STAT 516 Lec 05

One-way analysis of variance (ANOVA)

Author

Karl Gregory

Published

October 5, 2026

Rust inhibitors example

Data from Kutner et al. (2005).

Ten experimental units assigned to each of four brands of rust inhibitors.

link <- url("https://people.stat.sc.edu/gregorkb/data/KNNLrust.txt")
rust <- read.csv(link,col.names=c("score","brand","rep"),sep = "", header = FALSE)
head(rust)

  score brand rep
1  43.9     1   1
2  39.0     1   2
3  46.7     1   3
4  43.8     1   4
5  44.2     1   5
6  47.7     1   6

Do the brands differ in effectiveness? Is there a best brand?

Visually compare the treatment group means with boxplots or dotplots.

boxplot(score ~ brand, data = rust)

stripchart(score ~ brand, data = rust, vertical = T)

stripchart(score ~ brand, data = rust, vertical = T, method = "jitter")

aggregate(score ~ brand, mean, data = rust) # mean of each treatment group

  brand score
1     1 43.14
2     2 89.44
3     3 67.95
4     4 40.47

aggregate(score ~ brand, sd, data = rust) # standard deviation

  brand    score
1     1 3.000074
2     2 2.218207
3     3 2.168589
4     4 2.436322

aggregate(score ~ brand, length, data = rust) # number of replications

  brand score
1     1    10
2     2    10
3     3    10
4     4    10

 brand    1    2    3    4    5    6    7    8    9   10  mean    sd  n
     1 43.9 39.0 46.7 43.8 44.2 47.7 43.6 38.9 43.6 40.0 43.14 3.000 10
     2 89.8 87.1 92.7 90.6 87.7 92.4 86.1 88.1 90.8 89.1 89.44 2.218 10
     3 68.4 69.3 68.5 66.4 70.0 68.1 70.6 65.2 63.8 69.2 67.95 2.169 10
     4 36.2 45.2 40.7 40.5 39.3 40.3 43.2 38.7 40.9 39.7 40.47 2.436 10

Randomized experiments comparing treatments

Start with $N$ experimental units (EUs), e.g. subjects, mice, etc.

Randomly assign each EU to one of $a$ treatment groups.

Measure on each EU after treatment a response $Y$.

Compute the average of the responses in each treatment group

Questions we’d like to answer:

Is the response mean the same in all treatment groups?
If not, then which pairs of means are different?

One-way treatment effects model

Consider the model

\[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where

$Y_{ij}$ is the response for EU $j$ in treatment group $i$.
$\mu$ represents an overall or baseline mean.
$\tau_i$ is the treatment effect for treatment $i$.
The $\varepsilon_{ij}$ are independent $\text{Normal}(0,\sigma^2)$ error terms.
The $n_i$ are the numbers of replicates in the treatment groups.

Of central interest are the hypotheses \[ \text{$H_0$: $\tau_1 = \dots = \tau_a$ \quad versus \quad $H_1$: $\tau_i \neq \tau_{i'}$ for some $i \neq i'$.} \] If we reject $H_0$, we may wish to sort/compare the treatments.

Identifiability constraint in the treatment effects model

The model has $a+1$ parameters to describe $a$ treatment means.

To identify $\mu$, $\tau_1,\dots,\tau_a$ uniquely, we typically set $\tau_1 = 0$.

Alternative “cell means” setup

An alternate version of the model is

\[ Y_{ij} = \mu_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where

$Y_{ij}$ is the response for EU $j$ in treatment group $i$.
$\mu_i$ represents the mean of treatment group $i$.
The $\varepsilon_{ij}$ are error terms distributed as $\text{Normal}(0,\sigma^2)$.

In this version of the model the central hypotheses become \[ \text{$H_0$: $\mu_1 = \dots = \mu_a$ \quad versus \quad $H_1$: $\mu_i \neq \mu_i'$ for some $i \neq i'$}. \]

Goals in one-way ANOVA

In the one-way treatment effects model \[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where $\varepsilon_{ij}\overset{\operatorname{ind}}{\sim}\text{Normal}(0,\sigma^2)$, we wish to

Visualize the data.
Estimate the parameters $\mu,\tau_1,\dots,\tau_a$.
Estimate the error term variance $\sigma^2$.
Decompose the variation in the $Y_{ij}$ as signal plus noise.
Test whether there is any difference in treatment group means.
Sort/compare the treatment means if there is any difference.
Check whether the model assumptions are satisfied.

Treatment effect estimation in one-way ANOVA

For each $i = 1,\dots,a$ define the observed treatment group mean as \[ \bar Y_{i.} = \frac{1}{n_i}\sum_{j=1}^{n_i}Y_{ij}. \]
Then, setting $\tau_1 = 0$, estimate $\mu$ and $\tau_2,\dots,\tau_a$ as \[ \hat \mu = \bar Y_{1.} \quad \text{ and } \quad \hat \tau_i = \bar Y_{i.} - \bar Y_{1.} \quad \text{ for } i = 2,\dots,a. \]
So treatment group 1 is regarded as a baseline, where:
1. The baseline has estimated mean $\hat \mu$.
2. The estimates $\hat \tau_2,\dots,\hat \tau_a$ are deviations from the baseline.

One obtains the fitted values $\hat Y_{ij} = \hat \mu + \hat \tau_i = \bar Y_{i.}$ for $i = 1,\dots,a$.

Rust inhibitors example (cont)

Use lm() with as.factor() to fit the one-way ANOVA model.

# use as.factor() to designate brand as a "factor"
lm_out <- lm(score ~ as.factor(brand), data = rust)
lm_out


Call:
lm(formula = score ~ as.factor(brand), data = rust)

Coefficients:
      (Intercept)  as.factor(brand)2  as.factor(brand)3  as.factor(brand)4  
            43.14              46.30              24.81              -2.67

See how $\hat \mu,\hat \tau_2,\hat \tau_3,\hat\tau_4$ are related to $\bar Y_{1.},\bar Y_{2.},\bar Y_{3.},\bar Y_{4.}$.

# compute the group means
aggregate(score ~ brand,mean, data = rust)

  brand score
1     1 43.14
2     2 89.44
3     3 67.95
4     4 40.47

Estimation of the error term variance $\sigma^2$

As in linear regression, define the

fitted values $\hat Y_{ij}$ as $\hat Y_{ij} = \bar Y_{i.}$ for $j = 1,\dots,n_i$, and the
residuals $\hat \varepsilon_{ij}$ as $\hat \varepsilon_{ij} = Y_{ij} - \bar Y_{i.}$

for $j=1,\dots,n_i$, $i = 1,\dots,a$.

Then an unbiased estimator of $\sigma^2$ is given by \[ \hat \sigma^2 = \frac{1}{N - a} \sum_{i=1}^a\sum_{j=1}^{n_i} \hat \varepsilon_{ij}^2 = \frac{1}{N - a} \sum_{i=1}^a\sum_{j=1}^{n_i} (Y_{ij} - \bar {Y_{i.}})^2. \] Divide by $N-a$ since the $N$ residuals depend on $a$ estimated quantities

Rust inhibitors example (cont)

tab <- cbind(rust$brand,rust$score,lm_out$fitted.values,lm_out$residuals)
colnames(tab) <- c("brand","score","Fitted value","Residual")
head(tab,n = 13)

   brand score Fitted value Residual
1      1  43.9        43.14     0.76
2      1  39.0        43.14    -4.14
3      1  46.7        43.14     3.56
4      1  43.8        43.14     0.66
5      1  44.2        43.14     1.06
6      1  47.7        43.14     4.56
7      1  43.6        43.14     0.46
8      1  38.9        43.14    -4.24
9      1  43.6        43.14     0.46
10     1  40.0        43.14    -3.14
11     2  89.8        89.44     0.36
12     2  87.1        89.44    -2.34
13     2  92.7        89.44     3.26

sgsqhat <- sum(lm_out$residuals^2) / (nrow(rust) - 4)
sgsqhat

[1] 6.139833

The value of $\hat \sigma$ is printed in the summary() output:

summary(lm_out)


Call:
lm(formula = score ~ as.factor(brand), data = rust)

Residuals:
   Min     1Q Median     3Q    Max 
-4.270 -1.597  0.395  1.275  4.730 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        43.1400     0.7836  55.056   <2e-16 ***
as.factor(brand)2  46.3000     1.1081  41.782   <2e-16 ***
as.factor(brand)3  24.8100     1.1081  22.389   <2e-16 ***
as.factor(brand)4  -2.6700     1.1081  -2.409   0.0212 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.478 on 36 degrees of freedom
Multiple R-squared:  0.9863,    Adjusted R-squared:  0.9852 
F-statistic: 866.1 on 3 and 36 DF,  p-value: < 2.2e-16

Sums of squares in the one-way ANOVA model

As in linear regression we decompose the variation in the $Y_{ij}$ by defining:

Total sum of squares: $\operatorname{SS}_{\operatorname{Tot}} = \sum_{i=1}^a\sum_{j=1}^{n_i}(Y_{ij} - \bar Y_{..})^2$
Treatment sum of squares: $\operatorname{SS}_{\operatorname{Trt}} = \sum_{i=1}^a n_i(\bar Y_{i.} - \bar Y_{..})^2$
Error sum of squares: $\operatorname{SS}_{\operatorname{Error}} = \sum_{i=1}^a\sum_{j=1}^{n_i}(Y_{ij} - \bar Y_{i.})^2$

In the above, $\bar Y_{..}$ denotes the overall mean, defined as \[ \textstyle \bar Y_{..} = N^{-1}\sum_{i=1}^a\sum_{j=1}^{n_i} Y_{ij}, \quad\text{ where } N = n_1+\dots+n_a. \]

We have $\operatorname{SS}_{\operatorname{Tot}} = \operatorname{SS}_{\operatorname{Trt}} + \operatorname{SS}_{\operatorname{Error}}$.

Note that $\operatorname{SS}_{\operatorname{Trt}}$ is computed just like $\operatorname{SS}_{\operatorname{Reg}}$ in linear regression.

We again define $\displaystyle R^2 = \frac{\operatorname{SS}_{\operatorname{Trt}}}{\operatorname{SS}_{\operatorname{Tot}}}$.

Sampling distributions of our sums of squares

The SS, appropriately scaled, follow chi-square distributions:

$\displaystyle \operatorname{SS}_{\operatorname{Tot}}/\sigma^2 \sim \chi^2_{N-1}(\phi_{\operatorname{Tot}})$
$\displaystyle \operatorname{SS}_{\operatorname{Trt}}/\sigma^2 \sim \chi^2_{a-1}(\phi_{\operatorname{Trt}})$
$\displaystyle \operatorname{SS}_{\operatorname{Error}}/\sigma^2 \sim \chi^2_{N-a}$,

where $\phi_{\operatorname{Tot}}$ and $\phi_{\operatorname{Trt}}$ are noncentrality parameters.

The mean squares in the one-way ANOVA model

Dividing $\operatorname{SS}_{\operatorname{Trt}}$ and $\operatorname{SS}_{\operatorname{Error}}$ by their dfs, we define:

Treatment mean square: $\displaystyle \operatorname{MS}_{\operatorname{Trt}}= \frac{\operatorname{SS}_{\operatorname{Trt}}}{a-1}$
Error mean square: $\displaystyle \operatorname{MS}_{\operatorname{Error}} = \frac{\operatorname{SS}_{\operatorname{Error}}}{N-a}$

The ratio $\displaystyle F_{\operatorname{test}}= \frac{\operatorname{MS}_{\operatorname{Trt}}}{\operatorname{MS}_{\operatorname{Error}}}$ has an F distribution.

The Analysis of Variance (ANOVA) table

We often present the SS, df, and MS values in a table like this:

Source	Df	SS	MS	F value	p-value
Treatment	$a-1$	$\operatorname{SS}_{\operatorname{Trt}}$	$\operatorname{MS}_{\operatorname{Trt}}$	$F_{\operatorname{test}}$	$P(F > F_{\operatorname{test}})$
Error	$N-a$	$\operatorname{SS}_{\operatorname{Error}}$	$\operatorname{MS}_{\operatorname{Error}}$
Total	$N-1$	$\operatorname{SS}_{\operatorname{Tot}}$

In the table $\displaystyle F_{\operatorname{test}}= \frac{\operatorname{MS}_{\operatorname{Trt}}}{\operatorname{MS}_{\operatorname{Error}}}$.

The p-value is based on $F \sim F_{a-1,N-a}$.

Rust inhibitors example (cont)

Obtain the ANOVA table “manually”:

Y <- rust$score
Y.. <- mean(Y)
Yi. <- aggregate(score ~ brand, mean, data = rust)[,2]
n <- aggregate(score ~ brand, length, data = rust)[,2]
a <- length(Yi.)
N <- sum(n)

SStot <- sum((Y - Y..)**2)
SStrt <- sum(n*(Yi. - Y..)**2)
SSerror <- SStot - SStrt

MStrt <- SStrt / (a - 1)
MSerror <- SSerror / (N - a)

Ftest <- MStrt / MSerror
pval <- 1 - pf(Ftest, a-1, N - a)

Rust inhibitors example (cont)

 brand    1    2    3    4    5    6    7    8    9   10  mean    sd  n
     1 43.9 39.0 46.7 43.8 44.2 47.7 43.6 38.9 43.6 40.0 43.14 3.000 10
     2 89.8 87.1 92.7 90.6 87.7 92.4 86.1 88.1 90.8 89.1 89.44 2.218 10
     3 68.4 69.3 68.5 66.4 70.0 68.1 70.6 65.2 63.8 69.2 67.95 2.169 10
     4 36.2 45.2 40.7 40.5 39.3 40.3 43.2 38.7 40.9 39.7 40.47 2.436 10

Rust inhibitors example (cont)

Obtain the ANOVA table with the anova() function on the lm() output.

anova(lm_out)

Analysis of Variance Table

Response: score
                 Df Sum Sq Mean Sq F value    Pr(>F)    
as.factor(brand)  3  15954  5317.8  866.12 < 2.2e-16 ***
Residuals        36    221     6.1                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing whether there is any difference in treatment means

In the one-way treatment effects model we wish to test \[ \text{$H_0$: $\tau_1 = \dots = \tau_a$ \quad versus \quad $H_1$: $\tau_i \neq \tau_{i'}$ for some $i \neq i'$.} \]

We use the overall F test of significance:

Compute $\displaystyle F_{\operatorname{test}}= \frac{\operatorname{MS}_{\operatorname{Trt}}}{\operatorname{MS}_{\operatorname{Error}}}$
Reject $H_0$ at $\alpha$ if $F_{\operatorname{test}}> F_{a-1,N-a,\alpha}$.
Obtain p-value as $P(F > F_{\operatorname{test}})$, where $F \sim F_{a-1,N-a}$.

The value of $\displaystyle F_{\operatorname{test}}$ and the p-value are printed in the summary() output.

Interpretation of F statistic

Note that $F_{\operatorname{test}}$ is a ratio of the form $\displaystyle \frac{\text{Between treatment variation}}{\text{Within treatment variation}}$.

Exercise: For which data set will the F-statistic be largest/smallest?

Exercise: Compute $F_{\operatorname{test}}$ for the rust data using the summary info:

brand	replicates	mean	standard deviation
1	10	43.14	3.00
2	10	89.44	2.22
3	10	67.95	2.17
4	10	40.47	2.44

Hint: $\displaystyle \operatorname{SS}_{\operatorname{Error}}= \sum_{i=1}^a(n_i - 1)S_i^2$, where $\displaystyle S_i^2 = \frac{1}{n_i - 1}\sum_{j=1}^{n_i}(Y_{ij} - \bar Y_{i.})^2$

Some CI formulas (without familywise adjustment)

In the cell-means formulation of the model \[ Y_{ij} = \mu_i + \varepsilon_{ij}, \quad j = 1,\dots,n_i, \quad i = 1,\dots,a, \] where $\mu_i = \mu + \tau_i$, we have the following CI formulas:

Target	$(1-\alpha)100\%$ confidence interval
$\mu_i$	$\bar Y_{i.} \pm t_{N-a,\alpha/2}\hat \sigma \sqrt{\frac{1}{n_i}}$
$\mu_i - \mu_{i'}$	$\bar Y_{i.} - \bar Y_{i'.} \pm t_{N-a,\alpha/2}\hat \sigma \sqrt{\frac{1}{n_i} + \frac{1}{n_{i'}}}$

Rust inhibitors example (cont)

Compute 95% CIs for $\mu_1$ and $\mu_2 - \mu_1$.

alpha <- 0.05
lo1 <- y1bar - qt(1-alpha/2,N-a) * sqrt(sgsqhat) / sqrt(n1)
up1 <- y1bar + qt(1-alpha/2,N-a) * sqrt(sgsqhat) / sqrt(n1)
c(lo1,up1)

[1] 41.55084 44.72916

lo21 <- y2bar - y1bar - qt(1-alpha/2,N-a) * sqrt(sgsqhat) * sqrt(1/n1 + 1/n2)
up21 <- y2bar - y1bar + qt(1-alpha/2,N-a) * sqrt(sgsqhat) * sqrt(1/n1 + 1/n2)
c(lo21,up21)

[1] 44.05259 48.54741

Post-hoc comparisons of means

If we reject $H_0$: $\mu_1 = \dots = \mu_a$, then we may wish to compare means.
Call such comparisons post-hoc as we do them after the F-test.
We may wish to compare several pairs of means, which is like testing several hypotheses at once.
When several hypotheses are tested at once, the familywise Type I error rate is the probability that any Type I error is committed.
We discuss two methods for post-hoc comparisons of means which control the familywise Type I error rate.

Comparing all pairs of means

We want to build a CI for $\mu_i - \mu_{i'}$ for all pairs $i \neq i'$.
Suppose the design is balanced, i.e. $n_i = n$ for all $i = 1,\dots,a$.
If we build for all $i \neq i'$ the ordinary $(1-\alpha)\times 100\%$ CIs \[ \bar Y_{i.} - \bar Y_{i'.} \pm t_{a(n-1),\alpha/2} \hat \sigma \sqrt{2/n}, \] each one will cover its target with probability $1-\alpha$.
But now we want simultaneous coverage with probability $1-\alpha$, i.e. \[ P(\cap_{i \neq i'} \{\text{CI for $\mu_i - \mu_{i'}$ captures target}\}) = 1-\alpha. \]
Above probability is called the familywise coverage.

The venerable John Tukey

Multiple comparisons of means with Tukey’s HSD

Suppose the design is balanced, i.e. $n_i = n$ for all $i = 1,\dots,a$.
Suppose we could find the value $q_{a,a(n-1),\alpha}$ such that \[ P\left(\max_{i \neq i'} \left\{ \frac{|(\bar Y_{i.} - \bar Y_{i'.}) - (\mu_i - \mu_{i'}) |}{\hat \sigma /\sqrt{n}}\right\} \leq q_{a,a(n-1),\alpha} \right) = 1-\alpha. \]
Then with probability $1-\alpha$ the CIs \[ \bar Y_{i.} - \bar Y_{i'.} \pm q_{a,a(n-1),\alpha} \hat \sigma / \sqrt{n} \] will simultaneously cover the targets $\mu_i - \mu_{i'}$ for all $i \neq i'$. Show!
Tukey made tables of the values $q_{a,a(n-1),\alpha}$.
Can use the simultaneous intervals to sort/compare the means.

Table A.6 from Mohr, Wilson, and Freund (2021)

Rust inhibitors example (cont)

For the rust data we have $n = 10$ and $a = 4$.

At $\alpha = 0.05$ we have $q_{a,a(n-1),\alpha} = q_{4,36,0.05} \approx 3.85$ from table.

Obtain exact value with qtukey(.95,4,36) = 3.8087984.

Build the Tukey HSD CI for $\mu_2 - \mu_1$.

n <- 10
a <- 4
MSE <- sum(lm_out$residuals^2) / ( a*(n-1))
y1bar <- mean(rust$score[rust$brand == 1])
y2bar <- mean(rust$score[rust$brand == 2])
me <- qtukey(.95,a,a*(n-1)) * sqrt(MSE) / sqrt(10)
lo21 <- y2bar - y1bar - me
up21 <- y2bar - y1bar + me
c(lo21,up21)

[1] 43.31554 49.28446

Rust inhibitors example (cont)

Use TukeyHSD() on aov() output to obtain the simultaneous CIs.

# must use the aov() function instead of the lm() function
aov_out <- aov(score ~ as.factor(brand), data = rust)
TukeyHSD(aov_out)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = score ~ as.factor(brand), data = rust)

$`as.factor(brand)`
      diff        lwr         upr     p adj
2-1  46.30  43.315536  49.2844635 0.0000000
3-1  24.81  21.825536  27.7944635 0.0000000
4-1  -2.67  -5.654464   0.3144635 0.0933303
3-2 -21.49 -24.474464 -18.5055365 0.0000000
4-2 -48.97 -51.954464 -45.9855365 0.0000000
4-3 -27.48 -30.464464 -24.4955365 0.0000000

plot(TukeyHSD(aov_out))

Comparison of treatments with a baseline treatment

It may be that not all pairwise comparisons are of interest.
Then Tukey’s method is too conservative (CIs wider than necessary).
Say we want to compare all treatments to a “baseline” treatment.
Build CIs for $\mu_i - \mu_1$, $i = 2,\dots,a$, $1$ the baseline treatment.
This makes $a-1$ CIs instead of ${a\choose 2}$ CIs.
Can use Dunnett’s method, Dunnett (1964).

The equally venerable Charles Dunnett

Charles Dunnett, 1921 – 2007 (Canadian, served in WWII, photo taken in Belgium)

Dunnett’s method for comparisons with a baseline

Assume $n_i = n$ for all $i$ (balanced case).
Given a value $d_{n,a(n-1),\alpha}$ such that \[ P\left( \max_{2 \leq i \leq a} \Big| \frac{(\bar Y_{i.} - \bar Y_{1.}) - (\mu_i - \mu_1)}{\hat \sigma\sqrt{2/n}} \Big| \leq d_{n,a(n-1),\alpha}\right) = 1 - \alpha, \] with probability $1-\alpha$ the CIs \[ \bar Y_{i.} - \bar Y_{1.} \pm d_{n,a(n-1),\alpha} \hat \sigma\sqrt{2/n} \] will simultaneously cover the targets $\mu_i - \mu_1$ for all $i=2,\dots,a$.
Dunnett made tables of the values $d_{n,a(n-1),\alpha}$.
Cannot sort all the means after Dunnett’s.

Table A.5 from Mohr, Wilson, and Freund (2021)

Rust inhibitor data (cont)

For the rust data we have $n = 10$ and $a = 4$.

At $\alpha = 0.05$ we have $d_{a,a(n-1),\alpha} = d_{4,36,0.05}$.

Use value $2.44$ in the table (should be close).

Treat Brand 1 as the baseline and make comparisons with Dunnett’s.

# just show the comparison of treatment 2 to the baseline
y1bar <- mean(rust$score[rust$brand == 1])
y2bar <- mean(rust$score[rust$brand == 2])

me <- 2.44 * sqrt(MSE) * sqrt(2/10) # margin of error for Dunnett's

lo21 <- y2bar - y1bar - me
up21 <- y2bar - y1bar + me

c(y2bar - y1bar,lo21,up21)

[1] 46.30000 43.59615 49.00385

Rust inhibitor data (cont)

Use DunnettTest() from R package DescTools.

library(DescTools) # first time run install.packages("DescTools")

Warning: package 'DescTools' was built under R version 4.4.1

Dunnett_out <- DunnettTest(score ~ as.factor(brand), data = rust, control = "1")
Dunnett_out


  Dunnett's test for comparing several treatments with a control :  
    95% family-wise confidence level

$`1`
     diff    lwr.ci    upr.ci   pval    
2-1 46.30 43.582516 49.017484 <2e-16 ***
3-1 24.81 22.092516 27.527484 <2e-16 ***
4-1 -2.67 -5.387484  0.047484 0.0549 .  

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(Dunnett_out)

Dunnett’s vs Tukey’s

Tukey’s is for comparisons between all pairs of means.
Dunnett’s is for comparison of means with a baseline.
So Tukey’s must make greater adjustments to control the familywise Type I error.
Therefore Tukey intervals will be wider than Dunnett intervals.
Tukey’s allows you to sort the means, while Dunnett’s does not.
Both methods assume a balanced design, i.e. $n_i = n$ for all $i$. Modifications for unbalanced designs exist, but are not straightforward to implement in R.

Bonferroni correction

If building $B$ CIs you can ALWAYS use the Bonferroni correction:

Build each CI ordinarily, but use $\alpha/B$ instead of $\alpha$.
Ensures simultaneous coverage of all CIs with probability $\geq 1-\alpha$.
True prob of simultaneous coverage may be greater than $1-\alpha$
Bonferroni-corrected CIs will be wider than Dunnett’s and wider than Tukey’s if used for making those same comparisons.
Use when we do not know how to adjust for multiple comparisons.

Rust inhibitor data (cont)

Compare Brand 3 to 4 and Brand 1 to 3, using the Bonferroni correction to control the familywise error rate.

y1bar <- mean(rust$score[rust$brand == 1])
y3bar <- mean(rust$score[rust$brand == 3])
y4bar <- mean(rust$score[rust$brand == 4])
alpha <- 0.05
B <- 2
me <- qt(1 - (alpha/B)/2,a*(n-1)) * sqrt(MSE) * sqrt(2/n)
tab <- rbind(c(y3bar - y4bar - me,y3bar - y4bar + me),
             c(y1bar - y3bar - me,y1bar - y3bar + me))
rownames(tab) <- c("3-4","1-3")
colnames(tab) <- c("lower","upper")
tab

      lower   upper
3-4  24.888  30.072
1-3 -27.402 -22.218

Checking model assumptions

Validity of the foregoing analyses depends on these assumptions:

The responses are normally distributed around the treatment means (Check QQ plot of residuals).
The response has the same variance in all treatment groups (Check residuals vs fitted values plot).
The response values are independent of each other (No way to check; must trust experimental design).

Rust inhibitors example (cont)

plot(lm_out,which = 2)

Rust inhibitors example (cont)

plot(lm_out,which = 1, add.smooth = F) # we don't want the red line

Perception of slope example

Do axis re-scalings affect how we perceive an x-y relationship?

For a single data set with data pairs $(X_i,Y_i)$, with $X_i \sim \text{Normal}(0,1)$ and $Y_i = \text{Normal}(X_i,1)$ for $i = 1,\dots,50$, three scatterplot treatments were constructed:

“Control” used x and y plotting limits given by the range of the data.
“X” extended the x-limits by 1.5 in each direction.
“Y” extended the y-limits by 1.5 in each direction.

Each student in a class was randomly assigned a scatterplot and told to draw with a ruler the best-fitting line through the data. The slope of each student-drawn line was measured and recorded as the response.

Is the response mean the same in the three treatment groups?

An artifact from each treatment group:

slope <- c(1.23,1.80,1.81,1.29,2.89,1.58,0.99,1.24,
           1.26,1.57,1.27,1.19,1.82,1.76,1.91,1.25,
           1.09,1.29,1.12,1.51,2.13,1.16,0.62,1.04)
trt <- c("X","Y","X","X","Y","X","Y","C",
         "Y","C","C","C","Y","C","X","Y",
         "X","X","Y","C","Y","X","Y","C")

boxplot(slope ~ trt)

lm_slope <- lm(slope ~ as.factor(trt))
summary(lm_slope)


Call:
lm(formula = slope ~ as.factor(trt))

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9222 -0.2847 -0.1293  0.2628  1.3478 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1.36857    0.18161   7.536 2.12e-07 ***
as.factor(trt)X  0.05143    0.24868   0.207    0.838    
as.factor(trt)Y  0.17365    0.24215   0.717    0.481    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4805 on 21 degrees of freedom
Multiple R-squared:  0.02614,   Adjusted R-squared:  -0.06661 
F-statistic: 0.2818 on 2 and 21 DF,  p-value: 0.7572

plot(lm_slope,which = 2)

plot(lm_slope,which = 1)

Levene’s test for equality of variances

Checks if the mean magnitude of the residuals is equal across groups:

Obtain the residuals $\hat \varepsilon_{ij}$ from the one-way ANOVA model.
Treat the absolute values $|\hat \varepsilon_{ij}|$ of the residuals as new responses.
Test for equal means of the new responses with the F test.

So, do the ordinary F-test with the $|\hat \varepsilon_{ij}|$ as the responses.

Perception of slope example (cont)

Perform Levene’s test:

ehat <- lm_slope$residuals
lm_levene <- lm(abs(ehat) ~ as.factor(trt))
summary(lm_levene)


Call:
lm(formula = abs(ehat) ~ as.factor(trt))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.29136 -0.12769 -0.04980  0.08219  0.79864 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)      0.20980    0.09352   2.243   0.0358 *
as.factor(trt)X  0.05020    0.12805   0.392   0.6990  
as.factor(trt)Y  0.33934    0.12469   2.721   0.0128 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2474 on 21 degrees of freedom
Multiple R-squared:  0.303, Adjusted R-squared:  0.2367 
F-statistic: 4.565 on 2 and 21 DF,  p-value: 0.02258

Can also use the leveneTest() function in the R package car.

library(car)

Warning: package 'car' was built under R version 4.4.3

Loading required package: carData


Attaching package: 'car'

The following object is masked from 'package:DescTools':

    Recode

leveneTest(slope~as.factor(trt),center = mean)

Levene's Test for Homogeneity of Variance (center = mean)
      Df F value  Pr(>F)  
group  2  4.5652 0.02258 *
      21                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We conclude that the variances are not equal across treatment groups.

References

Dunnett, Charles W. 1964. “New Tables for Multiple Comparisons with a Control.” Biometrics 20 (3): 482–91.

Kutner, Michael H, Christopher J Nachtsheim, John Neter, and William Li. 2005. Applied Linear Statistical Models. McGraw-hill.

Mohr, Donna L, William J Wilson, and Rudolf J Freund. 2021. Statistical Methods. Academic Press.

Source	Df	SS	MS	F value	p-value
Treatment	\(a-1\)	\(\operatorname{SS}_{\operatorname{Trt}}\)	\(\operatorname{MS}_{\operatorname{Trt}}\)	\(F_{\operatorname{test}}\)	\(P(F > F_{\operatorname{test}})\)
Error	\(N-a\)	\(\operatorname{SS}_{\operatorname{Error}}\)	\(\operatorname{MS}_{\operatorname{Error}}\)
Total	\(N-1\)	\(\operatorname{SS}_{\operatorname{Tot}}\)

Target	\((1-\alpha)100\%\) confidence interval
\(\mu_i\)	\(\bar Y_{i.} \pm t_{N-a,\alpha/2}\hat \sigma \sqrt{\frac{1}{n_i}}\)
\(\mu_i - \mu_{i'}\)	\(\bar Y_{i.} - \bar Y_{i'.} \pm t_{N-a,\alpha/2}\hat \sigma \sqrt{\frac{1}{n_i} + \frac{1}{n_{i'}}}\)

Rust inhibitors example

Randomized experiments comparing treatments

One-way treatment effects model

Identifiability constraint in the treatment effects model

Alternative “cell means” setup

Goals in one-way ANOVA

Treatment effect estimation in one-way ANOVA

Rust inhibitors example (cont)

Estimation of the error term variance \(\sigma^2\)

Rust inhibitors example (cont)

Sums of squares in the one-way ANOVA model

Sampling distributions of our sums of squares

The mean squares in the one-way ANOVA model

The Analysis of Variance (ANOVA) table

Rust inhibitors example (cont)

Rust inhibitors example (cont)

Rust inhibitors example (cont)

Testing whether there is any difference in treatment means

Interpretation of F statistic

Some CI formulas (without familywise adjustment)

Rust inhibitors example (cont)

Post-hoc comparisons of means

Comparing all pairs of means

The venerable John Tukey

Multiple comparisons of means with Tukey’s HSD

Rust inhibitors example (cont)

Rust inhibitors example (cont)

Comparison of treatments with a baseline treatment

The equally venerable Charles Dunnett

Dunnett’s method for comparisons with a baseline

Rust inhibitor data (cont)

Rust inhibitor data (cont)

Dunnett’s vs Tukey’s

Bonferroni correction

Rust inhibitor data (cont)

Checking model assumptions

Rust inhibitors example (cont)

Rust inhibitors example (cont)

Perception of slope example

Levene’s test for equality of variances

Perception of slope example (cont)

References