13-power_analysis.Rmd

# Power Analysis

![](poweranalysis.jpg)

A big asset (and probably one of the reasons why meta-analysis can be helpful in practical research) of meta-analyses is that **they allow for large data to be combined** to attain a more precise pooled effect. **Lack of statistical power**, however, may still play an important role, even in meta-analysis. This is particularly true for the clinical field, where it is often the case that only **few studies are available for synthesis**. The median number of included studies in the *Cochrane Database for Systematic Reviews*, for example, is six [@borenstein2011]. This is even more grave once we consider that (1) many meta-analysts will also want to perform **subgroup analyses and meta-regression**, for which even more power is required, and (2) many meta-analyses have **high heterogeneity**, which reduces our precision, and thus our power.

Power is directly related to the **Type II error level** ($\beta$) we defined: $Power = 1- \beta$. It is common practice to set our **Type I error level** ($\alpha$) to $\alpha=0.05$, and thus to assume that the **Type I error is four times as grave as the Type II error** (i.e., falsely finding an effect while there is no effect in reality is four times as bad as not finding an effect while there is one in reality). The **Type II error** is therefore set at $\beta=0.20$, and the power should thus be $1-\beta=1-0.20=80\%$.

<br><br>

```{block,type='rmdinfo'}
**What assumptions should i make for my meta-analysis?**

While researchers conducting primary studies can **plan the size of their sample based on the effect size they want to find**, the situation is a little different in meta-analysis, where we can only work with the published material. However, we have some **control over the number of studies we want to include in our meta-analysis** (e.g., through more leniently or strictly defined inclusion criteria). Therefore, we can change our power to some extent by including more or less studies into the meta-analysis. There are **four things we have to make assumptions about when assessing the power of our meta-analysis a priori**.

* The **number of included or includable studies** $k$
* The **overall size of the studies we want to include** (are the studies in the field rather small or large?)
* The **effect size we want to determine**. This is particularly important, as we have to make assumptions about how **big an effect size has to be to still be clinically meaningful**. One study calculated that for interventions against depression, even effects as small as $SMD=0.24$ may still be meaningful for patients [@cuijpers2014threshold]. If we want to study **negative effects of an intervention** (e.g., death or symptom deterioration), even **very small effect sizes are extremely important and should be detected**.
* The **heterogeneity** of our studies' effect sizes, as this also affects the precision of our meta-analysis, and thus its potential to find significant effects.

Besides these parameters, it is also important to think about other analyses, such as the **subgroup analyses** we want to conduct. How many studies are there for each subgroup, and what effects do we want to find in the subgroups? This is particularly important if we **hypothesize that an intervention is not effective in a subgroup of patients**, because we do not want to falsely find a treatment to be ineffective simply because the power was insufficient.
```

```{block,type='rmdachtung'}
**Post-hoc power tests: the abuse of power**

Please note that power analyses should always be conducted **a priori**, meaning *before* you perform the meta-analysis.

Power analyses conducted *after* an analysis ("post hoc") are fundamentally flawed [@hoenig2001abuse], as they suffer from the so-called **"power approach paradox"**, in which an analysis yielding no significant effect is thought to show more evidence that the null hypothesis is true when the p-value is smaller, since then, the power to detect a true effect would be higher.

```

---

## Fixed-Effect Model {#fixed.power}

To determine the **power** of a meta-analysis **under the fixed-effect model**, we have to assume the **true value of a distribution when the alternative hypothesis is correct** (i.e., when there is an effect). For power analysis in a conventional study, this distribution is $Z$. Follwing Borenstein et al. [@borenstein2011], we will call the true value $\lambda$ here to make clear that we are dealing with a meta-analysis, and not a primary study. $\lambda$ is defined as:

$$\lambda=\frac{\delta}{\sqrt{V_{\delta}}}$$

Where $\delta$ is the **true effect size** and $V_{\delta}$ its variance.

$V_{\delta}$ can be calculated for meta-analysis using the fixed-effect model with this formula:


$$V_{\delta}=\frac{\frac{n_1+n_2}{n_1xn_2}+\frac{d^2}{2(n_1+n_2)}}{k}$$

Where $k$ are all the included studies, and $n_1$ and $n_2$ are the **average sample sizes in each trial arm we assume across our studies**.

Assuming a normal distribution and using $\lambda$, we can calculate the Power:

$$Power = 1- \beta$$
$$Power = 1- \Phi(c_{\alpha}-\lambda)+\Phi(-c_{\alpha}-\lambda) $$

