Evaluation of Residual-Based Fit Statistics for Item Response Theory Models in the Presence of Non-Responses

Abstract

Residual-based fit statistics, which compare observed item statistics (e.g., proportions) with model-implied probabilities, are widely used to evaluate model fit, item fit, and local dependence in item response theory (IRT) models. Despite the prevalence of item non-responses in empirical studies, their impact on these statistics has not been systematically examined. Existing software (package) often applies heuristic treatments (e.g., listwise or pairwise deletion), which can distort fit statistics because missing data further inflate discrepancies between observed and expected proportions. This study evaluates the appropriateness of such treatments through extensive simulation. Results show that deletion methods degrade the accuracy of fit testing: fit indices are inflated under both null and power conditions, with the bias worsening as missingness increases. In addition, the impact of missing data exceeds that of model misspecification. Practical recommendations and alternative methods are discussed to guide applied researchers.

Keywords

item response theory missing data fit statistics residual-based statistics model fit item fit local independence

Introduction

Item response theory (IRT) models are widely used in fields including large-scale educational assessments (e.g., PISA and NAEP; Mislevy et al., 1992; Organisation for Economic Co-operation and Development [OECD], 2024) and the measurement of patient-reported outcomes (e.g., Hansen et al., 2014; Yu et al., 2012). After obtaining item parameter estimates via maximum likelihood estimation (e.g., Bock & Aitkin, 1981), a critical next step is to evaluate model-data fit. This involves testing whether the model and its individual items fit the data well and whether the assumed ability structure (i.e., factor structure) violates the local (conditional) independence assumption. Researchers have proposed various test statistics to evaluate these distinct aspects of model fit; for a comprehensive review, see Sinharay and Monroe (2025).

Among the proposed methodologies, the family of residual-based fit statistics is arguably the most widely utilized in practice. These statistics function by quantifying the discrepancy between observed item response patterns (or their summary statistics) and the probabilities implied by the IRT model. This study focuses on three representative residual-based statistics: $M_{2}$ (Maydeu-Olivares & Joe, 2005, 2006) for overall model fit, $S - X^{2}$ (T. Kang & Chen, 2008; Orlando & Thissen, 2000, 2003) for item-level fit, and $LD - X^{2}$ (Chen & Thissen, 1997) for testing the local independence assumption. Given their implementation in popular software such as flexMIRT (Cai, 2024) and the mirt (Chalmers, 2012) package in R software (R Core Team, 2025), a working familiarity with them is presumed among applied researchers.

The theoretical foundation of these statistics rests on the assumption of a complete data matrix, a condition frequently unmet in applied research. The presence of item non-responses introduces ambiguity into the calculation of observed frequencies, posing a significant challenge for statistical inference. Without a formal framework for handling this ambiguity, researchers often resort to ad hoc procedures. However, the validity of these methods is highly questionable, as their use can introduce systematic bias, reduce statistical power, and distort sampling distributions (Graham, 2009; H. Kang, 2013; Myers, 2011; Newman, 2014; Pepinsky, 2018). The severity of these issues depends on the underlying missing data mechanism, which is typically categorized as missing completely at random (MCAR), missing at random (MAR), or missing not at random (MNAR) (Allison, 2009; Enders, 2022; Little & Rubin, 2019; Rubin, 1976). Consequently, applying fit statistics designed for complete data to incomplete datasets is a methodological concern that threatens the validity of any conclusions drawn about model fit.

This challenge is compounded by the default behaviors of popular IRT software (package). Programs such as flexMIRT and mirt often implement heuristic missing data treatments (e.g., listwise and pairwise deletion) as a practical necessity and due to the absence of an appropriate treatment. Although these methods are theoretically problematic, their actual impact on the performance of the $M_{2}$ , $S - X^{2}$ , and $LD - X^{2}$ statistics remain largely unexamined. This research gap leaves critical questions unanswered regarding how these ad hoc procedures inflate type I error rates or diminish statistical power. Therefore, this study uses an extensive simulation to evaluate these common missing data treatments when applied to overall model fit, item fit, and local dependence statistics under the MAR framework. Our goal is to provide evidence-based guidance for practitioners on interpreting these diagnostic tools when data are incomplete.

To ensure a clear evaluation, this study exclusively focuses on the MAR framework. While the simpler MCAR mechanism is also a possibility, it often represents a less realistic scenario in applied research and a less challenging statistical condition (Heitjan & Basu, 1996; Little, 1988). The MAR framework, where missingness systematically relates to other observed data, provides a more common and relevant case for investigation. Conversely, the MNAR mechanism was deliberately excluded because it fundamentally biases item parameter estimates, which would make it impossible to separate the effect of the missing data from the effect of biased parameters (Galimard et al., 2018; Waterbury, 2019). By adhering to the MAR assumption, for which robust estimation methods are available, this study can isolate and accurately assess the direct influence of missing data on model fit against the complete data benchmark.

Fit Statistics

This study focuses on the logistic version of the graded response model (GRM; Samejima, 1969). Note that GRM is chosen because of its popularity and for the sake of simplicity, while the discussed fit statistics are generally applicable to other model forms, such as three-parameter logistic (3PL), generalized partial credit, and nominal response models. Item responses are represented in an $N \times J$ matrix, where rows and columns denote respondents and items, respectively. For an item $j$ ( $j = 1, \dots, J$ ) with $K_{j}$ categories, the cumulative probability of responding at or above category $k$ ( $k = 0, \dots, K_{j} - 1$ ) is defined as

\Pr (Y_{j} \geq k | θ) = \frac{\exp (a_{j} θ + c_{jk})}{1 + \exp (a_{j} θ + c_{jk})},

(1)

where $a_{j}$ and $c_{jk}$ are the slope and intercept for item $j$ pertaining to category $k$ , with the boundary conditions $\Pr (Y_{j} \geq 0 | θ) = 1$ and $\Pr (Y_{j} \geq K_{j} | θ) = 0$ . The probability of responding in a specific category is the difference between adjacent cumulative probabilities

\Pr (Y_{j} = k | θ) = \Pr (Y_{j} \geq k | θ) - \Pr (Y_{j} \geq k + 1 | θ) .

(2)

When $K_{j} = 2$ , the GRM simplifies to the two-parameter logistic (2PL) model. A restricted form, analogous to the one-parameter logistic (1PL) model, constrains the slope parameters ( $a_{j}$ ) to be equal across all items. For simplicity, we refer to the GRM as a 1PL or 2PL model depending on this constraint. The following sections introduce the three fit statistics central to this study, $M_{2}$ for the overall model fit, $S - X^{2}$ for item fit, and $LD - X^{2}$ for local independence, and discuss their implementation for missing data in the flexMIRT and the mirt.

Model Fit: $M_{2}$

In a test consisting of many items, typical in IRT applications, Pearson’s $X^{2}$ and the likelihood ratio $G^{2}$ result in increased Type I errors because the resulting contingency table is highly sparse given fixed samples. $M_{2}$ statistics (Cai & Hansen, 2013; Cai & Monroe, 2014; Maydeu-Olivares & Joe, 2005, 2006) are a special case of limited information goodness-of-fit (LI–GOF) statistics (Joe & Maydeu-Olivares, 2010; Maydeu-Olivares & Joe, 2005, 2006; Reiser, 1996), which have been known to resolve this issue. Instead of using full item response pattern sample proportions and expected probabilities, $M_{2}$ statistics utilize linearly independent first- and second-order marginals to eliminate the sparseness of multidimensional contingency tables.

