Centering and Statistical Power in Multilevel Modeling

Determining the relationship to estimate and the data scale together

Multilevel modeling and mixed-effects models estimate the relationship of interest while separating variance attributable to groups, participants, and items. Centering establishes the reference point for estimating the relationship, and power analysis examines the combinations of sample sizes and observations needed to detect that relationship.

From group-mean centering to random slopes and cross-level interactions

A tutorial on two-level linear modeling presents a procedure for analyzing data in which individuals belong to groups using SPSS, Stata, R, and Mplus (Sommet & Morselli, 2021). The practice data consist of 175 members of 50 boy bands. They were drawn from the Internet Boy Band Database, which includes bands with songs that entered the US Billboard Hot 100 between 1980 and 2018. The outcome variable is a popularity score derived from log-transformed Instagram follower counts. The individual-level predictor, physical attractiveness, was measured by the number of appearances in internet rankings. The group-level variable, peak activity period, was coded as −0.5 for 1990–2000 and +0.5 for after 2000. The empty model (Y_{ij}=B_{00}+u_{0j}+e_{ij}) separates the grand mean, the group mean’s deviation from the grand mean, and the individual’s deviation from the group mean. Grand-mean centering an individual-level predictor preserves the slope in a model without interactions. Coefficients that do not distinguish within-group from between-group effects combine the two relationships. In the practice exercise, physical attractiveness was group-mean centered to estimate the relationship within bands, and the intercept was interpreted as the predicted popularity score of a member with average physical attractiveness in that band. The fixed slope represents the average effect across all groups. The random slope indicates how much each group’s effect deviates from that average. The intercept-slope covariance expresses the relationship between the intercept and slope deviations.

The ICC for the practice data was 0.82, meaning that 82% of the variance in popularity scores was attributable to differences between bands, and the DEFF, based on an average group size of 3.50 members, was 3.04. The DEFF is calculated as (1+(n-1)\times ICC). A provisional criterion suggested that ordinary regression might be used if the value is below 1.5. To assess random slopes before introducing the interaction, a model that added slope residual variance and intercept-slope covariance was compared with the preceding model. The likelihood ratio test yielded LR χ²(2) = 415.96 − 412.08 = 3.88, p = 0.144. The two terms were retained under the proposed model selection criterion of α = 0.20. In the final model, the coefficient for physical attractiveness was B10 = 0.03, 95% CI [−0.14, 0.21]. The coefficient for peak activity period was B01 = 1.59, 95% CI [0.94, 2.24]. The interaction was B11 = −0.39, 95% CI [−0.73, −0.04]. The simple slopes obtained by recoding the peak activity period were B = 0.23, 95% CI [0.05, 0.41] for 1990s bands and B = −0.16, 95% CI [−0.45, 0.14] for bands after 2000. The results are interpreted as a relationship between the two online indicators used as measures. The discussion of ML and REML explained that REML can be less biased in small samples. Models with different fixed-effects structures cannot be compared using REML.

In the simulation comparing software, each dataset contained 12,500 observations nested within 50 groups. The two main effects were set to 0, the ICC to 0.05, and the intercept-slope covariance to 0. A total of 9,000 datasets, 1,000 per condition, were analyzed across combinations of interaction sizes of 0.10, 0.05, and 0.00 and slope residual variances of 0.01, 0.005, and 0.001. Statistical power, verified using simglm in R, was 0.80 at α = 0.05 with an interaction of 0.10 and a slope residual variance of 0.01. The comparison metrics were the convergence rate within 100 iterations, the detection rate for slope residual variance in a likelihood ratio test at α = 0.20, and the error rate for the interaction test at α = 0.05. When the slope residual variance was 0.001, the convergence rate for SPSS was 88%–90%. The other software showed high convergence rates. The detection rate for residual variance was nearly identical across software. It was approximately 53% when the interaction was 0 and the slope residual variance was 0.001. For Stata, R, and Mplus, Type II error rates were 20%, 6%, and 0% across the residual variance conditions when the interaction was 0.10. They were 70%, 55%, and 23% when the interaction was 0.05. The Type II error rate for SPSS was higher, but the difference between software was not significant. With an interaction of 0, all software yielded Type I error rates of 3%, 3%, and 2% across the conditions. Because the detectability of small slope variations depends on the sample structure and parameter settings, decisions about including random slopes draw on both test results and theoretical expectations about differences in effects across groups. Open data, scripts, and supplementary materials provided guidance on extending the analysis to statistical power, effect sizes, three-level, cross-classified, and multiple-membership designs, and non-continuous outcome variables.

