SIMULATION STUDIES OF INFORMATION CRITERIA USED IN THE SELECTION OF MULTILEVEL STRUCTURAL EQUATION MODELS

ABSTRACT

This study examines the application of information criteria in the selection of multilevel structural equation models (MSEM), analyzing different estimation methods (ULS, GLS, and ML) and statistical criteria (AIC, BIC, BCC, and CAIC). The study defines MSEM and specifies fit functions for different estimation methods. Using Monte Carlo simulations, the analysis evaluates the performance of these criteria in model selection under various data heterogeneity and distribution scenarios. The results demonstrate that the choice of information criterion can significantly influence both variable selection and the interpretation of model results. The study concludes that a comprehensive understanding of the behavior of information criteria can help researchers select and interpret the most appropriate models for hierarchical data, recommending a combined approach tailored to the specific objectives of the analysis.

Keywords:
designs; outliers; hierarchical data

1 INTRODUCTION

Structural Equation Models (SEM) are widely recognized as a powerful methodology for analyzing complex relationships among latent variables in various research fields. However, as research increasingly moves beyond independent data structures to encompass multilevel frameworks, new challenges emerge that demand more advanced analytical approaches (Davidov et al., 2012; Rappaport et al., 2019). In this context, the application of multilevel structural equation models (MSEM) has gained prominence, allowing for the analysis of hierarchical data frequently encountered in disciplines such as education, healthcare, psychology, and agriculture.

A fundamental aspect of MSEM implementation is the selection of estimation methods, as it directly impacts the accuracy of model fit evaluation, whether assessed through goodness-of-fit criteria (de L. Oliveira & Maçada, 2012; Dogan, 2022) or statistical information criteria (Xia & Yang, 2019). This process involves theoretical approaches specific to the method used for estimating model parameters. Notable estimation methods include Ordinary Least Squares (OLS) and Maximum Likelihood (ML), as well as more recent Bayesian approaches, such as Bayesian Least Squares Estimation (BLSEM) (Esmenda & Barrios, 2018; Zhang et al., 2023a; Fan & Wang, 2023).

Notably, the choice of estimation method exerts a profound influence on the interpretation of results and the precision of inferences drawn from the model. Contradictory outcomes between different criteria are frequently observed, underscoring the complexity of model selection (Shi & Maydeu-Olivares, 2020).

Zhang et al. (2023b) highlight that the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) were developed as general model selection methods, grounded in statistics and information theory. However, they emphasize that, in many cases, the selection of an appropriate model depends on specific data conditions that satisfy certain regularity assumptions. Accordingly, the authors suggest that if the goal of model selection is future prediction, adaptive information criteria combining AIC and BIC should be employed. On the other hand, if the primary goal is variable selection, BIC is recommended to avoid the inclusion of variables unlikely to prove significant.

In light of the critiques and divergences found in the literature, this study aims to conduct a computational analysis of the key information criteria: Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Browne-Cudeck Criterion (BCC), and Consistent Akaike Information Criterion (CAIC), applied to the selection of multilevel structural equation models. The choice of this class of models as the subject of study is justified by the fact that, when handling hierarchical data, the selection of variables to include in the model becomes more constrained, alongside the risk of incorrectly specifying fixed parameters, which can influence the results of information criteria.

The expectation of this study is that the insights gained will guide researchers in the proper use and interpretation of these criteria, ultimately aiding in the selection of the most appropriate model.

2 METHODOLOGY

In line with the proposed objectives, the statistical methodology employed in this study is outlined in the following sections:

2.1 Definition of Structural Equation Models

Multilevel models correspond to a hierarchical or nested structure, characterized by distinct levels in which data are organized. Each sublevel represents a linear model that captures the relationship between variables within a given level and specifies how the variables at that level influence relationships across other levels.

The structure of the correlation and covariance matrices, when considering the hierarchical model, takes the form of a block-diagonal arrangement, as illustrated in the matrices below, corresponding to a two-level model (1) and a three-level model (2).

$Co{r}_1=\left[\begin{array}{cccccc}1& {\rho }_1& {\rho }_2& 0& 0& 0\\ {\rho }_1& 1& {\rho }_2& 0& 0& 0\\ {\rho }_2& {\rho }_2& 1& 0& 0& 0\\ 0& 0& 0& 1& {\rho }_1& {\rho }_1\\ 0& 0& 0& {\rho }_1& 1& {\rho }_1\\ 0& 0& 0& {\rho }_1& {\rho }_1& 1\end{array}\right]$ (1)

$Co{r}_2=\left[\begin{array}{ccccccccc}1& {\rho }_1& {\rho }_2& {\rho }_2& 0& 0& 0& 0& 0\\ {\rho }_1& 1& {\rho }_2& {\rho }_2& 0& 0& 0& 0& 0\\ {\rho }_2& {\rho }_2& 1& {\rho }_1& 0& 0& 0& 0& 0\\ {\rho }_2& {\rho }_2& {\rho }_1& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& {\rho }_1& {\rho }_1& 0& 0\\ 0& 0& 0& 0& {\rho }_1& 1& {\rho }_1& 0& 0\\ 0& 0& 0& 0& {\rho }_1& {\rho }_1& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& {\rho }_1\\ 0& 0& 0& 0& 0& 0& 0& {\rho }_1& 1\end{array}\right]$ (2)

