Multilevel SEM with random slopes in discrete data using the pairwise maximum likelihood

Pairwise maximum likelihood (PML) estimation is a promising method for multilevel models with discrete responses. Multilevel models take into account that units within a cluster tend to be more alike than units from different clusters. The pairwise likelihood is then obtained as the product of bivariate likelihoods for all within-cluster pairs of units and items. In this study, we investigate the PML estimation method with computationally intensive multilevel random intercept and random slope structural equation models (SEM) in discrete data. In pursuing this, we first reconsidered the general ‘wide format’ (WF) approach for SEM models and then extend the WF approach with random slopes. In a small simulation study we the determine accuracy and efficiency of the PML estimation method by varying the sample size (250, 500, 1000, 2000), response scales (two-point, four-point), and data-generating model (mediation model with three random slopes, factor model with one and two random slopes). Overall, results show that the PML estimation method is capable of estimating computationally intensive random intercept and random slopes multilevel models in the SEM framework with discrete data and many (six or more) latent variables with satisfactory accuracy and efficiency. However, the condition with 250 clusters combined with a two-point response scale shows more bias. K


I N TRODUC T ION
Structural equation modelling (SEM) is an effective modelling framework used for measuring latent variables through indicators and studying the associations between the latent variables and/or observed variables (see Bollen, 1989;Kline, 2015). The generality of SEM is reflected in the number of applications in the social sciences (e.g., Guo et al., 2008;Kline, 2015;MacCallum & Austin, 2000;Xiong et al., 2015). SEM was introduced for continuous data and later extended to handle more complex data.
Here, we will focus on the pairwise maximum likelihood (PML) estimation method for SEM to deal with two data complexities: discrete responses (i.e., binary coding or three-or four-point scales) and multilevel structures.
To estimate discrete data in the context of SEM, the PML estimation method was introduced by Jöreskog and Moustaki (2001). With PML estimation the product of bivariate (and sometimes univariate) likelihoods is calculated (see Jöreskog & Moustaki, 2001). Outside the SEM context, PML is part of a broader framework of composite maximum likelihood estimators (see Varin, 2008;Varin et al., 2011). Katsikatsou et al. (2012) showed that PML estimation produced low biases within the SEM framework. Compared to other frequentist SEM estimation methods for discrete data, the most prominent advantage of the PML estimation method is the possibility of computing models with a large number of latent variables. The marginal maximum likelihood (MML; see Bock & Aitkin, 1981) estimation method calculates a full likelihood and cannot estimate models with too many latent variables. Alternatively, the (weighted) least squares estimation method (Browne, 1984;Muthén et al., 1997) can be used to estimate discrete data. This estimation method also uses bivariate and univariate information and is in that respect comparable to the PML estimation method.
Estimating models with both discrete data and a multilevel structure complicates the estimation. In multilevel data, lower-level units are selected within higher-level units (i.e., clusters). As a consequence, units within a cluster tend to be more alike than units from different clusters. This dependency in the data introduces an additional source of variation that needs to be taken into account in analysing data (Hox et al., 2017). Multilevel data in SEM are usually analysed with the 'long format' (LF) approach (see McDonald & Goldstein, 1989;Muthén, 1990), where each row corresponds to the data of a single unit and multiple rows constitute a single cluster. In the 'wide format' (WF) or 'multivariate' approach, each data row is independent and corresponds to a single cluster (see Barendse & Rosseel, 2020, for a multilevel model framework with a random intercept). However, the number of columns can increase substantially with large clusters in the WF approach. The WF approach is therefore particularly useful for data with relatively small cluster sizes (e.g., Koomen et al., 2007;Lau et al., 2015;Mahlke et al., 2016;Moorman, 2016;NLSAH, 2005). This paper combines the challenges of multilevel structures including random slopes and discrete data using the PML estimation method. The pairwise likelihood in multilevel data is obtained as the product of bivariate likelihoods for within-cluster pairs of units and items. So far, the PML estimation method with random intercepts and random slopes has only been used for generalized mixed models for binary data (see Bellio & Varin, 2005;Cho & Rabe-Hesketh, 2011;Renard et al., 2004;Tibaldi et al., 2007). Random slopes show up in multilevel models if the effect of covariates is allowed to vary across clusters. To deal with the random slopes, casewise estimation is applied. Casewise estimation blurs away most of the differences between the WF and LF approaches. The aim of this study is to estimate complex multilevel random intercept and random slope SEM with at least six latent variables in the WF approach that cannot be estimated with other frequentist estimation methods. Computationally, these models cannot be estimated with intensive numerical (e.g., adaptive Gauss-Hermite quadrature) multilevel marginal maximum likelihood estimation methods as the number of latent variables that has to be integrated out is too high (MML; see Hedeker & Gibbons, 1994). In addition, these models cannot be estimated with multilevel weighted least squares (WLS; see Asparouhov & Muthén, 2007), as it does not allow for the casewise estimation that is essential for models with random slopes.
The paper is organized as follows. We first describe the general SEM framework and the PML estimation method for discrete data. Then we introduce a new procedure to analyse data in the WF

