Scenario1_notebook.Rmd

---
title: 'Scenario 1: Using an available well-powered design as a starting point'
output:
  html_document:
    df_print: paged
  pdf_document: default
---
Kumle, L., Vo, M. L-H., & Draschkow, D.

latest update: May 2021

 

This notebook aims at accompanying Scenario 1 in Kumle, Vo & Draschkow [(2021)](https://doi.org/10.3758/s13428-021-01546-0).  
Additionally to giving a brief general introduction to simulation-based power analyses, this notebook focuses on using power analysis as a tool for a priori sample size estimation. To do so, the packages [*mixedpower* ](https://github.com/DejanDraschkow/mixedpower) (Kumle, Vo & Draschkow, 2018) and [*simr*](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504) (Green & Macleod, 2016) are used and all simulations are based on (generalized) mixed models fitted with [*lme4*](https://arxiv.org/pdf/1406.5823.pdf) (Bates, Maechler, Bolker & Walker, 2015). All simulations shown in this notebook are applicable for linear mixed models as well as generalized linear mixed models.

 

See [Scenario 2](https://lkumle.github.io/power_notebooks/Scenario2_notebook.html) for power analyses varying factors other than *participants* and [Scenario 3](https://lkumle.github.io/power_notebooks/Scenario3_notebook.html) for starting a simulation-based power estimation from scratch.
 

***
***
## Introduction to simulation-based power analysis

Being able to estimate power is important as it is closely linked to the reliability and replicability of empirical findings. Classical solutions to power analysis work with analytical formulas. However, (generalized) linear mixed models ((G)LMMs) are often too complex to be solved analytically and therefore require a different approach. A flexible and more intuitive alternative to analytic power solutions are simulation-based power analyses (Brysbaert & Stevens, 2018; Thomas & Juanes, 1996). In simple terms, one basic question behind power analyses is: “Suppose there really is an effect of a certain size and I run my experiment one hundred times - how many times will I get a statistically significant result?” (Coppock, 2013). As it is possible to simulate the outcome of an experiment, power can be calculated based on the proportion of significant simulations to all simulations (Johnson et al., 2015; Thomas & Juanes, 1996).
 

The shared principle of all simulation-based power analyses solutions can therefore be broken down into the following steps:     
**1)** simulation of new data sets,   
**2)** analysis of each data set and test for statistical significance, and   
**3)** calculation of the proportion of significant to all simulations.  

![](/Users/leah/Dropbox/Power/manuscript/BRM/notebooks/Scenario_1/simulation_based_power.png)


