How to draw an interaction plot with confidence intervals?

My attempts:

1. I couldn't get confidence intervals in interaction.plot()

2. and on the other hand plotmeans() from package 'gplot' wouldn't display two graphs. Furthermore, I couldn't impose two plotmeans() graphs one on top of the other because by default the axis are different.

3. I had some success using plotCI() from package 'gplot' and superimposing two graphs but still the match of the axis wasn't perfect.

Any advice on how to make an interaction plot with confidence intervals? Either by one function, or advice on how to superimpose plotmeans() or plotCI() graphs.

code sample

br=structure(list(tangle = c(140L, 50L, 40L, 140L, 90L, 70L, 110L,
150L, 150L, 110L, 110L, 50L, 90L, 140L, 110L, 50L, 60L, 40L,
40L, 130L, 120L, 140L, 70L, 50L, 140L, 120L, 130L, 50L, 40L,
80L, 140L, 100L, 60L, 70L, 50L, 60L, 60L, 130L, 40L, 130L, 100L,
70L, 110L, 80L, 120L, 110L, 40L, 100L, 40L, 60L, 120L, 120L,
70L, 80L, 130L, 60L, 100L, 100L, 60L, 70L, 90L, 100L, 140L, 70L,
100L, 90L, 130L, 70L, 130L, 40L, 80L, 130L, 150L, 110L, 120L,
140L, 90L, 60L, 90L, 80L, 120L, 150L, 90L, 150L, 50L, 50L, 100L,
150L, 80L, 90L, 110L, 150L, 150L, 120L, 80L, 80L), gtangles = c(141L,
58L, 44L, 154L, 120L, 90L, 128L, 147L, 147L, 120L, 127L, 66L,
118L, 141L, 111L, 59L, 72L, 45L, 52L, 144L, 139L, 143L, 73L,
59L, 148L, 141L, 135L, 63L, 51L, 88L, 147L, 110L, 68L, 78L, 63L,
64L, 70L, 133L, 49L, 129L, 100L, 78L, 128L, 91L, 121L, 109L,
48L, 113L, 50L, 68L, 135L, 120L, 85L, 97L, 136L, 59L, 112L, 103L,
62L, 87L, 92L, 116L, 141L, 70L, 121L, 92L, 137L, 85L, 117L, 51L,
84L, 128L, 162L, 102L, 127L, 151L, 115L, 57L, 93L, 92L, 117L,
140L, 95L, 159L, 57L, 65L, 130L, 152L, 90L, 117L, 116L, 147L,
140L, 116L, 98L, 95L), up = c(-1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
-1L, -1L, 1L, 1L, 1L, 1L, -1L, -1L, -1L, -1L, 1L, 1L, -1L, -1L,
1L, 1L, -1L, 1L, 1L, -1L, 1L, 1L, 1L, 1L, 1L, -1L, -1L, 1L, 1L,
1L, 1L, -1L, -1L, 1L, 1L, -1L, -1L, -1L, -1L, -1L, -1L, -1L,
1L, -1L, -1L, -1L, -1L, -1L, 1L, -1L, 1L, 1L, -1L, -1L, -1L,
-1L, 1L, -1L, 1L, -1L, -1L, -1L, 1L, -1L, 1L, -1L, 1L, 1L, 1L,
-1L, -1L, -1L, -1L, -1L, -1L, 1L, -1L, 1L, 1L, -1L, -1L, 1L,
1L, 1L, -1L, 1L, 1L, 1L)), .Names = c("tangle", "gtangles", "up"
), class = "data.frame", row.names = c(NA, -96L))