where ρ₁ represents the intra-unit correlation and ρ₂ the intra-individual correlation.

Denoting the variables measured at the levels as B (level 1) and W (level 2), the covariance matrix S ₁(3) is partitioned into the covariances S _BB referring to the covariance of variables at level 1 and S _BW and S _WB which correspond to the covariances between variables at level 1 and variables at level 2.

Similarly, the extension to three levels is represented by the partitioned covariance matrix S ₂(4). The covariances of the variables within each group are shown on the diagonal as S _BB (level 1), S _WW (level 2), and S _ZZ (level 3). The off-diagonal elements represent the covariances between the groups.

$S_1=\left[\begin{array}{cc}{S}_{BB}& S_{BF}\\ S_{FB}& S_{FF}\end{array}\right]$ (3)

$S_2=\left[\begin{array}{ccc}{S}_{BB}& S_{BF}& S_{BZ}\\ S_{BF}& S_{FF}& S_{FZ}\\ S_{ZB}& S_{ZF}& S_{ZZ}\end{array}\right]$ (4)

The reproduction of the covariance matrix adjusted by an MSEM (Figure 1) begins with the specification of the relationships between variables, as illustrated in Figure 1. This figure represents a hierarchical structural model with two levels, including a latent factor at level 1 (within) and a latent factor at level 2 (between), assuming k=3 variables.

Figure 1
Hierarchical structural model with two levels and three variables.

In the path diagram of Figure 1, the rectangles represent three indicators at Level 1, while the unidirectional arrows indicate causal effects. The ellipses η_Wij and η_Bj represent the latent variables at the within (intra) and between (inter) levels, respectively. The smaller arrows next to the rectangles refer to the error term within the level, ε_W ; the ellipses labeled y at the intermediate level correspond to the indicator variable at the intermediate level, and the smaller arrows near the indicators at the intermediate level represent the between-level error term ε_B , also referred to as the random term (Davidov et al., 2012).

The interpretation of the path diagram in Figure 2 follows similarly. The circles labeled y at the first level represent three indicators at Level 1, with unidirectional arrows denoting causal effects. The ellipses η_Bj , η_Wij , and η_Zhij represent the latent variables at Levels 1, 2, and 3, respectively. The rectangles labeled y at the intermediate level refer to the indicator variable at Level 2, while the diamonds labeled y represent the indicator variable at Level 3.

Figure 2
Multilevel structural equation model for three levels and three variables.

Based on these specifications (Figures 1 and 2), the estimated covariance matrices were computed according to the equation v=Av+u (McDonald & Hartmann, 1992), which expresses the linear regressions represented by arrows in the diagram, A refers to the regression coefficient matrices that are empirically determinedf for the two-level MSEM (5) and three-level MSEM (6), while v and u represent the vectors composed of observed and latent variables formed by the parametric values used in the Monte Carlo simulation.

$A_1=\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& -0.5& 0\\ 0& 0& 0& 0& 0& 0& -0.7& 0\\ 0& 0& 0& 0& 0& 0& -0.2& 0\\ 0& 0& 0& 0& 0& 0& 0& -0.3\\ 0& 0& 0& 0& 0& 0& 0& -0.4\\ 0& 0& 0& 0& 0& 0& 0& -0.6\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$ (5)

