Multilevel Latent Profile Analysis With Covariates

Person-Centered Approach in the Single Versus Multilevel Context

A person-centered approach focuses on identifying subpopulations of individuals who show different patterns of values for certain criteria variables (Muthén & Muthén, 2000). In addition, the person-centered approach can also be used to identify profiles that differ from one another at the level of the relations between variables (i.e., mixture regression, factor mixture). Compared with a variable-centered approach, which applies a single model to the whole sample to estimate a single set of parameters, a person-centered approach is able to relax the assumption of homogeneity in the population through the estimation of a distinct set of parameters applicable to the various subpopulations present in the sample.

The concept of finite mixture models refers to statistical methods that probabilistically assign individuals into subpopulations (McLachlan & Peel, 2000). Latent class analysis (LCA) and latent profile analysis (LPA) are special cases of finite mixture models that aim to identify naturally occurring latent subpopulations that are homogeneous in relation to selected criteria variables (McLachlan & Peel, 2000). Although LCA is sometimes presented as a general statistical framework for mixture models, it typically refers to identifying subpopulations based on categorical variables, whereas LPA performs the same function using continuous variables (Oberski, 2016). In the present demonstration, we focus on LPA, which aims to identify latent subpopulations presenting different configurations on a set of indicators (McLachlan & Peel, 2000; Sterba, 2013). Consequently, the covariance structure of the observed variables is explained via differences in their mean values between the latent subpopulations rather than via the associations between continuous variables, as is done in variable-centered modeling.

In organizational research, LPA has previously been applied to identify, for example, different employee profiles in relation to the three symptoms of burnout, that is, exhaustion, cynicism, and reduced professional efficacy (for a review, see Mäkikangas & Kinnunen, 2016), or different mind-sets (i.e., affective, normative, and continuance) of workplace commitment (for a review, see Meyer & Morin, 2016). Because LPA assumes observations to be independent of each other, it is not a suitable method for use with hierarchical data sets where individuals are nested within groups (e.g., within units/departments/teams in an organization; Asparouhov & Muthén, 2008; Henry & Muthén, 2010). In addition to violation of the independency assumption of observations, ignoring the presence of nested data structures may lead to inaccuracy in the classification of profiles (e.g., to inaccuracy in the number of latent classes and to biased standard errors; Asparouhov & Muthén, 2008; Chen, Kwok, Luo, & Willson, 2010; Luo & Kwok, 2009; Meyers & Beretvas, 2006; Moerbeek, 2004; Park & Yu, 2016; Vermunt, 2003, 2008). Recent evidence indicates that multilevel modeling is advised in the case of nested data, even if no group-level effects are expected or evident (Bliese, Maltarich, & Hendricks, 2017).

To overcome the obstacles caused by hierarchical data structures, multilevel LPA (MLPA) is needed. The benefit of MLPA is that it enables the different levels of the data to be taken into account in modeling the phenomena of interest, and thus answer research questions which are beyond the scope of single-level LPA. In the organizational sciences, for instance, MLPA enables investigation of the role played by qualitatively different work environments. Examples of such novel questions would be, how do burnout profiles vary across work departments or organizations? Or how can work departments be clustered according to the different proportions of individual burnout profiles? In multilevel modeling, the associations between variables are modeled at two or more levels (e.g., at the level of the individual and at the level of the work department). In variable-centered multilevel modeling the associations of variables are modeled at different levels (e.g., between job characteristics and turnover intentions at different hierarchical levels of the data; Mauno, De Cuyper, Tolvanen, Kinnunen, & Mäkikangas, 2014). In contrast, MLPA explains the observed covariance structure through latent profiles and relates the information to latent profiles at different levels of the data (e.g., to predict the proportions of profiles at different levels).

