Specifying Location-Scale Models for Heterogeneous Variances as Multilevel SEMs

Abstract

Standard multilevel models focus on variables that predict the mean while the within-group variability is largely treated as a nuisance. Recent work has shown the advantage of including predictors for both the mean (the location submodel) and the variability (the scale submodel) within a single model. Constrained versions of the model can be fit in standard mixed effect model software, but the most general version with random effects in each of the location and scale submodels has been noted for being difficult to fit and estimate in software. However, the latest release of Mplus includes new capabilities that facilitate fitting the general version of the model as a multilevel structural equation model (SEM). This article introduces the general form of the model that includes location and scale random effects (called the location-scale model) and notes how it can be envisioned as a multilevel SEM. We provide a tutorial with example analyses and Mplus code for the model with two-level cross-sectional data and three-level repeated measures data and discuss how such a model has potential to extend recent developments in organizational science.

Keywords

multilevel model location-scale model Mplus heterogeneous variance model

Due to the hierarchical structure of organizations such that people are clustered within some higher-level unit, multilevel models are commonly applied in organizational sciences (Kozlowski & Klein, 2000). The benefit of multilevel models is that researchers can investigate relationships both within organizational groups and between organizational groups (e.g., Zhang, Zyphur, & Preacher, 2009). Researchers often take advantage of this axiom for modeling how the mean changes as a function of predictor variables, which is referred to as the location aspect of the model. Although location-only models are common, a more general version of the multilevel model allows researchers to also model how variances within groups change as a function of predictors (Bolger, Davis, & Rafaeli, 2003; Cleveland, Denby, & Liu, 2002). For example, it might be of interest to model which factors make a group more homogeneous or cohesive in addition to which factors that cause a group’s mean to increase or decrease.

Recent studies have noted that a joint modeling approach can be useful to simultaneously model how both mean and variances change as a function of predictors (e.g., Brunton-Smith, Sturgis, & Leckie, 2017; Kuppens & Yzerbyt, 2014; Lang, Bliese, & de Voogt, 2018; Lang, Bliese, & Runge, 2019). These types of models are broadly referred to as heterogeneous variance models because the groups have heterogeneous variances, conditional on predictors. For instance, Lang et al. (2018) advanced such a model to study consensus emergence in which groups become more similar over time and the rate at which this occurs varies based on predictors. Heterogeneous variance models are straightforward to apply in standard mixed effect software used for multilevel models like SAS, R, or Stata so long as researchers only want to include fixed effects in the scale submodel, where variance is the outcome (Jongerling, Laurenceau, & Hamaker, 2015). However, in the more likely case where predictors cannot completely account for between-group differences (Jongerling et al., 2015), the more general location-scale model is needed whereby random effects are included in the location submodel to capture each group’s mean (the location) and in the scale submodel to capture each group’s variance (the scale).

Walters, Hoffman, and Templin (2018) noted that there are few software options for location-scale models featuring both location and scale random effects. Software options that do exist tend to be low on user-friendliness. For instance, models with scale random effects can be fit in SAS Proc Nlimixed but require a custom model parameterization (Hedeker, Mermelstein, & Demirtas, 2008), which requires a fairly high level of statistical programming knowledge (Hedeker & Nordgren, 2013). The MIXREGLS freeware program is specifically designed for location-scale models, but it is not designed to accommodate random slopes (Hedeker & Nordgren, 2013, p. 30). Alternatively, the model can be manually programmed as a hierarchical model in general Bayesian software like JAGS or Stan (Rast, Hofer, & Sparks, 2012). Lester, Cullen-Lester, and Walters (2020) ¹ described how to fit location-scale models within these Bayesian packages in the R environment, but support for researchers not well-versed in the R environment is less widespread. So, despite potential advantages shown by articles advocating for modeling heterogeneity with random scale effects, application of the model in software has been a barrier to wider implementation.

Dynamic structural equation models (DSEMs) have recently been advanced in the statistical literature to incorporate multilevel, structural equation, and time-series modeling under one large framework (Asparouhov, Hamaker, & Muthén, 2018). Although this article is not about DSEM (see McNeish & Hamaker, 2020, for details on DSEM), the development of DSEM has spawned advances in multilevel SEM software capabilities. The benefit of these advances is that the recently released Version 8 of Mplus improved functionality of multilevel structural equation models (SEM) such that it can seamlessly incorporate random intercepts and slopes in the scale submodel, meaning that the general location-scale model can be fit within software familiar to organizational researchers with minimal programmatic effort, especially compared to existing alternatives.

In this article, we show how the location-scale model can be conceived as a multilevel SEM to take advantage of these improved software capabilities and to circumvent the computational and programmatic issues with location-scale models that can impede implementation. In this article, we introduce the basics of the location-scale model, discuss how the model can be conceived as a multilevel SEM, and provide Mplus code for fitting the model. We use a two-level cross-sectional running example when introducing the model before shifting to a three-level repeated measures example reminiscent of questions posed by organizational researchers. Ultimately, we hope that showing how the location-scale model can be conceived as a multilevel SEM and fit within Mplus will allow organizational researchers to more easily exploit the advantages of such models recently espoused in the literature.

The Standard Location-Only Multilevel Model

To overview the standard multilevel model, consider the study by Chen et al. (2007), who was interested in modeling individual empowerment (IP) of people nested within teams as a function of the quality of the leader-member exchanges (LMX), for which each person has a unique value (i.e., it is a within-group predictor). In hierarchical notation popularized by Raudenbush and Bryk (2002), the standard multilevel random coefficients model that disaggregates the within-group predictor for this example would be

\begin{array}{l} I P_{i j} = β_{0 j} + β_{1 j} (L M X_{i j} - {\bar{L M X}}_{j}) + r_{i j} \\ β_{0 j} = γ_{00} + γ_{01} {\bar{L M X}}_{j} + u_{0 j} \\ β_{1 j} = γ_{10} + u_{1 j} \end{array},

where $I P_{i j}$ is the individual empowerment outcome variable for the ith person in the jth team and $L M X_{i j}$ is the within-team predictor value for the ith person in the jth team (i.e., each member of the team has a potentially different LMX value). LMX is team-mean-centered in the first expression to avoid a conflated effect (i.e., by centering around the team-mean, there will be separate effects for a one-unit increase in a person’s LMX and a one-unit increase in a team’s LMX).² The intercept ( $β_{0 j}$ ) and the regression coefficient ( $β_{1 j}$ ) have j subscripts, meaning that the intercept and relation between LMX and IP could potentially be different for each team, so the β coefficients become outcomes in team-level equations to model the between-team differences. The $β_{0 j}$ equation states that the intercept in team j is equal to the average intercept across all teams ( $γ_{00}$ ) plus the effect of the team-average LMX ( $γ_{01}$ ) plus some team-specific deviation ( $u_{0 j}$ ). Similarly, the $β_{1 j}$ equation states that the relation between team-mean-centered LMX and IP in team j is equal to the average relation across all teams ( $γ_{10}$ ) plus some team-specific deviation ( $u_{1 j}$ ).

These average effects (denoted by γ) are referred to as fixed effects because they are constant across teams (i.e., each team-specific effect uses the fixed effect as a basis). The team-specific deviations capture how different the team-specific effect is from the fixed effect. These deviations are called random effects because they are not directly estimated in the model but instead are assumed to follow a multivariate normal distribution with a mean of zero and estimated covariance matrix: $u_{j} \sim N (0, τ)$ . In a standard multilevel model like the one presented in Equation 1, the u terms are location random effects because they dictate the coefficients that relate to the mean of the outcome for each team. The diagonal elements of $τ$ capture the variances of the effects across all teams. The effects are allowed to be different in each team, so this variance quantifies how different the effects are across the teams. The off-diagonal elements of $τ$ capture potential relations between team-specific effects (e.g., a positive value indicates that teams with higher intercepts tend to have larger slopes). After accounting for possible differences in the effects across teams, $r_{i j}$ is the within-team residual for the ith person in the jth team (i.e., $r_{i j} = I P_{i j} - I {\hat{P}}_{j}$ ). Similar to the random effects, $r_{i j}$ is not directly estimated but again assumed to come from a normal distribution with a mean of zero and estimated variance: $r_{i j} \sim N (0, σ^{2})$ .

