Using mlmpower Package to Conduct Multilevel Power Analysis

Illustration 1: Cross-Sectional Power Analysis

The first illustration demonstrates a power analysis for a cross-sectional application of the following multilevel model. To provide a substantive context, consider a prototypical education example where students are nested in schools. Intraclass correlations for achievement-related outcomes often range between .10 and .25 (Hedges & Hedberg, 2007; Hedges & Hedberg, 2013; Sellstrom & Bremberg, 2006; Spybrook et al., 2011; Stockford, 2009). To accommodate uncertainty about this important parameter, the power simulations investigate intraclass correlation values of .10 and .25.

The multilevel model for the illustration is

\[ \begin{split} Y_{ij} &= \left( \beta_{0} + b_{0j} \right) + \left( \beta_{1} + b_{1j} \right)\left( X_{1ij} - {\overline{X}}_{1j} \right) + \beta_{2}\left( X_{2ij} - {\overline{X}}_{2} \right) \\ &\phantom{=}+ \ \beta_{3}\left( Z_{1j} - {\overline{Z}}_{1} \right) + \ \beta_{4}\left( Z_{2j} - {\overline{Z}}_{2} \right) + \beta_{5}\left( X_{1ij} - {\overline{X}}_{1j} \right)\left( Z_{1j} - {\overline{Z}}_{1} \right) + \varepsilon_{ij}\\[1em] \mathbf{b}_{j}\ &\sim\ N\left( 0,\mathbf{\Sigma}_{\mathbf{b}} \right)\ \ \ \ \ \mathbf{\Sigma}_{\mathbf{b}}=\begin{pmatrix} \sigma_{b_{0}}^{2} & \\ \sigma_{b_{1},b_{0}} & \sigma_{b_{1}}^{2} \\ \end{pmatrix}\mathbf{\ \ \ \ \ }\varepsilon_{ij}\ \sim\ N(0,\sigma_{\varepsilon}^{2}) \end{split} \]

where \(Y_{ij}\) is the outcome score for observation i in cluster j, \(X_{1ij}\) is a focal within-cluster predictor, \({\overline{X}}_{1j}\) is the variable’s level-2 cluster mean, \(X_{2ij}\) is a grand mean centered level-1 covariate, \(Z_{1j}\) is a level-2 moderator score for cluster j, and \(Z_{2j}\) is a level-2 covariate. Turning to the residual terms, \(b_{0j}\) and \(b_{1j}\) are between-cluster random effects that capture residual variation in the cluster-specific intercepts and slopes, and \(\varepsilon_{ij}\) is a within-cluster residual. By assumption, the random effects follow a multivariate normal distribution with a between-cluster covariance matrix \(\mathbf{\Sigma}_{\mathbf{b}}\), and within-cluster residuals are normal with constant variance \(\sigma_{\varepsilon}^{2}\).

Importantly, group mean centering \(X_{1}\) yields a pure within-cluster predictor, whereas grand mean centering \(X_{2}\) gives a predictor with two sources of variation. Although it need not be true in practice, the within- and between-cluster parts of any level-1 predictor variables with non-zero intraclass correlations share a common slope coefficient. To simplify notation, Equation 1 can be rewritten as

\[Y_{ij} = \left( \beta_{0} + b_{0j} \right) + \left( \beta_{1} + b_{1j} \right)X_{1ij}^{w} + \beta_{2}\left( X_{2ij}^{w} + X_{2j}^{b} \right) + \beta_{3}Z_{1j}^{b} + \beta_{4}Z_{2j}^{b} + \beta_{5}X_{1ij}^{w}Z_{1j}^{b} + \varepsilon_{ij}\]

where the w and b superscripts reference within- and between-cluster deviation scores, respectively.

The key inputs to the model object are the unconditional intraclass correlation and effect size values. This illustration investigates power for intraclass correlation values of .10 and .25. The within- and between-cluster fixed effects are set to \(R_{w}^{2}\) = .065 and \(R_{b}^{2}\) = .065, respectively, the sum of which corresponds to Cohen’s (1988) medium effect size benchmark. Note that \(R_{b}^{2}\) cannot exceed the outcome’s intraclass correlation, as this would imply that between-cluster predictors explain more than 100% of the available variance. Following conventional wisdom that interaction effects tend to produce small effects (Aguinis, Beaty, Boik, & Pierce, 2005; Chaplin, 1991), the cross-level product is assigned \(R_{p}^{2}\) = .01. Finally, the random coefficient effect size is set to \(R_{rc}^{2}\) = .03 based on values from Table 2 in Enders, Keller, and Woller (2023).

By default, mlmpower assigns equal weights to all quantities contributing to a particular source of variance. To assign non-uniform weights, assume \(X_{1}\) and \(Z_{1}\) (the interacting variables) are the focal predictors and \(X_{2}\) and \(Z_{2}\) are covariates. To mimic a situation where the covariates explain a small amount of variation, the fixed effect \(R^{2}\) values are predominantly allocated to the focal predictors using weights of .80 and .20. A small covariate allocation could be appropriate for a set of background or sociodemographic characteristics that weakly predict the outcome. Researchers often need power estimates for individual partial regression slopes. Although the weights do not exactly carve \(R_{w}^{2}\) and \(R_{b}^{2}\) into additive components when predictors are correlated, adopting weak associations among the regressors allows us to infer that each focal predictor roughly accounts for 5% of the explained variation at each level (i.e., \(.80 \times R_{w}^{2}\) or \(R_{b}^{2} \approx .05\)).

The following code block shows the mlmpower model object for the example.

example1 <- (
    effect_size(
        icc          = c(.10, .25),
        within       = .065,
        between      = .065,
        product      = .01,
        random_slope = .03
    )
    + outcome('y', mean = 50, sd = 10)
    + within_predictor('x1', icc = 0, weight = .80)
    + within_predictor('x2', weight = .20)
    + between_predictor('z1', weight = .80)
    + between_predictor('z2', weight = .20)
    + product('x1','z1', weight = 1)
    + random_slope('x1', weight = 1)
)

The mlmpower package groups \(R^{2}\) values and intraclass correlations into a single object called effect_size(), with the five inputs separated by commas. The within argument corresponds to \(R_{w}^{2}\), the between argument aligns with \(R_{b}^{2}\), the product argument specifies \(R_{p}^{2}\), and the random_slope argument corresponds to \(R_{rc}^{2}\). The icc argument assigns a global intraclass correlation to all level-1 variables, and separate simulations are performed for each requested level.

Variable attributes are referenced by adding the following objects: outcome(), within_predictor(), between_predictor(), product(), and random_slope(). All five objects have an argument for the variable name, weight (weight =), mean (mean =), and standard deviation (sd =). Level-1 variables additionally have an intraclass correlation argument (icc =) that supersedes the global setting in effect_size().

The previous code block assigns explicit weights to all variables contributing to a given effect. The unit weights in the product() and random_slope() objects result from a single variable determining those sources of variation. Next, the within_predictor('x1', icc = 0, weight = .80) object overrides the global intraclass correlation setting, defining \(X_{1}\) as a pure within-cluster deviation variable with no between-cluster variation. Finally, with the exception of the outcome variable, which has a mean and standard deviation of 50 and 10, the script accepts the default settings for the means and variances of the predictors (0 and 1, respectively).

