knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
The CGPfunctions package includes functions that I find useful for teaching
statistics especially to novices (as well as an opportunity to sharpen my own R
skills). I only write functions when I have a real need -- no theory -- just
help for actually practicing the art. They typically are not "new" methods but
rather wrappers around either base R or other packages and are very task
focused. This vignette covers one function from the package that tries to help
users (especially students) do one thing well by pulling together pieces from a
variety of places in R
. Plot2WayANOVA
, which as the name implies conducts a
2 way ANOVA and plots the results. I always try and find the right balance
between keeping the number of dependencies to a minimum and not reinventing the
wheel and writing functions that others have done for me. The function makes use
of the following non base r packages.
ggplot2
as the work horse for all the actual plotting car
for it's ability to compute Type II sums of squares, we'll address why
that's important in more detail later in the scenario. We'll also make use of
it's leveneTest
.sjstats
which takes out ANOVA table and gives us other important
information such as the effect sizes ($\eta^2$ and $\omega^2$ ) through use of
its anova_stats
function. Prior to this version I had been using my own local
function but this runs rings around what I could do.broomExtra::glance
will also help us grab very important results like $R^2$
and display themDescTool::PostHocTest
for accomplishing post hoc testsThe ANOVA (Analysis of Variance) family of statistical techniques allow us to
compare mean differences of one outcome (dependent) variable across two or more
groups (levels) of one or more independent variables (factor). It is also true
that ANOVA is a special case of the GLM or regression models so as the number of
levels increase it might make more sense to try one of those approaches. The 2
Way ANOVA allows for comparisons of mean differences across 2 independent
variables factors
with a varying numbers of levels
in each factor
.
If you prefer a more regression based approach with a very similar plotted
result I highly recommend the interactions
package which I was unaware of
until just recently.
It is available through CRAN.
The Plot2WayANOVA
function conducts a classic analysis of variance (ANOVA) in
a sane and defensible, albeit opinionated, manner, not necessarily the only one.
It's real strength (I hope) lies in the fact that it is pulled together in
one function and more importantly allows you to visualize the results
concurrently with no additional work.
Imagine that you are interested in understanding whether a car's fuel efficiency
(mpg) varies based upon the type of transmission (automatic or manual) and the
number of cylinders the engine has. Let's imagine that the mtcars
data set is
actually a random sample of 32 cars from different manufacturers and use the
mean mpg
grouped by am
and cyl
to help inform our thinking. While we
expect variation across our sample we're interested in whether the differences
between the means by grouping of transmission type and cylinders is
significantly different than what we would expect in random variation across the
data.
In simplistic terms we want to know whether am
matters, cyl
matters or if it
depends on the interaction of the two. It's this interaction term that typically
confuses novices or is difficult to "see". That's where a good interaction graph
can hopefully play a key role, and that's what the Plot2WayANOVA
focuses on.
There's no lack or tools or capabilities in base R or in the many packages to do this task. What this function tries to do is pull together the disparate pieces with a set of sane defaults and a simple interface to work with it. At its simplest you would require the library and then enter this command:
Plot2WayANOVA(formula = mpg ~ am * cyl, dataframe = mtcars)
which lays our
question out in R's vernacular with a formula and a dataframe. Optionally we
can specify a different confidence level and choose a line or a bar graph. Over
time the function has gained a plethora of formatting options.
"Under the hood", however there's a lot of nice features at work.
am
and cyl
is equal. There's some debate as to how much
imbalance is permissible. But you'll be given fair warning if there is any.aov
or
Anova
you'll be presented with information about effect sizes like eta
squared $\eta^2$. They're calculated and appended as additional columns.
If you're unfamiliar with them and want to know more especially where the
numbers come from I recommend a good introductory stats text. I recommend
Learning Statistics with R
LSR see Table 14-1
on page 432.# Install from CRAN install.packages("CGPfunctions") # Or the development version from GitHub # install.packages("devtools") devtools::install_github("ibecav/CGPfunctions")
then load the library.
library(CGPfunctions)
The only two required parameters are a formula and a dataframe (like) object. If we run the function in its simplest form here's what we get.
Plot2WayANOVA(formula = mpg ~ am * cyl, dataframe = mtcars)
In the console you'll receive a series of messages chronicling your progress
and any diagnostic information. In this case am
and cyl
are being coerced to
factors and you're being prompted to make sure that's what is intended.
Next you receive a warning because you have a very unbalanced design. There are only two 8 cylinder cars with a manual transmission and twelve 8 cylinder cars with automatics. Whereas there are eight 4 cylinders with manual and only three that are automatics. Imbalance in our design worries us for two reasons. One is that it causes a lack of statistical power and creates some math challenges in deciding how to divide up the sums of the squared differences. This data set causes the more troublesome worry. Are the number of cylinders and manual versus automatic related systematically which would call our whole design into question. Make sure you can answer questions about which is at work here, or make sure you have a balanced design.
A table follows that is intended to summarize the findings. We'll discuss it more later when we examine the plot.
The overall measures table can be very handy for comparing numerous models. For
example how does the AIC number change if we were to eliminate am
?
The table of group means is useful for looking at summary by group. Want the best gas mileage? Buy a 4 cylinder manual car.
In our simple example the only statistically significant effect is for the main effect of number of cylinders. Accordingly the Scheffe test is run against the three types of cars with 4, 6, or 8 cylinders and we can see that with a difference of 7.3 mpg eight and four cylinder cars are statistically significant even as we control for the multiple simultaneous comparisons.
The next step is to test homogeneity of variance also known as (homoscedasticity). Since the math in our ANOVA rely on the assumption that the variance for the different groupings of cars is more or less equal, we need to check that assumption.
We'll use the Brown-Forsythe test. When you run the leveneTest
in R the default
is actually a Brown-Forsythe, to get a true Levene you must specify center =
mean
. Brown-Forsythe is actually more robust since it tests differences from
the median. Not surprisingly when we consult out table of group results we have
some reason for concern sine the standard deviations vary widely.
Finally, let’s address the assumption that our errors or residuals are normally distributed. We’re looking for evidence that our residuals are skewed or tailed or otherwise misshapen in a way that would influence our results. Surprisingly, there is actually quite a bit of controversy on this point since on the one hand we have strong reason to believe that our sample will be imperfect and that our population will not necessarily be “perfectly normal” either. Some argue that some simple plotting is all that is necessary looking for an unspecifiable amount of non normality that will trigger a search for the source. Other prefer a more formal approach using one or more statistical tests.
Plot2WayANOVA
runs the most common test of the normality assumption (there are
many) the Shapiro-Wilk test The statistics look good, no strong evidence in the
data we have.
The default settings for the resultant plot are deliberately minimalistic, allowing you to focus visually on the pattern of means and the connecting lines. If you're already used to looking at this sort of plot it is immediately apparent from the separation between the lines that the number of cylinders is having a significant impact on mileage. Automatic versus manual transmission seems to have less impact (judged by the relative lack of slope except for 4 cylinder models) and there does seem to be at least the start of an interaction between the two. (bear in mind this is a small data set and we are very unbalanced).
One other easy tip is warranted. Order matters and sometimes it is helpful to
run the command simply reversing the order of the independent variables to help
you better "see" the results, e.g. Plot2WayANOVA(formula = mpg ~ cyl * am,
dataframe = mtcars)
Note that if you want to "save" all these tables of data and information all
you need to do is store the results in an object as in MyResults <-
Plot2WayANOVA(formula = mpg ~ am * cyl, dataframe = mtcars)
then MyResults
can be accessed as a list.
Let's make some changes that are likely to be quite common:
Plot2WayANOVA(formula = mpg ~ cyl * am, dataframe = mtcars, confidence = .99, title = "MPG by cylinders and type transmission", xlab = "Cylinders", ylab = "Miles per gallon", mean.label = TRUE, mean.shape = 22, posthoc.method = "lsd", errorbar.display = "SEM" )
Please don't fail to notice how liberal Fisher's LSD is compared to Scheffe especially given we've demanded more confidence
Although the defaults are minimalistic to allow you to focus on the interaction pattern the function also has any number of optional ways of increasing complexity and showing more information.
Let's make some custom changes that are more uncommon:
# Create a new dataset library(dplyr) library(ggplot2) library(stringi) newmpg <- mpg %>% filter(cyl != 5) %>% mutate(am = stringi::stri_extract(trans, regex = "auto|manual")) Plot2WayANOVA(formula = hwy ~ am * cyl, dataframe = newmpg, ylab = "Highway mileage", xlab = "Transmission type", plottype = "line", offset.style = "wide", overlay.type = "box", mean.label = TRUE, mean.shape = 20, mean.size = 3, mean.label.size = 3, show.dots = TRUE, errorbar.display = "SD", ggtheme = ggplot2::theme_minimal(), ggplot.component = theme(axis.text.x = element_text(size=14, color="darkblue")) )
Many thanks to Danielle Navarro and the book * Learning Statistics with R.*
If you like CGPfunctions, please consider Filing a GitHub issue by leaving feedback here, or by contacting me at ibecav at gmail.com by email.
I hope you've found this useful. I am always open to comments, corrections and suggestions.
Chuck
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.