MLPA can be specified in many ways. Here, we focus on three specific types of MLPA. In the first model, Level 2 variations in the size of Level 1 profiles, the variation in the relative size of single-level profiles across the higher-level units (i.e., Level 2) is modeled. In the second model, Level 2 classes based on the relative frequency of the Level 1 profiles, latent classes of higher-level units (i.e., Level 2) based on the proportions of the single-level latent profiles are investigated. In the third model, Level 2 classes based on Level 1 ratings, latent classes of higher-level units are based on the variables observed at Level 1. In the following, we use the term profile to refer to latent subpopulations at the level of individuals (Level 1) and the term class to refer to subpopulations at the level of work departments (Level 2). Next, we explain in detail the steps to follow in constructing these models.

MLPA: Steps and Equations

Modeling MLPA follows six distinguishable steps: (a) definition of the problem and selection of the variables, (b) preliminary analysis of the data, (c) model specification, (d) model estimation for 1, 2,…, K latent profile (Level 1)/class (Level 2) solutions, (e) choosing the best-fitting model, and (f) investigating the role of potential covariates.

In the first step, the definition of the problem to be solved and selection of the variables are guided by the research target and questions. MLPA, like other person-centered methods, is typically conducted in an exploratory manner, that is, a priori hypotheses for the number or nature of the latent upper-level classes present in the data may not be given (Finch & French, 2014; see also Meyer & Morin, 2016; Morin & Wang, 2016). However, predictions of the characteristics of possible subpopulations can be theoretically anchored, keeping in mind that person-centered analysis is more model- than theory-driven (for exceptions, see Finch & Bronk, 2011; Wang, 2007). Thus, instead of a theory-driven definition of the problem, the theoretical rationale for person-centered studies is typically related to the following questions: Why is there a need for the identification of possible latent subpopulations in relation to the variables of interest? Why can we expect that the identification of these profiles will improve our understanding of the phenomenon under study? Do we expect these profiles to have a practical impact on different covariates?

After defining the problem and selecting the relevant variables, the second step, preliminary examination of the data, follows. This includes, for example, checking for possible outliers (e.g., by using Mahalanobis distance), ¹ evaluating the reliability of the variables, and conducting missing data analysis. In implementing MLPA, it is also necessary to calculate intraclass correlation coefficients (ICCs) to see whether and how strongly the responses given by the members of the same subpopulation are related (Bliese et al., 2017). A significant ICC(1), indicating the degree to which individuals within a designated unit (i.e., work department) are nonindependent (Bliese, 2000), is often regarded as an essential precursor for multilevel modeling. However, recent evidence indicates that multilevel models are appropriate even in the case of a low or nonsignificant ICC(1) (Bliese et al., 2017). In our demonstration, ICCs are calculated by dividing the variance between work departments with the total variance (variance within + variance between work departments), which corresponds to ICC(1).

The third step, specification of the model, needs to be well-aligned with the research question and with current knowledge on best practices. In this stage, it is decided which of the model parameters (e.g., means, variances) are to be freely estimated and which are to be fixed/constrained at different levels of the data (Asparouhov & Muthén, 2008; Vermunt, 2011). Simulation studies conducted for different mixture models (i.e., LPA, growth mixture modeling, mixture regression modeling) suggest that to achieve a less biased and proper representation of the data, the variances across profiles should be freely estimated (Diallo, Morin, & Lu, 2016; Enders & Tofighi, 2008; Kim, Lamont, et al., 2016; Peugh & Fan, 2013). However, despite this generic recommendation, it is not always possible to rely on the results of models in which both means and variances are freely estimated. The reason for this is that these models tend more frequently to converge on improper solutions, which is generally taken to indicate that the model might have been overparameterized and that a more parsimonious model should be preferred (Bauer & Curran, 2003; Chen, Bollen, Paxton, Curran, & Kirby, 2001; Diallo et al., 2016). It should be also noted that LPA relies on a conditional independence assumption (i.e., variables within latent profile are uncorrelated). Consequently, covariances within profiles should be estimated only with great caution and on the basis of a strong theoretical justification (see Meyer & Morin, 2016). In MLPA, raw data or grand-mean centering is recommended, as they both provide the same results. To aid interpretation of the results, either grand-mean centering could be applied or the variables could be standardized for the analysis. Group-mean centering, however, should not be used in MPLA because it can substantially change model interpretation (see also Bliese et al., 2017).

