2020-08-30-Getting Started with Meta-Analysis in R

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.

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

setwd("/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 29

head(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      3

gas<-subset(dat1,PARAM=="PN" | PARAM=="GS")  ## 
gas$PARAM<-droplevels(gas$PARAM)
dim(gas)

## [1] 80 29

head(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      6

Calculate 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

|
image.png

|

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)

image.png

Now, we can do the same for the gas exchange parameters

|
image.png

|

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] 30

length(unique(wt2$PAP_NO))  #How many unique studies?

## [1] 14
Residuals 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 ' ' 1
CO2 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")
image.png
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 ' ' 1
How 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

image.png

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

gas3<-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")
image.png
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 ' ' 1

gas1.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.03961

Let’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"))
image.png

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)
image.png
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

文章出处：https://rpubs.com/dylanjcraven/metaforr。我只是个搬砖工

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 199,711评论 5赞 468
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 83,932评论 2赞 376
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 146,770评论 0赞 330
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 53,799评论 1赞 271
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 62,697评论 5赞 359
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,069评论 1赞 276
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,535评论 3赞 390
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,200评论 0赞 254
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,353评论 1赞 294
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,290评论 2赞 317
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,331评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,020评论 3赞 315
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,610评论 3赞 303
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,694评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,927评论 1赞 255
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,330评论 2赞 346
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 41,904评论 2赞 341