ANOVA tables in R

I don’t know what fears keep you up at night, but for me it’s worrying that I might have copy-pasted the wrong values over from my output. No matter how carefully I check my work, there’s always the nagging suspicion that I could have confused the contrasts for two different factors, or missed a decimal point or a negative sign.

Although I’m usually overreacting, I think my paranoia isn’t completely misplaced — little errors are much too easy to make, and they can have horrifying consequences.

Through the years, I’ve learned that the only sure way to reduce human error is to give humans (including myself) as little opportunity to interfere in the process as possible. Happily, R integrates beautifully with output documents, allowing you to ask the computer to fill in the numbers in your tables and text for you, so you never have to wake up in a cold sweat panicking about a typo in your correlation matrix. These are called dynamic documents, and they’re awesome.

I almost never type out my results anymore; I let R do it for me. I wrote my entire dissertation in R Studio, in fact, using sweave to integrate my R code with LaTeX typesetting. I’m writing this blog post in R Studio as an R-markdown document; if you want to see the raw .rmd file for this post, it’s available on my github: ANOVA_tables.Rmd

Even if you’ve never used markdown or R-markdown before, you can jump right in and start getting APA output from R. In this tutorial, I’ll cover examples for one common model (an analysis of variance, or ANOVA) and show you how you can get APA style output automatically.

We’ll use one of the most basic functions for creating tables, kable, which is from one of the most user-friendly packages for combining R code and output together, knitr. (Side note to the fiber enthusiasts: Yes, you’re not imagining it — pretty much all of this stuff is playfully named after yarn, knitting, and textile references. Enjoy.) I recommend knitr and kable for people just getting into writing dynamic documents because they’re the easiest and most flexible tools, especially since they can be used to create Word documents (not just pdfs or html pages). Depending on your desired output, though, you may find other packages better suited to your needs. For example, if you’re creating pdf documents, you may prefer pander, xtable or stargazer, all of which are much more powerful and elegant. Although these are excellent packages, I find they don’t work consistently (or at all) for Word output, which is a deal breaker for a lot of people.

Quick links to content in this tutorial:

Running an ANOVA in R

Annotated output

Creating an APA style ANOVA table from R output

Inline APA-style output

Recommended references

This tutorial assumes…

Running an ANOVA in R

Set up

Since the kable function is part of the knitr package, you’ll need to load knitr before you can use it:

library(knitr)

We’ll use a data set that comes built into R: warpbreaks. Fittingly, it’s about yarn. It gives the number of breaks in yarn tested under conditions of low, medium, or high tension, for two types of wool (A and B). This data comes standard with R, so you already have it on your computer. You can read the help documentation about this data set by typing ?warpbreaks.

str(warpbreaks) # check out the structure of the data

We’ll run a 2x3 ANOVA to test if there are differences in the number of breaks based on the type of wool and the amount of tension.

Exploratory data analysis

Before running a model, you always want to plot the data, to check that your assumptions look okay. Here are a couple plots I might generate while analyzing these data:

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.3.2

# histograms, to check out the distribution within each group
ggplot(warpbreaks, aes(x=breaks)) + 
  geom_histogram(bins=10) + 
  facet_grid(wool ~ tension) + 
  theme_classic()

# boxplot, to highlight the group means
ggplot(warpbreaks, aes(y=breaks, x=tension, fill = wool)) + 
  geom_boxplot() + 
  theme_classic()

The box plot gives me an idea of what I might find in the ANOVA. It looks like there are differences between groups, with fewer breaks at higher tension, and perhaps fewer breaks in wool B vs. wool A at both low and high tension.

The distributions within each cell look pretty wonky, but that’s not particularly surprising given the small sample size (n=9):

xtabs(~ wool + tension, data = warpbreaks)

##     tension
## wool L M H
##    A 9 9 9
##    B 9 9 9

Running the model

One important consideration when running ANOVAs in R is the coding of factors (in this case, wool and tension). By default, R uses traditional dummy coding (also called “treatment” coding), which works great for regression-style output but can produce weird sums of squares estimates for ANOVA style output.

To be on the safe side, always use effects coding (contr.sum) or orthogonal contrast coding (e.g. contr.helmert, contr.poly) for factors when running an ANOVA. Here, I’m choosing to use effects coding for wool, and polynomial trend contrasts for tension.

model <- lm(breaks ~ wool * tension, 
            data = warpbreaks, 
            contrasts = list(wool = "contr.sum", tension = "contr.poly"))

Annotated ANOVA output

