library(dplyr)
library(tidyr)
library(ggplot2)
library(purrr)
library(lme4)
library(lmerTest)
library(broom.mixed)
library(faux)

This paper provides a tutorial and examples for more complex designs.

Building multilevel designs

You can build up a mixed effects model by adding random factors in a stepwise fashion using add_random(). The example below simulates 3 schools with 2, 3 and 4 classes nested in each school.

data <- add_random(school = 3) %>%
add_random(class = c(2, 3, 4), .nested_in = "school")
school class
s1 c1
s1 c2
s2 c3
s2 c4
s2 c5
s3 c6
s3 c7
s3 c8
s3 c9

A cross-classified design with 2 subjects and 3 items.

data <- add_random(subj = 2, item = 3)
subj item
s1 i1
s1 i2
s1 i3
s2 i1
s2 i2
s2 i3

Add within factors with add_within().

data <- add_random(subj = 2) %>%
condition = c("control", "test"))
subj time condition
s1 pre control
s1 pre test
s1 post control
s1 post test
s2 pre control
s2 pre test
s2 post control
s2 post test

Add between factors with add_between(). If you have more than one factor, they will be crossed. If you set .shuffle = TRUE, the factors will be added randomly (not in “order”), and separately added factors may end up confounded. This simulates true random allocation.

data <- add_random(subj = 4, item = 2) %>%
add_between("subj", cond = c("control", "test"), gen = c("X", "Z")) %>%
add_between("item", version = c("A", "B")) %>%
add_between(c("subj", "item"), shuffled = 1:4, .shuffle = TRUE) %>%
add_between(c("subj", "item"), not_shuffled = 1:4, .shuffle = FALSE)
subj item cond gen version shuffled not_shuffled
s1 i1 control X A 2 1
s1 i2 control X B 3 2
s2 i1 control Z A 2 3
s2 i2 control Z B 1 4
s3 i1 test X A 3 1
s3 i2 test X B 4 2
s4 i1 test Z A 4 3
s4 i2 test Z B 1 4

Unequal cells

If the levels of a factor don’t have equal probability, set the probability with .prob. If the sum of the values equals the number of groups, you will always get those exact numbers (e.g., always 8 control and 2 test; set shuffle = TRUE to get them in a random order). Otherwise, the values are sampled and the exact proportions will change for each simulation.

data <- add_random(subj = 10) %>%
add_between("subj", cond = c("control", "test"), .prob = c(8, 2)) %>%
add_between("subj", age = c("young", "old"), .prob = c(.8, .2))
subj cond age
s01 control old
s02 control young
s03 control old
s04 control young
s05 control young
s06 control young
s07 control young
s08 control old
s09 test young
s10 test young

If you have more than one between-subject factor, setting them together allows you to set joint proportions for each cell.

data <- add_random(subj = 10) %>%
cond = c("control", "test"),
age = c("young", "old"),
.prob = c(2, 3, 4, 1))
subj cond age
s01 control young
s02 control young
s03 control old
s04 control old
s05 control old
s06 test young
s07 test young
s08 test young
s09 test young
s10 test old

Counterbalanced designs

You can use the filter() function from dplyr to create counterbalanced designs. For example, if your items are grouped into A and B counterbalanced groups, add a between-subjects and a between-items factor with the same levels and filter only rows where the subject value matches the item value. In the example below, odd-numbered subjects only respond to odd-numbered items.

data <- add_random(subj = 4, item = 4) %>%
add_between("subj", subj_cb = c("odd", "even")) %>%
add_between("item", item_cb = c("odd", "even")) %>%
dplyr::filter(subj_cb == item_cb)
subj item subj_cb item_cb
s1 i1 odd odd
s1 i3 odd odd
s2 i2 even even
s2 i4 even even
s3 i1 odd odd
s3 i3 odd odd
s4 i2 even even
s4 i4 even even

Recoding

To set up data for analysis, you often need to recode categorical variables. Use the helper function add_contrast() for this. The code below creates anova-coded and treatment-coded versions of “cond” and the two variables needed to treatment-code the 3-level variable “type”. See the contrasts vignette for more details.

You can add_recode() for setting a manual contrast that isn’t available in add_contrast(), such as weighted contrasts.