| 329
MULTILEVEL SEM WITH RANDOM SLOPES approach in multilevel data with random slopes and show how PML estimation can be used within this framework if the data are discrete. Next, we perform a small simulation study using two complex multilevel models (i.e., a mediation model and a factor model) with random slopes to evaluate the PML estimation method in terms of the accuracy of the parameter estimates. Finally, we discuss the usefulness of the PML estimation method in multilevel data.

General SEM framework
Structural equation modelling can estimate a wide range of models. It can be considered as a combination of regression or path analyses and factor analysis. The general model for unit i with continuous data can be described by a measurement model, that relates measured variables to latent variables, and a structural model, that relates latent variables to one another: where y i is a pdimensional vector of observed variables for unit i, is a pdimensional vector of intercepts, Λ is a p × l matrix of factor loadings relating the observed variables to the ldimensional latent variables i , i is a pdimensional vector of measurement errors or residuals, is an ldimensional vector of latent factor means and intercepts, B is an l × l matrix of regression coefficients among the latent factors, and i is an ldimensional vector of residuals for unit i. Assuming Cov( , ) = 0, Cov( , ) = 0, E( i ) = 0, E( ) = 0, diag(B ) = 0, and that (I −B) is invertible, where I is an l × l identity matrix, we can find expressions for the p × p covariance matrix Σ i and the ldimensional vector of the mean structure i (Equation 4) of y i : where the variances and covariances of and are denoted by Ψ and Θ, respectively. The model parameter vector includes the free parameters in Λ, B, Θ, Ψ, , and . Before estimating models, the scale for the latent variables (i.e., ) needs to be defined, for example by fixing the first factor loading of each latent variable to unity.

