I often need to conduct power and sample size calculations for multilevel regression, typically for a fixed effect. I use the simr package which works great, but I can't figure out how to test the standardized, instead of unstandardized regression coefficient in a simulated multilevel regression, or otherwise get it to test a direct measure of effect size.
For context, let’s say I want to conduct power calculations for the fixed predictor in the following simple random slope model
(please note that my version of lme4 requires the predictor entered both as a separate fixed effect and as the random slope, and my simr version requires me to specify the fixed effect. I've noticed most people use versions with different notation so just FYI)
model<-lmer(y ~ var1+(var1|id))
Based on previous research, I believe my fixed effect will be the size of .20, random intercept .60, residual .80 and random slope variance .10.
I assume I'll be able to collect a sample of 60 clusters and 10 observations per cluster.
I create the data:
id1<-1:60 index1<-1:10 id<-rep(id1, each=10) index<-rep(index1, times=60) var1<-rnorm(600, 3, 1) simdf<-data.frame(id, index, var1) #I cluster-mean-center the predictor library(misty) simdf$var1_c<-center(simdf$var1, type=”CWC”, cluster=simdf$id) #I create the parameters: fixed<-c(3.0, 0.2) random<-matrix(c(0.6, 0.05, 0.05, 0.1),2) resid<-0.8 #I construct the simulated model library(lme4) library(simr) modelsim<-makeLmer(y ~ var1_c+(var1_c|id), data=simdf, fixef=fixed, VarCorr=random, sigma=resid) #The output: Linear mixed model fit by REML ['lmerMod'] Formula: y ~ var1_c + (var1_c | id) Data: simdf REML criterion at convergence: 1615.7 Scaled residuals: Min 1Q Median 3Q Max -2.8692 -0.6599 0.0073 0.5803 3.7198 Random effects: Groups Name Variance Std.Dev. Corr id (Intercept) 0.60 0.7746 var1_c 0.10 0.3162 0.20 Residual 0.64 0.8000 Number of obs: 600, groups: id, 60 Fixed effects: Estimate Std. Error t value (Intercept) 3.0000 0.1052 28.518 var1_c 0.2000 0.0552 3.623 Correlation of Fixed Effects: (Intr) var1_c 0.144 #I conduct the power analysis: powerSim(modelsim, test=fixed("var1_c"), nsim=100, seed=1234) #Output Power for predictor 'var1_c', (95% confidence interval):= 94.00% (87.40, 97.77) Test: Kenward Roger (package pbkrtest) Effect size for var1_c is 0.20 Based on 100 simulations, (0 warnings, 0 errors) alpha = 0.05, nrow = 600
However, the fixed effect tested is based on the unstandardized regression coefficient which is not a direct measure of effect size. As I want to investigate within-person relations between the predictor and dependent, I can't standardize the predictor, it needs to be cluster-mean centered (or could I standardize it within cluster? I don't think I've ever seen anyone do that but maybe I'm just ignorant. Would that help?).
Do you have any suggestions? I can of course standardize the coefficients after running the model, and then report what the standardized beta ended up being, but this does not feel like the right approach and then I can't determine effect size in advance.