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
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:

Meta-analysis summary plot

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.

Advertisements

6 thoughts on “R code for Meta-analysis summary plots

  1. Nice rarefaction curves! And I like the way you build the other plot up line by line in your example code. Any chance of some example code for plotting a rarefaction curve?
    Thanks 🙂

  2. Very cool! I find ggplot2 also works beautifully with the metafor package for meta-analysis. Simple extract coefficient and CI information and go! Although I often use geom_pointrange.

    1. I agree. I tend to use metafor as well, but just didn’t want to include it in this example as I just wanted to focus on the ggplot2 side of things. Given what you’ve said (and that stats posts seem to be very popular on here) I might knock up a post on using metafor for meta-analysis, I know at least a few people that would find it useful.

  3. Thank you so much! I am pretty new in meta-analysis and found this very useful. Now the groups for meta-analysis shows up alphabetically, is there a way to change that? Thank you!!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s