Estimation of Standard Error, Linking Error, and Total Error for Robust and Nonrobust Linking Methods in the Two-Parameter Logistic Model

Abstract

The two-parameter logistic (2PL) item response theory model is a statistical model for analyzing multivariate binary data. In this article, two groups are brought onto a common metric using the 2PL model using linking methods. The linking methods of mean–mean linking, mean–geometric–mean linking, and Haebara linking are investigated in nonrobust and robust specifications in the presence of differential item functioning (DIF). M-estimation theory is applied to derive linking errors for the studied linking methods. However, estimated linking errors are prone to sampling error in estimated item parameters, thus resulting in artificially increased the linking error estimates in finite samples. For this reason, a bias-corrected linking error estimate is proposed. The usefulness of the modified linking error estimate is demonstrated in a simulation study. It is shown that a simultaneous assessment of the standard error and linking error in a total error must be conducted to obtain valid statistical inference. In the computation of the total error, using the bias-corrected linking error estimate instead of the usually employed linking error provides more accurate coverage rates.

1. Introduction

Item response theory (IRT) models [1,2,3] are multivariate statistical models for multivariate binary random variables. These models are frequently used to model cognitive testing data from educational or psychological applications. For example, IRT models are operationally utilized in educational large-scale assessments [4,5], like the programme for international student assessment (PISA; [6]) study.

In this article, we only treat unidimensional IRT models [7]. Let $\mathbf{X}=(X_1,\dots ,X_I)$ be the vector of I dichotomous random variables $X_i\in \{0,1\}$ (also referred to as items or (scored) item responses). A unidimensional IRT model [8] is a statistical model for the probability distribution $\mathsf{P}(\mathbf{X}=\mathit{x})$ for $\mathit{x}=(x_1,\dots ,x_I)\in {\{0,1\}}^I$, where

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

(1)

where $\varphi$ denotes the density of the normal distribution, with the mean $\mathsf{\mu }$ and the standard deviation $\mathsf{\sigma }$. The distribution parameters of the latent variable $\mathsf{\theta }$ (also referred to as the factor variable, trait, or ability) are contained in the vector $\mathbf{\delta }=(\mathsf{\mu },\mathsf{\sigma })$. The vector $\mathbf{\gamma }=({\mathbf{\gamma }}_1,\dots ,{\mathbf{\gamma }}_I)$ contains all the estimated item parameters of item response functions (IRFs); where $P_i(\mathsf{\theta };{\mathbf{\gamma }}_i)=\mathsf{P}\left(X_i=1\right|\mathsf{\theta })$ ($i=1,\dots ,I$). The two-parameter logistic (2PL) model [9] possesses the following IRF:

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

(2)

Using the item discrimination $a_i$ and item difficulty $b_i$, $\Psi \left(x\right)={(1+exp(-x))}^{-1}$ denotes the logistic distribution function. The 2PL model could also be estimated for non-normal distributions [10,11,12,13,14,15]. In this case, an identification constraint is typically applied to an item $X_{i_0}$ such that $a_{i_0}=1$ and $b_{i_0}=0$.

The Rasch model [16] is obtained from the 2PL model as a particular case in which all item discriminations equal 1. Some researchers believe that the Rasch model offers particular measurement (i.e., metrological) properties in contrast to the 2PL model (e.g., [17,18,19]). However, in our view, the Rasch model has only one advantage over the 2PL model, which is that a conditional maximum likelihood estimation is applicable [20]. Moreover, the Rasch model possesses the unweighted sum score as a sufficient statistic for $\mathsf{\theta }$, which offers many interpretational advantages [21,22,23,24,25,26]. There is a disbelief that group comparisons can only be conducted with the Rasch model because it has a so-called property of separability, which entails specific objective comparisons [27,28]. However, this reasoning is incorrect and can be disproved by empirical data [29]. In fact, any IRT model with invariant item parameters across groups allows for invariant group comparisons [7,30], although proponents of the Rasch model frequently claim otherwise [31,32].

If independent and identically distributed observations ${\mathit{x}}_1,\dots ,{\mathit{x}}_N$ of N persons from the distribution of the random variable $\mathit{X}$ are available, and the unknown model parameters of the IRT model in (1) can be estimated using the marginal maximum likelihood (MML) using an expectation–maximization algorithm [33,34].

IRT models are frequently used to compare the performance of two groups in a test (i.e., on a set of items) regarding the factor variable $\mathsf{\theta }$ in the IRT model (1). In the following, we only discuss the 2PL model. Two primary approaches can be distinguished [35]. First, concurrent calibration can be applied in which a joint IRT model is estimated in the two groups by assuming common (i.e., invariant) item discriminations and item difficulties in the two groups. While the mean and the standard deviation of $\mathsf{\theta }$ are fixed in the first group for identification reasons, the mean $\mathsf{\mu }$ and the standard deviation $\mathsf{\sigma }$ can be identified for the second group. Hence, these two parameters summarize the group difference regarding the factor variable $\mathsf{\theta }$. Second, the 2PL model can be separately estimated in each of the two groups. This approach allows items to function differently across groups, which is a property that is referred to as differential item functioning (DIF; see [36,37,38]). In the second step, the differences in item parameters are used to determine the group difference regarding the $\mathsf{\theta }$ variable by means of a linking method [39,40,41]. The occurrence of DIF causes additional variability in the estimated mean $\mathsf{\mu }$ and standard deviation $\mathsf{\sigma }$ [42,43,44,45,46]. Therefore, the estimated distribution parameters $(\mathsf{\mu },\mathsf{\sigma })$ depend on the choice of selected items, even for infinite sample sizes of persons. This variability is quantified in the linking error [47,48,49,50,51,52].

An anonymous reviewer pointed out that a simultaneous estimation, while allowing for group-specific DIF effects, would also be possible [53,54,55,56]. In fact, the two-step procedure of separate scaling with subsequent linking can be equivalently formulated as a one-step simultaneous estimation (i.e., concurrent calibration) with nonlinear constraints on group-specific item parameters [57].

This article investigates the computation of the total uncertainty of linking methods in a general treatment based on M-estimation theory [58]. A new bias-corrected linking error is derived, which results in the better performance of the coverage rates of constructed confidence intervals for linking estimates. In particular, it turned out that the newly proposed bias-corrected linking error estimate has a smaller bias for the linking error than the estimators currently employed in the literature.

The rest of the article is organized as follows. Section 2 formalizes the linking methods in the statistical language of estimating equations (i.e., M-estimation). Section 3 presents examples of linking methods that are subsequently investigated in two simulation studies. In Section 4, M-estimation theory is applied to compute the linking error, standard error, total error, and the newly proposed bias-corrected linking error. Then, Section 5 and Section 6 present the findings from two simulation studies. Finally, the article closes with a discussion in Section 7.

