How to determine sample size needed for sufficient power to compare between levels of a fixed effect in a mixed model?

Question

I ran a pilot study for a within-subjects experiment where the fixed effect of interest is categorical, with 4 levels. The research question involves comparing which levels of this IV can statistically be considered similar versus different.

Here's a simplified example. Let's say that for Levels A, B, C, D, the mean reaction times in those conditions are 1000ms, 1100ms, 1150ms, and 1175ms. By rotating the reference group we determined that: (1) levels A B and C are all statistically different from each other, (2) Level D is different from A and B, but (3) C and D are equivalent.

As mentioned, I have pilot data, and am also happy to conduct a simulation-based analysis. But all the tools I've found seem to be telling me the sample size needed to detect an effect of the fixed effect in general, which is clearly not sufficient for my actual goals. For example, simR is convinced I only need roughly 12 participants for 80% power, despite the varied effect sizes and the small trial size (30) per condition.

I have also tried running simR based on just subsets of the data (for example, just comparing Levels 2 and 3 or 3 and 4), but of course this is incredibly flawed (e.g., Levels 3 and 4 suddenly become significantly different) and besides, I understand that I'd need a larger sample to account for the additional levels in the real analysis. So there's no real utility to this approach.

Please note that it is a strong convention in my field that I use a mixed effect model, and that is also what is documented in the project's pre-registered analysis plan. If possible, please do not just tell me that it would be easier to switch to [insert your favorite analysis method here].

Answer 1

I'm not sure if this solves the issue, but with simr package you can test the power for different sample sizes either for the "whole" categorical predictor or for specific contrasts. For instance like this

library(simr)
library(lme4)

# Make a toy dataset with 60 participants in 4 different conditions 
# and 30 trials per participant

subj <- factor(1:60)
cat <- letters[1:4]
subj_id <- rep(subj, each=30)
cat_pred <- rep(rep(cat, each=30), times=15)
simdf <- data.frame(subj_id, cat_pred)

# Specify fixed and random effects, you can use the ones from your 
# pilot data, I just plug in some values

b <- c(3.0, 1.2, 0.01, 0.9) # Intercept is set to 3.0, a-b contrast to 
  # 1.2, a-c contrast to 0.01 etc.
random <- 0.6 
  # Assumes that participant-level variance in outcome is 0.6
res <- 1.0 # residual variance in outcome is set to 1

# compose the mock model

modelsim <- makeLmer(y ~ (1|subj_id) + cat_pred, fixef=b, 
                     VarCorr=random, sigma=res, data=simdf)
summary(modelsim)

# To compute the power for the categorical predictor as a whole 
# (i.e., if there are any differences between levels), use

powercat <- powerSim(modelsim, test=fixed("cat_pred", "lr")) 
 # uses likelihood ratio test
powercat

# But if you want to check the power for, for instance a-b contrast, 
# use

powera_b <- powerSim(modelsim, test=fixed("cat_predb", "z"))
powera_b

You can then change the fixed effects estimates to get different comparisons, and number of participants to test different sample sizes. You can also test several different sample sizes at once using the powerCurve function, but I'm bad at using it so I'll leave that to the tutorials :), there are several good ones, e.g. https://humburg.github.io/Power-Analysis/simr_power_analysis.html .

EDITED TO ADD: you may want to initially add the command nsim=50 into the powerSim function to check that it works with smaller number of simulations, default is 1000 which takes a long time.