Discrete data with the PML estimation method
To deal with discrete data, the PML estimation method assumes an underlying normally distributed continuous latent response variable. A variable y ig for individual i on item g with C g response scales stems from an underlying continuous variable y * ig with a normal distribution N y * ig | 0, 2 g and τ g,c values that refer to thresholds for categories c g = 1, 2, … , C g , with g,0 = − ∞ and g,C = + ∞. Instead of calculating a full likelihood, PML estimation breaks down the complex likelihood. The log-likelihood of the PML estimation method for an individual is then calculated as the sum of p ⋆ = p( p − 1)/2 components, each component being the bivariate log-likelihood of two variables (i.e., g and h): (1) The total log-likelihood of the data is the sum of all individual contributions in a random sample Y = {y 1 , y 2 , …, y I } of size I and equals The exact form of f y ig , y ih ; in Equation 7 for discrete indicators g and h equals with where gh is the model implied correlation between y * ig and y * ih , and Φ(τ 1 , τ 2 ; ) is the bivariate cumulative normal distribution with correlation evaluated at the point (τ 1 , τ 2 ). Before estimating the model with the PML method, the metric for y * needs to be determined. Two popular ways are to fix the total variance (socalled delta parameterization: Θ = Δ −2 − diag(Σ*), where Σ*= Λ (I -B) -1 Ψ(I -B) -1T Λ T and Δ the scaling factors) or to fix the residual variance (so-called theta parameterization: Δ −2 = diag(Σ*) + Θ). Other ways of scaling are also possible (see Lee et al., 1990Lee et al., , 1992. Research has shown that the PML estimation method provides accurate and efficient results in the SEM context (see Jöreskog & Moustaki, 2001;Katsikatsou et al., 2012) as well as in the broader framework of composite maximum likelihood estimators (see Lindsay, 1988;Varin, 2008). In theory, Equation 6 is general and can deal with any type of data (continuous or discrete, and combinations thereof). Barendse and Rosseel (2020) estimated SEM with the PML estimation method for a mixture of binary and continuous data, by estimating Pearson, tetrachoric, polychoric, and polyserial correlations. Standard errors and missing-data procedures for the PML estimation method have been developed by Katsikatsou et al. (2012) and Katsikatsou and Moustaki (2017).

Multilevel SEM in the WF approach
Before investigating the PML estimation method for multilevel data in the WF approach, we will first explain the WF framework for continuous data. With continuous multilevel data the LF and WF approaches obtain identical results (see Mehta & Neale, 2005). A random intercept, random slope regression model in the LF approach can be estimated in the WF approach as a single-level confirmatory factor analysis in the SEM framework (see Bryk & Raudenbush, 1987;Chou et al., 1998;MacCallum et al., 1997;McArdle & Epstein, 1987;Mehta & Neale, 2005;Meredith & Tisak, 1990). The random slope allows the relationship between the covariate and the dependent variable (y) to be different for each cluster. The mean of the random slope represents the fixed effect of the covariate and the variance of the random slope reflects the variability. Figure 1a and 1b show the data layout, a graphical representation of the model, and the formula of a random intercept and random slope multilevel linear regression model for covariate x in the LF approach and in the WF approach, respectively. In Figure 1b for the WF approach there is no single covariance matrix for the entire sample as with conventional SEM, but a cluster-dependent covariance matrix (Mehta & Neale, 2005), where Λ j is now cluster-specific. To estimate this model in the WF approach, we have to use casewise calculations to compute a likelihood with cluster-specific vectors, not unlike the full-information maximum likelihood approach that was introduced in SEM for handling missing data (Arbuckle, 1996). Consequently, a multilevel random intercept and random slope model can only be estimated with software that allows for casewise estimation to deal with a cluster-specific model implied covariance and mean structure.
F I G U R E 1 A random intercept and slope model with the corresponding data, figure and formula in both (a) the LF approach and (b) WF approach. Note: Triangles indicate for the mean or intercepts. In (1a), x ′ ij equals the fixed effects matrix and z ′ ij indicates the random effects matrix. In (1b), the rounded boxes represent the within-and between-components of the original variable y for three units per cluster. y corresponds to the unit-specific version of y at the between-level, which equals a factor with loadings fixed at unity. rs corresponds to a random slope at the between-level with indicators fitted to the values of x j . Residual variances (not shown here) contain equality restrictions to ensure only one residual variance is estimated.

1a) 1b)
A multilevel regression model can be extended to a multilevel SEM. Multilevel random intercept SEM in the LF approach is described by different authors (see Bryk & Raudenbush, 1987;McDonald, 1993;McDonald & Goldstein, 1989;Muthén, 1989Muthén, , 1990Schmidt, 1969). Mehta and Neale (2005), Curran (2003), and Bauer (2003) were among the first authors to describe random intercept SEM in the WF approach. Barendse and Rosseel (2020) described a general SEM WF framework for both continuous and discrete data and proposed steps to perform random intercept models. However, a complication arises in the case of a covariate that exists at the within-level and the between-level (see Lüdtke et al., 2008). In a more general sense, this also holds for every regression coefficient that exists at the within-level and the between-level. To obtain the correct between-level effect, one has to disentangle the between effect by estimating the between-level effect and subtract this from the within-level effect using a definition variable (see Lüdtke et al., 2008). To avoid using definition variables to calculate the between effect, we adjusted the estimation of models in the WF approach in a way that is more intuitive than the parameterization presented in Barendse and Rosseel (2020). In this new approach we disentangle the between-and within-levels by creating latent variables, which are then used to build models at the within-level and the between-level.

M E T HODS
We first reformulate the WF approach of Barendse and Rosseel (2020) to allow for a more intuitive model parameterization in the presence of a covariate or a regression coefficient at both the within and the between-level (see Lüdtke et al., 2008). Then we explain how to apply the PML estimation method with discrete data in the WF approach.