data <- add_random(subj = 2, item = 3) %>%
add_between("subj", cond = c("A", "B")) %>%
add_between("item", type = c("X", "Y", "Z")) %>%
# you can change the default column names
add_recode("type", "type.w", X = -1.3, Y = 0, Z = 0.92)
subj item cond type cond.B-A cond.B-A.tr Y Z type.w
s1 i1 A X -0.5 0 0 0 -1.30
s1 i2 A Y -0.5 0 1 0 0.00
s1 i3 A Z -0.5 0 0 1 0.92
s2 i1 B X 0.5 1 0 0 -1.30
s2 i2 B Y 0.5 1 1 0 0.00
s2 i3 B Z 0.5 1 0 1 0.92

To simulate multilevel data, you need to add random intercepts and slopes for each random factor (or combination of random factors). These are randomly sampled each time you simulate a new sample, so you can only characterise them by their standard deviation. You can use the function add_ranef() to add random effects with specified SDs. If you add more than one random effect for the same group (e.g., a random intercept and a random slope), you can specify their correlation with .cors (you can specify correlations in the same way as for rnorm_multi()).

The code below sets up a simple cross-classified design where 2 subjects are crossed with 2 items. It adds a between-subject, within-item factor of “version”. Then it adds random effects by subject, item, and their interaction, as well as an error term (sigma).

data <- add_random(subj = 4, item = 2) %>%
# add by-item random intercept and slope
add_ranef("item", u0i = 1.5, u1i = 1.5, .cors = 0.3) %>%
add_ranef(c("subj", "item"), u0si = 1.7) %>%
# add error term (by observation)
add_ranef(sigma = 2.2)
subj item version u0s u0i u1i u0si sigma
s1 i1 1 1.175 -2.549 0.952 -1.531 1.503
s1 i2 1 1.175 -0.811 -1.595 -0.256 1.518
s2 i1 2 -2.014 -2.549 0.952 -1.407 1.174
s2 i2 2 -2.014 -0.811 -1.595 3.376 -0.409
s3 i1 1 1.329 -2.549 0.952 0.075 0.842
s3 i2 1 1.329 -0.811 -1.595 -0.687 0.828
s4 i1 2 0.195 -2.549 0.952 -0.804 2.538
s4 i2 2 0.195 -0.811 -1.595 -0.705 3.465

Simulating data

Now you can define your data-generating parameters and put everything together to simulate a dataset. In this example,subject are crossed with items, and there is a single treatment-coded between-subject, within-item fixed factor of condition with levels “control” and “test”. The intercept and effect of condition are both set to 0. The SDs of the random intercepts and slopes are all set to 1 (you will need pilot data to estimate realistic values for your design), the correlation between the random intercept and slope by items is set to 0. The SD of the error term is set to 2.

# define parameters
subj_n = 10  # number of subjects
item_n = 10  # number of items
b0 = 0       # intercept
b1 = 0       # fixed effect of condition
u0s_sd = 1   # random intercept SD for subjects
u0i_sd = 1   # random intercept SD for items
u1i_sd = 1   # random b1 slope SD for items
r01i = 0     # correlation between random effects 0 and 1 for items
sigma_sd = 2 # error SD

# set up data structure
data <- add_random(subj = subj_n, item = item_n) %>%
# add and recode categorical variables
add_between("subj", cond = c("control", "test")) %>%
add_recode("cond", "cond.t", control = 0, test = 1) %>%
add_ranef("item", u0i = u0i_sd, u1i = u1i_sd, .cors = r01i) %>%
# calculate DV
mutate(dv = b0 + u0s + u0i + (b1 + u1i) * cond.t + sigma)
subj item cond cond.t u0s u0i u1i sigma dv
s01 i01 control 0 0.589 2.135 0.684 -2.351 0.373
s01 i02 control 0 0.589 -0.214 -0.252 -0.284 0.090
s01 i03 control 0 0.589 0.947 -0.151 1.802 3.337
s01 i04 control 0 0.589 0.652 -0.875 0.716 1.957
s01 i05 control 0 0.589 -0.729 -1.977 -1.837 -1.977
s01 i06 control 0 0.589 -1.050 0.238 -2.644 -3.105

You can analyse these data with lme4::lmer().

m <- lmer(dv ~ cond.t + (1 | subj) + (1 + cond.t | item), data = data)

