Correcting the Bias of the Root Mean Squared Error of Approximation Under Missing Data

Summary of This Article

In this article, we will first discuss the problem that missing data introduces into the realm of model fit indices in SEM, focusing on the RMSEA. We then introduce the mathematical basis for a bias correction of the RMESA. Next, we present data simulations and plots to illustrate the problem and demonstrate to what extent the proposed solutions are better suited to capture model misfit under missing data than the uncorrected RMSEA. We discuss our motivation for a correction, broader implications of the problem of missing data in SEM, limitations of our correction, and future directions in the research of model misfit indices.

Mathematical Statement of the Problem

The RMSEA is an estimate based on the quantity $\sqrt{\frac{F}{df}}$, where $\mathit{df}$ represents the difference of degrees of freedom between a fitted model and a saturated model, and $F$ is a population quantity of the discrepancy between a normal population distribution and the best-fitting normal distribution under a given model. $F$ is thus independent of sample size and, conventionally, $\frac{{\chi }^2}{N}$ is used as estimator of this discrepancy. The resulting population value, $\sqrt{\frac{F}{df}}$, is thus a measure of absolute misfit per person and per degree of freedom. The fraction $\sqrt{\frac{{\chi }^2}{N\cdot df}}$ is an estimate of this population quantity for model misfit, which forms the basis of RMSEA.

For complete data, the RMSEA for a proposed model is defined as

where $\mathit{df}$ is the difference in degrees of freedom between the proposed model and a saturated normal model, $N$ the number of participants, and ${\chi }^2$ the difference in the minus two log-likelihood between the proposed and saturated models for the given data set. Without model misspecification, the ${\chi }^2$ index will be distributed as a central ${\chi }^2$ distribution with $\mathit{df}$ degrees of freedom; therefore, the fraction will have an expected value of zero, with decreasing sampling variance as $N$ increases. If the model is misspecified, ${\chi }^2$ follows a non-central ${\chi }^2$ distribution with a non-centrality parameter proportional to $N$, and thus RMSEA will converge towards a fixed non-zero number in expectation. The degree to which RMSEA is larger than zero is therefore a measure of model misfit.

As a data set has an increasing proportion of missing data (e.g., assuming missing completely at random; MCAR), the minus two log likelihoods of both the proposed and saturated models will in expectation decrease proportionally to that amount. In consequence, the ${\chi }^2$ – the difference between the two – will also decrease proportionally to the proportion of missing data, while $\mathit{df}$ and $N$ remain constant. Hence, RMSEA values will become smaller in value and, hence, more optimistic for increasing proportions of missing data. Note that when there is no model misspecification, the proportion of missing data does not bias RMSEA.

Defining the Correction

How can we best correct the computation of RMSEA such that it is invariant under the proportion of missing data for misspecified models? We require that the bias-correction must leave RMSEA unchanged if there are no missing data. Second, we require the bias-correction to yield asymptotically identical RMSEA values under misspecification and different levels of missing data under the assumption of MCAR and MAR. In the following, we propose a method to correct the bias incurred in the RMSEA by missing data.

To account for the downward bias in the ${\chi }^2$ as estimator of model misfit, the ${\chi }^2$ value can be replaced by an estimate of its expected value if no missing data were present. We obtain this bias correction by computing the expected divergence of a single person with complete observations and the best-fitting model and multiply this by the effective sample size, that is, the sample size multiplied by the proportion of observed values, to obtain the expected discrepancy if there had been no missing data.

Under MCAR or MAR missingness, the model-implied covariance matrix from a model with missing data is an estimate of the population covariance matrix, and by extension, of the model-implied covariance matrix if no missing data were present. Using the Kullback-Leibler divergence (Kullback & Leibler, 1951), also known as relative entropy, between the model-implied distribution and the saturated model’s distribution, we compute the expected value of the ${\chi }^2$ value of a single, fully observed person. Hence, the negative two log-likelihood of the complete data can be estimated as $2N$ times the Kullback–Leibler (Kullback & Leibler, 1951) divergence between the two distributions, just as $2N$ is multiplied by the fit function in covariance modeling to generate ${\chi }^2$ values. $2KL$ is given by

where $\text{Σ}$ and $\mu$ are the model distribution covariance matrices and mean, $S$ and $m$ the saturated model distribution covariance matrix and mean, and $k$ the number of variables (i.e., the dimensionality of the distributions).