However, accuracy of the power estimate strongly depends on the accuracy of our simulation - informing the parameters of the simulation therefore is a critical step. A more detailed discussion can be found in Kumle, Vo & Draschkow [(2021)](https://doi.org/10.3758/s13428-021-01546-0).


***
***

### Import Data

All simulations in this notebook will be based on an available and preceding well-powered design. This is possibly the most desired solution since we can utilize a (G)LMM fitted on real and independent empirical data to inform the simulation parameters.
   
Data used for this tutorial originates from Yan, Zhou, Shu, Yusupu, Miao, Kruegel & Kliegl (2014). All analyses and data were obtained from http://read.psych.uni-potsdam.de where the authors made them available.  
Yan et all. (2014) tested 48 subjects, each of whom read 120 sentences investigating the effect of different factors on the first landing position (FLP, i.e. the position in a sentence your eyes first land on).


Suppose the goal is to conduct a study replicating and further investigating the effect of morphological complexity (i.e. number of suffixes) and word length – we would need to conduct a power analysis in order to inform the sample size of the follow-up study.

To focus on the actual power analysis procedure, we already preprocessed the data set according to the author's analysis and the data can be downloaded [here](https://github.com/lkumle/analyses_power_tutorial/blob/master/Scenario%201%20%26%202/Yan_et_al.Rdata). Moreover, several variables have been renamed for more clarity.
```{r include = F}
# get data
load("~/Dropbox/Power/manuscript/BRM/analyses_BRM/Yan_et_al.RData") # data set is called "YanData"
data <- YanData
data["word_length"] <- data$wl.c
data["complexity"] <- data$sn.c
data["subject"] <- data$nsub
data["sentence"] <- data$nsen
YanData <- data
```

```{r eval = F}
# get data
load("Yan_et_al.RData") # data set is called "YanData"
```

&nbsp;  


First, we need a corresponding model fitted with *lme4*.
```{r message=FALSE, warning=FALSE}
library(lme4)
#  LMM including word length and word complexity as fixed effects and random intercepts for subject and sentence
FLPmodel <- lmer(flp ~ word_length * complexity + (1|subject) + (1|sentence),
                 data = YanData)

summary(FLPmodel, corr = F) # let's have a look!
```
The FLPmodel includes word length ( ß = 1.511) and morphological complexity (ß = - 0.075) as well as their interaction (ß = 0.116) as fixed effects (see Table 1). Moreover, we included random intercepts for the random effects of subjects and sentence (i.e. stimuli) making this model a typical example with crossed random effects as described by Baayen et al. (2008).  

***
***
### Power analysis using mixedpower    

While it is possible to implement power analyses from scratch, using preprogrammed packaged saves us a lot of time and effort. To start with, we will make use of the [*mixedpower*](https://github.com/DejanDraschkow/mixedpower) package (Kumle, Vo & Draschkow, 2018). An overview of all functions included in mixedpower can be found in its [documentation](https://lkumle.github.io/power_notebooks/intro/Introduction_to_mixedpower.pdf) and more theoretical guidance on the specific parameters can be found in Kumle, Vo & Draschkow [(2021)](https://doi.org/10.3758/s13428-021-01546-0).

```{r eval=FALSE}
# install mixedpower
if (!require("devtools")) {
    install.packages("devtools", dependencies = TRUE)}
devtools::install_github("DejanDraschkow/mixedpower") # mixedpower is hosted on GitHub

# load library
library(mixedpower)
```
&nbsp;  

#### **Varying sample size**
Since we want to explore power for different sample sizes, which corresponds to one of our random variables (i.e. subject), we will make use of the mixedpower() function.

As can be seen in the [documentation](https://lkumle.github.io/power_notebooks/intro/Introduction_to_mixedpower.pdf), using mixedpower() requires to specify various parameters. To get a quick overview of those parameters, *?mixedpower* can be used as well. While those parameters can be directly specified as arguments in the function, we chose to make it more explicit and specify them beforehand by assigning them to variables.   


First, we will provide some general information about the model we have chosen to simulate power for.
```{r eval=FALSE}

# ------------------------------------------ #
# INFORMATION ABOUT MODEL USED FOR SIMULATION

model <- FLPmodel # which model do we want to simulate power for?
data <- YanData # data used to fit the model
fixed_effects <- c("word_length", "complexity") # all fixed effects specified in FLPmodel
simvar <- "subject" # which random effect do we want to vary in the simulation?

```
&nbsp;  

Next, we will set the parameters determining the details of our simulation.
```{r eval=FALSE}

# ------------------------------------------ #
# SIMULATION PARAMETERS
steps <- c(20, 30, 40, 50, 60) # which sample sizes do we want to look at?
critical_value <- 2 # which t/z value do we want to use to test for significance?
n_sim <- 1000 # how many single simulations should be used to estimate power?
```
&nbsp;  

Finally, we will combine this information in the mixedpower()-function and run the power analysis. This might take a while (depending on the complexity and size of the model as well as on the number of cores on the machine used for simulation)
… time to take a short break! Running on a machine with 6 cores, the following simulation will take approximately 10-15 minutes …
```{r eval=FALSE}

# ------------------------------------------ #
# RUN SIMULATION
power_FLP <- mixedpower(model = FLPmodel, data = YanData,
                        fixed_effects = c("word_length", "complexity"),
                        simvar = "subject", steps = c(20,30,40,50,60),
                        critical_value = 2, n_sim = 1000)

```
&nbsp;   
```{r include=FALSE}
# load output in the background
load("~/Dropbox/Power/manuscript/BRM/analyses_BRM/Scenario 1/mixedpower_scenario1_SESOI.RData")
power_FLP <- power_FLP[1:3,]
```



```{r}
# let's have a first look at the results:
power_FLP
```
&nbsp;  
&nbsp;   

#### **SESOI**

Comparing the parameters we set above with the documentation or help page of mixedpower reveals that we skipped two parameters - namely *SESOI = F* and *databased = T*.
Both parameters determine the effect sizes used for simulation and the default is to include a *databased* simulation (i.e. using the beta coefficients found in the specified model) and no *SESOI* simulation (i.e. specifiying a **Smalles effect of interest** which is used to simulate power). However, we can provide a vector with new beta coefficients and hand it to mixedpower() using the SESOI argument. Theoretical guidance on this issue can be found in the accompanying tutorial paper (Kumle, Vo & Draschkow [(2021)](https://doi.org/10.3758/s13428-021-01546-0)).

So let’s include a SESOI simulation as well:

```{r eval=FALSE}

# ------------------------------------------ #
# INCLUDE SESOI SIMULATION
SESOI <- c(3.66, 0.75, -0.065, 0.09) # specify SESOI

power_SESOI <- mixedpower(model = FLPmodel, data = YanData,
                              fixed_effects = c("word_length", "complexity"),
                              simvar = "subject", steps = c(20,30,40,50,60),
                              critical_value = 2, n_sim = 1000,
                              SESOI = SESOI, databased = T)

```
Note that we set *databased = T* again. Since we already estimated databased power above, we could save some time and not run it again this time (this will result it slightly different estimates due to the nature of our simulations). However, we will run it again to later plot both estimations at once.


```{r include=FALSE}
# load output in the background
load("~/Dropbox/Power/manuscript/BRM/analyses_BRM/Scenario 1/mixedpower_scenario1_SESOI.RData")
power_SESOI <- power_FLP
```



```{r}
# let's have a first look at the results:
power_SESOI
```
&nbsp;  
&nbsp;  

#### **Plot and interpret results**

Being able to interpret and combine the results of power simulations is important for a power analysis to be a useful assistant in experimental planning. Before we try to interpret the results of our power analysis, we will plot it using the plotting function **multiplotPower()** included in the mixedpower package.
```{r eval=FALSE}

# ------------------------------------------ #
# PLOT RESULTS
multiplotPower(power_SESOI)

```

![](/Users/leah/Dropbox/Power/manuscript/BRM/analyses_BRM/Scenario 1/plot_mixedpower_S1.png)

Finally, we need to transfer the plots and results into a decision regarding the sample size of our planned study. General guidance can be found in the accompanying tutorial paper. In this notebook, we will adapt a strategy resulting in > 80% power for all effects included in the model. Using the *databased* estimate, we would agree to test 50 participants. However, using the *SESOI* estimate, we would test more than 60 participants as we do not reach a power of >80% for all effects with a sample size of 60. Moreover, additional SESOI-simulations with larger sample sizes would be necessary to reach a sample size that fulfils our power criterion.

***
***


### Power analysis using simr  
&nbsp;  

There are other resources to perform simulation-based power analyses, one of them being [**simr**](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504) (Green & Macleod, 2016). The tutorial mainly focusses on mixedpower, as this package has advantages in computational speed and allows for more factors to be varied simultaneously. However, simr provides some exceptionally easy to implement functions that are worth looking at.
&nbsp;  

Similar to the mixedpower()-function, simr provides a function called *powerCurve()* which allows to explore power for different sample sizes. However, all functions included in simr simulate power for one effect at a time - we therefore need to specify which effect we want to test.
To get a first overview, we would like to explore power for 20, 30, 40 and 60 participants. As this includes a sample size larger than in the data set used for simulation, we need to *extend* the data set before we hand it to the simulation.

```{r eval=FALSE}
library(simr)
# -------------------------------------------- #
# EXTEND DATA SET
FLPmodel_extended <- extend(FLPmodel, along="subject", n = 60) #extend data set to include 60 participants
```

All functions included in simr simulate power for one effect at a time - we therefore need to specify which effect we want to test. Moreove, since computational speed in simr is limited for more complex models and data, we will only include 100 simulations.

```{r eval=FALSE}
library(simr)
# -------------------------------------------- #
# POWER ANALYSIS WITH SIMR
powerC <- powerCurve(fit = FLPmodel_extended, test = fixed("complexity"), along = "subject",
                      breaks = c(20, 30, 40, 60))
```

```{r include=FALSE}
# load output in the background
load("~/Dropbox/Power/manuscript/BRM/analyses_BRM/Scenario 1/simr_powerC_scenario1_snc.RData")
library(simr)
```

```{r}
# let's have a look at the results
print(powerC)
```
We see that results are similar to the results achieved with mixedpower although they differ slightly. This is to be expected, as the simulations in simr are based on only 100 repetitions and therefore are more prone to variation.
&nbsp;   

#### **SESOI in simr**  

Simr also comes with an option to change the effect sizes in accordance with a SESOI approach. All that has to be done is to modify the effect sizes in the model before handing it to the powerCurve()-function.

```{r eval=FALSE}

# -------------------------------------------- #
# SESOI WITH SIMR
fixef(FLPmodel)["complexity"] <- -0.065
```

After this step, we can just repeat the power analysis shown above.   

***
***
Continue with [Notebook 2](https://lkumle.github.io/power_notebooks/Scenario2_notebook.html) for power analyses varying factors other than *sample size* and [Notebook 3](https://lkumle.github.io/power_notebooks/Scenario3_notebook.html) for starting a simulation-based power estimation from scratch.   

***
***
### References     

Baayen, R. H., Davidson, D. J., & Bates, D. M. (2008). Mixed-effects modeling with crossed random effects for subjects and items. Journal of Memory and Language, 59(4), 390–412. https://doi.org/10.1016/j.jml.2007.12.005

Bates, D., Maechler, M., Bolker, B., & Walker, S. (2014). Fitting Linear Mixed-Effects Models using lme4. Journal of Statistical Software, 67(1). https://doi.org/10.18637/jss.v067.i01

Brysbaert, M., & Stevens, M. (2018). Power Analysis and Effect Size in Mixed Effects Models: A Tutorial. Journal of Cognition, 1(1), 1???20. https://doi.org/10.5334/joc.10

Coppock, A. (2013). 10 Things to Know About Statistical Power. Retrieved September 20,
2018, from http://egap.org/methods-guides/10-things-you-need-know-about-statistical-power

Green, P., & Macleod, C. J. (2016). SIMR: An R package for power analysis of generalized linear mixed models by simulation. Methods in Ecology and Evolution, 7(4), 493-498. [https://doi.org/10.1111/2041-210X.12504](https://doi.org/10.1111/2041-210X.12504)

Johnson, P. C. D., Barry, S. J. E., Ferguson, H. M., & Müller, P. (2015). Power analysis for
generalized linear mixed models in ecology and evolution. Methods in Ecology and
Evolution, 6(2), 133–142. https://doi.org/10.1111/2041-210X.12306

Kumle, L., Võ, M.LH. & Draschkow, D. Estimating power in (generalized) linear mixed models: An open introduction and tutorial in R. Behav Res (2021). https://doi.org/10.3758/s13428-021-01546-0

Kumle, L., Vo, M. L-H., & Draschkow, D. (in preparation). Estimating power in linear and generalized linear mixed models: an open introduction and tutorial in R.

Thomas, L., & Juanes, F. (1996). The importance of statistical power analysis: An example from Animal Behaviour. Animal Behaviour, 52(4), 856-859. https://doi.org/10.1006/anbe.1996.0232

Yan, M., Zhou, W., Shu, H., Yusupu, R., Miao, D., Kruegel, A., & Kliegl, R. (2014). Eye movements guided by morphological structure: Evidence from the Uighur language. Cognition, 132(2), 181???215. https://doi.org/10.1016/j.cognition.2014.03.008

***
***