Graphic Depiction

Because the subscripts and notation can complicate the underlying process of the multilevel model, consider the following graphic depiction for simplified intuition behind the model. Figure 1 shows hypothetical data from four teams with the outcome variable (IP) on the vertical axis and the within-team predictor variable (LMX) on the horizontal axis. The top panel shows a scatter plot of all the data, and the bottom panel compares the team-specific line for Team 1 to the average line for all four teams. The dashed line is the team-specific regression that would be formed from including the random intercept and random slope (the u terms in Equation 1): The intercept is a little higher and the slope a little steeper than the fixed effect regression solid line.

Figure 1.

The top panel shows a scatterplot of four teams of data. The bottom panel shows the team-specific regression line for Team 1 (dashed line) and the average regression line for all four teams (solid line).

The random effects capture how far the team-specific line is from the average line. For this reason, the random effects are sometimes referred to as the “team-level residuals”—the total distance between the data from Team 1 and the average regression line can be accounted for by the fact that a different, team-specific regression line is more appropriate for the Team 1 data. Knowing that these data come from Team 1 reduces the within-team residual (i.e., the team-specific line is better for prediction than the average line, so the data points are much closer to the dashed line than to the solid line). Note that this variation remains unexplained: We do not know why Team 1 has a higher intercept and steeper slope. The random effects are not explaining any variation but rather they are simply partitioning the variance into different sources. The data points are not perfectly on the team-specific line, so the distance between the dashed line and the data is captured by within-team residual ( $r_{i j}$ in Equation 1).

Explaining Team-Level Variation

The random effects help to partition the residual by level, but this variation remains unexplained. However, team-level predictors can be added to the model to explain team-level variance. The elements of $τ$ tell us how much variation there is between teams, and one goal of multilevel models is to explain why these teams are different. Imagine that the data contain a measure of leadership climate (LC) to assess attributes of the team (i.e., this variable is constant across all members within a team but can vary across teams). This variable can be included in the team-level equations (i.e., the equations with β as the outcome in Equation 1) with what Raudenbush and Bryk (2002) referred to as a means-and-slopes as outcomes model:

\begin{array}{l} I P_{i j} = β_{0 j} + β_{1 j} (L M X_{i j} - {\bar{L M X}}_{j}) + r_{i j} \\ β_{0 j} = γ_{00} + γ_{01} {\bar{L M X}}_{j} + γ_{02} (L C_{j} - \bar{L C}) + u_{0 j} \\ β_{1 j} = γ_{10} + γ_{11} (L C_{j} - \bar{L C}) + u_{1 j} . \end{array}

Equation 2 states that the intercept in team j changes by $γ_{01}$ units for every one-unit increase in grand-mean-centered LC, shifting the source of differing team intercepts away from the unexplained random effect ( $u_{0 j}$ ) and assigning it to a particular predictor (LC). A similar explanation applies to the team-specific slope of LMX represented by $β_{1 j}$ . The slope for LMX in team j changes by $γ_{11}$ units for every one-unit increase in LC and similarly shifts the source of team-specific differences from a random effect ( $u_{1 j}$ ) to the predictor, LC.

Example Analysis

Data Description

Similar to Aguinis, Gottfredson, and Culpepper (2013), we used data patterned after the study by Chen et al. (2007) on IP (predicted by LMX within teams and LC at the team level) that was used illustratively in the preceding sections. The data feature 630 people clustered within 105 teams such that each team has 6 members. All variables are treated as continuous, and 0 on the outcome represents average IP across the entire data set. Data and code are provided on the author’s Open Science Framework page, https://osf.io/8peqb.

To these data, we fit a location-only Bayesian multilevel model in Mplus 8.3. Frequentist estimation methods can have difficulty accommodating the location-scale model (but not the standard model; Hedeker et al., 2009; Rast et al., 2012), so we fit the location-only model in a Bayesian framework to compare the results from this section to the next section. By default, Mplus uses Bayesian methods for computational convenience via noninformative priors rather than for philosophical reasons (Asparouhov et al., 2018), so the results are interpreted very similarly to a frequentist analysis. Basics of Bayesian estimation are discussed in the Appendix for readers who may not be well oriented to these methods. The estimation was conducted with two chains, a minimum of 1,000 iterations, and the iterations were stopped when the potential scale reduction (PSR; Gelman & Rubin, 1992) fell below 1.05 (PSR is similar to convergence criteria in maximum likelihood to determine when the estimates are considered to be “stable enough”).

Mplus Code

Annotated Mplus code for this model is as follows:

VARIABLE:

NAMES ARE lc ip lmx team;

USEVAR ARE lc ip lmx lmx_b;

! The Cluster ID is the Team variable;

CLUSTER = team;

! LMX varies within Teams;

WITHIN = lmx;

! LC only varies between Teams;

! LMX_B is the team mean of LMX;

BETWEEN = lc lmx_b;

DEFINE:

! team-mean center LMX;

CENTER lmx (GROUPMEAN);

! Grand-mean center LC;

CENTER lc (GRANDMEAN);

! Create new variable for tean mean of LMX;

LMX_B=Cluster_Mean(LMX);

! Grand-mean center team mean of LMX;

CENTER LMX_B (GRANDMEAN);

ANALYSIS:

! Two-Level multilevel model with random slopes;

TYPE = TWOLEVEL RANDOM;

! Use Bayesian MCMC estimation;

ESTIMATOR = BAYES;

! Run at least 1000 iteration of the MCMC algorithm because allowing convergence;

BITERATIONS =(1000);

MODEL:

! The within-team equation;

%WITHIN%

! Ip is predicted by lmx, this slope is random and varies across teams;

beta1| ip ON lmx;

! Estimate the within-team residual variance for ip;

ip;

! team-level equations;

%BETWEEN%

! Estimate the fixed effect for the effect of ip ON lmx;

[beta1];

! Estimate the fixed effect for the intercept of ip;

[ip];

! Estimate the random effect variance of the ip on lmx slope;

beta1 (a);

! estimate the random effect variance for the intercept of ip;

Ip (b);

! The slope of ip ON lmx is predicted by lc;

beta1 ON lc;

! the intercept of ip is predicted by lc and the team mean of lmx;

ip ON lc lmx_b;

! covariance between random effects;

ip with beta1 (c);

MODEL CONSTRAINTS:

! track new parameter called corr;

new(corr);

! the correlation is the covariance divided by square roots of the variances;

corr = c / sqrt(a*b);

! Provides standardized coefficients and R-square estimates;

OUTPUT: STANDARDIZED;

With multilevel SEM in Mplus, a vertical pipe is included to signify that paths are team-specific, meaning that they have random effects and the path will appear team-level model (i.e., under the %BETWEEN% heading). The label for the random effect appears before the vertical pipe, and researchers can assign any name they wish to this label because it does not correspond to a variable in the data set. Here, we use Greek letters from Equation 2. For instance, “beta1| ip ON lmx” means that (team-mean-centered) LMX predicts IP and that this relationship varies across teams. Then, in the team-level equation, the fixed effect is specified with [beta1], which is akin to $γ_{10}$ in Equation 2. Team-level predictors are included via regression by placing the random slope label as the dependent variable: beta1 on lc corresponds to $γ_{11}$ in Equation 2.³

These guidelines generalize to all predictors, but intercepts function slightly differently in Mplus. The outcome variable by itself in the within-team equation under the %WITHIN% heading (i.e., “ip;”) corresponds to the within-team residual variance, $σ^{2}$ . There is no “beta0|” label under the %WITHIN% heading because Mplus uses the outcome variable to specify the intercept fixed effect (when in brackets under the %BETWEEN% heading), the team-level intercept variance (the outcome variable label by itself under the %BETWEEN% heading), and the residual variance (the outcome variable label by itself under the %WITHIN% heading).

