Confidence Interval Estimation for RMSD and MD Item Fit Statistics

1. Introduction

Item response theory (IRT) models [1,2,3] are multivariate statistical models for analyzing vectors of discrete random variables. IRT models are widely applied in the social sciences, particularly in educational large-scale assessment (LSA; [4,5]) studies in which cognitive tasks are administered to persons.

This article focuses on dichotomous (i.e., binary) random variables. Let $\mathit{X}=(X_1,\dots ,X_I)$ denote a vector of I random variables $X_i$ for $i=1,\dots ,I$, commonly referred to as items or scored item responses. A unidimensional IRT model [6] is a parametric model for the probability distribution $\mathsf{P}(\mathit{X}=\mathit{x})$ for $\mathit{x}=(x_1,\dots ,x_I)\in {\{0,1\}}^I$ given by

$$\begin{array}{c}\mathsf{P}(\mathit{X}=\mathit{x};\mathit{\delta },\mathit{\gamma })\hfill \\ {\displaystyle =\int \prod _{i=1}^I\left[{P_i(\theta ;{\mathit{\gamma }}_i)}^{x_i}{\left(1-P_i(\theta ;{\mathit{\gamma }}_i)\right)}^{1-x_i}\right]\varphi (\theta ;\mathit{\delta })\mathrm{d}\theta ,}\hfill \end{array}$$

(1)

where $\varphi$ denotes the density of the normal distribution with mean $\mu$ and standard deviation (SD) $\sigma$, collected in the distribution parameter $\mathit{\delta }=(\mu ,\sigma )$ for the unidimensional latent variable $\theta$, often called the trait or ability variable. The vector $\mathit{\gamma }=({\mathit{\gamma }}_1,\dots ,{\mathit{\gamma }}_I)$ contains the item parameters for the item response functions (IRFs) $P_i(\theta ;{\mathit{\gamma }}_i)=\mathsf{P}(X_i=1\mid \theta )$ for $i=1,\dots ,I$. The IRF of the two-parameter logistic (2PL) model [7] is

$$P_i(\theta ;{\mathit{\gamma }}_i)=\mathsf{\Psi }\left(a_i(\theta -b_i)\right),$$

(2)

where $a_i$ and $b_i$ are the item discrimination and difficulty parameters, respectively, and $\mathsf{\Psi }\left(x\right)={(1+exp(-x))}^{-1}$ is the logistic function. For the 2PL model, the item parameter vector is ${\mathit{\gamma }}_i=(a_i,b_i)$.

For a sample of N individuals with independently and identically distributed observations ${\mathit{x}}_1,\dots ,{\mathit{x}}_N$ from the distribution of the random variable $\mathit{X}$, the model parameters of the IRT model in (1) can be consistently estimated using marginal maximum likelihood estimation (MML; [8,9,10]). Note that identification constraints are required when estimating item parameters and distribution parameters simultaneously [11].

In LSA applications of IRT models, such as programme for international student assessment (PISA; [12]), a parametric model is imposed for the IRFs in (1). These studies involve a large number of countries, and item parameters are generally assumed to be identical across countries, implying parameter invariance [13,14,15]. In practice, however, this assumption may be violated, as certain items can systematically favor or disadvantage specific countries. This phenomenon is referred to as differential item functioning (DIF; [16,17,18]).

In LSA applications of the IRT model, the imposed invariance of item parameters across countries implies that the assumed IRF represents a misspecification of the true IRT model (1) for a given country (or group). Under this misspecified IRT model, the multivariate random vector $\mathit{X}$ is represented as

$$\begin{array}{c}\mathsf{P}(\mathit{X}=\mathit{x};\mathit{\delta },\mathit{\gamma })\hfill \\ {\displaystyle \simeq \int \prod _{i=1}^I\left[{P_i^{*}(\theta ;{\mathit{\gamma }}_i^{*})}^{x_i}{\left(1-P_i^{*}(\theta ;{\mathit{\gamma }}_i^{*})\right)}^{1-x_i}\right]\varphi (\theta ;\mathit{\delta })\mathrm{d}\theta ,}\hfill \end{array}$$

(3)