Reformulation of the WF approach
Based on the steps introduced by Barendse and Rosseel (2020), we formulate adjusted steps that first disentangle the between-and within-level of the model and then explicitly build models at both levels.
1. Rearrange the data in such a way that each row corresponds to a single cluster. 2. Disentangle each continuous endogenous variable into a between part (e.g., ) and a within part (e.g., y w ) by introducing new variables with factor loadings fixed to unity. For the between variables, representing the random intercepts of the model, one has to construct a latent variable where the indicators correspond to the unit-specific observations of that variable. 3. Construct a model with the newly introduced within-level variables involving the variables that belong to a single unit in a cluster and repeat this model as many times as the maximum cluster size. 4. Put equality constraints on all parameters across units in a cluster in the within part of the model. For example, in a one-factor model, equality constraints are necessary on the factor loadings, factor variances, and error variances. If variables are both at the within and the between-level, the intercepts at the within-level should be fixed to zero. 5. Construct a model at the between-level with the newly constructed between-level latent variables.
The solid lines in Figures 2 and 3 show a random intercept mediation model and a random intercept factor model, respectively. Appendix A shows the lavaan syntax for a multilevel mediation model (i.e., the model shown in Figure 2) and Appendix C shows the lavaan syntax for a multilevel factor model (i.e., the model shown in Figure 3).
This multilevel random intercept SEM can be extended with a random slope. Rockwood (2020) gives a detailed description of multilevel models in the LF approach with random slopes (e.g., Figure 6 in this paper shows a one-factor model in the LF approach with a random slope). In the WF approach | 333 MULTILEVEL SEM WITH RANDOM SLOPES we can also extend multilevel random intercept models with one or more random slope(s), by adding an additional step: 6. For each regression coefficient at the within-level that we allow to vary across clusters, we construct a latent variable ( rs ) at the between-level with the dependent variable as indicators and factor loadings that are fixed to the values of the covariate. The latent variable rs represents the random slope in the model (see the dashed lines in Figures 2 and 3, and Appendices B and D for additional syntax to include a random slope).

Multilevel discrete data with the PML estimation method in the WF approach
With multilevel data, the pairwise likelihood is obtained as the product of bivariate likelihoods for within-cluster pairs of units (see Renard et al., 2004). Equations 6-9 are applicable for multilevel data by substituting individual i for unit j, where each separate row represents a single cluster instead of a single unit (see Figure 1). If the PML estimation method is used with multilevel discrete data, only the within-level needs to be adjusted. Instead of estimating error variances at the withinlevel, we estimate thresholds with equality constraints. rs with the dependent variable as indicators and factor loadings that are fixed to the values of the covariate x j . Residual variances (not shown here) also contain equality restrictions at the within-level. Appendix A shows the corresponding lavaan syntax for the random intercept model and Appendix B shows the additional lavaan syntax to add a random slope. Appendix E shows the alternative lavaan syntax with discrete variables. generalized (random intercept, random slope) multilevel regression models with a binary outcome. For fitting multilevel SEM with discrete data, we need to adjust our formulated step 2: 2a. For each discrete endogenous variable at the within-level, construct a latent variable (i.e., y * ), where the indicators correspond to the unit-specific observation of that variable. The thresholds are estimated with equality constraints on thresholds across units within a cluster (see step 4). Figure 4 shows a generalized multilevel regression model with both a random intercept and a random slope and four response options. This model is similar to the one presented in Figure 1b, but now contains four response options instead of continuous data. Appendix E shows the lavaan syntax for a multilevel mediation model with discrete data and Appendix F shows the lavaan syntax for a multilevel factor model with discrete data. The syntax in Appendix B or D can be added to include a random slope.
F I G U R E 3 A multilevel factor model in the WF approach with three items and three units within each cluster. Note: The rounded boxes represent the within-and between-components of the original variables y 1.1 to y 3.3 for three units per cluster. correspond to the unit-specific observation of that variable at the between-level, where labels a, b, c reflect the between-level factor loadings. For identification purposes, one factor loading or the variance of fb must be fixed to unity. The structures below y 1.1 to y 3.3 account for the within-level latent variables (i.e., ) with parameter labels d, e, f for the factor loadings. For identification purposes, one factor loading of each within factor ( fw ) or the variance of fw must be fixed to unity. The dashed lines indicate a random slope rs with the dependent variable as indicators and factor loadings that are fixed to the values of the covariate x j . The parameter h is an indicator of the correlation between the between-level factor and the random slope. Identical parameter labels indicate equality constraints in the model. Residual variances (not shown here) also contain equality restrictions. Appendix C shows the corresponding lavaan syntax for the random intercept model, and Appendix D shows the additional lavaan syntax to add a random slope. Appendix F shows the alternative lavaan syntax with discrete variables.