To conclude, model specification and tests of alternative models should be conducted after careful reasoning and described in detail. As stated earlier, we present three possible ways to specify MLPA, each of which, depending on the research question, can be implemented separately. However, two of these (Level 2 variations in the size of the Level 1 profiles and Level 2 classes based on the relative frequency of the Level 1 profiles) are based on single-level profiles, and therefore it is recommended to conduct single-level LPA before MLPA.

The model for Level 2 variations in the size of the Level 1 profiles is performed to investigate whether the proportional distributions of the single-level profiles vary across the upper-level units. The model equation with two levels in which latent profiles j = 1, 2,…, J (across individuals i in k = 1,., K variables) appear at Level 1 and the proportions of the latent profiles vary randomly across Level 2 units u is

Level 1 : y i j k u = μ j k + ε i k ε i k ∼ N ( 0 , σ k 2 )

Probability of latent class j P ( y i j k u ) = p j u is associated to unit via latent variables

C 1 u , C 2 u , … , C J u ;

Level 2 : C j u ∼ N ( μ j , σ j 2 ) , j = 1 , … , J − 1 .

Mean μ_j and variance σ j 2 of the latent categorical profile variable C_ju are estimated in Level 2. The mean proportion p_j of latent profile j is obtained via the value of C_ju, where p j = e μ j 1 + ∑ j = 1 J − 1 e μ j , j = 1,.., J-1 and p J = 1 1 + ∑ j = 1 J − 1 e μ j . Consequently, this multilevel mixture model examines whether the probability that an employee belongs to a specific job characteristics profile varies significantly across work departments. In other words, whether in some work departments there is a larger probability that an employee belongs to a certain profile than in other work departments. The model can be extended by regressing C_ju on the Level 1 and/or Level 2 variables. The expected value of μ_cj of C_ju is then dependent on the Level 1 and/or Level 2 variables and thus also explains part of the variance σ c j 2 .

The model for Level 2 classes based on the relative frequency of the Level 1 profiles investigates whether Level 2 latent classes can be identified based on the proportions of the Level 1 latent profiles. Consequently, the equation for a two-level model in which latent profiles appear in Level 1 and the proportions of the Level 1 latent profiles are regressed on the Level 2 latent classes is

Level 1 : y i j k u = μ j k u + ε i j k u ε i j k u ∼ N ( 0 , σ j k 2 )

Probability of latent class j P ( y i j k u ) = p j u is associated with unit via latent variables

C 1 u , C 2 u , … , C J u ;

Level 2 : C j u | b = i n t e r c e p t j + β j b × C B b , j = 1 , … , J − 1 , b = 1 , … , B − 1 .

The value of C j u | b depends on which of the between latent profiles the C_Bb unit u belongs to. The proportion of latent profiles in latent class C_Bb is obtained via a regression equation where p j | b = e i n t e r c e p t j + β j b × C B b 1 + ∑ j = 1 J − 1 e i n t e r c e p t j + β j b × C B b , j = 1,2,., J-1 and p J | b = 1 1 + ∑ j = 1 J − 1 e i n t e r c e p t j + β j b × C B b . Similarly, the proportion of C_Bb in the whole data is p b = e C B j b 1 + ∑ k = 1 K − 1 e C B j b , k = 1,., K-1 and p K = 1 1 + ∑ k = 1 K − 1 e C B j b when b = B. To concretize, this multilevel mixture model extends the single-level LPA to the multilevel context (in our example, nesting employees in work departments) by specifying latent profiles for both employees (Level 1) and work departments (Level 2). The model can be extended by regressing C j u | b on the Level 1 and/or Level 2 variables and regressing C_Bb on the Level 2 variables.

The equation for the Level 2 classes based on Level 1 ratings with two levels in which latent classes j = 1,2,…, J (across individuals i in k = 1,., K variables) appear at Level 2 is

Level 1: y i j k u = μ j k u + ε i j k u ε i j k u ∼ N ( 0 , σ j k 2 )

