Post-doctoral researcher
German Centre for Integrative Biodiversity Research
Personal website

Brief introduction

Recent advances in meta analysis, such as mixed effects and structural equation models, allow ecologists to tackle increasingly more complex questions. Visualizing these results, however, is not straightforward. In a previous post, I described how to set up the data set and how to analyse it using a meta-analytical mixed-effects model. Here, I will briefly show how to graph results of bi-variate random effects model using ggplot2 and metafor.

I strongly recommend that users thoroughly read Wolfgang Viechtbauer’s web site on meta-analysis in R, which has a lot of useful information about common problems in meta-analysis.

For graphing results from ‘normal’ mixed effects models check out a) ben bolker’s FAQ and b) Lionel Hertzog’s blog

Getting started

Install and then load packages:

#install.packages(c("metafor","ggplot2"))

require(metafor) # for fitting meta-regression models
## Loading required package: metafor
## Loading 'metafor' package (version 1.9-5). For an overview 
## and introduction to the package please type: help(metafor).
require(ggplot2) # for visualising results
## Loading required package: ggplot2

For this tutorial, we’ll be using data from the Curtis database on plant responses to elevated C02. To save space and not to replicate code from the previous post, I’ve uploaded the data to Dropbox

Load the data

dat<-read.delim("Curtis_PN_GS_data.csv",sep=",",header=T)
str(dat)
## 'data.frame':    27 obs. of  7 variables:
##  $ PAP_NO   : int  121 121 121 121 168 2035 2035 2047 2065 2122 ...
##  $ GENUS    : Factor w/ 18 levels "ACER","ALNUS",..: 1 8 16 16 14 7 16 9 12 2 ...
##  $ SPECIES  : Factor w/ 23 levels "ALBA","ATLANTICA",..: 20 5 17 19 9 23 1 4 22 8 ...
##  $ GS_LRR   : num  0 0.0541 0.0715 -0.2231 -0.8535 ...
##  $ PN_LRR   : num  0.541 0.389 0.731 0.188 0.545 ...
##  $ GS_LRRvar: num  0.286 0.333 0.333 0.333 0.333 ...
##  $ PN_LRRvar: num  0.286 0.333 0.333 0.333 0.333 ...

Fit candidate models

gas.1<-rma.mv(yi=PN_LRR,V=PN_LRRvar, mods=~GS_LRR,random= ~PAP_NO,method="REML",digits=4,data=dat) # Study as random effect
gas.2<-rma.mv(yi=PN_LRR,V=PN_LRRvar, mods=~GS_LRR,random= ~GENUS|PAP_NO,method="REML",digits=4,data=dat) #SPP nested in Study as random effect


aa<-fitstats(gas.1)[5,1]
bb<-fitstats(gas.2)[5,1]

AICC<-data.frame(c(aa,bb))
mod<-data.frame(c("a","b"))
AICC<-cbind(mod,AICC)
AICC<-AICC[order(AICC[,2]),]
colnames(AICC)[1]<-"Mod"
colnames(AICC)[2]<-"AICc"
AICC$deltaAICc<-AICC$AICc-AICC$AICc[1]
AICC
##   Mod     AICc deltaAICc
## 1   a 32.10565  0.000000
## 2   b 34.96280  2.857143
# gas.1 has the most parsimonious random effects structure (lowest AICc)

Model summary

summary(gas.1)
## 
## Multivariate Meta-Analysis Model (k = 27; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc  
## -12.4814   24.9628   30.9628   34.6194   32.1057  
## 
## Variance Components: 
## 
##             estim    sqrt  nlvls  fixed  factor
## sigma^2    0.0000  0.0000     14     no  PAP_NO
## 
## Test for Residual Heterogeneity: 
## QE(df = 25) = 5.6089, p-val = 1.0000
## 
## Test of Moderators (coefficient(s) 2): 
## QM(df = 1) = 3.0809, p-val = 0.0792
## 
## Model Results:
## 
##          estimate      se    zval    pval    ci.lb   ci.ub     
## intrcpt    0.4242  0.1211  3.5042  0.0005   0.1869  0.6615  ***
## GS_LRR     0.4404  0.2509  1.7552  0.0792  -0.0514  0.9322    .
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Set up the prediction matrix

