Phylometa was created by Marc Lajeunesse at University of South Florida, and is described in his 2009 AmNat paper. Phylometa can be downloaded free here.
Save phylometa_fxn.R (get here) to your working directory. Then use the block of code below to call the functions within phylometa_fxn.R.
The program Phylometa needs to be in the working directory you are calling from.
Let me know what doesn't work, and what improvements can be made; I'm sure there are many!
---This code below is also available here on Github.
########Directions
#Place phylometa software to your working directory
#Put your phylogeny, in format required by phylometa, in your working directory #Put your meta-analysis dataset, in format required by phylometa, in your working directory #Set working directory #Use below functions #Beware: only use a moderator variable with up to 6 groups ########Install packages install.packages(c("plyr","ggplot2")) library(plyr) library(ggplot2) ########Set the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY] setwd("/Users/Scott/Documents/phylometa") #Call and run functions (used below) in the working directory [NOTE:CHANGE TO YOUR WORKING DIRECTORY] source("/Users/Scott/Documents/phylometa")
###########################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('"phyloMeta_v1-2_beta.exe" phylogeny.txt metadata_2g.txt'),intern=T) ####Process phylometa output #E.g. myoutput <- phylometa.process(phylometa.run,groups) ####Get output from phylometa.run phylometa.output(myoutput) #Prints all five tables phylometa.output.table(myoutput,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(myoutput,2)) #Tradiational effect size table phylog_effsizes <- data.frame(analysis,phylometa.output.table(myoutput,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(myoutput,2),phylometa.output.table(myoutput,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()
Below is an example output figure from the code. This example is from an analysis using 5 groups (i.e., 5 levels in the explanatory variable).
No comments:
Post a Comment