library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.1 ✓ dplyr 1.0.6
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(simr)## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4##
## Attaching package: 'simr'## The following object is masked from 'package:lme4':
##
## getData## The following object is masked from 'package:stringr':
##
## fixedlibrary(lme4)
# Create data
subj <- factor(1:20) # 20 students per class
class_id <- factor(1:50) # 50 classes
time <- 0:1 # 2 time points
group <- c("control", "intervention") # two groups
subj_full <- rep(subj, 15)
class_full <- rep(rep(class_id, each=10), 3)
time_full <- rep(time, each=50)
group_full <- rep(rep(group, each=5), 15)
covars <- data.frame(id=subj_full, class=class_full, treat=group_full, time=factor(time_full))
covars <- as_tibble(covars)
covars## # A tibble: 1,500 x 4
## id class treat time
## <fct> <fct> <chr> <fct>
## 1 1 1 control 0
## 2 2 1 control 0
## 3 3 1 control 0
## 4 4 1 control 0
## 5 5 1 control 0
## 6 6 1 intervention 0
## 7 7 1 intervention 0
## 8 8 1 intervention 0
## 9 9 1 intervention 0
## 10 10 1 intervention 0
## # … with 1,490 more rows# Set betas and other parameters
fixed <- c(0, 0, 0.05, 0.20) # Intercept and slopes for time and intervention and time:intervention)
rand <- list(0.16) # equivalent to an ICC of 20%
res <- .80 # residual variance
# model
model <- makeLmer(y ~ treat*time + (1|class), fixef=fixed, VarCorr=rand, sigma=res, data=covars)
model## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ treat * time + (1 | class)
## Data: covars
## REML criterion at convergence: 3690.646
## Random effects:
## Groups Name Std.Dev.
## class (Intercept) 0.4
## Residual 0.8
## Number of obs: 1500, groups: class, 50
## Fixed Effects:
## (Intercept) treatintervention time1
## 0.00 0.00 0.05
## treatintervention:time1
## 0.20# simulate!
sim_treat <- powerSim(model, nsim=100, test = fcompare(y~time))## Simulating: | |Simulating: |= |Simulating: |== |Simulating: |=== |Simulating: |==== |Simulating: |===== |Simulating: |====== |Simulating: |======= |Simulating: |======== |Simulating: |========= |Simulating: |========== |Simulating: |=========== |Simulating: |============ |Simulating: |============= |Simulating: |============== |Simulating: |=============== |Simulating: |================ |Simulating: |================= |Simulating: |================== |Simulating: |=================== |Simulating: |==================== |Simulating: |===================== |Simulating: |====================== |Simulating: |======================= |Simulating: |======================== |Simulating: |========================= |Simulating: |========================== |Simulating: |=========================== |Simulating: |============================ |Simulating: |============================= |Simulating: |============================== |Simulating: |=============================== |Simulating: |================================ |Simulating: |================================= |Simulating: |================================== |Simulating: |=================================== |Simulating: |==================================== |Simulating: |===================================== |Simulating: |====================================== |Simulating: |======================================= |Simulating: |======================================== |Simulating: |========================================= |Simulating: |========================================== |Simulating: |=========================================== |Simulating: |============================================ |Simulating: |============================================= |Simulating: |============================================== |Simulating: |=============================================== |Simulating: |================================================ |Simulating: |================================================= |Simulating: |================================================== |Simulating: |=================================================== |Simulating: |==================================================== |Simulating: |===================================================== |Simulating: |====================================================== |Simulating: |======================================================= |Simulating: |======================================================== |Simulating: |========================================================= |Simulating: |========================================================== |Simulating: |=========================================================== |Simulating: |============================================================ |Simulating: |============================================================= |Simulating: |============================================================== |Simulating: |=============================================================== |Simulating: |================================================================ |Simulating: |================================================================= |Simulating: |==================================================================|sim_treat## Power for model comparison, (95% confidence interval):
## 90.00% (82.38, 95.10)
##
## Test: Likelihood ratio
## Comparison to y ~ time + [re]
##
## Based on 100 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 1500
##
## Time elapsed: 0 h 0 m 12 sresources