Evaluation of Model Fit in Structural Equation Models with Ordinal Missing Data: A Comparison of the D₂ and MI2S Methods

ABSTRACT

Social science research often utilizes measurement instruments that generate ordinal data (e.g., Likert scales). Many empirical studies also face the challenge of missing data, which can be addressed by performing multiple imputation followed by analyses of the imputed datasets. However, when missing data exist on ordinal variables, there has been limited research on how to evaluate model fit of structural equation models for ordinal variables. Recent studies suggest that two multiple-imputation-based approaches show great promise: The D₂ procedure, and the Multiple Imputation Two-step (MI2S) approach, though the two have not been systematically compared. This study extends previous research by comparing the D₂ with the MI2S fit statistics in a wider range of conditions than previous studies. Our findings revealed a number of factors that can influence the performance of these test statistics.

KEYWORDS:

• missing data
• model fit
• multiple imputation
• ordinal data

Acknowledgements

This work was completed in part with resources provided by the Research Computing Data Core at the University of Houston.

Supplemental Material

Supplemental data for this article can be accessed on the publisher’s website.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Notes

¹ Following Millsap and Tein (2004) and Millsap (2011), standardized threshold parameters refer to the threshold parameters for a continuous latent response variable that follows a standard normal distribution.

² Note that the MI2S approach applied to ordinal data does not utilize the standard errors or the mean- or mean-and-variance-corrected fit statistics produced by the SEM software packages to compute the MI2S fit statistics ${\tilde{T}}_B$ and ${\tilde{T}}_{YB}$. Thus, using the different variants of the ULS estimator (“ULS”, “ULSM”, “ULSMV”, or “ULSMVS”) in lavaan in Step 2 of the MI2S approach does not influence the calculated MI2S fit statistics ${\tilde{T}}_B$ and ${\tilde{T}}_{YB}$.

³ Liu and Sriutaisuk (2020) investigated the performance of the D₂ method to combine mean- or mean-and-variance-corrected model fit statistics of the variants of ULS and DWLS estimators (ULSM, WLSM, ULSMV, and WLSMV), and found that using the D₂ procedure with the robust test statistics produced by ULSMV or WLSMV had similar performance in their simulations, and generated D₂ statistics with better calibrated empirical means and variances than those produced by ULSM or WLSM.

⁴ To do this, we utilized the calculated D2() function in the R package semTools (Jorgensen et al., 2019).

⁵ Within each imputed dataset, we obtained the thresholds and polychoric correlations using the lavCor() function in lavaan, and used vcov() in lavaan to obtain the asymptotic covariance matrix of these thresholds and polychoric correlations. The thresholds, polychoric correlations, and the asymptotic covariance matrices were further combined in R.

⁶ Issues occurred in the form of model nonconvergence, negative estimated residual variances of the latent response variables underlying ordinal items, and none positive-definite covariance matrix of the latent factors.

⁷ Issues occurred in the form of model nonconvergence, non-invertible Information matrix, negative estimated residual variances of the latent response variables underlying ordinal items, and none positive-definite covariance matrix of the latent factors.

⁸ Except for the threshold for the C = 2 conditions, where the population value was 0. For the threshold in the C = 2 conditions we instead computed the standardized bias, which was always below 40%, the cutoff for practically significant bias suggested in Collins et al. (2001).

⁹ As can be seen from the results described below, as the proportion of incomplete analysis variables increased from 33% to 67%, while the model size remains the same (three-factor CFA for 15 items), the performance of the D₂ and MI2S test statistics worsened. However, situations with 50% incomplete analysis variables involved a change in model size in addition to a change in the proportion of incomplete analysis variables. We conducted additional analyses to manipulate the proportion of incomplete items, without the confounding effect of model size. Specifically, in conditions with incomplete M and Y items, we fitted a Correct Model 3 (a two-factor CFA model with one factor for the five M items, and another factor for the five Y items). Correct Model 3 has 10 analysis variables and 100% incomplete analysis items, and can be directly compared to Correct Model 2 (a two-factor CFA model for the complete X and the incomplete Y items) which also has 10 analysis variables but only 50% incomplete items in conditions with incomplete M and Y items. We summarized in Supplemental Material 2 the comparison of the performance of the test statistics for Correct Models 2 versus 3 in conditions with incomplete M and Y items. As can be seen from Supplemental Material 2, Correct Model 3 (two-factor CFA for 10 items with 100% incomplete analysis variables) is always associated with a greater risk of and a greater extent of Type 1 error inflation of D₂ and ${\tilde{T}}_{YB}$, compared to Correct Model 2 (two-factor CFA for 10 items with 50% incomplete analysis variables). Thus, we found the same common trend in comparing situations with 33% versus 67% incomplete analysis variables (both using three-factor CFA for 15 items), and in comparing situations with 50% versus 100% incomplete variables (both using two-factor CFA for 10 items). That is, the proportion of incomplete analysis variables in the analysis model has a strong impact on the performance of the D₂ and MI2S test statistics.

¹⁰ In addition to a high-low loading configuration, these three conditions with 50% missing data on 33% of the analysis variables in which D₂ showed inflated Type 1 error rates also had a smaller sample size of 250 or 500, and fewer response categories (C = 2 or 3). One of these three conditions (C = 2, N = 250, high-low loading configuration, 50% missingness on Y items only) had a maximum relative bias of parameter estimates from the standard multiple imputation approach (for the calculation of D₂) that was over 10% (13.4%).

¹¹ This N = 250 condition (C = 2, N = 250, high-low loading configuration, 50% missingness on Y items only) had a maximum relative bias of parameter estimates from the MI2S approach that was close to 10% (9.6%).

¹² At N = 250 with 50% missing data on the incomplete variables, conditions with incomplete M and Y items in general had greater challenges in getting the MCMC chains for multiple imputation to converge, and in getting the analysis models fitted to the imputed datasets to converge to a proper solution, compared to conditions with incomplete Y items only. This likely caused the slight differences in the performance of the test statistics in Correct Model 2 of X and Y items at N = 250, across data generation conditions with incomplete M and Y items and conditions with incomplete Y items only.

¹³ The statistical power information, along with the numbers of usable replications and the average usable m in retained replications, can be found in Supplemental Material 5.

¹⁴ We also conducted a small-scale additional simulation to explore if the same pattern of results would be observed if the data were missing completely at random (MCAR). The results are summarized in Supplemental Material 6. The patterns of results under MCAR are consistent with the patterns of results in corresponding conditions under MAR reported in the text, suggesting that our findings based on MAR should hold when the missing data mechanism is MCAR.