Table 1 shows the estimates from the location-only multilevel model. Despite the Bayesian estimation, the interpretation is essentially the same as a frequentist multilevel model (with the exception that posterior summaries and credible intervals [CI] are substituted for point estimates and confidence intervals).

Table 1.

Posterior Summaries and 95% Credible Intervals From the Standard Multilevel Model.

Effect	Notation	Posterior	95%
Effect	Notation	Median	Credible Interval
Fixed Effects
Intercept	$γ_{00}$	−0.14	[−0.37, 0.11]
LC (GMC)	$γ_{01}$	0.23	[−0.05, 0.49]
LMX (TM)	$γ_{02}$	1.26	[0.74, 1.74]
LMX (TMC)	$γ_{10}$	0.84	[0.67, 1.01]
LMX × LC	$γ_{11}$	0.32	[0.12, 0.53]
Variance Components
Intercept	$τ_{00}$	0.98	[0.66, 1.46]
LMX	$τ_{11}$	0.38	[0.19, 0.68]
Corr(Intercept, LMX)	$C o r r (u_{0 j}, u_{1 j})$	0.06	[−0.31, 0.46]
Within-team residual	$σ^{2}$	2.34	[2.05, 2.66]

Note: LMX = leader-member exchange; LC = leadership climate; GMC = grand-mean centered; TM = team mean; TMC = team-mean centered.

It does not appear that LC affects IP because 0 appears in the credible interval ( $γ_{01} = 0.23, 95 % CI = [- 0.05, 0.49]$ ). LC does appear to affect the slope of LMX ( $γ_{11} = 0.32, 95 % CI = [0.12, 0.53]$ ). Figure 2 shows a region of significance plot created from $γ_{10}$ , $γ_{01}$ , and $γ_{11}$ as recommended by Hoffman (2015) to assess where the effect of LMX on IP occurs. Using the procedure described by Preacher, Curran, and Bauer (2006), the simple slope is positive and nonnull for LC values above –1.502.⁴ Conditional on predictors in the model, teams appear to vary in their IP intercept ( $τ_{00} = 1.31, 95 % CI = [0.92, 1.84]$ ), and the slope of LMXvaries across teams ( $τ_{11} = 0.43, 95 % CI = [0.23, 0.77]$ ). The within-team residual variance is constrained across all teams and is estimated to be 2.33 ( $95 % CI = [2.05, 2.66]$ ).

Figure 2.

Region of significance plot for cross-level interaction of LMX and LC. The black line is the simple slope at each value of LC; the dashed lines are the confidence limits for the simple slope. The simple slope of LMX is null in the gray shaded area (LC values below –1.502) and nonnull elsewhere. LMC = leader-member exchange; LC = leadership climate.

Weakness of the Standard Location-Only Multilevel Model

In the previous sections, we noted how the β coefficients in the within-team equation (the first equation with IP as the outcome) become outcomes in the team-level equations where they can take random effects or team-level predictors to account for differences across teams. This follows from the colloquial rule that any term on the right side of the within-team equation can drop down to become an outcome in a team-level equation. The location-only multilevel model follows this adage with respect to the coefficients associated with location; however, it does not allow the scale to vary by team (i.e., the within-team residual variance is constant across teams).

Specifically, the within-team residual variance $σ^{2}$ is also a within-team quantity, but it is treated differently compared to the other within-team quantities (Hoffman, 2007). That is, although the within-team coefficients become outcomes at the team level, the within-team residual variance is constant across all teams even though it is conceivable that this residual variance differs across teams. Figure 3 shows the same hypothetical data from Figure 1 but now shows two team-specific lines. Data from the black team show relatively little variation around the team-specific regression line, whereas data from the gray team show much larger variation, suggesting that the variances in these teams are not equivalent.

Figure 3.

Comparison of variation around the team-specific regression lines for two separate hypothetical teams. The black team adheres much more closely to the team-specific line, whereas the gray team has much greater dispersion. In the location-only multilevel model, the variance around the team-specific lines is constrained and would miss the opportunity to model the differing variability in these teams.

In the location-only multilevel model, the within-team residual variance is captured by a single parameter $σ^{2}$ that does not have a team-specific j subscript. In addition to applying a possibly unwarranted and frequently untested assumption to the model, there is a missed opportunity for researchers to model systematic reasons why the gray team is more variable than the black team, which can be as illuminating as modeling why the team-specific intercepts and slopes differ. For example, does a characteristic of the team leader lead to more heterogeneity in the gray team? Do different demographic compositions of teams lead to different levels of empowerment at the same values of LMX?

Commonly, these types of questions are addressed by fitting a second model with a summary metric of variability as the outcome, such as the standard deviation of the dependent variable in each team or the coefficient of variation (Lester et al., 2020). Although this approach captures the spirit of the intended analysis, the two-step approach has limitations. First, using a descriptive summary measure such as the raw unconditional standard deviation of the dependent variables does not allow researchers to account for variance that may be explained by other predictors in the model. Second, if an approach that accounts for other predictors is used (e.g., saving team-specific values from a multilevel model), a two-step approach does not acknowledge uncertainty in the variability estimates for each team, which will inflate Type I error rates in the subsequent regression, especially for teams with fewer members (Lester et al., 2020). Third, a two-step approach omits the ability to correlate the random effects across submodels (e.g., scale submodel random effects associated with variability and location submodel random effects associated with the mean). Fourth, a two-step approach results in zero within-team variability in the outcome, meaning that within-team predictors would be excluded from the analysis (Leckie et al., 2014).

Heterogeneous variance models have been proposed to encompass the analysis of both the location and scale submodels within a single model (e.g., Hedeker & Mermelstein, 2007; Lang et al., 2018). Heterogeneous variance models deterministically model the residual variance as a function of predictors with fixed effects to address these types of questions. If there is unexplained variance in the scale submodel and a random effect is not included, simulation studies have shown that Type I error rates will be inflated (Leckie et al., 2014; Walters et al., 2018). However, the model can be extended to include random effects in both the location and the scale submodels with what is deemed a location-scale model (Hedeker et al., 2008). The next section and the remainder of this article focus on this model.

The Location-Scale Multilevel Model

The location-scale multilevel model looks very similar to the location-only multilevel model presented in Equation 1 but with one small change. Specifically, the location-scale version of the model in Equation 1 is written as:

\begin{array}{l} I P_{i j} = β_{0 j} + β_{1 j} (L M X_{i j} - \bar{L M X_{j}}) + r_{i j} \\ β_{0 j} = γ_{00} + γ_{01} \bar{L M X_{j}} + u_{0 j} \\ β_{1 j} = γ_{10} + u_{1 j}, \end{array}

where⁵

r_{i j} \sim N (0, σ_{j}^{2}) .

The subtle difference is easy to miss if reading quickly: The within-team residual variance $σ^{2}$ in Equation 4 now has a j subscript, meaning that it is allowed to vary across teams. Recall that any portion of the within-team equation that has a j subscript becomes an outcome in a team-level equation, so the team-specific residual variance now becomes an outcome. This means the location-scale multilevel model includes another equation, called the scale submodel:

σ_{j}^{2} = exp (ω_{0} + u_{2 j}),

where $ω_{0}$ is the average within-team residual across all teams and $u_{2 j}$ is a scale random effect that captures the deviation of the jth team’s within-team residual from the average, which is assumed to follow a normal distribution with a mean of zero with an estimated variance. This extends the u vector to three elements (random location intercept, random location slope, random scale intercept) such that

[\begin{matrix} u_{0 j} \\ u_{1 j} \\ u_{2 j} \end{matrix}] \sim N ([\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}], [\begin{matrix} τ_{00} \\ τ_{10} & τ_{11} \\ τ_{20} & τ_{21} & τ_{22} \end{matrix}]) .