2. Linking Method

In this section, we formally define linking methods as a particular M-estimation [58,59] problem. We refer to Section 3 for examples of linking methods. Throughout this paper, it is assumed that ${\widehat{\mathbf{\gamma }}}_i=({\widehat{a}}_{i1},{\widehat{b}}_{i1},{\widehat{a}}_{i2},{\widehat{b}}_{i2})$ contains the estimated item parameters from the 2PL model for item i in both groups. One can assume that there exists a true item parameter ${\mathbf{\gamma }}_i=(a_{i1},b_{i1},a_{i2},b_{i2})$ in the population (i.e., for an infinite sample size). The quantity ${\widehat{\mathbf{\gamma }}}_i$ is a consistent estimate for ${\mathbf{\gamma }}_i$ (under the scheme $N\to \infty$). The goal of a linking method consists of estimating the mean $\mathsf{\mu }$ and the standard deviation $\mathsf{\sigma }$ in the second group based on all the estimated item parameters ${\widehat{\gamma }}_i$.

2.1. One-Step Linking Method

In a one-step method linking method, the parameter estimate $\widehat{\mathbf{\delta }}=(\widehat{\mathsf{\mu }},\widehat{\mathsf{\sigma }})$ is obtained as the minimizer of a nonlinear optimization function H that is a sum of terms, in which each term refers to a single item. Formally, we define

$$\widehat{\mathbf{\delta }}=\underset{\mathbf{\delta }}{arg\; min}H\left(\mathbf{\delta }\right)\mathrm{with}H\left(\mathbf{\delta }\right)=\sum _{i=1}^{I}h(\mathbf{\delta };{\widehat{\mathbf{\gamma }}}_i).$$

(3)

Now, we define the partial derivatives $H_{\mathsf{\mu }}=(\partial H)/(\partial \mathsf{\mu })$, $H_{\mathsf{\sigma }}=(\partial H)/(\partial \mathsf{\sigma })$, $h_{\mathsf{\mu }}=(\partial h)/(\partial \mathsf{\mu })$, $h_{\mathsf{\sigma }}=(\partial h)/(\partial \mathsf{\sigma })$ and ${\mathit{H}}_{\mathbf{\delta }}=(H_{\mathsf{\mu }},H_{\mathsf{\sigma }})$, and ${\mathit{h}}_{\mathbf{\delta }}=(h_{\mathsf{\mu }},h_{\mathsf{\sigma }})$ for brevity. The parameter estimate $\widehat{\mathbf{\delta }}$ solves the following estimating equation:

$${\mathit{H}}_{\mathbf{\delta }}\left(\mathbf{\delta }\right)=\sum _{i=1}^I{\mathit{h}}_{\mathbf{\delta }}(\mathbf{\delta };{\widehat{\mathbf{\gamma }}}_i)=\mathbf{0}.$$

(4)

M-estimation theory is applied for computing a variance estimate for $\widehat{\mathbf{\delta }}$ under the assumption that $I\to \infty$ and ${\left\{{\widehat{\mathbf{\gamma }}}_i\right\}}_{i=1,\dots ,I}$ are independent realizations from a distribution [58]. In this sense, the uncertainty of $\widehat{\mathbf{\delta }}$ regarding the selected set of items is somehow quantified.

2.2. Two-Step Linking Method

In a two-step linking method, the standard deviation $\mathsf{\sigma }$ is estimated in the first step. Afterward, the mean $\mathsf{\mu }$ is estimated in the second step. The estimate $\widehat{\mathsf{\sigma }}$ of the standard deviation $\mathsf{\sigma }$ is determined as a root of the nonlinear equation (that is, it is additive in items).

$$H_{\mathsf{\sigma }}\left(\mathsf{\sigma }\right)=\sum _{i=1}^I{h}_{\mathsf{\sigma }}(\mathsf{\sigma };{\widehat{\mathbf{\gamma }}}_i)=0.$$

(5)

In the second step, the estimate $\widehat{\mathsf{\mu }}$ is obtained as the root of the estimating equation

$$H_{\mathsf{\mu }}(\mathsf{\mu },\widehat{\mathsf{\sigma }})=\sum _{i=1}^I{h}_{\mathsf{\mu }}(\mathsf{\mu },\widehat{\mathsf{\sigma }};{\widehat{\mathbf{\gamma }}}_i)=0.$$

(6)

Note that the two-step method can be alternatively considered as a one-step linking method, because the estimate $\widehat{\mathbf{\delta }}=(\widehat{\mathsf{\mu }},\widehat{\mathsf{\sigma }})$ solves the stacked estimating equation (see [58]).

$${\mathit{H}}_{\mathbf{\delta }}\left(\mathbf{\delta }\right)=\mathbf{0},\mathrm{where}{\mathit{H}}_{\mathbf{\delta }}\left(\mathbf{\delta }\right)=\left(\begin{array}{c}{H}_{\mathsf{\mu }}(\mathsf{\mu },\mathsf{\sigma })\\ H_{\mathsf{\sigma }}\left(\mathsf{\sigma }\right)\end{array}\right)\mathrm{for}\mathbf{\delta }=(\mathsf{\mu },\mathsf{\sigma }).$$

(7)

Hence, one-step and two-step linking methods can be simultaneously analyzed using M-estimation theory.

2.3. Statistical Inference

This article is concerned with estimating the uncertainty of the estimated mean and standard deviation that are contained in the vector $\mathbf{\delta }$. The uncertainty is due to the sampling (or selection) of persons (i.e., under the scheme $N\to \infty$) and the selection (or sampling/modeling the randomness in group comparisons) of items (i.e., under the scheme $I\to \infty$).

Most of the linking literature that treats the uncertainty in $\widehat{\mathbf{\gamma }}$ by computing a standard error (SE) of $\widehat{\mathbf{\delta }}$ for a fixed number of items [60,61,62]. In this case, the variability in $\widehat{\mathbf{\gamma }}$ exists because there is sampling variability in the estimated item parameters ${\widehat{\mathbf{\gamma }}}_i$.

We have argued that the estimated item parameters ${\widehat{\mathbf{\gamma }}}_i$ have a population analogue ${\mathbf{\gamma }}_i$ in an infinite sample size. If there is a model error (i.e., DIF), the estimated linking parameter $\widehat{\mathbf{\delta }}$ depends on the chosen set of items, even for an infinite sample size. This variability due to item selection, which is referred to as the linking error (LE; see [48,63]) and the variance estimation, operates under the scheme where $I\to \infty$.

