A-Priori Power Analysis for Multilevel Models

A Tutorial Based on Arend & Schäfer (2019)

Overview

In the previous tutorials, we demonstrated how to conduct post-hoc power analyses for the available multilevel data. We can also do a power analysis prior to having data for analysis. In this tutorial, we construct a multilevel model and conduct an a-priori power analysis using the simr package in R. The process described here can be used to obtain power estimates for a variety of common use cases, such as determining power for multilevel and growth models for grant proposals. It is adapted from the paper:

  • Arend, M. G., & Schäfer, T. (2019). Statistical power in two-level models: A tutorial based on Monte Carlo simulation. Psychological Methods, 24(1), 1-19. https://doi.org/10.1037/met0000195.

Preliminaries

Loading libraries used in this script.

library(psych)   #for describing data
library(lme4)    #for multilevel models
library(simr)    #for power analysis simulations
library(ggplot2) #for data visualization
library(dplyr)   #for data manipulation
library(plyr)    #for data manipulation including lists and arrays

Model Formulation

To obtain power estimates we must formulate a model and make determinations about effect size and other parameters. Our use case is based on a standard simple multilevel model of the form

\[y_{ti} = \beta_{0i} + \beta_{1i}x_{ti} + e_{ti} \tag{Level-1}\] where

\[\beta_{0i} = \gamma_{00} + \gamma_{01}Z_i + u_{0i} \\ \beta_{1i} = \gamma_{10} + \gamma_{11}Z_i + u_{1i} \tag{Level-2}\] with \[e_{it} \sim N(0,\sigma^{2})\] and \[[u_{0i},u_{1i}] \sim N(0,\Psi)\] where \[\Psi = \left[\begin{array} \t \tau_{00} & \tau_{01} \\ \tau_{01} & \tau_{11} \\ \end{array}\right]\]

After algebraic substitution, the main equation becomes: \[y_{ti} = \gamma_{00} + \gamma_{01}Z_i + \gamma_{10}x_{ti} + \gamma_{11}x_{ti}Z_i + [u_{0i} + u_{1i}x_{ti}] + e_{ti}\]

This equation more closely matches the R code layout that will be used for data generation and model fitting in the makeLmer() function used later.

We will need to specify assumptions about each and every parameter, including Level 1 (L1) and Level 2 (L2) sample sizes,

  • \(\gamma_{00}\) = fixed intercept
  • \(\gamma_{01}\) = L2 direct effect
  • \(\gamma_{10}\) = L1 direct effect
  • \(\gamma_{11}\) = cross-level interaction (CLI) effect,

as well as the Intraclass Correlation Coefficient, or ICC (relative balance of variance across L1 and L2). We do the parameter specification such that all the parameters and variances are in a “standardized” form. Following Arend and Schäfer (2019), we fix the L1 variance = 1.00 and scale everything relative to that standardization.

Specifying Standardized Input Parameters

We specify input parameters that underpin the power calculation, including significance level (alpha.S), cluster size (Size.clus), sample sizes (N.clus), standardized effect size estimates (L1_DE_standardized, L2_DE_standardized, CLI_E_standardized), and ICC (ICC). As proposed by Arend and Schäfer (2019), the fixed intercept \(\gamma_{00}\) will be set to 0.

#Seed value for replication
myseed <- 1234
#Significance level alpha at .05 
alpha.S <- .05

#L1 sample size n = 20
Size.clus <- 20
#L2 sample size N = 30
N.clus <- 40

#L1 direct effect, gamma_10.std = .30
L1_DE_standardized <- .30
#L2 direct effect, gamma_01.std = .30
L2_DE_standardized <- .30
#CLI effect, gamma_11.std  = .50
CLI_E_standardized <- .50

#ICC, rho.std = .30
ICC <- .30
#Std. random slope, tau_11.std = .09
rand.sl <- .09

For the selection of the random slope (rand.sl), we took guidance from Arend & Schaefer (2019). They state (pp. 6-7),

