setwd("C:/Work Files/Collaboration/Andi")Andi R&R power analysis
I am using two packages to create a simulation in which to test the power of a 2-level model with level-1: students (n=450) and level=2: schools (n=10). The simulated model assumes approximately 45 students per school.
library(simr)Warning: package 'simr' was built under R version 4.4.2Loading required package: lme4Warning: package 'lme4' was built under R version 4.4.2Loading required package: Matrix
Attaching package: 'simr'The following object is masked from 'package:lme4':
getDatalibrary(lme4)
library(lmerTest)
Attaching package: 'lmerTest'The following object is masked from 'package:lme4':
lmerThe following object is masked from 'package:stats':
stepFirst I am generating a simulated data set that will serve as the population model for the Monte Carlo power simulation. This data set has 1-outcome, 1-student predictor, and 1-school level effect, all of which are on a standard scale ~N(0,1). I chose a moderate student level effect of \(\beta = 0.3\). I can also change the random effect standard deviations to make them narrower or wider depending on how you want to model the precision of the estimates.
set.seed(123) # Reproducibility
# Parameters
n_students <- 450 # Level 1 sample size
n_schools <- 10 # Level 2 sample size
students_per_school <- n_students / n_schools # Average students per school
# Generate school-level (random intercepts)
school_id <- rep(1:n_schools, each = students_per_school)
# Level-1 covariates (e.g., student scores)
student_covariate <- rnorm(n_students, mean = 0, sd = 1) # Predictor
# Fixed effect of student covariate (moderate effect size ~ 0.3)
beta <- 0.3
# Random effects (school-level variability)
school_effect <- rnorm(n_schools, mean = 0, sd = .8)
random_intercept <- school_effect[school_id]
# Generate outcome variable (Y)
epsilon <- rnorm(n_students, mean = 0, sd = 1) # Residual error
y <- 2 + random_intercept + beta * student_covariate + epsilon
# Combine into a data frame
data <- data.frame(
school_id = factor(school_id),
student_covariate = student_covariate,
y = y
)This block of code provides a summary of the simulated data confirming the population model. The summary(model) output uses the Satterthwaite method for d.f. and the anova output uses the Kenward-Roger adjusted d.f.
# Model: Random intercept for schools
model <- lmer(y ~ student_covariate + (1 | school_id), data = data)
# Summarize the model
summary(model)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: y ~ student_covariate + (1 | school_id)
Data: data
REML criterion at convergence: 1313
Scaled residuals:
Min 1Q Median 3Q Max
-2.66062 -0.65156 0.02398 0.66852 2.70787
Random effects:
Groups Name Variance Std.Dev.
school_id (Intercept) 0.9938 0.9969
Residual 0.9889 0.9944
Number of obs: 450, groups: school_id, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.22712 0.31871 8.99964 6.988 6.41e-05 ***
student_covariate 0.27133 0.04827 439.10443 5.621 3.38e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
studnt_cvrt -0.001anova(model, type=2, ddf = "Kenward-Roger")Type II Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
student_covariate 31.24 31.24 1 439.1 31.592 3.386e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This last block of code uses the coefficients from the above simulated population data to estimate the power across a 1,000 random draws from the population data. There is ample power here
# Extend the model for power simulation
model_ext <- extend(model, along = "school_id", n = n_schools)
# Power analysis for the random effect at the school level
power_results_school <- powerSim(model_ext, test = ranef, nsim = 1000)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: |==================================================================|Warning in observedPowerWarning(sim): This appears to be an "observed power"
calculation# Power analysis for the fixed effect of student_covariate
power_results_students <- powerSim(model_ext, test = fixed("student_covariate"), nsim = 1000)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: |==================================================================|Warning in observedPowerWarning(sim): This appears to be an "observed power"
calculation# Display results
print(power_results_school)Power [user defined], (95% confidence interval):
0.00% ( 0.00, 0.37)
Test: [user defined function]
Based on 1000 simulations, (0 warnings, 1000 errors)
alpha = 0.05, nrow = 450
Time elapsed: 0 h 0 m 49 s
nb: result might be an observed power calculationprint(power_results_students)Power for predictor 'student_covariate', (95% confidence interval):
99.90% (99.44, 100.00)
Test: unknown test
Effect size for student_covariate is 0.27
Based on 1000 simulations, (0 warnings, 0 errors)
alpha = 0.05, nrow = 450
Time elapsed: 0 h 1 m 51 s
nb: result might be an observed power calculation