Note that $2\cdot KL$ is identical to $\frac{{\chi }^2}{N}$ for complete data sets and, thus, the bias correction will not change RMSEA for complete data sets. To ensure an unbiased estimation of the average discrepancy for missingness with rate $1-p$, we need to discount the sample size, $N$, by the proportion of observed cases, $p$. The resulting term $2pN\cdot KL$ has the same expectation as ${\chi }^2$ for all rates of missingness if $S=\text{Σ}$, that is, if the model is correctly specified. So we obtain a bias correction if we replace our estimator of $F$ with $2KL$. For $S\ne \text{Σ}$ in the population, $2pN\cdot KL-\mathit{df}$ has a non-zero expectation, that is, asymptotically proportional to $pN$. Adding the maximum operator to avoid negative values under the square root, we obtained the bias-corrected RMSEA:

Note that even though the fraction term is guaranteed to have a constant expectation for different missingness rates, the expression under the square root (which includes the $max$ operator) may show an increasing expectation with $p$ if the distribution of the fraction term includes negative results. However, the actual center of the distribution is not changed; the effect is just due to the negative values being moved to zero. The median and mode of the distribution stay constant.

Simulation Design

To demonstrate the effect of the proposed bias correction, we performed two simulations. The first uses a small SEM to demonstrate the effect best without interference, while the second uses a larger Latent Growth Curve Model taken from a real substantive study.

In the first simulation, the generating model was a bivariate normal distribution with zero mean, variance one, and a correlation which we varied in three steps, no correlation (r = 0.0), a typical weak association one may encounter in psychological research (r = 0.125), and a very strong correlation (r = 0.9) to investigate the behavior close to the limit. N = 1000 participants were simulated. We then analyzed the data with a model in which we fixed the covariance to be zero, using a saturated mean structure (this is akin to a covariance-only model). Note that this miss-specified model is identical to the null model if one would test against a zero covariance between the two manifest variables.

For a larger substantive example, we simulated misfit in a linear Latent Growth Curve Model (LGCM; McArdle, 1988). The linear LGCM models change across repeated observations by specifying an intercept and a slope component describing the overall level and change as well as individual differences in each. Our simulation was inspired by estimates that are typical for cognitive ageing research (Ghisletta et al., 2020). The model contained five manifest variables. We set the intercept mean to 20 and the intercept variance to 30, the slope mean to –4 and the slope variance to 5 and the residual error to 14. This corresponds to a moderate growth curve reliability (Figure 1). We used the same model for data generation, but to add some typical misspecification, we added a quadratic growth component with mean = –5 and variance = 5. $N=9000$ participants were generated with this model.

In both simulations, we then created missing data by removing all values except for the first variable in some participants. The choice of participants with missing data was done either by an MCAR or MAR process: For the MCAR process, each participant had the same probability $p_{miss}$ to miss all values but the first. For the MAR process, the probability $p_{miss;i}$ of the $i$ th participant to show missingness depended on the first value $x_i$ which was never missing and the overall missingness rate $p_{miss}$ by a squashed sigmoid function:

Participants were chosen with this probability until exactly $p_{miss}$ of the participant showed missingness. Throughout the simulations, a single data set was created in each condition, and missing data were then simulated in the same data set for each missingness scenario. This data was then analyzed with the misspecified model.

Overall, we simulated 1000 trials in both simulations and recorded the RMSEA with and without the correction. We used Onyx (von Oertzen, Brandmaier, & Tsang, 2015) and lavaan (Rosseel, 2012) for the simulations.

Description of Plots

Results are presented visually by condition, with one plot per type of missingness (MCAR vs. MAR) and observed covariance between variables. The horizontal ( $x$) axes of each plot represent the percent of data missing, while the vertical ( $y$) axes of each plot represent RMSEA. To aid in the interpretation of these simulations, we follow general guidelines in the literature to define RMSEA values less than 0.05 to indicate good fit, greater than 0.1 to indicate poor fit, and values in between 0.05 and 0.1 to indicate mediocre fit (Browne & Cudeck, 1993). A line for RMSEA of 0.05 is included where necessary.

Simulation results revealed consistent patterns for RMSEA calculations. RMSEA values were most strongly related to model misspecification (i.e., the strength of the correlation between variables that was constrained to zero in the misspecified model), but also showed effects of missingness patterns (MCAR vs. MAR). In the bivariate model, for both covariance values of 0.0 and 0.125, MAR and MCAR data generally yielded similar patterns.