SI M U L ATION ST U DY
To investigate computationally intensive multilevel SEM models with discrete data in the WF approach using the PML estimation method, we ran a simulation study with a multilevel mediation model and two multilevel factor models (i.e., one with one random slope and one with two random slopes). All the models have at least six latent variables (i.e., random slopes or factors) in the model. Within each type of model, we vary response scales (two-point, four-point) and the number of clusters (250,500,1000,2000). The cluster size is always fixed at three. In a fully crossed design, these factors yield 3 × 2 × 4 = 24 different conditions. The performance of the PML method is evaluated by calculating the relative bias. The relative bias is calculated for the estimated parameters as % bias = (θ −θ)/θ × 100 and for the standard errors as % bias = (SE-SD)/SD × 100, where SE is the mean of the estimated standard errors across replications and the SD refers to the standard deviation of the parameter estimates across replications. Figures 5 and 6 show the data-generation model in the LF approach for a mediation model and a factor model, respectively. The mediation model is generated with three random slopes and the factor model is generated with one or two random slopes. The dashed lines represent the random slopes (i.e.,

Data generation
). The solid lines of models in the LF approach are identical to those in the WF approach presented in Figures 2 and 3. In all models we generated more variance at the within-level than at the between-level to mimic real data examples (see Snijders & Bosker, 1999). For example, using the parameter values of the basic factor model from Figure 6, the intraclass correlation for all indicators equals ) refers to a random intercept, ( y w 1 ) refers to the within representation of the variables, refers to the thresholds, and y * refers to the underlying latent response variables. Identical parameter labels indicate equality constraints. Figure 4 is similar to Figure 1 with the new parameterization and discrete data.
(1 2 × .5 + .2)/((1 2 × .5 + .2) + (1 2 × 1 + 1)) = .260. The multilevel models are generated with three variables with three units in each cluster to limit the estimation time in the simulation study. Continuous data were drawn from a multivariate normal distribution. Discrete data with two-point response scales are obtained by discretizing the continuous data such that the population proportions equal approximately .50 per category and four-point response scales are obtained by discretiszing the data such that the population proportions equal approximately .16, .34, .34, and .16. We have chosen the mean of the random slope (see the triangles in Figures 5 and 6) to be equal to the variance of the random slope. The residual variances equal an identity matrix to resemble the theta parameterization. For each condition 500 data sets are generated.

Estimation
The R package lavaan (version 0.6-12: Rosseel, 2012) is used for the calculations. Casewise estimation is necessary to estimate multilevel models with random slopes. To perform casewise calculations in this study, we wrote custom scripts to estimate the likelihood for each cluster separately. 1 Based on the theory of Katsikatsou et al. (2012), we used a closed-form solution to calculate the PML standard errors casewise (i.e., cluster-specific). In the estimated factor models, we fixed the first factor loading to unity to set the metric of the factor. 1 All R scripts (e.g., data-generation scripts and scripts to analyse the data) are available at https://osf.io/346vj/.

Results
After applying each of the 12000 data sets, we found that all models converged. Inspection of the results shows that especially models with 250 clusters and two-point scales show unreasonable parameter estimates and standard errors. We removed all results that included at least one parameter or standard error at more than four standard deviations from the mean 2 . This resulted in 3.90% deletions for the mediation model with three random slopes, 2.78% deletions for the factor model with one random slope, and 4.38% deletions for the factor model with two random slopes. Below we describe the relative bias (expressed as percentages) of the estimated mediation model and the factor models for the within-and between-level separately. Notice that the coefficients displayed in the plots refer to Equations 3 and 4.
2 All raw results with and without outliers are available at https://osf.io/346vj/.

Within-level
Bias of the parameter estimates Figure 7 shows the relative bias of the three regression coefficients in the mediation model at the withinlevel. The relative bias is lower across all conditions with a larger number of clusters and a four-point response scale. In particular, the condition with 250 clusters with a two-point response scale influenced the overall results. Raw results show that all the parameter estimates across all conditions with a fourpoint scale have less than 1.3% bias and that all conditions with a two-point scale and 1000 or 2000 clusters have less than 1.5% bias.

Bias standard errors
To study the efficiency of the estimates, we calculated the relative bias of the standard errors. The relative bias of the standard errors associated with the within-level regression coefficients is shown in Figure 8. High percentages of relative bias were found in conditions with 250 clusters and a two-point response scale (about 70% to 275%). To put this in context, an absolute bias between the standard errors and the standard deviations of the parameter estimates that equals .09 resulted in about 100% bias. Conditions with a four-point response scale with at least 500 clusters and conditions with a two-point response scale with at least 1000 clusters were more efficient, with less than 6% bias and less than 8% bias, respectively.