APA style ANOVA tables generally include the sums of squares, degrees of freedom, F statistic, and p value for each effect. You can get all of those calculations with the Anova function from the car package. It’s important to use the Anova function rather than the summary.aov function in base R because Anova allows you to control the type of sums of squares you want to calculate, whereas summary.aov only uses Type 1 (generally not what you want, especially if you have an unblanced design and/or any missing data).

library(car)
sstable <- Anova(model, type = 3) # Type III sums of squares is typical in social science research (it's the default in SPSS)

sstable 

## Anova Table (Type III tests)
## 
## Response: breaks
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept)   42785  1 357.4672 < 2.2e-16 ***
## wool            451  1   3.7653 0.0582130 .  
## tension        2034  2   8.4980 0.0006926 ***
## wool:tension   1003  2   4.1891 0.0210442 *  
## Residuals      5745 48                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The above code runs the Anova function on the model I saved before, using Type III sums of squares, and saves the resulting table as a new object called sstable.

In sstable, you can see a row for each predictor in the model, including the intercept, and the error term (Residuals) at the bottom. The wool:tension term is the interaction between wool and tension (R uses : to specify interaction terms, and * as shorthand for an interaction term with both main effects). There are two levels of wool (A and B), so you’ll see 1 degree of freedom for that effect. There are three levels of tension (low, medium, and high), so that has 2 degrees of freedom. The interaction has the df for both terms multiplied together, i.e. 1 * 2 = 2. The degrees of freedom for the residual are based on the total number of observations in the data (N=54) minus the number of groups, i.e. 54-6=48.

The F-statistic for each effect is the *SS*/*df* for that effect divided by the *SS*/*df* for the residual. The Pr(>F) gives the p value for that test, i.e. the probability of observing an F ratio greater than that given the null hypothesis is true.

Contrast estimates

summary.aov(model, split = list(tension=list(L=1, Q=2)))

##                   Df Sum Sq Mean Sq F value   Pr(>F)    
## wool               1    451   450.7   3.765 0.058213 .  
## tension            2   2034  1017.1   8.498 0.000693 ***
##   tension: L       1   1951  1950.7  16.298 0.000194 ***
##   tension: Q       1     84    83.6   0.698 0.407537    
## wool:tension       2   1003   501.4   4.189 0.021044 *  
##   wool:tension: L  1    251   250.7   2.095 0.154327    
##   wool:tension: Q  1    752   752.1   6.284 0.015626 *  
## Residuals         48   5745   119.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Note that if you run summary(model) instead, you’ll get the default regression-style output, which is the same information, but represented as regression coefficients with standard errors and t-tests instead of sums of squares estimates with F ratios.)

What we see here is a significant linear trend in the main effect for tension — from the box plot we made before, we know that it’s a negative linear trend. There’s no overall quadratic trend in tension, but there is a quadratic trend in the interaction between wool and tension. That means the estimate of the quadratic trend contrast is different for wool A compared to wool B. Referencing the box plot again, you can see that the direction of the quadratic trend appears to differ in wool A compared to wool B (wool A looks like a positive quadratic trend, with an upward swoop, whereas wool B looks like a negative quadratic trend, with an upside-down U swoop). There is no linear trend for the interaction, however, so that the estimate of the linear trend in wool A is not significantly different from the estimate of the linear trend in wool B.

Remember that summary.aov is using Type I sums of squares, so the estimates for some effects may not be what we want. In this example, the design is balanced and there are no missing data, so the SS estimates using Type I and Type III work out to be the same, but in your own data there may be a difference. Note that our orthogonal contrasts here are simple comparisons between means, and aren’t affected by the type of SS used. If you are concerned about Type of SS, you may want to grab the contrast estimates from this output and put them into your other sstable object. Here’s how you could do that:

Note which rows in the output correspond to the contrasts you want. In this case, it’s rows 3 and 4 for the contrasts on the main effect of tension, and rows 6 and 7 for the contrasts on the interaction. I select those rows with c(3, 4, 6, 7). I’m also selecting and reordering the columns in the output, so they’ll match what we have in sstable. I select the 2nd column (Sum Sq), then the first (Df), then the fourth (F value), then the fifth (Pr(>F)) with c(2, 1, 4, 5).

Remember that you can use [ , ] to select particular combinations of rows and columns from a given matrix or dataframe. Just put the rows you want as the first argument, and the columns as the second, i.e. [r, c]. If you leave either the rows or the columns blank, it will return all (so [r, ] will return row r and all columns).

# this pulls out just the specified rows
contrasts <- summary.aov(model, split = list(tension=list(L=1, Q=2)))[[1]][c(3, 4, 6, 7), c(2, 1, 4, 5)]

contrasts

