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.

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.

This tutorial assumes…

That you already have a basic understanding of what an
ANOVA is and
when you might use it.

That you’re not brand new to
R.
If you are, the descriptions may still be useful to you, but you may
run into problems replicating the analysis on your own computer or
editing the code to suit your needs.

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()

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):

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.

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).

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.

(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).

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:

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.

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:

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.

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.