Conduct phylogenetic meta-analysis from R calling Phylometa.

phylometa_inR_run.R

#####################################################################
# Created by Scott Chamberlain
# Ecology and Evolutionary Biology Dept., Rice University
# Houston, TX 77005, USA
# myrmecocystus@gmail.com
# http://schamber.wordpress.com/
#####################################################################

########Directions
#Put your phylogeny, in format required by phylometa (see Phylometa manual), in your working directory
#Put your meta-analysis dataset, in format required by phylometa (see Phylometa manual), in your working directory
#Set working directory
#Use below functions
#Beware: only use a moderator variable with up to 4 groups 

########Install packages
install.packages(c("plyr","ggplot2")) # if not already installed
library(plyr); library(ggplot2)

########Set the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY]
setwd("Z:/Mac/R_stuff/Blog_etc/Phylometa_from_R")

#Call and run functions (used below) in the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY]
source("Z:/Mac/R_stuff/Blog_etc/Phylometa_from_R/phylometa_fxn.r")

###########################Functions to to a phylogenetic meta-analysis
#Define number of groups in moderator variable
groups <- 2

####Run phylometa. Change file names as needed
phylometa.run <- system(paste('"new_phyloMeta_1.2b.exe" phylogeny.txt metadata_2g.txt'),intern=T) 

####Process phylometa output
myoutput2g <- phylometa.process(phylometa.run,groups)

####Get output from phylometa.run
phylometa.output(myoutput2g) #Prints all five tables
phylometa.output.table(myoutput2g,2) #Prints the table you specify, from 1 to 5, in this example, table 2 is output

###################################################
#########Plot effect sizes. These are various ways to look at the data. Go through them to see what they do. Output pdf's are in your working directory
#Make table for plotting
analysis <- c(rep("fixed",groups+1),rep("random",groups+1))
trad_effsizes <- data.frame(analysis,phylometa.output.table(myoutput2g,2)) #Tradiational effect size table
phylog_effsizes <- data.frame(analysis,phylometa.output.table(myoutput2g,4)) #Phylogenetic effect size table

#The arrange method
limits <- aes(ymax = effsize + (CI_high-effsize), ymin = effsize - (effsize-CI_low))
dodge <- position_dodge(width=0.3)
plot01 <- ggplot(trad_effsizes,aes(y=effsize,x=analysis,colour=Group)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank(),title="Traditional meta-analysis") + labs(x="Group",y="Effect size") + geom_errorbar(limits, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2)
plot02 <- ggplot(phylog_effsizes,aes(y=effsize,x=analysis,colour=Group)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank(),title="Phylogenetic meta-analysis") + labs(x="Group",y="Effect size") + geom_errorbar(limits, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2)
pdf("plots_effsizes_arrange.pdf",width = 8, height = 11)
arrange_(plot01,plot02,ncol=1)
dev.off() 

#used in the two plotting methods below
bothanalyses<-data.frame(tradphy=c(rep("Traditional",(groups*2)+2),rep("Phylogenetic",(groups*2)+2)),fixrand=rep(analysis,2),rbind.fill(phylometa.output.table(myoutput2g,2),phylometa.output.table(myoutput2g,4))) #Table of both trad and phylo
limits2 <- aes(ymax = effsize + (CI_high-effsize), ymin = effsize - (effsize-CI_low))
dodge <- position_dodge(width=0.3)

#The grid/wrap method, version 1
plot03 <- ggplot(bothanalyses,aes(y=effsize,x=tradphy,colour=Group)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank()) + labs(x="Group",y="Effect size") + geom_errorbar(limits2, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2) + facet_grid(.~fixrand)
pdf("plots_effsizes_wrap1.pdf")
plot03
dev.off() 

#The grid/wrap method, version 2 (excuse the sloppy x-axis labels)
plot04 <- ggplot(bothanalyses,aes(y=effsize,x=Group,colour=tradphy)) + geom_point(size=3,position=dodge) + theme_bw() + opts(panel.grid.major = theme_blank(),panel.grid.minor=theme_blank()) + labs(x="Group",y="Effect size") + geom_errorbar(limits2, width=0.2, position=dodge) + geom_hline(yintercept=0,linetype=2) + facet_grid(.~fixrand)
pdf("plots_effsizes_wrap2.pdf")
plot04
dev.off()
Created by Pretty R at inside-R.org