“… When no sufficiently informed estimator for the size of the random slope is available, choosing a medium effect size (tau_11.std = .09) is a good option, because small random slope variances would yield anti-conservative estimates and large slope variances would yield over-conservative estimates of the power for CLI effects.”

Creating Variables for Power Simulation

Next, we create a dataset with one L1-predictor, \(x_{ti}\) (standardized based on within-cluster variance), and one L2-predictor, \(Z_{i}\) (standardized based on total variance).

#creating t subscript variable
occasion <- rep(1:Size.clus)
#creating i subscript variable
person <- as.factor(1:N.clus)

#creating occasions x persons data frame
dat <- cbind(expand.grid("occasion"=occasion, "id"=person))

#creating standardized time-varying x variable (L1 predictor) from the occasion subscript 
#(e.g., standardizing time)
dat$x <- scale(dat$occasion)

#creating standardized time-invariant Z variable (L2 predictor) from the i subscript 
#(e.g., standardizing the id index in a numeric variable)
dat$Z <- scale(as.numeric(dat$id))

#printing the data
head(dat)
##   occasion id          x         Z
## 1        1  1 -1.6464789 -1.688221
## 2        2  1 -1.4731654 -1.688221
## 3        3  1 -1.2998518 -1.688221
## 4        4  1 -1.1265382 -1.688221
## 5        5  1 -0.9532246 -1.688221
## 6        6  1 -0.7799111 -1.688221

Adapting the Standardized Parameters

To standardize the variances, we calculate all the model parameters so that they are scaled in relation to the L1 variance. The resulting parameters will become inputs for the model-based data generation.

#L1 variance component (sigma^2) 
varL1 <- 1
#L2 variance component (tau_00) 
varL2 <- ICC/(1-ICC)
#Random slope variance (tau_11) 
varRS <- rand.sl*varL1

#L1 direct effect (gamma_10) 
L1_DE <- L1_DE_standardized*sqrt(varL1)
#L2 direct effect (gamma_01) 
L2_DE <- L2_DE_standardized*sqrt(varL2)
#CLI effect (gamma_11) 
CLI_E <- CLI_E_standardized*sqrt(varRS)

Creating Conditional Variances

Using these unconditional variance components, we derive the conditional variance components for the population model.

#Creating conditional variances
#L1 variance (sigma_{Y|X}^2) 
s <- sqrt((varL1)*(1-(L1_DE_standardized^2)))
#L2 variance (tau_00Y|Z) 
V1 <-  varL2*(1-(L2_DE_standardized^2))
#Random slope variance (tau_11Y|XZ) 
rand_sl.con <- varRS*(1-(CLI_E_standardized^2))

Creating a Population Model for Simulation

We collect the parameters into a vector, b, and a matrix, V2.

#Vector of fixed effects (fixed intercept, L1-, L2- direct and CLI effect)
b <- c(0, L1_DE, L2_DE, CLI_E)
#Random effects covariance matrix (with covariances set to 0)
V2 <- matrix(c(V1,0,0,rand_sl.con), 2)

We then use the calculated parameters to generate data in accordance with the hypothesized model using the makeLmer() function.

The model can have a fixed (1|id) or random (x|id) slope (when the slope is fixed, VarCorr = V1 must be used). We include the random slope in this tutorial.

#Making model  
#x = gamma_10; Z = gamma_01; x:Z = gamma_11
model <- simr::makeLmer(y ~ x + Z + x:Z + (x|id),
                        fixef = b,
                        VarCorr = V2,
                        sigma = s,
                        data = dat)