The multilevel model parameters also require three sets of correlations: within-cluster correlations among level-1 predictors, between-cluster correlations among level-2 predictors, and random effect correlations. These correlations are optional, but they can specified using the correlations() object. The earlier code block omits this object, thereby accepting the default specification. When the correlations() object is omitted, the mlmpower package iteratively samples all correlations from a uniform distribution between .10 and .30, such that the resulting power estimates average over a distribution of possible associations. This default range spans Cohen’s (1988) small to medium effect size benchmarks, and it brackets common correlation values from published research (Bosco, Aguinis, Singh, Field, & Pierce, 2015; Funder & Ozer, 2019; Gignac & Szodorai, 2016).

The default specifications could be invoked by explicitly appending the correlations() object to the earlier script, as follows.

example1 <- (
    effect_size(
        icc          = c(.10, .25),
        within       = .065,
        between      = .065,
        product      = .01,
        random_slope = .03
    )
    + outcome('y', mean = 50, sd = 10)
    + within_predictor('x1', icc = 0, weight = .80)
    + within_predictor('x2', weight = .20)
    + between_predictor('z1', weight = .80)
    + between_predictor('z2', weight = .20)
    + product('x1','z1', weight = 1)
    + random_slope('x1', weight = 1)
    + correlations(
        within  = random(0.1, 0.3),
        between = random(0.1, 0.3),
        randeff = random(0.1, 0.3)
    )
)

Researchers can modify the upper and lower limits of each correlation range, or they can specify a constant correlation by inputting a scalar value. For example, specifying randeff = 0 would define \(\mathbf{\Sigma}_{\mathbf{b}}\) as a diagonal matrix. The illustrative simulations presented in Enders et al. (2023) suggest that predictor and random effect correlations tend not to matter very much.

It may be instructive to inspect the population parameters prior to running the simulation. Executing summary(example1) returns tabular summaries of the multilevel model parameters.

summary(example1)
#> ! Parameters are solved based on average correlation.
#> $`icc = 0.1`
#> $r2
#>                                           Proportion
#> Variance Within-Cluster Fixed Effects    0.065000000
#> Variance Cross-Level Interaction Effects 0.010000000
#> Variance Random Slopes                   0.030000000
#> Within-Cluster Error Variance            0.795000000
#> Variance Between-Cluster Fixed Effects   0.071476494
#> Incremental Variance Level-2 Predictors  0.065000000
#> Between Variance Level-1 Covariates      0.006476494
#> Variance Random Intercepts               0.028523506
#> 
#> $phi_b
#>            x2_b         z1         z2
#> x2_b 0.10000000 0.06324555 0.06324555
#> z1   0.06324555 1.00000000 0.20000000
#> z2   0.06324555 0.20000000 1.00000000
#> 
#> $phi_p
#>         x1_w*z1
#> x1_w*z1       1
#> 
#> $phi_w
#>           x1_w      x2_w
#> x1_w 1.0000000 0.1897367
#> x2_w 0.1897367 0.9000000
#> 
#> $mean_Z
#> z1 z2 
#>  0  0 
#> 
#> $mean_X
#> x1 x2 
#>  0  0 
#> 
#> $var_e
#>           Value
#> Res. Var.  79.5
#> 
#> $tau
#>           Icept      x1_w
#> Icept 2.8523506 0.5850488
#> x1_w  0.5850488 3.0000000
#> 
#> $gammas
#>              Value
#> Icept   50.0000000
#> x1_w     2.3646137
#> x2_w     0.6231304
#> x1_w*z1  1.0000000
#> x2_b     0.6231304
#> z1       2.4308622
#> z2       0.6077155
#> 
#> $mean_Y
#>  y 
#> 50 
#> 
#> 
#> $`icc = 0.25`
#> $r2
#>                                           Proportion
#> Variance Within-Cluster Fixed Effects    0.065000000
#> Variance Cross-Level Interaction Effects 0.010000000
#> Variance Random Slopes                   0.030000000
#> Within-Cluster Error Variance            0.645000000
#> Variance Between-Cluster Fixed Effects   0.074006354
#> Incremental Variance Level-2 Predictors  0.065000000
#> Between Variance Level-1 Covariates      0.009006354
#> Variance Random Intercepts               0.175993646
#> 
#> $phi_b
#>      x2_b  z1  z2
#> x2_b 0.25 0.1 0.1
#> z1   0.10 1.0 0.2
#> z2   0.10 0.2 1.0
#> 
#> $phi_p
#>         x1_w*z1
#> x1_w*z1       1
#> 
#> $phi_w
#>           x1_w      x2_w
#> x1_w 1.0000000 0.1732051
#> x2_w 0.1732051 0.7500000
#> 
#> $mean_Z
#> z1 z2 
#>  0  0 
#> 
#> $mean_X
#> x1 x2 
#>  0  0 
#> 
#> $var_e
#>           Value
#> Res. Var.  64.5
#> 
#> $tau
#>           Icept     x1_w
#> Icept 17.599365 1.453246
#> x1_w   1.453246 3.000000
#> 
#> $gammas
#>              Value
#> Icept   50.0000000
#> x1_w     2.3646137
#> x2_w     0.6826052
#> x1_w*z1  1.0000000
#> x2_b     0.6826052
#> z1       2.4308622
#> z2       0.6077155
#> 
#> $mean_Y
#>  y 
#> 50

Having specified the target model, you next use the power_analysis() function to conduct simulations. The function requires four inputs: the model argument specifies the parameter value object (e.g., example1), the replications input specifies the number of artificial data sets, n_between is a vector of level-2 sample sizes, and n_within is a vector of level-1 sample sizes (i.e., the number of observations per cluster). The code block below pairs four level-2 sample size values (\(J\) = 30, 60, 90, and 120) with three level-1 sample sizes (\(n_{j}\) = 10, 20, or 30 observations per cluster), and it requests 2,000 artificial data sets for each combination.

# Set seed for replicable results
set.seed(2318971)

# Run Power Analysis
powersim1 <-
    power_analysis(
        model = example1,
        replications = 2000,
        n_between = c(30, 60, 90, 120),
        n_within = c(10, 20, 30)
    )

The package uses lme4 (Bates et al., 2021) for model fitting, and it defaults to an alpha level of .05 for all significance tests. Significance tests of random slope variation use a likelihood ratio test with a mixture chi-square reference distribution (i.e., a chi-bar distribution; Snijders & Bosker, 2012, p. 99), as implemented in the varTestnlme package (Baey & Kuhn, 2022). Executing summary(powersim1) returns the tabular summaries of the simulation results shown