summary(m)
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: dv ~ cond.t + (1 | subj) + (1 + cond.t | item)
#>    Data: data
#>
#> REML criterion at convergence: 421
#>
#> Scaled residuals:
#>      Min       1Q   Median       3Q      Max
#> -1.94554 -0.66137  0.03661  0.66655  2.39501
#>
#> Random effects:
#>  Groups   Name        Variance Std.Dev. Corr
#>  subj     (Intercept) 0.3239   0.5691
#>  item     (Intercept) 0.7649   0.8746
#>           cond.t      1.1253   1.0608   -0.66
#>  Residual             3.2817   1.8115
#> Number of obs: 100, groups:  subj, 10; item, 10
#>
#> Fixed effects:
#>             Estimate Std. Error       df t value Pr(>|t|)
#> (Intercept) -0.05392    0.45486  9.66296  -0.119    0.908
#> cond.t      -0.58770    0.61102  9.08258  -0.962    0.361
#>
#> Correlation of Fixed Effects:
#>        (Intr)
#> cond.t -0.690

Power simulation

Include this code in a function so you can easily change the Ns, fixed effects, and random effects to any values. Then run a mixed effects model on the data and return the model.

sim <- function(subj_n = 10, item_n = 10,
b0 = 0, b1 = 0,         # fixed effects
u0s_sd = 1, u0i_sd = 1, # random intercepts
u1i_sd = 1, r01i = 0,   # random slope and cor
sigma_sd = 2,           # error term
... # helps the function work with pmap() below
) {

# set up data structure
data <- add_random(subj = subj_n, item = item_n) %>%
# add and recode categorical variables
add_between("subj", cond = c("control", "test")) %>%
add_recode("cond", "cond.t", control = 0, test = 1) %>%
add_ranef("item", u0i = u0i_sd, u1i = u1i_sd, .cors = r01i) %>%
# calculate DV
mutate(dv = b0 + u0s + u0i + (b1 + u1i) * cond.t + sigma)

# run mixed effect model and return relevant values
m <- lmer(dv ~ cond.t + (1 | subj) + (1 + cond.t | item), data = data)

broom.mixed::tidy(m)
}

Check the function. Here, we’re simulating a fixed effect of condition (b1) of 0.5 for 50 subjects and 40 items, with a correlation between the random intercept and slope for items of 0.2, and the default values for all other parameters that we set above.

sim(subj_n = 50, item_n = 40, b1 = 0.5, r01i = 0.2)
#> # A tibble: 7 x 8
#>   effect   group   term               estimate std.error statistic    df p.value
#>   <chr>    <chr>   <chr>                 <dbl>     <dbl>     <dbl> <dbl>   <dbl>
#> 1 fixed    <NA>    (Intercept)         -0.0715     0.286    -0.250  71.4   0.803
#> 2 fixed    <NA>    cond.t               0.499      0.395     1.27   68.1   0.210
#> 3 ran_pars subj    sd__(Intercept)      1.19      NA        NA      NA    NA
#> 4 ran_pars item    sd__(Intercept)      0.913     NA        NA      NA    NA
#> 5 ran_pars item    cor__(Intercept).…   0.256     NA        NA      NA    NA
#> 6 ran_pars item    sd__cond.t           1.16      NA        NA      NA    NA
#> 7 ran_pars Residu… sd__Observation      2.08      NA        NA      NA    NA

Run the simulation repeatedly. I’m only running it 50 times per parameter combination here so my demo doesn’t take forever to run, but once you are done testing a range of parameters, you probably want to run the final simulation 100-1000 times. You might get a few warnings like “Model failed to converge” or “singular boundary”. You don’t need to worry too much if you only get a few of these. If most of your models have warnings, the simulation parameters are likely to be off.

x <- crossing(
rep = 1:50, # number of replicates
subj_n = c(50, 100), # range of subject N
item_n = 25, # fixed item N
b1 = c(0.25, 0.5, 0.75), # range of effects
r01i = 0.2 # fixed correlation
) %>%
mutate(analysis = pmap(., sim)) %>%
unnest(analysis)

Filter and/or group the resulting table and calculate the proportion of p.values above your alpha threshold to get power for the fixed effects.

# calculate power for alpha = 0.05
filter(x, effect == "fixed", term == "cond.t") %>%
group_by(b1, subj_n) %>%
summarise(power = mean(p.value < .05),
.groups = "drop") %>%
ggplot(aes(b1, subj_n, fill = power)) +
geom_tile() +
geom_text(aes(label = sprintf("%.2f", power)), color = "white", size = 10) +
scale_fill_viridis_c(limits = c(0, 1))