Note that the location and scale random effects are modeled within a single matrix, which allows for direct modeling of covariances between random effects in different submodels.

This generalization allows the within-team residual to differ, meaning that it would detect that the black and gray teams from Figure 3 vary in different amounts around their team-specific regression lines. If variables are thought to be related to the scale, Equation 5 can be appended to include predictor variables to explain why some teams are more variable than others. If LC was thought to also explain why the within-team residual variances differ across teams, Equation 5 could be expanded to

σ_{j}^{2} = exp [ω_{0} + ω_{1} (L C_{j} - \bar{L C}) + u_{2 j}],

where $ω_{1}$ would represent how the within-team residual variance changes as a function of LC (e.g., perhaps the black team in Figure 3 has low within-team residual variance because their LC is strong, so IP is less erratic conditional on LMX). Unlike heterogeneous variance models, the between-team variance in the scale submodel need not be a deterministic function of predictors, and between-team differences are allowed to be unexplained via the presence of $u_{2 j}$ .

An important point in Equations 5 and 7 is that the within-team residual variance is modeled as the exponential of the linear model on the right (i.e., it is a log-linear model). Alternatively, one could apply the natural logarithm to each side such that the natural log of the within-team residual variance follows a linear model. This is necessary because variances cannot be negative, and using a log-linear model precludes negative values. This will have important consequences for interpretation of results, which we cover in the next section.

Example Analysis

Using the same data as the previous example, we fit the same location-scale model defined by Equations 3, 4, 6, and 7 in Mplus Version 8.3 using the same estimation options as the location-only model. The Mplus code for this model is as follows (with differences from the previous code bolded for emphasis):

VARIABLE:

NAMES ARE lc ip lmx team;

USEVAR ARE lc ip lmx lmx_b;

! The Cluster ID is the Team variable;

CLUSTER = team;

! LMX varies within Teams;

WITHIN = lmx;

! LC only varies between Teams;

! LMX_B is the team mean of LMX;

BETWEEN = lc lmx_b;

DEFINE:

! team-mean center LMX;

CENTER lmx (GROUPMEAN);

! Grand-mean center LC;

CENTER lc (GRANDMEAN);

! Create new variable for tean mean of LMX;

LMX_B=Cluster_Mean(LMX);

! Grand-mean center team mean of LMX;

CENTER LMX_B (GRANDMEAN);

ANALYSIS:

! Two-Level multilevel model with random slopes;

TYPE = TWOLEVEL RANDOM;

! Using Bayesian MCMC estimation;

ESTIMATOR = BAYES;

! Run at least 1000 iteration of the MCMC algorithm because allowing convergence;

BITERATIONS =(1000);

MODEL:

! The within-team equation;

%WITHIN%

! Ip is predicted by lmx, this slope is random;

beta1| ip ON lmx;

! place a random effect on the within-team residual variance;

Logv|ip;

! team-level equations;

%BETWEEN%

! Estimate the fixed effect for the effect of ip ON lmx;

[beta1];

! Estimate the fixed effect for the intercept of ip;

[ip];

! Estimate the average residual variance ON THE LOG SCALE;

[logv];

! Estimate the random effect variance of the ip on lmx effect;

beta1 (a);

! estimate the random effect variance for the intercept of ip;

Ip (b);

! Estimate the random effect variance for the residual variance ON THE LOG SCALE;

Logv (c);

! The slope of ip ON lmx is predicted by lc;

beta1 ON lc;

! the intercept of ip is predicted by lc and the team mean of lmx;

ip ON lc lmx_b;

! the residual variance is predicted by lc with a LOG-LINEAR MODEL;

logv on lc;

! covariance between intercept and slope random effects;

ip with beta1 (c1);

! covariance between intercept and residual variance random effects;

ip with logv (c2);

! covariance between slope and residual variance random effects;

beta1 with logv (c3);

MODEL CONSTRAINTS:

! new parameter for intercept and slope random effect correlation;

new(corr_I_lmx);

! new parameter for intercept and residual variance random effect correlation;

new(corr_I_logv);

! new parameter for slope and residual variance random effect correlation;

new(corr_lmx_logv);

! the correlation formula for intercept and slope random location effect correlation;

corr_I_lmx = c1 / sqrt(a*b);

! the correlation formula for intercept location and intercept scale random effect correlation;

corr_I_logv = c2 / sqrt(b*c);

! the correlation formula for slope location and intercept scale random effect correlation;

corr_lmx_logv = c3/ sqrt(a*c);

! Provides standardized coefficients and R-square estimates;

OUTPUT: STANDARDIZED;

The main difference between the code for the location-only multilevel model and the location-scale multilevel model takes place in the within-team equation. A random effect is placed on the within-team residual variance, logv|ip. Again, the label does not have to be “logv,” and researchers can choose whichever label they wish. Mplus will automatically model within-team residual variances that vary across teams with a log-linear model, so we use “logv” as a reminder that the estimates using this variable as an outcome are on a different scale than the other estimates. Otherwise, the remainder of the code follows the same guidelines provided in the previous section. The intercept of Equation 7 ( $ω_{0}$ ) is represented by [logv] in a team-level equation, the effect of LC in Equation 7 ( $ω_{1}$ ) is represented by logv on lc, and including logv by itself in a team-level equation estimates the within-team residual variance random effect variance across teams ( $τ_{22}$ in Equation 6).

Table 2 shows the results from the location-scale multilevel model. The location portion of the model features similar estimates to Table 1 with similar conclusions, so we will not plot the interaction effect. The notable difference in Table 2 occurs at the bottom, under the “Scale Fixed Effects” heading. The estimate of intercept for Equation 7 (the expected residual variance for a team at the grand-mean of LC) is .76, on the log scale. To transform this back onto the scale of the outcome as in Table 1, we need to exponeniate: $exp (0.76) = 2.14$ , meaning that the raw residual variance is 2.14, on average.

Table 2.

Posterior Summaries and 95% Credible Intervals From the Location-Scale Multilevel Model.

Effect	Notation	Posterior	95%
Effect	Notation	Median	Credible Interval
Location Fixed Effects
Intercept	$γ_{00}$	−0.15	[−0.37, 0.10]
LC (GMC)	$γ_{01}$	0.26	[−0.02, 0.53]
LMX (TM)	$γ_{02}$	1.32	[0.80, 1.81]
LMX (TMC)	$γ_{10}$	0.85	[0.66, 1.03]
LMX × LC	$γ_{11}$	0.31	[0.10, 0.52]
Location Random Effect Covariance Structure
Intercept	$τ_{00}$	1.01	[0.69, 1.49]
LMX	$τ_{11}$	0.39	[0.17, 0.73]
Corr(Intercept, LMX)	$C o r r (u_{0 j}, u_{1 j})$	−0.04	[−0.42, 0.28]
Scale Fixed Effects
Ln (Intercept)	$ω_{0}$	0.76	[0.57, 0.93]
Ln [LC (GMC)]	$ω_{1}$	−0.21	[−0.42, −0.05]
Scale Random Effect Covariance Structure
Team-level scale variance	$τ_{22}$	0.23	[0.09, 0.49]
Corr(Intercept, Scale Variance)	$C o r r (u_{0 j}, u_{2 j})$	0.14	[−0.28, 0.51]
Corr(LMX, Scale Variance)	$C o r r (u_{1 j}, u_{2 j})$	−0.38	[−0.76, 0.16]

Note: LMX = leader-member exchange; LC = leadership climate; GMC = grand-mean centered; TM = team mean; TMC = team-mean centered.