Bias of the parameter estimates
The relative bias of the regression coefficients across all conditions in the mediation model is shown in Figure 9. Almost all conditions had very low percentages of relative bias. We only observed a little more . Figure 10 shows the relative bias of the parameters related to the variances and covariances in the mediation model. For the sake of efficiency, we only show the variance of z at the between-level (i.e., z). Results show that the random slopes have more bias than the regression coefficients. Figure 10 shows low percentages of bias across all conditions in the estimated covariances and higher percentages of bias in the variances. Overall, we observe similar patterns at the between-level to those we observed at the within-level, and more accuracy across all conditions with more clusters and a four-point response scale. Bias standard errors Figures 11 and 12 show the relative bias of the standard errors for all parameters at the between-level. Once again, high percentages of bias were found in conditions with fewer clusters and two-point response scales. The bias related to the variances is higher than the bias related to the regression coefficients. We identified similar patterns of high and low percentages of relative bias to those we found at the within-level. In particular, a sample size of 250 with a two-point response scale showed very high percentages of relative bias across all parameters (about 65% to 230%). All other conditions show reasonable percentages of relative bias.

Factor models with random slopes
As the result tendencies of the factor model with one random slope are very similar to the results with two random slopes, we will only show the results of the factor model with two random slopes. In general, the bias in the factor model with one random slope is a little lower across both parameter estimates and standard errors. The results of the factor model with one random slope are shown in the online supplementary materials.

Within-level
Bias of the parameter estimates Figure 13 shows the relative bias of one of the factor loadings (we chose the third factor loading) and the variance of the within-level factor in a factor model with two random slopes. The factor variance ( fw ) shows higher percentages of bias than the factor loading (λ w ). Overall, the percentages of relative bias were quite low. Raw results show that the condition with 250 clusters with a two-point response scale showed a maximum of 17% relative bias. All other conditions have a maximum of 5.5% relative bias.

Bias standard errors
The relative bias of the within-level parameter estimates is shown in Figure 14. Notice that we used a different scale on the yaxis as the relative percentages of bias were much lower in the factor models than those observed in the mediation model. The relative bias found across all conditions with a four-point response scale was <5%. In conditions with a two-point response scale a maximum of 15% relative bias was observed.

Bias of the parameter estimates
The relative bias of the parameter estimates is shown in Figures 15 and 16. The variances and covariances show higher percentages of bias than the other parameters. In general, we observe a similar trend with more bias in conditions with fewer clusters and fewer response scales. Except in conditions with a cluster size of 250 and a two-point response scale, we observe a maximum of 12% bias.

Bias standard errors
Figures 17 and 18 show the relative bias of the standard errors related to the parameter estimates at the between-level. The factor loading and the variance of the between-level show a pattern in which more relative bias is related to a smaller number of clusters. The bias of the standard errors related to the other parameter estimates is lower. Taking into account that the percentages of bias in the standard errors increases quite quickly, the percentages of relative bias were quite low. In addition, we observe here that calculating bias across all conditions slightly increases the bias of the standard errors. For example, the raw results with 1000 and 2000 clusters show no more than 10% relative bias.