Where $c_{\alpha}$ is the critical value of a $Z$-distribution. $\Phi$ is the **standard normal density function**, which we we need to calcuate the power using this equation:
$$\Phi(Z)=\frac{1}{\sqrt {2\pi}}e^{-\frac{Z^2}{2}}$$

Luckily, you don't have too think about these statistical details too much, as we have prepared a **function** for you with which you can easily conduct a **power analysis** using the fixed-effect model yourself. The function is called `power.analysis.function` and its code can be found below.

Again, R doesn't know this function yet, so we have to let R learn it by **copying and pasting** the code underneath **in its entirety** into the **console** on the bottom left pane of RStudio, and then hit **Enter ⏎**.

```{r}
power.analysis.fixed<-function(d,k,n1,n2,p){

  n1<-n1
  n2<-n2
  d<-d
  k<-k
  p<-p
  title<-c("Power for a fixed-effect meta-analysis:")

  v.d<-((n1+n2)/(n1*n2))+((d*d)/(2*(n1+n2)))
  v.m<-v.d/k
  lambda<-(d/sqrt(v.m))
  plevel<-1-(p/2)
  zval<-qnorm(p=plevel, 0,1)
  power<-1-(pnorm(zval-lambda))+(pnorm(-zval-lambda))
  return(power)
}
```

---

**For this function, we have to specify the following parameters:**

```{r,echo=FALSE}
library(kableExtra)
Package<-c("d","k","n1","n2","p")
Description<-c("The effect size we want to be able to detect","The number of studies we will likely be able to include into our analysis","The average number of participats we assume in the intervention arms of our included studies","The average number of participats we assume in the control arms of our included studies","The Type I error rate (p-level). It is common to use 'p=0.05'")
m<-data.frame(Package,Description)
names<-c("Parameter", "Description")
colnames(m)<-names
kable(m)
```

Now, let's give an example. I assume that an effect of $d=0.30$ is likely and meaningful for the field of my meta-analysis. I also assume that on average, the studies in my analysis will be rather small, with 25 participants in each trial arm, and that there will be 10 studies in my analysis. I will set the $\alpha$-level to 0.05, as is convention.

```{r, eval=FALSE}
power.analysis.fixed(d=0.30,k=10,n1=25,n2=25,p=0.05)
```

The output of the function is:

```{r,echo=FALSE}
power.analysis.fixed(d=0.30,k=10,n1=25,n2=25,p=0.05)
```

Meaning that my power is 91%. This is more than the desired 80%, so given that my assumptions are remotely true, my meta-analysis will have **sufficient power using the fixed-effect model to detect a clinically relevant effect if it exists**.

So, if i assume an effect of $d = 0.30$ in this example, i am lucky. If we play around with the effect size a little, however, while holding the other paramters constant, this can look very different.

```{r,fig.width=5,fig.align='center',echo=FALSE,message=FALSE}
library(ggplot2)
library(reshape)
k <- seq(0, 50, length=1000)
pow.vals01<-lapply(k,function(k) power.analysis.fixed(d=0.10,k=k,n1=25,n2=25,p=0.05))
pow.vals02<-lapply(k,function(k) power.analysis.fixed(d=0.20,k=k,n1=25,n2=25,p=0.05))
pow.vals03<-lapply(k,function(k) power.analysis.fixed(d=0.30,k=k,n1=25,n2=25,p=0.05))
pow.vals01<-as.numeric(pow.vals01)
pow.vals02<-as.numeric(pow.vals02)
pow.vals03<-as.numeric(pow.vals03)
data<-data.frame(k,pow.vals01,pow.vals02,pow.vals03)
ggplot()+
  geom_line(data = data, aes(x = k, y = pow.vals01), color = "blue",size=2) +
  geom_line(data = data, aes(x = k, y = pow.vals02), color = "red",size=2) +
  geom_line(data = data, aes(x = k, y = pow.vals03), color = "green",size=2) +
  xlab('Number of Studies') +
  ylab('Power')+
  scale_y_continuous(labels = scales::percent)+
  theme(
        axis.line= element_line(color = "black",size = 1,linetype = "solid"),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(linetype="solid",
                                         colour ="black"),
        legend.title = element_blank(),
        legend.key.size = unit(0.75,"cm"),
        legend.text=element_text(size=14))+
annotate("text", x = 2, y = 0.7, label = "d = 0.3",size=5)+
  annotate("text", x = 14, y = 0.6, label = "d = 0.2",size=5)+
  annotate("text", x = 20, y = 0.24, label = "d = 0.1",size=5)+
  geom_hline(yintercept=0.8,linetype="dashed")
```

