Getting Started with Meta-Analysis in R
Dylan Craven (http://dylancraven.weebly.com/)
Overview
- Organize data
- Calculate size effects
- Fit meta-regression models
- Test for publication bias
This is a gentle introduction to meta-analysis in R for ecologists, but is by no means exhaustive. I would encourage those interested to consult recent books on meta-analysis in ecology (e.g. Gurevitch et al. 2013) as well as the J. of Statistical Software article that accompanies the ‘metafor’ package (??metafor) as well as the package site: http://www.metafor-project.org/doku.php
We’ll focus on data from the Curtis (1996, 1998) database on plant responses to elevated C02, first in terms of plant total weight (“TOTWT”) and then we will look at how the relationship between photosynthesis (“PN”) and stomatal conductance (“GS”) changes.
- First, download the file (‘Curtis_CO2_database.csv’ to your local directory)
- Available via the NCEAS reposository of meta-analysis data: https://www.nceas.ucsb.edu/meta/publications.html#d_t_t
Organize data for analysis
Here, we will select relevant subsets of the data set, one for looking at plant biomass and the other for physiological characteristics.
Side note
reshape2 is an excellent package that makes it easy to prepare data for analysis
- ‘dcast’ and ‘melt’ are the main functions
- Look here for a nice tutorial: http://goo.gl/FfzOUL
require(metafor) # calculate size effects and perform meta analysis## Loading required package: metafor
## Loading required package: Formula
## Loading required package: Matrix
##
## Loading 'metafor' package (version 1.9-3). For an overview
## and introduction to the package please type: help(metafor).require(ggplot2) # pretty graphs## Loading required package: ggplot2setwd("/home/dylan/Documents/Post_Docs/Meta_Analysis_Workshop")
dat<-read.delim("Curtis_CO2_database.csv",sep=",",header=T)
dat1<-subset(dat,XTRT=="NONE")
dat1$XTRT <- droplevels(dat1$XTRT)
wt<-subset(dat1,PARAM=="TOTWT")
wt$PARAM<-droplevels(wt$PARAM)
dim(wt)## [1] 30 29head(wt)## OBSNO PAP_NO PARAM P_UNIT GENUS SPECIES DIV1 DIV2 AMBC ELEV
## 27 27 121 TOTWT g ACER RUBRUM WOODY ANGIO 350 700
## 32 32 121 TOTWT g QUERCUS PRINUS WOODY ANGIO 350 700
## 35 35 121 TOTWT g MALUS DOMESTICA WOODY ANGIO 350 700
## 38 38 121 TOTWT g ACER SACCHARINUM WOODY ANGIO 350 700
## 44 44 159 TOTWT g CASTANEA SATIVA WOODY ANGIO 350 700
## 63 63 183 TOTWT g CITRUS SINENSIS WOODY ANGIO 395 795
## CO2_UNIT TIME POT METHOD STOCK XTRT LEVEL QUANT SOURCE X_AMB SE_AMB
## 27 ppm 59 2.6 GH SEED NONE . . T4 1.93 0.2469
## 32 ppm 70 2.6 GH SEED NONE . . T4 6.62 0.7294
## 35 ppm 64 2.6 GH SEED NONE . . T4 4.10 0.6285
## 38 ppm 50 2.6 GH SEED NONE . . F2 6.42 1.1700
## 44 ppm 730 GRND GC SAP NONE . . T1 127.30 27.4000
## 63 ppm 365 9 GH SAP NONE . . T1 1140.60 47.9000
## SD_AMB CV._AMB N_AMB X_ELEV SE_ELEV SD_ELEV CV._ELEV N_ELEV
## 27 0.552 30.03 5 2.99 0.3828 0.856 30.06 5
## 32 1.631 25.87 5 5.91 0.7790 1.742 30.95 5
## 35 1.257 32.57 4 4.61 0.7035 1.407 32.43 4
## 38 2.026 34.19 3 10.78 0.5200 1.163 11.33 5
## 44 47.458 40.39 3 153.50 15.7000 27.193 19.19 3
## 63 82.965 7.88 3 1439.00 82.0000 142.028 10.69 3gas<-subset(dat1,PARAM=="PN" | PARAM=="GS") ##
gas$PARAM<-droplevels(gas$PARAM)
dim(gas)## [1] 80 29head(gas)## OBSNO PAP_NO PARAM P_UNIT GENUS SPECIES DIV1 DIV2 AMBC ELEV
## 28 28 121 GS molH2O/m2/s QUERCUS ROBUR WOODY ANGIO 350 700
## 29 29 121 PN umolCO2/m2/s QUERCUS ROBUR WOODY ANGIO 350 700
## 30 30 121 GS molH2O/m2/s QUERCUS PRINUS WOODY ANGIO 350 700
## 31 31 121 PN umolCO2/m2/s QUERCUS PRINUS WOODY ANGIO 350 700
## 33 33 121 GS molH2O/m2/s MALUS DOMESTICA WOODY ANGIO 350 700
## 34 34 121 PN umolCO2/m2/s MALUS DOMESTICA WOODY ANGIO 350 700
## CO2_UNIT TIME POT METHOD STOCK XTRT LEVEL QUANT SOURCE X_AMB SE_AMB
## 28 ppm 79 2.6 GH SEED NONE . . T2 0.100 0.015
## 29 ppm 79 2.6 GH SEED NONE . . T2 5.800 0.620
## 30 ppm 76 2.6 GH SEED NONE . . T2 0.054 0.007
## 31 ppm 76 2.6 GH SEED NONE . . T2 2.600 0.430
## 33 ppm 77 2.6 GH SEED NONE . . T2 0.180 0.045
## 34 ppm 77 2.6 GH SEED NONE . . T3a 8.200 0.460
## SD_AMB CV._AMB N_AMB X_ELEV SE_ELEV SD_ELEV CV._ELEV N_ELEV
## 28 0.0367 38.23 6 0.080 0.016 0.0392 51.04 6
## 29 1.5187 27.28 6 7.000 0.200 0.4899 7.29 6
## 30 0.0171 32.99 6 0.058 0.012 0.0294 52.80 6
## 31 1.0533 42.20 6 5.400 0.580 1.4207 27.41 6
## 33 0.1102 63.77 6 0.190 0.032 0.0784 42.98 6
## 34 1.1268 14.31 6 12.100 0.650 1.5922 13.71 6Calculate size effects
- We will use the log response ratio (LRR = log (Treatment/Control)) as our size effect (‘ROM’ in Metafor)
- Other options include Hedge’s d, Fisher’s r-to-z score, etc.
wt2<-escalc(measure="ROM",m1i=X_ELEV,m2i=X_AMB,sd1i=X_ELEV,sd2i=X_AMB,n1i=N_ELEV,n2i=N_AMB, data=wt,var.names=c("LRR","LRR_var"),digits=4)
gas2<-escalc(measure="ROM",m1i=X_ELEV,m2i=X_AMB,sd1i=X_ELEV,sd2i=X_AMB,n1i=N_ELEV,n2i=N_AMB, data=gas,var.names=c("LRR","LRR_var"),digits=4)Let’s look at these size effects first as a histogram
 |
Then as a forest plot
A forest plot is an easy way to visualise individual effect sizes
Extra options
- The size of each observation can be weighted by the number of observations
- You can plot this after your analysis, and then included the mean ES (and associated 95% CIs)
wt2<-wt2[order(wt2$GENUS),]
forest(wt2$LRR,wt2$LRR_var,slab=wt2$GENUS,pch=19)
Now, we can do the same for the gas exchange parameters
 |
A couple of extra notes on effect sizes in R
- Metafor can directly estimate ES for each study/observation, as well as sample variance
- compute.es (another package) can convert between ES, e.g. from Hedge’s D to Fishers r-to-z, which is really handy if you want to be consistent and use the same ES as part of the same analysis.
How can we test whether CO2 enrichment changes plant biomass?
length(unique(wt2$OBSNO)) #How many unique observations?## [1] 30length(unique(wt2$PAP_NO)) #How many unique studies?## [1] 14Residuals look pretty normal
summary(null)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.8993 27.7986 31.7986 34.5332 32.2601
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0000 0.0000 14 no PAP_NO
##
## Test for Heterogeneity:
## Q(df = 29) = 6.9145, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.3695 0.0938 3.9375 <.0001 0.1856 0.5534 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CO2 has a significant and positive effect on plant biomass
Let’s test for differences in effect sizes between broadloaf and evergreen species
mod.1<-rma.mv(yi=LRR,V=LRR_var, mods=~DIV2,random= ~1|PAP_NO,struct="CS",method="REML",digits=4,data=wt2)
qqnorm(residuals(mod.1,type="pearson"),main="QQ plot: residuals")
qqline(residuals(mod.1,type="pearson"),col="red")
summary(mod.1)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.4911 26.9823 32.9823 36.9789 33.9823
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0000 0.0000 14 no PAP_NO
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 6.5851, p-val = 1.0000
##
## Test of Moderators (coefficient(s) 2):
## QM(df = 1) = 0.3294, p-val = 0.5660
##
## Model Results:
##
## se zval pval ci.lb ci.ub
## intrcpt 0.3909 0.1010 3.8710 0.0001 0.1930 0.5888 ***
## DIV2GYMNO -0.1569 0.2733 -0.5739 0.5660 -0.6926 0.3789
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1How different are angiosperms and gymnosperms in their response to CO2 ?
- Model results show no differences
- We can visualize these differences by plotting estimated size effects for both groups by fitting a model without an intercept and plotting them
out_int<-rma.mv(yi=LRR,V=LRR_var, mods=~DIV2-1,random= ~1|PAP_NO,struct="CS",method="REML",digits=4,data=wt2)
y<-summary(out_int)$b
ci_l<-summary(out_int)$ci.lb
ci_h<-summary(out_int)$ci.ub
fg1<-data.frame(cbind(y,ci_l,ci_h))
colnames(fg1)[1]<-"y"
colnames(fg1)[2]<-"ci_l"
colnames(fg1)[3]<-"ci_h"
fg1$Sperm<-c("Angiosperm","Gymnosperm")
fg1$Sperm<-as.factor(fg1$Sperm)
fg1## y ci_l ci_h Sperm
## DIV2ANGIO 0.3909 0.1930 0.5888 Angiosperm
## DIV2GYMNO 0.2340 -0.2638 0.7318 Gymnosperm
Do photosynthesis and stomatal conductance respond similarly to CO2 enrichment?
We will compare ES of both physiological parameters using a meta-regression model
First, we need to re-arrange the data (using reshape2)
require(reshape2)## Loading required package: reshape2gas3<-dcast(gas2, PAP_NO+GENUS+SPECIES+DIV2~PARAM,value.var="LRR",mean)
gas3<-gas3[!is.nan(gas3[,5]),]
gas3<-gas3[!is.nan(gas3[,6]),]
colnames(gas3)[5]<-"PN_LRR"
colnames(gas3)[6]<-"GS_LRR"
gas4<-dcast(gas2, PAP_NO+GENUS+SPECIES+DIV2~PARAM,value.var="LRR_var",mean)
gas4<-gas4[!is.nan(gas4[,5]),]
gas4<-gas4[!is.nan(gas4[,6]),]
colnames(gas4)[5]<-"PN_var"
colnames(gas4)[6]<-"GS_var"
gass<-merge(gas3,gas4,by.y=c("PAP_NO","GENUS","SPECIES","DIV2"))Fit a model and perform likelihood ratio test
gas.1<-rma.mv(yi=PN_LRR,V=PN_var, mods=~GS_LRR,random= ~1|PAP_NO,struct="CS",method="REML",digits=4,data=gass)
qqnorm(residuals(gas.1,type="pearson"),main="QQ plot: residuals")
qqline(residuals(gas.1,type="pearson"),col="red")
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:
##
## 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 ' ' 1gas1.1<-rma.mv(yi=PN_LRR,V=PN_var, mods=~GS_LRR,random= ~1|PAP_NO,struct="CS",method="ML",digits=4,data=gass)
gas1.2<-rma.mv(yi=PN_LRR,V=PN_var, mods=~1,random= ~1|PAP_NO,struct="CS",method="ML",digits=4,data=gass)
lrt <- as.numeric(2*(logLik(gas1.1) - logLik(gas1.2)))
pvalue<-0.5 * (1 - pchisq(lrt, 1))
pvalue## [1] 0.03961Let’s plot this!
### Use model to make predictions
preds<-predict(gas.1,levels=0, addx=TRUE)
preds<-do.call(cbind.data.frame, preds)
colnames(preds)[8]<-"Gs"
#### fit the line over the real data points
ggplot(preds,aes(x=Gs,y=pred))+
stat_smooth(method="lm",formula=y~x,fullrange=T,se=FALSE,size=1,colour="red")+
geom_point(data=gass,aes(x=GS_LRR,y=PN_LRR),position=position_jitter(width=0.2),pch=19,colour="blue",size=3)+
labs(x="Gs response to CO2 ",y ="Photosynthesis response to CO2") +
theme(axis.title.x=element_text(colour="black",face="bold",size=15,vjust=-0.5),
axis.title.y=element_text(colour="black",face="bold",size=15,vjust=1),
axis.text.y=element_text(colour="black",face="bold",size=12),
axis.text.x=element_text(colour="black",face="bold",size=12),
panel.background =element_rect(fill="transparent",colour="black"),
panel.border=element_rect(fill=NA,colour="black"))
Publication Bias
We can show this graphically using funnel plots, and then test for asymmetry using a modification of Egger’s regression
funnel(gas.1,ylim=c(1:4,by=2),yaxis="seinv",level=c(90, 95, 99),ylab="Precision (1/SE)" ,shade=c("white", "gray", "darkgray"), refline=0)
regtest(gas.1,model="rma",predictor="sei")##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: z = -0.9416, p = 0.3464- The contour funnel plot is fairly symmetric around 0 (good!)
- Egger’s regression test for asymmetry confirms the lack of a publication bias
- Check out our lab website for fascinating research on aboveground-belowground interaction ecology: http://www.eisilab.uni-jena.de/
- We are located now @iDiv (Leipzig, Germany): http://www.idiv.de