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? –  Adam SA Apr 20 '11 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 '11 at 20:05
@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 '11 at 18:36