plotmeans2 <- function(br, alph) {
dt=br;   tmp   <- split(br$gtangles, br$tangle);
means <- sapply(tmp, mean);  stdev <- sqrt(sapply(tmp, var));
n <- sapply(tmp,length);
ciw   <- qt(alph, n) * stdev / sqrt(n)
plotCI(x=means, uiw=ciw, col="black", barcol="blue", lwd=1,ylim=c(40,150),  xlim=c(1,12));
par(new=TRUE) dt= subset(br,up==1);
tmp   <- split(dt$gtangles, dt$tangle);
means <- sapply(tmp, mean);
stdev <- sqrt(sapply(tmp, var));
n <- sapply(tmp,length);
ciw  <- qt(0.95, n) * stdev / sqrt(n)
plotCI(x=means, uiw=ciw, type='l',col="black", barcol="red", lwd=1,ylim=c(40,150), xlim=c(1,12),pch='+');
abline(v=6);abline(h=90);abline(30,10); par(new=TRUE);
dt=subset(br,up==-1);
tmp <- split(dt$gtangles, dt$tangle);
means <- sapply(tmp, mean);
stdev <- sqrt(sapply(tmp, var));
n <- sapply(tmp,length);
ciw <- qt(0.95, n) * stdev / sqrt(n)
plotCI(x=means, uiw=ciw, type='l', col="black", barcol="blue",   lwd=1,ylim=c(40,150), xlim=c(1,12),pch='-');abline(v=6);abline(h=90);
abline(30,10);
}

plotmeans2(br,.95)


If you're willing to use ggplot, you can try the following code.

With a continuous predictor

library(ggplot2)
gp <- ggplot(data=br, aes(x=tangle, y=gtangles))
gp + geom_point() + stat_smooth(method="lm", fullrange=T) + facet_grid(. ~ up)


for a facetted interaction plot

For a standard interaction plot (like the one produced by interaction.plot()), you just have to remove the facetting.

gp <- ggplot(data=br, aes(x=tangle, y=gtangles, colour=factor(up)))
gp + geom_point() + stat_smooth(method="lm")


With a discrete predictor

Using the ToothGrowth dataset (see help(ToothGrowth)),

ToothGrowth$dose.cat <- factor(ToothGrowth$dose, labels=paste("d", 1:3, sep=""))
df <- with(ToothGrowth , aggregate(len, list(supp=supp, dose=dose.cat), mean))
df\$se <- with(ToothGrowth , aggregate(len, list(supp=supp, dose=dose.cat),
function(x) sd(x)/sqrt(10)))[,3]

opar <- theme_update(panel.grid.major = theme_blank(),
panel.grid.minor = theme_blank(),
panel.background = theme_rect(colour = "black"))
gp <- ggplot(df, aes(x=dose, y=x, colour=supp, group=supp))
gp + geom_line(aes(linetype=supp), size=.6) +
geom_point(aes(shape=supp), size=3) +
geom_errorbar(aes(ymax=x+se, ymin=x-se), width=.1)
theme_set(opar)


• Thank you so much for the detailed response. I wanted to ask, is there a way to make vertical confidence intervals at each level of the independent variable? Is there a way to remove the background and revert to 'old style' graph? Apr 20, 2011 at 16:08
• @Adam I updated my response with the case of 2 categorical variables + a continuous response variable -- hope this is what you meant. I also added code to show how to customize ggplot theme. Generally, you can say gp + theme_bw() to just remove the grey background; here, I also removed the grid.
– chl
Apr 20, 2011 at 20:05

There's also Fox and Hong's effects package in R. See the J. Stat. Soft. papers here and here for examples with confidence intervals and generating R code.

It's not quite as pretty as a ggplot solution, but quite a bit more general, and a lifesaver for moderately complex GLMs.

• (+1) I must admit I prefer this approach :-)
– chl
Apr 13, 2011 at 20:39
• @chl and/or Conjugate, can you say more about why you prefer this approach? It would help people like me decide which method to invest time in. Dec 9, 2011 at 23:43
• @MichaelBishop Essentially because it wraps up a lot of tricky things (plotting on link vs. response scale, displaying 95% CI for GLMMM, marginalization against interaction terms, etc.) that would be hard to handle in few R commands (and personally, I very much like lattice graphics :)
– chl
Dec 10, 2011 at 18:36