The total error (TE) includes both sources of uncertainty: the standard error due to randomness due to persons and the linking error due to randomness in items [47,50,51,64]. However, it has been argued that the ordinary linking error estimate is partly affected by the sampling error [65]. In this article, a bias-corrected linking error resulting in a bias-corrected total error estimate is examined to try to reduce the portion of the estimated linking error variance that is due to the sampling error. We outline the statistical underpinnings of the estimators in Section 4.

3. Robust and Nonrobust Linking Methods

In this section, we discuss the most frequently employed linking methods for the 2PL model.

3.1. Mean–Mean Linking (MM)

The mean–mean (MM) linking method is a two-step linking method [40]. The standard deviation $\mathsf{\sigma }$ is estimated in the first step as the ratio of the means of the item discriminations in the two groups.

$$\widehat{\mathsf{\sigma }}=\frac{{\displaystyle \frac{1}{I}\sum _{i=1}^I{\widehat{a}}_{2i}}}{{\displaystyle \frac{1}{I}\sum _{i=1}^I{\widehat{a}}_{1i}}}.$$

(8)

In the second step, the mean $\mathsf{\mu }$ is estimated by

$$\widehat{\mathsf{\mu }}=-\widehat{\mathsf{\sigma }}\frac{1}{I}\sum _{i=1}^I{\widehat{b}}_{i2}+\frac{1}{I}\sum _{i=1}^I{\widehat{b}}_{i1}.$$

(9)

The estimate $\widehat{\mathbf{\delta }}=(\widehat{\mathsf{\mu }},\widehat{\mathsf{\sigma }})$ can be written as the solution $\mathbf{\delta }=(\mathsf{\mu },\mathsf{\sigma })$ of the estimating equations

$$\sum _{i=1}^I\left(\mathsf{\sigma }{\widehat{b}}_{i2}-{\widehat{b}}_{i1}+\mathsf{\mu }\right)=0\mathrm{and}\sum _{i=1}^I\left(\mathsf{\sigma }{\widehat{a}}_{i1}-{\widehat{a}}_{i2}\right)=0.$$

(10)

The MM linking method is considered nonrobust to outliers in the item parameter differences between groups, because the mean is an outlier-sensitive location measure.

3.2. Mean–Geometric–Mean Linking (MGM and RMGM)

The mean–geometric–mean (MGM) linking [40] (also referred to as log–mean–mean linking in [66]) is another two-step linking method that estimates $\mathsf{\sigma }$ as the ratio of the geometric means of the item discriminations. Using a general loss function $\mathsf{\rho }$, the standard deviation $\mathsf{\sigma }$ is estimated in the first step as

$$\widehat{\mathsf{\sigma }}=\underset{\mathsf{\sigma }}{arg\; min}\sum _{i=1}^I\mathsf{\rho }\left(log{\widehat{a}}_{i2}-log\left(\mathsf{\sigma }\right)-log{\widehat{a}}_{i1}\right),$$

(11)

where $\mathsf{\rho }\left(x\right)={(x^2+\mathsf{\epsilon })}^{p/2}$ is a differentiable approximation of the $L_p$ loss function $x\mapsto {\left|x\right|}^p$ for $p\in (0,2]$ and a sufficiently small $\epsilon >0$ (see [67]). In this article, we use $\epsilon =0.01$ in the simulation studies [68]. The ordinary mean is obtained using the $L_p$ loss function for $p=2$. The median as a location measure can be approximately obtained for $p=1$. The choice $p=0.5$ is advantageous for asymmetrically distributed error distributions for DIF effects [69] and appears in invariance alignment [68,70]. In the second step, the mean $\mathsf{\mu }$ is estimated by

$$\widehat{\mathsf{\mu }}=\underset{\mathsf{\mu }}{arg\; min}\sum _{i=1}^I\mathsf{\rho }\left(\widehat{\mathsf{\sigma }}{\widehat{b}}_{i2}-{\widehat{b}}_{i1}+\mathsf{\mu }\right).$$

(12)

By defining ${\mathsf{\rho }}_1=\mathrm{d}\mathsf{\rho }/\mathrm{d}x$, the solution in MGM linking is given as the root of the estimating equations as follows:

$$\sum _{i=1}^I{\mathsf{\rho }}_1\left(log{\widehat{a}}_{i2}-log\left(\mathsf{\sigma }\right)-log{\widehat{a}}_{i1}\right)=0\mathrm{and}\sum _{i=1}^I{\mathsf{\rho }}_1\left(\mathsf{\sigma }{\widehat{b}}_{i2}-{\widehat{b}}_{i1}+\mathsf{\mu }\right)=0.$$

(13)

MGM linking in a nonrobust variant is defined with $p=2$, while robust MGM (abbreviated as RMGM) linking is defined by choosing $p=0.5$ in the loss function $\mathsf{\rho }$.

3.3. Symmetric Haebara Linking (HAE and RHAE)

Haebara linking is a one-step linking method [71] that relies on aligning IRFs to determine $\mathsf{\mu }$ and $\mathsf{\sigma }$. Using a discrete grid of $\mathsf{\theta }$ points ${\mathsf{\theta }}_t$ ($t=1,\dots ,T$) and weights ${\omega }_t$, the mean $\mathsf{\mu }$ and the standard deviation $\mathsf{\sigma }$ are estimated by minimizing a weighted distance between the IRFs, which is shown as follows:

$$H(\mathsf{\mu },\mathsf{\sigma })=\sum _{i=1}^I\sum _{t=1}^T{\omega }_t\mathsf{\rho }\left(P(\mathsf{\sigma }{\mathsf{\theta }}_t+\mathsf{\mu };{\widehat{a}}_{i1},{\widehat{b}}_{i1})-P({\mathsf{\theta }}_t;{\widehat{a}}_{i2},{\widehat{b}}_{i2})\right),$$

(14)

where $P(\mathsf{\theta };a,b)=\Psi \left(a\right(\mathsf{\theta }-b\left)\right)$, and $\mathsf{\rho }$ again define a differentiable approximation of the $L_p$ loss function. For example, the $\mathsf{\theta }$ grid can be equidistantly chosen between $-6$ and 6, and ${\omega }_t$ could be proportional to the values of the density of a normal distribution with a standard deviation of 2. Some researchers alternatively prefer to set all weights ${\omega }_t$ equal to 1. The nonrobust loss function $\mathsf{\rho }\left(x\right)=x^2$ has been originally proposed by Haebara [71] and is widely used. The robust loss function with $p=1$ has been proposed in Refs. [72,73]. The general robust $L_p$ loss function in Haebara linking for $p\le 1$ was utilized in [74]. Corresponding estimating equations to (14) are obtained by computing the partial derivatives of H with respect to $\mathsf{\mu }$ and $\mathsf{\sigma }$.

