07-metaregression.Rmd

# Meta-Regression {#metareg}

![](regressionbild.jpg)


Conceptually, **Meta-Regression** does not differ much from a **subgroup analysis**. In fact, subgroup analyses with more than two groups are nothing more than a meta-regression with categorial covariates. Meta-regression with **continuous moderators** reveals whether values of this continuous variable are linearly associated with the effect size.



```{block,type='rmdinfo'}
**The idea behind meta-regression**

In a conventional regression, we specify a model predicting the dependent variable $y$, across $_i$ participants, based on their values on $p$ predictor variables, $x_{i1} \dots x_{ip}$. The residual error is referred to as $\epsilon_i$. A standard regression equation therefore looks like this:

$$y_i=\beta_0 + \beta_1x_{1i} + ...+\beta_px_{pi} + \epsilon_i$$

In a meta-regression, we want to estimate the **effect size** $\theta$ of several studies $_k$, as a function of between-studies moderators. There are two sources of heterogeneity: sampling error, $\epsilon_k$, and between-studies heterogeneity, $\zeta_k$ so our regression looks like this:

$$\theta_k = \beta_0 + \beta_1x_{1k} + ... + \beta_px_{pk} + \epsilon_k + \zeta_k$$

You might have seen that when estimating the effect size $\theta_k$ of a study $k$ in our regression model, there are two **extra terms in the equation**, $\epsilon_k$ and $\zeta_k$. The same terms can also be found in the equation for the random-effects-model in [Chapter 4.2](#random). The two terms signify two types of **independent errors** which cause our regression prediction to be **imperfect**. The first one, $\epsilon_k$, is the sampling error through which the effect size of the study deviates from its "true" effect. The second one, $\zeta_k$, denotes that even the true effect size of the study is only sampled from **an overarching distribution of effect sizes** (see the Chapter on the [Random-Effects-Model](#random) for more details). In a **fixed-effect-model**, we assume that all studies actually share the **same true effect size** and that the **between-study heterogeneity** $\tau^2 = 0$. In this case, we do not consider $\zeta_k$ in our equation, but only $\epsilon_k$.

As the equation above has includes **fixed effects** (the $\beta$ coefficients) as well as **random effects** ($\zeta_k$), the model used in meta-regression is often called **a mixed-effects-model**. Mathematically, this model is identical to the **mixed-effects-model** we described in [Chapter 7](#subgroup) where we explained how **subgroup analyses** work.

Indeed **subgroup analyses** are nothing else than a **meta-regression** with a **categorical predictor**. For meta-regression, these subgroups are then **dummy-coded**, e.g.

$$ D_k =  \{\begin{array}{c}0:ACT \\1:CBT \end{array}$$

$$\hat \theta_k = \theta + \beta x_{k} + D_k \gamma + \epsilon_k + \zeta_k$$

In this case, we assume the same **regression line**, which is simply "shifted" **up or down for the different subgroups** $D_k$.

```

```{r, echo=FALSE, fig.width=6,fig.cap="Visualisation of a Meta-Regression with dummy-coded categorial predictors",fig.align='center'}
library(metaforest)
df <- curry
library(png)
library(grid)
img <- readPNG("dummy.PNG")
grid.raster(img)
```

<br><br>

```{block,type='rmdinfo'}
**Assessing the fit of a regression model**

To evaluate the **statistical significance of a predictor**, we conduct a **t-test** of its $\beta$-weight.

$$ t=\frac{\beta}{SE_{\beta}}$$

Which provides a $p$-value telling us if a variable significantly predicts effect size differences in our regression model.

If we fit a regression model, our aim is to find a model **which explains as much as possible of the current variability in effect sizes** we find in our data.

In conventional regression, $R^2$ is commonly used to quantify the percentage of variance in the data explained by the model, as a percentage of total variance (0-100%). As this measure is commonly used, and many researchers know how to to interpret it, we can also calculate a $R^2$ analog for meta-regression using this formula:

$$R_2=\frac{\hat\tau^2_{REM}-\hat\tau^2_{MEM}}{\hat\tau^2_{REM}}$$

Where $\hat\tau^2_{REM}$ is the estimated total heterogenetiy based on the random-effects-model and $\hat\tau^2_{REM}$ the total heterogeneity of our mixed-effects regression model.

NOTE however, that $R^2$ refers to variance explained in the observed data. The more predictors you add, the better your model will explain the observed data. But this can decrease the generalizability of you model, through a process called **overfitting**: You're capturing noise in your dataset, not true effects that exist in the real world.
```

