Motivation

Today we examine methods that summarize quantitative data from across published studies. From each published study an effect size and weight are calculated. These are then summarized using GLM, where the weights are the inverse of study variances.

First let’s read in some data:

library(metafor)
## Loading required package: Matrix
## Loading 'metafor' package (version 2.4-0). For an overview 
## and introduction to the package please type: help(metafor).
library(RRPP)
mydata<-read.csv("Data/Lab-14-gur_hedge.csv",header=T)
habitat<-mydata[,(1)]
effectd<-as.matrix(mydata[,(2)])                                                             
var<-as.matrix(mydata[,(4)])

1: Fixed Effects Meta-Analysis

Here is the simplest meta-analytic model: a fixed effects m-a with no structure (i.e., all studies belong to a single group)

No Structure Fixed Effects M-A

#no structure
ma.no<-rma.uni(yi=effectd,vi=var,data=mydata,method="FE")
summary(ma.no)
## 
## Fixed-Effects Model (k = 43)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -62.7549   85.9775  127.5097  129.2709  127.6073   
## 
## I^2 (total heterogeneity / total variability):   51.15%
## H^2 (total variability / sampling variability):  2.05
## 
## Test for Heterogeneity:
## Q(df = 42) = 85.9775, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval   ci.lb   ci.ub 
##   1.0099  0.0838  12.0536  <.0001  0.8457  1.1741  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(ma.no)  #plot of study effect sizes

Categorical Structure Fixed Effects M-A

The next model is equivalent to meta-analytic anova:

ma.cat<-rma.uni(yi=effectd,vi=var,mods= ~habitat, data=mydata,method="FE")
summary(ma.cat)
## 
## Fixed-Effects with Moderators Model (k = 43)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -54.5152   69.4982  115.0304  120.3140  115.6458   
## 
## I^2 (residual heterogeneity / unaccounted variability): 42.44%
## H^2 (unaccounted variability / sampling variability):   1.74
## 
## Test for Residual Heterogeneity:
## QE(df = 40) = 69.4982, p-val = 0.0026
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 16.4793, p-val = 0.0003
## 
## Model Results:
## 
##                     estimate      se     zval    pval    ci.lb    ci.ub 
## intrcpt               4.1072  0.8857   4.6373  <.0001   2.3713   5.8431  *** 
## habitatMarine        -3.3087  0.8942  -3.7001  0.0002  -5.0614  -1.5561  *** 
## habitatTerrestrial   -2.9655  0.8931  -3.3203  0.0009  -4.7160  -1.2150  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Group means
ma.cat$b[1] #lentic
## [1] 4.107203
ma.cat$b[1]+ma.cat$b[2] #marine
## [1] 0.7984534
ma.cat$b[1]+ma.cat$b[3] #terrestrial
## [1] 1.14171

2: Random Effects Meta-Analysis

These models incorporate a random effect to the study variance, much like random effects in a standard anova.

No Structure Random Effects M-A

rma.uni(yi=effectd,vi=var,data=mydata)
## 
## Random-Effects Model (k = 43; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.3013 (SE = 0.1410)
## tau (square root of estimated tau^2 value):      0.5489
## I^2 (total heterogeneity / total variability):   49.60%
## H^2 (total variability / sampling variability):  1.98
## 
## Test for Heterogeneity:
## Q(df = 42) = 85.9775, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.9346  0.1266  7.3804  <.0001  0.6864  1.1828  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Categorical Structure Random Effects M-A

This model is equivalent to a random effects m-a anova:

ma.catr<-rma.uni(yi=effectd,vi=var,mods= ~habitat, data=mydata)
summary(ma.catr)
## 
## Mixed-Effects Model (k = 43; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -48.1754   96.3508  104.3508  111.1063  105.4937   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.2438 (SE = 0.1284)
## tau (square root of estimated tau^2 value):             0.4937
## I^2 (residual heterogeneity / unaccounted variability): 44.44%
## H^2 (unaccounted variability / sampling variability):   1.80
## R^2 (amount of heterogeneity accounted for):            19.09%
## 
## Test for Residual Heterogeneity:
## QE(df = 40) = 69.4982, p-val = 0.0026
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 13.7892, p-val = 0.0010
## 
## Model Results:
## 
##                     estimate      se     zval    pval    ci.lb    ci.ub 
## intrcpt               4.1174  0.9534   4.3186  <.0001   2.2488   5.9861  *** 
## habitatMarine        -3.4205  0.9687  -3.5310  0.0004  -5.3191  -1.5219  *** 
## habitatTerrestrial   -3.0368  0.9688  -3.1347  0.0017  -4.9355  -1.1381   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Group means
ma.catr$b[1] #lentic
## [1] 4.117416
ma.catr$b[1]+ma.catr$b[2] #marine
## [1] 0.6968947
ma.catr$b[1]+ma.catr$b[3] #terrestrial
## [1] 1.080635