The linking function H in (14) is referred to as asymmetric Haebara linking, because it aligns the IRF of the first group with the IRF of the second group. This kind of asymmetry induces non-negligibly biased estimates for the standard deviation and for the mean to a smaller extent [66]. To this end, symmetric Haebara linking [75] has been proposed to simultaneously align the IRFs of both groups and to reduce the bias [66]. The linking function of symmetric Haebara linking is defined by

$$H(\mathsf{\mu },\mathsf{\sigma })=\sum _{i=1}^I\sum _{t=1}^T{\omega }_t\mathsf{\rho }\left(P(\mathsf{\sigma }{\mathsf{\theta }}_t+\mathsf{\mu };a_{i1},b_{i1})-P({\mathsf{\theta }}_t;a_{i2},b_{i2})\right)+\sum _{i=1}^I\sum _{t=1}^T{\omega }_t\mathsf{\rho }\left(P({\mathsf{\theta }}_t;a_{i1},b_{i1})-P({\mathsf{\sigma }}^{-1}{\mathsf{\theta }}_t-{\mathsf{\sigma }}^{-1}\mathsf{\mu };a_{i2},b_{i2})\right).$$

(15)

The estimate $\widehat{\mathbf{\delta }}=(\widehat{\mathsf{\mu }},\widehat{\mathsf{\sigma }})$ is obtained by minimizing H in (15) with respect to $\mathbf{\delta }=(\mathsf{\mu },\mathsf{\sigma })$. Nonrobust symmetric Haebara (HAE) linking is obtained with the choice $p=2$ of the $L_p$ loss function, while robust symmetric Haebara (RHAE) linking is obtained with the choice $p=0.5$.

4. Estimation of Standard Error, Linking Error, and Total Error

In this section, we derive the standard error, linking error, and total error for linking estimates in the framework of M-estimation theory [58,59,76]. The treatment in this section is an extension of the material presented in [65]. Assume the item parameter estimate ${\widehat{\mathbf{\gamma }}}_i=({\widehat{a}}_{i1},{\widehat{b}}_{i1},{\widehat{a}}_{i2},{\widehat{b}}_{i2})$ with an estimated variance matrix ${\mathit{V}}_{{\widehat{\mathbf{\gamma }}}_i}$. Moreover, let $\widehat{\mathbf{\gamma }}=({\widehat{\mathbf{\gamma }}}_1,\dots ,{\widehat{\mathbf{\gamma }}}_I)$ be the vector of all the item parameters with an estimated variance matrix ${\mathit{V}}_{\widehat{\mathbf{\gamma }}}$. The corresponding population analogs of the estimators are denoted by ${\mathbf{\gamma }}_i$ and $\mathbf{\gamma }$ and are effectively the estimates for an infinite sample size. As described in Section 2, the linking method provides an estimate $\widehat{\mathbf{\delta }}=(\widehat{\mathsf{\mu }},\widehat{\mathsf{\sigma }})$ for the population parameter $\mathbf{\delta }=(\mathsf{\mu },\mathsf{\sigma })$ as a root of the estimating equation.

$${\mathit{H}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }})=\sum _{i=1}^I{\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)=\mathbf{0}.$$

(16)

4.1. Linking Error

First, we derive the linking error of the estimate $\widehat{\mathbf{\delta }}$ of the estimating equation (16) that quantifies the uncertainty in the estimate due to the selection (or randomness) of items. As is usual in M-estimation theory, we carry out a Taylor approximation of ${\mathit{H}}_{\mathbf{\delta }}$ around the true parameter $\mathbf{\delta }$ (sometimes also referred to as the pseudotrue parameter; [58])

$${\mathit{H}}_{\mathbf{\delta }}\left(\widehat{\mathbf{\delta }}\right)={\mathit{H}}_{\mathbf{\delta }}\left(\mathbf{\delta }\right)+{\mathit{H}}_{\mathbf{\delta }\mathbf{\delta }}\left(\mathbf{\delta }\right)(\widehat{\mathbf{\delta }}-\mathbf{\delta })=\mathbf{0},$$

(17)

where ${\mathit{H}}_{\mathbf{\delta }\mathbf{\delta }}$ denotes the vector of partial derivatives of ${\mathit{H}}_{\mathbf{\delta }}$ with respect to $\mathbf{\delta }$. Moreover, it holds that ${\mathit{H}}_{\mathbf{\delta }}\left(\mathbf{\delta }\right)=0$ because of the definition of the true parameter [58]. Hence, we obtain from (17) the following:

$$\widehat{\mathbf{\delta }}-\mathbf{\delta }=-{\mathit{H}}_{\mathbf{\delta }\mathbf{\delta }}{\left(\mathbf{\delta }\right)}^{-1}{\mathit{H}}_{\mathbf{\delta }}\left(\widehat{\mathbf{\delta }}\right).$$

(18)

M-estimation theory provides the variance matrix of $\widehat{\mathbf{\delta }}$ as the sandwich variance estimate:

$${\mathit{V}}_{\mathrm{LE}}=\mathsf{Var}\left(\widehat{\mathbf{\delta }}\right)={\mathit{A}}^{-1}\mathit{B}{\mathit{A}}^{-\top },\mathrm{where}$$

(19)

$$\mathit{A}={\mathit{H}}_{\mathbf{\delta }\mathbf{\delta }}\left(\mathbf{\delta }\right)=\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\mathsf{\delta }}(\mathbf{\delta };{\widehat{\mathbf{\gamma }}}_i)\mathrm{and}$$

(20)

$$\mathit{B}=\mathsf{Var}\left({\mathit{H}}_{\mathbf{\delta }}\left(\widehat{\mathbf{\delta }}\right)\right)=\sum _{i=1}^I\mathsf{Var}\left({\mathit{h}}_{\mathbf{\delta }}(\mathbf{\delta };{\widehat{\mathbf{\gamma }}}_i)\right).$$

(21)

In (21), we used the approximate independence of the item parameters across items. In M-estimation, the matrix $\mathit{A}$ is called the bread matrix, and $\mathit{B}$ is the meat matrix. The unknown quantities in (19) can be estimated by

$$\widehat{\mathit{A}}=\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\mathsf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)\mathrm{and}$$

(22)

$$\widehat{\mathit{B}}=\sum _{i=1}^I{\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i){\mathit{h}}_{\mathbf{\delta }}{(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)}^{\top }.$$