The estimated effect of LC is −0.21, which again if we exponentiate is $exp (- 0.21) = 0.81$ . Broadly, the negative coefficient on the log scale means that the within-team residual variance decreases when LC increases. In more detail, effects in log-linear models are multiplicative, not additive. This means that the interpretation is: For a one-unit change above the grand-mean of LC, the residual variance is expected be 0.81 times its intercept value ( $2.14 \times 0.81 = 1.73$ ). Rather than multiplying as in linear models, the exponent of the effect is changed to note a change in a predictor. So for a 1.5-unit increase above the grand-mean of LC, the effect would be ${0.81}^{1.5} = 0.73$ times the intercept value ( $2.16 \times 0.73 = 1.58$ ); for a one-unit decrease, the effect would be ${0.81}^{- 1} = 1.25$ times the intercept value ( $2.16 \times 1.25 = 2.70$ ). The scale intercept random effect variance ( $τ_{22} = 0.19, 95 % CI = [0.04, 0.40]$ ) suggests that there is variation in the within-team residual variance across teams that is not explained by LC.

Despite the unorthodox interpretation, the substantive interpretation of the model is useful beyond that provided by the results in Table 1. Specifically, the negative effect of LC in the scale submodel indicates that as LC increases, there is less variability on IP within-teams. To put this into context, Figure 4 displays the team-specific regression lines for two teams in the data. One team is high on LC (Team 34 in black), and one team is low on LC (Team 69 in gray). Based on the location-scale multilevel model, the members from Team 34 (with high LC) are much closer to one another after conditioning on LMX (i.e., the black dots very closely follow the dashed black line). Contrarily, the members of Team 69 are much less consistent, and they vary much more around the team-specific regression line (i.e., there is more within-team heterogeneity).

Figure 4.

Comparison of a homogeneous team (Team 34 in black) and a heterogeneous team (Team 69 in gray). The leadership climate in Team 34 was 1.50 SD above the overall mean, whereas the Team 69 leadership climate was 1.57 SD below the mean.

Example With Three-Level Repeated Measures Data

To demonstrate how specifying the model in the multilevel SEM framework can help address difficulties that currently exist in organizational research, consider the case study where repeated measures are nested within people who are further nested within groups. One such example appears in Lang et al. (2018), who described a three-level model with this structure for consensus emergence. One of their examples features data measuring job satisfaction three times for 471 US Army soldiers nested within 34 units. The data are publicly available in the multilevel R package and are posted on the OSF page for this article. Lang et al. (2018, 2019) noted that software limitations can hinder inclusion of scale random effects and that even then, computation is usually feasible only with large sample sizes. This is certainly true in mixed effect framework and with frequentist estimation, but this issue can be circumvented in a Bayesian multilevel SEM framework. In this section, we first show how to reproduce the three-level model from Lang et al. (2018) as a Bayesian multilevel SEM. Then, we show how multiple scale random effects can be included and correlated with location random effects in a three-level model, even with a Level 3 sample size of only 34 units.

First, scale random effects in multilevel SEM are not yet available for pure three-level models in Mplus, but they are available for two-level models. Asparouhov and Muthén (2016) noted that hierarchies can be reduced by one if repeated measures at the lowest level are transposed from the long format (one row per repeated measure) to the wide format (one column per repeated measure). Doing so yields a two-level hierarchy of soldiers within units, which is more amenable to the multilevel SEM framework.

Next, Lang et al. (2018) accurately noted that Mplus does not permit predictors in a MODEL CONSTRAINT statement with two-level models, which seemingly prevents scale submodel predictors and random effects. However, the MODEL CONSTRAINT statement can be bypassed via phantom variables. Phantom variables originated decades ago as a way to bypass software limitations (Rindskopf,1984). Although less common in the modern software environment, phantom variables retain their utility in contexts where statistical developments outpace software developments. Essentially, a phantom variable is a latent variable with a single indicator whose loading is constrained to 1 and whose residual variance is constrained to 0. A phantom variable retains the exact same properties as its indicator and is used simply to transform a parameter that is not permitted in modeling statements into a variable that is permitted to appear in modeling statements. In this context, phantom variables make the within-person residual variance directly available as an outcome rather than requiring researchers to use a MODEL CONSTRAINT statement.

We will first reproduce Model 3c from Table 9 in Lang et al. (2018). The model can be written

\begin{array}{l} J o b S a t i s f a c t i o n_{t i j} = π_{0 ij} + π_{1 ij} T i m e_{t} + e_{t i j} \\ π_{0 i j} = β_{00 j} + r_{0 i j} \\ π_{1 i j} = β_{10 j} \\ β_{00 j} = γ_{000} + γ_{010} R e a d i n e s s_{j} + u_{00 j} \\ β_{10 j} = γ_{100} + γ_{110} R e a d i n e s s_{j} + u_{10 j} \\ e_{t i j} \sim N (0, σ_{j}^{2}) \\ σ_{j}^{2} = exp [ω_{0} + ω_{1} T i m e_{t} + ω_{2} R e a d i n e s s_{j} + ω_{3} (T i m e_{t} \times R e a d i n e s s_{j})] \\ r_{0 i j} \sim N (0, τ_{00}) \\ [\begin{matrix} u_{00 j} \\ u_{10 j} \end{matrix}] \sim M V N ([\begin{matrix} 0 \\ 0 \end{matrix}], [\begin{matrix} ν_{00} \\ ν_{01} & ν_{11} \end{matrix}]), \end{array}

where readiness is a z-scored predictor of combat readiness of each unit. The data are nested within three levels, so a t subscript is added to index repeated measures nested within the ith person in the jth unit. The Mplus code for the model is:

VARIABLE:

Names Are Unit ready js0-js2;

USEVARIABLES = js0-js2 readyz;

Missing are all (9999);

! Readiness Z-score is has not variability within units;

Between=readyz;

! Unit is the Unit-level ID variable;

CLUSTER=Unit;

! Create Z-score for the readiness;

DEFINE: readyz=(ready-3)/.309303;

ANALYSIS:

! Two-Level multilevel model with random slopes;

TYPE = TWOLEVEL RANDOM;

! Use Bayesian MCMC estimation;

ESTIMATOR = BAYES;

! Run at least 10000 iteration of the MCMC algorithm before allowing convergence;

BITERATIONS =(10000);

MODEL:

%WITHIN%

! Person-level component of growth model;

Pi0 by js0@1 js1@1 js2@1;

! Phantom variable for Time 1 residual variance;

logv1 | js0;

! Phantom variable for Time 2 residual variance;

logv2 | js1;

! Phantom variable for Time 3 residual variance;

logv3 | js2;

! Person-level variance for the intercept;

Pi0;

! Person-level intercept mean is 0;

[Pi0@0];

%BETWEEN%

! Unit-level component of growth model;

Beta0 Beta1 | js0@0 js1@1 js2@2;

! Fixed effects for intercept and slope;

[Beta0 Beta1] (gamma000 gamma100);

! Unit-level variance of growth factors;

Beta0 (LI);

Beta1 (LS);

Beta0 WITH Beta1 (LCov);

! Growth model on the phantom variables to create scale submodel;

Omega0 Omega1 | logv1@0 logv2@1 logv3@2;

! Scale Model Intercept;

[Omega0](d0);

! Scale Model Time Effect;

[Omega1](d1);

! Convariance of Scale submodel growth factors. These are 0 in Lang et al. but can be estimated;

Omega0@0 (SI);

Omega1@0 (SS);

Omega0 with Omega1@0 (SCov);

! Location, Scale Covariance. These are 0 in Lang et al. but can be estimated;

Omega0 with Beta0@0 (LiSiCov);

Omega0 with Beta1@0 (LSSiCov);

Omega1 with Beta0@0 (LiSsCov);

Omega1 with Beta1@0 (LsSsCov);

! variance of phantom variables is 0;

logv1@0;

logv2@0;

logv3@0;

! Unit-level variance of outcome is already accounted for phantom variables;

js0@0 js1@0 js2@0;

! Growth factors and within-person variance predicted by Z-score of readiness;

Beta0 Beta1 on readyz (Gamma010 Gamma110);

Omega0 Omega1 on readyz (Omega2 Omega3);

MODEL CONSTRAINT:

! Put estimates onto SD scale to match Lang et al.;

new(delta1); new(delta2);new(delta3); new(sigma2);

