# Power analysis for simple LMER model in R from scratch

I'd like to compute power for some rather simple lmer models. I am aware of some packages like simr but the tutorials are not easy to follow because they demonstrate much more complicated cases than mine. When I try to adapt the cases to fit my models, I continually get results like 0.00% power.

So I'd like to compute power from scratch. In the past, when I've compute power for simple regression or t-tests, I've created a code that essentially:

1. simulates data
2. runs n tests
3. Sums up the number of p-values < .05 and divides by n tests.

However, no matter how I change the parameters, it always returns the same power values- usually between 68-70%.

I'm wondering: why is this happening? What have I done wrong or misunderstood? How can I correct this?

I adapted that code for a simple lmer model here (actually using lme because it was easier to save the p-vaues):

CORRECT CODE BASED ON FIX FROM COMMENTS

options(scipen = 999, digits = 3)
#Power simulation for mixed effects models

library(nlme)
library(dplyr)
library(tidyr)
library(purrr)
library(magrittr)

# 1. First, determine your parameters. The below numbers will create 4 dfs with different Ns
#(eg N=5, 10, 50, & 200) and cross each level with the other level of the other parameters specified.

params <- expand.grid(
N = c(20, 30, 50), #this is the N for each sample df
B1 = c(.05, .1, .5),
B2 = c(.3, .5, .7),#the effect size B, a linear coefficient of x predicting y
items = c(2, 5, 7)) #the number of repeating within-subject time points, items, etc
n_samples <- 50 #how many simulations you want to run

#2. Then, run this entire chunk with cmd + shift + enter

#Don't modify anything below
###############################################################
#create a sample
sim_sample <- function(N, B1, B2, items){
tid <- rep(1:items, N) #time id, not really needed but maybe useful
subj <- rep(1:N, each=items) #subject id
x = rnorm(n = items*N) #a fixed effect predictor
m = rnorm(n = items*N)
y <- B1*x + B2*m + rnorm(n = items*N) #y is defined as being predicted by x + random error
smaller.df <- data.frame(tid, subj, x, m, y) #combine them into a df

lmeout <- summary(lme(y ~ x + m, random=~1|subj, data=smaller.df))$tTable %>% data.frame() %>% tibble::rownames_to_column() #build a model, extract coefs, put coefs into a table } #another function to simulate all the parameters specified above sim_all_params <- function(params) { sim_results <- purrr::pmap(params, sim_sample) params %>% mutate(sim_results = sim_results) %>% tidyr::unnest(sim_results) } #execute the simulation and re-run it the pre-specified amount of times #sim_all_params(params) full_results <- purrr::rerun(.n = n_samples, sim_all_params(params)) %>% bind_rows(.id = 'sample_num') %>% as.data.frame(full_results) #save the full results in a table table <- full_results %>% dplyr::group_by(rowname, N, B1, items) %>% dplyr::summarise(power = sum(p.value < .05)/n()) #compute power for each row #print the table table  It is set up to produce a table of power at various combinations of effect size and sample size. Similar to the "power curve" functions you might see elsewhere. • You calculate the mean power for all three coefficients of your model, intercept, x and m. This gives roughly power$(0.05 + 1 + 1)/3 \approx 0.7\$
– Ute
Commented Jul 1, 2023 at 21:13
• +1 for a nice code template (when fixed, and purr::rerun got replaced) :-)
– Ute
Commented Jul 2, 2023 at 9:54
• Are you comfortable with adapting the code to newer version of purr?
– Ute
Commented Jul 2, 2023 at 14:06
• @Ute I am not sure exactly what has changed- apparently I haven't updated my R in a while :) Commented Jul 2, 2023 at 16:47

You need a small modification to your code - as it is now, it collects $$p$$-values for all three coefficients of your model, corresponding to Intercept, x and m. With the chosen values of B1 and B2, you have power 1 for x and m. For the Intercept, the power is $$0.05$$, so the results you get are reasonable.
1. add the names of your model coefficients to your output table. You can do this with function rownames_to_column from package tibble; in the last code line of function sim_sample, replace data.frame() by data.frame() %>% tibble::rownames_to_column()
2. In your main program, add rowname as the first grouping variable when you calculate the final table, that is, replace the fourth last line dplyr::group_by(N, B1, items) by  dplyr::group_by(rowname, N, B1, items)