(23)

Hence, an estimate of the variance matrix ${\mathit{V}}_{\mathrm{LE}}$ is given by

$${\widehat{\mathit{V}}}_{\mathrm{LE}}=\frac{I}{I-1}·{\widehat{\mathit{A}}}^{-1}\widehat{\mathit{B}}{\widehat{\mathit{A}}}^{-\top }.$$

(24)

The factor $I/(I-1)$ in (24) is included to correct for finite-sample bias [65,77,78,79].

4.2. Standard Error

We now compute the standard error of $\widehat{\mathbf{\delta }}$ due to the sampling of persons (see [80,81,82,83,84,85]). A Taylor approximation of ${\mathit{h}}_{\mathbf{\delta }}$ around $(\mathbf{\delta },{\mathbf{\gamma }}_i)$ is carried out and results in

$${\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)={\mathit{h}}_{\mathbf{\delta }}(\mathbf{\delta },{\mathbf{\gamma }}_i)+{\mathit{h}}_{\mathsf{\delta }\gamma }(\mathbf{\delta },{\mathbf{\gamma }}_i)({\widehat{\mathbf{\gamma }}}_i-{\mathbf{\gamma }}_i)+{\mathit{h}}_{\mathsf{\delta }\mathsf{\delta }}(\mathbf{\delta },{\mathbf{\gamma }}_i)(\widehat{\mathbf{\delta }}-\mathbf{\delta }).$$

(25)

We can now use

$$\sum _{i=1}^I{\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)=\mathbf{0}\mathrm{and}\sum _{i=1}^I{\mathit{h}}_{\mathbf{\delta }}(\mathbf{\delta },{\mathbf{\gamma }}_i)=\mathbf{0},\mathrm{and}\mathrm{arrive}\mathrm{at}$$

(26)

$$\widehat{\mathbf{\delta }}-\mathbf{\delta }={\left(\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\mathsf{\delta }}(\mathbf{\delta },{\mathbf{\gamma }}_i)\right)}^{-1}\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\gamma }(\mathbf{\delta },{\mathbf{\gamma }}_i)({\widehat{\mathbf{\gamma }}}_i-{\mathbf{\gamma }}_i)={\mathit{A}}^{-1}\mathit{C}(\widehat{\mathbf{\gamma }}-\mathbf{\gamma }),\mathrm{where}$$

(27)

$$\mathit{C}=\left(\begin{array}{ccc}{\mathit{h}}_{\mathsf{\delta }\gamma }(\mathbf{\delta };{\mathbf{\gamma }}_1)& \cdots & {\mathit{h}}_{\mathsf{\delta }\gamma }(\mathbf{\delta };{\mathbf{\gamma }}_I)\end{array}\right).$$

(28)

This allows us to compute the variance matrix in $\widehat{\mathbf{\delta }}$ due to the sampling error as follows:

$${\mathit{V}}_{\mathrm{SE}}={\mathit{A}}^{-1}\mathit{D}{\mathit{A}}^{-\top }\mathrm{with}\mathit{D}=\mathit{C}{\mathit{V}}_{\widehat{\mathbf{\gamma }}}{\mathit{C}}^{\top }.$$

(29)

The unknown quantities in (29) can be estimated using $\widehat{\mathit{A}}$ in (22), and we obtain the following:

$$\widehat{\mathit{D}}=\widehat{\mathit{C}}{\mathit{V}}_{\widehat{\mathbf{\gamma }}}{\widehat{\mathit{C}}}^{\top }\mathrm{with}\widehat{\mathbf{C}}=\left(\begin{array}{ccc}{\mathit{h}}_{\mathsf{\delta }\gamma }(\widehat{\mathbf{\delta }},{\widehat{\mathbf{\gamma }}}_1)& \dots & {\mathit{h}}_{\mathsf{\delta }\gamma }(\widehat{\mathbf{\delta }},{\widehat{\mathbf{\gamma }}}_I)\end{array}\right).$$

(30)

4.3. Total Error and Bias-Corrected Linking Error

We now compute the total uncertainty in $\widehat{\mathbf{\delta }}$ (i.e., the total error). The variance as the total error has been defined as the sum of the variances due to the sampling error and linking error, and it is written as follows (see [50,65]):

$${\mathit{V}}_{\mathrm{TE}}={\mathit{V}}_{\mathrm{SE}}+{\mathit{V}}_{\mathrm{LE}}\mathrm{and}{\widehat{\mathit{V}}}_{\mathrm{TE}}={\widehat{\mathit{V}}}_{\mathrm{SE}}+{\widehat{\mathit{V}}}_{\mathrm{LE}}.$$

(31)

We now derive a bias-corrected estimate of the linking error variance matrix ${\mathit{V}}_{\mathrm{LE}}$, which, in turn, allows us to compute a bias-corrected variance matrix for the linking error. The estimated meat matrix in the variance matrix for the linking error is given as

$$\widehat{\mathit{B}}=\sum _{i=1}^I{\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i){\mathit{h}}_{\mathbf{\delta }}{(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)}^{\top }.$$

(32)

However, the linking error should only be computed based on the true item parameters ${\mathbf{\gamma }}_i$ instead of the estimated item parameters ${\widehat{\mathbf{\gamma }}}_i$, which appear in (32). A Taylor approximation provides

$${\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\widehat{\mathbf{\gamma }}}_i)={\mathit{h}}_{\mathbf{\delta }}(\widehat{\mathbf{\delta }};{\mathbf{\gamma }}_i)+{\mathit{h}}_{\mathsf{\delta }\gamma }(\widehat{\mathbf{\delta }},{\mathbf{\gamma }}_i)({\widehat{\mathbf{\gamma }}}_i-{\mathbf{\gamma }}_i).$$

(33)

Hence, the inflated variance contribution in $\widehat{\mathit{B}}$ due to the sampling error can be determined as

$$\mathsf{Var}\left(\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\gamma }(\widehat{\mathbf{\delta }},{\mathbf{\gamma }}_i)({\widehat{\mathbf{\gamma }}}_i-{\mathbf{\gamma }}_i)\right)=\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\gamma }(\widehat{\mathbf{\delta }},{\mathbf{\gamma }}_i){\mathit{V}}_{{\widehat{\mathbf{\gamma }}}_i}{\mathit{h}}_{\mathsf{\delta }\gamma }{(\widehat{\mathbf{\delta }},{\mathbf{\gamma }}_i)}^{\top },$$

(34)

where we used the approximate independence of the item parameters across items. As a result, we compute a bias-corrected meat matrix as follows:

