03-pooling_effect_sizes.Rmd

# Pooling Effect Sizes {#pool}

![](pooling.jpg)

Now, let's get to the core of every Meta-Analysis: **pooling your effect sizes** to get one overall effect size estimate of the studies.

```{block,type='rmdinfo'}
When pooling effect sizes in Meta-Analysis, there are two basic approaches: the **Fixed-Effect-Model**, or the **Random-Effects-Model** [@borenstein2011]. The fixed effect model assumes that one true effect size exists in the population; the random effects model assumes that the true effect varies (and is normally distributed). There is an extensive debate on which model fits best in which context [@fleiss1993review], with no clear consensus in sight. Although it has been recommended to **only resort to the Random-Effects-Pooling model** in clinical psychology and the health sciences [@cuijpers2016meta], we will describe how to conduct both in R here.

Both of these models only require an **effect size**, and a **dispersion (variance)** estimate for each study, of which the inverse is taken. This is why the methods are often called **generic inverse-variance methods**.
```

We will describe in-depth how to conduct meta-analyses in R with **continuous variables** (such as effect sizes), as these are the most common ones in psychology and the health science field. Later on, we will present briefly how to do meta-analyses with **binary outcome** data too, which might be important if you're focusing on prevention trials.

For these meta-analyses, we'll use the `metafor` package [@viechtbauer2010conducting]. In [Section 2.1](#RStudio), we showed how to install the package. We load a different package: The `metaforest` package, which in turn loads the `metafor` package, and contains the example data for this tutorial.

```{r, warning=FALSE,message=FALSE}
library(metaforest)
```

```{r, echo=FALSE,warning=FALSE,message=FALSE}
library(tidyverse)
library(knitr)
```

<br><br>

---

## Fixed-Effects-Model {#fixedef}


```{block, type='rmdinfo'}
**The idea behind the fixed-effects-model**

The fixed-effects-model assumes that all observed effect sizes stem from a single *true* population effect [@borenstein2011]. To calculate the overall effect, we therefore average all effect sizes, but give studies with greater precision a higher weight. Precision relates to the fact that studies with a larger **N** will provide more accurate estimates of the true population effect, as reflected by a smaller **Standard Error** of the effect size estimate.

For this weighing, we use the **inverse of the variance** $1/\hat\sigma^2_k$ of each study $k$. We then calculate a weighted average of all studies, our fixed effect size estimator $\hat\theta_F$:
```

\begin{equation}
\hat\theta_F = \frac{\sum\limits_{k=1}^K \hat\theta_k/ \hat\sigma^2_k}{\sum\limits_{k=1}^K 1/\hat\sigma^2_k}
\end{equation}


First, let's assume you already have a dataset with the **calucated effects and SE** for each study. The `curry` data set will do (you've used it before to practice loading data from Excel).

```{r,echo=FALSE}
df<-curry
```

This dataset has **continuous outcome data**. As our effect sizes are already calculated, we can directly use the `rma` function. For this function, we can specify loads of parameters, all of which you can accessed by typing `?rma` in your console once the `metafor` package is loaded, or selecting the function and pressing F1.
  
**Here is a table with the most important parameters for our code:**

```{r,echo=FALSE}
i<-c("yi","vi","method")
ii<-c("A vector with the effect sizes",
      "A vector with the sampling variances",
      "A character string, indicating what type of meta-analysis to run. FE runs a fixed-effect model")
ms<-data.frame(i,ii)
names<-c("Parameter", "Function")
colnames(ms)<-names
kable(ms)
```

Let's conduct our first fixed-effects-model Meta-Analysis. We we will give the results of this analysis the simple name `m`.

```{r}
m <- rma(yi = df$d,     # The d-column of the df, which contains Cohen's d
         vi = df$vi,    # The vi-column of the df, which contains the variances
         method = "FE") # Run a fixed-effect model
m

```

We now see a summary of the results of our Meta-Analysis, including

* The total **number of included studies** (k)
* The **overall effect** (in our case, *g* = 0.2059) and its confidence interval and p-value
* The *Q*-test of heterogeneity

The object `m` is a **list**: A collection of different types of data. A data.frame is a special kind of list, where every entry must have the same length. Using the `$` command, we can look at parts of the list directly. For example, the estimates of heterogeneity:

```{r}
m$I2
m$tau2
```

There are many ways to export results from R. The function `sink()` redirects
output, which would normally be printed to the console, to a text file:

```{r,eval=FALSE}
sink("results.txt")
print(m)
sink()
```


## Random-Effects-Model {#random}

We can only use the fixed-effect-model when we can assume that **all included studies tap into one true effect size**. In practice this is hardly ever the case: interventions may vary in certain characteristics, the sample used in each study might be slightly different, or its methods. A more appropriate assumption in these cases might be that the true effect size follows a normal distribution.