Level 2: μ j k u = α j k

Probability p_j of latent class j is received via latent variable CBj = μ_j p j = e μ j 1 + ∑ j = 1 J − 1 e μ j , j = 1,., J-1 and p J = 1 1 + ∑ j = 1 J − 1 e μ j . Consequently, this model identifies the higher-level classes based on the Level 1 ratings of certain criteria variables.

To sum up, both models, Level 2 classes based on the relative frequency of the Level 1 profiles and Level 2 classes based on Level 1 ratings, can be used where the interest is in classifying upper-level units (in this case, work departments). The main difference between these models is that Level 2 classes based on the relative frequency of the Level 1 profiles investigates latent classes based on within-person profiles, whereas Level 2 classes based on Level 1 ratings examines upper-level classes on the basis of variables observed on the single level. Consequently, model selection depends on whether the research target is to investigate latent subpopulations simultaneously at both levels of the data, or focus solely on upper-level latent classes. In turn, in the model Level 2 variations in the size of the Level 1 profiles, upper-level latent subpopulations are not investigated. Instead, the main research question is to investigate intercontextual variability in the Level 1 profiles, that is, whether the rate of individual (i.e., single-level) experiences vary by higher-level units.

In the fourth step, the estimation of latent profiles (Level 1)/classes (Level 2) for alternative models is conducted, starting from one profile/class solution and continuing to K/B latent profile/class solutions. Because mixed normal distributions of observations are assumed, and typically data appear to be missing, the full information maximum likelihood estimator with robust standard error estimation is needed (Asparouhov & Muthén, 2008). While trying to find the maximum likelihood value, the algorithm may prematurely stop and return a suboptimal set of parameter values, that is, lead to local maxima. To overcome the local maxima problem and find the global maximum (i.e., optimal solution), it is necessary to use several sets of starting values for parameters (Asparouhov & Muthén, 2008; McLachlan & Peel, 2000; see also Li, Harring, & Macready, 2014). It should also be noted that different MLPA models may require different starting value sets.

In the fifth step, the final model is selected. Selection of the optimal number of profiles should be based on the inspection of several statistical information criteria (e.g., Akaike information criterion [AIC], consistent AIC [CAIC], Bayesian information criterion [BIC], and the sample-size-adjusted BIC [aBIC]) and tests (e.g., Lo-Mendell-Rubin likelihood ratio test [LMR] and the bootstrap likelihood-ratio test [BLRT]), the statistical adequacy of the solution (e.g., convergence and absence of negative variance estimates), and inspection of the substantive meaning and theoretical conformity of the solution (McLachlan & Peel, 2000; see also Morin, Morizot, Boudrias, & Madore, 2011). Simulation studies indicate that the BIC, adjusted BIC, consistent AIC, and BLRT tests appear to be particularly effective in choosing the model that best covers the sample’s true parameters on the single level (Henson, Reise, & Kim, 2007; Nylund, Asparouhov, & Muthén, 2007; Tein, Coxe, & Cham, 2013; Tofighi & Enders, 2008; Tolvanen, 2007; Yu & Park, 2014). For deciding the number of upper-level latent classes, the BIC has proved to be the most efficient criterion (Finch & French, 2014; Yu & Park, 2014; Zhang, Zhang, Zhang, & Jiou, 2014). Graphical presentations of alternative solutions also aid the selection of the final model. Moreover, as recommended by Meyer and Morin (2016), the final model results, based on both standardized and raw scores, are also presented graphically.

In the sixth and final step, it is possible to extend MLPA to include single- and upper-level covariates and to investigate their cross-level interaction. According to the latest knowledge, the covariates should not be incorporated into the models until after the final unconditional solution has been identified (Diallo & Lu, 2017; Diallo, Morin, & Lu, 2017; Kim, Vermunt, Bakk, Jaki, & Van Horn, 2016; Nylund-Gibson & Masyn, 2016). Furthermore, to avoid the problem of profile switching, the starts values from the final unconditional MLPA model should be used.