print(model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + Z + x:Z + (x | id)
##    Data: dat
## REML criterion at convergence: 2354.852
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  id       (Intercept) 0.6245       
##           x           0.2598   0.00
##  Residual             0.9539       
## Number of obs: 800, groups:  id, 40
## Fixed Effects:
## (Intercept)            x            Z          x:Z  
##      0.0000       0.3000       0.1964       0.1500

We can plot the predicted individual trajectories.

#Plotting data
dat$y <- predict(model)
dat %>%
  ggplot(aes(y=y, x=occasion, group=id, color=Z)) +
  geom_line(alpha=.5) +
  scale_color_gradient(low="red", high="blue") +
  scale_x_continuous(limits=c(1,20),
                     breaks=c(1,5,10,15,20), 
                     name="Occasion") +
  scale_y_continuous(name="Model Predictions") +
  theme_bw()

Cool! As described in the model results, we see that predicted values tend to increase over time (x, or “Occasion” in the plot = 0.3000), with higher Z values associated with higher predicted values (Z, represented by color in the plot = 0.1964) specifically on later occasions (x:Z = 0.1500).

Now let’s implement the power analysis.

Implementing the Power Analysis

We will estimate power for x based on 10 repetitions (for more repetitions, change nsim to > 1000).

Note that the code method="kr" here specifies the Kenward-Roger approximation for correcting standard errors. See ?tests for available options/details about the method argument.

#Power for L1_DE_standardized
sim.ef_gamma_10 <- simr::powerSim(model, 
                                  test = simr::fixed("x", method="kr"),
                                  seed = myseed, #set for replication
                                  alpha = alpha.S, 
                                  nsim=10)
## Simulating: |                                                                  |Simulating: |======                                                            |Simulating: |=============                                                     |Simulating: |===================                                               |Simulating: |==========================                                        |Simulating: |=================================                                 |Simulating: |=======================================                           |Simulating: |==============================================                    |Simulating: |====================================================              |Simulating: |===========================================================       |Simulating: |==================================================================|
print(sim.ef_gamma_10)
## Power for predictor 'x', (95% confidence interval):
##       100.0% (69.15, 100.0)
## 
## Test: Kenward Roger (package pbkrtest)
##       Effect size for x is 0.30
## 
## Based on 10 simulations, (10 warnings, 0 errors)
## alpha = 0.05, nrow = 800
## 
## Time elapsed: 0 h 0 m 0 s

Next, we estimate power for Z based on 10 repetitions.

#Power for L2_DE_standardized
sim.ef_gamma_01 <- simr::powerSim(model, 
                                  test = simr::fixed("Z", method="kr"),
                                  seed = myseed, #set for replication
                                  alpha = alpha.S, 
                                  nsim = 10)
## Simulating: |                                                                  |Simulating: |======                                                            |Simulating: |=============                                                     |Simulating: |===================                                               |Simulating: |==========================                                        |Simulating: |=================================                                 |Simulating: |=======================================                           |Simulating: |==============================================                    |Simulating: |====================================================              |Simulating: |===========================================================       |Simulating: |==================================================================|
print(sim.ef_gamma_01)
## Power for predictor 'Z', (95% confidence interval):
##       40.00% (12.16, 73.76)
## 
## Test: Kenward Roger (package pbkrtest)
##       Effect size for Z is 0.20
## 
## Based on 10 simulations, (10 warnings, 0 errors)
## alpha = 0.05, nrow = 800
## 
## Time elapsed: 0 h 0 m 0 s

Finally, we estimate power for x:Z based on 10 repetitions.

#Power for CLI_E_standardized
sim.ef_gamma_11 <-  simr::powerSim(model, 
                                  test = simr::fixed("x:Z", method="kr"),
                                  seed = myseed, #set for replication
                                  alpha = alpha.S, 
                                  nsim = 10)
## Simulating: |                                                                  |Simulating: |======                                                            |Simulating: |=============                                                     |Simulating: |===================                                               |Simulating: |==========================                                        |Simulating: |=================================                                 |Simulating: |=======================================                           |Simulating: |==============================================                    |Simulating: |====================================================              |Simulating: |===========================================================       |Simulating: |==================================================================|
print(sim.ef_gamma_11)
## Power for predictor 'x:Z', (95% confidence interval):
##       90.00% (55.50, 99.75)
## 
## Test: Kenward Roger (package pbkrtest)
##       Effect size for x:Z is 0.15
## 
## Based on 10 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 800
## 
## Time elapsed: 0 h 0 m 0 s

Plotting Power Curves

We can also plot a power curve with different effect sizes on the x-axis and power on the y-axis. The following code chunk defines a function called powerCurveData() that tests a set of effect sizes. You can specify the series of effect sizes that should be used. Run this entire code chunk at the same time.

powerCurveData = function(
  #alpha  
  alpha.S = c(.05),
  #L1 sample size n = 20 
  Size.clus = c(20),
  #L2 sample size N = 40 
  N.clus = c(40),
  #L1 direct effect, gamma_10.std = .30 
  L1_DE_standardized = c(.30),
  #L2 direct effect, gamma_01.std = .30 
  L2_DE_standardized = c(.30),
  #CLI effect, gamma_11.std = .50 
  CLI_E_standardized = c(.50),
  #ICC, rho.std = .30 
  ICC = c(.30),
  #Std. random slope, tau_11.std = .09 
  rand.sl = c(.09),
  #number of simulations
  numbersim = 10,
  #setting seed
  myseed = 1234 ) {
  
  # Create a dataframe where simPower results will be stored (and then used to make a plot)
  powerData <- cbind(expand.grid("alpha.S" = alpha.S,
                                 "Size.clus" = Size.clus,
                                 "N.clus" = N.clus,
                                 "L1_DE_standardized" = L1_DE_standardized, 
                                 "L2_DE_standardized" = L2_DE_standardized,
                                 "CLI_E_standardized" = CLI_E_standardized,
                                 "ICC" = ICC,
                                 "rand.sl" = rand.sl,
                                 "numbersim" = numbersim,
                                 "myseed" = myseed))
  
  #Generating Power Estimates
  for(i in 1:nrow(powerData)) {
    #  i <- 1
    # Create data (any of these could be parameters in the function too)
    #creating t subscript indicator
    occasion <- rep(1:powerData$Size.clus[i])
    #creating i subscript indicator
    person <- as.factor(1:powerData$N.clus[i])
    
    #creating occasions x persons data frame
    dat <- cbind(expand.grid("occasion"=occasion, "id"=person))
    
    #creating standardized time-varying x variable (L1 predictor)
    #using shortcut to make x from the occasion subscript (e.g., standardizing time)
    dat$x <- scale(dat$occasion)
    
    #creating standardized time-invariant Z variable (L2 predictor)
    #using shortcut to make Z from the i subscript (e.g., standardizing the id index in a numeric variable)
    dat$Z <- scale(as.numeric(dat$id))
    
    #L1 variance component (sigma^2) 
    varL1 <- 1
    #L2 variance component (tau_00) 
    varL2 <- powerData$ICC[i]/(1-powerData$ICC[i])
    #Random slope variance (tau_11) 
    varRS <- powerData$rand.sl[i]*varL1
    
    #L1 direct effect (gamma_10) 
    L1_DE <- powerData$L1_DE_standardized[i]*sqrt(varL1)
    #L2 direct effect (gamma_01) 
    L2_DE <- powerData$L2_DE_standardized[i]*sqrt(varL2)
    #CLI effect (gamma_11) 
    CLI_E <- powerData$CLI_E_standardized[i]*sqrt(varRS)
    
    #Creating conditional variances
    #L1 variance (sigma_{Y|X}^2) 
    s <- sqrt((varL1)*(1-(powerData$L1_DE_standardized[i]^2)))
    #L2 variance (tau_00Y|Z) 
    V1 <-  varL2*(1-(powerData$L2_DE_standardized[i]^2))
    #Random slope variance (tau_11Y|XZ) 
    rand_sl.con <- varRS*(1-(powerData$CLI_E_standardized[i]^2))
    
    #Creating fixed and random effects
    #Vector of fixed effects (fixed intercept, L1-, L2- direct and CLI effect).
    b <- c(0, L1_DE, L2_DE, CLI_E)
    #Random effects covariance matrix (with covariances set to 0).
    V2 <- matrix(c(V1,0,0,rand_sl.con), 2)
    
    #Making model  
    #x = gamma_10; Z = gamma_01; x:Z = gamma_11
    model <- simr::makeLmer(y ~ x + Z + x:Z + (x|id),
                            fixef = b, 
                            VarCorr = V2, 
                            sigma = s, 
                            data = dat)
    
    #Power for L1_DE_standardized
    sim.ef_gamma_10 <- simr::powerSim(model, 
                                      test = simr::fixed("x", method="kr"),
                                      seed = powerData$myseed[i], #set for replication
                                      alpha = powerData$alpha.S[i], 
                                      nsim = numbersim, 
                                      progress = FALSE) #options
    #Add results to dataframe
    powerData[i,"Power_L1_DE"] <- summary(sim.ef_gamma_10)$mean #Power estimate
    powerData[i,"Power_L1_DE_lower"] <- summary(sim.ef_gamma_10)$lower #Power estimate lower CI
    powerData[i,"Power_L1_DE_upper"] <- summary(sim.ef_gamma_10)$upper #Power estimate upper CI
    
    #Power for L2_DE_standardized
    sim.ef_gamma_01 <- simr::powerSim(model, 
                                      test = simr::fixed("Z", method="kr"),
                                      seed = powerData$myseed[i], #set for replication
                                      alpha = powerData$alpha.S[i], 
                                      nsim = numbersim,
                                      progress = FALSE) #options
    #Add results to dataframe
    powerData[i,"Power_L2_DE"] <- summary(sim.ef_gamma_01)$mean #Power estimate
    powerData[i,"Power_L2_DE_lower"] <- summary(sim.ef_gamma_01)$lower #Power estimate lower CI
    powerData[i,"Power_L2_DE_upper"] <- summary(sim.ef_gamma_01)$upper #Power estimate upper CI
    
    #Power for CLI_E_standardized
    sim.ef_gamma_11 <- simr::powerSim(model, 
                                      test = simr::fixed("x:Z", method="kr"),
                                      seed = powerData$myseed[i], #set for replication
                                      alpha = powerData$alpha.S[i], 
                                      nsim = numbersim,
                                      progress = FALSE) #options
    #Add results to dataframe
    powerData[i,"Power_CLI_E"] <- summary(sim.ef_gamma_11)$mean #Power estimate
    powerData[i,"Power_CLI_E_lower"] <- summary(sim.ef_gamma_11)$lower #Power estimate lower CI
    powerData[i,"Power_CLI_E_upper"] <- summary(sim.ef_gamma_11)$upper #Power estimate upper CI
    
  }
  return(powerData)
}

We next run a power analysis for each of the three effects. Here we assume that 191 participants (N.clus = c(191)) are repeatedly measured for three times each (Size.clus = c(3)).

powerCurve_L1 <- powerCurveData(Size.clus = c(3), N.clus = c(191), 
                                L1_DE_standardized = c(.05,.1,.2,.3,.4,.5),
                                L2_DE_standardized = c(.0),
                                CLI_E_standardized = c(.0))

powerCurve_L2 <- powerCurveData(Size.clus = c(3), N.clus = c(191), 
                                L1_DE_standardized = c(.0),  
                                L2_DE_standardized = c(.05,.1,.2,.3,.4,.5),
                                CLI_E_standardized = c(.0))

powerCurve_CLI <- powerCurveData(Size.clus = c(3), N.clus = c(191), 
                                 L1_DE_standardized = c(.0),  
                                 L2_DE_standardized = c(.0),
                                 CLI_E_standardized = c(.1,.2,.3,.4,.5,.6,.7,.8))

Using the power analyses conducted above, we’ll make power curve plots for each of the three effects. We’ll start with the power curve plot for L1_DE_standardized.

powerCurve_L1 %>%
  ggplot(aes(y=Power_L1_DE, x=L1_DE_standardized)) +
  #power goal
  geom_hline(yintercept=0.8, color="gray40", lty = 2) +
  #power curve
  geom_line(color="black") +
  geom_point() +
  geom_errorbar(aes(ymin=Power_L1_DE_lower, ymax=Power_L1_DE_upper)) +
  scale_x_continuous(breaks=powerCurve_L1$L1_DE_standardized,
                     name="Effect Size") +
  scale_y_continuous(limits=c(0,1.01), 
                     breaks=c(0.0,0.2,0.4,0.6,0.8,1.0), 
                     name="Power") +  
  theme_classic() +
  theme(axis.title=element_text(size=14),
       axis.text=element_text(size=12,colour="black"),
       plot.title=element_text(size=14, hjust=.5)) +
  ggtitle("Power Curve for L1_DE_standardized")

We see that the standardized effect size of L1 direct effect (\(\gamma_{10}\)) [x] should be larger than 0.20 for adequate power (using the typical threshold of .80).

We’ll next make the power curve plot for L2_DE_standardized.

powerCurve_L2 %>%
  ggplot(aes(y=Power_L2_DE, x=L2_DE_standardized)) +
  #power goal
  geom_hline(yintercept=0.8, color="gray40", lty = 2) +
  #power curve
  geom_line(color="black") +
  geom_point() +
  geom_errorbar(aes(ymin=Power_L2_DE_lower, ymax=Power_L2_DE_upper)) +
  scale_x_continuous(breaks=powerCurve_L2$L2_DE_standardized,
                     name="Effect Size") +
  scale_y_continuous(limits=c(0,1.01), 
                     breaks=c(0.0,0.2,0.4,0.6,0.8,1.0), 
                     name="Power") +  
  theme_classic() +
  theme(axis.title=element_text(size=14),
       axis.text=element_text(size=12,colour="black"),
       plot.title=element_text(size=14, hjust=.5)) +
  ggtitle("Power Curve for L2_DE_standardized")

We see that the standardized effect size of L2 direct effect (\(\gamma_{01}\)) [Z] should be larger than 0.30 for adequate power (again, using the threshold of .80).

Finally, we’ll make the power curve plot for CLI_E_standardized.

powerCurve_CLI %>%
  ggplot(aes(y=Power_CLI_E, x=CLI_E_standardized)) +
  #power goal
  geom_hline(yintercept=0.8, color="gray40", lty = 2) +
  #power curve
  geom_line(color="black") +
  geom_point() +
  geom_errorbar(aes(ymin=Power_CLI_E_lower, ymax=Power_CLI_E_upper)) +
  scale_x_continuous(breaks=powerCurve_CLI$CLI_E_standardized,
                     name="Effect Size") +
  scale_y_continuous(limits=c(0,1.01), 
                     breaks=c(0.0,0.2,0.4,0.6,0.8,1.0), 
                     name="Power") +  
  theme_classic() +
  theme(axis.title=element_text(size=14),
       axis.text=element_text(size=12,colour="black"),
       plot.title=element_text(size=14, hjust=.5)) +
  ggtitle("Power Curve for CLI_E_standardized")

We see that the standardized effect size of cross-level interaction (CLI) effect (\(\gamma_{11}\)) [x:Z] should be larger than 0.50 for adequate power (based on a .80 threshold).

Conclusion

In this tutorial, we illustrated how to specify a two-level model and conduct a power analysis in R using the simr package. Prior to having data for analysis, researchers can assess power based on a set of input parameters (e.g., L1/L2 direct effect, CLI effect, ICC, random slope) that specify the population model of interest. Customizations could include:

  • Data generation (e.g., number of occasions, number of persons, covariates)

  • Effect of interest for which power is calculated (currently labeled “Effect Size” in the plots above, which corresponds to L1 direct effect, L2 direct effect, and CLI effect)

  • Fixed and random effects parameters (e.g., ICC, random slope)

  • Power curves based on different alpha values (e.g., \(\alpha = .005\))

Thanks for powering through with us!