Designing mixed-effects model statistical power through combinations of participants and items

A tutorial on statistical power for mixed-effects models uses simulation to compare sample configurations for detecting fixed effects of interest, with the model and random-effects structure specified in advance (Kumle et al., 2021). Statistical power is 1−β, the probability of rejecting the null hypothesis when a true effect exists. It is estimated as the proportion of significant results obtained by repeatedly generating data, refitting the model, and testing the effect. R packages lme4, mixedpower, and simr were used for the analysis, along with open data and Notebooks. The first scenario derives parameters from a model fitted to existing data. The example was an eye-tracking study in which 48 participants read 120 sentences. The outcome variable was the first landing position (FLP), the location where the eye first lands on a word. Word length, morphological complexity, and their interaction were specified as fixed effects. Random intercepts for participants and sentences were specified as crossed random effects. The coefficients of the model with centered continuous predictors were β = 1.511 for word length, β = −0.075 for morphological complexity, and β = 0.116 for the interaction. Participant sample sizes of 20, 30, 40, 50, and 60 were compared using mixedpower(). New response values were generated using simulate.merMod() from lme4. A t-value of 2 served as an approximate threshold corresponding to a 5% significance level. As the number of participants increased, statistical power for morphological complexity and the interaction increased. Power for word length, which had a large effect, remained high across the sample sizes examined. All coefficients were reduced by 15% to examine the Smallest Effect Size of Interest (SESOI) and uncertainty in effect sizes. Power for morphological complexity and the interaction decreased. The reduction was chosen for illustration.

In the second scenario, the numbers of participants and sentences were varied jointly. With the number of participants held at 48, varying the number of stimuli from the original 120 sentences showed that statistical power increased as sentences were added. R2power() enabled comparisons across combinations of 30, 48, and 60 participants with varying sentence counts. The third scenario constructs artificial data and a model based on prior assumptions. A binomial GLMM was used to analyze the effects of native English speaker status and word frequency on lexical decision accuracy. The initial data had a fully crossed structure in which 20 participants viewed all 100 words. Equal numbers of participants were assigned to native and non-native speaker groups, coded as −0.5 and 0.5. Word frequency was generated from a normal distribution with a mean of 5 and a standard deviation of 1 based on the lexdec data from languageR, and then centered. The model included random intercepts for participants and words. Fixed-effect coefficients, random-effect variances, and a binomial distribution were specified in makeGlmer() from simr. Parameters were obtained from a similar model fitted to the lexdec data. The initial power estimate for the native language status effect, obtained with powerSim(), was 16.7%. Reducing the standard deviation of word frequency to 0.5 produced a difference of up to 24% in power estimates. Analyses that varied group ratios demonstrated that unbalanced designs could reduce statistical power. Setting all effect sizes to 0 allowed the Type I error rate to be estimated using the same procedure.

In the third example, 180 participants were needed to achieve 80% power for all effects with 100 words, but if the effects were 15% smaller, 180 participants were insufficient. The effect of native language status determined the required sample size. The analysis distinguishes between each effect reaching the target power and the probability of multiple effects being significant simultaneously. At least 1,000 Monte Carlo iterations were recommended for the final analysis. Increasing the number of iterations reduces the variance of the estimates but lengthens execution time. The estimates depend on assumptions about the data-generating model, effect sizes, random-effect variances, and predictor distributions. mixedpower is suited to comparing multiple effects and sample conditions through parallel computing. simr allows detailed configuration of artificial models, confidence intervals for power estimates, and significance testing methods. The tutorial presented an approach in which power analysis is conducted a priori using data and models independent of the planned research. This approach avoids using post-hoc power, which has a monotonic relationship with p-values, to interpret results already obtained.

Reading notes

In HRD and medical education research, modeling data structures in which learners belong to organizations, perform tasks, or are observed repeatedly connects choices about centering with sample design. A study in which learners belong to hospitals and perform multiple tasks can compare how power to detect an effect of interest changes as the numbers of hospitals, learners, and tasks are adjusted. Varying the random slope variance or the expected distribution of predictors identifies conditions that affect detectability even when recruitment numbers remain the same.

References

• Sommet, N., & Morselli, D. (2021). Keep calm and learn multilevel linear modeling: A three-step procedure using SPSS, Stata, R, and Mplus. International Review of Social Psychology, 34(1), 24. https://doi.org/10.5334/irsp.555
• 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(6), 2528–2543. https://doi.org/10.3758/s13428-021-01546-0