summary(powersim1)
#> 
#> ── Power Analysis Specification ──
#> 
#> Replications = 2000
#> N_Within = 10, 20, and 30
#> N_Between = 30, 60, 90, and 120
#> 
#> ── lme4 Model Syntax
#> 
#> y ~ 1 + cwc(x1) + cgm(x2) + cwc(x1):cgm(z1) + cgm(z1) + cgm(z2) + (1 + cwc(x1)
#> | `_id`)
#> 
#> ℹ cwc() = centering within cluster
#> ℹ cgm() = centering with grand mean
#> 
#> ── Effect Sizes
#> 
#> • WITHIN = 0.065
#> • BETWEEN = 0.065
#> • RANDOM SLOPE = 0.03
#> • PRODUCT = 0.01
#> 
#> ── Power Analysis Results ──
#> 
#> ── global icc = 0.1, n_within = 10, n_between = 30
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         0.94 ± 0.01
#> Fixed:  cgm(x2)         0.23 ± 0.02
#> Fixed:  cgm(z1)         0.92 ± 0.01
#> Fixed:  cgm(z2)         0.15 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.29 ± 0.02
#> Random: Slopes Test     0.17 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 10, n_between = 30
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         0.95 ± 0.01
#> Fixed:  cgm(x2)         0.26 ± 0.02
#> Fixed:  cgm(z1)         0.66 ± 0.02
#> Fixed:  cgm(z2)         0.11 ± 0.01
#> Fixed:  cwc(x1):cgm(z1) 0.36 ± 0.02
#> Random: Slopes Test     0.23 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 20, n_between = 30
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         0.99 ± 0.00
#> Fixed:  cgm(x2)         0.35 ± 0.02
#> Fixed:  cgm(z1)         0.98 ± 0.01
#> Fixed:  cgm(z2)         0.22 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.48 ± 0.02
#> Random: Slopes Test     0.46 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 20, n_between = 30
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         0.99 ± 0.00
#> Fixed:  cgm(x2)         0.42 ± 0.02
#> Fixed:  cgm(z1)         0.75 ± 0.02
#> Fixed:  cgm(z2)         0.12 ± 0.01
#> Fixed:  cwc(x1):cgm(z1) 0.54 ± 0.02
#> Random: Slopes Test     0.59 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 30
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         0.99 ± 0.00
#> Fixed:  cgm(x2)         0.53 ± 0.02
#> Fixed:  cgm(z1)         0.99 ± 0.00
#> Fixed:  cgm(z2)         0.26 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.56 ± 0.02
#> Random: Slopes Test     0.73 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 30
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.60 ± 0.02
#> Fixed:  cgm(z1)         0.75 ± 0.02
#> Fixed:  cgm(z2)         0.12 ± 0.01
#> Fixed:  cwc(x1):cgm(z1) 0.59 ± 0.02
#> Random: Slopes Test     0.83 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 10, n_between = 60
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.40 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.25 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.56 ± 0.02
#> Random: Slopes Test     0.31 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 10, n_between = 60
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.43 ± 0.02
#> Fixed:  cgm(z1)         0.94 ± 0.01
#> Fixed:  cgm(z2)         0.15 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.63 ± 0.02
#> Random: Slopes Test     0.46 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 20, n_between = 60
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.63 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.40 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.78 ± 0.02
#> Random: Slopes Test     0.80 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 20, n_between = 60
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.70 ± 0.02
#> Fixed:  cgm(z1)         0.96 ± 0.01
#> Fixed:  cgm(z2)         0.17 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.83 ± 0.02
#> Random: Slopes Test     0.90 ± 0.01
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 60
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.81 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.45 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.86 ± 0.02
#> Random: Slopes Test     0.95 ± 0.01
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 60
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.87 ± 0.01
#> Fixed:  cgm(z1)         0.97 ± 0.01
#> Fixed:  cgm(z2)         0.19 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.90 ± 0.01
#> Random: Slopes Test     0.99 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 10, n_between = 90
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.54 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.35 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.74 ± 0.02
#> Random: Slopes Test     0.46 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 10, n_between = 90
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.60 ± 0.02
#> Fixed:  cgm(z1)         0.99 ± 0.00
#> Fixed:  cgm(z2)         0.21 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.79 ± 0.02
#> Random: Slopes Test     0.64 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 20, n_between = 90
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.81 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.53 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.94 ± 0.01
#> Random: Slopes Test     0.91 ± 0.01
#> 
#> ── global icc = 0.25, n_within = 20, n_between = 90
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.86 ± 0.02
#> Fixed:  cgm(z1)         0.99 ± 0.00
#> Fixed:  cgm(z2)         0.23 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.95 ± 0.01
#> Random: Slopes Test     0.98 ± 0.01
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 90
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.93 ± 0.01
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.65 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.97 ± 0.01
#> Random: Slopes Test     0.99 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 90
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.96 ± 0.01
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.25 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.98 ± 0.01
#> Random: Slopes Test     1.00 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 10, n_between = 120
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.67 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.47 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.86 ± 0.02
#> Random: Slopes Test     0.58 ± 0.02
#> 
#> ── global icc = 0.25, n_within = 10, n_between = 120
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.69 ± 0.02
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.26 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.91 ± 0.01
#> Random: Slopes Test     0.74 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 20, n_between = 120
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.91 ± 0.01
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.66 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.97 ± 0.01
#> Random: Slopes Test     0.97 ± 0.01
#> 
#> ── global icc = 0.25, n_within = 20, n_between = 120
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.94 ± 0.01
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.30 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.98 ± 0.01
#> Random: Slopes Test     1.00 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 120
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.98 ± 0.01
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.77 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.99 ± 0.00
#> Random: Slopes Test     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 120
#>                         Power      
#> Fixed:  (Intercept)     1.00 ± 0.00
#> Fixed:  cwc(x1)         1.00 ± 0.00
#> Fixed:  cgm(x2)         0.99 ± 0.00
#> Fixed:  cgm(z1)         1.00 ± 0.00
#> Fixed:  cgm(z2)         0.31 ± 0.02
#> Fixed:  cwc(x1):cgm(z1) 0.99 ± 0.00
#> Random: Slopes Test     1.00 ± 0.00
#> 
#> ℹ Margin of error (MOE) computed as ± 1.96 * standard error

Illustration 2: Growth Curve Power Analysis

The second illustration demonstrates a power analysis for a longitudinal growth curve model with a pair of cross-level interactions involving a binary level-2 moderator. Intraclass correlations for longitudinal and intensive repeated measures data often reach values of .40 or higher (Arend & Schäfer, 2019; Bolger & Laurenceau, 2013; Singer & Willett, 2003). To accommodate uncertainty about this important parameter, the simulation investigates intraclass correlation values of .40 and .60.

The multilevel model for the illustration is

