Interaction plot from cell means

February 24, 2010
By

(This article was first published on Psychological Statistics, and kindly contributed to R-bloggers)

I needed to produce a few a interaction plots for my book in R and, while the interaction.plot() function is useful it has a couple of drawbacks. First, the default output isn't very pretty. Second, it works from the raw data, whereas I often need plots from cell means. For teaching purposes it is quite common to produce plots without raw data (for hypothetical data or from published examples).

My first attempts at the plots involved setting them up element by element. Just going over some examples I decided to turn the basic plot (for a 2 x 2 ANOVA) into a simple function. Nothing fancy, just a regular interaction plot in black and white that I think is prettier than the SPSS, Excel or R defaults. At some point I may have a go turning it into a general I x J ANOVA plot (or maybe even add CIs, but I'll probably do that from raw data if I ever get round to it).


plot.2by2 <- function(A1B1,A1B2, A2B1, A2B2, group.names, legend = TRUE, leg.loc=NULL, factor.labels=c('Factor A', 'Factor B'), swap = FALSE, ylab= NULL, main = NULL){
group.means <- c(A1B1, A2B1, A1B2, A2B2)
if(missing(ylab)) ylab <- expression(italic(DV))
if(swap==TRUE) {
group.names <- list(group.names[[2]], group.names[[1]]) ; group.means <- c(A1B1, A1B2, A2B1, A2B2); factor.labels <- c(factor.labels[2], factor.labels[1])
}
plot(group.means, pch=NA, ylim=c(min(group.means)*.95, max(group.means)*1.025), xlim=c(0.8,2.2), ylab=ylab, xaxt='n', xlab=factor.labels[1], main=main)
points(group.means[1:2], pch = 21)
points(group.means[3:4], pch = 19)
axis(side = 1, at = c(1:2), labels = group.names[[1]])
lines(group.means[1:2], lwd = .6, lty = 2)
lines(group.means[3:4], lwd = .6)
if(missing(leg.loc)) leg.loc <- c(1,max(group.means))
if(legend ==TRUE) legend(leg.loc[1], leg.loc[2],legend = group.names[[2]],  title = factor.labels[2], lty = c(3,1))
}


Call the function by entering the four cell means in conventional order: A1B1, A1B2 and so on where A1B1 is the mean of level 1 of factor A at level 1 of factor B. You also need a two item list containing text strings of the two level names of each factor. For instance:

lev.names <- list(c('A1', 'A2'), c('B1', 'B2'))
plot.2by2(5,15,10,20, lev.names)
You can swap the axes by adding the argument swap = TRUE:

plot.2by2(5,15,10,20, lev.names, swap = TRUE)



The default factor names are 'Factor A' and 'Factor B', but these are over-ridden in the call:

plot.2by2(5,15,10,20,lev.names, swap = TRUE, factor.labels= c('Factor 1','Factor 2'))


You can also change the y-axis label with ylab or add a main title with main.  The legend can be dropped (legend = FALSE) if you don't want one or need it to be located outside the plot. To move the legend just specify coordinates with an argument such as leg.loc = c(1,10). You can also edit the source code directly. Here is an example with title and meaningful labels:


group.names <- list(c('placebo','drug'), c('male', 'female'))


plot.2by2(10,10,15,20, group.names, factor.labels=c('Drug', 'Sex'), swap = FALSE)

As this just uses basic plotting functions in R you can also manipulate the plot in other ways: adding lines with segments(), adding text with text() changing graphical parameters with par() and so on. Depending on your platform it is also easy to extract the plot as a .pdf or .jpg file. On a mac I save it as a .pdf file and open it in preview which allows me to save it as .png, .gif or whatever I need.

To leave a comment for the author, please follow the link and comment on his blog: Psychological Statistics.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.