Having provided an overview of the steps of MLPA, we now demonstrate the method using a study example in which job characteristics profiles based on the JDCS model are analyzed with hierarchical data on university employees.

Data Used in the Demonstration

Steps and Results of the MPLA Demonstration

Step 2: Preliminary Analysis of The Data

Using Mahalanobis distance and corresponding p values, the preliminary analyses revealed the presence of deviant cases. At the p < .05 probability level, 5% of the cases were expected to be statistically significant at random. In our case, the percentage of statistically significant cases was 5.9%. However, in the cases with the smallest p values, further inspection revealed that the observed values of job characteristics were highly similar, suggesting the existence of a subpopulation. Therefore, all these deviant cases were included in the analysis (see also Note 1). The preliminary analyses showed further that the reliabilities of the observed variables were acceptably high (i.e., > .73; see the Measures section). ICCs (together with descriptive information on the variables) are presented in Table 1. As shown, the ICCs were significant for each of the four job characteristics variables, thus demonstrating the existence of meaningful variation between the Level 2 units (i.e., work departments). Significant ICCs mean that as well as varying from one individual to another, perceived job characteristics also varied from one work department to another.

OPEN IN VIEWER

Table 1.

Means (M), Standard Deviations (SD), Intraclass Correlations (ICC), and Correlations Between Study Variables at Level 1 (below the diagonal, N = 1,958 employees) and at Level 2 (above the diagonal, N = 78 departments).

Variable	M	SD	ICC	1	2	3	4	5	6	7	8
1. Workload	3.44	0.18	.05**	—	–.59	–.48	–.07
2. Job control	3.52	0.27	.14**		—	.11	–.05
3. Supervisory support	3.54	0.19	.03*			—	.73
4. Colleague support	3.90	0.19	.04**				—
5. Workloada	—	0.76						—
6. Job controla	—	0.68						–.33	—
7. Supervisory supporta	—	1.15						–.18	.30	—
8. Colleague supporta	—	0.89						–.17	.27	.53	—

^aLevel 1 variable.

*p < .05, **p < .01.

Step 5: Choosing and Interpreting the Best-Fitting Model

Single-Level LPA

The LPAs were first estimated with unequal means and variances across the profiles as recommended in the literature (Diallo et al., 2016; Enders & Tofighi, 2008; Kim, Lamont, et al., 2016; Peugh & Fan, 2013; see the appendix for Mplus script; Model 1). This analysis produced a two-profile solution. The solutions with three or more profiles did not converge. Therefore, we next estimated the more parsimonious model, that is, a model where the latent profiles differed in the mean values of job characteristics (see the appendix for Mplus script; Model 1). Figure 1 presents the elbow plot of the information criteria (BIC, aBIC, AIC, and CAIC) for the different profile solutions of the single-level LPA based solely on the mean differences between the profiles. As illustrated, between the one- and two-profile solutions, all the information criteria decreased strongly, but remained at the same level up to the six-profile solution. After a six-profile solution, the information criteria again decreased compared with the previous solutions. When the k and k-1 profile solutions were compared, the p value of the BLRT test was < .001; consequently, in all the comparisons with the k-1 profile solutions, the null hypothesis was rejected.

OPEN IN VIEWER

Figure 1.

Elbow plot of information criteria for different Level 1 profile solutions.

Typically, the change (i.e., decrease) in ICs depicted via an elbow plot should be U-shaped, while the point after which the slope flattens out indicates the optimal number of profiles in the data (see, e.g., Morin, Meyer, McInerney, Marsh, & Ganotice, 2015; Wang, Morin, Ryan, & Liu, 2016). The shape of the elbow plot obtained here may reflect the fact that the observed variables were not strictly continuous. Inspection of the profile solutions indicated that the six-profile solution and solutions with more than six profiles contained profiles with very small memberships (i.e., 1-2% of cases). We further investigated the two- to five-solution profiles graphically. It became evident that in the two-profile solution, the profiles differed from each other qualitatively (i.e., the two profiles had distinct shapes). The three- to five-profile solutions contained similarly shaped profiles with small level differences in workload, job control and facets of social support. Thus, these profile-solutions did not reveal new job types, as depicted in the JDCS model. Also, the classification quality of the two-profile solution was high: the average latent profile probabilities for the most likely latent profile membership (AvePP) were .90 and .94. In the three- to five-profile solutions, the AvePPs were much lower, the lowest varying between .74 and .83.