delta1=d1/2;

delta2=Omega2/2;

delta3=Omega3/2;

sigma2=exp(d0);

In the %WITHIN% portion of the code, three phantom latent variables are created for the within-person residual variances (labeled logv1, logv2, and logv3). This is done to trick Mplus into treating the residual variances as variables rather than parameters. In the %BETWEEN% portion of the code, we then introduce VI and VS as growth factors for the phantom variables. [VI] is the log of the residual variance at Time = 0 ( $ω_{0}$ in Equation 8), and [VS] is the effect of time on the log of the residual variance ( $ω_{1}$ in Equation 8). Including scale submodel predictors is then accomplished by adding them as predictors of the phantom variable growth model intercept (for main effects) or the phantom variable growth model slope (for interactions with time). Lang et al. (2018) put their scale submodel predictor effects on the scale of the standard deviation rather than the variance, so the MODEL CONSTRAINT statement divides the scale submodel fixed effects by 2 to make the estimates comparable to their results. Table 3 shows the estimates between the multilevel SEM specification and the estimates reported in Table 9 of Lang et al. (2018). Results show only slight differences due to different estimation methods.

Table 3.

Comparison of Estimates From the Model 3c in Lang, Bliese, and Runge (2019), the Reproduction of That Model as a Bayesian Multilevel SEM, and the Extension of the Model That Allows for Scale Random Effects Estimated as a Multilevel SEM.

	Lang et al.Model 3c	ML-SEMModel 3c	UnconstrainedML-SEM
Location Fixed Effects
Intercept, $γ_{000}$	3.22	3.23	3.25
Time, $γ_{100}$	0.06	0.05	0.05
Readiness, $γ_{010}$	0.26	0.26	0.25
Time × Readiness, $γ_{110}$	–0.06	–0.06	–0.06
Location Random Effect Covariances^a
Unit Intercept Variance, $ν_{00}$	0.00	0.02	0.07
Unit Time Variance, $ν_{11}$	0.00	0.01	0.01
Unit Time, Intercept Covariance, $ν_{10}$	0.00	–0.01	–0.02
Person Intercept Variance, $τ_{00}$	0.37	0.38	0.39
Scale Fixed Effects
Intercept, $exp (ω_{0})$	0.51	0.52	0.42
Time, $ω_{1}$	–0.09	–0.11	–0.11
Readiness, $ω_{2}$	–0.15	–0.13	–0.14
Time × Readiness, $ω_{3}$	0.09	0.08	0.07
Scale Random Effect Variances
Intercept Variance, $ν_{22}$	—	—	0.36
Time Variance, $ν_{33}$	—	—	0.11
Intercept, Time Covariance, $ν_{32}$	—	—	–0.06
Location-Scale Random Effect Covariances
Location Intercept, Scale Intercept, $ν_{20}$	—	—	–0.08
Location Intercept, Scale Time, $ν_{30}$	—	—	–0.01
Location Time, Scale Intercept, $ν_{21}$	—	—	0.02
Location Time, Scale Time, $ν_{31}$	—	—	0.01

Note: Model 3c = the model originally presented in Table 9 of Lang Bliese, and de Voogt (2018). Estimates from models estimated with Bayesian methods are the median of the posterior distribution for each parameter. ML-SEM = multilevel structural equation model.

^a Covariances were reported in Lang et al. (2018), so we provide covariances estimates for ease of comparison rather than correlation as we present in Tables 1 and 2.

Whereas scale submodel random effects and random effect correlation are cumbersome and complex to obtain in a mixed effects framework, they can be easily obtained with multilevel SEM by simply unconstraining certain parameters. We can add scale submodel random intercepts and random slopes by unconstraining the variance of VI and VS.⁶ The scale submodel random effects can be correlated with each other and also with the location submodel random effects. The full code is included on the Open Science Framework Page, and results are shown in the last column of Table 3. This would augment Equation 8 such that

\begin{array}{l} σ_{j}^{2} = exp [(ω_{0} + u_{20 j}) + (ω_{1} + u_{30 j}) T i m e_{t} + ω_{2} R e a d i n e s s_{j} + ω_{3} (T i m e_{t} \times R e a d i n e s s_{j})] \\ [\begin{matrix} u_{00 j} \\ u_{10 j} \\ u_{20 j} \\ u_{30 j} \end{matrix}] \sim M V N ([\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \end{matrix}], [\begin{matrix} ν_{00} \\ ν_{10} & ν_{11} \\ ν_{20} & ν_{21} & ν_{22} \\ ν_{30} & ν_{31} & ν_{32} & ν_{33} \end{matrix}]) . \end{array}

Whereas Lang et al. (2019) noted the hurdles for fitting this type of model in the mixed effect framework with frequentist estimation, allowing for all of these aspects presents no issues in the Bayesian multilevel SEM framework, even with only 34 units. In fact, the model converged in only 31 seconds, vastly eclipsing the demanding nonlinear estimation that would be required if trying to fit the model in another suitable software environment such as SAS Proc Nlmixed. By allowing these effects in in the model, Table 3 shows that the random effect variance of the scale intercept and scale effect of time are fairly large ( $ν_{22} = .38, ν_{33} = .11$ ). Substantively, this indicates that there are unexplained reasons beyond readiness that account for changes in variability over time. Without these random effects, it would be assumed that readiness perfectly accounts for the reasons why variability changes over time.

Implications for Organizational Research

By its nature, organizational research has a vested interest not just in averages of organizations but also in their cohesiveness. In addition to the two examples provided in this article, other examples for organizational research where cross-sectional questions about variance is prominent include consensus in ratings of job competencies (Lievens et al., 2010), heterogeneity in occupational knowledge and skills for people with similar jobs (Sitzmann, Ployhart, & Kim, 2019), perceptions of team efficacy (DeRue et al., 2010), and variation in policy implementation (Pak & Kim, 2018). When considering variability repeated measures data, possibilities include variability in decision making (Melkonyan & Safra, 2016), variability in personality across situations (Dalal et al., 2015), variability in job performance over time (Barnes, Reb, & Ang, 2012; Reb & Greguras, 2010), or fluctuations in personality dynamics and its effect on personnel selection (Sosnowska, Hofmans, & Lievens, in press).

Despite the frequency of questions about cohesion and heterogeneity in the organizational literature, Lang et al. (2018, 2019) noted how there are few methodological tools that have been extended to study aspects of cohesiveness or similarity within organizational groups. Heterogeneous variance models are undoubtedly an improvement and serve as an excellent starting point to include a scale submodel. However, these models can be further expanded with the inclusion of random scale effects to quantify unexplained variance in the scale submodel and protect against inflated Type I error rates for fixed effects that can arise from misspecifying the random effect structure.

Adopting location-scale models in a multilevel SEM framework would seem to have particular advantages for addressing the possible limitations noted with applications of heterogeneous variance models in the recent literature. For instance, Lang et al. (2019) noted random scale effects as a future direction with the qualification that such models often cannot include random location slopes on account of their computational encumbrances (e.g., the MIXREGLS program does not permit random location slopes). As we show in our examples, location-scale models in the multilevel SEM framework can accommodate multiple scale submodel random effects and can correlate them with location scale effects without issue. Furthermore, the computational burden was essentially nonexistent given that the models converged in less than 1 minute.

Lang et al. (2019) also discussed the possibility of merging heterogeneous variance models with other common organizational techniques such as mediation models but noted “it is not clear how the approach [mediation] could be extended to emergence models” (p. 15). To contextualize this in the context of organizational research, perhaps the effect of leadership climate on task efficiency is mediated by the within-team variability of individual empowerment such that teams with a strong leadership climate are more homogeneous on individual empowerment, which leads to improved efficiency because team members work more cohesively. Whereas mediation would be a difficult leap for a heterogeneous variance model in a mixed effect framework and would likely require a multistage model whose estimates would require Croon’s correction (Croon, 2002) to be unbiased (Devlieger & Rosseel, 2017; Kelcey, Cox, & Dong, 2019), mediation is another straightforward extension if modeling the random scale effects in the multivariate multilevel SEM framework (Hoffman, 2019). McNeish and Hamaker (2020) noted such advantages of multilevel SEM over mixed effects models, noting that multilevel SEM allows any hypothesized structural relations to be modeled between random effects. That is, if the scale submodel has a random effect, it can serve as a mediator or predictor of other variables of interest. Schultzberg and Muthén (2018) included an example of a mediation model with a scale random effect in their study (Model 9, p. 427).