where $P_i^{*}$ denotes the assumed IRF and $P_i$ represents the true IRF in the data-generating model for the group. The item parameters are denoted by ${\mathit{\gamma }}_i^{*}$, reflecting that they introduce a model misspecification. In LSA applications, it is generally assumed that replacing $P_i$ with $P_i^{*}$ results in only minimal distortion of the estimated distribution parameters $\mathit{\delta }=(\mu ,\sigma )$.

The assessment of the adequacy of parametric IRFs (i.e., item fit; [19,20]) is a central topic in psychometrics. The discrepancy between the true IRF $P_i$ and the assumed parametric IRF $P_i^{*}$ should be quantified using an appropriate item fit statistic, which can also serve as an effect size measure relatively independent of sample size. Ideally, a clear definition of the corresponding population value of the item fit statistic for an infinite sample should exist. It is also desirable that statistical inference, such as the computation of confidence intervals (CI) for these effect sizes, be available. Of primary interest to applied researchers is the detection of misfitting items i, for which the assumed IRFs $P_i^{*}$ deviate substantially from $P_i$.

A lot of item fit statistics has been proposed in the psychometric literature [20]. The present study targets on the root mean square deviation (RMSD; [21,22,23,24,25]) and the related mean deviation (MD) statistic. The motivation for a thorough investigation of these statistics lies in their operational use in current PISA studies.

The following two subsections outline the definition and estimation of the RMSD and MD item fit statistics.

1.1. RMSD Statistic

The RMSD statistic for item i defined at the population level is defined as (see [26,27])

$${\mathsf{RMSD}}_i=\sqrt{\int {(P_i\left(\theta \right)-P_i^{*}\left(\theta \right))}^{2}w\left(\theta \right)\mathrm{d}\theta },$$

(4)

where w denotes the estimated group density function, $P_i$ represents the true IRF, and $P_i^{*}$ denotes the assumed IRF in the IRT model. The RMSD statistic captures the weighted squared discrepancy between $P_i$ and $P_i^{*}$, and equals 0 for an item that fits perfectly.

Typically, the integral in (4) is approximated using a grid of T points ${\theta }_t$ for $t=1,\dots ,T$ of the $\theta$ variable [28], giving

$${\mathsf{RMSD}}_i=\sqrt{{\displaystyle \sum _{t=1}^T}{(p_{it}-p_{it}^{*})}^2{w}_t},$$

(5)

where $p_{it}=P_i\left({\theta }_t\right)$ and $p_{it}^{*}=P_i^{*}\left({\theta }_t\right)$. The distribution weights $w_t$ are assumed to be proportional to $w\left({\theta }_t\right)$ and sum to 1, i.e., ${\sum }_{t=1}^T{w}_t=1$.

Using the notation ${\mathit{p}}_i=(p_{i1},\dots ,p_{iT})$, ${\mathit{p}}_i^{*}=(p_{i1}^{*},\dots ,p_{iT}^{*})$, $\mathit{w}=(w_1,\dots ,w_T)$, and $\mathit{W}=\mathrm{diag}\left(\mathit{w}\right)$, the RMSD definition in (5) can be written compactly as