$${\widehat{\mathit{B}}}_{\mathrm{bc}}=\widehat{\mathit{B}}-\tilde{\mathit{D}}\mathrm{with}\tilde{\mathit{D}}=\sum _{i=1}^I{\mathit{h}}_{\mathsf{\delta }\gamma }(\widehat{\mathbf{\delta }},{\widehat{\mathbf{\gamma }}}_i){\mathit{V}}_{{\widehat{\mathbf{\gamma }}}_i}{\mathit{h}}_{\mathsf{\delta }\gamma }{(\widehat{\mathbf{\delta }},{\widehat{\mathbf{\gamma }}}_i)}^{\top }$$

(35)

Note that the correction term $\tilde{\mathit{D}}$ in (35) corresponds to the matrix $\widehat{\mathit{D}}$ in (29) in the variance due to standard errors if the item parameters ${\widehat{\mathbf{\gamma }}}_i$ were uncorrelated across items. Next, a bias-corrected variance matrix due to linking error is given as

$${\widehat{\mathit{V}}}_{\mathrm{LE},\mathrm{bc}}=\frac{I}{I-1}·{\widehat{\mathit{A}}}^{-1}{\widehat{\mathit{B}}}_{\mathrm{bc}}{\widehat{\mathit{A}}}^{-\top }$$

(36)

and the variance matrix for the total error is given by

$${\widehat{\mathit{V}}}_{\mathrm{TE},\mathrm{bc}}={\widehat{\mathit{V}}}_{\mathrm{SE}}+{\widehat{\mathit{V}}}_{\mathrm{LE},\mathrm{bc}}.$$

(37)

To sum up, the variance matrix referring to the total error can be written as

$${\widehat{\mathit{V}}}_{\mathrm{TE}}={\widehat{\mathit{A}}}^{-1}\left(\frac{I}{I-1}\widehat{\mathit{B}}+\widehat{\mathit{D}}\right){\widehat{\mathit{A}}}^{-\top }\mathrm{and}{\widehat{\mathit{V}}}_{\mathrm{TE},\mathrm{bc}}={\widehat{\mathit{A}}}^{-1}\left(\frac{I}{I-1}\widehat{\mathit{B}}+\widehat{\mathit{D}}-\frac{I}{I-1}\tilde{\mathit{D}}\right){\widehat{\mathit{A}}}^{-\top }.$$

(38)

To obtain standard errors, linking errors, bias-corrected linking errors, total errors and bias-corrected total errors, the square root of the diagonal elements of the corresponding matrices can be taken. In cases of negative variances for bias-corrected estimates, a corresponding linking error estimate is set to zero.

5. Simulation Study 1: Assessing Linking Error in Infinite Sample Size

In the first Simulation Study 1, the validity of the variance estimates based on M-estimation (see Section 4) for the estimated mean $\widehat{\mathsf{\mu }}$ and the estimated standard deviation $\widehat{\mathsf{\sigma }}$ was investigated for an infinite sample size of persons.

5.1. Method

In this study, only item parameters were simulated in each replication. No item responses were simulated, because the case of an infinite sample size N was investigated. The 2PL model was used as the IRF in the IRT model. There were two groups. For identification reasons, the mean and the standard deviation of the factor variable $\mathsf{\theta }$ in the first group were set to 0 to 1, respectively. The mean $\mathsf{\mu }$ and the standard deviation $\mathsf{\sigma }$ of the second group parametrized the group differences. Throughout the simulation, the choices $\mathsf{\mu }=0.3$ and $\mathsf{\sigma }=1.2$ were made.

We simulated item parameters for $I=10$, 20, and 40 items. The group-specific item parameters $a_{ig}$ and $b_{ig}$ for $i=1,\dots ,I$ and $g=1,2$ relied on base item parameters that were fixed in the simulation and a random DIF effect that was simulated in each replication of the simulation study. The base item discriminations $a_{i0}$ in the case of $I=10$ items were chosen as 0.73, 1.25, 1.20, 1.47, 0.97, 1.38, 1.05, 1.14, 1.15, and 0.67. The base item discriminations $b_{i0}$ were chosen as −1.31, 1.44, −1.20, 0.10, 0.10, −0.74, 1.48, −0.61, 0.82, and −0.07. The same item parameters were also chosen in [65]. For item numbers as multiples of 10, we duplicated the item parameters of the 10 items accordingly. The group-specific item difficulties $b_{ig}$$(g=1,2)$ were simulated as

$$b_{i1}=b_{i0}-e_i/2\mathrm{and}{b}_{i2}=b_{i0}+e_i/2,$$

(39)

where $e_i$ is a random DIF effect. Note that $e_i=b_{i2}-b_{i1}$ parametrizes a uniform DIF effect as the difference in group-specific item difficulties. Group-specific item discriminations $a_{ig}$ $(g=1,2)$ were simulated as

$$a_{i1}=a_{i0}exp(-f_i/2)\mathrm{and}{a}_{i2}=a_{i0}exp(f_i/2),$$

(40)

where $f_i$ is another random DIF effect. The nonuniform DIF effect $f_i$ can be computed as the difference of logarithms of item discriminations (i.e., $f_i=log{a}_{i2}-log{a}_{i1}$). In the simulation, we assumed that $e_i$ and $f_i$ were uncorrelated. Both DIF effects had zero means and had standard deviation $\tau$ for $e_i$ and $0.3\times \tau$ for $f_i$. Two distributions of $e_i$ and $f_i$ were specified: a normal distribution or a scaled t distribution with three degrees of freedom (denoted as $t_3$). In this simulation study, we varied the DIF standard deviation for DIF effects for $e_i$ item difficulties as 0.25 and 0.50. According to the definition, the respective DIF standard deviation for DIF effects $f_i$ for (logarithmized) item discriminations were 0.075, and 0.15, respectively.

Five different linking methods were utilized to compute to estimate the mean difference $\widehat{\mathsf{\mu }}$ and the standard deviation $\widehat{\mathsf{\sigma }}$: mean–mean (MM) linking, mean–geometric–mean (MGM) linking, robust mean–geometric–mean (RGM) linking, symmetric Haebara (HAE) linking, and robust symmetric Haebara (RHAE) linking. The linking methods rely on estimated item discriminations ${\widehat{a}}_{ig}$ and item difficulties ${\widehat{b}}_{ig}$ ($i=1,\dots ,I$, $g=1,2$). In an infinite sample size, these identified item parameters are given as

$${\widehat{a}}_{i1}=a_{i1},{\widehat{b}}_{i1}=b_{i1},{\widehat{a}}_{i2}={\mathsf{\sigma }}^{-1}{a}_{i2},\mathrm{and}{\widehat{b}}_{i2}={\mathsf{\sigma }}^{-1}(b_{i2}-\mathsf{\mu }).$$