Limitations

Bayesian methods are noted for their better performance with small samples (Muthén & Asparouhov, 2012); however, the computational Bayes approach taken by Mplus can be biased with smaller samples (van de Schoot et al., 2015). The issue stems from the use of diffuse, improper priors under the default settings, which are known to artificially inflate the variance estimates when the number of teams is about 50 or less (Hox, van de Schoot, & Matthijsse, 2012). This problem can be avoided by overriding the defaults with informative priors in Mplus (McNeish, 2016), but effectively doing so usually requires some knowledge of Bayesian statistics. If details on Bayesian statistics are outside of a researcher’s capabilities, a model with fewer than 50 groups run with Mplus’s defaults will reduce the ability to determine that predictors have nonnull effects; the fixed effects coefficients will still be accurately estimated, however.

It is also important to note that the within-group residual variance is based on the within-group sample size. As the team size or number of repeated measures decreases, estimates of the residual variance become less reliable (Estabrook, Grimm, & Bowles, 2012). Although location-scale models can be fit effectively when the within-unit sample size is quite low (Williams, Zimprich, & Rast, 2019), ideally the data would have 20 or more within-group observations (Rast & Zimprich, 2010). Estabrook et al. (2012) noted that small within-group sample sizes do not preclude location-scale modeling; power (or its Bayesian analog) of predictors in the scale submodel will be reduced, however. For two-level models, if groups have varying sizes and there are some very small groups, these small groups need not be deleted because the model will be able to borrow information from the larger groups. Three-level models with multivariate repeated measures have difficulty with varying times of observation or a different number of repeated measures per person (e.g., McNeish & Matta, 2018)

Concluding Remarks

The organizational literature has recently acknowledged the benefit of statistical models with heterogeneous variances across groups. Most proposed models have fixed effects in the scale submodel but have tended to avoid random scale effects due to (a) lack of implementation in standard mixed effect software, (b) existing programs requiring users to be rather skilled at statistical programming, or (c) limitations on the number of random effects in existing programs. Specifying models as Bayesian multilevel SEMs in the latest release of Mplus can alleviate many of these issues in that (a) Mplus is familiar to a large segment of organizational researchers, (b) the computation is expedient, (c) there are less stringent limitations with respect to the number of random effects that can be accommodated, and (d) the model can be embedded into larger models like mediation.

Appendix: Overview of Bayesian Estimation in Mplus

We will not fully review Bayesian statistics because they are used in Mplus more so as a computational tool rather than for its philosophical principles (Asparouhov, Hamaker, & Muthén, 2018). However, we will cover some relevant differences in terminology for readers who are not well versed in Bayesian statistics to facilitate proper interpretation the output. Readers interested in a deeper understanding of Bayesian statistics are referred to Kruschke, Aguinis, and Joo, (2012); van de Schoot et al. (2014); or Zyphur, Oswald, and Rupp (2015) for pedagogical introductions.

A major difference between maximum likelihood and Bayesian Markov chain Monte Carlo (MCMC) is that maximum likelihood yields a single point estimate for each parameter, whereas MCMC yields an entire distribution of possible values for each parameter. This distribution is referred to as the posterior distribution. This posterior distribution is formed by updating researchers’ prior beliefs for each parameter (prior distributions) with information from the data (the likelihood, the same likelihood used in maximum likelihood).

If researchers do not have specific prior beliefs, Mplus will assign parameters uninformative priors such that all admissible values are equally plausible. Retaining the default prior distributions in Mplus yields results that are very close to maximum likelihood, especially at larger samples (Muthén & Asparouhov, 2012). Instead of a point estimate for each parameter, in MCMC, the values of the posterior distribution are summarized with a measure of central tendency (i.e., mean, median, or mode) to provide a single representative estimate for the parameter. By providing a posterior distribution instead of point estimates, MCMC provides more intuitive analogs of standard errors or confidence intervals that do not rely on assumptions or asymptotic theory. Instead, the MCMC analog of the frequentist confidence interval is the Bayesian credible interval, for which the bounds of a 95% interval are the 2.5 and 97.5 percentiles of the posterior distribution. With MCMC, there are no null hypotheses for any effects; consequently, there are no p values.⁷ To determine whether an estimate is null in the population, a common test is whether 0 is within the 95% credible interval for the parameter of interest.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Daniel McNeish

Notes

References

Aguinis

Gottfredson

R. K.

Culpepper

S. A.

(2013). Best-practice recommendations for estimating cross-level interaction effects using multilevel modeling. Journal of Management, 39, 1490–1528.

Asparouhov

Hamaker

E. L.

Muthén

(2018). Dynamic structural equation models. Structural Equation Modeling, 25, 359–388.

Asparouhov

Muthén

. (2016). General random effect latent variable modeling: Random subjects, items, contexts, and parameters. In Harring

J. R.

Stapleton

L. M.

Beretvas

S. N.

(Eds.), Advances in multilevel modeling for educational research: Addressing practical issues found in real-world applications (pp. 163–192). Information Age.

Barnes

C. M.

Reb

Ang

(2012). More than just the mean: Moving to a dynamic view of performance-based compensation. Journal of Applied Psychology, 97, 711–718.

Bolger

Davis

Rafaeli

(2003). Diary methods: Capturing life as it is lived. Annual Review of Psychology, 54, 579–616.

Brunton-Smith

Sturgis

Leckie

. (2017). Detecting and understanding interviewer effects on survey data by using a cross-classified mixed effects location–scale model. Journal of the Royal Statistical Society: Series A, 180, 551–568.

Chen

Kirkman

B. L.

Kanfer

Allen

Rosen

(2007). A multilevel study of leadership, empowerment, and performance in teams. Journal of Applied Psychology, 92, 331–346.

Cleveland

W. S.

Denby

Liu

(2002). Random location and scale effects: Model building methods for a general class of models. Computing Science and Statistics, 32, 3–10.

Croon

. (2002). Using predicted latent scores in general latent structure models. In Marcoulides

Moustaki

(Eds.), Latent variable and latent structure modeling (pp. 195–223). Lawrence Erlbaum.

10.

Dalal

R. S.

Meyer

R. D.

Bradshaw

R. P.

Green

J. P.

Kelly

E. D.

Zhu

(2015). Personality strength and situational influences on behavior: A conceptual review and research agenda. Journal of Management, 41, 261–287.

11.

DeRue

D. S.

Hollenbeck

Ilgen

Feltz

(2010). Efficacy dispersion in teams: Moving beyond agreement and aggregation. Personnel Psychology, 63, 1–40.

12.

Devlieger

Rosseel

(2017). Factor score path analysis. Methodology, 13, 31–38.

13.

Estabrook

Grimm

K. J.

Bowles

R. P.

(2012). A Monte Carlo simulation study of the reliability of intraindividual variability. Psychology and Aging, 27, 560–576.

14.

Gelman

Rubin

D. B.

