forked from lkumle/power_notebooks
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathScenario2_notebook.Rmd
247 lines (180 loc) · 13 KB
/
Scenario2_notebook.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
---
title: 'Scenario 2: Simulating different units (random variables)'
output:
html_document:
df_print: paged
pdf_document: default
---
Kumle, L., Vo, M. L-H., & Draschkow, D.
latest update: March 2020
This notebook accompanies Scenario 2 in Kumle, Vo & Draschkow [(2021)](https://doi.org/10.3758/s13428-021-01546-0).
For a general introduction to simulation-based power analyses as well as using simulations to explore power for different sample sizes, see [Scenario 1](https://lkumle.github.io/power_notebooks/Scenario1_notebook.html).
In this notebook, we will cover scenarios in which we wish to simulate power for random variables other than subject. 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 3](https://lkumle.github.io/power_notebooks/Scenario3_notebook.html) for starting a simulation-based power estimation from scratch.
***
***
### Import Data
As in [Scenario 1](https://lkumle.github.io/power_notebooks/Scenario1_notebook.html), we will base our simulation on data published by Yan, Zhou, Shu, Yusupu, Miao, Kruegel & Kliegl (2014) retrieved from http://read.psych.uni-potsdam.de. A preprocessed data set can be downloaded [here](https://github.com/lkumle/analyses_power_tutorial/blob/master/Scenario%201%20%26%202/Yan_et_al.Rdata).
```{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"
```
Yan et all. (2014) tested 48 subjects, each of whom read 120 sentences investigating the effect of different factors on the first landing position while reading (FLP, i.e. the position in a sentence your eyes first land on) - a more detailed introduction to the data set can be found in [Scenario 1](https://lkumle.github.io/power_notebooks/Scenario1_notebook.html) where we varied the number of subjects. The levels of the second random parameter (i.e. stimuli/items), however, have been kept constant until now.
Again, we are interested in the effect of morphological complexity (i.e. number of suffixes) and word length and therefore will include them as fixed effects in our 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 corresponding 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).
***
***
### Varying the number of stimuli
This time, we are not interested in changing the number of subjects but the number of sentences around the original number of 120.
#### **mixedpower**
Using the mixedpower package, all that has to be changed compared to Scenario 1 is the specification of **simvar** and its **steps**. For clarity, we will again specify all parameters before running the actual power analysis.
```{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)
```
First, we will provide some general information about the model we have chosen to simulate power for. Note how we specified "sentence" as the random variable we want to vary.
```{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 <- "sentence" # which random variable do we want to vary in the simulation?
```
Next, we will set the parameters determining the details of our simulation. As the original data set contains 120 sentences, we will explore power for a range of sentences around that number.
Similar to Scenario 1, we will also include the SESOI (i.e. smallest effect of interest) specification. An introduction to the SESOI can be found in Notebook 1 and the tutorial paper Kumle, Vo & Draschkow [(2021)](https://doi.org/10.3758/s13428-021-01546-0).
```{r eval=FALSE}
# ------------------------------------------ #
# SIMULATION PARAMETERS
steps <- c(100, 120, 140, 160, 180) # 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?
# ------------------------------------------ #
# INCLUDE SESOI SIMULATION
SESOI <- c(3.66, 0.75, -0.065, 0.09) # specify SESOI
```
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_sentences <- mixedpower(model = FLPmodel, data = YanData,
fixed_effects = c("word_length", "complexity"),
simvar = "subject", steps = c(100, 120, 140, 160, 180),
critical_value = 2, n_sim = 1000,
SESOI = SESOI, databased = T)
```
```{r include=FALSE}
# load output in the background
load("~/Dropbox/Power/manuscript/BRM/analyses_BRM/Scenario 2/R2_scenario2_48subjects.Rdata")
power_sentences <- power48_sentences
```
Let's have a look at the results and plot them:
```{r}
# let's have a first look at the results:
power_sentences
```
```{r eval=FALSE}
# ------------------------------------------ #
# PLOT RESULTS
multiplotPower(power_sentences)
```
***
#### **simr**
The same analyses can be done with the [*simr*](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504) package (Green & Macleod, 2016) - with the difference that we have to test one effect at a time. Similar to the analysis with mixedpower, we just specify “sentence” as the variable we want to simulate *along* and extend the model along this random effect first.
```{r eval=FALSE}
library(simr)
library(simr)
# -------------------------------------------- #
# EXTEND DATA SET
FLPmodel_extended <- extend(FLPmodel, along="sentence", n = 180) #extend data set to include 180 sentences
# -------------------------------------------- #
# POWER ANALYSIS WITH SIMR
powerC_sentences <- powerCurve(fit = FLPmodel_extended, test = fixed("complexity"), along = "sentence",
breaks = c(100, 120, 140, 160, 180))
```
A more detailed tutorial for simr can be found [here](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504) and the inclusion of the simr package in this notebook simply aims at introducing the package to give a more complete overview of existing resources.
***
***
### Varying the levels multiple random variables simultaneously
So far, we only varied the levels of one random variable while we kept the other one constant. However, being able to vary both simultaneously would be desirable as it allows to estimate power for different combinations of levels.
To do so, we will make use of the **R2power()-function** in mixedpower which allows to do exactly that (check ?R2power or the [documentation](https://lkumle.github.io/power_notebooks/intro/Introduction_to_mixedpower.pdf) of mixedpower for more details).
The principle behind the simulations in R2power stays the same as in the other approaches in mixedpower introduced earlier. First, new data sets are simulated before they are used to refit the model entered into the simulation. Second, all fixed effects (and interactions) in the refitted models are checked for significance before power is estimated as the proportion of significant to all simulation.
However, the simulation of new data sets now follows a two-step simulation process where the levels of one random variable are simulated before the levels of the second random variable are adapted. To do so, we need to specify the desired levels and names of *both* random variables.
Imagine wanting to explore power for a range of sentences - but not for 48 subjects like in the data published by Yan et al. (2014) but for 30 subjects. Consequentially, we would keep sentences as our “simvar” and add subject as the second varying random variable. All other parameters will stay the same as in the simulation run with mixedpower() above.
```{r eval=FALSE}
# ------------------------------------------ #
# SPECIFY R2 PARAMETERS
R2var <- "subject" # specify "subject" as the second random effect we want to vary
R2level <- 30 # which level should "subject" have?
# ------------------------------------------ #
# RUN SIMULATION
R2power <- R2power(model = FLPmodel, data = YanData,
fixed_effects = c("word_length", "complexity"),
simvar = "sentence", steps = c(100,120,140,160,180),
R2var = "subject", R2level = 30, critical_value = 2,
n_sim = 1000, SESOI = SESOI, databased = T)
```
```{r include=FALSE}
# load output in the background
load("~/Dropbox/Power/manuscript/BRM/analyses_BRM/Scenario 2/R2power_Notebook2.Rdata")
R2power <- R2power_Notebook2
```
```{r}
# let's have a look at the results
R2power
```
```{r eval=FALSE}
# ------------------------------------------ #
# PLOT RESULTS
multiplotPower(R2power)
```
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. To reach this threshold, we could either test 48 subjects each of whom reads 160 sentences. In a scenario where we only have resources for 30 subjects even 180 sentences would not be enough to power our study with > 80% power for all effects if we base our decision on the SESOI effect (as can be seen in the results of the R2power simulation). Moreover, we could run more simulations for even more combinations and base our decision on them. Repeating this process with different combinations (e.g. different R2levels) enables an extensive overview of power and the factors that influence it, which can be especially useful in cases in which one of the random effects cannot be increased, i.e. a ceiling on the amount of participants or stimuli.
***
***
Continue with [Notebook 3](https://lkumle.github.io/power_notebooks/Scenario3_notebook.html) for starting a simulation-based power estimation from scratch.
***
***
### References
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
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)
Kumle, L., Vo, M. L-H., & Draschkow, D. (2018). Mixedpower: a library for
estimating simulation-based power for mixed models in R. https://doi.org/10.5281/zenodo.1341047
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
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
***