In a test with $J$ items, the number of cells is $C = Π_{j = 1}^{J} K_{j}$ for the full item response patterns. Let $p$ and $π (γ)$ be vectors of population and model-implied probabilities for $C$ full item response patterns, respectively, where $γ$ is a $q \times 1$ IRT model parameter vector. For $M_{2}$ statistics, considered are the linearly independent first- and second-order marginal moments, denoted by $p_{2}$ and $π_{2} (γ)$ with the length of $d = \sum_{j = 1}^{J} (K_{j} - 1) + \sum_{j' = 2}^{J} \sum_{j = 1}^{j' - 1} (K_{j} - 1) (K_{j^{'}} - 1)$ (Cai & Hansen, 2013; Maydeu-Olivares & Joe, 2006). $p_{2}$ and $π_{2} (γ)$ are obtained by a linear transformation matrix $L_{2}$ , such that $p_{2} = L_{2} p$ and $π_{2} (γ) = L_{2} π (γ)$ . An example of $J = 3$ and $K_{j} = 2$ is presented in Section S1 of the supplemental material.

When maximum likelihood estimates (MLEs) $\hat{γ}$ are used for a sample size $N$ , Maydeu-Olivares and Joe (2005, 2006) showed that

\sqrt{N} {\hat{e}}_{2} = \sqrt{N} ({\hat{p}}_{2} - π_{2} (\hat{γ})) \to^{D} N_{d} (0, Ω_{2} (γ)),

(3)

where ${\hat{p}}_{2}$ is the vector of the sample proportions corresponding to $p$ ; the covariance matrix $Ω_{2} (γ) = Ξ_{2} (γ) - Δ_{2} (γ) I^{- 1} (γ) Δ_{2}^{T} (γ)$ ; $Ξ_{2} (γ) = L_{2} Ξ (γ) L_{2}^{T}$ ; $Ξ (γ) = diag (π (γ)) - π (γ) π^{T} (γ)$ , the multinomial covariance matrix; $Δ_{2} (γ) = \partial π_{2} (γ) / \partial γ^{T} = L_{2} \partial π (γ) / \partial γ^{T}$ , the Jacobian matrix; and $I (γ)$ is the Fisher information matrix.

Finally, the $M_{2}$ statistics are defined by

M_{2} = N {\hat{e}}_{2}^{T} Δ_{2}^{(c)} (\hat{γ}) {[Δ_{2}^{(c) T} (\hat{γ}) Ξ_{2} (\hat{γ}) Δ_{2}^{(c)} (\hat{γ})]}^{- 1} Δ_{2}^{(c) T} (\hat{γ}) {\hat{e}}_{2},

(4)

where $Δ_{2}^{(c)} (\hat{γ})$ is an orthogonal complement to $Δ_{2} (\hat{γ})$ evaluated at the MLEs. Under the null hypothesis, the $M_{2}$ statistic is $χ^{2}$ -distributed with the degrees of freedom $df = d - q$ . For polytomous items, as the resulting $p_{2}$ may still be sparse, Cai and Hansen (2013) and Cai and Monroe (2014) proposed $M_{2}^{*}$ and $C_{2}$ , respectively, to further reduce the first- and/or second-order margins by utilizing the results of Joe and Maydeu-Olivares (2010). We use $M_{2}^{*}$ statistics for the simulation study as we employ both dichotomous and polytomous items, but call it $M_{2}$ because the statistic is an extension of the original $M_{2}$ .

Item Fit: $S - X^{2}$

Orlando and Thissen (2000) proposed $S - X^{2}$ , formulated as a Pearson’s $X^{2}$ -type chi-square test, to assess the fit of dichotomous items, and T. Kang and Chen (2008) later extended this approach to polytomous items. The $S - X^{2}$ statistics for a specific item is calculated by summing the $X^{2}$ values for that item across groups of respondents with the same observed total score. To elaborate, suppose there are $S + 1$ possible total scores, and for each score $s$ ( $s = 0, \dots, S$ ), there are $N_{s}$ respondents. The $S - X^{2}$ statistic for item $j$ is defined as

S - X_{j}^{2} = \sum_{s = 1}^{S - 1} \sum_{k = 0}^{K_{j} - 1} N_{s} \frac{{(p_{jsk} - π_{jsk} (\hat{γ}))}^{2}}{π_{jsk} (\hat{γ})},

(5)

where $p_{jsk}$ and $π_{jsk} (\hat{γ})$ are the observed proportion and the expected (IRT model-implied) proportion of respondents who answered with category $k$ for item $j$ within the group of respondents with a sum score of $s$ , respectively.

The expected proportion, $π_{jsk} (\hat{γ})$ , is obtained by the Lord–Wingersky recursion formula (Lord & Wingersky, 1984; Thissen et al., 1995). This recursive algorithm is used to compute the likelihood of obtaining a certain total score, both with and without item $j$ . This method has also been extended to multidimensional IRT models by Cai (2015) and Huang and Cai (2021). $π_{jsk} (\hat{γ})$ is given by

π_{jsk} (\hat{γ}) = \frac{\int \Pr (Y_{j} = k | θ) π_{s - k}^{j -} (θ) g (θ) d θ}{\int π_{s} (θ) g (θ) d θ},

(6)

where $\Pr (Y_{j} = k | θ)$ is the item response function as defined in Equation (2), $π_{s - k}^{j -} (θ)$ is the sum-score likelihood at $s - k$ on all items except item $j$ for a given ability level $θ$ , and $π_{s} (θ)$ is the sum-score likelihood at $s$ on all items for the same ability level $θ$ , and the distribution of $θ$ , $g (θ)$ , is assumed to be a standard normal distribution, $N (0, 1^{2})$ .

One should note that some observed ( $p_{jsk}$ ) and expected ( $π_{jsk}$ ) proportions may not be defined for extremely low and high sum scores. For example, if $K_{j} = 3$ (i.e., $k = 0, 1, 2$ ), the sum scores of 0 and 1 have zero expectations to $k = 1, 2$ (the sum score of 0) and $k = 2$ (the sum score of 1), respectively. Analogously, the sum scores of $S$ and $S - 1$ have zero expectations to $k = 0, 1$ (the sum score of $S$ ) and $k = 0$ (the sum score of $S - 1$ ). To address the issue of extremely small $π_{jsk}$ , T. Kang and Chen (2008) suggested collapsing adjacent response categories for sum score $s$ . As a result, $S - X^{2}$ has different degrees of freedom ( $df$ ) across items, which can be defined for item $j$ as $df = S_{j} - q_{j}$ , where $S_{j}$ is the number of the available sum scores (and categories) and $q_{j}$ is the number of estimated item parameters for that item $j$ .

Recently, Han et al. (2023) argued that the $S - X^{2}$ does not follow a $χ^{2}$ distribution and therefore proposed a modification to the original formula. However, we do not consider this new statistic in our study for two primary reasons. First, practitioners are likely unfamiliar with this modification because, to our knowledge, it has not yet been implemented in the two software packages of interest. Second, previous simulation studies have shown that the original $S - X^{2}$ recovers nominal rejection rates well at $α = 0.05$ when the sample size is sufficiently large (e.g., $N > 1, 000$ ; Chon et al., 2010; T. Kang & Chen, 2008; Orlando & Thissen, 2000). Rather than a rigorous evaluation of the statistic’s distributional properties, our study seeks to determine whether systematic differences in the performance of the original $S - X^{2}$ occur due to the presence of missing data.

Local Dependence Index: $LD - X^{2}$

The most critical assumption in IRT is the local independence assumption, under which item responses are independent given latent variables (in this study, unidimensional $θ$ ). Two forms of local independence exist, referred to as the strong and weak local independence assumptions, which are closely related to the estimation methods: full and limited information estimations, respectively (Liu & Maydeu-Olivares, 2013; Liu & Thissen, 2012, 2014; Maydeu-Olivares, 2013). In strong local independence, the joint probability of whole items must be the product of individual conditional probabilities, given $θ$ . On the other hand, weak local independence is satisfied if the residual covariance, after partialling $θ$ out, is zero for all item pairs. While strong local independence is the basis of standard full-information maximum likelihood estimation (Bock & Aitkin, 1981) in IRT, it is not practically feasible to examine the full item response patterns. Therefore, the tools for the test of local independence typically rely on examining item pairs.

Chen and Thissen (1997)’s $LD - X^{2}$ utilizes Pearson’s $X^{2}$ for item pairs. For items $j$ and $j^{'}$ ( $j \neq j^{'}$ ) with $K_{j}$ and $K_{j^{'}}$ categories, respectively, two $K_{j} \times K_{j^{'}}$ matrices are constructed for the observed and expected item response probabilities. Using the notations defined earlier, we denote the cells in the two matrices $p_{k k^{'}}$ and $π_{k k^{'}}$ for categories $k$ and $k^{'}$ of items $j$ and $j^{'}$ , respectively. $π_{k k^{'}}$ is given by

π_{k k^{'}} = \int \Pr (Y_{j} = k | θ) \Pr (Y_{j^{'}} = k^{'} | θ) g (θ) d θ,

(7)

at MLEs $\hat{γ}$ . The $LD - X^{2}$ index for the item pair is defined as

LD - X_{j j^{'}}^{2} = N \sum_{k = 0}^{K_{j} - 1} \sum_{k^{'} = 0}^{K_{j^{'}} - 1} \frac{{(p_{k k^{'}} - π_{k k^{'}})}^{2}}{π_{k k^{'}}} .

(8)

The theoretical distribution of $LD - X^{2}$ is not known, but Chen and Thissen (1997) suggested $χ^{2}$ with $df = (2 - 1) \times (2 - 1) = 1$ as a reference distribution for dichotomous items. For polytomous items, $df$ can be considered $(K_{j} - 1) \times (K_{j^{'}} - 1)$ (e.g., Liu & Thissen, 2014). However, the empirical distribution of $LD - X^{2}$ under the null seems to vary as a function of test length and sample size. In Chen and Thissen's (1997) simulation study, the mean and standard deviation of $LD - X^{2}$ gradually increase as test length increases. In Liu and Thissen (2012), decreasing sample size tends to result in larger values. Therefore, it is important to examine the overall patterns of $LD - X^{2}$ values instead of focusing on whether or not the values exceed thresholds at type I error rates in $χ^{2}$ with $df = (K_{j} - 1) \times (K_{j^{'}} - 1)$ , which is also recommended in the flexMIRT manual. In this study, rather than testing the null, we mainly focus on how $LD - X^{2}$ changes as missing data occur.

Missing Data Treatment in flexMIRT and mirt

The treatment of missing data for IRT model fit testing in flexMIRT and mirt is based on heuristic methods, namely listwise and pairwise deletion. In listwise deletion, respondents with any missing responses are removed entirely before computing the fit statistics. If $N^{*}$ out of $N$ respondents have complete data, only the data from those $N^{*}$ respondents are used. In contrast, pairwise deletion uses all respondents but excludes those with missing data only from the specific calculations where their data is needed. For example, in a three-item test where item 3 has missing data, the bivariate proportions for items 1 and 2 are based on all $N$ respondents, while the proportions for items 2 and 3 are based on a smaller subset of respondents.

Table 1 summarizes the default deletion method used in each program. While the general procedure in flexMIRT for statistics like $M_{2}$ and $LD - X^{2}$ amounts to pairwise deletion, there is a key distinction for $S - X^{2}$ . For the $S - X^{2}$ statistic, the actual computations for the observed probabilities effectively amount to listwise deletion. A direct consequence is that $S - X^{2}$ cannot be calculated in study designs where no respondent provides a complete set of item responses (e.g., incomplete randomized blocks design). In addition, flexMIRT has an adjustment for $M_{2}$ and produces the standardized $LD - X^{2}$ . For the $M_{2}$ statistic, each element of the covariance matrix ( $Ξ_{2}$ ) is multiplied by the ratio of the total sample size to the available sample size for each element. The standardized $LD - X^{2}$ is computed by $(X^{2} - df) / \sqrt{2 df}$ since the expected mean and variance are $df$ and $2 df$ if $LD - X^{2}$ is $χ^{2}$ -distributed. However, mirt produces raw $LD - X^{2}$ statistics.

Table 1.

Default Missing Data Treatment Methods for Fit Statistics in flexMIRT and mirt

Fit statistic	flexMIRT	mirt
$M_{2}$	Pairwise	Listwise
$S - X^{2}$	Listwise	Listwise
$LD - X^{2}$	Pairwise	Pairwise

Note. Clarification regarding the flexMIRT procedures was obtained from Vector Psychometric Group via the official software support portal (August 5, 2025).

Deletion methods can introduce significant bias into residual-based statistics because their accuracy depends on correctly computed observed proportions, assuming parameter estimation is unbiased and consistent. To illustrate this, we selected datasets from one condition of the simulation study (2PL model, $N = 3, 000$ , $K = 2$ , for details refer to Complete and Missing Data Generation in Monte Carlo Simulation). Using these datasets, we calculated univariate proportions for an item with missing data and bivariate proportions for an item pair with and without missing data after applying deletion methods. It should be noted that because this illustration considers at most two items at a time, listwise and pairwise deletion methods yield identical results.

Table 2 presents the model-implied probabilities from the 2PL model alongside the sample proportions from the three missing data conditions, averaged over 1,000 replications. As shown in the table, while the proportions from the complete data closely matched their expected values, those calculated using deletion methods diverged substantially as the missing data rate increased. This distortion creates misleading residuals, implying that these naive deletion methods will lead to artificially inflated Type I error rates for fit statistics, even when the model is correct. This effect is exacerbated by larger sample sizes, as the sample size $N$ is a multiplier in the formulae for most fit statistics. The main simulation study will demonstrate these consequences in detail.

Table 2.

Average Observed Proportions Compared to Expected Proportions

Proportion	Univariate	Bivariate
Expected	0.832	0.374
Missing 0%	0.831	0.373
Missing 15%	0.861	0.413
Missing 30%	0.902	0.468

Note: “Expected” indicates the marginal probability using true 2PL parameters. “Missing X%” indicates proportions calculated using pairwise deletion. Values are averaged across 1,000 replications.

Monte Carlo Simulation

A Monte Carlo simulation was conducted to evaluate the performance of residual-based fit statistics under various missing data conditions. This section details the simulation design and analysis procedures. Our study considered both correctly specified (null) and misspecified (power) conditions, tailored to the known strengths of each fit statistic. For example, $S - X^{2}$ statistic is primarily used to detect a failure in model form, such as when a 1PL model is incorrectly fitted to 2PL data. In contrast, the $LD - X^{2}$ statistic is sensitive to a failure in dimensionality, which occurs when a unidimensional model is fitted to multidimensional (e.g., bi-factor) data. As a global fit index, the $M_{2}$ statistic was evaluated under both types of misspecification.

Complete and Missing Data Generation

Complete data were first generated for a 21-item test for two sample sizes: $N = 1, 000$ and $N = 3, 000$ . Six data-generating models were used, spanning unidimensional and bi-factor structures with two and three response categories ( $K = 2, 3$ ). For the unidimensional conditions, the 2PL model was used for dichotomous items ( $K = 2$ ), while the GRM was used for three-category items ( $K = 3$ ). The bi-factor models, used exclusively to simulate a failure in dimensionality, consisted of three specific factors (items 1–7, 8–14, and 15–21) and included both mild and moderate levels of misspecification defined by the size of the specific factor slopes relative to the general factor slopes. A total of 1,000 datasets were generated for each condition using the mirt. The specific procedures for generating the slope and difficulty parameters are detailed below, and the full list of true parameters is available in Table S2.1 in the supplemental material.

The (true) item parameters for data generation were drawn from specific distributions using a slope-difficulty formulation, $a_{j} (θ - b_{j})$ , and subsequently transformed into intercepts where $c_{j} = - a_{j} b_{j}$ . General factor slopes ( $a_{j}$ ) were drawn from a uniform distribution, $a_{j} ~ U (0.50, 1.50)$ . For the bi-factor models, specific factor slopes ( $a_{js}$ ) were generated relative to the variance of the general slopes ( $a_{j}^{2}$ ) to create two levels of misspecification. In the mild condition, specific slopes were generated to account for 10% to 25% of the general slope variance, with values drawn from the range $0.1 a_{j}^{2} \leq a_{sj}^{2} \leq 0.25 a_{j}^{2}$ . In the moderate condition, specific slopes were generated to account for 40% to 50% of the general slope variance, with values drawn from the range $0.4 a_{j}^{2} \leq a_{sj}^{2} \leq 0.5 a_{j}^{2}$ . Item difficulties were generated differently according to the number of response categories. For dichotomous items ( $K = 2$ ), item difficulties ( $b_{j}$ ) were drawn from a uniform distribution, $b_{j} ~ U (- 2, 2)$ . For three-category items ( $K = 3$ ), a sequential process was used to ensure ordered difficulties. The first difficulty parameter ( $b_{j 1}$ ) was drawn from $U (- 2, 0)$ , and the second was generated by adding a positive increment: $b_{j 2} = b_{j 1} + δ_{j}$ , where $δ_{j} ~ U (1, 2)$ .

Missing data were introduced using an MAR mechanism, where lower-performing examinees are more likely to have missing responses. The simulation manipulated the missing data rate at two levels: a small rate (affecting approximately 15% of respondents) and a large rate (30%). Missingness was deliberately targeted for items 7, 14, and 21 to represent easy, moderate, and difficult items with similar item slopes, respectively. The procedure for generating missing data, following Finch (2008) and Chung and Cai (2019), involved several steps. First, sum scores were calculated using the 18 non-target items (all items except 7, 14, and 21). Based on the scores, respondents were grouped into four quartiles (the sum scores by 25% intervals). Higher probabilities of missingness were then assigned to the lower-scoring quartiles. For the small rate, the chances of having a missing response for each quartile (from the lowest to highest score) were 0.375, 0.150, 0.056, and 0.019, respectively. For the large rate, these probabilities were doubled (i.e., 0.750, 0.300, 0.113, and 0.038). When a respondent was selected for missingness, all three target items were deleted, resulting in a monotone missing pattern.

Fitted Models and Results of Interest

The generated datasets were analyzed using both flexMIRT and mirt to assess the performance of the fit statistics, with results from the complete data conditions serving as a benchmark. The analysis focused on investigating changes in both fit indices and empirical rejection rates for the null and power conditions across missing rates. Table 3 summarizes the specific data-generating and fitted models used for each null and misspecified condition. The null conditions were considered only for the datasets that were generated and analyzed by the 2PL model, and all three fit statistics ( $M_{2}$ , $S - X^{2}$ , and $LD - X^{2}$ ) were examined. Statistical power was evaluated against two misspecification scenarios. The first, a failure in model form, involved fitting a 1PL model to 2PL-generated data, for which the performance of $M_{2}$ and $S - X^{2}$ was assessed. The second, a failure in dimensionality, involved fitting a unidimensional model to bi-factor data, for which the performance of $M_{2}$ and $LD - X^{2}$ was assessed.

Table 3.

Summary of Data-Generating and Fitted Models for Null and Power Conditions

		Data generating			Fitted
#	$K$	Dimension	Type	Specific slopes	Dimension	Type	Remark
1. Null Rejection Rates
1	2	Uni	2PL		Uni	2PL
2	3	Uni	2PL		Uni	2PL
2. Power—Failure in Model Forms
3	2	Uni	2PL		Uni	1PL	# 1
4	3	Uni	2PL		Uni	1PL	# 2
3. Power—Failure in Dimensionality
5	2	Bi	2PL	Mild	Uni	2PL	# 1
6	2	Bi	2PL	Moderate	Uni	2PL	# 1
7	3	Bi	2PL	Mild	Uni	2PL	# 2
8	3	Bi	2PL	Moderate	Uni	2PL	# 2

Note. In the “Remark” column, the numbers indicate the paired null conditions for each misspecified model.

For $M_{2}$ statistics, three key metrics were summarized: the average ratio of the $M_{2}$ statistic to its degrees of freedom ( ${\bar{M}}_{2} / df$ ), the empirical rejection rate at $α = 0.05$ , and the average of the root mean squared error of approximation (RMSEA).

For $S - X^{2}$ statistics, we considered the average ratio of $S - X^{2}$ to its $df$ ( $\bar{S - X^{2}} / df$ ) and the empirical rejection rates at $α = 0.05$ , as the degrees of freedom for $S - X^{2}$ differed across items and raw score distributions in each replication. Items without missing data were categorized into “middle slope” (central 50%) and “extreme slope” (outer 50%) groups based on their true slope parameters because Orlando and Thissen (2000) found that items with extreme slopes were more likely to exhibit greater statistical power when 1PL models were fit. Note that the target missing items (i.e., items 7, 14, and 21) can be categorized into “middle slope,” as found in Table S2.1.

Finally, to evaluate local dependence, the $LD - X^{2}$ statistics were summarized for meaningful item pair groups. In the null condition, pairs were categorized based on the missing data status of the items involved: (1) pairs without missing data, (2) pairs with one item containing missing data, and (3) pairs with both items containing missing data. For the power condition (i.e., failure in dimensionality), the analysis was further stratified by the specific dimensions ( $s = 1, 2, 3$ ) to which items belonged, allowing for an examination of the joint impact of missing data and violations of local dependence. Note that the metric differs by software. For flexMIRT, the average of the absolute values of the standardized $LD - X^{2}$ was presented, while for mirt, the average of the $LD - X^{2} / df$ ratio was examined.

The ratio-based measures in ${\bar{M}}_{2} / df$ , $\bar{S - X^{2}} / df$ , and $LD - X^{2} / df$ (for mirt) are expected to be close to 1 when a model is correctly specified, as the expected value of a $χ^{2}$ variable is its $df$ . In contrast, the standardized $LD - X^{2}$ in flexMIRT is expected to be close to 0 under the null condition. If the missing data handling methods implemented in the two software packages are effective, the trends observed between complete and missing data conditions should be similar in both the null and power conditions.

Results

Result 1: $M_{2}$

For the null and power conditions, Figure 1 presents the trends for the mean $M_{2} / df$ ratio ${\bar{M}}_{2} / df$ , $p$ -value at $α = 0.05$ , and averaged RMSEA across the simulation conditions. The corresponding values are presented in Tables S2.2, S2.3, and S2.4 in the supplemental material for the null, first power, and second power conditions, respectively.

Null

As shown, when item responses were complete, the indices from the two packages were identical; ${\bar{M}}_{2} / df$ values were close to 1 (Figure 1A), and the empirical rejection rates were close to $α = 0.05$ level (Figure 1B). However, in the presence of missing data, the two software packages yielded different results due to their distinct missing data handling methods. In flexMIRT, which uses a pairwise deletion with an adjustment to $Ξ_{2}$ , the results varied by $K$ . For the $K = 3$ conditions, all indices of misfit increased as the missing data rate increased. In contrast, for the $K = 2$ condition, the test became more conservative at the 30% missing rate than at the 15% rate. Referring to Table S2.2, for example, with $N = 1, 000$ and $K = 2$ , the empirical rejection ( $p$ -value) rate was 0.816 for the 15% missing condition but dropped to 0.005 for the 30% missing condition, a value much lower than both the nominal alpha level and the rate for complete data (0.037).

In contrast, mirt, which implements a listwise deletion method, showed a consistent pattern where all misfit indices increased as the missing data rate increased. Given that mirt does not implement any additional adjustment, we speculate that the distinct behavior observed in flexMIRT for $K = 2$ is attributable to the adjustment in $Ξ_{2}$ . We also observed that a larger sample size exacerbated the negative impact of missing data. In $K = 3$ , both software packages generally showed larger ${\bar{M}}_{2} / df$ , $p$ -value, and RMSEA values when $N = 3, 000$ compared to $N = 1, 000$ . This trend can be explained by the formulation of the $M_{2}$ statistic (Equation 4) and the trends in Table 2. As missing rates increase, the observed response proportions deviate further from the expected proportions; these deviations are then amplified by larger sample sizes, leading to greater misfit.

Power

In both power conditions (failure in model form and failure in dimensionality), the impact of missing item responses followed a similar trend to that observed in the null condition. The presence of missing data appeared to increase the statistical power of the hypothesis tests, as misfit indices generally increased with the rate of missingness. However, this interpretation requires caution. For instance, in the first model misspecification scenario (Figure 1; Table S2.3), the empirical rejection rate and mean RMSEA were 0.336 and 0.010, respectively, for the complete data condition where $N = 1, 000$ and $K = 3$ . In this case, the statistical power was not high and the model fit was acceptable according to the RMSEA, despite the misspecified item slopes. In contrast, when the missing data rate was 30%, the null hypothesis was always rejected ( $p = 1.000$ ) and the mean RMSEA increased significantly to 0.050. This suggests that the presence of missing data can inflate misfit indices, potentially leading researchers to overestimate the true impact of the underlying model misspecification.

Result 2: $S - X^{2}$

For the null and power conditions, Figures 2 and 3 present the trends for the $\bar{S - X^{2}} / df$ ratio and the empirical rejection $p$ -value rate at $α = 0.05$ for $K = 2$ and $K = 3$ , respectively. In the supplemental material, the corresponding values for the null and power conditions are presented for each program in Tables S2.5, S2.6, S2.7, and S2.8.

Null

When the data were complete (0% of missing), the two programs produced almost identical results in both $\bar{S - X^{2}} / df$ ’s and $p$ -values; $\bar{S - X^{2}} / df$ ’s were close to 1, and $p$ -values were close to 0.05 in all other manipulating conditions. Numerical differences between the two programs in this case can be considered as the differences in how to collapse the sum scores with extremely small expected frequencies.

In contrast, the two programs resulted in different trends for items with or without missing responses in both $K = 2$ and $K = 3$ . For items without missing responses, flexMIRT resulted in higher $\bar{S - X^{2}} / df$ ’s and $p$ -values in general as missing rates increased, while mirt produced $\bar{S - X^{2}} / df$ ’s and $p$ -values close to 1 and 0.05, regardless of missing rates, respectively. The result of flexMIRT was exacerbated when both sample size and missing rate were large. $K = 2$ tended to show higher $\bar{S - X^{2}} / df$ and $p$ -values with greater variability than $K = 3$ . In terms of $p$ -value for $K = 2$ , for example, items with extreme slopes showed the range of 0.154 to 0.608 in flexMIRT with $N = 3, 000$ and 30% of missing rates (Figure 2C and Table S2.5). On the other hand, the $p$ -values ranged from 0.053 to 0.084 in mirt for the corresponding condition (Figure 2C and Table S2.6), which are close to the nominal $α = 0.05$ , though a little bit inflated.

Recall that items with missing values had similar slopes to items within the “middle slope” category. For those items, $\bar{S - X^{2}} / df$ and $p$ -value of flexMIRT retained about 1 and 0.05 without regard to missing rates and item difficulty levels. Across conditions, the $p$ -values ranged from 0.042 to 0.066, from 0.045 to 0.075, and from 0.045 to 0.067 for items 7, 14, and 21, representing easy, moderately difficult, and difficult items, respectively (Figures 2C and 3C, and Table S2.5). In opposition, $\bar{S - X^{2}} / df$ and $p$ -value of mirt deviated from 1 and 0.05, as missing rate increased. This was highlighted in items with low and moderate difficulties (items 7 and 14). $p$ -values of these items increased as a function of missing rates, in which the increment was much greater when $N = 3, 000$ and $K = 2$ —from 0.043 to 0.517 (item 7), from 0.050 to 0.515 (item 14), and $K = 3$ —from 0.049 to 0.199 (item 7), and from 0.052 to 0.128 (item 14) (Figures 2C and 3C and Table S2.6). On the other hand, item 21, with high difficulty maintained $\bar{S - X^{2}} / df$ and $p$ -value approximately at 1 and 0.05, when $N = 1, 000$ . When N = 3,000 and as missing rates increased, the $p$ -values of item 21 were from 0.051 to 0.136 for $K = 2$ and from 0.046 to 0.086 for $K = 3$ (Figures 2C and 3C, and Table S2.6), of which the increment was not as substantial as items 7 and 14.

Although both programs utilize a listwise deletion method for computing $S - X^{2}$ in the presence of missing responses, the reason for the different trends observed between items with and without missing data is not immediately clear. We speculate that this discrepancy arised from the different score table collapsing rules implemented in flexMIRT and mirt. These differing rules, evidenced by the small numerical differences observed even in the complete data case, likely produced more substantial divergences when listwise deletion is applied to incomplete data.

Power: Failure in Model Form

Like the null condition, the two programs provided almost identical results for both types of items for the complete data. The power was greater when the sample size was large, and especially for “extreme slope” items in terms of both $\bar{S - X^{2}} / df$ ’s and $p$ -values. For instance, with $K = 2$ , “extreme slope” items’ $p$ -value in flexMIRT ranged from 0.171 to 0.975 when $N = 1, 000$ , while it ranged from 0.641 to 1.000 when $N = 3, 000$ . In contrast, that of “middle slope” items ranged from 0.066 to 0.223 and from 0.068 to 0.775 when $N = 1, 000, 3, 000$ , respectively. When missing data occurred, the two programs showed different trends for items without missing responses and with middle slope and items with missing responses, although generally more powerful with more sample sizes.

For items without missing responses and with middle slopes, flexMIRT resulted in a similar power at 15% of missing to the complete data case, but the $\bar{S - X^{2}} / df$ and $p$ -value substantially increased at 30% of missing. On the other hand, mirt had decreasing $\bar{S - X^{2}} / df$ ’s and $p$ -values in terms of both median and range, as the missing rate increased. For items without missing responses and with extreme slopes, both programs had similar ranges of $\bar{S - X^{2}} / df$ and $p$ -value regardless of missing rates. It means that the item fit test is effective, even though missing data exist. For items with missing responses, regardless of difficulty level, the two programs showed different trends. In flexMIRT, the power gradually decreased for all items as missing rates increased. In mirt, on the other hand, the power gradually increased for all items as missing rates increased.

Result 3: $LD - X^{2}$

As a reminder, flexMIRT produces the standardized $LD - X^{2}$ , computed as $(X^{2} - df) / \sqrt{2 df}$ , while mirt produces the raw $LD - X^{2}$ statistic. Therefore, the values presented in the following results are the averages of the absolute standardized $LD - X^{2}$ for flexMIRT and the $LD - X^{2} / df$ ratio for mirt, calculated for each type of item pair. Our interpretation thus focuses on the changes in these $LD - X^{2}$ indices across different missing rates.

Null

Figure 4 shows the trends of the $LD - X^{2}$ statistics from both packages under the null condition, with corresponding values presented in Table S2.9. The trends were identical for both software packages with respect to the number of item categories, sample size, and missing data rate. As expected, when the data were complete, all types of item pairs showed equivalent $LD - X^{2}$ indices. We observed that although the indices were slightly larger for $K = 3$ , there were no substantial differences across various sample sizes or item-pair types.

However, as the rate of missingness increased, the $LD - X^{2}$ indices also increased, except for pairs composed of two items without missing data; these pairs retained the same values as in the complete data case. For pairs with at least one item containing missing data, the $LD - X^{2}$ indices increased as a function of both sample size and missing rate. This could lead a researcher unaware of these effects to falsely conclude that such item pairs exhibit local dependence, even when they are not locally dependent. When comparing the two types of item pairs involving missingness (“No missing”—“Missing” vs. “Missing”—“Missing”), the indices tended to be slightly larger for the “No missing”—“Missing” pairs, but these differences were not substantial.

Power: Failure in Dimensionality

Figures 5 and 6 present the trends for the $LD - X^{2}$ indices under the mild and moderate misspecification conditions, respectively; the corresponding values are presented in Tables S2.10 through S2.13 in the supplemental material. First, we confirmed that the trends were substantially similar across both programs. When the data were complete, the calculated $LD - X^{2}$ indices were larger for item pairs within the same specific dimension than for pairs between different specific dimensions. In the moderate misspecification condition for mirt (Table S2.13), for instance, the mean $LD - X^{2} / df$ for within-dimension pairs ranged from 2.47 to 3.81, while it ranged from 0.86 to 1.48 for between-dimension pairs when $N = 1, 000$ and $K = 2$ . This indicates local dependence, where the additional association between items is not appropriately captured by the fitted IRT model.

However, the impact of missing data was consistent with what was observed in the null condition. As missing rates and sample sizes increased, the $LD - X^{2}$ indices for the item pairs involving at least one item with missing data were greatly amplified compared to those for pairs with complete responses. Our observations suggest this impact of missing data is much greater than the impact of the underlying local dependence. For example, under the moderate misspecification in flexMIRT with $N = 1, 000$ , $K = 2$ , and a 15% missing rate, the $LD - X^{2}$ indices for pairs involving missing data were large (ranging from 10.74 to 14.76). In contrast, the indices for pairs without missing data under the same condition were substantially smaller (ranging from 0.74 to 4.63).

Particularly when the missing rates were high (e.g., 30%), the underlying patterns of local dependence became more difficult to detect. For instance, in mirt with $N = 1, 000$ and $K = 3$ , the distinction between within- and between-dimension pairs blurred. The $LD - X^{2} / df$ ranges for within-dimension pairs (e.g., 7.60–9.11 in the mild condition) and between-dimension pairs (e.g., 7.97–9.36 in the mild condition) became highly overlapping. This is in sharp contrast to the clear separation observed in the complete data conditions. These results suggest that detecting true patterns of local dependence is highly challenging when missing data are present. As shown in the examples above, the $LD - X^{2}$ indices no longer clearly distinguish between within- and between-dimension item pairs, unlike in the complete data case. Even if one attempts to account for the impact of missing data by separately examining pairs with and without missingness, that strategy is unlikely to be effective here, as the inflation of the indices due to missing data appears to be much larger than the effect of the underlying local dependence.

Conclusion

This study investigated the performance of three widely used residual-based fit statistics in IRT— $M_{2}$ for overall model fit, $S - X^{2}$ for item-level fit, and $LD - X^{2}$ for local dependence—when confronted with item non-responses, handled by heuristic deletion methods. We conclude this research by providing a summary of the results, suggestions for practitioners, and briefly introducing alternative methods to overcome the impact of missing item responses.

Summary of Findings

Across all fit statistics, we observed that heuristic listwise and pairwise deletion methods, as implemented in flexMIRT and mirt, significantly distort the statistics and can lead to erroneous conclusions about model-data fit. The core issue is that these deletion methods alter sample proportions, causing them to diverge from what would be expected under complete data and thereby creating misleading residuals even when the underlying IRT model is correctly specified. Ironically, this effect was often amplified by larger sample sizes, where greater statistical precision is ensured in parameter estimation. A large sample size magnified the discrepancy between the distorted observed proportions and the model-implied probabilities. Critically, the results showed that the impact of missing data on the fit indices was generally much greater than the impact of an actual model misspecification, as observed in the $M_{2}$ and $LD - X^{2}$ results.

The performance of the fit statistics, particularly for $M_{2}$ and $S - X^{2}$ , also varied by software package, a difference attributable to their distinct missing data treatments and statistical adjustments. For the $M_{2}$ statistic, the listwise deletion in mirt led to a consistent increase in misfit indices as the missing rate grew, whereas the adjusted pairwise deletion in flexMIRT produced inconsistent results, even becoming conservative under certain conditions for dichotomous items. The software-related differences for the $S - X^{2}$ statistic were particularly complex. Although both programs use listwise deletion for $S - X^{2}$ , they likely implement different rules for collapsing score tables, leading to divergent outcomes. For items without missing data, flexMIRT showed inflated misfit as the overall rate of missingness in the dataset increased, while mirt’s results for these items remained stable. Conversely, for items with missing data, this pattern reversed: flexMIRT’s statistics remained stable, while mirt’s deviated from their expected values.

Implications for Practice

Based on these findings, the following suggestions are offered to practitioners. First and foremost, researchers must exercise extreme caution when interpreting the output of $M_{2}$ , $S - X^{2}$ , and $LD - X^{2}$ from standard software (or packages) in the presence of item non-responses, as high values of misfit may be more indicative of the impact of missing data than of a poorly specified model. Practitioners should therefore be trained to recognize specific patterns of distortion; for instance, if $LD - X^{2}$ values are dramatically higher for item pairs involving missingness compared to pairs with complete data; this should be considered a likely artifact of the deletion method rather than evidence of local dependence. Finally, researchers must be aware of the specific missing data procedures their software employs by default (e.g., listwise vs. pairwise deletion), as this study demonstrates that these choices lead to different results and patterns of bias.

Alternative Approaches

Given the severe limitations of heuristic deletion methods, practitioners are strongly encouraged to move away from these approaches for fit evaluation in favor of alternatives where the impact of missing data can be minimized. Employing principled missing data techniques, such as multiple imputation (Rubin, 1987), can alleviate the distortion caused by missing data. The availability of general and flexible software for multiple imputation, including mice (van Buuren & Groothuis-Oudshoorn, 2011) and Blimp (Keller & Enders, 2023), has facilitated the use of multiple imputation. A theoretical advantage of multiple imputation is that sample proportions are less likely to be distorted because they are calculated from complete, albeit imputed, datasets. In categorical data factor analysis literature, studies have shown that for global model fit, multiple imputation successfully recovers nominal rejection rates under null conditions while also providing moderate statistical power, though this power is often slightly lower than that of its complete-data counterpart (Chung & Cai, 2019; Liu & Sriutaisuk, 2020; Sriutaisuk et al., 2025). Even though this is not a direct result about $M_{2}$ , we expect similar results to be produced for $M_{2}$ because $M_{2}$ can be understood as a variant of Browne (1984).

A promising alternative that does not rely on residuals is the Jackknife slope index, proposed by Edwards et al. (2018) to detect local dependence. Because this procedure relies on standard item parameter estimates and their standard errors, it can be readily implemented using output from flexMIRT and mirt. Instead of calculating residuals, the Jackknife slope index utilizes changes in item slopes after removing a target item, weighted by the standard errors of the estimates. A large weighted change in the slope of item $j$ after the removal of item $k$ suggests that items $j$ and $k$ are likely to be locally dependent. Because parameter and standard error estimation under missing data is well-established in IRT, the index is not expected to be substantially influenced by the presence of non-responses. As an illustrative example, Figures S3.1–S3.3 in the supplemental material present the bivariate distributions of Jackknife slope index values from a single dataset used in the simulation. Although this is only one example, the figures highlight the potential of the Jackknife slope index for properly detecting local dependence in the presence of missing data; regardless of the missing data rate, the Jackknife slope index values are highly consistent. We emphasize that this stability is attributable to the robust estimation of item parameters and their corresponding standard errors.

Concluding Remark

It is important to acknowledge the limitations of this study, which also suggest avenues for future research. Our simulation conditions, while comprehensive, did not include mixed-item formats (e.g., a combination of dichotomous and polytomous items), which are common in practice. The analysis could also be expanded; for instance, future studies on the $S - X^{2}$ could investigate scenarios where model misfit is specific to only a subset of items (e.g., some items are fit to a 2PL model while others are fit to a 1PL model). Finally, this study did not include an empirical data example, as the true data-generating model is unknown in real-world settings, making it difficult to definitively evaluate the source of observed misfit.

These limitations point to several important directions for future research. A critical next step is the systematic evaluation of the principled alternatives discussed. Although this study highlighted their theoretical advantages, comprehensive simulations are needed to confirm the effectiveness of multiple imputation for recovering global model fit and to validate the performance of non-residual-based methods such as the Jackknife slope index. Furthermore, this work highlights the need for research into the development of new fit statistics that are inherently robust to missing data.

In conclusion, the findings from our study demand that practitioners interpret results from incomplete data with extreme caution, as apparent misfit is often an artifact of the method rather than a true reflection of model inadequacy. These findings advocate for a broader adoption of principled techniques, such as multiple imputation for overall fit and the non-residual-based Jackknife slope index for local dependence. Furthermore, this work should stimulate future research and prompt software developers to implement safer default procedures for practitioners.

Supplemental Material

sj-pdf-1-epm-10.1177_00131644251393444 – Supplemental material for Evaluation of Residual-Based Fit Statistics for Item Response Theory Models in the Presence of Non-Responses

Supplemental material, sj-pdf-1-epm-10.1177_00131644251393444 for Evaluation of Residual-Based Fit Statistics for Item Response Theory Models in the Presence of Non-Responses by Minho Lee and Juyoung Jung in Educational and Psychological Measurement

Footnotes

Acknowledgements

The authors thank Dr. Hwanggyu Lim for providing valuable comments that improved an earlier draft.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Minho Lee

Juyoung Jung

Supplemental Material

Supplemental material for this article is available online.

References

Allison

P. D.

(2009). Missing data. In Milsap

R. E.

Maydeu-Olivares

(Eds.), The SAGE handbook of quantitative methods in psychology (pp. 72–89). Sage.

Bock

R. D.

Aitkin

(1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46(4), 443–459.

Browne

M. W.

(1984). Asymptotically distribution-free methods for the analysis of covariance structures. British Journal of Mathematical and Statistical Psychology, 37(1), 62–83.

Cai

(2015). Lord–Wingersky algorithm version 2.0 for hierarchical item factor models with applications in test scoring, scale alignment, and model fit testing. Psychometrika, 80(2), 535–559.

Cai

(2024). flexMIRT: Flexible multilevel multidimensional item analysis and test scoring (Version 3.7) [Computer software]. Vector Psychometric Group, LLC.

Cai

Hansen

(2013). Limited-information goodness-of-fit testing of hierarchical item factor models. British Journal of Mathematical and Statistical Psychology, 66(2), 245–276.

Cai

Monroe

(2014). A new statistic for evaluating item response theory models for ordinal data (CRESST Report No. 839). University of California, National Center for Research on Evaluation, Standards, and Student Testing (CRESST).

Chalmers

R. P.

(2012). Mirt: A multidimensional item response theory package for the R Environment. Journal of Statistical Software, 48(6), 1–29.

Chen

W.-H.

Thissen

(1997). Local dependence indexes for item pairs using item response theory. Journal of Educational and Behavioral Statistics, 22(3), 265–289.

10.

Chon

K. H.

Lee

W.-C.

Dunbar

S. B.

(2010). A comparison of item fit statistics for mixed IRT models. Journal of Educational Measurement, 47(3), 318–338.

11.

Chung

Cai

(2019). Alternative multiple imputation inference for categorical structural equation modeling. Multivariate Behavioral Research, 54(3), 323–337.

12.

Edwards

M. C.

Houts

C. R.

Cai

(2018). A diagnostic procedure to detect departures from local independence in item response theory models. Psychological Methods, 23(1), 138–149.

13.

Enders

C. K.

(2022). Applied missing data analysis (2nd ed.). Guilford Press.

14.

Finch

(2008). Estimation of item response theory parameters in the presence of missing data. Journal of Educational Measurement, 45(3), 225–245.

15.

Galimard

J.-E.

Chevret

Curis

Resche-Rigon

(2018). Heckman imputation models for binary or continuous MNAR outcomes and MAR predictors. BMC Medical Research Methodology, 18(1), 90.

16.

Graham

J. W.

(2009). Missing data analysis: Making it work in the real world. Annual Review of Psychology, 60, 549–576.

17.

Han

Sinharay

Johnson

M. S.

Liu

(2023). The standardized S-X2 statistic for assessing item fit. Applied Psychological Measurement, 47(1), 3–18.

18.

Hansen

Cai

Stucky

B. D.

Tucker

J. S.

Shadel

W. G.

Edelen

M. O.

(2014). Methodology for developing and evaluating the PROMIS smoking item banks. Nicotine & Tobacco Research, 16(Suppl. 3), S175–S189.

19.

Heitjan

D. F.

Basu

(1996). Distinguishing “missing at random” and “missing completely at random.” The American Statistician, 50(3), 207–213.

20.

Huang

Cai

(2021). Lord–Wingersky algorithm version 2.5 with applications. Psychometrika, 86(4), 973–993.

21.

Joe

Maydeu-Olivares

(2010). A general family of limited information goodness-of-fit statistics for multinomial data. Psychometrika, 75(3), 393–419.

22.

Kang

(2013). The prevention and handling of the missing data. Korean Journal of Anesthesiology, 64(5), 402–406.

23.

Kang

Chen

T. T.

(2008). Performance of the generalized S-X2 item fit index for polytomous IRT models. Journal of Educational Measurement, 45(4), 391–406.

24.

Keller

B. T.

Enders

C. K.

(2023). Blimp user’s guide (version 3). https://www.appliedmissingdata.com/blimp

25.

Little

R. J. A.

(1988). A test of missing completely at random for multivariate data with missing values. Journal of the American Statistical Association, 83(404), 1198–1202.

26.

Little

R. J. A.

Rubin

D. B.

(2019). Statistical analysis with missing data (3rd ed.). Wiley.

27.

Liu

Maydeu-Olivares

(2013). Local dependence diagnostics in IRT modeling of binary data. Educational and Psychological Measurement, 73(2), 254–274.

28.

Liu

Sriutaisuk

(2020). Evaluation of model fit in structural equation models with ordinal missing data: An examination of the D2 method. Structural Equation Modeling: A Multidisciplinary Journal, 27(4), 561–583.

29.

Liu

Thissen

(2012). Identifying local dependence with a score test statistic based on the bifactor logistic model. Applied Psychological Measurement, 36(8), 670–688.

30.

Liu

Thissen

(2014). Comparing score tests and other local dependence diagnostics for the graded response model. British Journal of Mathematical and Statistical Psychology, 67(3), 496–513.

31.

Lord

F. M.

Wingersky

M. S.

(1984). Comparison of IRT true-score and equipercentile observed-score “equatings.” Applied Psychological Measurement, 8(4), 453–461.

32.

Maydeu-Olivares

(2013). Goodness-of-fit assessment of item response theory models. Measurement: Interdisciplinary Research and Perspectives, 11(3), 71–101.

33.

Maydeu-Olivares

Joe

(2005). Limited- and full-information estimation and goodness-of-fit testing in 2n contingency tables. Journal of the American Statistical Association, 100(471), 1009–1020.

34.

Maydeu-Olivares

Joe

(2006). Limited information goodness-of-fit testing in multidimensional contingency tables. Psychometrika, 71(4), 713–732.

35.

Mislevy

R. J.

Johnson

E. G.

Muraki

(1992). Chapter 3: Scaling procedures in NAEP. Journal of Educational Statistics, 17(2), 131–154.

36.

Myers

T. A.

(2011). Goodbye, listwise deletion: Presenting hot deck imputation as an easy and effective tool for handling missing data. Communication Methods and Measures, 5(4), 297–310.

37.

Newman

D. A.

(2014). Missing data: Five practical guidelines. Organizational Research Methods, 17(4), 372–411.

38.

Organisation for Economic Co-operation and Development. (2024). PISA 2022 technical report. PISA, OECD Publishing.

39.

Orlando

Thissen

(2000). Likelihood-based item-fit indices for dichotomous item response theory models. Applied Psychological Measurement, 24(1), 50–64.

40.

Orlando

Thissen

(2003). Further investigation of the performance of S—X2: An item fit index for use with dichotomous item response theory models. Applied Psychological Measurement, 27(4), 289–298.

41.

Pepinsky

T. B.

(2018). A note on listwise deletion versus multiple imputation. Political Analysis, 26(4), 480–488.

42.

R Core Team. (2025). R: A language and environment for statistical computing (Version 4.5.1) [Computer software]. R Foundation for Statistical Computing. https://www.R-project.org

43.

Reiser

(1996). Analysis of residuals for the multinomial item response model. Psychometrika, 61(3), 509–528.

44.

Rubin

D. B.

(1976). Inference and missing data. Biometrika, 63(3), 581–592.

45.

Rubin

D. B.

(1987). Multiple imputation for nonresponse in surveys. Wiley-Interscience.

46.

Samejima

(1969). Estimation of latent ability using a response pattern of graded scores. Psychometrika, 34(S1), 1–97.

47.

Sinharay

Monroe

(2025). Assessment of fit of item response theory models: A critical review of the status quo and some future directions. British Journal of Mathematical and Statistical Psychology.

48.

Sriutaisuk

Liu

Chung

Kim

(2025). Evaluating imputation-based fit statistics in structural equation modeling with ordinal data: The MI2S approach. Educational and Psychological Measurement, 85(1), 82–113.

49.

Thissen

Pommerich

Billeaud

Williams

V. S. L.

(1995). Item response theory for scores on tests including polytomous items with ordered responses. Applied Psychological Measurement, 19(1), 39–49.

50.

van Buuren

Groothuis-Oudshoorn

. (2011). Mice: Multivariate imputation by chained equations in R. Journal of Statistical Software, 45(3), 1–67.

51.

Waterbury

G. T.

(2019). Missing data and the Rasch model: The effects of missing data mechanisms on item parameter estimation. Journal of Applied Measurement, 20(2), 154–166.

52.

Buysse

D. J.

Germain

Moul

D. E.

Stover

Dodds

N. E.

Johnston

K. L.

Pilkonis

P. A.

(2012). Development of short forms from the PROMISTM sleep disturbance and sleep-related impairment item banks. Behavioral Sleep Medicine, 10(1), 6–24.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.84 MB

Evaluation of Residual-Based Fit Statistics for Item Response Theory Models in the Presence of Non-Responses

Abstract

Keywords

Introduction

Fit Statistics

Model Fit: M 2

Item Fit: S - X 2

Local Dependence Index: LD - X 2

Missing Data Treatment in flexMIRT and mirt

Monte Carlo Simulation

Complete and Missing Data Generation

Fitted Models and Results of Interest

Results

Result 1: M 2

Null

Power

Result 2: S - X 2

Null

Power: Failure in Model Form

Result 3: LD - X 2

Null

Power: Failure in Dimensionality

Conclusion

Summary of Findings

Implications for Practice

Alternative Approaches

Concluding Remark

Supplemental Material

sj-pdf-1-epm-10.1177_00131644251393444 – Supplemental material for Evaluation of Residual-Based Fit Statistics for Item Response Theory Models in the Presence of Non-Responses

Footnotes

Acknowledgements

Declaration of Conflicting Interests

Funding

ORCID iDs

Supplemental Material

References

Supplementary Material

Model Fit: $M_{2}$

Item Fit: $S - X^{2}$

Local Dependence Index: $LD - X^{2}$

Result 1: $M_{2}$

Result 2: $S - X^{2}$

Result 3: $LD - X^{2}$