2
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ 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$ $\endgroup$
    – Ute
    Commented Jul 1, 2023 at 21:13
  • 1
    $\begingroup$ +1 for a nice code template (when fixed, and purr::rerun got replaced) :-) $\endgroup$
    – Ute
    Commented Jul 2, 2023 at 9:54
  • $\begingroup$ Are you comfortable with adapting the code to newer version of purr? $\endgroup$
    – Ute
    Commented Jul 2, 2023 at 14:06
  • $\begingroup$ @Ute I am not sure exactly what has changed- apparently I haven't updated my R in a while :) $\endgroup$ Commented Jul 2, 2023 at 16:47

1 Answer 1

1
$\begingroup$

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.

A quick fix of your code:

  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)

$\endgroup$
2
  • $\begingroup$ Thanks-sometimes you just need another set of eyes! I will update my code in the original question so hopefully it can be useful to others. $\endgroup$ Commented Jul 2, 2023 at 9:26
  • 2
    $\begingroup$ yes, that is right - it is hard to spot glitches in one's own code. Your code seems otherwise very suitable for modification. I got some "deprecated" message for one of the tidyverse functions, maybe you will fix that, too ? $\endgroup$
    – Ute
    Commented Jul 2, 2023 at 9:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.