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:

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
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)

#I construct the simulated model

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:
var1_c 0.144

#I conduct  the power analysis:

powerSim(modelsim, test=fixed("var1_c"), nsim=100, seed=1234)


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.



Your Answer

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