\[ \begin{split} Y_{ij} &= \left( \beta_{0} + b_{0j} \right) + \left( \beta_{1} + b_{1j} \right)X_{1ij}^{w} + \beta_{2}\left( X_{2ij}^{w} + X_{2j}^{b} \right) + \beta_{3}\left( X_{3ij}^{w} + X_{3j}^{b} \right)\\ &\phantom{=}+ \ \beta_{4}Z_{1j}^{b} + \ \beta_{5}Z_{2j}^{b} + \beta_{6}Z_{3j}^{b} + \beta_{7}X_{1ij}^{w}Z_{1j}^{b} + \varepsilon_{ij} \\ \mathbf{b}_{j}\ &\sim\ N\left( 0,\mathbf{\Sigma}_{\mathbf{b}} \right)\ \ \ \ \ \mathbf{\Sigma}_{\mathbf{b}}=\begin{pmatrix} \sigma_{b_{0}}^{2} & \\ \sigma_{b_{1},b_{0}} & \sigma_{b_{1}}^{2} \\ \end{pmatrix}\mathbf{\ \ \ \ \ }\varepsilon_{ij}\ \sim\ N(0,\sigma_{\varepsilon}^{2}) \end{split} \]

Following established notation (see Illustration 1), the w and b superscripts reference within- and between-cluster deviation scores, respectively. The explanatory variables include a time score predictor with a random coefficient, a pair of time-varying covariates, a binary level-2 moderator, a pair of level-2 covariates, and a cross-level (group-by-time) interaction. For the purposes of weighting, we designated \(X_{1}\) and \(Z_{1}\) (the interacting variables) as focal predictors, \(X_{2}\) and \(X_{3}\) as level-1 (time-varying) covariates, and \(Z_{2}\) and \(Z_{3}\) as level-2 covariates. To mimic the scaling of a typical temporal index, we assume the major time increments are coded \(X_{1}^{w}\) = (0, 1, 2, 3, 4).

A brief discussion of \(Z_{1}\) (the binary moderator) is warranted before continuing. First, like other level-2 variables, this variable’s population mean must be 0 (i.e., \(Z_{1}\) is centered in the population) in order to maintain the orthogonality of the cross-level interaction terms. Grand mean centering a level-2 binary predictor creates an ANOVA-like contrast code, such that \(\beta_{0}\) is the grand mean and \(\beta_{4}\) is the predicted group mean difference when \(X_{1}^{w}\) (the time score predictor) equals 0. In this case, a code of 0 corresponds to the baseline assessment. Because \(\beta_{4}\) is a conditional effect that represents the group mean difference at the first occasion, \(Z_{1}\)’s contribution to the between-cluster effect size depends on whether we view this regressor as a naturally occurring classification or an intervention assignment indicator. In the former case, we might expect a group mean difference at baseline, and the presence of such a difference would require a non-zero weight. In contrast, random assignment to conditions would eliminate a group mean difference at baseline, and \(Z_{1}\)’s weight would equal 0. Note that this conclusion changes if the time scores are coded differently. For illustration purposes, we assume \(Z_{1}\) is an intervention assignment indicator.

The key inputs to the model object are the unconditional intraclass correlation and effect size values. This illustration investigates power for intraclass correlation values of .40 and .60. The within- and between-cluster fixed effects are set to \(R_{w}^{2}\) = .13 and \(R_{b}^{2}\) = .065, respectively; the former corresponds to Cohen’s (1988) medium effect size benchmark. Note that \(R_{b}^{2}\) cannot exceed the outcome’s intraclass correlation, as this would imply that between-cluster predictors explain more than 100% of the available variance. Following conventional wisdom that interaction effects tend to produce small effects (Aguinis et al., 2005; Chaplin, 1991), \(R_{p}^{2}\) = .05 is assigned to the pair of cross-level product terms. Finally, the random coefficient effect size is set to \(R_{rc}^{2}\) = .03 based on values from Table 2 in Enders et al. (2023).

To refresh, we designated \(X_{1}\) and \(Z_{1}\) (the interacting variables) as focal predictors, \(X_{2}\) and \(X_{3}\) as level-1 (time-varying) covariates, and \(Z_{2}\) and \(Z_{3}\) as level-2 covariates. To mimic a situation where the linear change predominates the level-1 model, we used weights of .50, .25, .25 to allocate the within-cluster \(R^{2}\) to \(X_{1}\), \(X_{2}\), and \(X_{3}\). At level-2, we used weights equal to 0, .50, .50 to allocate the between-cluster \(R^{2}\) to \(Z_{1}\), \(Z_{2}\), and \(Z_{3}\). As noted previously, \(Z_{1}\)’s slope represents the group mean difference at baseline, and we are assuming that random assignment nullifies this effect. Finally, \(R_{p}^{2}\) and \(R_{rc}^{2}\) do not require weights because a single variable determines each source of variation.

To illustrate how to modify default correlation settings, we sampled within-cluster predictor correlations between the range of .20 and .40 under the assumption that the time scores could have stronger associations with other time-varying predictors. We similarly sampled the random effect correlations between the range of .30 and .50 to mimic a scenario where higher baseline scores (i.e., random intercepts) are associated with higher (more positive) growth rates. Finally, we adopted the default correlation range for the between-cluster predictors. The simulations from the previous cross-sectional example suggest that the correlations are somewhat arbitrary and would not have a material impact on power estimates.

The following code block shows the mlmpower model object for this example.

example2 <- (
    effect_size(
        icc          = c(.40, .60),
        within       = .13,
        between      = .065,
        product      = .03,
        random_slope = .10
    )
    + outcome('y', mean = 50, sd = 10)
    + within_time_predictor('x1', weight = .50, values = 0:4)
    + within_predictor('x2', weight = .25)
    + within_predictor('x3', weight = .25)
    + between_binary_predictor('z1', proportion = .50, weight = 0)
    + between_predictor('z2', weight = .50)
    + between_predictor('z3', weight = .50)
    + product('x1','z1', weight = 1)
    + random_slope('x1', weight = 1)
    + correlations(
        within  = random(.20, .40),
        between = random(.10, .30),
        randeff = random(.30, .50)
    )
)

Variable attributes are referenced by adding the following base objects: outcome(), within_predictor(), between_predictor(), product(), and random_slope(). All five objects have an argument for the variable name, weight (weight =), mean (mean =), and standard deviation (sd =). Level-1 variables additionally have an intraclass correlation argument (icc =) that supersedes the global setting in effect_size().

This illustration additionally uses the within_time_predictor() object to specify a set of fixed time scores for \(X_{1}\), and it uses between_binary_predictor() object to define the level-2 moderator \(Z_{1}\) as a binary predictor. In addition to a name and weight, the within_time_predictor() object requires a vector of time scores as an argument (values =). The between_binary_predictor() object requires a name, weight, and the highest category proportion (proportion =).

First, the within_time_predictor() object specifies \(X_{1}\) as a temporal predictor with the fixed set of time scores described earlier. The values = 0:4 argument specifies an integer sequence, but unequally spaced increments can also be specified using a vector as input (e.g., values = c(0,1,3,6)). This object does not require a mean or standard deviation argument, as these quantities are determined from the time scores. Additionally, the variable’s intraclass correlation is automatically fixed to 0 because the time scores are constant across level-2 units. Next the between_binary_predictor('z1', proportion = .50, weight = 0) object specifies \(Z_{1}\) (the moderator variable) as a binary predictor with a 50/50 split. This object does not require a mean or standard deviation argument, as the category proportions determine these quantities. Finally, within the exception of the outcome, the code block uses default values for all means and standard deviations (0 and 1, respectively). The means and standard deviations of the time scores and binary predictor are automatically determined by the user inputs.

