R Companion: Two-way Anova

The following commands will install these packages if they are not already installed:

if(!require(FSA)){install.packages("FSA")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(car)){install.packages("car")}
if(!require(multcomp)){install.packages("multcomp")}
if(!require(emmeans)){install.packages("emmeans")}
if(!require(nlme)){install.packages("nlme")}
if(!require(lme4)){install.packages("lme4")}
if(!require(lmerTest)){install.packages("lmerTest")}
if(!require(rcompanion)){install.packages("rcompanion")}

When to use it

Null hypotheses

How the test works

Assumptions

See the Handbook for information on these topics.

Examples

The rattlesnake example is shown at the end of the “How to do the test” section.

How to do the test

For notes on linear models and conducting anova, see the “How to do the test” section in the One-way anova chapter of this book. For two-way anova with robust regression, see the chapter on Two-way Anova with Robust Estimation.

Two-way anova example

### --------------------------------------------------------------
### Two-way anova, SAS example, pp. 179–180
### --------------------------------------------------------------

Input = ("
id Sex     Genotype Activity
1 male    ff        1.884
2 male    ff        2.283
3 male    fs        2.396
4 female ff        2.838
5 male    fs        2.956
6 female ff        4.216
7 female ss        3.620
8 female ff        2.889
9 female fs        3.550
10 male    fs        3.105
11 female fs        4.556
12 female fs        3.087
13 male    ff        4.939
14 male    ff        3.486
15 female ss        3.079
16 male    fs        2.649
17 female fs        1.943
19 female ff        4.198
20 female ff        2.473
22 female ff        2.033
24 female fs        2.200
25 female fs        2.157
26 male    ss        2.801
28 male    ss        3.421
29 female ff        1.811
30 female fs        4.281
32 female fs        4.772
34 female ss        3.586
36 female ff        3.944
38 female ss        2.669
39 female ss        3.050
41 male    ss        4.275
43 female ss        2.963
46 female ss        3.236
48 female ss        3.673
49 male    ss        3.110
")

Data = read.table(textConnection(Input),header=TRUE)

Means and summary statistics by group

library(FSA)

Sum = Summarize(Activity ~ Sex + Genotype,
data = Data)

Sum

     Sex Genotype n    mean        sd   min    Q1 median    Q3   max
1 female       ff 8 3.05025 0.9599032 1.811 2.363 2.864 4.008 4.216
2   male       ff 4 3.14800 1.3745115 1.884 2.183 2.884 3.849 4.939
3 female       fs 8 3.31825 1.1445388 1.943 2.189 3.318 4.350 4.772
4   male       fs 4 2.77650 0.3168433 2.396 2.586 2.802 2.993 3.105
5 female       ss 8 3.23450 0.3617754 2.669 3.028 3.158 3.594 3.673
6   male       ss 4 3.40175 0.6348109 2.801 3.033 3.266 3.634 4.275

### Add standard error

Sum$se = Sum$sd/sqrt(Sum$n)

Sum

     Sex Genotype n    mean        sd   min    Q1 median    Q3   max        se
1 female       ff 8 3.05025 0.9599032 1.811 2.363 2.864 4.008 4.216 0.3393770
2   male       ff 4 3.14800 1.3745115 1.884 2.183 2.884 3.849 4.939 0.6872558
3 female       fs 8 3.31825 1.1445388 1.943 2.189 3.318 4.350 4.772 0.4046556
4   male       fs 4 2.77650 0.3168433 2.396 2.586 2.802 2.993 3.105 0.1584216
5 female       ss 8 3.23450 0.3617754 2.669 3.028 3.158 3.594 3.673 0.1279069
6   male       ss 4 3.40175 0.6348109 2.801 3.033 3.266 3.634 4.275 0.3174054

Interaction plot using summary statistics

library(ggplot2)

pd = position_dodge(.2)

ggplot(Sum, aes(x=Genotype,
                y=mean,
                color=Sex)) +
    geom_errorbar(aes(ymin=mean-se,
                      ymax=mean+se),
                   width=.2, linewidth=0.7, position=pd) +
    geom_point(shape=15, size=4, position=pd) +
    theme_bw() +
    theme(axis.title.y = element_text(vjust= 1.8),
          axis.title.x = element_text(vjust= -0.5),
          axis.title = element_text(face = "bold")) +

scale_color_manual(values = c("black", "blue"))+

ylab("Activity")

Interaction plot for a two-way anova. Square points represent means for groups, and error bars indicate standard errors of the mean.

Simple box plot of main effect and interaction

boxplot(Activity ~ Genotype,
        data = Data,
        xlab = "Genotype",
        ylab = "MPI Activity",
        col = "white")

Rplot

boxplot(Activity ~ Genotype:Sex,
        data = Data,
        xlab = "Genotype x Sex",
        ylab = "MPI Activity",
        col = "white")

Rplot03

Fit the linear model and conduct ANOVA

model = lm(Activity ~ Sex + Genotype + Sex:Genotype,
           data=Data)

library(car)

Anova(model,
      type="II")   ### Type II sum of squares

   ### If you use type="III", you need the following line before the analysis
   ### options(contrasts = c("contr.sum", "contr.poly"))

              Sum Sq Df F value Pr(>F)
Sex           0.0681 1 0.0861 0.7712
Genotype      0.2772 2 0.1754 0.8400
Sex:Genotype 0.8146 2 0.5153 0.6025
Residuals    23.7138 30

anova(model) # Produces type I sum of squares

             Df Sum Sq Mean Sq F value Pr(>F)
Sex           1 0.0681 0.06808 0.0861 0.7712
Genotype      2 0.2772 0.13862 0.1754 0.8400
Sex:Genotype 2 0.8146 0.40732 0.5153 0.6025
Residuals    30 23.7138 0.79046

summary(model) # Produces r-square, overall p-value, parameter estimates

Multiple R-squared: 0.04663, Adjusted R-squared: -0.1123

F-statistic: 0.2935 on 5 and 30 DF, p-value: 0.9128

Checking assumptions of the model

hist(residuals(model),
col="darkgray")

RPlot

A histogram of residuals from a linear model. The distribution of these residuals should be approximately normal.

plot(fitted(model),
residuals(model))

RPlot

A plot of residuals vs. predicted values. The residuals should be unbiased and homoscedastic. For an illustration of these properties, see the Introduction to Parametric Tests in SAEPER (rcompanion.org/handbook/I_01.html).

Additional model checking plots

plot(model)

Post-hoc comparison of least-square means

For notes on least-square means, see the “Post-hoc comparison of least-square means” section in the Nested anova chapter in this book.

One advantage of the using the emmeans package for post-hoc tests is that it can produce comparisons for interaction effects.

In general, if the interaction effect is significant, you will want to look at comparisons of means for the interactions. If the interaction effect is not significant but a main effect is, it is appropriate to look at comparisons among the means for that main effect. In this case, because no effect of Sex, Genotype, or Sex:Genotype was significant, we would not actually perform any mean separation test.

Mean separations for main factor with emmeans

library(emmeans)

marginal = emmeans(model, ~ Genotype)

pairs(marginal, adjust="tukey")

contrast estimate    SE df t.ratio p.value
ff - fs    0.0517 0.385 30   0.134 0.9901
ff - ss   -0.2190 0.385 30 -0.569 0.8376
fs - ss   -0.2707 0.385 30 -0.703 0.7634

Results are averaged over the levels of: Sex
P value adjustment: tukey method for comparing a family of 3 estimates

library(multcomp)

cld(marginal, Letters=letters, adjust="tukey")

Genotype emmean    SE df lower.CL upper.CL .group
fs         3.05 0.272 30     2.36     3.74 a
ff         3.10 0.272 30     2.41     3.79 a
ss         3.32 0.272 30     2.63     4.01 a

Results are averaged over the levels of: Sex
Confidence level used: 0.95
Conf-level adjustment: sidak method for 3 estimates
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05

    ### Means sharing a letter in .group are not significantly different

Mean separations for interaction effect with emmeans

library(emmeans)

marginal = emmeans(model, ~ Sex:Genotype)

pairs(marginal, adjust="tukey")

### Results not shown

library(multcomp)

cld(marginal, Letters=letters, adjust="tukey")

Sex Genotype emmean SE df lower.CL upper.CL .group

male fs 2.78 0.445 30 1.52 4.03 a

female ff 3.05 0.314 30 2.17 3.94 a

male ff 3.15 0.445 30 1.90 4.40 a

female ss 3.23 0.314 30 2.35 4.12 a

female fs 3.32 0.314 30 2.43 4.20 a

male ss 3.40 0.445 30 2.15 4.65 a

Confidence level used: 0.95

Conf-level adjustment: sidak method for 6 estimates

P value adjustment: tukey method for comparing a family of 6 estimates

significance level used: alpha = 0.05

### Note that means are listed from low to high,
### not in the same order as Summarize

Graphing the results

Simple bar plot with categories and no error bars

The means for the Sex x Genotype interaction can be extracted from a data frame created with the Summarize function in the FSA package, and organized into a matrix, so that they can be plotted in a simple bar plot.

library(FSA)

Sum = Summarize(Activity ~ Sex + Genotype,
data = Data)

Matriz = xtabs(mean ~ Sex + Genotype, data=Sum)

Matriz

        Genotype
Sex           ff      fs      ss
female 3.05025 3.31825 3.23450
male   3.14800 2.77650 3.40175

row.names(Matriz) = c("Female", "Male")

Matriz

        Genotype
Sex           ff      fs      ss
Female 3.05025 3.31825 3.23450
Male   3.14800 2.77650 3.40175

barplot(Matriz,
        beside=TRUE,
        legend=TRUE,
        ylim=c(0, 5),
        xlab="Genotype",
        ylab="MPI Activity")

RPlot

Bar plot with error bars with ggplot2

Again, a data frame is created with the Summarize function. Error bars indicate standard error of the means (se in the data frame).

library(FSA)

Sum = Summarize(Activity ~ Sex + Genotype,
data = Data)

Sum

### Add standard error

Sum$se = Sum$sd/sqrt(Sum$n)

Sum

### Plot adapted from:
###   shinyapps.stat.ubc.ca/r-graph-catalog/

library(ggplot2)

ggplot(Sum,
   aes(x = Genotype,
       y = mean,
       fill = Sex,
       ymax=mean+se,
       ymin=mean-se)) +
       geom_bar(stat="identity", position = "dodge", width = 0.7) +
       geom_bar(stat="identity", position = "dodge",
                colour = "black", width = 0.7,
                show.legend = FALSE) +
       scale_y_continuous(breaks = seq(0, 4, 0.5),
                limits = c(0, 4),
                expand = c(0, 0)) +
       scale_fill_manual(name = "Count type" ,
                 values = c('grey80', 'grey30'),
                 labels = c("Female",
                            "Male")) +
       geom_errorbar(position=position_dodge(width=0.7),
                     width=0.0, size=0.5, color="black") +
       labs(x = "\nGenotype",
            y = "MPI Activity\n") +
       ## ggtitle("Main title") +
       theme_bw() +
       theme(panel.grid.major.x = element_blank(),
             panel.grid.major.y = element_line(colour = "grey50"),
             plot.title = element_text(size = rel(1.5),
             face = "bold", vjust = 1.5),
             axis.title = element_text(face = "bold"),
             legend.position = "top",
             legend.title = element_blank(),
             legend.key.size = unit(0.4, "cm"),
             legend.key = element_rect(fill = "black"),
             axis.title.y = element_text(vjust= 1.8),
             axis.title.x = element_text(vjust= -0.5))

RPlot

Bar plot for a two-way anova. Bar heights represent means for groups, and error bars indicate standard errors of the mean.

Rattlesnake example – two-way anova without replication, repeated measures

This example could be interpreted as two-way anova without replication or as a one-way repeated measures experiment. Below it is analyzed as a two-way fixed effects model using the lm function, and as a mixed effects model using the nlme package and lme4 package. In addition, a traditional repeated measures anova is conducted using the aov() function in the native stats package.

### --------------------------------------------------------------
### Two-way anova, rattlesnake example, pp. 177–178
### --------------------------------------------------------------

Data = read.table(header=TRUE, stringsAsFactors=TRUE, text="

Day Snake Openings
1    D1        85
1    D3       107
1    D5        61
1    D8        22
1    D11       40
1    D12       65
2    D1        58
2    D3        51
2    D5        60
2    D8        41
2    D11       45
2    D12       27
3    D1        15
3    D3        30
3    D5        68
3    D8        63
3    D11       28
3    D12        3
4    D1        57
4    D3        12
4    D5        36
4    D8        21
4    D11       10
4    D12       16
")

### Treat Day as a factor variable

Data$Day = as.factor(Data$Day)

str(Data)

'data.frame': 24 obs. of 3 variables:

$ Day : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 2 2 2 2 ...

$ Snake : Factor w/ 6 levels "D1","D11","D12",..: 1 4 5 6 2 3 1 4 5 6 ...

$ Openings: int 85 107 61 22 40 65 58 51 60 41 ...

Using two-way fixed effects model