example regression plots with pretty CI shades

regplots.R

require(RColorBrewer)

# RColorBrewer provides the palettes developped by Cynthia Brewer, you can check them
# out at http://colorbrewer2.org/ .

# make some mock data with interaction

x1 <- factor(rep(c('a','b','c'),each=33))
x2 <- rnorm(99,10,2)

y <- c(
		3*x2[1:33]+rnorm(33,-10,1),
		1*x2[34:66]+rnorm(33,5,1),
		-0.5*x2[67:99]+rnorm(33,10,1)
		)
		
dat <- data.frame(y,x1,x2)
rm(y,x1,x2)
		
plot(y~x2, data=dat,col=factor(x1))

# fit a linear model

m1 <- lm(y~x1*x2, data=dat)

summary(m1)

# get model estimates

#tip: put the categorical last so it'll be more organized

newdat = expand.grid(
	x2 = seq(min(dat$x2),max(dat$x2),length=50),
	x1 = levels(dat$x1),
	y = 0 )


head(newdat)

pred = predict(m1,newdat,type='terms', interval='confidence')

# to get the predicted values
# the estimates are in pred$fit, but you need to add the global intercept to it.
# the global intercept is stored in attr(pred$fit,'constant')

# from the newdat we created, you can see that the first 1/3 of the data is for level "a",
# the second third for "b", and the final third for "c".
 
estim <- rowSums(pred$fit)+attr(pred$fit,"constant")

# remember these data refere to the "newdat" we created, so we need to plot accordingly

x.estim <- rep(seq(min(dat$x2),max(dat$x2),length=50),3)

# now for the confidence intervals
# to plot, we will use polygon() the tricky part here is that the polygon is traced along
# the sequence you provide! So you need to tell it to:
# start from the "top left", go to the "top right", then to the "bottom right", and finally
# back to the "bottom left" !!!! Or else it will look like a mess.
#
# Also, for the linear part you could just use the limits with segments(), but not here.
# this is the reason we did the predict with a bunch of points: to capture the "curviness"
# of the confidence interval.

CI.a <- c (
			rowSums(pred$upr[1:50,]), #top part
			rowSums(pred$lwr[50:1,])  #NOTE REVERSE ORDER! will trace in the right order
			) + attr(pred$upr,'constant')

CI.b <- c (
			rowSums(pred$upr[51:100,]), #top part
			rowSums(pred$lwr[100:51,])  #NOTE REVERSE ORDER! will trace in the right order
			) + attr(pred$upr,'constant')

CI.c <- c (
			rowSums(pred$upr[101:150,]), #top part
			rowSums(pred$lwr[150:101,])  #NOTE REVERSE ORDER! will trace in the right order
			) + attr(pred$upr,'constant')


CI.x <- c(
		seq(min(dat$x2),max(dat$x2),length=50),
		seq(max(dat$x2),min(dat$x2),length=50) #again, reverse order
		)

######
#PLOT#
######

# let's get some nice colors from colorbrewer

cols = brewer.pal(3,'Set1')

# but we want transparency, so we need to convert it to RGB and to proportion

cols = col2rgb(cols)/255

# start with a blank plot, adjust axes and labels according to need

plot(y~x2, data=dat,xlim=c(min(dat$x2)-1,max(dat$x2)+1), ylim=c(min(dat$y)-2.5,max(dat$y+2.5)), type='n')

# you want to plot the intervals first so the shading won't obscure the lines

polygon(CI.a~CI.x, border=NA, col=rgb(cols[1,1],cols[2,1],cols[3,1], alpha=0.4))
polygon(CI.b~CI.x, border=NA, col=rgb(cols[1,2],cols[2,2],cols[3,2], alpha=0.4))
polygon(CI.c~CI.x,border=NA, col=rgb(cols[1,3],cols[2,3],cols[3,3], alpha=0.4))

# you need to plot the 3 lines separately. 
# in a linear model you only really need the info for the extreme values of x2, which you
# can plot with segments(), but to make this more broadly applicable, I will use lines()
# especially if you're creating vector graphics, you'd like to use segments (minimize size)

lines(estim[1:50]~x.estim[1:50])
lines(estim[51:100]~x.estim[51:100], lty=2)
lines(estim[101:150]~x.estim[101:150], lty=4)

# and finally, the legend

legend('topleft', legend = c('pairs','males','mixed'), bty='n',
	pch=22, pt.cex=3, pt.lwd=0, lty=c(1,2,4), seg.len=1.2,
			pt.bg= c(rgb(cols[1,1],cols[2,1],cols[3,1], alpha=0.4),
					 rgb(cols[1,2],cols[2,2],cols[3,2], alpha=0.4),
					 rgb(cols[1,3],cols[2,3],cols[3,3], alpha=0.4)
					)
)