This step creates a data frame with predicted values from the random effects model,
which we then use to graph the predicted line.

preds<-predict(gas.1,levels=0, addx=TRUE,transf=exp) # let's transfer LRR to RR
preds<-do.call(cbind.data.frame, preds)

preds$GS_RR<-exp(preds[,8]) # convert log response ratios to response ratios

dat$GS_RR <- exp(dat$GS_LRR) # convert log response ratios to response ratios
dat$PN_RR <- exp(dat$PN_LRR) # convert log response ratios to response ratios

Make the graph!

One of the nice features about the ‘predict’ function in metafor is that it allows for plotting either confidence intervals or prediction intervals

g0 <- ggplot(preds,aes(x=GS_RR,y=pred))+ 
  stat_smooth(method="lm",formula=y~x,fullrange=T,se=FALSE,size=1,colour="black")+
  stat_smooth(data=preds,aes(x=GS_RR,y=ci.lb),method="lm",formula=y~x,fullrange=T,  
              se=FALSE,size=1,linetype=3,colour="black")+
  stat_smooth(data=preds,aes(x=GS_RR,y=ci.ub),method="lm",formula=y~x,fullrange=T,
              se=FALSE,size=1,linetype=3,colour="black")+
  
  geom_point(data=dat,aes(x=GS_RR,y=PN_RR),colour="black",fill="white",
             position="jitter",alpha=0.4,shape= 20,size=3)+
  
  scale_x_continuous(breaks=c(0,0.5,1,1.5,2,2.5,3))+scale_y_continuous(breaks=c(0,0.5,1,1.5,2,2.5,3))+
  
  labs(x=expression(bold(paste("Response of Gs to elevated CO"[2]))),
  y =expression(bold(paste("Response of A to elevated CO"[2]))))

Additional thoughts

One can also graph separate lines and adjust other parameters , e.g.,point size, shape, and color,
for each level of a predictor variable

Here’s an example:


g0 <- ggplot(preds,aes(x=Gs,y=pred,fill=ANGIO_CONIF,linetype=ANGIO_CONIF))+ 
 
  stat_smooth(method="lm",formula=y~x,fullrange=T,se=FALSE,colour="red",size=1)+
  geom_point(data=dat,aes(x=Gs,y=RR,fill=ANGIO_CONIF),colour="red",position="jitter",
        alpha=0.4,shape= 21,cex=3)+
  # this will change the color with which each point is filled based on 'ANGIO_CONIF'  
  # and the outline of the circle will be red
  
  scale_fill_manual(name="Angiosperm or Coniferous",
                values=c("Angiosperm"="red","Coniferous"="white")) +
  # this will specify the color used to fill the points
  
  scale_linetype_manual(name="Angiosperm or Coniferous",values=c("Angiosperm"=2,"Coniferous"=1))+
  # this will specify the linetype used to fill the points

To select colors that will look good on a printed page, check out Colorbrewer2

To adjust shapes and linetypes, check out this the R cookbook

Adjust ggplot2 theme parameters

Bivar<-g0+ theme(axis.title.x=element_text(colour="black",face="bold",size=10,vjust=-.5),
                    axis.title.y=element_text(colour="black",face="bold",size=10,vjust=2),
                    axis.text.y=element_text(colour="black",face="bold",size=8),
                    axis.text.x=element_text(colour="black",face="bold",size=8),
                    panel.background =element_rect(fill="transparent",colour="black"),
                    panel.border=element_rect(fill=NA,colour="black"))

The size and resolution can be adjusted to make publication-quality figures

png(filename="PS_GS_response.png", 
    type="cairo",
    units="in", 
    width=3.5, 
    height=3.5, 
    pointsize=2, 
    res=500)
Bivar

dev.off()
## png 
##   2

Now admire your hardwork!

alt text