(41)

In each of the 2 (type of distribution) × 2 (DIF standard deviation $\tau$) × 3 (number of items I) = 12 cells of the simulation, 5000 replications were conducted. We computed bias and root mean square error (RMSE) for the estimated mean $\widehat{\mathsf{\mu }}$ and the estimated standard deviation $\widehat{\mathsf{\sigma }}$. A relative percentage RMSE was computed as the ratio of the RMSE values of a particular linking method and the chosen reference method of MGM linking. We also assessed the coverage rate for $\widehat{\mathsf{\mu }}$ and $\widehat{\mathsf{\sigma }}$ at the 95% confidence level based on the normal distribution as the percentage of events in which an estimated confidence interval contained the true value $\mathsf{\mu }=0.3$ or $\mathsf{\sigma }=1.2$, respectively.

The R software (Version 4.3.0; [86]) was used for the entire analysis in this simulation study. We wrote an R function linking_2groups_dich() that allows for the computation of estimates and their standard errors of any user-defined linking method. The function and replication material for this Simulation Study 1 can be found at https://osf.io/6bp3t (accessed on 29 April 2024).

5.2. Results

It turned out that all five linking methods resulted in approximately unbiased estimates for the mean $\widehat{\mathsf{\mu }}$ and the standard deviation $\widehat{\mathsf{\sigma }}$. Table 1 displays the relative RMSE for the estimates as a function of the DIF effect distribution type, the DIF standard deviation $\tau$, and the number of items I. The MM linking method performed similarly to the MGM method with respect to the mean $\mathsf{\mu }$ (i.e., the MM had relative RMSE values close to 100; that is, those of the MGM) but was slightly less efficient with respect to the standard deviation $\mathsf{\sigma }$, particularly for DIF effects that followed the $t_3$ distribution. The RMGM linking method had substantial efficiency losses for normally distributed DIF effects (i.e., it had RMSE values much larger than 100). On the other hand, the RMGM method provided large efficiency gains compared to the MGM method for the heavy-tailed $t_3$ distribution for $\mathsf{\mu }$ and $\mathsf{\sigma }$ (i.e., it had RMSE values much smaller than 100). Note that the RHAE method outperformed the HAE method for DIF effects with a $t_3$ distribution when analyzing the RMSE for the estimated mean $\widehat{\mathsf{\mu }}$. However, RHAE (or HAE) should not be preferred over RMGM (or MGM) with respect to the RMSE of the estimated standard deviation $\widehat{\mathsf{\sigma }}$, because it had RMSE values much larger than 100.

Table 2 reports the coverage rate for $\widehat{\mathsf{\mu }}$ and $\widehat{\mathsf{\sigma }}$ for the five linking methods. Coverage rates that are within the interval [92.5,97.5] indicate acceptable performance. Overall, the coverage rates based on M-estimation theory performed satisfactorily for at least 20 items, except for the RMGM linking method, which resulted in undercoverage for $\widehat{\mathsf{\mu }}$ (i.e., it had coverage rates much lower than 92.5). In this case, the estimated linking errors were, on average, too small compared to the standard deviation of the $\widehat{\mathsf{\mu }}$ estimates across replications in this simulation study. As expected, the coverage rates improved with an increasing number of items I.

6. Simulation Study 2: Assessing Total Error in Finite Sample Sizes

In the second Simulation Study 2, we investigated the statistical performance of the linking estimation methods in finite samples. In this case, there was uncertainty regarding the sampling of persons, thus resulting in standard error results and randomness in the DIF effects, which resulted in linking error results. Both sources of errors can be summarized as the total error.

6.1. Method

The item responses were simulated according to the 2PL model for a test with $I=10$, 20, or 40 items. The same item parameters as in Simulation Study 1 (see Section 5.1) were used. The factor variable $\mathsf{\theta }$ was assumed to be normally distributed in both groups. Like in Simulation Study 1, the mean and the standard deviation for the factor variable $\mathsf{\theta }$ in the first group were fixed at 0 and 1, respectively. The variable $\mathsf{\theta }$ in the second group had a mean $\mathsf{\mu }=0.3$ and a standard deviation $\mathsf{\sigma }=1.2$. We used the normal distribution and the scaled $t_3$ distribution for DIF effects and varied the DIF standard deviations $\tau$ as 0.25 and 0.5. Moreover, we simulated a condition of no DIF effects (i.e., $\tau =0$). The sample sizes $N=500$, 1000, and 2000 were chosen in order to mimic sample sizes that are typically in the applications of the 2PL model [8].

In contrast to Simulation Study 1, the item parameters of the 2PL model were separately estimated for the two groups in the first step using MML estimation. In the second step, the performance of the five linking methods MM, MGM, RMGM, HAE, and RHAE was studied. The estimated mean $\widehat{\mathsf{\mu }}$ and the estimated standard deviation $\widehat{\mathsf{\sigma }}$ for the five methods were compared regarding the bias, RMSE, and relative RMSE. As in Simulation Study 1, MGM linking was used as the reference method for computing the RMSE.

In total, 5000 replications were conducted in each of the 5 (type of distribution for DIF effects combined with DIF standard deviation $\mathsf{\sigma }$) × 3 (number of items I) × 3 (sample size N) = 45 cells of the simulation.

In the analysis, we computed the median of the linking error estimate $\mathrm{LE}$ based on (24) and the median of the bias-corrected linking error estimate ${\mathrm{LE}}_{\mathrm{bc}}$ based on (36). Moreover, we compared the coverage rates for the estimates $\widehat{\mathsf{\mu }}$ and $\widehat{\mathsf{\sigma }}$ based on the standard error $\mathrm{SE}$, as well as the (uncorrected) total error $\mathrm{TE}$, based on (31), and the bias-corrected total error ${\mathrm{TE}}_{\mathrm{bc}}$ was based on (37).

The R software [86] was used for the entire analysis in this simulation study. The 2PL model was estimated with the sirt::xxirt() function in the R package sirt [87]. As in Simulation Study 1, the R function linking_2groups_dich() was used for computing the estimates $\widehat{\mathsf{\mu }}$ and $\widehat{\mathsf{\sigma }}$ and their standard errors for the five linking methods. Replication material for this Simulation Study 2 can be found at https://osf.io/6bp3t (accessed on 29 April 2024).

6.2. Results