##                    Sum Sq Df F value    Pr(>F)    
##   tension: L      1950.69  1 16.2979 0.0001938 ***
##   tension: Q        83.56  1  0.6982 0.4075366    
##   wool:tension: L  250.69  1  2.0945 0.1543266    
##   wool:tension: Q  752.08  1  6.2836 0.0156262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now use rbind to create a sort of Frankenstein table, splicing the contrasts estimate rows in with the other rows of sstable.

# select the rows to combine
maineffects <- sstable[c(1,2,3), ]
me_contrasts <- contrasts[c(1,2), ]
interaction <- sstable[4, ]
int_contrasts <- contrasts[c(3,4), ]
resid <- sstable[5, ]

# bind the rows together in the desired order
sstable <- rbind(maineffects, me_contrasts, interaction, int_contrasts, resid)

sstable # ta-da!

## Anova Table (Type III tests)
## 
## Response: breaks
##                   Sum Sq Df  F value    Pr(>F)    
## (Intercept)        42785  1 357.4672 < 2.2e-16 ***
## wool                 451  1   3.7653 0.0582130 .  
## tension             2034  2   8.4980 0.0006926 ***
##   tension: L        1951  1  16.2979 0.0001938 ***
##   tension: Q          84  1   0.6982 0.4075366    
## wool:tension        1003  2   4.1891 0.0210442 *  
##   wool:tension: L    251  1   2.0945 0.1543266    
##   wool:tension: Q    752  1   6.2836 0.0156262 *  
## Residuals           5745 48                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates of effect size

A popular measure of effect size for ANOVAs (and other linear models) is partial eta-squared. It’s the sums of squares for each effect divided by the error SS. The following code adds a column to the sstable object with partial eta-squared estimates for each effect:

sstable$pes <- c(sstable$'Sum Sq'[-nrow(sstable)], NA)/(sstable$'Sum Sq' + sstable$'Sum Sq'[nrow(sstable)]) # SS for each effect divided by the last SS (SS_residual)

sstable

## Anova Table (Type III tests)
## 
## Response: breaks
##                   Sum Sq Df  F value  Pr(>F)     pes
## (Intercept)        42785  1 357.4672 0.00000 0.88162
## wool                 451  1   3.7653 0.05821 0.07274
## tension             2034  2   8.4980 0.00069 0.26149
##   tension: L        1951  1  16.2979 0.00019 0.25348
##   tension: Q          84  1   0.6982 0.40754 0.01434
## wool:tension        1003  2   4.1891 0.02104 0.14861
##   wool:tension: L    251  1   2.0945 0.15433 0.04181
##   wool:tension: Q    752  1   6.2836 0.01563 0.11576
## Residuals           5745 48

Okay great! There’s your output, but you don’t want to just copy-paste that mess into your manuscript. Let’s get R to generate a nice, clean table we can use in Word.

Creating an APA style ANOVA table from R output

kable(sstable, digits = 3) # the digits argument controls rounding
Sum Sq Df F value Pr(>F) pes
(Intercept) 42785.185 1 357.467 0.000 0.882
wool 450.667 1 3.765 0.058 0.073
tension 2034.259 2 8.498 0.001 0.261
tension: L 1950.694 1 16.298 0.000 0.253
tension: Q 83.565 1 0.698 0.408 0.014
wool:tension 1002.778 2 4.189 0.021 0.149
wool:tension: L 250.694 1 2.095 0.154 0.042
wool:tension: Q 752.083 1 6.284 0.016 0.116
Residuals 5745.111 48 NA NA NA

Wait, what? That was so easy!

Yes, yes it was. :)

In a lot of cases, that will be all you need to get a workable ANOVA table in your document. Just for fun, though, let’s play around with customizing it a little.

Hide missing values

By default, kable displays missing values in a table as NA, but in this case we’d rather have them just be blank. You can control that with the options command:

options(knitr.kable.NA = '') # this will hide missing values in the kable table

kable(sstable, digits = 3)
Sum Sq Df F value Pr(>F) pes
(Intercept) 42785.185 1 357.467 0.000 0.882
wool 450.667 1 3.765 0.058 0.073
tension 2034.259 2 8.498 0.001 0.261
tension: L 1950.694 1 16.298 0.000 0.253
tension: Q 83.565 1 0.698 0.408 0.014
wool:tension 1002.778 2 4.189 0.021 0.149
wool:tension: L 250.694 1 2.095 0.154 0.042
wool:tension: Q 752.083 1 6.284 0.016 0.116
Residuals 5745.111 48

Add a table caption

You can add a title for the table if you change the format to “pandoc”. Depending on your final document output (pdf, html, Word, etc.), you can get automatic table numbering this way as well, which saves much time and many headaches.