Results in the Latent Growth Curve Model

Results of the LGCM simulation are shown in Figure 5. The corrected RMSEA value stays mostly constant up to a rate where 80% of the participants show missingness both under MAR and under MCAR. In contrast, both for MAR and MCAR the uncorrected RMSEA reduces from 0.32 initially to 0.15 for very high missingness rates.

Results Summary

Among models with some misspecification, we see a a downward bias with the uncorrected RMSEA calculation and largely constant values for the $KL$ corrected RMSEA. This shows up among both MCAR and MAR data simulations among models with both low (cov = 0.125) and high (cov = 0.9) model misspecification.

For the models having no misspecification, a pessimistic trend for KL-corrected values is observed when the amount of missing data is larger but not when it is small. In this case, the uncorrected RMSEA is unbiased and the KL-correction shows a slight pessimistic bias.

Limitations

Although the simulations strongly support the notion that $KL$-corrected RMSEA values are unaffected by missing data, further investigation is needed for more general models. The model used in our simulation was minimal and was designed as a proof of concept that demonstrates the behavior of our proposed corrections In most cases, such models are oversimplified and do not adequately represent more complicated datasets. Increasing the complexity of data simulation models and reassessing RMSEA behavior will allow for broader confirmation of the effectiveness of the $KL$ correction as a modification for the RMSEA fit index.

A possible objection to implementing the $KL$ correction for RMSEA lies in the realization that, in presence of missing data, corrected RMSEA values indicate worse model fit than uncorrected values. Though our rationale for an RMSEA correction is rooted in theory and practice, it may be hard to accept that RMSEA values that have been reported thus far have been impling artificially inflated goodness-of-fit for models estimated from missing data. Accepting the $KL$ correction will pressure researchers to reassess model fit, and there may be resistance if $KL$ corrected RMSEA values point to a different conclusion than uncorrected RMSEA values had indicated previously. Considering that the typically used cut-off values are based on rules-of-thumb derived from observed RMSEA values without using the correction, it may be a viable strategy to re-think the typically employed cut-off values. Changing the RMSEA calculation under missing data at large, with higher RMSEA values considered acceptable, may reduce resistance to such a change without necessarily accepting considerably worse models as acceptable but rather restore fairness among models derived from completely and partially observed data. RMSEA values may then become comparable across different levels of missingness, and preference will no longer be given to models where missing data rates are high.

Summary and Future Directions

The original formulation of RMSEA was derived under the assumption of complete observations. We found that a naive computation of the RMSEA under missing data results in overoptimistic model fit (i.e., RMSEA values that are too low). Hence, we must conclude that most reported model fits using RMSEA are too optimistic. We propose that the RMSEA can be corrected by replacing the ${\chi }^2$ with an estimator based on the $KL$ divergence in the RMSEA equation. This leads to an intuitive extension of the original idea behind the RMSEA to settings of missing data.

This paradigm can and should be extended to other fit indices as well, and be incorporated into the updating of SEM fit statistics to account for novel estimation approaches. Other fit indices are affected by missing data, and the $KL$ correction should be explored and applied to CFI, TLI, and other indices. As the $KL$ correction is an estimate, estimation accuracy needs to be considered when creating confidence intervals around $KL$-corrected RMSEA. Beyond these corrections, we must accept and further explore ways of evaluating SEM model fit that do not depend on the assumption of a complete saturated model, especially as we extend SEM into planned missingness, design and definition variables, and other important advances that take SEM from a covariance method to flexible GLM-based modeling.

Informing researchers about the benefits of an RMSEA correction and its ability to handle missing data may encourage the use of SEM in the realistic scenario of missing data. In order to do this, it is important to provide the scientific community easy access to the correction. A feasible solution is writing a package for R which includes wrapper tools that will extract data from existing functions and calculate RMSEA corrected values. In this line, we attach a small R script to this article that derives $KL$-corrected RMSEA values for a given estimation result. The script is available in the Supplementary Materials. $\Omega$ nyx is currently the only software providing both uncorrected and corrected RMSEA. We hope that eventually other software developer teams will follow and integrate the $KL$ correction into existing data software tools to encourage and support widespread use of corrected RMSEA estimates.