The two-profile solution is shown in Figures 2 (standardized scores) and 3 (raw scores). The first profile contained 64% of the participants, and was characterized by high social support, high levels of job control and low levels of workload. Consequently, this profile was labeled “collective low-strain job.” The second profile comprised 36% of the participants, and was characterized by high levels of workload, low levels of job control and low levels of social support from both supervisor and colleagues. Thus, using the terminology of the JDCS model, this profile was labeled “isolated high-strain job” (Karasek & Theorell, 1990). Previous person-centered research has also identified these low- and high-strain profiles (Mauno et al., 2016), and facets of social support further specified them as collective and isolated, respectively.

OPEN IN VIEWER

Figure 2.

Level 1 two-profile solution with equal variances across profiles based on standardized scores.

OPEN IN VIEWER

Figure 3.

Level 1 two-profile solution with equal variances across profiles based on raw scores.

Simulation studies suggest that LPAs with unequal variances lead to a less biased and a more accurate representation than models with only freely estimated means (Diallo et al., 2016; Enders & Tofighi, 2008; Kim, Lamont, et al., 2016; Peugh & Fan, 2013). As the estimation with equal variances resulted in a two-profile solution, as shown above, we also interpreted a two-profile solution based on unequal variances across the profiles. According to the BIC, the model with unequal variances (19151.664) was statistically better than the model with equal variances (21210.085). The two-profile solution based on unequal variances is shown in Figure 4. The form of the profiles is the same as in the two-profile solution with equal variances across profiles, but the size of the profiles is noticeably different (64% vs. 41% and 36% vs. 59%). Below, the MLPAs are illustrated based on the two-profile solution with unequal variances across the profiles. ²

OPEN IN VIEWER

Figure 4.

Level 1 two-profile solution with unequal variances across profiles based on raw scores.

MLPA

Of the alternative MLPA models, Level 2 variations in the size of the Level 1 profiles was first tested (see the appendix for Mplus script; Model 2). This was done by setting the starting values according to the single-level LPA and the random starting values as zero, thereby ensuring that the Level 1 profiles would remain unchanged (see Mplus script for Model 2). The results of the random effect MLPA showed that the size of the latent profile varied significantly between work departments (variance estimate = 0.322, SE = 0.162, p = .047), suggesting that the typicality of the “collective low-strain job” and “isolated high-strain job” profiles varied from one department to another. Overall, these results indicate that, depending on the university department, the proportional distributions of the two individual job characteristics profiles are different.

The model of Level 2 classes based on the relative frequency of the Level 1 profiles was then estimated (see the appendix for Mplus script; Model 3). We compared the different Level 2 latent class solutions using the BIC criterion (Finch & French, 2014; Yu & Park, 2014; Zhang et al., 2014). The corresponding values for the one- to three-class solutions were 19139.26, 19137.84, and 19473.53. Consequently, the two-class solution was deemed sufficient. These classes are presented graphically in Figure 5. The results showed that the first Level 2 latent class included departments in which the great majority of employees were characterized by the profile “isolated high-strain job” (8% of all work departments); these departments were thus labeled “isolated high-strain departments.” The second Level 2 latent class contained work departments which had relatively equal proportions of “isolated high-strain job” and “collective low-strain job” profiles (92% of all work departments); the Latent class 2 was thus labeled as “mixed strain departments.” In this class solution (two Level 1 profiles, two Level 2 latent classes), the multinomial regression coefficient obtained from regressing the Level 1 latent profiles on the Level 2 latent classes was 2.314 (SE = 0.815), p = .002. This regression coefficient indicates that in Level 2 latent class 1, the probability for the “isolated high strain job” profile is .936 and for the “collective low strain job” profile .064. In Level 2 latent class 2, the corresponding probabilities are .590 and .410.