As you can see from this plot, sufficient power (see the **dashed line**) is soon reached for $d=0.30$, even if only few studies are included. If i assume a smaller effect size of $d=0.10$, however, **even 50 studies will not be sufficient to find a true effect**. 

---

## Random-Effects Model

For power analyses under the **random-effects model**, the formula to calculate the variance of my true mean effect looks slightly different:

$$V_{\delta}^*=\frac{V_Y+\tau^2}{k}$$

We see that again, $tau^2$ has to be included to take the **between-study heterogeneity into account** (see [Chapter 4.2](#random) for more details). However, i **do not know the between-study heterogeneity of my analysis** before i perform it, so what value should i assume?

According to Hedges and Pigott [@hedges2004power], the follwing formulae may be used to calculate the power in the random-effect model assuming **small, moderate or high heterogeneity**:


**Small heterogeneity:**

$$V_{\delta}^*=1.33\times\frac{V_Y}{k}$$

**Moderate heterogeneity:**

$$V_{\delta}^*=1.67\times\frac{V_Y}{k}$$

**Large heterogeneity:**

$$V_{\delta}^*=2\times\frac{V_Y}{k}$$

Again, you don't have to worry about the statistical details here. We have put the entire calculations into the `power.analysis.random` function, which can be found below.

Again, R doesn't know this function yet, so we have to let R learn it by **copying and pasting** the code underneath **in its entirety** into the **console** on the bottom left pane of RStudio, and then hit **Enter ⏎**.

```{r}
power.analysis.random<-function(d,k,n1,n2,p,heterogeneity){

  n1<-n1
  n2<-n2
  d<-d
  k<-k
  p<-p
  heterogeneity<-heterogeneity
  
  if(heterogeneity=="low"){

  v.d<-((n1+n2)/(n1*n2))+((d*d)/(2*(n1+n2)))
  v.m<-v.d/k
  v.m<-1.33*v.m
  lambda<-(d/sqrt(v.m))
  plevel<-1-(p/2)
  zval<-qnorm(p=plevel, 0,1)
  power<-1-(pnorm(zval-lambda))+(pnorm(-zval-lambda))
  return(power)
  }
  
  if(heterogeneity=="moderate"){
    
      v.d<-((n1+n2)/(n1*n2))+((d*d)/(2*(n1+n2)))
  v.m<-v.d/k
  v.m<-1.67*v.m
  lambda<-(d/sqrt(v.m))
  plevel<-1-(p/2)
  zval<-qnorm(p=plevel, 0,1)
  power<-1-(pnorm(zval-lambda))+(pnorm(-zval-lambda))
  return(power)
  }
    
    if(heterogeneity=="high"){
    
      v.d<-((n1+n2)/(n1*n2))+((d*d)/(2*(n1+n2)))
  v.m<-v.d/k
  v.m<-2*v.m
  lambda<-(d/sqrt(v.m))
  plevel<-1-(p/2)
  zval<-qnorm(p=plevel, 0,1)
  power<-1-(pnorm(zval-lambda))+(pnorm(-zval-lambda))
  return(power)
  }
    
}
```

Now you are set and ready to use the function. I will assume the same parameters from [before](#fixed.power), but this time i will also have to specify the `heterogeneity` in the function, which can take the values `"low"`, `"moderate"` and `"high"`. I will choose `"moderate"` for this example.

```{r,eval=FALSE}
power.analysis.random(d=0.30,k=10,n1=25,n2=25,p=0.05,
                      heterogeneity = "moderate")

```

The output i get is:

```{r,echo=FALSE}
power.analysis.random(d=0.30,k=10,n1=25,n2=25,p=0.05,heterogeneity = "moderate")

```

Interestingly, we see that this value is 73%, which is smaller than the value of 91% which was calculated using the **fixed-effect model**. The value is also below 80%, meaning that i would not have optimal power to find the desired effect of $d=0.30$ to be statistically significant if it exists.

This has to do with the **larger heterogeneity** i assume in this simulation, which decreases the precision of my effect size estimate, and thus increases my need for statistical power.

**The graph below visualizes this relationship:**

```{r,fig.width=5,fig.align='center',echo=FALSE,fig.cap="Power in the random-effects-model. Darker colors indicate higher heterogeneity",message=FALSE}
library(ggplot2)
library(reshape)

k <- seq(0, 50, length=1000)
pow.vals01<-lapply(k,function(k) power.analysis.random(d=0.10,k=k,n1=25,n2=25,p=0.05,heterogeneity = "moderate"))
pow.vals02<-lapply(k,function(k) power.analysis.random(d=0.20,k=k,n1=25,n2=25,p=0.05,heterogeneity = "moderate"))
pow.vals03<-lapply(k,function(k) power.analysis.random(d=0.30,k=k,n1=25,n2=25,p=0.05,heterogeneity = "moderate"))
pow.vals01<-as.numeric(pow.vals01)
pow.vals02<-as.numeric(pow.vals02)
pow.vals03<-as.numeric(pow.vals03)
data1<-data.frame(k,pow.vals01,pow.vals02,pow.vals03)

k <- seq(0, 50, length=1000)
pow.vals01<-lapply(k,function(k) power.analysis.random(d=0.10,k=k,n1=25,n2=25,p=0.05,heterogeneity = "low"))
pow.vals02<-lapply(k,function(k) power.analysis.random(d=0.20,k=k,n1=25,n2=25,p=0.05,heterogeneity = "low"))
pow.vals03<-lapply(k,function(k) power.analysis.random(d=0.30,k=k,n1=25,n2=25,p=0.05,heterogeneity = "low"))
pow.vals01<-as.numeric(pow.vals01)
pow.vals02<-as.numeric(pow.vals02)
pow.vals03<-as.numeric(pow.vals03)
data2<-data.frame(k,pow.vals01,pow.vals02,pow.vals03)

k <- seq(0, 50, length=1000)
pow.vals01<-lapply(k,function(k) power.analysis.random(d=0.10,k=k,n1=25,n2=25,p=0.05,heterogeneity = "high"))
pow.vals02<-lapply(k,function(k) power.analysis.random(d=0.20,k=k,n1=25,n2=25,p=0.05,heterogeneity = "high"))
pow.vals03<-lapply(k,function(k) power.analysis.random(d=0.30,k=k,n1=25,n2=25,p=0.05,heterogeneity = "high"))
pow.vals01<-as.numeric(pow.vals01)
pow.vals02<-as.numeric(pow.vals02)
pow.vals03<-as.numeric(pow.vals03)
data3<-data.frame(k,pow.vals01,pow.vals02,pow.vals03)


ggplot()+
  geom_line(data = data1, aes(x = k, y = pow.vals01), color = "blue",size=2) +
  geom_line(data = data1, aes(x = k, y = pow.vals02), color = "#ff0000",size=2) +
  geom_line(data = data1, aes(x = k, y = pow.vals03), color = "#00cc00",size=2) +
  geom_line(data = data2, aes(x = k, y = pow.vals01), color = "lightblue",size=2) +
  geom_line(data = data2, aes(x = k, y = pow.vals02), color = "#ff4d4d",size=2) +
  geom_line(data = data2, aes(x = k, y = pow.vals03), color = "#1aff1a",size=2) +
  geom_line(data = data3, aes(x = k, y = pow.vals01), color = "darkblue",size=2) +
  geom_line(data = data3, aes(x = k, y = pow.vals02), color = "#b30000",size=2) +
  geom_line(data = data3, aes(x = k, y = pow.vals03), color = "#008000",size=2) +
  xlab('Number of Studies') +
  ylab('Power')+
  scale_y_continuous(labels = scales::percent)+
  theme(
        axis.line= element_line(color = "black",size = 1,linetype = "solid"),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.background = element_rect(linetype="solid",
                                         colour ="black"),
        legend.title = element_blank(),
        legend.key.size = unit(0.75,"cm"),
        legend.text=element_text(size=14))+
annotate("text", x = 3, y = 0.7, label = "d = 0.3",size=5)+
  annotate("text", x = 25, y = 0.6, label = "d = 0.2",size=5)+
  annotate("text", x = 20, y = 0.13, label = "d = 0.1",size=5)+
  geom_hline(yintercept=0.8,linetype="dashed")
```

---

## Power Calculator Tool

If you're feeling lazy, or if you want to quickly **check for the power of your meta-analysis under varying assumptions**, you might need a tool which makes it easier for you to calculate the power without having to run the R functions we described before each time.

We therefore built a online **Power Calculator Tool**, which you can find below. The calculations are based on the formulae and functions we described in the previous chapters.

[View in full page mode](https://mathiasharrer.shinyapps.io/power_calculator_meta_analysis/)

```{r,echo=FALSE, eval = FALSE}
knitr::include_app("https://mathiasharrer.shinyapps.io/power_calculator_meta_analysis/", height = "1500px")
```

---