ANOVA Examples in R

• ANOVA Examples
  • One way ANOVA
    • Effect Size $\omega^2$
  • Two-Way ANOVA

ANOVA Examples

Examples adapted from the following work:

Pagano, R. R. (2002). Understanding statistics in the behavioral sciences (6th ed.). Belmont, CA: Wadsworth.

One way ANOVA

From pages 357 - 360.

The null hypothesis:

$H_0: mu_1 = mu_2 = mu_3$
$H_1:$ At least two means differ.

The data in grouped, tabular format:

Create the data:

x  <- c(2,3,7,2,6,10,8,7,5,10,10,13,14,13,15)
g  <- c(rep("group1", 5), rep("group2", 5), rep("group3", 5))
xg <- data.frame(x,g)

Examine the data:

    x      g
1   2 group1
2   3 group1
3   7 group1
4   2 group1
5   6 group1
6  10 group2
7   8 group2
8   7 group2
9   5 group2
10 10 group2
11 10 group3
12 13 group3
13 14 group3
14 13 group3
15 15 group3

Run the test:

fit.1 <- aov(x ~ g, data = xg)
summary(fit.1)
            Df Sum Sq Mean Sq F value   Pr(>F)    
g            2  203.3   101.7   22.59 8.54e-05 ***
Residuals   12   54.0     4.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Effect Size $\omega^2$

Calculate size of effect (p. 367):

$\omega^2 = \frac{SS_B - (k - 1)s_w^2}{SS_T + s_w^2}$

Where:

• $SS_B$: Sum of Squares Between (Sum Sq of g in R summary output)
• $k$: Number of groups
• $s_w^2$: Within-groups variance (Mean Sq of Residuals in R summary output)
• $SS_T$: Sum of Squares Total (Not shown in R summary output; add Sum Sq values)

Found this function to compute $\omega^2$ on StackOverflow:

https://stats.stackexchange.com/questions/2962/omega-squared-for-measure-of-effect-in-r

omega_sq <- function(aovm){
  sum_stats <- summary(aovm)[[1]]
  SSm <- sum_stats[["Sum Sq"]][1]
  SSr <- sum_stats[["Sum Sq"]][2]
  DFm <- sum_stats[["Df"]][1]
  MSr <- sum_stats[["Mean Sq"]][2]
  W2 <- (SSm-DFm*MSr)/(SSm+SSr+MSr)
  return(W2)
}

Then run the function on the model:

omega_sq(fit.1)
[1] 0.7422024

Two-Way ANOVA

There are three null hypotheses in a two-way ANOVA:

• $H_0$ Main Effect A (Rows): $\mu_{a1} = \mu_{a2}$
• $H_0$ Main Effect B (Columns): $\mu_{b1} = \mu_{b2} = \mu_{b3}$
• $H_0$ Interaction Effect A:B: $\mu_{a_1b_1} = \mu_{a_1b_3} = \mu_{a_2b_1} = \mu_{a_2b_2} = \mu_{a_2b_3}$

The data in grouped, tabular format:

|---------------|--------------------------------|
|               |  Exercise                      |
|---------------|--------------------------------|
|  Time of Day  |  Light  |  Moderate  |  Heavy  |
|---------------|---------|------------|---------|
|  Morning      |  6.5    |  7.4       |  8.0    |
|               |  7.3    |  6.8       |  7.7    |
|               |  6.6    |  6.7       |  7.1    |
|               |  7.4    |  7.3       |  7.6    |
|               |  7.2    |  7.6       |  6.6    |
|               |  6.8    |  7.4       |  7.2    |
|---------------|---------|------------|---------|
|  Evening      |  7.1    |  7.4       |  8.2    |
|               |  7.9    |  8.1       |  8.5    |
|               |  8.2    |  8.2       |  9.5    |
|               |  7.7    |  8.0       |  8.7    |
|               |  7.5    |  7.6       |  9.6    |
|               |  7.6    |  8.0       |  9.4    |
|---------------|---------|------------|---------|

First, load the data:

hours <- c(6.5, 7.3, 6.6, 7.4, 7.2, 6.8, 7.4, 6.8, 6.7,
           7.3, 7.6, 7.4, 8.0, 7.7, 7.1, 7.6, 6.6, 7.2,
           7.1, 7.9, 8.2, 7.7, 7.5, 7.6, 7.4, 8.1, 8.2,
           8.0, 7.6, 8.0, 8.2, 8.5, 9.5, 8.7, 9.6, 9.4)
tod   <- c(rep("morning", 18), rep("evening", 18))
exer  <- rep(c(rep("light", 6), rep("moderate", 6), rep("heavy", 6)), 2)
an.ex <- data.frame(hours, tod, exer) ; rm(hours, tod, exer)

These are ordered factors:

an.ex$tod  <- ordered(an.ex$tod, levels = c("morning", "evening"))
an.ex$exer <- ordered(an.ex$exer, levels = c("light", "moderate", "heavy"))

Examine the data frame (truncated here):

hours     tod     exer
1    6.5 morning    light
2    7.3 morning    light
...
7    7.4 morning moderate
8    6.8 morning moderate
...
35   9.6 evening    heavy
36   9.4 evening    heavy

Run the model:

fit.1 <- aov(hours ~ tod + exer + tod:exer, data = an.ex)
summary(fit.1)
            Df Sum Sq Mean Sq F value   Pr(>F)    
tod          1  9.000   9.000  48.416 9.93e-08 ***
exer         2  4.754   2.377  12.787 9.63e-05 ***
tod:exer     2  1.712   0.856   4.604    0.018 *  
Residuals   30  5.577   0.186                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Alternate method:

library(car) # alternate way
fit.2 <- lm(hours ~ tod + exer + tod:exer, data = an.ex)
Anova(fit.2, type = "II")