All five linking methods were approximately unbiased in all conditions of the simulation study. Table 3 presents the relative RMSE as a function of the different DIF distribution types, the DIF effect standard deviation $\tau$, the number of items I, and the sample size N. As expected from the literature, the HAE method was the most efficient linking method in the condition of no DIF (i.e., $\tau =0$). Across all conditions, the MM method had comparable performance to the MGM method regarding $\widehat{\mathsf{\mu }}$, but it was slightly more efficient for $\widehat{\mathsf{\sigma }}$. Efficiency gains of the RMGM method were only realized for the heavy-tailed $t_3$ distribution in large sample sizes. In these situations, the RHAE method outperformed the RMGM method for the estimated mean but not for the estimated standard deviation.

Table 4 presents the median of the estimated linking error $\mathrm{LE}$ and bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$ for the estimated mean $\widehat{\mathsf{\mu }}$. The results for the DIF effects with a $t_3$ distribution and a standard deviation of $\tau =0.5$ were omitted for space reasons. In Table 4, we have also reported the estimated linking error for an infinite sample size (i.e., $N=$ Inf) that was obtained from Simulation Study 1. It can be seen that the estimated linking errors $\mathrm{LE}$ and ${\mathrm{LE}}_{\mathrm{bc}}$ (almost always) converged to the $\mathrm{LE}$ for an infinite sample size with an increasing sample size.

It turned out that the estimated linking error $\mathrm{LE}$ was positively biased, while the bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$ was negatively biased (to a lesser extent). In particular, the median estimated linking error of 0.061 for the MM method for 10 items and a small sample size $N=500$ was substantially larger than the true value of 0 in the condition of no DIF (i.e., $\tau =0$). On the other hand, the median of the estimated bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$ was 0 in all situations in which no DIF was simulated in the item parameters. Overall, one could conclude that the bias in both linking error types can be reduced with an increasing sample size and an increasing number of items. For all linking methods except for the RMGM method, the uncorrected linking error $\mathrm{LE}$ had worse performance compared to the bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$. Hence, ${\mathrm{LE}}_{\mathrm{bc}}$ could be the preferred choice for a reported linking error.

Table 5 reports the median values of the estimated linking error $\mathrm{LE}$ and bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$ for the estimated standard deviation $\widehat{\mathsf{\sigma }}$. Overall, we observed a similar pattern of findings as in the case of the estimated mean $\widehat{\mathsf{\mu }}$. Again, the bias-corrected linking error estimates for the RMGM method were unsatisfactory. The bias in the uncorrected linking error estimates was slightly larger for $\widehat{\mathsf{\sigma }}$ than for $\widehat{\mathsf{\mu }}$. A simple idea might be to use the mean of both of the linking error estimates $\mathrm{LE}$ and ${\mathrm{LE}}_{\mathrm{bc}}$ as another linking error to improve the performance of the linking error estimate.

In Table 6, the coverage rate for the estimated mean $\widehat{\mathsf{\mu }}$ is displayed. In the no DIF condition $\tau =0$, the coverage rates based on the standard error $\mathrm{SE}$ performed satisfactorily. In the presence of DIF, the uncertainty in the estimated means of $\widehat{\mathsf{\mu }}$ was underestimated when the standard error was used when computing confidence intervals, thus resulting in substantial undercoverage. The confidence intervals based on the total error $\mathrm{TE}$ tended to have slightly increased coverage rates. In such situations, the coverage rates for the confidence intervals based on the bias-corrected linking error ${\mathrm{TE}}_{\mathrm{bc}}$ were slightly better. Generally, the RMGM linking method did not have adequate coverage rates in many situations.

Finally, Table 7 reports the coverage rates for the estimated standard deviation $\widehat{\mathsf{\sigma }}$. The bias-corrected total error outperformed the uncorrected total error regarding coverage rates. In many situations, the coverage rates based on the total error were too high. However, the RMGM method had substantial overcoverage in many conditions for confidence intervals based on $\mathrm{TE}$ and ${\mathrm{TE}}_{\mathrm{bc}}$, particularly for fewer items or smaller sample sizes.

7. Discussion

In this article, we simultaneously treated standard errors and linking errors for linking methods in the 2PL model. We proposed a bias-corrected linking error estimate, which delivers a bias-corrected total error estimate. This bias-corrected total error outperformed the usually employed total error that is given as the simple variances due to standard error and the usual uncorrected linking error. In a simulation study, it turned out that the confidence intervals for the linking parameters based on the bias-corrected total error outperformed those based on the usual total error regarding coverage rates. Moreover, the bias-corrected linking error estimate was less biased than the uncorrected linking error estimate.

As with any simulation study, our study had several limitations. First, our study only treated the 2PL model for dichotomous item responses. However, the performance of the linking estimators and their variance estimates could also be investigated for the simpler Rasch model for dichotomous item responses [16] or the generalized partial credit model for polytomous item responses [88]. Furthermore, the theory in this article could also be adapted to the chain linking of multiple groups [83,84]. In addition, the distribution types of the DIF effects in the simulation studies were restricted to the symmetric normal distribution and the t distributions with three degrees of freedom. Future research could focus on alternative and asymmetric distributions such as mixture, uniform, or discrete distributions. Moreover, the factor variable $\mathsf{\theta }$ was assumed to be normally distributed in both simulation studies. The 2PL model could also be estimated with non-normal $\mathsf{\theta }$ distributions [10,13], which could be investigated in future studies. Next, followup research could focus on linking with smaller sample sizes, such as $N=250$, as well as the case of unbalanced group sizes. Furthermore, we only employed 10, 20, or 40 items in the two simulation studies. Future research could also investigate a larger number of items. We do not think that linking should be conducted with an even smaller number of items, such as $I=5$ items, because the group comparisons will likely become unstable in the presence of DIF, and the representativity of the link items might be questioned (but, see [89]). Also, the extent of nonuniform DIF was not independently manipulated from the extent of uniform DIF in the two simulation studies. Finally, the performance of our proposed error estimates could also be applied to mis-specified IRT models. For example, the 2PL model could be employed for linking if the item response data would be generated from the logistic positive exponential model [90,91] or the monotonic polynomial IRT model [92,93]. All of these limitations could be addressed in extensive future research.

As a final side note, I would like to add that a comparison of two groups regarding the distribution of the factor variable $\mathsf{\theta }$ could also be conducted using concurrent calibration by assuming invariant (i.e., the same) item parameters across the groups. Some researchers argue that linking uncertainty is reduced by assuming invariant item parameters (see [94,95]). I think that this belief is unjustified. The fact that there could be variability due to item selection does not disappear, because the variability in the model parameters is not represented in the statistical model. The computation of the linking errors under the assumption of invariant item parameters in the statistical model has been carved out in Ref. [96].