The multilevel model parameters also require three sets of correlations: within-cluster correlations among level-1 predictors, between-cluster correlations among level-2 predictors, and random effect correlations. These correlations are optional, but they can specified using the correlations() object. This example samples within-cluster predictor correlation values between .20 and .40 (e.g., to mimic a situation where the time scores have salient correlations with other repeated measures predictors). Random effect correlation values are similarly sampled between the range of .30 and .50 to mimic a scenario where higher baseline scores (i.e., random intercepts) are associated with higher (more positive) growth rates. Finally, the script specifies the default setting for between-cluster correlations, which is to iteratively sample correlation values between .10 and .30.

It may be instructive to inspect the population parameters prior to running the simulation. Executing summary(example2) returns tabular summaries of the multilevel model parameters.

summary(example2)
#> ! Parameters are solved based on average correlation.
#> $`icc = 0.4`
#> $r2
#>                                          Proportion
#> Variance Within-Cluster Fixed Effects     0.1300000
#> Variance Cross-Level Interaction Effects  0.0300000
#> Variance Random Slopes                    0.1000000
#> Within-Cluster Error Variance             0.3400000
#> Variance Between-Cluster Fixed Effects    0.1247412
#> Incremental Variance Level-2 Predictors   0.0650000
#> Between Variance Level-1 Covariates       0.0597412
#> Variance Random Intercepts                0.2752588
#> 
#> $phi_b
#>                 x2_b       x3_b  z1_binary         z2         z3
#> x2_b      0.40000000 0.12000000 0.06324555 0.12649111 0.12649111
#> x3_b      0.12000000 0.40000000 0.06324555 0.12649111 0.12649111
#> z1_binary 0.06324555 0.06324555 0.25000000 0.07978846 0.07978846
#> z2        0.12649111 0.12649111 0.07978846 1.00000000 0.20000000
#> z3        0.12649111 0.12649111 0.07978846 0.20000000 1.00000000
#> 
#> $phi_p
#>                x1_w*z1_binary
#> x1_w*z1_binary            0.5
#> 
#> $phi_w
#>           x1_w      x2_w      x3_w
#> x1_w 2.0000000 0.3286335 0.3286335
#> x2_w 0.3286335 0.6000000 0.1800000
#> x3_w 0.3286335 0.1800000 0.6000000
#> 
#> $mean_Z
#>  z1  z2  z3 
#> 0.5 0.0 0.0 
#> 
#> $mean_X
#> x1 x2 x3 
#>  2  0  0 
#> 
#> $var_e
#>           Value
#> Res. Var.    34
#> 
#> $tau
#>          Icept     x1_w
#> Icept 2.208745 1.329284
#> x1_w  1.329284 5.000000
#> 
#> $gammas
#>                    Value
#> Icept          46.600654
#> x1_w            1.699673
#> x2_w            1.551582
#> x3_w            1.551582
#> x1_w*z1_binary  2.449490
#> x2_b            1.551582
#> x3_b            1.551582
#> z2              1.737198
#> z3              1.737198
#> 
#> $mean_Y
#>  y 
#> 50 
#> 
#> 
#> $`icc = 0.6`
#> $r2
#>                                          Proportion
#> Variance Within-Cluster Fixed Effects     0.1300000
#> Variance Cross-Level Interaction Effects  0.0300000
#> Variance Random Slopes                    0.1000000
#> Within-Cluster Error Variance             0.1400000
#> Variance Between-Cluster Fixed Effects    0.1696753
#> Incremental Variance Level-2 Predictors   0.0650000
#> Between Variance Level-1 Covariates       0.1046753
#> Variance Random Intercepts                0.4303247
#> 
#> $phi_b
#>                 x2_b       x3_b  z1_binary         z2         z3
#> x2_b      0.60000000 0.18000000 0.07745967 0.15491933 0.15491933
#> x3_b      0.18000000 0.60000000 0.07745967 0.15491933 0.15491933
#> z1_binary 0.07745967 0.07745967 0.25000000 0.07978846 0.07978846
#> z2        0.15491933 0.15491933 0.07978846 1.00000000 0.20000000
#> z3        0.15491933 0.15491933 0.07978846 0.20000000 1.00000000
#> 
#> $phi_p
#>                x1_w*z1_binary
#> x1_w*z1_binary            0.5
#> 
#> $phi_w
#>           x1_w      x2_w      x3_w
#> x1_w 2.0000000 0.2683282 0.2683282
#> x2_w 0.2683282 0.4000000 0.1200000
#> x3_w 0.2683282 0.1200000 0.4000000
#> 
#> $mean_Z
#>  z1  z2  z3 
#> 0.5 0.0 0.0 
#> 
#> $mean_X
#> x1 x2 x3 
#>  2  0  0 
#> 
#> $var_e
#>           Value
#> Res. Var.    14
#> 
#> $tau
#>           Icept     x1_w
#> Icept 11.108290 2.981045
#> x1_w   2.981045 5.000000
#> 
#> $gammas
#>                    Value
#> Icept          46.600654
#> x1_w            1.699673
#> x2_w            1.900292
#> x3_w            1.900292
#> x1_w*z1_binary  2.449490
#> x2_b            1.900292
#> x3_b            1.900292
#> z2              1.737198
#> z3              1.737198
#> 
#> $mean_Y
#>  y 
#> 50

Having specified the target model, you next use the power_analysis() function to conduct simulations. The function requires four inputs: the model argument specifies the parameter value object (e.g., example2), the replications input specifies the number of artificial data sets, n_between is a vector of level-2 sample sizes, and n_within is a vector of level-1 sample sizes (i.e., the number of observations per cluster). The code block below specifies six level-2 sample size conditions (\(J =\) 50, 60, 70, 80, 90, and 100), each with five repeated measurements and 2,000 replications.

# Set seed for replicable results
set.seed(12379)

# Run Power Analysis
powersim2 <-
    power_analysis(
        model = example2,
        replications = 2000,
        n_between = c(50, 60, 70, 80, 90, 100),
        n_within = 5
    )

The package uses lme4 (Bates et al., 2021) for model fitting, and it defaults to an alpha level of .05 for all significance tests. Significance tests of random slope variation use a likelihood ratio test with a mixture chi-square reference distribution (i.e., a chi-bar distribution; Snijders & Bosker, 2012, p. 99), as implemented in the varTestnlme package (Baey & Kuhn, 2022). Executing summary(powersim2) returns the tabular summaries of the simulation results shown below.

