R code for Meta-analysis summary plots

I think meta-analysis is great. I am aware that some out there are a bit less complementary about it, but I like it. Sure, it can be used stupidly and give spurious results, but then so can ANOVAs or GLMs if you don’t know what you’re doing.

I also love R and have recently been converted to using ggplot2, one of the fantastic packages put together by Hadley Wickham. For me ggplot2 is the most aesthetically pleasing way to plot your data in R. It also has the ability to draw error bars easily, which in base graphics can be a bit of a faff. It’s endlessly customisable and can produce things of beautiful simplicity. I was quite proud of the rarefaction curves I produced the other day for example:

Some rarefaction curves we came up with the other day

Anyway, I recently did a meta-analysis and found the best way to plot the results was in ggplot2. It wasn’t really easy to get into the form commonly used for summary plots so I thought I’d stick up some sample code here in the hope that it might help some people stuck with the same problems. The data I used here is all made up but the code for the graphs should prove useful.

#this is a file to create the standard meta-analysis 
#summary plot using ggplot2

#first load in/create some data, this should be in the form of 
#summarised means and upper and lower confidence intervals

#these are the mean values
means<-as.numeric(c(-1,0,1,0.5,-0.5)) 
#and this is a bit of code to come up with randomly sized CIs
ci<-c(rnorm(1,mean=-1,sd=3),rnorm(1,mean=0,sd=3),rnorm(1,mean=1,sd=3),rnorm(1,mean=0.5,sd=3),rnorm(1,mean=-0.5,sd=3))
ci2<-sqrt(ci^2)
#this is the upper 95% confidence interval
upper<-as.numeric(means+ci2) 
#this is the lower 95% confidence interval
lower<-as.numeric(means-ci2) 
# and this is a buch of different taxonomic groups for my imaginary meta-analysis
treatment<-c("Birds","Mammals","Amphibians", "Insects", "Plants") 
#stick all the data together
meta_test<-data.frame(cbind(as.numeric(means),as.numeric(upper),as.numeric(lower),treatment)) 

#now install and run ggplot from the library
install.packages("ggplot2")
library(ggplot2)

#the easiest way for me to make a ggplot figure is to build things up
#a bit at a time. That way if you only need to change one bit you can do it
#without hunting through a heap of code

#this defines the elements to go in the plot, both the x and y and upper and lower CIs
a<-ggplot(meta_test,aes(x=treatment,y=means,ymax=upper,ymin=lower,size=2))
#this defines the plot type
b<-a+geom_pointrange()
#this flips the co-ordinates so your x axis becomes your y and vice versa
c<-b+coord_flip()+scale_area(range=c(1.5))
#this puts in a dotted line at the point of group difference
d<-c+geom_hline(aes(x=0), lty=2,size=1)
#all of this gets rid of the grey grid and legends
e<-d+opts(panel.grid.minor=theme_blank(), panel.grid.major=theme_blank())+opts(axis.ticks = theme_blank(), axis.text.x = theme_blank())+ theme_bw()+opts(legend.position = "none")
#this sets x and y axis titles
f<-e+ xlab('Taxonomic group') +ylab ('Change following removal of herbivores')
#this sets axis label size
g<-f+opts(axis.text.x = theme_text(size = 16, colour = 'black')) +opts(axis.text.y = theme_text(size = 16, colour = 'black'))
#this sets axis title size and there is your finished summary plot!
g+opts(axis.title.x = theme_text(size = 20, colour = 'black'))+opts(axis.title.y = theme_text(size = 20, colour = 'black'))

Created by Pretty R at inside-R.org

At the end you should have something that looks a bit like this:

though your error bars will obviously be a different width.

Hope someone finds this useful. Drop me a comment if you did.

EDIT!

I have been asked below to supply the code for the rarefaction curves. The data for the curves was produced using a program other than R (shock horror – do they even exist?). I think my friend used estimateS. Anyway once you have done that you can get the data in R and do something like this to it:

#script for drawing rarefraction curves#

#read in data#
leah<-read.csv("C:/Users/Phil/Documents/My Dropbox/rare2.csv")

#load ggplot2 librar (or install it if you don't have it yet!)
library(ggplot2)

#this sorts the orders of the seperate panels for each plot
leah$Order<-sort(rep((1:20),100))
leah$sampleNo2<-reorder(leah$sampleNo,leah$Order)

#this tells ggplot what data to plot and the limits for the 'ribbon' around the curve
a<-ggplot(leah,aes(x=Samplesize,y=Mean,ymin=Mean-(1.96*SE),ymax=Mean+(1.96*SE)))

#this plots your line along with the ribbon and different 'facets' or panels and sets the y axis limits
b<-a+geom_line(shape=16,size=0.5)+geom_ribbon(colour=NA,fill="blue",alpha=0.2)+facet_wrap(~sampleNo2,scales="free_x")+ylim(0,80)+opts(panel.grid.major = theme_line(colour = "white"))
#this changes the colour of your plots to white and gets rid of gridlines
c<-b+ theme_bw()+opts(panel.grid.major = theme_line(colour =NA))
#this changes the angle of your x axis labels
d<-c+opts(axis.text.x=theme_text(angle=90, hjust=1))
#this changes the size of your x axis labels
e<-d+opts(axis.title.x = theme_text(size = 20, colour = 'black'))+opts(axis.title.y = theme_text(angle=90,size = 20, colour = 'black'))
#this changes you x axis title names
f<-e+ylab ('Mean richness')+xlab ('Sample size')
#this plots the final plot
f

#and this saves it
setwd("C:/Documents and Settings/PMART/My Documents/Dropbox/")
ggsave("rarefracation.png",height=6.75,width=9,dpi=300)

Created by Pretty R at inside-R.org

The data we used can be found here.