$A_2=\left[\begin{array}{cccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.3& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -0.4& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.8& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.7\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -0.5\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.2\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$ (6)

Thus, assuming a as the number of variables defined in v and k as the number of observed variables, in order for the covariance matrix imposed by the structural model (Figures 1 and 2) to be fully specified, the estimation of the structural model with hierarchical levels is given by (7):

${\Sigma }_{\theta k\times k}=J{\left(I_a-A\right)}^{-1}P{\left[{\left(I_a-A\right)}^{-1}\right]}^t{J}^t$ (7)

where P is the covariance matrix of the elements for two levels (8) and three levels (9) of the vector u, for two levels and three levels (10), assuming the parametric values used in the Monte Carlo simulation. J is the matrix whose first k columns form the identity matrix I _k .

$P_1=\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$ (8)

$P_2=\left[\begin{array}{cccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$ (9)

2.2 Specification of the Information Criteria Used in MSEM Selection

Initially, the information criteria and respective fit functions were defined according to each estimation method. Specifically, the following methods were analyzed: Maximum Likelihood (ML), Generalized Least Squares (GLS), and Unweighted Least Squares (ULS), as outlined in Table 1 (Bollen, 1989; Schermelleh-Engel et al., 2003; Thakkar, 2020).

Table 1
Specification of fit functions for Maximum Likelihood (F _ML ), Generalized Least Squares (F _GLS ), and Unweighted Least Squares (F _ULS ) methods.

Method	Fit Function
Maximum Likelihood (F ML )	$F_{ML}=\log \left|{\Sigma }_{\theta }\right|+tr\left(S\times \sum _{\theta }^{-1}\right)-\log \left|S\right|-\left(p+q\right)$, where the matrices are Σ and $\sum _{\theta }^{-1}$, in which |·| represent, respectively the determinant and the trace associated with these matrices, with dimensions p+q considered as observed variables.
Generalized Least Squares (F GLS )	$F_{GLS}=\frac{1}{2}tr{\left(\left[S-{\Sigma }_{\theta }\right]W^{-1}\right)}^2$ where $W_{k\times k}^{-1}$ is the weight matrix for the residual matrix, with the most common form being W -1=S -1.
Unweighted Least Squares (F ULS )	$F_{ULS}=\frac{1}{2}tr\left({\left[S-{\Sigma }_{\theta }\right]}^2\right)$

The statistics that define the information criteria were constructed following the expression that defines the minimum sample discrepancy <mml:math><mml:mover><mml:mi>C</mml:mi><mml:mo>^</mml:mo></mml:mover></mml:math>, as shown in Equation (10), for each discrepancy function described in Table 1.

${\hat{C}}_{\left(MV\right)}=N\times F_{\left(MV\right)};{\hat{C}}_{\left(GLS\right)}=N\times F_{\left(GLS\right)};{\hat{C}}_{\left(ULS\right)}=N\times F_{\left(ULS\right)}$ (10)

The information criteria listed in Table 2 were defined for each estimation method, following the discrepancy functions ${\hat{C}}_{\left(MV\right)}$, ${\hat{C}}_{\left(GLS\right)}$, and ${\hat{C}}_{\left(ULS\right)}$, in accordance with Equation (10).

Table 2
Specification of information criteria based on estimation methods: Maximum Likelihood $\left(\hat{C}={\hat{C}}_{\left(MV\right)}\right)$, Generalized Least Squares $\left(\hat{C}={\hat{C}}_{\left(GLS\right)}\right)$, and Unweighted Least Squares $\left(\hat{C}={\hat{C}}_{\left(ULS\right)}\right)$.

Criteria	Description	Author
$AIC=\hat{C}+2t$	t: the number of free parameters.	Akaike (1974)
$BIC=\hat{C}+t\times \ln \left(N_k\right)$	k: the number of observed variables.	Schwarz (1978)
	N: sample size
$BCC=\frac{\hat{C}+2t}{k\left({\displaystyle \frac{k+3}{N-k-2}}\right)}$	k: number of observed variables.	Browne & Cudeck (1992)
	N: sample size, b=N-1
$CAIC=\hat{C}+t\left(\ln \left(N\right)+1\right)$	t: the number of free parameters.	Akaike (1974)

It is recommended to select the model with the lowest AIC, BIC, BCC, and CAIC values. It is important to consider multiple information criteria simultaneously for a more comprehensive evaluation of model selection.

2.3 Evaluation Scenarios Used in the Monte Carlo Procedure

The differentiation concerning the dispersion of the sample covariance matrix (S) and the model-adjusted covariance matrix (Σ_θ) was expressed through the heterogeneity factor GH, obtained by the relative efficiency between the determinants (det) of the covariance matrices involved, as shown in Equation (11).

${\left(\frac{det\left|S\right|}{det\left|{\Sigma }_{\theta }\right|}\right)}^{\frac{1}{k}}=G_H;implyingthatS=\sqrt{G_H{\Sigma }_{\theta }}$ (11)

Thus, when considering low values for G _H , such as 2 and 4, we can observe a similarity in the dispersion of the covariance matrices between S and Σ_θ. In this case, the sample values will be more concentrated around the mean vector, and the information criteria will tend to not reject the proposed model. On the other hand, with higher values of G _H , such as 6 and 8, there is greater evidence that the criteria will tend to favor rejecting the model.

In the context of violating the assumption of multivariate normality, a simulation was conducted using multivariate normal distributions contaminated by outliers generated from symmetric distributions with excess kurtosis. In this context, the Student’s t and Lognormal distributions were considered, with mixture probabilities set at α=0.05, 0.10, and 0.15, as suggested by Cirillo & Barroso (2016) and Santos & Cirillo (2021).

The combination of multivariate distributions, mixture rates, and heterogeneity factors defined the scenarios used for reproducing the results. A script was developed in R software version 4.2.3 (R Core Team, 2023) to study the criteria (see Table 2), using 1,000 Monte Carlo simulations. In each simulation, the mean of the estimates for the proposed criteria was calculated.

3 RESULTS AND DISCUSSION

3.1 Performance of Information Criteria in Assessing Model Fit Quality, Considering Degrees of Heterogeneity for Two and Three Levels Under the Assumption of Multivariate Normality

The results presented in Table 3 correspond to the estimates of the information criteria AIC, BIC, BCC, and CAIC, under the assumption of multivariate normality. The results vary in terms of model fit quality, as indicated by the degree of heterogeneity (G _H ) between the sample covariance matrices and the covariances obtained from the model fit.

Table 3
Information Criteria for MSEM Estimated Using the ML, GLS, and ULS Methods, Obtained in the Simulation Study of Multivariate Normal Distribution with Different Degrees of Heterogeneity.

	G H	Two Levels			Three Levels
		ML	GLS	ULS	ML	GLS	ULS
AIC	2	344.00	230.50	519.00	266.00	144.11	291.00
	4	934.00	547.99	1,564.00	526.00	427.12	931.00
	6	1,414.00	1,073.14	2,619.00	966.00	561.14	1,931.00
	8	2,189.00	1,288.43	3,989.00	1,231.00	783.62	2,331.00
BIC	2	360.62	247.12	535.62	292.14	170.25	317.14
	4	950.62	564.61	1,580.00	552.14	453.28	957.14
	6	1,430.62	1,089.76	2,635.62	992.14	587.28	1,951.14
	8	2,205.62	1,305.05	4,005.62	1,257.14	809.76	2,357.14
BCC	2	344.00	230.50	519.00	266.01	144.11	291.01
	4	934.00	548.00	1,564.00	526.01	427.13	931.01
	6	1,414.00	1,073.15	2,619.00	966.01	561.15	1,931.01
	8	2,189.00	1,288.44	3,989.00	1,231.01	783.63	2,331.01
CAIC	2	245.43	230.50	534.03	288.55	166.66	313.55
	4	949.03	563.03	1,579.03	548.55	449.68	953.55
	6	1,429.03	1,088.18	2,634.03	988.55	583.69	1,953.55
	8	2,204.03	1,303.47	4,004.03	1,253.55	806.16	2,353.55

In the simulation of the ideal scenario, where the structural models for two and three levels (Figures 1 and 2) are fitted to data that meet the assumption of normality, it is important to note that lower values of the information criteria suggest a better model fit to the data. Accordingly, the criteria generated by model fitting using the Generalized Least Squares (GLS) method demonstrated superior performance across the simulated scenarios. This method produced the lowest values for its criteria compared to models fitted using the Maximum Likelihood (ML) and Unweighted Least Squares (ULS) methods, highlighting the effectiveness of GLS as an estimation method.

This study revealed notable discrepancies in the information criteria, contingent upon the estimation method employed for model fitting. The criteria values showed a marked increase, especially when comparing the fit achieved through the GLS method to that of the ULS method.

As anticipated, the information criteria values rose with increasing heterogeneity between the model and sample covariance matrices, reflecting a higher degree of model misfit. This indicates a deterioration in model fit, as shown by higher criterion values. It was noted that the increase in the criteria values generated by the model fit using the GLS method was not as pronounced as with the ML and ULS estimation methods. This suggests greater stability and robustness of the GLS method in relation to model specifications when compared to the ML and ULS methods.

Maximum Likelihood (ML) estimation in traditional structural equation modeling involves minimizing the deviations between the expected covariance matrix and the sample covariance matrix derived from observed data. In the context of multilevel data, the model’s log-likelihood criteria are adjusted at both levels. Specifically, this is described by equations for the aggregated -2 log likelihood across the levels of analysis (Liang & Bentler, 2004; Rappaport et al., 2019). These equations indicate that the overall model fit combines all levels of analysis, which can obscure the source of potential model misspecifications. Moreover, the overall model fit is more sensitive to misspecifications at level 1 than at the intermediate level (Hsu et al., 2015; Rappaport et al., 2019; Ryu & West, 2009; Yuan & Bentler, 2007). One approach to address this issue is to estimate the model fit separately at each level of a multilevel structural model (Rappaport et al., 2019; Ryu & West, 2009; Yuan & Bentler, 2007).

Unweighted Least Squares (ULS) is an option for estimating models with non-normal data. It is a consistent estimator (Bollen, 1989; Kyriazos & Poga-Kyriazou, 2023). However, ULS assumes that all observed variables are measured on the same scale (i.e., it is not scale-invariant). Despite this limitation, ULS remains a viable option when the data are not normally distributed (Kline, 2016; Lei & Wu, 2012).

ULS is a commonly used estimation method in traditional Structural Equation Modeling (SEM), but its use in the context of Multilevel Structural Equation Modeling (MSEM) is less frequent compared to other methods, such as Maximum Likelihood (ML) and Generalized Least Squares (GLS). ULS does not account for the multilevel structure of the data and assumes homoscedasticity (Hox et al., 2017), which may have contributed to the poorer performance of the simulated criteria.

One of the advantages of ULS is its computational simplicity and speed in estimation compared to ML and GLS, which can be useful for large multilevel datasets. However, its use should be exercised with caution to avoid misleading conclusions.

Generalized Least Squares (GLS) is an alternative to ML, assuming normally distributed and continuous data, or data with mild non-normality and non-extreme kurtosis. GLS is less computationally demanding than ML, while still producing comparable model fit quality, especially with large samples (Browne, 1974; Brown, 2015; Yuan & Chan, 2016).

GLS estimates are considered consistent, unbiased, and asymptotically normally distributed, meaning that the distribution of its parameter estimates approaches normality as the sample size increases (Wang & Wang, 2020). Additionally, GLS is efficient (Lei & Wu, 2012). It allows for the consideration of complex covariance structures between variables and across hierarchical levels, resulting in more precise parameter estimates for the model. GLS offers a more flexible, robust, and efficient approach for modeling multilevel data, particularly in cases where the assumptions of traditional methods may not be adequately met (Leeuw & Meijer, 2008).

Studies on the impact of estimation methods on SEM information criteria, such as RMSEA, CFI, and SRMR, have been conducted and yielded similar results. Given the same type and level of model misspecification, the choice of estimation method had a significant effect on the population values of RMSEA and CFI. Specifically, for RMSEA, population values based on ML may be higher or lower than those calculated using least squares methods (i.e., unweighted least squares and diagonally weighted least squares - DWLS) (Shi & Maydeu-Olivares, 2020; Xia & Yang, 2019). According to (Xia & Yang, 2018, 2019; Shi & Maydeu-Olivares, 2020), with respect to CFI, the CFIML values were generally lower than those obtained using ULS or DWLS, suggesting a worse model fit.

3.2 Performance of Information Criteria in Evaluating Model Fit Quality, Considering Degrees of Heterogeneity for Two and Three Levels with Violation of the Multivariate Normality Assumption

Table 4 and Figures 3 and 4 present the influence of outlier observations on the information criteria, considering the degree of heterogeneity (G _H ) between the sample and structural covariance matrices, as well as the percentage of outlier observations (α) generated by the multivariate Student’s t distribution for two and three levels of MSEM models.

Figure 3
Simulated BIC Values for Two- and Three-Level Models, with Outliers Generated by the Multivariate t-Student Distribution, Varying the Mixture Rate α and the Degree of Heterogeneity.

Figure 4
Simulated BCC Values for Two- and Three-Level Models, with Outliers Generated by the Multivariate t-Student Distribution, Varying the Mixture Rate α and the Degree of Heterogeneity.

Table 4
Information Criteria for MSEM Estimated Using the ML, GLS, and ULS Methods, Obtained in the Simulation Study Considering Outliers Generated by the Multivariate t-Student Distribution with Mixture Rate α and Different Degrees of Heterogeneity.

TWO LEVELS
G H	α	AIC ML	CAIC ML	AIC GLS	CAIC GLS	AIC ULS	CAIC ULS
2	0.05	274.00	289.03	179.86	194.89	389.00	404.03
4		849.00	864.03	657.62	672.65	1,429.00	1,444.03
6		1,434.00	1,449.03	1,052.74	1,067.77	2,654.00	2,669.03
8		2,044.00	2,059.03	1,196.70	1,211.73	3,824.00	3,839.03
2	0.10	214.00	229.03	145.13	160.17	304.00	319.03
4		974.00	989.03	648.89	663.93	1,599.00	1,654.03
6		1,544.00	1,559.03	1,001.19	1,016.23	2,764.00	2,779.03
8		2,189.00	2,204.03	1,200.57	1,215.60	4,184.00	4,199.03
2	0.15	259.00	274.03	169.79	184.83	374.00	389.03
4		929.00	944.03	636.66	651.69	1,639.00	1,654.03
6		1,474.00	1,489.03	1985.29	1,000.33	2,649.00	2,664.03
8		2,374.00	2,389.03	1,456.13	1,471.16	4,604.01	4,619.03
THREE LEVELS
2	0.05	261.00	283.55	201.07	223.62	276.00	298.55
4		571.00	593.55	335.22	357.78	996.00	1,018.55
6		1,176.00	1,198.55	739.22	761.77	2,091.00	2,113.55
8		1,336.00	1,358.55	772.15	794.71	2,686.00	2,708.55
2	0.10	216.00	283.55	216.93	239.48	331.00	353.55
4		646.00	668.55	431.58	454.13	1,2760.00	1,298.55
6		1,221.00	1,243.55	972.05	994.60	2,251.00	2,273.55
8		1,326.00	1,348.00	783.85	806.39	2,591.00	2,613.55
2	0.15	271.00	293.55	150.60	172.71	356.00	378.55
4		691.00	713.55	456.11	478.66	1,111.00	1,133.55
6		1,166.00	1,188.55	804.06	826.61	2,386.00	2,408.55
8		1,421.00	1,443.55	1,081.17	1,103.72	2,771.00	2,793.55

Based on the results, it was generally observed that, as expected, the lowest values for the criteria were achieved when a low degree of heterogeneity and low outlier rates (α=5%) were adopted, regardless of the estimatio(G _H =2)n method used. This implies a better model fit.

The worst outcomes, reflected in the highest information criterion values, occurred with a high degree of heterogeneity (G _H =8) and elevated outlier rates (α=15%). These results underscore the importance of carefully managing data quality in high-heterogeneity environments to avoid significant model misfit.

Intermediate degrees of heterogeneity G _H =4 and 6 did not show a significant increase in the criteria values when the outlier contamination rate increased from 10% to 15%. In some instances, a slight decrease in criterion values was observed, particularly in the two-level. In some instances, a slight decrease in criterion values was observed, particularly in the two-level structural model, potentially indicating that smaller models are less sensitive to moderate increases in outlier concentration. This suggests the importance of refining the criteria to incorporate the effect of levels, as the disturbances violating the normality assumption caused by the factors evaluated in this study affected the criteria estimates depending on the levels adopted in the model, which may influence the interpretation of the structural model’s fit quality.

Regarding the best estimation method used to fit the models, results similar to those from the multivariate normal distribution cases were found. The GLS method stood out, while ULS produced the worst results, showing a rapid increase in the values of the information criteria as the degree of heterogeneity increased.

Although ULS offers computational simplicity and robustness for non-normal data (Mindrila, 2010), it may not be the best choice in multilevel modeling situations due to its inability to directly handle the hierarchical structure of the data and its inadequacy in modeling multilevel variability. Therefore, it is crucial to carefully consider the characteristics of your data and the objectives of your analysis when selecting the estimation method in the context of MSEM.

According to Jung & Takane (2008), ULS frequently generates imprecise solutions, such as negative estimates or limitations in unique variances, due to its handling of measurement errors in observed variables. Unique variance represents disturbance in SEM, reflecting random error due to unreliability or measurement error. It is the portion of the total variance of a latent or observed variable that is not explained by the independent variables included in the model. In MSEM, unique variance can be viewed as the portion of a variable’s variance not attributed to higher levels of the hierarchy. In multilevel models, it is crucial to distinguish between unique variance within each level of the hierarchy and variance shared between levels. Additionally, the use of ULS tends to introduce more bias and reduce the accuracy of parameter estimates (Forero et al., 2009).

Yuan & Zhong (2013) state that outliers tend to provide a seemingly good fit to the model. However, many criteria based on maximum likelihood estimation indicate that the model is poor. Yuan & Bentler (2001) highlight a key deficiency: even if the model is correctly specified, a small proportion of outliers can distort the covariance matrix estimates imposed by the structural model, undermining the ability to make reliable assertions about the model’s fit quality.

Considering the outliers generated by the Lognormal distribution, characterized by skewness and excess kurtosis, it was observed, through Table 5 and Figures 5 and 6, that once again the criteria generated by model fitting using the GLS estimation method stood out. For this method, with a high degree of heterogeneity (G _H =8) and outlier concentrations of 10% and 15%, the criteria remained similar when considering the two-level structural model.

Figure 5
Simulated BIC Values for Two- and Three-Level Models, with Outliers Generated by the Multivariate Lognormal Distribution, Varying the Mixture Rate α and the Degree of Heterogeneity.

Figure 6
Simulated BCC Values for Two- and Three-Level Models, with Outliers Generated by the Multivariate Lognormal Distribution, Varying the Mixture Rate α and the Degree of Heterogeneity.

Table 5
Information Criteria for MSEM Estimated Using the ML, GLS, and ULS Methods, Obtained in the Simulation Study Considering Outliers Generated by the Multivariate t-Student Distribution with Mixture Rate α and Different Degrees of Heterogeneity.

TWO LEVELS
G H	α	AIC ML	CAIC ML	AIC GLS	CAIC GLS	AIC ULS	CAIC ULS
2	0.05	274.00	289.03	143.23	158.26	399.00	414.03
4		754.00	769.03	591.91	606.95	1,354.00	1,369.03
6		1,819.00	1,834.03	1,238.53	1,253.57	3,399.00	3,414.00
8		1,929.00	1,944.03	1,202.67	1,217.71	3,814.00	3,829.03
2	0.10	189.00	204.03	101.42	116.46	254.00	269.03
4		999.00	1,014.03	654.30	669.33	1,764.00	1,779.03
6		1,524.00	1,539.03	953.81	968.84	2,939.00	2,954.03
8		2,174.00	2,189.03	1,202.82	1,217.86	3,884.00	3,899.03
2	0.15	274.00	289.03	215.25	230.28	384.00	399.03
4		994.00	1,009.03	643.85	658.88	1,769.00	1,784.03
6		1,589.00	1,604.03	997.39	1,012.42	2,979.00	2,994.03
8		2,224.00	2,239.03	1,433.38	1,448.42	4,224.00	4,239.03
THREE LEVELS
2	0.05	271.00	293.55	191.18	213.73	421.00	443.55
4		726.00	748.55	474.22	497.48	1,291.00	1,313.55
6		911.00	933.55	563.69	586.24	1,746.00	1,768.55
8		1,466.00	1,488.55	926.23	948.78	3,136.00	3,158.55
2	0.10	206.00	228.55	103.68	126.23	421.00	443.55
4		661.00	683.55	556.02	578.57	1,091.00	1,113.55
6		991.00	1,013.55	480,70	503.25	1,996.00	2,018.55
8		1,301.00	1,323.00	908.72	931.27	2,621.00	2,643.55
2	0.15	296.00	318.55	160.95	183.50	441.00	463.55
4		766.00	788.55	428.82	451.37	1,391.00	1,413.55
6		1,036.00	1,058.55	745.14	767.69	1,846.00	1,868.55
8		1,531.00	1,553.55	990.08	1,012.64	2,901.00	2,923.55

Despite showing the best results, according to Savalei & Bentler (2006) and Brown (2015), GLS is appropriate only when variables do not exhibit excess kurtosis. Therefore, the results for GLS may not be reliable in cases of high kurtosis.

In this specific case of the Lognormal distribution, the criteria generated by model fitting using the ML method should be considered as the best-fitting case. This is because the values found are quite robust to non-normality and proved to be less biased than those obtained through model fitting using other estimation methods.

Especially in cases with a low level of outlier contamination (α=5%), Maximum Likelihood (ML) can also be applied to variables that deviate slightly from normality (Bollen, 1989; Kyriazos & Poga-Kyriazou, 2023), particularly when the primary focus is on parameter estimates (Raykov & Marcoulides, 2006). However, using ML with extremely non-normal data could lead to an artificial increase in χ² and an excessive rejection of the mo Lower degrees of heterogeneity (G _H =2 and 4) produced similar criterion values across both low and high outlier concentrations (α=5% and 15%), irrespective of the estimation method. This may indicate that the criteria are inconclusive in determining whether the model is well-fitted.

For the three-level model, the best results were observed with an intermediate contamination rate (α=10%), with the exception of the case with G _H =6.

When the data are not normally distributed, the estimates of factor loadings can be biased or less efficient, which may lead to inaccurate conclusions about the relationships between latent and observed variables. As a result, a substantial body of literature has advocated for the use of regularization methods to address this issue (Jacobucci et al., 2016; Yuan & Bentler, 2017; Yuan & Chan, 2016; Zulkifli et al., 2023).

4 CONCLUSIONS

The comparative study of information criteria applied to multilevel structural equation models (MSEM) revealed valuable insights into the performance of different estimation methods across scenarios with varying heterogeneity and violations of normality.

The GLS method consistently outperformed other estimation methods across most scenarios, particularly when the data met the assumption of multivariate normality. This superior performance is likely due to the method’s ability to accommodate complex covariance structures while maintaining computational efficiency.

In contrast, the Unweighted Least Squares (ULS) method showed the poorest performance, as evidenced by significantly higher information criteria, particularly in high-heterogeneity scenarios. This finding suggests that ULS is ill-suited for multilevel models with substantial heterogeneity and should be applied with caution when data deviate significantly from normality, despite its computational simplicity. As it may result in less accurate and potentially biased estimates.

In summary, the combined application of AIC, BIC, BCC, and CAIC criteria provides a robust framework for guiding researchers in selecting and interpreting the most appropriate models for hierarchical data. Future analyses should consider the specific characteristics of the dataset, such as heterogeneity and outlier concentration, to ensure accurate and meaningful conclusions.

References

• AKAIKE H. 1974. A new look at statistical model identification. IEEE Transactions on Automatic Control, 19(6): 716-723.
• BOLLEN KA. 1989. Structural Equations with Latent Variables. New York: John Wiley & Sons.
• BROWN TA. 2015. Confirmatory Factor Analysis for Applied Research. 2nd ed. New York, NY: Guilford.
• BROWNE MW. 1974. Generalized least squares estimators in the analysis of covariance structures. South African Statistical Journal, 8(1): 1-24.
• BROWNE MW & CUDECK R. 1992. Alternative ways of assessing model fit. Sociological Methods & Research, 21(2): 230-258.
• CIRILLO MA & BARROSO LP. 2016. Effect of outliers on the GFI quality adjustment index in structural equation model and proposal of alternative indices. Communications in Statistics - Simulation and Computation, 46(3): 1895-1905.
• DAVIDOV E, DÜLMER H, SCHLÜTER E, SCHMIDT P & MEULEMAN B. 2012. Using a multilevel structural equation modeling approach to explain cross-cultural measurement noninvariance. Journal of Cross-Cultural Psychology, 43(4): 558-575.
• DE L OLIVEIRA D & MAÇADA ACG. 2012. Valor da TI na perspectiva das capacidades internas: uma análise de desempenho multinível. Revista de Administração e Negócios da Amazônia, 4(3): 140-155.
• DOGAN I. 2022. A simulation study comparing model fit measures of structural equation modeling with multivariate contaminated normal distribution. Communications in Statistics - Simulation and Computation, 51(5): 2526-2536.
• ESMENDA MJ & BARRIOS EB. 2018. Robust estimation of a multilevel model with structural change. Communications in Statistics - Simulation and Computation, 47(4): 1014-1027.
• FAN TH & WANG YF. 2023. Bayesian model selection for structural equation models for myopia data. Communications in Statistics - Simulation and Computation, 52(11): 5680-5694.
• FORERO CG, MAYDEU-OLIVARES A & GALLARDO-PUJOL D. 2009. Factor analysis with ordinal indicators: A Monte Carlo study comparing DWLS and ULS estimation. Structural Equation Modeling, 16(4): 625-641.
• HOX J, MOERBEEK M & VAN DE SCHOOT R. 2017. Multilevel Analysis: Techniques and Applications. 3rd ed. New York, NY: Routledge.
• HSU HY, KWOK O, LIN JH & ACOSTA S. 2015. Detecting misspecified multilevel structural equation models with common fit indices: A Monte Carlo study. Multivariate Behavioral Research, 50: 197-215.
• JACOBUCCI R, GRIMM KJ & MCARDLE JJ. 2016. Regularized structural equation modeling. Structural Equation Modeling, 23(4): 555-566.
• JUNG S & TAKANE Y. 2008. Regularized common factor analysis. In: New Trends in Psychometrics. pp. 141-149.
• KLINE RB. 2016. Principles and Practice of Structural Equation Modeling. 4th ed. New York, NY: Guilford Press.
• KYRIAZOS T & POGA-KYRIAZOU M. 2023. Applied psychometrics: Estimator considerations in commonly encountered conditions in CFA, SEM, and EFA practice. Psychology, 14: 799-828.
• LEEUW JD & MEIJER E (Eds.). 2008. Handbook of Multilevel Analysis. Springer.
• LEI PW & WU Q. 2012. Estimation in structural equation modeling. In: HOYLE RH (Ed.), Handbook of Structural Equation Modeling. pp. 164-180. Guilford.
• LIANG J & BENTLER PM. 2004. An EM algorithm for fitting two-level structural equation models. Psychometrika, 69: 101-122.
• MCDONALD RP & HARTMANN WM. 1992. A procedure for obtaining initial values of parameters in the RAM model. Multivariate Behavioral Research, 27(1): 57-76.
• MINDRILA D. 2010. Maximum likelihood (ML) and diagonally weighted least squares (DWLS) estimation procedures: A comparison of estimation bias with ordinal and multivariate non-normal data. International Journal of Digital Society, 1(1): 60-66.
• R CORE TEAM. 2023. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. Available at: https://www.R-project.org/
  » https://www.R-project.org/
• RAPPAPORT LM, AMSTADTER AB & NEALE MC. 2019. Model fit estimation for multilevel structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 0: 1-12.
• RAYKOV T & MARCOULIDES GA. 2006. A First Course in Structural Equation Modeling. Mahwah, NJ: Erlbaum.
• RYU E & WEST SG. 2009. Level-specific evaluation of model fit in multilevel structural equation modeling. Structural Equation Modeling, 16: 583-601.
• SANTOS PM & CIRILLO MA. 2021. Construction of the average variance extracted index for construct validation in structural equation models with adaptive regressions. Communications in Statistics-Simulation and Computation, 50(4): 1639-1650.
• SAVALEI V & BENTLER PM. 2006. Structural equation modeling. In: The Corsini Encyclopedia of Psychology, vol. 330. p. 36.
• SCHERMELLEH-ENGEL K, MOOSBRUGGER H & MÜLLER H. 2003. Evaluating the fit of structural equation models: Tests of significance and descriptive goodness-of-fit measures. Methods of Psychological Research, 8(2): 23-74.
• SCHWARZ G. 1978. Estimating the dimension of a model. Annals of Statistics, 6: 461-464.
• SHI D & MAYDEU-OLIVARES A. 2020. The effect of estimation methods on SEM fit indices. Educational and Psychological Measurement, 80(3): 421-445.
• THAKKAR JJ. 2020. Structural Equation Modelling: Application for Research and Practice.
• WANG J & WANG X. 2020. Structural Equation Modeling. 2nd ed. Wiley.
• XIA Y & YANG Y. 2018. The influence of number of categories and threshold ordered categorical data. Multivariate Behavioral Research, 53(5): 731-755.
• XIA Y & YANG Y. 2019. RMSEA, CFI, and TLI in structural equation modeling with ordered categorical data: The story they tell depends on the estimation methods. Behavior Research Methods, 51(1): 409-428.
• YUAN KH & BENTLER PM. 2001. Effects of ’outliers’ on estimators and tests in covariance structure analysis. British Journal of Mathematical and Statistical Psychology, 54(1): 161-175.
• YUAN KH & BENTLER PM. 2007. Multilevel covariance structure analysis by fitting multiple single-level models. Sociological Methodology, 37: 53-82.
• YUAN KH & BENTLER PM. 2017. Improving the convergence rate and speed of Fisher-scoring algorithm: Ridge and antiridge methods in structural equation modeling. Annals of the Institute of Statistical Mathematics, 69: 1-27.
• YUAN KH & CHAN W. 2016. Structural equation modeling with unknown population distributions: Ridge generalized least squares. Structural Equation Modeling: A Multidisciplinary Journal, 23(2): 163-179.
• YUAN KH & ZHONG X. 2013. Robustness of fit indices to outliers and leverage observations in structural equation modeling. Psychological Methods, 18(2): 121-136.
• ZHANG J, LU J, XU X & TAO J. 2023a. Bayesian multilevel multidimensional item response modeling approach for multiple latent variables in a hierarchical structure. Communications in Statistics - Simulation and Computation, 52(7): 2822-2842.
• ZHANG J, YANG Y & DING J. 2023b. Information criteria for model selection. Wiley Interdisciplinary Reviews: Computational Statistics, 15(5): e1607.
• ZULKIFLI NR, AIMRAN N & DENI SM. 2023. The performance of unweighted least squares and regularized unweighted least squares in estimating factor loadings in structural equation modeling. International Journal of Data and Network Science, 7(3): 1017-1024.

• Data Availability
  All data used in this study were generated through simulation, as described in the methodology section. The code used for data generation is available upon request from the authors.
• Funding
  The authors would like to thank the Research Supporting Foundation of Minas Gerais State (Fundação de Amparo à Pesquisa do Estado de Minas Gerais - FAPEMIG), for the financial support through research in the project n° BPD-00448/22 and National Council for Scientific and Technological Development (Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq) for the financial support through research in project 304939/2021-8.

Edited by

• Editor responsible for the review
  Editor-in-Chief: Annibal Parracho Sant’Anna.

Data availability

All data used in this study were generated through simulation, as described in the methodology section. The code used for data generation is available upon request from the authors.

Publication Dates

• Publication in this collection
  08 Aug 2025
• Date of issue
  2025

History

• Received
  18 Mar 2025
• Accepted
  20 June 2025