<br><br>

---

## Meta-regression in R

Meta-regressions can be conducted in R using the `rma` function in `metafor`. To show the similarity between `subgroup` analysis and `meta-regression`, consider the code for our regression-specified subgroup analysis again:

```{r, eval = FALSE}
m_dummy <- rma(yi = d, vi = vi, mods = ~ Country, data = df)
```

This syntax used **"country"** as a dummy-coded variable. R will always do this for `factor` variables. It is possible to drop the intercept and estimate the means for all groups instead by specifying `~Country-1`, and it is even possible to change the default dummy coding if you want to compare a set of groups with another set of groups: [a tutorial](https://rstudio-pubs-static.s3.amazonaws.com/65059_586f394d8eb84f84b1baaf56ffb6b47f.html).

<br><br>

**Continuous variables**

Imagine you want to check if the **proportion of male participants** is associated with effect size. The variable `sex` contains this information. You can use this predictor in a meta-regression:

```{r,echo=FALSE}
# load("Meta_Analysis_Data.RData")
# madata<-Meta_Analysis_Data
# madata$pub_year<-c(2001,2002,2011,2013,2013,2014,1999,2018,2001,2002,2011,2013,2013,2014,1999,2018,2003,2005)
# madata$pub_year<-as.numeric(madata$pub_year)
# m.pubyear<-metagen(TE,seTE,studlab = paste(Author),comb.fixed = FALSE,data=madata)
```

```{r}
m_reg <- rma(yi = d,
             vi = vi,
             mods = ~sex,
             data = df)
m_reg
```

As you can see from the output, `sex` was now included as a **predictor**, but it is not significantly associated with the effect size ($p=.1416$).

<br><br>

---


## Plotting regressions

To plot our meta-regression output, we can make a **bubble plot** using ggplot. This is essentially a weighted scatter plot, where the size of the scatter is mapped to the inverse SE of each effect size, which means the area of the scatter is proportional to the inverse variance:

```{r}
library(ggplot2)
df$weights <- 1/sqrt(df$vi)
# Specify basic plot, mapping sex to the x-axis, effect size 'd' to the y-axis,
# and 'weights' to the weight parameter.
ggplot(df, aes(x = sex, y = d, size = weights)) +
  geom_point(shape = 1) + # Add scatter
  geom_abline(intercept = m_reg$b[1], slope = m_reg$b[2]) + # Add regression line
  theme_bw() + # Apply black and white theme
  theme(legend.position = "none") # Remove legend
```

<br><br>

---

## Multiple Meta-Regression

![](mvreg.jpg)

Previously, we only considered the scenario in which we use **one predictor** $\beta_1x_1$ in our meta-regression. When we add more than one predictor, we're using **multiple meta-regression**. In multiple meta-regression we use several moderators (variables) to predict the outcome (effect sizes). When we look back at the **general meta-regression** formula we defined before, we actually see that the formula already provides us with this feature through the $\beta_nx_{pk}$ part. Here, the parameter $p$ denotes that we can include $p$ more predictors/variables into our meta-regression, making it a multiple meta-regression.

Imagine, for example, that we expect the effect size to be determined by `sex`, but also by the type of outcome variable that was measured, `outcomecode`. We could include both predictors in a mutliple meta-regression as follows:

```{r}
m_multi <- rma(yi = d,
               vi = vi,
               mods = ~sex+outcomecode,
               data = df)
m_multi
```

We will discuss a few important pitfalls in multiple meta-regression and how we can build multiple meta-regression models which are robust and trustworthy. But first, let's cover another important feature of multiple meta-regression: **interactions**.


## Interactions

So far, in our multiple meta-regression model, we only considered the case where we have multiple predictor variables $x_1,x_2, ... x_p$, and along with their predictor estimates $\beta_p$, **add them together** to calculate our estimate of the true effect size $\hat \theta_k$ for each study $k$. In multiple meta-regression models, however, we can not only model such **additive relationships**. We can also model so-called **interactions**. Interactions mean that the **relationship** between one **predictor variable** (e.g., $x_1$) and the **estimated effect size** is different for different values of another predictor variable (e.g. $x_2$).

Imagine a scenario where we want to model two predictors and their relationship to the effect size: the **sex** ($x_1$) of participants and the **donor type** ($x_2$) of participants in a study (either Anxious or Typical).

As we described before, we can now imagine a meta-regression model in which we combine these two predictors $x_1$ and $x_2$ and assume an additive relationship. We can do this by simply adding them:

$$\theta_k = \beta_0 + \beta_1x_{1k} + \beta_2x_{2k} + \epsilon_k + \zeta_k$$

We could now ask ourselves if the relationship of participants' sex **depends** on the type of participants (Anxious or Typical; $x_2$). For example, maybe sex predicts effect size more strongly in Anxious populations? To examine such questions, we can add an **interaction term** to our meta-regression model. This interaction term lets predictions of $x_1$ vary for different values of $x_2$ (and vice versa). We can denote this additional interactional relationship in our model by introducing a third predictor, $\beta_3$, which captures this interaction $x_{1k}x_{2k}$ we want to test in our model:

$$\theta_k = \beta_0 + \beta_1x_{1k} + \beta_2x_{2k} + \beta_3x_{1k}x_{2k}+ \epsilon_k + \zeta_k$$

The R-syntax for this interaction is:

```{r}
m_int <- rma(yi = d,
             vi = vi,
             mods = ~sex*donorcode,
             data = df)
m_int
```

<br></br>

## Common pitfalls {#pitfalls} 

As we have mentioned before, multiple meta-regression, while **very useful when applied properly**, comes with **certain caveats** we have to know and consider when fitting a model. Indeed, some argue that (multiple) meta-regression is **often improperly used and interpreted in practice**, leading to a low validity of many meta-regression models [@higgins2004controlling]. Thus, there are some points we have to keep in mind when fitting multiple meta-regression models, which we will describe in the following.

```{block,type='rmdachtung'}
**Overfitting: seeing a signal when there is none**

To better understand the risks of (multiple) meta-regression models, we have to understand the concept of **overfitting**. Overfitting occurs when we build a statistical model which fits the data **too closely**. In essence, this means that we build a statistical model which can predict the data at hand very well, but performs bad at predicting future data it has never seen before. This happens if our model assumes that some variation in our data **stems from a true "signal" in our data, when in fact we only model random noise** [@iniesta2016machine]. As a result, our statistical model produces **false positive results**: it sees relationships where there are none.

![Illustration of an overfitted model vs. model with a good fit](overfitting.png)

Regression methods, which usually utilize minimization or maximization procedures such as Ordinary Least Squares or Maximum Likelihood estimatation, can be prone to overfitting [@gigerenzer2004mindless; @gigerenzer2008rationality]. Unfortunately, the risk of building a **non-robust model, which produces false-positive results**, is **even higher** once we go from conventional regression to **meta-regression** [@higgins2004controlling]. There are several reasons for this:

1.  In Meta-Regression, our **sample of data is mostly small**, as we can only use the synthesized data of all analyzed studies $k$.
2.  As our meta-analysis aims to be a **comprehensive overview of all evidence**, we have no **additional data** on which we can "test" how well our regression model can predict high or low effect sizes.
3.  In meta-regressions, we have to deal with the potential presence of effect size heterogeneity ([see Chapter 6](#heterogeneity)). Imagine a case in which we have two studies with different effect sizes and non-overlapping confidence intervals. Every variable which has different values for the different studies might be a potential explanation for effect size difference we find, while it seems straightforward that **most of such explanations are spurious** [@higgins2004controlling].
4.  Meta-regression as such, and multiple meta-regression in particular, make it very easy to **"play around" with predictors**. We can test numerous meta-regression models, include many more predictors or remove them in an attempt to explain the heterogeneity in our data. Such an approach is of course tempting, and often found in practice, because we, as meta-analysts, want to find a significant model which explains why effect sizes differ [@higgins2002statistical]. However, such behavior **massively increases the risk of spurious findings** [@higgins2004controlling], because we can change parts of our model indefinitely until we find a significant model, which is then very likely to be overfitted (i.e., it mostly models statistical noise).

**Some guidelines have been proposed to avoid an excessive false positive rate when building meta-regression models:**

- Minimize the number of investigated predictors a priori. Predictor selection should be based on **predefined scientific or theoretical questions** we want to answer in our meta-regression.
- When evaluating the fit of a meta-regression model, we prefer models which achieve a good fit with less predictors. Use fit indices that balance fit and parsimony, such as the *Akaike* and *Bayesian information criterion*, to determine which model to retain if you compare several models.
- When the number of studies is low (which is very likely to be the case), and we want to compute the significance of a predictor, you can use the Knapp-Hartung adjustment to obtain more reliable estimates [@higgins2002statistical], by specifying `test = "knha` when calling `rma()`.

```


<br></br>


---