OPEN IN VIEWER

Figure 5.

Level 2 classes based on the relative frequency of the Level 1 profiles.

Finally, a model of Level 2 classes based on Level 1 ratings was estimated (see the appendix for Mplus script; Model 4). In the one-class solution, the BIC was value 22028.68 and in the two-class solution 22037.56. Thus, it seems that no Level 2 latent classes were present in the data when the MLPA was based on Level 1 observations.

Discussion

This study presented and illustrated the use of MLPA, a method that offers an interesting extension of single-level LPA as it allows estimation of the distribution of Level 1 profiles across Level 2 units and the estimation of Level 2 classes based on either single-level profiles or variables. These methods have been successfully applied in various disciplines. For instance, they have been used to analyze categorical variables in the school context in the educational sciences (e.g., Allison, Adlaf, Irving, Schoueri-Mychasiw, & Rehm, 2016), and regional health behaviors in the health sciences (e.g., Erickson, Lenk, Toomey, Nelson, & Jones-Webb, 2016), and even for clustering scientific disciplines according to research output profiles (Mutz, Bornmann, & Daniel, 2013). However, in the organizational sciences the use of multilevel mixture modeling continues to be uncommon, despite data sets that are often hierarchical; for example, employees are nested in work units and teams are nested in organizations. Moreover, in the organizational sciences the variables analyzed are typically continuous rather than categorical, and thus MLPA is a highly advantageous modeling technique in this field of inquiry.

The consequences of not taking the nested structure of the data into account, such as inaccuracy in parameter estimation and biased standard errors due to using the Level 1 model alone, have been widely discussed and demonstrated (Asparouhov & Muthén, 2008; Chen et al., 2010; Luo & Kwok, 2009; Meyers & Beretvas, 2006; Moerbeek, 2004; Park & Yu, 2016; Vermunt, 2003, 2008). Thus, it is essential to recognize whether the data have been collected in a nested context, and to use multilevel modeling even in the case of low ICCs (Bliese et al., 2017). Conversely, aside from avoiding inaccuracy in profile estimation, the existence of a higher-level nested structure offers interesting novel analytical possibilities, as this study has demonstrated.

First, we showed in our demonstration that the distribution of profiles of individual latent job characteristics varied between university departments. That is, in some work departments the “collective low-strain job” profile was more typical, whereas in others the “isolated-high strain job” profile was more typical. Second, we were able to demonstrate in some detail how work departments were clustered according to the proportions of single-level profiles. In our data demonstration, two Level 2 latent classes were identified: “collective low-strain” work departments, composing most (92%) of the departments studied and “isolated high-strain” departments, accounting for the remaining 8%. To conclude, although job characteristics describe individual experiences of job control, demands, and support, a high-strain job profile was more typical in some work departments than others, and the identification of these departments was achieved through the application of MLPA.

Overall, our results indicate that the sample of a university employees studied typically belonged to a high-strain job profile. A novel finding is that the facets of social support did not reveal the additional job types described in the JDCS model (Karasek & Theorell, 1990). Thus, it seems that the original JDC model could be sufficient to characterize and differentiate types of jobs (Karasek, 1979), at least among highly educated university employees. Moreover, and perhaps more interesting, through the application of a multilevel mixture approach, we were able to identify latent classes of “healthy” and “unhealthy” work departments. The unhealthy “isolated high-strain” departments contained nearly all the employees with a high-strain job. This suggests that psychosocial job characteristics are not only individually experienced but may also be linked to collective stress and thus shared with colleagues.

Implications and Recommendations