3: Additional Analyses

As we discussed, there are many additional ways of exploring meta-analytic data: phylogenetic meta-analysis, cumulative meta-analysis, funnel plots, fail safe numbers, etc. Many of these may be found in the ‘meta’, ‘rmeta’ and ‘metafor’ package. Below are several examples.

#Funnel plot for outliers
funnel(ma.no)

# fail safe number
fsn(yi=effectd,vi=var,data=mydata)
## 
## Fail-safe N Calculation Using the Rosenthal Approach
## 
## Observed Significance Level: <.0001
## Target Significance Level:   0.05
## 
## Fail-safe N: 1775
fsn(yi=effectd,vi=var,data=mydata,type="Rosenberg")
## 
## Fail-safe N Calculation Using the Rosenberg Approach
## 
## Average Effect Size:        1.0099
## Observed Significance Level: <.0001
## Target Significance Level:   0.05
## 
## Fail-safe N: 1584

4: Using GLS ‘by hand’ in RRPP

Recognizing that meta-analysis is simply a GLS model where the weight matrix is a diagonal matrix containing study weights, one can use RRPP to perform the same task:

rdf <- rrpp.data.frame(effectd = effectd, habitat=habitat, var.e = diag(as.vector(var)))

anova(lm.rrpp(effectd~1,Cov = rdf$var.e, data = rdf, print.progress = FALSE))
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Generalized Least-Squares (via OLS projection) 
## No model effects; simple summary provided
## 
##           df     SS     MS
## Residuals 42 85.978 2.0471
summary(ma.no)  #same SS!
## 
## Fixed-Effects Model (k = 43)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -62.7549   85.9775  127.5097  129.2709  127.6073   
## 
## I^2 (total heterogeneity / total variability):   51.15%
## H^2 (total variability / sampling variability):  2.05
## 
## Test for Heterogeneity:
## Q(df = 42) = 85.9775, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval   ci.lb   ci.ub 
##   1.0099  0.0838  12.0536  <.0001  0.8457  1.1741  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res<- lm.rrpp(effectd~habitat,Cov = rdf$var.e, data = rdf, print.progress = FALSE)
anova(res)
## 
## Analysis of Variance, using Residual Randomization
## Permutation procedure: Randomization of null model residuals 
## Number of permutations: 1000 
## Estimation method: Generalized Least-Squares (via OLS projection) 
## Sums of Squares and Cross-products: Type I 
## Effect sizes (Z) based on F distributions
## 
##           Df     SS     MS     Rsq      F      Z Pr(>F)  
## habitat    2 16.479 8.2397 0.19167 4.7424 2.1782  0.014 *
## Residuals 40 69.498 1.7375 0.80833                       
## Total     42 85.978                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call: lm.rrpp(f1 = effectd ~ habitat, data = rdf, Cov = rdf$var.e,  
##     print.progress = FALSE)
summary(ma.cat) #same Q tests!
## 
## Fixed-Effects with Moderators Model (k = 43)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -54.5152   69.4982  115.0304  120.3140  115.6458   
## 
## I^2 (residual heterogeneity / unaccounted variability): 42.44%
## H^2 (unaccounted variability / sampling variability):   1.74
## 
## Test for Residual Heterogeneity:
## QE(df = 40) = 69.4982, p-val = 0.0026
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 16.4793, p-val = 0.0003
## 
## Model Results:
## 
##                     estimate      se     zval    pval    ci.lb    ci.ub 
## intrcpt               4.1072  0.8857   4.6373  <.0001   2.3713   5.8431  *** 
## habitatMarine        -3.3087  0.8942  -3.7001  0.0002  -5.0614  -1.5561  *** 
## habitatTerrestrial   -2.9655  0.8931  -3.3203  0.0009  -4.7160  -1.2150  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$LM$gls.coefficients #same group effect sizes
##        (Intercept)      habitatMarine habitatTerrestrial 
##           4.107203          -3.308749          -2.965493

Note that for both models, one obtais identical results. This confirms that conceptually, meta-analysis is simply a GLS model containing a weight matrix which is a diagonal matrix with study weights on the diagonal.