(1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457–472.

15.

Hedeker

Mermelstein

R. J

. (2007). Mixed-effects regression models with heterogeneous variance: Analyzing ecological momentary assessment data of smoking. In Little

T. D.

Bovaird

J. A.

Card

N. A.

(Eds.), Modeling ecological and contextual effects in longitudinal studies of human development (pp. 183–206). Erlbaum.

16.

Hedeker

Mermelstein

R. J.

Berbaum

M. L.

Campbell

R. T.

(2009). Modeling mood variation associated with smoking: An application of a heterogeneous mixed-effects model for analysis of ecological momentary assessment (EMA) data. Addiction, 104, 297–307.

17.

Hedeker

Mermelstein

R. J.

Demirtas

(2008). An application of a mixed-effects location scale model for analysis of ecological momentary assessment (EMA) data. Biometrics, 64, 627–634.

18.

Hedeker

Nordgren

(2013). MIXREGLS: A program for mixed-effects location scale analysis. Journal of Statistical Software, 52, 1–38.

19.

Hoffman

(2007). Multilevel models for examining individual differences in within-person variation and covariation over time. Multivariate Behavioral Research, 42, 609–629.

20.

Hoffman

. (2015). Longitudinal analysis: Modeling within person fluctuation and change. Routledge

21.

Hoffman

(2019). On the interpretation of multivariate multilevel model parameters across different combinations of model specification and estimation. Advances in Methods and Practices in Psychological Science, 2, 288–311.

22.

Hox

J. J.

van de Schoot

Matthijsse

(2012). How few countries will do? Comparative survey analysis from a Bayesian perspective. Survey Research Methods, 6, 87–93.

23.

Jongerling

Laurenceau

J. P.

Hamaker

E. L.

(2015). A multilevel AR (1) model: Allowing for inter-individual differences in trait-scores, inertia, and innovation variance. Multivariate Behavioral Research, 50, 334–349.

24.

Kelcey

Cox

Dong

(2019). Croon’s bias-corrected factor score path analysis for small-to moderate-sample multilevel structural equation models. Organizational Research Methods. Advance online publication. https://doi.org/10.1177/1094428119879758

25.

Kozlowski

S. W.

Klein

K. J.

(2000). A multilevel approach to theory and research in organizations: Contextual, temporal, and emergent processes. In Klein

Kozlowski

(Eds.), Multilevel theory, research, and methods in organizations: Foundations, extensions, and new directions (pp. 3–90). Jossey-Bass.

26.

Kruschke

J. K.

Aguinis

Joo

(2012). The time has come: Bayesian methods for data analysis in the organizational sciences. Organizational Research Methods, 15, 722–752.

27.

Kuppens

Yzerbyt

V. Y.

(2014). When are emotions related to group-based appraisals? A comparison between group-based emotions and general group emotions. Personality and Social Psychology Bulletin, 40, 1574–1588.

28.

Lang

J. W.

Bliese

P. D.

de Voogt

(2018). Modeling consensus emergence in groups using longitudinal multilevel methods. Personnel Psychology, 71, 255–281.

29.

Lang

J. W.

Bliese

P. D.

Runge

J. M.

(2019). Detecting consensus emergence in organizational multilevel data: Power simulations. Organizational Research Methods. Advance online publication. https://doi.org/10.1177/1094428119873950

30.

Leckie

French

Charlton

Browne

(2014). Modeling heterogeneous variance-covariance components in two-level models. Journal of Educational and Behavioral Statistics, 39, 307–332.

31.

Lester

H. F.

Cullen-Lester

K. L.

Walters

R. W.

(2020). From nuisance to novel research questions: Using multilevel models to predict heterogeneous variances. Organizational Research Methods. Advance online publication. https://doi.org/10.1177/1094428119887434

32.

Lievens

Sanchez

J. I.

Bartram

Brown

(2010). Lack of consensus among competency ratings of the same occupation: Noise or substance? Journal of Applied Psychology, 95, 562–571.

33.

Lüdtke

Marsh

H. W.

Robitzsch

Trautwein

Asparouhov

Muthén

(2008). The multilevel latent covariate model: A new, more reliable approach to group-level effects in contextual studies. Psychological Methods, 13, 203–229.

34.

McNeish

(2016). On using Bayesian methods to address small sample problems. Structural Equation Modeling, 23, 750–773.

35.

McNeish

Hamaker

E. L.

(2020). A primer on two-level dynamic structural equation models for intensive longitudinal data. Psychological Methods. Advance online publication. https://doi.org/10.1037/met0000250

36.

McNeish

Matta

(2018). Differentiating between mixed-effects and latent-curve approaches to growth modeling. Behavior Research Methods, 50, 1398–1414.

37.

Melkonyan

Safra

(2016). Intrinsic variability in group and individual decision making. Management Science, 62, 2651–2667.

38.

Muthén

(2010). Bayesian analysis in Mplus: A brief introduction. Los Angeles, CA: Author. Retrieved from. https://www.statmodel.com/download/IntroBayesVersion%203.pdf

39.

Muthén

Asparouhov

(2012). Bayesian structural equation modeling: A more flexible representation of substantive theory. Psychological Methods, 17, 313–335.

40.

Pak

Kim

(2018). Team manager’s implementation, high performance work systems intensity, and performance: A multilevel investigation. Journal of Management, 44, 2690–2715.

41.

Preacher

K. J.

Curran

P. J.

Bauer

D. J.

(2006). Computational tools for probing interactions in multiple linear regression, multilevel modeling, and latent curve analysis. Journal of Educational and Behavioral Statistics, 31, 437–448.

42.

Rast

Hofer

S. M.

Sparks

(2012). Modeling individual differences in within-person variation of negative and positive affect in a mixed effects location scale model using BUGS/JAGS. Multivariate Behavioral Research, 47, 177–200.

43.

Rast

Zimprich

(2010). Individual differences in a positional learning task across the adult lifespan. Learning and Individual Differences, 20, 1–7.

44.

Raudenbush

S. W.

Bryk

A. S.

(2002). Hierarchical linear models: Applications and data analysis methods. Sage.

45.

Reb

Greguras

G. J.

(2010). Understanding performance ratings: Dynamic performance, attributions, and rating purpose. Journal of Applied Psychology, 95, 213–220.

46.

Rindskopf

(1984). Using phantom and imaginary latent variables to parameterize constraints in linear structural models. Psychometrika, 49, 37–47.

47.

Schultzberg

Muthén

(2018). Number of subjects and time points needed for multilevel time-series analysis: A simulation study of dynamic structural equation modeling. Structural Equation Modeling, 25, 495–515.

48.

Sitzmann

Ployhart

R. E.

Kim

(2019). A process model linking occupational strength to attitudes and behaviors: The explanatory role of occupational personality heterogeneity. Journal of Applied Psychology, 104, 247–269.

49.

Sosnowska

Hofmans

Lievens

. (in press). Assessing personality dynamics in personnel selection. In Rauthmann

J. F.

(Ed.), The handbook of personality dynamics and processes. Elsevier.

50.

Van De Schoot

Broere

J. J.

Perryck

K. H.

Zondervan-Zwijnenburg

Van Loey

N. E.

(2015). Analyzing small data sets using Bayesian estimation: The case of posttraumatic stress symptoms following mechanical ventilation in burn survivors. European Journal of Psychotraumatology, 6, Article 25216. https://doi.org/10.3402/ejpt.v6.25216

51.

Van de Schoot

Kaplan

Denissen

Asendorpf

J. B.

Neyer

F. J.

Van Aken

M. A.

(2014). A gentle introduction to Bayesian analysis: Applications to developmental research. Child Development, 85, 842–860.

52.

Walters

R. W.

Hoffman

Templin

(2018). The power to detect and predict individual differences in intra-individual variability using the mixed-effects location-scale model. Multivariate Behavioral Research, 53, 360–374.

53.

Williams

D. R.

Zimprich

D. R.

Rast

(2019). A Bayesian nonlinear mixed-effects location scale model for learning. Behavior Research Methods, 51, 1968–1986.

54.

Zhang

Zyphur

M. J.

Preacher

K. J.

(2009). Testing multilevel mediation using hierarchical linear models: Problems and solutions. Organizational Research Methods, 12, 695–719.

55.

Zyphur

M. J.

Oswald

F. L.

Rupp

D. E.

(2015). Rendezvous Overdue: Bayes Analysis Meets Organizational Research. Journal of Management, 41, 387–389.