summary(powersim2)
#> 
#> ── Power Analysis Specification ──
#> 
#> Replications = 2000
#> N_Within = 5
#> N_Between = 50, 60, 70, 80, 90, and 100
#> 
#> ── lme4 Model Syntax
#> 
#> y ~ 1 + x1 + cgm(x2) + cgm(x3) + x1:cgm(z1) + cgm(z1) + cgm(z2) + cgm(z3) + (1
#> + x1 | `_id`)
#> 
#> ℹ cwc() = centering within cluster
#> ℹ cgm() = centering with grand mean
#> 
#> ── Effect Sizes
#> 
#> • WITHIN = 0.13
#> • BETWEEN = 0.065
#> • RANDOM SLOPE = 0.1
#> • PRODUCT = 0.03
#> 
#> ── Power Analysis Results ──
#> 
#> ── global icc = 0.4, n_within = 5, n_between = 50
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          0.97 ± 0.01
#> Fixed:  cgm(x2)     0.84 ± 0.02
#> Fixed:  cgm(x3)     0.84 ± 0.02
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.69 ± 0.02
#> Fixed:  cgm(z3)     0.67 ± 0.02
#> Fixed:  x1:cgm(z1)  0.82 ± 0.02
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.6, n_within = 5, n_between = 50
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          0.99 ± 0.00
#> Fixed:  cgm(x2)     0.98 ± 0.01
#> Fixed:  cgm(x3)     0.99 ± 0.00
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.61 ± 0.02
#> Fixed:  cgm(z3)     0.62 ± 0.02
#> Fixed:  x1:cgm(z1)  0.92 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.4, n_within = 5, n_between = 60
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          0.98 ± 0.01
#> Fixed:  cgm(x2)     0.90 ± 0.01
#> Fixed:  cgm(x3)     0.89 ± 0.01
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.75 ± 0.02
#> Fixed:  cgm(z3)     0.76 ± 0.02
#> Fixed:  x1:cgm(z1)  0.89 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.6, n_within = 5, n_between = 60
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          0.99 ± 0.00
#> Fixed:  cgm(x2)     1.00 ± 0.00
#> Fixed:  cgm(x3)     1.00 ± 0.00
#> Fixed:  cgm(z1)     0.06 ± 0.01
#> Fixed:  cgm(z2)     0.72 ± 0.02
#> Fixed:  cgm(z3)     0.69 ± 0.02
#> Fixed:  x1:cgm(z1)  0.96 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.4, n_within = 5, n_between = 70
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          0.99 ± 0.00
#> Fixed:  cgm(x2)     0.94 ± 0.01
#> Fixed:  cgm(x3)     0.94 ± 0.01
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.82 ± 0.02
#> Fixed:  cgm(z3)     0.81 ± 0.02
#> Fixed:  x1:cgm(z1)  0.94 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.6, n_within = 5, n_between = 70
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     1.00 ± 0.00
#> Fixed:  cgm(x3)     1.00 ± 0.00
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.79 ± 0.02
#> Fixed:  cgm(z3)     0.79 ± 0.02
#> Fixed:  x1:cgm(z1)  0.98 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.4, n_within = 5, n_between = 80
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     0.96 ± 0.01
#> Fixed:  cgm(x3)     0.97 ± 0.01
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.88 ± 0.01
#> Fixed:  cgm(z3)     0.86 ± 0.02
#> Fixed:  x1:cgm(z1)  0.96 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.6, n_within = 5, n_between = 80
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     1.00 ± 0.00
#> Fixed:  cgm(x3)     1.00 ± 0.00
#> Fixed:  cgm(z1)     0.04 ± 0.01
#> Fixed:  cgm(z2)     0.83 ± 0.02
#> Fixed:  cgm(z3)     0.84 ± 0.02
#> Fixed:  x1:cgm(z1)  0.99 ± 0.00
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.4, n_within = 5, n_between = 90
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     0.98 ± 0.01
#> Fixed:  cgm(x3)     0.97 ± 0.01
#> Fixed:  cgm(z1)     0.06 ± 0.01
#> Fixed:  cgm(z2)     0.90 ± 0.01
#> Fixed:  cgm(z3)     0.91 ± 0.01
#> Fixed:  x1:cgm(z1)  0.98 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.6, n_within = 5, n_between = 90
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     1.00 ± 0.00
#> Fixed:  cgm(x3)     1.00 ± 0.00
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.87 ± 0.01
#> Fixed:  cgm(z3)     0.88 ± 0.01
#> Fixed:  x1:cgm(z1)  1.00 ± 0.00
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.4, n_within = 5, n_between = 100
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     0.99 ± 0.00
#> Fixed:  cgm(x3)     0.99 ± 0.00
#> Fixed:  cgm(z1)     0.04 ± 0.01
#> Fixed:  cgm(z2)     0.95 ± 0.01
#> Fixed:  cgm(z3)     0.93 ± 0.01
#> Fixed:  x1:cgm(z1)  0.98 ± 0.01
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ── global icc = 0.6, n_within = 5, n_between = 100
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  x1          1.00 ± 0.00
#> Fixed:  cgm(x2)     1.00 ± 0.00
#> Fixed:  cgm(x3)     1.00 ± 0.00
#> Fixed:  cgm(z1)     0.05 ± 0.01
#> Fixed:  cgm(z2)     0.92 ± 0.01
#> Fixed:  cgm(z3)     0.90 ± 0.01
#> Fixed:  x1:cgm(z1)  1.00 ± 0.00
#> Random: Slopes Test 1.00 ± 0.00
#> 
#> ℹ Margin of error (MOE) computed as ± 1.96 * standard error

Illustration 3 Vignette: Cluster-Randomized Design

The third vignette demonstrates a power analysis for a cluster-randomized design (Raudenbush, 1997) where level-2 units are randomly assigned to one of two experimental conditions. To provide a substantive context, consider a prototypical education example where students are nested in schools, and schools are the unit of randomization. Intraclass correlations for achievement-related outcomes often range between .10 and .25 (Hedges & Hedberg, 2007; Hedges & Hedberg, 2013; Sellström & Bremberg, 2006; Spybrook et al., 2011; Stockford, 2009). To accommodate uncertainty about this important parameter, the simulation investigates intraclass correlation values of .10 and .25.

The multilevel model for the illustration is

\[ \begin{split} Y_{ij} &= \beta_{0} + \beta_{1}\left( X_{1ij}^{w} + X_{1j}^{b} \right) + \beta_{2}\left( X_{2ij}^{w} + X_{2j}^{b} \right)\\ &\phantom{=}+ \ \beta_{3}\left( X_{3ij}^{w} + X_{3j}^{b} \right) + \beta_{4}\left( X_{4ij}^{w} + X_{4j}^{b} \right) + \ \beta_{5}Z_{1j}^{b} + b_{0j} + \varepsilon_{ij} \end{split} \]

where \(X_{1}\) to \(X_{4}\) are grand mean centered level-1 covariates, and \(Z_{1}\) is a binary intervention assignment indicator. Following established notation, the w and b superscripts reference within- and between-cluster deviation scores, respectively. The notation conveys that all level-1 predictors contain within- and between-cluster variation. By default, all predictors are centered, including the binary dummy code. Grand mean centering a level-2 binary predictor creates an ANOVA-like contrast code, such that \(\beta_{0}\) is the grand mean and \(\beta_{5}\) is the predicted mean difference, adjusted for the covariates. Turning to the residual terms, \(b_{0j}\) is a between-cluster random effect that captures residual variation in the cluster-specific intercepts, and \(\varepsilon_{ij}\) is a within-cluster residual. By assumption, both residuals are normal with constant variance.

