I want to conduct a power simulation for my research with crossed random effects model in the R environment. Specifically, in the pilot study, each participant randomly rates 10 stimuli out of a stimulus pool consisting of 100 stimuli. Each stimulus has 20 random raters. I fit a linear mixed effects model with random intercepts for both participants and stimuli (I didn't include any random slopes for some specific reasons).

I want to do a power simulation to estimate how many participants and stimuli should be used for the formal study. Moreover, I also want to take into account the ICC values and want to manipulate the number of raters for each stimulus. I didn't find any tutorials for my case. I have proposed the following three steps as a solution but I'm not sure whether it is the right way to do so:

  1. Simulate predictor values based on between- and within-stimulus correlations (I guess there is no way to simulate based on both between- and within-stimulus correlations and between- and within-participant correlations so I chose only between- and within-stimulus correlations).
  2. Simulate random intercepts for participants and stimuli.
  3. Calculate the outcomes from the regression model.

Besides the question of whether the steps are right, I also have a question about Step 1. The function sim.multilevel from the psych package might be useful to implement Step 1, but after checking the package manual I didn't understand how to fit the function for my case.

  • $\begingroup$ I have figured out how to use function sim.multilevel from the psych package to simulate finish step 1 although I found that between-stimulus correlations are far from being accurate even using a very large sample. But maybe I can use it anyway. $\endgroup$
    – Bo WANG
    Jun 24, 2023 at 9:45
  • $\begingroup$ I now have a new question because I forgot to mention that I also want to see the interaction between an individual differences variable (i.e., personality) and one of the predictors. The problem is that personality is a between-participant variable unlike other predictors. Therefore, I cannot use the function sim.multilevel to simulate all the predictors including the personality. However, personality also correlates with the predictors which means I cannot use the function rnorm to simulate itself as if it is unrelated with other predictors. $\endgroup$
    – Bo WANG
    Jun 24, 2023 at 9:49
  • $\begingroup$ My solution is that I can conduct another regression model where the personality is regressed on all the participant-level means of all other predictors (although this model doesn't make sense). Then after I use the function sim.multilevel simulate the data of the predictors except for the personality. I can calculate the personality data from the participant-level means of the simulated values of all other predictors. $\endgroup$
    – Bo WANG
    Jun 24, 2023 at 9:52

1 Answer 1



Generally speaking, there is a good body of literature that says performing power analysis on pilot data is not a great idea (see Lakens, 2022 cited below for more information).

enter image description here

While you haven't gone off and collected that full dataset already, the main issue is that you don't want to generate power analyses from pilot data or draw inferences about how to simulate power based off what you have found already. Rather, it should be derived from population-level data from previous studies.

The second thing to consider here is how to conduct power in the case you are in a position to. There are generally two methods you can employ. The first is a simulation-based power analysis, which should cover everything you are looking for in your question. There is a four-step guide at this website that is pretty comprehensive in terms of achieving this yourself (the fourth part specifically covers mixed models, but I advise going through the whole thing). Because this approach requires specificity, you will need to figure out what you will include in this simulation. If that is a lot to figure out, I know that the packages simr and mixedpower are useful for figuring out this information.


  • Green, P., & MacLeod, C. J. (2016). SIMR: an R package for power analysis of generalized linear mixed models by simulation. Methods in Ecology and Evolution, 7(4), 493–498. https://doi.org/10.1111/2041-210X.12504
  • Kain, M. P., Bolker, B. M., & McCoy, M. W. (2015). A practical guide and power analysis for GLMMs: Detecting among treatment variation in random effects. PeerJ, 3, e1226. https://doi.org/10.7717/peerj.1226
  • Kumle, L., Võ, M. L.-H., & Draschkow, D. (2021). Estimating power in (generalized) linear mixed models: An open introduction and tutorial in R. Behavior Research Methods, 53, 2528–2543. https://doi.org/10.3758/s13428-021-01546-0
  • Lakens, D. (2022). Sample size justification. Collabra: Psychology, 8(1), 33267.
  • $\begingroup$ Thank you so much! I will follow Lakens (2020) to use a more conservative effect size (e.g., using the lower limit of the 60% CI of the target effect size) for simulation. I understand that I can use the function extend from the package simr to extend the sample for simulation. The problem is that (if it is correct) the extend function can only keep the data structure (e.g., how many raters for each stimulus) and cannot change it. That's why I chose to define custom steps as a function to do the simulation. $\endgroup$
    – Bo WANG
    Jun 24, 2023 at 9:38
  • $\begingroup$ I would go for the simulation approach in the tutorial I linked then. $\endgroup$ Jun 24, 2023 at 9:49

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.