```{block,type='rmdinfo'}
**The Idea behind the Random-Effects-Model**

In the Random-Effects-Model, we want to account for our assumption that the population effect size is normally distributed [@schwarzer2015meta]. The random-effects-model works under the so-called **assumption of exchangeability**. 

This means that in Random-Effects-Model Meta-Analyses, we not only assume that effects of individual studies deviate from the true intervention effect of all studies due to sampling error, but that there is another source of variance introduced by the fact that the studies do not stem from one single population, but are drawn from a "universe" of populations. We therefore assume that there is not only one true effect size, but **a distribution of true effect sizes**. We therefore want to estimate the mean and variance of this distribution of true effect sizes.

The fixed-effect-model assumes that when the observed effect size $\hat\theta_k$ of an individual study $k$ deviates from the true effect size $\theta_F$, the only reason for this is that the estimate is burdened by (sampling) error $\epsilon_k$.

$$\hat\theta_k = \theta_F + \epsilon_k$$


While the random-effects-model assumes that, in addition, there is **a second source of error** $\zeta_k$.This second source of error is introduced by the fact that even the true effect size $\theta_k$ of our study $k$ is also only part of an over-arching distribution of true effect sizes with the mean $\mu$ [@borenstein2011]. 

```

```{r, echo=FALSE, fig.width=6,fig.height=4,fig.align='center'}
library(png)
library(grid)
img <- readPNG("density.png")
grid.raster(img)
```

<p style="text-align: center;">
*An illustration of  parameters of the random-effects-model*
</p>
 
```{block,type='rmdinfo'}
The formula for the random-effects-model therefore looks like this:

$$\hat\theta_k = \mu + \epsilon_k + \zeta_k$$

When calculating a random-effects-model meta-analysis, where therefore also have to take the error $\zeta_k$ into account. To do this, we have to **estimate the variance of the distribution of true effect sizes**, which is denoted by $\tau^{2}$, or *tau^2^*. There are several estimators for $\tau^{2}$, all of which are implemented in `metafor`.
```

```{block,type='rmdachtung'}
Even though it is **conventional** to use random-effects-model meta-analyses in psychological outcome research, applying this model is **not undisputed**. The random-effects-model pays **more attention to small studies** when pooling the overall effect in a meta-analysis [@schwarzer2015meta]. Yet, small studies in particular are often fraught with **bias** (see [Chapter 8.1](#smallstudyeffects)). This is why some have argued that the fixed-effects-model should be nearly always preferred [@poole1999; @furukawa2003].
```

### Estimators for *tau^2^* in the random-effects-model {#tau2}

Operationally, conducting a random-effects-model meta-analysis in R is not so different from conducting a fixed-effects-model meta-analyis. Yet, we do have choose an estimator for $\tau^{2}$. Here are the estimators implemented in `metafor`, which we can choose using the `method` argument when calling the function.

```{r,echo=FALSE}
i<-c("DL","PM","REML","ML","HS","SJ","HE","EB")
ii<-c("DerSimonian-Laird","Paule-Mandel","Restricted Maximum-Likelihood","Maximum-likelihood","Hunter-Schmidt","Sidik-Jonkman","Hedges","Empirical Bayes")
ms<-data.frame(i,ii)
names<-c("Code", "Estimator")
colnames(ms)<-names
kable(ms)
```

```{block,type='rmdinfo'}
**Which estimator should I use?**

All of these estimators derive $\tau^{2}$ using a slightly different approach, leading to somewhat different pooled effect size estimates and confidence intervals. If one of these approaches is more or less biased often depends on the context, and parameters such as the number of studies $k$, the number of participants $n$ in each study, how much $n$ varies from study to study, and how big $\tau^{2}$ is. 

An overview paper by Veroniki and colleagues [@veroniki2016methods] provides an excellent summary on current evidence which estimator might be more or less biased in which situation. The article is openly accessible, and you can read it [here](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4950030/). This paper suggests that the **Restricted Maximum-Likelihood** estimator performs best, and it is the default estimator in metafor.
```

### Conducting the analysis {#random.precalc}

Random-effects meta-analyses are very easy to code in R. Compared to the fixed-effects-model [Chapter 5.1](#fixedef), we can simply remove the `method = "FE"` argument, if we want to use the default REML estimator:

```{r}
m_re <- rma(yi = df$d,     # The d-column of the df, which contains Cohen's d
            vi = df$vi)    # The vi-column of the df, which contains the variances
m_re

```


The output shows that our estimated effect is $g=0.2393$, and the 95% confidence interval stretches from $g=0.16$ to $0.32$ (rounded). The estimated heterogeneity is $\tau^2 = 0.06$. The percentage of variation across effect sizes that is due to heterogeneity rather than change is estimated at $I^2 = 67.77\%$.

The summary effect estimate is different (and larger) than the one we found in the fixed-effects-model meta-analysis in [Chapter 5.1](#fixedef) ($g=0.21$).