kable(sstable, digits = 3, format = "pandoc", caption = "ANOVA table")
ANOVA table
Sum Sq Df F value Pr(>F) pes
(Intercept) 42785.185 1 357.467 0.000 0.882
wool 450.667 1 3.765 0.058 0.073
tension 2034.259 2 8.498 0.001 0.261
tension: L 1950.694 1 16.298 0.000 0.253
tension: Q 83.565 1 0.698 0.408 0.014
wool:tension 1002.778 2 4.189 0.021 0.149
wool:tension: L 250.694 1 2.095 0.154 0.042
wool:tension: Q 752.083 1 6.284 0.016 0.116
Residuals 5745.111 48

Modify column and row names

Often, the automatic row names and column names aren’t quite what you want. If so, you’ll need to modify them for the sstable object itself, and then run kable on the updated object.

colnames(sstable) <- c("SS", "df", "$F$", "$p$", "partial $\\eta^2$")

rownames(sstable) <- c("(Intercept)", "Wool", "Tension", "Tension: Linear Trend", "Tension: Quadratic Trend", "Wool x Tension", "Wool x Tension: Linear Trend", "Wool x Tension: Quadratic Trend", "Residuals")

kable(sstable, digits = 3, format = "pandoc", caption = "ANOVA table")
ANOVA table
SS df F p partial η2
(Intercept) 42785.185 1 357.467 0.000 0.882
Wool 450.667 1 3.765 0.058 0.073
Tension 2034.259 2 8.498 0.001 0.261
Tension: Linear Trend 1950.694 1 16.298 0.000 0.253
Tension: Quadratic Trend 83.565 1 0.698 0.408 0.014
Wool x Tension 1002.778 2 4.189 0.021 0.149
Wool x Tension: Linear Trend 250.694 1 2.095 0.154 0.042
Wool x Tension: Quadratic Trend 752.083 1 6.284 0.016 0.116
Residuals 5745.111 48

Omit the intercept row

For many models, the intercept is not of any theoretical interest, and you may want to omit it from the output. If you just want to drop one row (or column), the easiest approach is to indicate that row’s number and put a minus sign before it:

kable(sstable[-1, ], digits = 3, format = "pandoc", caption = "ANOVA table")
ANOVA table
SS df F p partial η2
Wool 450.667 1 3.765 0.058 0.073
Tension 2034.259 2 8.498 0.001 0.261
Tension: Linear Trend 1950.694 1 16.298 0.000 0.253
Tension: Quadratic Trend 83.565 1 0.698 0.408 0.014
Wool x Tension 1002.778 2 4.189 0.021 0.149
Wool x Tension: Linear Trend 250.694 1 2.095 0.154 0.042
Wool x Tension: Quadratic Trend 752.083 1 6.284 0.016 0.116
Residuals 5745.111 48

Inline APA-style output

You can also knit R output right into your typed sentences! To include inline R output, use back-ticks, like this:

Here's a sentence, and I want to let you know that the total number of cats I own is `r length(my_cats)`.

The back-ticks mark out the code to run, and the r after the first back-tick tells knitr that it’s R code (if you feel the need, you can incorporate code from pretty much any language you like, not just R). Assuming there’s a vector called my_cats, when we knit the document, that line of code will be evaluated and the result (the number of items in the vector my_cats) will be printed right in that sentence.

Let’s work this into our ANOVA reporting.

Here’s an example write-up of this ANOVA, using inline code to plug in the stats. Since the inline code can be a little hard to read, I like to save all of the variables I want to use inline with convenient names first.

fstat <- unname(summary(model)$fstatistic[1])
df_model <- unname(summary(model)$fstatistic[2])
df_res <- unname(summary(model)$fstatistic[3])
rsq <- summary(model)$r.squared
p <- pf(fstat, df_model, df_res, lower.tail = FALSE)

Then I can plug them in as needed in my writing.

A 2x2 factorial ANOVA revealed that tension (low, medium, or high) and wool type (A or B) predict a significant amount of variane in number of breaks, $R^2$=`r round(rsq, 2)`, $F(`r round(df_model, 0)`, `r round(df_res, 0)`)=`r round(fstat, 2)`$, $p=`r ifelse(round(p, 3) == 0, "<.001", round(p, 3))`$.

Here’s how that code renders when knit: A 2x2 factorial ANOVA revealed that tension (low, medium, or high) and wool type (A or B) predict a significant amount of variane in number of breaks, R2=0.38, F(5, 48)=5.83, *p* = <.001.

References and Further Reading

Coming soon. :)