$${\mathsf{RMSD}}_i=\sqrt{{({\mathit{p}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{W}({\mathit{p}}_i-{\mathit{p}}_i^{*})}.$$

(6)

As the RMSD statistic in (6) contains the unknown group-specific IRF ${\mathit{p}}_i$, the expression cannot be used directly for estimation. An estimate ${\widehat{\mathit{p}}}_i$ of the group-specific IRF ${\mathit{p}}_i$ is therefore required. This estimation relies on individual posterior distributions, which arise as a by-product of MML estimation in the IRT model (see [29,30]).

Let $h_{nt}=\mathsf{P}({\theta }_t\mid {\mathit{x}}_n)$ denote the estimated posterior probability for person n at grid point ${\theta }_t$. It is computed as

$$h_{nt}={\displaystyle \frac{{\displaystyle w_t\prod _{i=1}^I{\left(p_{it}^{*}\right)}^{x_{ni}}{(1-p_{it}^{*})}^{1-x_{ni}}}}{{\displaystyle \sum _{u=1}^T{w}_u\prod _{i=1}^I{\left(p_{iu}^{*}\right)}^{x_{ni}}{(1-p_{iu}^{*})}^{1-x_{ni}}}}}.$$

(7)

Using this posterior distribution, an estimate for item response probabilities $p_{it}=P_i\left({\theta }_t\right)$ can be obtained by

$${\widehat{p}}_{it}={\displaystyle \frac{{\displaystyle \sum _{n=1}^N{h}_{nt}{x}_{ni}}}{{\displaystyle \sum _{n=1}^N{h}_{nt}}}}.$$

(8)

Define ${\widehat{\mathit{p}}}_i=({\widehat{p}}_{i1},\dots ,{\widehat{p}}_{iT})$. It follows that

$$\mathsf{E}\left({\widehat{\mathit{p}}}_i\right)={\mathit{p}}_i.$$

(9)

Substituting ${\mathit{p}}_i$ with ${\widehat{\mathit{p}}}_i$ in (6) yields the sample RMSD statistic

$${\widehat{\mathsf{RMSD}}}_i=\sqrt{{({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{W}({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*})}.$$

(10)

It is important to emphasize that sampling variability in ${\widehat{\mathit{p}}}_i$ induces corresponding sampling variability in the estimated RMSD statistic ${\widehat{\mathsf{RMSD}}}_i$. Because RMSD item fit statistics are used to identify potentially misfitting items, point estimates of the RMSD statistic should be supplemented with inferential tools such as standard errors, CIs, and tests assessing whether the population RMSD differs from 0 (but see [31]). However, the existing literature on this topic is extremely limited.

Finally, the sample-based RMSD statistic depends on sample size through its expected value, because

$$\begin{array}{c}\mathsf{E}{\left({\widehat{\mathsf{RMSD}}}_i\right)}^2\hfill \\ ={\left({\mathsf{RMSD}}_i\right)}^2+\mathsf{E}\left[{({\widehat{\mathit{p}}}_i-{\mathit{p}}_i)}^{\top }\mathit{W}({\widehat{\mathit{p}}}_i-{\mathit{p}}_i)\right].\hfill \end{array}$$

(11)

Note that this relies on the property that ${\widehat{\mathit{p}}}_i$ is an unbiased estimate of ${\mathit{p}}_i$ (see (9)). The second term in (11) introduces bias determined by the sampling variance of ${\widehat{\mathit{p}}}_i$ (see [32,33]). This bias in ${\left({\widehat{\mathsf{RMSD}}}_i\right)}^2$ carries over to the RMSD estimate ${\widehat{\mathsf{RMSD}}}_i$. As the sampling variance of ${\widehat{\mathit{p}}}_i$ approaches zero with increasing sample size, the sample-based statistic ${\widehat{\mathsf{RMSD}}}_i$ converges to the population value ${\mathsf{RMSD}}_i$. Consequently, statistical inference for ${\widehat{\mathsf{RMSD}}}_i$ at a fixed sample size pertains to inference about the pseudo-true value $\mathsf{E}\left({\widehat{\mathsf{RMSD}}}_i\right)$, which exceeds ${\mathsf{RMSD}}_i$. The pseudo-true RMSD value corresponds to the value of this statistic at a fixed sample size for a given dataset. Note that the proportion of misfitting items affects the pseudo-true value, which is larger for smaller sample sizes [32].

1.2. MD Statistic

The MD statistic for item i at the population level is defined as (see [33])

$${\mathsf{MD}}_i=\int (P_i\left(\theta \right)-P_i^{*}\left(\theta \right))w\left(\theta \right)\mathrm{d}\theta .$$

(12)

In contrast to the RMSD statistic, the MD statistic is signed and may take positive or negative values. It is particularly informative for detecting uniform DIF, where the differences $P_i-P_i^{*}$ are uniformly positive or uniformly negative; such patterns are captured directly by ${\mathsf{MD}}_i$. However, the MD statistic is less suited to identifying misfit in the functional form of IRFs.

Using the discrete grid of $\theta$ values, the MD statistic becomes

$${\mathsf{MD}}_i={({\mathit{p}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{w}.$$

(13)

The corresponding sample-based version is

$${\widehat{\mathsf{MD}}}_i={({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{w}.$$

(14)

Using (9), it follows that

$$\mathsf{E}\left({\widehat{\mathsf{MD}}}_i\right)={\mathsf{MD}}_i.$$

(15)

Hence, the expected value of the MD statistic is independent of sample size, which may be viewed as an advantage over the RMSD statistic.

1.3. Goal of this Paper

Although the RMSD and MD item fit statistics are widely used in LSA applications [12,23,34], methodological work on statistical inference for these statistics remains limited. In particular, the availability of CIs would be a valuable addition for identifying misfitting items. Preliminary developments have been presented in [32], but that work relies on a coarse approximation for quantifying uncertainty in the RMSD statistic. An alternative approach in the same study employs a nonparametric bootstrap [32], which is computationally intensive because it requires repeated estimation of the IRT scaling model for each bootstrap sample.

To address these limitations, the present article develops CI estimation procedures based on asymptotic theory and a less computationally demanding parametric bootstrap. The proposed methods are evaluated through a simulation study.

The remainder of the article is structured as follows. Section 2 outlines the estimation procedures for constructing CIs for the RMSD and MD item fit statistics. Section 3 presents a simulation study evaluating the performance of these CI estimators under various conditions. Finally, Section 4 and Section 5 present the discussion and the conclusion.

2. Confidence Intervals for RMSD and MD Statistics

In this section, CI estimation methods for the RMSD and MD item fit statistics are introduced based on asymptotic theory and parametric bootstrap.

The key observation is that the sampling variability of the sample-based RMSD and MD statistics arises from the sampling variability of the group-specific IRF estimates ${\widehat{\mathit{p}}}_i=({\widehat{p}}_{i1},\dots ,{\widehat{p}}_{iT})$. Consequently, uncertainty quantification for these item fit statistics requires the variance matrix ${\mathit{V}}_i=\mathsf{Var}\left({\widehat{\mathit{p}}}_i\right)$.

An estimate of ${\mathit{V}}_i$ can be obtained using M-estimation theory [35,36,37]. From (8), the estimates ${\widehat{p}}_{it}$ satisfy

$$\sum _{n=1}^N{h}_{nt}(x_{ni}-p_{it})=0\mathrm{for}t=1,\dots ,T.$$

(16)

Hence, the variance matrix ${\mathit{V}}_i$ follows from standard M-estimation arguments [35]. The explicit expressions for the entries of ${\mathit{V}}_i$ are as follows.

For the diagonal entries $v_{i,tt}=\mathsf{Var}\left({\widehat{p}}_{it}\right)$, it holds that

$$v_{i,tt}={\displaystyle \frac{{\displaystyle \sum _{n=1}^N{h}_{nt}^2{(x_{ni}-p_{it})}^2}}{{\displaystyle {\left(\sum _{n=1}^N{h}_{nt}\right)}^2}}}.$$

(17)

For the off-diagonal entries $v_{i,st}=\mathsf{Cov}({\widehat{p}}_{is},{\widehat{p}}_{it})$ with $s\ne t$, it holds that

$$v_{i,st}={\displaystyle \frac{{\displaystyle \sum _{n=1}^N{h}_{ns}{h}_{nt}(x_{ni}-p_{is})(x_{ni}-p_{it})}}{{\displaystyle \left(\sum _{n=1}^N{h}_{ns}\right)\left(\sum _{n=1}^N{h}_{nt}\right)}}}.$$

(18)

Note that these results were already presented in [33,38].

M-estimation theory implies that the estimates ${\widehat{\mathit{p}}}_i$ are approximately multivariate normally distributed with mean vector ${\mathit{p}}_i$ and variance matrix ${\mathit{V}}_i$.

The CI derivations in the next two subsections are based on the estimated variance matrix ${\mathit{V}}_i$.

2.1. Confidence Intervals Based on Asymptotic Theory

CI estimates for the MD and RMSD statistics can now be derived using asymptotic theory. The derivation for the MD item fit statistic is presented first.

Using the definition of the sample-based MD statistic ${\widehat{\mathsf{MD}}}_i$ in (14), its variance is

$$\mathsf{Var}\left({\widehat{\mathsf{MD}}}_i\right)={\mathit{w}}^{\top }{\mathit{V}}_i\mathit{w}.$$

(19)

The corresponding standard error is therefore

$$\mathsf{SE}\left({\widehat{\mathsf{MD}}}_i\right)=\sqrt{{\mathit{w}}^{\top }{\mathit{V}}_i\mathit{w}}.$$

(20)

A 95% CI for ${\widehat{\mathsf{MD}}}_i$ is obtained as

$$\left[{\widehat{\mathsf{MD}}}_i-1.96\mathsf{SE}\left({\widehat{\mathsf{MD}}}_i\right),{\widehat{\mathsf{MD}}}_i+1.96\mathsf{SE}\left({\widehat{\mathsf{MD}}}_i\right)\right].$$

(21)

A CI estimate for the sample-based RMSD statistic ${\widehat{\mathsf{RMSD}}}_i$ is now derived. Let the squared RMSD statistic be denoted as the mean square deviation (MSD). Because ${\widehat{\mathsf{MSD}}}_i={\left({\widehat{\mathsf{RMSD}}}_i\right)}^2$ is a quadratic form of a multivariate normal vector, its distribution does not have a convenient closed-form expression [39]. However, its variance can be approximated as

$$\mathsf{Var}\left({\widehat{\mathsf{MSD}}}_i\right)=4{({\mathit{p}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{W}{\mathit{V}}_i\mathit{W}({\mathit{p}}_i-{\mathit{p}}_i^{*}).$$

(22)

A computable estimate is obtained by substituting ${\mathit{p}}_i$ with ${\widehat{\mathit{p}}}_i$, giving

$$\mathsf{Var}\left({\widehat{\mathsf{MSD}}}_i\right)=4{({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{W}{\mathit{V}}_i\mathit{W}({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*}).$$

(23)

Since the RMSD statistic is the square root of the MSD statistic, applying the delta method [36] to $f\left(x\right)=\sqrt{x}$ yields (see [33])

$$\mathsf{Var}\left({\widehat{\mathsf{RMSD}}}_i\right)={\displaystyle \frac{\mathsf{Var}\left({\widehat{\mathsf{MSD}}}_i\right)}{4{\widehat{\mathsf{MSD}}}_i}}.$$

(24)

The standard error of ${\widehat{\mathsf{RMSD}}}_i$ is therefore

$$\mathsf{SE}\left({\widehat{\mathsf{RMSD}}}_i\right)={\displaystyle \frac{\sqrt{{({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{W}{\mathit{V}}_i\mathit{W}({\widehat{\mathit{p}}}_i-{\mathit{p}}_i^{*})}}{{\widehat{\mathsf{RMSD}}}_i}}.$$

(25)

A 95% CI for the sample-based RMSD statistic is given by

$$\begin{array}{c}[{\widehat{\mathsf{RMSD}}}_i-1.96\mathsf{SE}\left({\widehat{\mathsf{RMSD}}}_i\right),\\ {\widehat{\mathsf{RMSD}}}_i+1.96\mathsf{SE}\left({\widehat{\mathsf{RMSD}}}_i\right)]\end{array}$$

(26)

2.2. Confidence Intervals Based on Parametric Bootstrap

The distribution of the RMSD and MD statistics is approximated numerically via a parametric bootstrap. Note that the IRF estimates ${\widehat{\mathit{p}}}_i$ are approximately multivariate normally distributed with mean vector ${\mathit{p}}_i$ and variance matrix ${\mathit{V}}_i$. Because the item fit statistics are functions of ${\widehat{\mathit{p}}}_i$, the parametric bootstrap generates repeated draws ${\tilde{\mathit{p}}}_i$ from a multivariate normal distribution with mean ${\widehat{\mathit{p}}}_i$ and variance matrix ${\mathit{V}}_i$ and evaluates the statistics for each draw. The resulting empirical distribution of the computed values serves for uncertainty quantification and hypothesis testing (see also [40] for a related approach).

To this end, a vector ${\tilde{\mathit{e}}}_i$ is simulated from a multivariate normal distribution with zero mean and variance matrix ${\mathit{V}}_i$. The bootstrap-based estimated IRF ${\tilde{\mathit{p}}}_i$ is then computed as

$${\tilde{\mathit{p}}}_i={\widehat{\mathit{p}}}_i+{\tilde{\mathit{e}}}_i.$$

(27)

The RMSD statistic for a bootstrap replicate, denoted by ${\widetilde{\mathsf{RMSD}}}_i$, is defined as

$${\widetilde{\mathsf{RMSD}}}_i=\sqrt{{({\tilde{\mathit{p}}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{W}({\tilde{\mathit{p}}}_i-{\mathit{p}}_i^{*})}.$$

(28)

The corresponding MD statistic based on bootstrap is

$${\widetilde{\mathsf{MD}}}_i={({\tilde{\mathit{p}}}_i-{\mathit{p}}_i^{*})}^{\top }\mathit{w}={\widehat{\mathsf{MD}}}_i+{\tilde{\mathit{e}}}_i^{\top }\mathit{w}.$$

(29)

The principal advantage of the numerical bootstrap approach is that the resulting distributions of the RMSD and MD statistics are not required to follow a normal distribution. This article applies two alternative CI estimation procedures based on bootstrap.

First, a CI is constructed under a normality assumption, using the SD of the RMSD or MD estimates across bootstrap samples as the standard error in CI computation. Because the normality assumption may be inappropriate for the RMSD statistic, a bootstrap percentile CI is also employed, with the bounds of the 95% CI defined by the 2.5% and 97.5% empirical percentiles.

Although no substantial differences are expected between the CI estimation approaches for the parametric bootstrap applied to the linear MD statistic, meaningful differences may arise for the nonlinear RMSD statistic. These differences are anticipated because the RMSD statistic is more likely to violate the normality assumption.

3. Simulation Study

In this Simulation Study, the performance of alternative CI estimation methods for the RMSD and MD item fit statistics is examined.

3.1. Method

The 2PL model was used for data generaation and for the IRT scaling model. The mean and SD of the normally distributed factor variable $\theta$ were set to 0 and 1, respectively.

The simulation study utilized $I=40$ items to represent conditions involving sufficiently long tests, which facilitates DIF detection. Ten base items with fixed item discriminations $a_i$ and item difficulties $b_i$ were defined and then duplicated four times to form a test comprising 40 items. All item discriminations $a_i$ were fixed at 1. The item difficulties $b_i$ of the ten base items were set to −1.80, −1.40, −1.00, −0.60, −0.20, 0.20, 0.60, 1.00, 1.40, and 1.80. These ten base item parameters were subsequently duplicated four times.

In total, two of the 40 items were simulated to exhibit uniform DIF. Items j and $10-j+1$ for $j=1,\dots ,5$ were specified to display uniform DIF in their difficulties $b_i$ with values $\delta$ and $-\delta$, respectively. The DIF effect size $\delta$ was set to $-0.6$ and $0.6$, representing large DIF magnitudes [18,27,41].

The sample sizes N were set to 125, 250, 500, 1000, 2000, and 4000, reflecting typical applications of item fit statistics in small-scale or large-scale assessment studies using the 2PL model [4,5,42].

In each of the 6 (sample size N) × 5 (selected DIF items) × 2 (DIF effect size $\delta$) $=60$ simulation conditions, 3000 replications were conducted. Item parameters in the 2PL IRT scaling model were fixed at the parameters of the ten base items. However, the presence of DIF was ignored in the scaling step to allow DIF detection by the MD and RMSD item fit statistics. In the scaling step, only the mean $\mu$ and the SD $\sigma$ were estimated, while the item parameters remained fixed. Based on the individual posterior distributions obtained in the scaling model, the MD and RMSD statistics were computed.

Means, SDs, and skewness were computed for the item fit statistics as descriptive summaries. CIs were obtained using asymptotic theory (ASY; see Section 2.1) and a parametric bootstrap for the residuals ${\tilde{\mathit{e}}}_i$ (see Section 2.2). The parametric bootstrap does not require drawing new person samples or re-estimating the scaling model; it only recomputes the fit statistics based on simulated residuals ${\tilde{\mathit{e}}}_i$. Because the parametric bootstrap of residuals relies on repeated draws from a multivariate normal distribution, it is subject to Monte Carlo error, which decreases only with larger numbers of bootstrap samples.

As an alternative, quasi Monte Carlo methods based on a deterministic distribution that minimizes Monte Carlo error were used. Specifically, 1000 multidimensional data points covering a multidimensional uniform distribution with independent components were generated using a Sobol sequence [43], implemented by the function qrng::sobol() in R (Version 4.4.1; [44]) within the qrng package (Version 0.0-10; [45]). These vectors were transformed via the inverse standard normal distribution function to approximate a multivariate normal distribution. This approach enables a fully numerical evaluation of the distribution of the MD and RMSD statistics in the parametric bootstrap procedure.

CIs under the parametric bootstrap were computed using two approaches: assuming a normal distribution for the bootstrap samples (BNO) and using empirical percentiles (BPE). A confidence level of 95% was applied throughout this Simulation Study.

To evaluate the performance of the alternative CI estimation methods, we computed coverage rates, power rates and Type I error rates. The coverage rate was computed as the percentage rate of events in which the pseudo-true parameter of the MD or RMSD item fit statistic lies within the computed CI. A pseudo-true parameter was computed as the average estimate of the fit statistic in a corresponding simulation condition because it is known that the expected value of the RMSD statistic depends on sample size [32]. Coverage rates within the range of 91% and 98% are considered acceptable.

Moreover, the percentage of cases was computed in which the lower bound of the CI included the value 0.05 for RMSD, and in which the CI for MD included 0.05 or −0.05, serving as a test of minimum effect size or close fit [46]. An RMSD value of 0.05 has been shown to be an effective cutoff for detecting DIF [47]. An exact test of an RMSD value of 0 is not feasible with the BPE method because empirical percentiles are always greater than 0. In addition, testing exact model fit (i.e., a population fit statistic of 0) is of limited practical relevance. For DIF items, the resulting percentage represents the power for rejecting close fit. For non-DIF items, these percentages correspond to Type I error rates.

All analyses for this Simulation Study were conducted in R (Version 4.4.1; [44]) using the sirt::xxirt() function in the sirt package (Version 4.2-133; [48]) to fit the 2PL scaling model. Dedicated R functions were written for computing the RMSD and MD statistics and their associated CI estimation methods. Replication materials are available at https://osf.io/jdf38 (accessed on 12 December 2025).

3.2. Results

Figure 1 displays boxplots of the MD statistic for the third item under conditions with DIF and without DIF. The sample distributions were approximately symmetric, and variability decreased as sample size increased. Moreover, the average MD statistics were relatively independent of sample size.

Table 1 presents descriptive statistics, coverage rates, and power rates for the MD item fit statistic for DIF items as a function of item difficulty $b_i$, DIF effect ${\delta }_i$ and sample size N. The mean of the MD statistic was independent of sample size, whereas the SD decreased as sample size increased. In addition, the MD statistic had smaller absolute values for extreme item difficulties and the largest absolute values for item difficulties near 0. The distribution of the MD statistic was approximately symmetric, as indicated by skewness values close to 0.

The coverage rates for the ASY, BNO, and BPE CI estimation methods were all close to the nominal 95% level, indicating satisfactory performance of the MD item fit statistic for DIF items. Power rates increased with larger sample sizes, and the three CI estimation methods exhibited very similar power characteristics.

Table 2 presents descriptive statistics, coverage rates, and Type I error rates for the MD item fit statistic for non-DIF items as a function of item difficulty $b_i$ and sample size N. The mean of the MD statistic was essentially zero across conditions, indicating that the statistic accurately reflected the absence of DIF. Coverage rates for the three CI estimation methods were also close to the nominal 95% level. The Type I error rates were effectively 0%, indicating that false detections of DIF were virtually absent when items did not exhibit DIF in the data-generating model.

Figure 2 displays boxplots of the RMSD statistic for the third item under conditions with and without DIF. The sample distributions showed slight skewness, particularly for smaller sample sizes. Moreover, the average RMSD was larger for smaller sample sizes and converged to a fixed value as the sample size increased.

Table 3 presents descriptive statistics, coverage rates, and power rates for the RMSD item fit statistic for DIF items as a function of item difficulty $b_i$, DIF effect ${\delta }_i$ and sample size N. The mean of the RMSD statistic decreased with increasing sample size, indicating a small-sample bias. As with the MD statistic, the SD of the RMSD statistic decreased as sample size increased. The mean of the RMSD statistic was lower for items with extreme item difficulties. Moreover, notable skewness was observed in small samples, whereas the RMSD distribution for DIF items was approximately symmetric in large samples.

The coverage rates were acceptable for the CI based on asymptotic theory (ASY). In small samples, the percentile-based bootstrap CI (BPE) provided better coverage than the normal distribution-based bootstrap CI (BNO). However, in small samples, the power to detect DIF for items that exhibited DIF in the data-generating model was substantially higher for BPE compared to ASY and BNO.

Table 4 presents descriptive statistics, coverage rates, and Type I error rates for the RMSD item fit statistic for non-DIF items as a function of item difficulty $b_i$ and sample size N. As with the RMSD statistic for DIF items, the mean RMSD statistic decreased with increasing sample size. Substantial skewness was observed in the RMSD distribution, calling into question the validity of CI estimation methods that rely on normality assumptions.

Coverage rates were satisfactory for BPE, whereas inflated coverage was observed for the ASY and BNO methods, both of which depend on the normal distribution assumption. Type I error rates were close to 0, with the exception of BPE in the smallest sample size of $N=125$.

4. Discussion

In this article, CI estimation methods for the MD and RMSD statistics were compared, including those based on asymptotic theory (ASY) and parametric bootstrap. The parametric bootstrap was computationally efficient because only the residuals directly used in the MD and RMSD statistics were resimulated, without requiring repeated estimation of the 2PL scaling model. Such repeated estimation would be substantially more demanding when using large numbers of bootstrap samples. The parametric bootstrap yielded empirical distributions of the resampled MD and RMSD statistics, enabling the computation of bootstrap CIs based on the normal distribution (BNO) and empirical percentiles (BPE).

Overall, the three CI estimation methods exhibited highly similar performance with respect to coverage rates, power rates, and Type I error rates for the MD item fit statistic. This result aligns with the observation that the MD statistic is approximately normally distributed, with no indication of skewness.

The situation differed when considering CI estimation for the RMSD item fit statistic. This statistic exhibited strongly skewed distributions, particularly in small sample sizes. As a consequence of this skewness, the ASY and BNO methods produced unacceptable coverage rates and lower power for rejecting a close fit of RMSD; that is, they were less effective in detecting cases where the sample RMSD was significantly larger than a population RMSD value of 0.05, a threshold indicating close fit and thus an item with no or only minor DIF (or misfit).

As with any simulation study, several limitations apply that might suggest directions for future research on constructing CIs for the MD and RMSD item fit statistics. First, only DIF in item difficulties was examined; DIF in item discriminations was not considered. It is expected that DIF in item discriminations is more difficult to detect with the RMSD statistic, and the MD statistic is generally less suitable for detecting misfitting items in this case. Second, the data-generating model was limited by the absence of variation in item discriminations; that is, all item discriminations were fixed at 1 in the simulation study. This setting is unrealistic in practical datasets but was chosen for convenience in order not to confound the effects of item difficulty on the performance of the fit statistics and their estimated CIs with item discriminations, since the performance of the procedures can be expected to depend on item discriminations. Third, DIF was the sole source of misfit under investigation, although misfit in the functional form of the item response function is also relevant. This situation differs from the assessment of item fit across groups. Misspecifications in the item response functions are frequently assessed within a single group, implying that the uncertainty in estimated item parameters, which also determines the expected item response functions in the RMSD and MD statistics, must be taken into account for CI estimation. M-estimation theory and the parametric bootstrap must be adapted in this case. Fourth, only dichotomous item responses were studied, despite the fact that real-world assessments often include polytomous responses. No conceptual differences are apparent for CI computation in the case of polytomous item responses. However, ambiguity remains regarding how the RMSD statistic can be extended from the dichotomous to the polytomous case.

In educational large-scale assessments, item response data frequently involved student sampling weights and clustered samples. The CI computation techniques presented in this paper relied on the variance matrix of empirical item response probabilities calculated using M-estimation theory (see Section 2). The variance easily accommodated sampling weights by multiplying the case-wise terms by these weights. Moreover, M-estimation theory also allowed variance estimation under stratified cluster sampling, which could be incorporated into the formulas. Alternatively, the required variance matrix for CI computation could be obtained using resampling techniques such as jackknife, balanced repeated replication, or the nonparametric bootstrap.

5. Conclusion

It has been shown that confidence intervals based on both asymptotic theory and the parametric bootstrap work well for the MD item fit statistic, which is approximately symmetric in samples. However, because the RMSD statistic has a skewed distribution, the parametric bootstrap with a percentile-based confidence interval outperforms the confidence interval based on asymptotic theory, which relies on the normal distribution.