DISCUS SION
In this study we combined the challenges of multilevel-and discrete-data SEM using the PML estimation method. In pursuing this, we first reformulated the general framework of multilevel SEM in the WF approach of Barendse and Rosseel (2020) into a parameterization that is more intuitive and avoids using definition variables in case covariates or regression coefficients exist at both the within-and between-level. In this approach, the between-and within-levels are disentangled via latent variables that are then used to specify models at each level separately. Then we extended the reformulated WF approach for multilevel models with random slopes. In this study, we specifically focused on estimating structural equation models with discrete data and many latent variables (six or more, including random slopes) as these are computationally too intensive to estimate with the multilevel marginal maximum likelihood estimation method. The least squares estimation methods cannot be used either because they do not allow for random slopes due to relying on summary statistics. The PML estimation method, on the other hand, is well suited for these kinds of models, as PML allows for random slopes and calculates the products of bivariate likelihoods that are computationally easy to handle.
We conducted a small simulation study to examine whether the PML estimation method can deal with complex multilevel data including random slopes. Results across all conditions show that the PML estimation method can deal with three-variable mediation models with three random slopes and threevariable factor models with one or two random slopes. The mediation model showed higher percentages of relative bias than the factor models. This may have to do with the number of random slopes. The simulation study also showed that a factor model with one random slope can easily be extended to a model with two random slopes, with only slightly less accurate and efficient results. In general, a larger number of clusters and four-point response scales leads to low percentages of relative bias. The condition with 250 clusters and a two-point response scale showed higher percentages of relative bias. More specifically, the most prominent bias of the standard errors is observed in the mediation model. Overall, we observed more bias at the within-level than at the between-level. Additional research to estimate the same multilevel model with scripts that do not make use of equality constraints at the within-level revealed that the equality constraints did not cause the additional bias. The bias at the within-level may be improved by increasing the number of units in each cluster. Generally speaking, we conclude that the PML estimation method is quite accurate and efficient when there are enough data to estimate the parameters. Only the condition with 250 clusters and a two-point response scale lacked sufficient accuracy and efficiency. Comparing our results to those obtained in Barendse and Rosseel (2020), we conclude that models with random slopes need a larger number of cluster to be accurate and efficient. The random slopes bring an extra level of variability to the model. More research is needed to investigate whether this can be improved by increasing the number of units in the cluster or using techniques such as bounded estimation (De Jonckere & Rosseel, 2022). In this study we used a different parameterization than in the study of Barendse and Rosseel (2020). This parameterization does not influence the results in most of the models. However, in models with regression coefficients that exist at both the within-and betweenlevel, the newly described parameterization is highly recommended as it causes fewer difficulties in model estimation (i.e., it avoids calculation the between-level with the within-level using definition variables). Compared to the previous approach, we estimate more latent variables in the newly described approach. However, the PML estimation is not affected by this as it can handle many latent variables.
Adding random slopes complicates the estimation of a multilevel structural equation model, but it is still a sophisticated method to deal with within-level covariates. We recommend that each within-level covariate is first investigated via a random slope. Ignoring a random slope creates a misspecification in the model. This misspecification can result in incorrect standard errors and incorrect fit statistics. In principle, there are different ways of incorporating covariates that can be combined when fitting multilevel models. In case a random slope does not contain enough variance, one can continue with a fixed covariate. If one is not interested in the effect of the covariate, it is also possible to treat the variable as an exogenous covariate and regress out the covariates first and perform all other calculations on the residual correlations (Katsikatsou, 2017). When the covariate is a dichotomy, it is also possible to investigate the effect of the covariate in a multigroup analysis. Depending on the model, one can chose and combine different ways to incorporate covariates.
Even though PML estimation can estimate complex multilevel models in the WF approach, there are a number of limitations. As the PML with multilevel data is obtained as the product of bivariate likelihoods for within-cluster pairs of units and variables, the PML in the WF approach is slower with large cluster sizes. Further research should investigate whether we can improve the estimation by, for example, deleting the bivariate pairs that do not contribute much to the likelihood or applying a two-step approach where the thresholds are fixed in the second step to reduce the number of estimated parameters.
Notwithstanding the difficulties of the PML estimation method in the WF approach, the PML estimation method seems promising and suited to fit complex multilevel structural equation models in the WF approach. The suggested two-level models with PML estimation methods can in theory also be extended to more than two levels. In addition, PML is computationally efficient enough to deal with many latent variables at both the within-and between-level. That does not mean that estimation is fast. Indeed, it currently takes a couple of hours to estimate a model. Another advantage of the PML estimation method with multilevel data is that the number of columns (variables × units) can be larger than the number of rows (number of clusters). In this simulation study we used three units in each cluster to limit the estimation time, but models with larger cluster sizes can be estimated. Finally, the PML estimation method can deal with all types of data -discrete, continuous, and combinations thereof.
In this paper we only focused on frequentist estimation methods and did not take into account the multilevel Bayesian estimation method (e.g., Fox, 2010), which is a full-information method that is able to estimate random slopes. For further research, it would be interesting to compare the PML estimation method to the multilevel Bayesian estimation method, which can also estimate complex multilevel models.

C ON F L IC T OF I N T ER E S T
All authors declare no conflict of interest.

DATA AVA I L A BI L I T Y S TAT E M E N T
Data has been shared (following the link https://osf.io/346vj/).

SU PP ORT I NG I N FOR M AT ION
Additional supporting information can be found online in the Supporting Information section at the end of this article.