The key inputs to the model object are the unconditional intraclass correlation and effect size values. This illustration investigates power for intraclass correlation values of .10 and .25. To mimic a scenario with a strong covariate set (e.g., one that includes a pretest measure of the outcome), the within-cluster effect size is set to \(R_{w}^{2}\) = .18 (a value roughly midway between Cohen’s small and medium benchmarks). Most of this variation was allocated to \(X_{1}\) (the pretest) by assigning it a weight of .70, and the remaining three predictors had their within-cluster weights set to .10. We previously argued that the allocation of the weights among covariates is arbitrary because these predictors are not the focus. To demonstrate that point, we conducted a second simulation that weighted the four level-1 covariates equally.

Turning to the level-2 predictor, researchers often prefer the Cohen’s (1988) d effect size when working with binary explanatory variables. To illustrate, consider power for d = .40, the midway point between Cohen’s small and medium benchmarks. Substituting this value into the conversion equation below returns \(R_{b}^{2}\) = .038.

\[R^{2} = \frac{d^{2}}{d^{2} + 4}\]

The following code block shows the mlmpower model object for the example.

example3 <- (
    effect_size(
        icc     = c(.10, .25),
        within  = .18,
        between = .038,
    )
    + outcome('y')
    + within_predictor('x1', weight = .70)
    + within_predictor('x2', weight = .10)
    + within_predictor('x3', weight = .10)
    + within_predictor('x4', weight = .10)
    + between_binary_predictor('z1', proportion = .50, weight = 1)
)

The mlmpower package groups \(R^{2}\) values and intraclass correlations into a single object called effect_size(), with three inputs separated by commas. The within argument corresponds to \(R_{w}^{2}\), the between argument aligns with \(R_{b}^{2}\), and the icc argument assigns a global intraclass correlation to all level-1 variables. Separate simulations are performed for each requested level.

Variable attributes for this example require three objects: outcome(), within_predictor(), and between_binary_predictor(). The first two objects have an argument for the variable name, weight (weight =), mean (mean =), standard deviation (sd =), and an intraclass correlation argument (icc =) that supersedes the global setting in effect_size(). The between_binary_predictor() object requires a name, weight, and the highest category proportion (proportion =).

The previous code block assigns explicit weights to all variables contributing to a given effect, as described previously. The unit weight in the between_binary_predictor() object reflects the fact that this \(Z_{1}\) solely determines \(R_{b}^{2}\). Finally, the proportion argument assigns a 50/50 split to the level-2 groups.

The multilevel model parameters also require two sets of correlations: within-cluster correlations among level-1 predictors, and between-cluster correlations among level-2 predictors (in this case, the cluster means and the intervention assignment indicator). These correlations are optional, but they can specified using the correlations() object. The earlier code block omits this object, thereby accepting the default specification. When the correlations() object is omitted, the mlmpower package iteratively samples all correlations from a uniform distribution between .10 and .30, such that the resulting power estimates average over a distribution of possible associations. This default range spans Cohen’s (1988) small to medium effect size benchmarks, and it also brackets common correlations from published research (Bosco et al., 2015; Funder & Ozer, 2019; Gignac & Szodorai, 2016). The default specifications could be invoked by explicitly appending the correlations() object to the earlier script, as follows.

example3 <- (
    example3
    + correlations(
        within  = random(0.1, 0.3),
        between = random(0.1, 0.3)
    )
)

Researchers can modify the upper and lower limits of each correlation range, or they can specify a constant correlation by inputting a scalar value. The illustrative simulations presented in Enders et al. (2023) suggest that predictor correlations tend not to matter very much.

It may be instructive to inspect the population parameters prior to running the simulation. Executing summary(example3) returns tabular summaries of the multilevel model parameters.

summary(example3)
#> ! Parameters are solved based on average correlation.
#> $`icc = 0.1`
#> $r2
#>                                          Proportion
#> Variance Within-Cluster Fixed Effects    0.18000000
#> Variance Cross-Level Interaction Effects 0.00000000
#> Variance Random Slopes                   0.00000000
#> Within-Cluster Error Variance            0.72000000
#> Variance Between-Cluster Fixed Effects   0.07703223
#> Incremental Variance Level-2 Predictors  0.03800000
#> Between Variance Level-1 Covariates      0.03903223
#> Variance Random Intercepts               0.02296777
#> 
#> $phi_b
#>                 x1_b       x2_b       x3_b       x4_b  z1_binary
#> x1_b      0.10000000 0.02000000 0.02000000 0.02000000 0.03162278
#> x2_b      0.02000000 0.10000000 0.02000000 0.02000000 0.03162278
#> x3_b      0.02000000 0.02000000 0.10000000 0.02000000 0.03162278
#> x4_b      0.02000000 0.02000000 0.02000000 0.10000000 0.03162278
#> z1_binary 0.03162278 0.03162278 0.03162278 0.03162278 0.25000000
#> 
#> $phi_p
#> <0 x 0 matrix>
#> 
#> $phi_w
#>      x1_w x2_w x3_w x4_w
#> x1_w 0.90 0.18 0.18 0.18
#> x2_w 0.18 0.90 0.18 0.18
#> x3_w 0.18 0.18 0.90 0.18
#> x4_w 0.18 0.18 0.18 0.90
#> 
#> $mean_Z
#>  z1 
#> 0.5 
#> 
#> $mean_X
#> x1 x2 x3 x4 
#>  0  0  0  0 
#> 
#> $var_e
#>           Value
#> Res. Var.    18
#> 
#> $tau
#>           Icept
#> Icept 0.5741943
#> 
#> $gammas
#>                Value
#> Icept     10.0000000
#> x1_w       1.9943101
#> x2_w       0.2849014
#> x3_w       0.2849014
#> x4_w       0.2849014
#> x1_b       1.9943101
#> x2_b       0.2849014
#> x3_b       0.2849014
#> x4_b       0.2849014
#> z1_binary  2.0548047
#> 
#> $mean_Y
#>  y 
#> 10 
#> 
#> 
#> $`icc = 0.25`
#> $r2
#>                                          Proportion
#> Variance Within-Cluster Fixed Effects     0.1800000
#> Variance Cross-Level Interaction Effects  0.0000000
#> Variance Random Slopes                    0.0000000
#> Within-Cluster Error Variance             0.5700000
#> Variance Between-Cluster Fixed Effects    0.1278739
#> Incremental Variance Level-2 Predictors   0.0380000
#> Between Variance Level-1 Covariates       0.0898739
#> Variance Random Intercepts                0.1221261
#> 
#> $phi_b
#>           x1_b x2_b x3_b x4_b z1_binary
#> x1_b      0.25 0.05 0.05 0.05      0.05
#> x2_b      0.05 0.25 0.05 0.05      0.05
#> x3_b      0.05 0.05 0.25 0.05      0.05
#> x4_b      0.05 0.05 0.05 0.25      0.05
#> z1_binary 0.05 0.05 0.05 0.05      0.25
#> 
#> $phi_p
#> <0 x 0 matrix>
#> 
#> $phi_w
#>      x1_w x2_w x3_w x4_w
#> x1_w 0.75 0.15 0.15 0.15
#> x2_w 0.15 0.75 0.15 0.15
#> x3_w 0.15 0.15 0.75 0.15
#> x4_w 0.15 0.15 0.15 0.75
#> 
#> $mean_Z
#>  z1 
#> 0.5 
#> 
#> $mean_X
#> x1 x2 x3 x4 
#>  0  0  0  0 
#> 
#> $var_e
#>           Value
#> Res. Var. 14.25
#> 
#> $tau
#>          Icept
#> Icept 3.053152
#> 
#> $gammas
#>                Value
#> Icept     10.0000000
#> x1_w       2.1846572
#> x2_w       0.3120939
#> x3_w       0.3120939
#> x4_w       0.3120939
#> x1_b       2.1846572
#> x2_b       0.3120939
#> x3_b       0.3120939
#> x4_b       0.3120939
#> z1_binary  2.0548047
#> 
#> $mean_Y
#>  y 
#> 10