MLPA importantly extends latent profile identification beyond the level of the individual employee, thereby enlarging the applicability of a person-centered methodology. Moreover, MLPA also enriches the outcomes of multilevel modeling. For example, in our case example, rather than simply estimating the amount of variance in job characteristics between work departments, it was possible to identify qualitatively different work departments. No hypotheses were set for the Level 2 latent classes, and thus they were identified in a purely data-driven way. In the context of multilevel person-centered research it would be advisable to use and combine both individual and organizational theories. Use of multiple theoretical frameworks also enables the possibility to conduct person-centered studies in a theory-driven way, as demonstrated by Wang (2007). Multilevel mixture modeling also has the potential to answer the generalization problem of person-centered research (see Meyer & Morin, 2016; Mäkikangas & Kinnunen, 2016), that is, whether the profiles found are unit-specific or more universal. Via the identification of Level 2 latent classes, it is possible to discover which profile type is predominant (versus underrepresented) across the higher-level units.

Besides offering theoretically innovative and interesting novel research targets, multilevel mixture modeling also provides plenty of modeling extension possibilities. In our MLPA demonstration, it was assumed that the Level 1 profiles are the same across the work departments, and vary only in their relative proportions. However, it is recommended that the Level 1 profile solution be verified when estimating the alternative MLPA models. In the longitudinal context, multilevel mixture models can be used to investigate Level 1 and Level 2 growth trajectories simultaneously (Chen et al., 2010). It is also worth noticing that all the multilevel mixture models can be based on latent factors instead of the sum scores used in this demonstration. However, too much complexity in the model easily induces convergence problems, and hence it might be reasonable to start the modeling procedure with a simpler model and increase the level of complexity within the boundaries of the tested model and sample size.

Overall, multilevel mixture models are highly recommended for analyzing nested data sets where the research interest is in studying latent classes in within-person profiles or their intercontextual variability (see also Henry & Muthén, 2010; Vermunt, 2003, 2008, 2011). One challenge in LPA—as in other types of mixture modeling—is to find the best fit solution, that is, the one showing the highest likelihood. It is typical of this kind of analysis that the iteration process stops at the local maximum likelihood value. Therefore, given the aim of finding the correct number of latent profiles, the use of several random sets of starting values is recommended (see also Meyer & Morin, 2016).

Although MLPA models offer indisputable benefits in modeling latent subpopulations in nested data sets, they also have some limitations. First, as models of these kinds are relatively complex and computationally heavy, they require a large sample size at both levels: a large number of observations at Level 1 and a large number of Level 2 clusters. Large sample sizes are especially required when the latent subpopulations are multifaceted and less distinct (Park & Yu, 2017). For example, the present data included 78 departments, of which only 8% were in the “isolated high-strain” class. An insufficient sample size in mixture modeling can lead to convergence issues, improper solutions, and the inability to identify small but meaningful subpopulations (Berlin, Williams, & Parra, 2014). However, a sufficient sample size is difficult to determine, as it depends, for example, on the size of the model, the distribution of the variables, the amount of missing data, the reliability of the variables, and the strength of the relations among the variables (Muthén & Muthén, 2002). Consequently, further statistical research is needed on the calculation of statistical power and adequate sample size for different MLPA models.

Second, a simplification of the model may be needed, such as in the case of a large number of Level 1 latent profiles whose random means correlate highly at Level 2. In this case, using a common factor to model the different random means at Level 2 may be a solution, as recommended by Henry and Muthén (2010). Third, modeling the Level 2 latent subpopulations and investigating the effects of the Level 2 covariates requires careful planning and theoretical reasoning prior to the data collection. For example, in the present demonstration, the Level 2 covariates were nonsignificant and their number was rather limited.

In sum, multilevel mixture modeling offers a fruitful approach to the analysis of hierarchical data in person-centered research, which is nowadays more commonly undertaken in the organizational sciences. Given the importance attributed to understanding human behavior in organizational settings, this method allows latent classes on different hierarchical levels (e.g., individual, team, and organization) in organizations to be identified at the same time, thereby also disentangling the diverse effects of covariates expressed on the levels of interest. Via studying the effects of various contextual factors on lower- and upper-level classes, this method has potential for broader application in organizational research. Owing to their suitability for studying, in particular, the effects of contextual-level predictors on individual typologies of behavior, we believe that these advanced statistical methods will encourage the formulation of new research questions, lead to new theorizing, and eventually also lead to new answers to issues pertaining to the employee-organization relationship.