30 Meta Analysis
Meta-analyses are increasingly used in ecology, evolution and the environmental sciences to find general patterns among many studies, settle controversies among conflicting studies and to generate new hypotheses. These tutorials provide an introduction to running meta-analyses with the R package metafor.
-
Meta-analyses 1: Introduction and calculating effect sizes
-
Meta-analyses 2: Fixed effect and random effect models
- Meta-analyses 3: More complex models
30.1 Introduction
30.1.1 Background
What is a meta-analysis?
A meta-analysis is a quantitative summary of studies on the same topic.
Why do we want to perform a meta-analysis?
- Finding generalities
- Increasing power and precision
- Exploring differences between studies
- Settling controversies from conflicting studies (testing hypotheses)
- Generating new hypotheses
The process of meta-analysis
How many steps involved in meta-analysis?
One answer is 5 steps
- Formulating questions & hypothesis or finding a topic
- Literature search & paper collection
- Data extraction & coding
- Meta-analysis & publication bias tests
- Reporting & publication
We only consider the step iv in this tutorial. You will need to learn the other steps elsewhere. To get you started, we recently wrote an overview paper which divides the process of meta-analysis into 10 questions (Nakagawa et al. 2017). The 10 questions will guide you through judging the quality of a meta-analysis.
- Is the search systematic and transparently documented?
- What question and what effect size?
- Is non-independence taken into account?
- Which meta-analytic model?
- Is the level of consistency among studies reported?
- Are the causes of variation among studies investigated?
- Are effects interpreted in terms of biological importance?
- Has publication bias been considered?
- Are results really robust and unbiased?
- Is the current state (and lack) of knowledge summarized?
30.1.2 Metafor for meta-analysis
I think the R package metafor (Viechtbauer 2010) is the most comprehensive meta-analytic software and the author Wolfgang Viechtbauer, who, I have to say, has the coolest name among my friends, is still actively developing it.
First, install and load the metafor package.
Have a look at the data set named dat.curtis1998, from the metadat package that comes with metafor. To see the other data sets included in that package, try help(package="metadat").
dat <- metadat::dat.curtis1998
str(dat)'data.frame': 102 obs. of 20 variables:
$ id : int 21 22 27 32 35 38 44 63 86 87 ...
$ paper : int 44 44 121 121 121 121 159 183 209 209 ...
$ genus : chr "ALNUS" "ALNUS" "ACER" "QUERCUS" ...
$ species : chr "RUBRA" "RUBRA" "RUBRUM" "PRINUS" ...
$ fungrp : chr "N2FIX" "N2FIX" "ANGIO" "ANGIO" ...
$ co2.ambi: num 350 350 350 350 350 350 350 395 350 350 ...
$ co2.elev: num 650 650 700 700 700 700 700 795 700 700 ...
$ units : chr "ul/l" "ul/l" "ppm" "ppm" ...
$ time : int 47 47 59 70 64 50 730 365 365 365 ...
$ pot : chr "0.5" "0.5" "2.6" "2.6" ...
$ method : chr "GC" "GC" "GH" "GH" ...
$ stock : chr "SEED" "SEED" "SEED" "SEED" ...
$ xtrt : chr "FERT" "FERT" "NONE" "NONE" ...
$ level : chr "HIGH" "CONTROL" "." "." ...
$ m1i : num 6.82 2.6 2.99 5.91 4.61 ...
$ sd1i : num 1.77 0.667 0.856 1.742 1.407 ...
$ n1i : int 3 5 5 5 4 5 3 3 20 16 ...
$ m2i : num 3.94 2.25 1.93 6.62 4.1 ...
$ sd2i : num 1.116 0.328 0.552 1.631 1.257 ...
$ n2i : int 5 5 5 5 4 3 3 3 20 16 ...This data set is from the paper by Curtis and Wang (1998). They looked at the effect of increased CO\(_2\) on plant traits (mainly changes in biomass). So we have information on control group (1) and experimental group (2) (m = mean, sd = standard deviation) along with species information and experimental details. In meta-analysis, these variables are often referred to as ‘moderators’ (we will get to this a bit later).
30.1.3 Calculating ‘standardized’ effect sizes
To compare the effect of increased CO\(_2\) across multiple studies, we first need to calculate an effect size for each study - a metric that quantifies the difference between our control and experimental groups.
There are several ‘standardized’ effect sizes, which are unit-less. When we have two groups to compare, we use two types of effect size statistics. The first one is standardized mean difference (SMD or also known as Cohen’s \(d\) or Hedge’s \(d\) or \(g\); some subtle differences between them, but we do not worry about them for now):
\[\begin{equation} \mathrm{SMD}=\frac{\bar{x}_{E}-\bar{x}_{C}}{\sqrt{\frac{(n_{C}-1)sd^2_{C}+(n_{E}-1)sd^2_{E}}{n_{C}+n_{E}-2}}} \end{equation}\] where \(\bar{x}_{C}\) and \(\bar{x}_{E}\) are the means of the control and experimental group, respectively, \(sd\) is sample standard deviation (\(sd^2\) is sample variance) and \(n\) is sample size.
And its sample error variance is:
\[\begin{equation} se^2_{\mathrm{SMD}}= \frac{n_{C}+n_{E}}{n_{C}n_{E}}+\frac{\mathrm{SMD}^2}{2(n_{C}+n_{E})} \end{equation}\]
The square root of this is referred to as ‘standard error’ (or the standard deviation of the estimate – confused?). The inverse of this (\(1/se^2\)) is used as ‘weight’, but things are bit more complicated than this as we will find out below.
Another common index is called ‘response ratio’, which is usually presented in its natural logarithm form (lnRR):
\[\begin{equation} \mathrm{lnRR}=\ln\left({\frac{\bar{x}_{E}}{\bar{x}_{C}}}\right) \end{equation}\]
And its sampling error variance is:
\[\begin{equation} se^2_\mathrm{lnRR}=\frac{sd^2_{C}}{n_{C}\bar{x}^2_{C}}+\frac{sd^2_{E}}{n_{E}\bar{x}^2_{E}} \end{equation}\]
Let’s get these using the function escalc in metafor. To obtain the standardised mean difference, we use:
# SMD
SMD <- escalc(measure = "SMD", n1i = dat$n1i, n2i = dat$n2i, m1i = dat$m1i, m2i = dat$m2i,
sd1i = dat$sd1i, sd2i = dat$sd2i)where n1i and n2i are the sample sizes, m1i and m2i are the means, and sd1i and sd2i the standard deviations from each study.
The object created now has an effect size (yi) and its variance (vi) for each study
yi vi
1 1.8222 0.7408
2 0.5922 0.4175
3 1.3286 0.4883
4 -0.3798 0.4072
5 0.3321 0.5069
6 2.5137 0.9282 To obtain the response ratio (log transformed ratio of means), we would use:
lnRR <- escalc(measure = "ROM", n1i = dat$n1i, n2i = dat$n2i, m1i = dat$m1i, m2 = dat$m2i,
sd1i = dat$sd1i, sd2i = dat$sd2i)The original paper used lnRR so we will use it, but you may want to repeat analysis below using SMD to see whether results are consistent.
Add the effect sizes to the original data set with cbind or bind_cols from the package dplyr
You should see yi (effect size) and vi (sampling variance) are added.
'data.frame': 102 obs. of 22 variables:
$ id : int 21 22 27 32 35 38 44 63 86 87 ...
$ paper : int 44 44 121 121 121 121 159 183 209 209 ...
$ genus : chr "ALNUS" "ALNUS" "ACER" "QUERCUS" ...
$ species : chr "RUBRA" "RUBRA" "RUBRUM" "PRINUS" ...
$ fungrp : chr "N2FIX" "N2FIX" "ANGIO" "ANGIO" ...
$ co2.ambi: num 350 350 350 350 350 350 350 395 350 350 ...
$ co2.elev: num 650 650 700 700 700 700 700 795 700 700 ...
$ units : chr "ul/l" "ul/l" "ppm" "ppm" ...
$ time : int 47 47 59 70 64 50 730 365 365 365 ...
$ pot : chr "0.5" "0.5" "2.6" "2.6" ...
$ method : chr "GC" "GC" "GH" "GH" ...
$ stock : chr "SEED" "SEED" "SEED" "SEED" ...
$ xtrt : chr "FERT" "FERT" "NONE" "NONE" ...
$ level : chr "HIGH" "CONTROL" "." "." ...
$ m1i : num 6.82 2.6 2.99 5.91 4.61 ...
$ sd1i : num 1.77 0.667 0.856 1.742 1.407 ...
$ n1i : int 3 5 5 5 4 5 3 3 20 16 ...
$ m2i : num 3.94 2.25 1.93 6.62 4.1 ...
$ sd2i : num 1.116 0.328 0.552 1.631 1.257 ...
$ n2i : int 5 5 5 5 4 3 3 3 20 16 ...
$ yi : num 0.55 0.147 0.438 -0.111 0.117 ...
..- attr(*, "ni")= int [1:102] 8 10 10 10 8 8 6 6 40 32 ...
..- attr(*, "measure")= chr "ROM"
$ vi : num 0.0385 0.0175 0.0328 0.0295 0.0468 ...Visualising effect size. We can visualize point estimates (effect size) and their 95% confidence intervals, CIs (based on sampling error variance) by using the forest function, which draws a forest plot for us.
forest(dat$yi, dat$vi)
The problem you see is that when there are many studies, a forest plot does not really work (unless you have very large screen!). Let’s look at just the first 12 studies.
forest(dat$yi[1:12], dat$vi[1:12])
We can calculate many different kinds of effect sizes with escalc; other common effect size statistics include \(Zr\) (Fisher’s z-transformed correlation). By the way, along with my colleagues, I have proposed a new standardized effect size called lnCVR (the log of coefficient of variation ratio – mouthful!), which compares the variability of two groups rather than means. See whether you can calculate it with these data. Actually, the development version of metafor, let you do this with escalc– github page. lnCVR is called “CVR” in escalc. Actually, if you re-analysis this data with lnCVR, you may be able to publish a paper! Nobody has done it yet. Do it tonight!
Once you have calculated effect sizes, move on to the next page: Meta-analysis 2: fixed-effect and random-effect models
30.1.4 Further help (references)
Go to the metafor package’s website. There you find many worked examples.
Curtis, P. S., and X. Z. Wang. 1998. A meta-analysis of elevated CO2 effects on woody plant mass, form, and physiology. Oecologia 113:299-313.
Nakagawa, S., R. Poulin, K. Mengersen, K. Reinhold, L. Engqvist, M. Lagisz, and A. M. Senior. 2015. Meta-analysis of variation: ecological and evolutionary applications and beyond. Methods in Ecology and Evolution 6:143-152.
Viechtbauer, W. 2010. Conducting meta-analyses in R with the metafor package. Journal of Statistical Software 36:1-48.
Authors: Shinichi Nakagawa and Malgorzata (Losia) Lagisz
Year: 2016
30.2 Fixed-effect and Random-effect Models
First read our introduction to meta analysis with instructions for calculating effect sizes. Now, we consider the statistical models to estimate mean effect sizes and the influence of moderator variables.
30.2.1 Fixed-effect and random-effects models
There are two main models of meta-analysis: 1) the fixed-effect model and 2) the random-effect model (actually we say 3 types, but the third one is the extension of the second model). Because fixed effects mean something different in another context, this naming is a bit confusing. So people now call the first model ‘common-effect’ model. You can see a visual representation of these two models (from Figure 4 in Nakagawa et al. 2017).
A math representation for the common-effect model is:
\[\begin{equation} y_i=b_0+e_i, \\ e_i\sim \mathcal{N}(0,v_i),\ \end{equation}\]\end{equation}
where \(y_i\) is the \(i\)th effect size (from \(i\)th study), \(b_0\) is the overall mean (or meta-analytic mean), \(e_i\) is a deviation from the overall mean and it is normally distributed with \(v_i\), which is the study specific sampling variance (you can get this from formulas above). Note that weights for this model are \(1/v_i\). As you can see, this is a very simple model; basically, we are estimating the mean considering weights.
Importantly, this model is assuming that there is a common mean for all studies included. Is this realistic? Probably not, especially, if you are combining different species!
The random-effects model can be written as:
\[\begin{equation}
y_i=b_0+s_i+e_i,
\\
s_i\sim \mathcal{N}(0,\tau^2),\
\\
e_i\sim \mathcal{N}(0,v_i),\
\end{equation}\]e_i(0,v_i),
\end{equation}
where \(s_i\) is a study-specific deviation from the overall mean for \(i\)th study, it is normally distributed with the between-study variance, often referred to as \(\tau^2\), and the others are the same as above. Unlike the common-effect model, a random-effect model assumes that different studies have different means. Note that weights for this model are \(1/(\tau^2+v_i)\). We revisit this point, as it turns out to be quite important.
30.2.2 Running a common-effect model
Let’s use the function rma from metafor to run a common-effect model.
common_m <- rma(yi = yi, vi = vi, method = "FE", data = dat)We specify the effect size (yi), its variance (vi), the method (“FE” for fixed-effect) and the data frame (dat).
To see the output, use summary on the model object:
summary(common_m)
Fixed-Effects Model (k = 102)
logLik deviance AIC BIC AICc
-244.4167 765.9358 490.8334 493.4583 490.8734
I^2 (total heterogeneity / total variability): 86.81%
H^2 (total variability / sampling variability): 7.58
Test for Heterogeneity:
Q(df = 101) = 765.9358, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2087 0.0054 38.3072 <.0001 0.1980 0.2193 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1OK, it was easy! The overall mean is statistically significant and it’s around 0.2. What does 0.2 mean? Let’s convert this back into the original scale - a response ratio of the control and experimental means.
exp(0.2)[1] 1.221403This means a plant trait (e.g. mass) was 22% larger in the experimental group (RR\(=\bar{x}_{E}/\bar{x}_{C}\)), which seems like a pretty large effect to me (remember we need to interpret our results in a biological meaningful way). Note that Jensen’s inequality states that this is a bit wrong, but we will not get into this today.
30.2.3 Running a random-effects model
Now, we move onto the random-effects model - a more realistic model. Again, we use the rma function, but this time change the method to REML which is the default and the best method for the random-effect meta-analysis.
Random-Effects Model (k = 102; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
8.0220 -16.0439 -12.0439 -6.8137 -11.9215
tau^2 (estimated amount of total heterogeneity): 0.0258 (SE = 0.0052)
tau (square root of estimated tau^2 value): 0.1605
I^2 (total heterogeneity / total variability): 88.73%
H^2 (total variability / sampling variability): 8.87
Test for Heterogeneity:
Q(df = 101) = 765.9358, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2540 0.0197 12.9017 <.0001 0.2154 0.2925 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Compare the overall mean from this model with the common-effect model. Oh, the overall mean of the random-effects model is actually bigger than that of the fixed-effect model. OK, that sometimes happens (we will find out that this probably is an over-estimation later). We expect 95% CI is wider (i.e. more realistic) in this random-effects model as this model has a better assumption than the common-effect model.
30.2.4 Understanding heterogeneity
The output from the random-effects model has more things than the common-effect model. We have tau^2 (\(\tau^2\)) and I^2 (\(I^2\)), two very common measures of heterogeneity (note that H^2, or \(H^2\) is a transformation of \(I^2\)).
Heterogeneity is variation in effect sizes, which is not accounted by the sampling error variance.
In other words, real variation in the data. I think \(I^2\) is a quite important index as it can tell the percentage of real variation in your meta-analytic data. It can be given by (which you saw within the figure of the random-effects model above):
\[\begin{equation} I^2=\frac{\tau^2}{(\tau^2+\bar{v})}, \end{equation}\]
where \(\bar{v}\) is a representative value of \(v_i\) (or think \(\bar{v}\) as the average of \(v_i\) although it is not quite it). Note that the denominator is the whole variance which exists in the data. The benchmark values for \(I^2\) are 25, 50 and 75% for low, moderate and high heterogeneity, respectively (Higgins et al. 2003).
Our \(I^2\) value is 88.9% so very high and this value is, as you expect, statistically significant when tested with the \(Q\) value (which follows a \(\chi^2\) distribution defined by the df value, in this case \(df = 101\)).
We recently did a meta-analysis of meta-analyses (a secondary meta-analysis) looking at what is the average value of \(I^2\) in the field of ecology and evolution (Senior et al. 2016). The average value was 92%! So this indicates that we should really be fitting the random-effects model rather than the common-effect model because the latter assumes heterogeneity to be zero or \(\tau^2=0\) and \(I^2 = 0\). Or is it really? We find this out later.
30.2.5 Meta-regression (the random-effects model)
The existence of heterogeneity sets a scene for meta-regression. This means that we now put predictors (‘moderators’ in the meta-analytic terminology) into our model to explain heterogeneity (equivalent to normal regression models).
Fitting meta-regression is pretty straightforward too. In this example, let’s fit three moderators that were collected by the authors: 1) time (how long the experiment was), 2) method (different ways of increasing CO\(_2\)), and 3) fungroup (functional group, i.e., angiosperm, gymnosperm or N\(_2\) fixer).
We use rma again, but now have a model statement.
Mixed-Effects Model (k = 102; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
6.2121 -12.4241 1.5759 19.5263 2.8486
tau^2 (estimated amount of residual heterogeneity): 0.0262 (SE = 0.0055)
tau (square root of estimated tau^2 value): 0.1618
I^2 (residual heterogeneity / unaccounted variability): 86.93%
H^2 (unaccounted variability / sampling variability): 7.65
R^2 (amount of heterogeneity accounted for): 0.00%
Test for Residual Heterogeneity:
QE(df = 96) = 655.1885, p-val < .0001
Test of Moderators (coefficients 2:6):
QM(df = 5) = 3.0606, p-val = 0.6906
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3033 0.0512 5.9216 <.0001 0.2029 0.4037 ***
time -0.0001 0.0001 -1.0830 0.2788 -0.0002 0.0000
methodGH -0.0374 0.0563 -0.6645 0.5064 -0.1478 0.0729
methodOTC 0.0311 0.0895 0.3471 0.7285 -0.1444 0.2066
fungrpGYMNO -0.0443 0.0600 -0.7388 0.4600 -0.1620 0.0733
fungrpN2FIX 0.0090 0.1692 0.0534 0.9574 -0.3226 0.3406
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Well, they do not explain anything! Look at the \(R^2\) value and, as you expect, the test for moderators (again the Q value) say they are not significant. A terrible model! So we give up here (in a real meta-analysis, you need to do this more systematically, preferably based on your priori hypotheses).
30.2.6 Checking for publication bias
OK, it seems like an CO\(_2\) increase promotes plant growth (which may not be surprising), but we are assuming the data set we have does not suffer from publication bias.
Publication bias in its simplest form is that significant results are more likely to be published than non-significant results.
But there are several methods, which people have been using. The two commonest methods, often used as a set, are: 1) funnel plot, which one uses to detect a funnel asymmetry (a sign of publication bias), and 2) Egger’s regression test with which you test funnel asymmetry statistically.
What am I talking about by ‘funnel asymmetry’? If you plot effect sizes and its standard error in a decreasing order (\(se\); see the figure below), we suppose to see a funnel up-side-down with effect sizes with low \(se\) (or high weights) are more tightly aggregated than effect sizes with high \(se\) (or low weights). But this funnel shape occurs only when we do not have publication bias. If there is, we should see funnel asymmetry. This is because studies with small sample sizes (i.e. high \(se\), which leads to non-significance) are less likely to be published. We will see an example of this funnel asymmetry in our data set.
To create a funnel plot and run Egger’s test:
funnel(random_m)
# Note that the orignal Egger's test is regtest(random_m, model='lm')
regtest(random_m)
Regression Test for Funnel Plot Asymmetry
Model: mixed-effects meta-regression model
Predictor: standard error
Test for Funnel Plot Asymmetry: z = 3.0911, p = 0.0020
Limit Estimate (as sei -> 0): b = 0.1610 (CI: 0.0918, 0.2301)The funnel plot and Egger’s test both suggest funnel asymmetry. But we need to be careful. Funnel asymmetry cannot be caused not only by publication bias but also by heterogeneity (one or more undetected moderators are distorting a funnel shape). Given we have a lot of unexplained variance (i.e. heterogeneity), we are not sure which is causing this asymmetry. But we have some evidence for publication bias. A relating point is that if meta-regression explains a lot of heterogeneity, you should use that meta-regression model in the function regrest (in our case regrest(metareg), but I did not do it as our model metareg did not explain anything).
There is an alternative method, which has a cool name, the trim-and-fill method. You will see the reason why. We can use this method by using the function trimfill.
# Note that we are using the default estimator ('L0'), but there are two others
# available
tf_m <- trimfill(random_m)
tf_m
Estimated number of missing studies on the left side: 15 (SE = 6.6035)
Random-Effects Model (k = 117; tau^2 estimator: REML)
tau^2 (estimated amount of total heterogeneity): 0.0435 (SE = 0.0077)
tau (square root of estimated tau^2 value): 0.2085
I^2 (total heterogeneity / total variability): 92.24%
H^2 (total variability / sampling variability): 12.88
Test for Heterogeneity:
Q(df = 116) = 899.7193, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2084 0.0228 9.1479 <.0001 0.1637 0.2530 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1funnel(tf_m)
As you can see this method uses the asymmetry to add more points and provide a revised overall mean, which is smaller than that of the original random-effect model. Although this effect is still significant, this method could turn a significant overall mean into a non-significant one. But rather than taking this is a real estimate of the overall mean, we need to see this as a part of sensitivity analysis.
There are more methods for publication bias tests, none of which are perfect, but we do need to do some of these tests (see more for Nakagawa et al. 2017 and references therein).
30.2.7 Another model (I did not know), but important!
I recently learned that there is another meta-analytic model (neither common nor random), which is more robust to publication bias. Remember what I said above.
- The random-effect model is a more realistic model so we should use this model.
- The random-effect model has the weight of \(1/(\tau^2+v_i)\).
- \(I^2\) is very large in ecological and evolutionary meta-data sets so this means \(\tau^2\) is very large, meaning \(1/(\tau^2+v_i) \approx 1/\tau^2\) as \(v_i\) is negligible
- Then, all data points having the same weight \(1/\tau^2\) and basically it becomes an unweighted model!
This is problematic because under publication bias like we saw it just above, effect sizes with small sample sizes (very high \(v_i\)) should get small weights compared to effect sizes with large sample sizes (low \(v_i\)). When there exists publication bias, the random-effect model could provide us with a biased estimate… What could we do? Can we use the weight of \(1/v_i\) like the common-effect model, but use the structure of the random-effects model? It turns out you can!
Random-Effects Model (k = 102; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
5.3743 -10.7486 -6.7486 -1.5184 -6.6262
tau^2 (estimated amount of total heterogeneity): 0.0258 (SE = 0.0052)
tau (square root of estimated tau^2 value): 0.1605
I^2 (total heterogeneity / total variability): 88.73%
H^2 (total variability / sampling variability): 8.87
Test for Heterogeneity:
Q(df = 101) = 765.9358, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2087 0.0469 4.4451 <.0001 0.1167 0.3007 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Actually, this point estimate is the exactly the same as that of the common-effect model (this is not just a coincidence as this model calculates the overall mean like the common-effect model). But it is important to notice this model has a wider 97% CI than the common-effect model and this level of uncertainty is more realistic (just like the random-effects model). Basically, this model has got the best of both worlds!
This model is based on the paper by Henmi and Copas (2010). Although it is more robust to publication bias, I do not know any meta-analyses in the field of ecology and evolution, which uses this model. Probably we should be using this model rather than a random-effect model? I think so - let’s call this model the ‘robust’ model. I think we should use this robust model for meta-regression too. I acknowledge Wolfgang telling me about this model to me.
If you are happy with all this, move on the more complex models
30.2.8 Further help (references)
Any questions? Or email me (s(-dot-)nakagawa(-at-)unsw(-dot-)edu(-dot-)au). Also visit our website
Go to the metafor package’s website. There you find many worked examples.
Henmi, M., and J. B. Copas. 2010. Confidence intervals for random effects meta-analysis and robustness to publication bias. Statistics in Medicine 29:2969-2983.
Higgins, J. P. T., S. G. Thompson, J. J. Deeks, and D. G. Altman. 2003. Measuring inconsistency in meta-analyses. British Medical Journal 327:557-560.
Nakagawa, S., D. W. A. Noble, A. M. Senior, and M. Lagisz. 2017. Meta-evaluation of meta-analysis: ten appraisal questions for biologists. BMC Biology 15:18.
Senior, A. M., C. E. Grueber, T. Kamiya, M. Lagisz, K. O’Dwyer, E. S. A. Santos, and S. Nakagawa. 2016. Heterogeneity in ecological and evolutionary meta-analyses: its magnitude and implications. Ecology 97:3293-3299.
Authors: Shinichi Nakagawa and Malgorzata (Losia) Lagisz
Year: 2016
30.3 Multilevel Models
First read our introduction to meta analysis and become familiar with commonly used statistical models for meta-analysis. Now, we consider more complex models.
30.3.1 More complex models
Up to now, we have completely ignored non-independence among effect sizes. We assumed that we have one effect size from one study (or paper). But in reality, a paper usually contains multiple effect sizes. These effect sizes from the same studies are not independent of each other. Let’s consider a model where we have several effect sizes from single studies like our data set. First, I give you a math representation:
\[\begin{equation}
y_{ij}=b_0+s_i+u_{ij}+e_{ij}
\\
s_i\sim \mathcal{N}(0,\tau^2)\
\\
u_{ij}\sim \mathcal{N}(0,\sigma^2)\
\\
e_{ij}\sim \mathcal{N}(0,v_{ij})\
\end{equation}\]\ e_{ij}(0,v_{ij})
\end{equation}
where \(u_{ij}\) is a deviation from \(s_i\) (the within-study effect; the \(j\)th effect size from the \(i\)th study), it is normally distributed with \(\sigma^2\) (other notations are comparable as above).
We can visualize this (again Figure 4 from Nakagawa et al. 2017). And we can see why this model is called, a ‘multilevel’ meta-analytic model, an extension of the random-effects model.
We can fit an equivalent model using the function rma.mv. We need to add paper (the between-study effect) and id (different effect sizes; the within-study effect) in our data set to the multilevel model as random factors.
Multivariate Meta-Analysis Model (k = 102; method: REML)
logLik Deviance AIC BIC AICc
8.0845 -16.1690 -10.1690 -2.3237 -9.9216
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0014 0.0377 29 no paper
sigma^2.2 0.0244 0.1562 102 no id
Test for Heterogeneity:
Q(df = 101) = 765.9358, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2564 0.0219 11.6818 <.0001 0.2134 0.2994 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1OK, this does not have \(I^2\). We have actually proposed a multilevel-model version of \(I^2\) (Nakagawa & Santos 2012), which, in this case, can be written as:
\[\begin{equation} I^2=\frac{\tau^2+\sigma^2}{(\tau^2+\sigma^2+\bar{v})}, \end{equation}\]
Note that we do have \(\tau^2\) and \(\sigma^2\), which are sigma^2.1 and sigma^2.2 in the output above, respectively. Using this formula, we have the total heterogeneity \(I^2\) of 88.93% (the \(\bar{v} = 0.0033\) for our data set; see Nakagawa & Santos 2012 for how to get this). As you might expect, this value is nearly identical to what we got from the random-effect model (88.9%). But this model is better as we are explicitly dealing with non-independence arising from effect sizes from the same studies (although it turns out the problem is not completely solved…).
As you could probably imagine, we can add more levels to this multilevel models. For example, we could add genus in the data set, as related species are probably more similar to each other. But it is better to model this taxonomic non-independence using phylogeny (which is the topic of the next section). Note that we can also run meta-regression using rma.mv; more complex models (different versions of a multilevel model) are explained in Nakagawa & Santos (2012).
30.3.2 Phylogenetic meta-analysis
As I just mentioned, we have, so far, also ignored another type of non-independence, i.e. phylogenetic relatedness. Chuck Darwin unfortunately (or fortunately) found out all species on earth are related so we need to deal with this issue.
Actually, we just need to model phylogenetic non-independence by adding the degree of relatedness among species as a correlation matrix. The term, phylogeny or phylo, which we will create below, can be added as a random factor to a multilevel model.
This means that we need a phylogenetic tree for our data set. For this data set, we have prepared a tree for you to download here. But it is not that difficult to get a tree together by using the package called rotl.
First install and load the package ape that we will use to import the tree file and visualise the phylogeny.
We can then make a correlation matrix (a relatedness matrix among species). I’ve skipped explanations of these operations - other than saying we have a correlation matrix to fit to the model
tree <- compute.brlen(tree)
cor <- vcv(tree, cor = T)We need a bit more preparation as we do have not a column which has the whole species names (we call it phylo). Also, it turns out that we need to correct some typos in genus and species columns in our data.
library(Hmisc)
phylo <- tolower(paste(dat$genus, dat$species, sep = "_"))
# note: 'populusx_euramericana' should be same as 'populus_euramericana'
phylo <- gsub("populusx_euramericana", "populus_euramericana", phylo)
# these two species are the two different names of the same species
phylo <- gsub("nothofagus_fusca", "fuscospora_fusca", phylo)
phylo <- capitalize(phylo)
dat[, "phylo"] <- phyloWe now have our phylogenetic correlation cor and a column with species names phylo, and can run our meta-analysis again with a phylogenetic effect.
Multivariate Meta-Analysis Model (k = 102; method: REML)
logLik Deviance AIC BIC AICc
8.0845 -16.1690 -8.1690 2.2915 -7.7524
Variance Components:
estim sqrt nlvls fixed factor R
sigma^2.1 0.0000 0.0000 36 no phylo yes
sigma^2.2 0.0014 0.0377 29 no paper no
sigma^2.3 0.0244 0.1562 102 no id no
Test for Heterogeneity:
Q(df = 101) = 765.9358, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2564 0.0219 11.6818 <.0001 0.2134 0.2994 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All this effort, there is no variation due to phylogeny! So we do not need this phylogeny term (i.e. phylo).
Also, this multilevel model can be considered as a comparative phylogenetic method. There are quite a few things you need to know and to be careful of, which I cannot cover (e.g. we assumed that the Brownian-motion model of evolution in the model above - what does this even mean?). But Will and I have written a nice ‘primer’ so please read that primer – Cornwell & Nakagawa (2017)
Unfortunately, there are other types of non-independence from what covered here. We summaries all types in our recent paper – Noble et al. (2017). So read this as well if you are interested.
30.3.3 Complex robust models (the last section!)
We have a multilevel version of the robust model too. It is easy to fit using the function rma.mv (we do not include phylo as it did not explain any variance).
Multivariate Meta-Analysis Model (k = 102; method: REML)
logLik Deviance AIC BIC AICc
5.7208 -11.4415 -5.4415 2.4039 -5.1941
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0014 0.0377 29 no paper
sigma^2.2 0.0244 0.1562 102 no id
Test for Heterogeneity:
Q(df = 101) = 765.9358, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.2087 0.0479 4.3600 <.0001 0.1149 0.3025 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1I think this is the model we have been looking for, i.e. our final model. At least for this data set.
30.4 Further help (references)
Any questions? Or email me (s(-dot-)nakagawa(-at-)unsw(-dot-)edu(-dot-)au). Also visit our website
Go to the metafor package’s website. There you find many worked examples.
Cornwell, W., and S. Nakagawa. 2017. Phylogenetic comparative methods. Current Biology 27:R333-R336.
Nakagawa, S., D. W. A. Noble, A. M. Senior, and M. Lagisz. 2017. Meta-evaluation of meta-analysis: ten appraisal questions for biologists. BMC Biology 15:18.
Nakagawa, S., and E. S. A. Santos. 2012. Methodological issues and advances in biological meta-analysis. Evolutionary Ecology 26:1253-1274.
Noble, D. W. A., M. Lagisz, R. E. O’Dea, and S. Nakagawa. 2017. Nonindependence and sensitivity analyses in ecological and evolutionary meta-analyses. Molecular Ecology 26:2410-2425.
Authors: Shinichi Nakagawa and Malgorzata (Losia) Lagisz
Year: 2016