Having specified the target model, you next use the power_analysis() function to conduct the simulations. The function requires four inputs: the model argument specifies the parameter value object (e.g., example3), the replications input specifies the number of artificial data sets, n_between is a vector of level-2 sample sizes, and n_within is a vector of level-1 sample sizes (i.e., the number of observations per cluster). The code block below pairs four level-2 sample size values (\(J =\) 30, 60, 90, and 120) with two level-1 sample sizes (\(n_{j} =\) 15 or 30 observations per cluster), and it requests 2,000 artificial data sets for each combination.

The package uses lme4 (Bates et al., 2021) for model fitting, and it defaults to an alpha level of .05 for all significance tests. Executing summary(powersim) returns the tabular summaries of the simulation results shown below.

# Set seed for replicable results
set.seed(981723)

# Run Power Analysis
powersim3 <-
    power_analysis(
        model = example3,
        replications = 2000,
        n_between = c(30, 60, 90, 120),
        n_within = c(15, 30)
    )

The package uses lme4 (Bates et al., 2021) for model fitting, and it defaults to an alpha level of .05 for all significance tests. Executing summary(powersim3) returns the tabular summaries of the simulation results shown below.

summary(powersim3)
#> 
#> ── Power Analysis Specification ──
#> 
#> Replications = 2000
#> N_Within = 15 and 30
#> N_Between = 30, 60, 90, and 120
#> 
#> ── lme4 Model Syntax
#> 
#> y ~ 1 + cgm(x1) + cgm(x2) + cgm(x3) + cgm(x4) + cgm(z1) + (1 | `_id`)
#> 
#> ℹ cwc() = centering within cluster
#> ℹ cgm() = centering with grand mean
#> 
#> ── Effect Sizes
#> 
#> • WITHIN = 0.18
#> • BETWEEN = 0.038
#> • RANDOM SLOPE = 0
#> • PRODUCT = 0
#> 
#> ── Power Analysis Results ──
#> 
#> ── global icc = 0.1, n_within = 15, n_between = 30
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.27 ± 0.02
#> Fixed:  cgm(x3)     0.25 ± 0.02
#> Fixed:  cgm(x4)     0.26 ± 0.02
#> Fixed:  cgm(z1)     0.98 ± 0.01
#> 
#> ── global icc = 0.25, n_within = 15, n_between = 30
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.30 ± 0.02
#> Fixed:  cgm(x3)     0.31 ± 0.02
#> Fixed:  cgm(x4)     0.31 ± 0.02
#> Fixed:  cgm(z1)     0.76 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 30
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.45 ± 0.02
#> Fixed:  cgm(x3)     0.46 ± 0.02
#> Fixed:  cgm(x4)     0.48 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 30
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.54 ± 0.02
#> Fixed:  cgm(x3)     0.54 ± 0.02
#> Fixed:  cgm(x4)     0.53 ± 0.02
#> Fixed:  cgm(z1)     0.79 ± 0.02
#> 
#> ── global icc = 0.1, n_within = 15, n_between = 60
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.46 ± 0.02
#> Fixed:  cgm(x3)     0.48 ± 0.02
#> Fixed:  cgm(x4)     0.47 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 15, n_between = 60
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.54 ± 0.02
#> Fixed:  cgm(x3)     0.53 ± 0.02
#> Fixed:  cgm(x4)     0.55 ± 0.02
#> Fixed:  cgm(z1)     0.98 ± 0.01
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 60
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.76 ± 0.02
#> Fixed:  cgm(x3)     0.75 ± 0.02
#> Fixed:  cgm(x4)     0.75 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 60
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.83 ± 0.02
#> Fixed:  cgm(x3)     0.83 ± 0.02
#> Fixed:  cgm(x4)     0.84 ± 0.02
#> Fixed:  cgm(z1)     0.99 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 15, n_between = 90
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.64 ± 0.02
#> Fixed:  cgm(x3)     0.63 ± 0.02
#> Fixed:  cgm(x4)     0.63 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 15, n_between = 90
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.72 ± 0.02
#> Fixed:  cgm(x3)     0.71 ± 0.02
#> Fixed:  cgm(x4)     0.71 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 90
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.90 ± 0.01
#> Fixed:  cgm(x3)     0.90 ± 0.01
#> Fixed:  cgm(x4)     0.89 ± 0.01
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 90
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.95 ± 0.01
#> Fixed:  cgm(x3)     0.95 ± 0.01
#> Fixed:  cgm(x4)     0.95 ± 0.01
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 15, n_between = 120
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.75 ± 0.02
#> Fixed:  cgm(x3)     0.75 ± 0.02
#> Fixed:  cgm(x4)     0.77 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 15, n_between = 120
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.84 ± 0.02
#> Fixed:  cgm(x3)     0.83 ± 0.02
#> Fixed:  cgm(x4)     0.84 ± 0.02
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.1, n_within = 30, n_between = 120
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.96 ± 0.01
#> Fixed:  cgm(x3)     0.96 ± 0.01
#> Fixed:  cgm(x4)     0.97 ± 0.01
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ── global icc = 0.25, n_within = 30, n_between = 120
#>                     Power      
#> Fixed:  (Intercept) 1.00 ± 0.00
#> Fixed:  cgm(x1)     1.00 ± 0.00
#> Fixed:  cgm(x2)     0.98 ± 0.01
#> Fixed:  cgm(x3)     0.99 ± 0.01
#> Fixed:  cgm(x4)     0.98 ± 0.01
#> Fixed:  cgm(z1)     1.00 ± 0.00
#> 
#> ℹ Margin of error (MOE) computed as ± 1.96 * standard error

Using `mlmpower` Package to Conduct Multilevel Power Analysis

Brian T. Keller

Illustration 1: Cross-Sectional Power Analysis

Illustration 2: Growth Curve Power Analysis

Illustration 3 Vignette: Cluster-Randomized Design