Linking Error Estimation in Stocking–Lord Linking

Abstract

Stocking–Lord (SL) linking is a widely used linking method based on item response theory (IRT). This article examines the variability in SL linking parameter estimates within the two-parameter logistic (2PL) model. The uncertainty in SL linking arises from the sampling variability (standard error) and item selection (linking error), which can induce variability due to random differential item functioning (DIF). Three linking error estimation approaches are compared in this paper: the conventional jackknife linking error method, a newly developed approximate jackknife linking error method, and a Taylor approximation-based estimate. Simulation studies showed that the approximate jackknife method closely aligns with the traditional jackknife linking error method and outperforms the linking error estimation approach based on Taylor approximation. The adequacy of coverage rates for SL linking parameter estimates was also assessed using estimates of the total error. Results from a simulation study demonstrate that the bias-corrected total error provides superior coverage rates compared to both the conventional total error and the standard error, which does not account for item-related uncertainty due to random DIF.

1. Introduction

Item response theory (IRT) models [1] are multivariate statistical models for analyzing multivariate binary (i.e., dichotomous) random variables. Let $\mathit{X}=(X_1,\dots ,X_I)$ denote a vector of I binary random variables $X_i\in \{0,1\}$, often referred to as items or (scored) item responses. A unidimensional IRT model [1] is a statistical model for the probability distribution $\mathsf{P}(\mathit{X}=\mathit{x})$

$${\displaystyle \mathsf{P}(\mathit{X}=\mathit{x};\mathit{\delta },\mathit{\gamma })=\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 ;\mu ,\sigma )\mathrm{d}\theta }$$

(1)

for $\mathit{x}=(x_1,\dots ,x_I)\in {\{0,1\}}^I$. Here, $\varphi$ denotes the density of the normal distribution, depending on the mean $\mu$ and the standard deviation (SD) $\sigma$; the distribution parameters of the latent variable $\theta$, often referred to as a trait or ability variable, are contained in the vector $\mathit{\delta }=(\mu ,\sigma )$. The vector $\mathit{\gamma }=({\mathit{\gamma }}_1,\dots ,{\mathit{\gamma }}_I)$ includes the item parameters for the item response functions (IRFs) $P_i(\theta ;{\mathit{\gamma }}_i)=\mathsf{P}(X_i=1|\theta )$ (for $i=1,\dots ,I$). The IRF of the two-parameter logistic (2PL) model [2] is provided by

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

(2)

where $a_i$ is the item discrimination and $b_i$ is the item difficulty $b_i$, with $\Psi \left(x\right)={(1+exp(-x))}^{-1}$ as the logistic distribution function. For independently and identically distributed observations ${\mathit{x}}_1,\dots ,{\mathit{x}}_N$ across N persons from the distribution of the random variable $\mathit{X}$, the unknown model parameters in the IRT model (1) can be consistently estimated through marginal maximum likelihood estimation (MML; [3]). It is worth noting that identification constraints in the IRT model are often required, as item and distribution parameters would otherwise not be statistically identified.

IRT models are commonly used to compare the performance of two groups on a test, examining differences in the factor variable $\theta$ within the IRT model (1). This article focuses on linking methods [4] in the 2PL model. In the first step of the linking method, the 2PL model can be estimated separately for each group, which accommodates the possibility of items functioning differently across groups, known as differential item functioning (DIF; [5,6]). In the second step, differences in item parameters are utilized to calculate group differences in $\theta$ through a linking method [4,7].

Random DIF [8] introduces additional variability in the estimated mean $\mu$ and the SD $\sigma$ of the second group in a linking method, making the distribution parameters $\mathit{\delta }=(\mu ,\sigma )$ sensitive to item selection even with an infinite number of persons. This variability, known as the linking error, can be quantified to assess the stability of linking parameter estimates across selected items [9,10,11,12,13].

In a separate estimation of the 2PL model, as the first step of the linking procedure, group-specific item discriminations ${\widehat{a}}_{ig}$ and item difficulties ${\widehat{b}}_{ig}$ are obtained ($i=1,\dots ,I$, $g=1,2$). We arrange all estimated item parameters into vectors ${\widehat{\mathit{a}}}_g$ and ${\widehat{\mathit{b}}}_g$ for the two groups $g=1,2$. As presented in [14], the Stocking–Lord (SL) linking method utilizes the linking function H, which is defined as

$$H(\mathit{\delta },{\widehat{\mathit{a}}}_1,{\widehat{\mathit{b}}}_1,{\widehat{\mathit{a}}}_2,{\widehat{\mathit{b}}}_2)=\sum _{t=1}^T{\omega }_t{\left[\sum _{i=1}^I\Psi \left({\widehat{a}}_{i1}(\sigma {\theta }_t+\mu -{\widehat{b}}_{i1})\right)-\sum _{i=1}^I\Psi \left({\widehat{a}}_{i2}({\theta }_t-{\widehat{b}}_{i2})\right)\right]}^2.$$

(3)

for $\mathit{\delta }=(\mu ,\sigma )$. The weights ${\omega }_t$ are predetermined, and an equidistant grid for the ability variable $\theta$ is established spanning from ${\theta }_1$ to ${\theta }_T$ for ${\theta }_t$ with $t=1,\dots ,T$. The weights ${\omega }_t$ may either be uniform (equal to one) or proportional to a discretized density of a normal distribution with an SD larger than 1 (e.g., an SD of 2). The linking function (3) is referred to as asymmetric Stocking–Lord linking because it aligns the test characteristic function (TCF) from the first group with that of the second group.

The distribution parameters $\mu$ and $\sigma$ can be estimated by minimizing the linking function H

$$\widehat{\mathit{\delta }}=(\widehat{\mu },\widehat{\sigma })=\underset{\mathit{\delta }}{arg\; min}H(\mathit{\delta },{\widehat{\mathit{a}}}_1,{\widehat{\mathit{b}}}_1,{\widehat{\mathit{a}}}_2,{\widehat{\mathit{b}}}_2).$$

(4)

The estimating equations $\mu$ and $\sigma$ are derived by calculating the partial derivatives of H with respect to $\mu$ and $\sigma$. The estimate $\widehat{\mathit{\delta }}$ matches the true parameter $\mathit{\delta }$ in the absence of sampling errors and random DIF. This property validates the use of the linking function H in SL linking.

This article is devoted to the computation of the linking error in SL linking. To th best of our knowledge, there have been no attempts in the literature to derive a linking error estimation approach for SL linking. This article fills this gap by deriving an estimate based on M-estimation theory [15]. Moreover, linking error estimation can alternatively be carried out by applying the jackknife approach to items [10]. The jackknife linking error approach can be computationally intensive for the iterative SL linking method for a moderate to large number of items. For this reason, an approximation of the jackknife linking error is proposed in this article. Finally, bias-corrected linking error estimates are derived which have improved behavior in small sample sizes. The adequacy of the different error estimation approaches is compared through empirical coverage rates of $\widehat{\mu }$ and $\widehat{\sigma }$ obtained from SL linking.

The rest of this article is organized as follows: Section 2 discusses the estimation of the standard error, linking error, and total error for SL linking; in Section 3, we present findings from a simulation study conducted with an infinite number of persons; Section 4 reports results from a simulation study that examines finite sample sizes; finally, the article closes with a discussion in Section 5 and conclusions in Section 6.

2. Assessing the Standard Error, Linking Error, and Total Error in Stocking–Lord Linking

The SL linking method can generally be formulated as the estimation problem

$$\widehat{\mathit{\delta }}=\underset{\mathit{\delta }}{arg\; min}H(\mathit{\delta },\widehat{\mathit{\gamma }}),$$

(5)

where H is the linking function used in SL linking (see (3)), $\mathit{\delta }=(\mu ,\sigma )$ is the vector of linking parameters $\mu$ and $\sigma$ that represent the mean and the SD of the second group, the vector $\widehat{\mathit{\delta }}=(\widehat{\mu },\widehat{\sigma })$ contains the estimated linking parameters obtained from SL linking, and the vector $\widehat{\mathit{\gamma }}$ contains all estimated item parameters from separate estimation of the 2PL model in the first step of the linking method. By using partial derivatives ${\mathit{H}}_{\mathit{\delta }}=(\partial H)/(\partial \mathit{\delta })$, the parameter estimate $\widehat{\mathit{\delta }}$ from (5) fulfills the nonlinear equation

$${\mathit{H}}_{\mathit{\delta }}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=\mathbf{0}.$$

(6)

For a general linking function H, the true parameter vector $\delta$ must satisfy the population estimation equation ${\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\mathit{\gamma })=0$, where the item parameters $\mathit{\gamma }$ are assumed to be invariant across the two groups and free from sampling errors.

This article addresses estimation of the uncertainty of the estimated mean $\widehat{\mu }$ and SD $\widehat{\sigma }$ in the parameter vector $\mathit{\delta }$. This uncertainty arises from two main sources during group comparisons: the sampling of persons (as represented in the scenario where $N\to \infty$), and the selection or modeling of items (in the context of $I\to \infty$). The variability in linking parameter estimates due to the selection of items is a consequence of the presence of random DIF.

Much of the existing literature on linking addresses uncertainty in $\widehat{\mathit{\gamma }}$ by calculating the standard error (SE) of $\widehat{\mathit{\delta }}$ for a fixed number of items [9,16,17,18]. In this situation, the variability in $\widehat{\mathit{\gamma }}$ is attributed to the sampling variability present in the estimated item parameters ${\widehat{\mathit{\gamma }}}_i$. Variance estimation for $\widehat{\mathit{\delta }}$ operates under the assumption that $N\to \infty$.

Estimated item parameters ${\widehat{\mathit{\gamma }}}_i$ for item i correspond to a population analog ${\mathit{\gamma }}_i$ when considering an infinite sample size of persons. However, in the presence of model error (e.g., random DIF), the estimated linking parameter $\widehat{\mathit{\delta }}$ is influenced by the chosen set of items even with an infinite sample size of persons N. This variability stemming from item selection due to DIF is referred to as the linking error (LE, [10,19]). Variance estimation operates under the assumption of $I\to \infty$.

It has been shown that linking error estimates also contain portions of sampling error [20]. In detail, the variance estimate ${\mathit{V}}_{\mathrm{LE}}$ for the linking error can be decomposed into

$${\mathit{V}}_{\mathrm{LE}}={\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}}}+f\left(\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)\right),$$

(7)

where $f\left(\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)\right)$ is the expected value of the variance contribution due to the sampling variance $\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)$ in estimated item parameters $\widehat{\mathit{\gamma }}$. If an estimate $\widehat{f}\left(\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)\right)$ for this term were available, a bias-corrected variance estimate ${\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}}}$ for the linking error could be obtained from (7) as

$${\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}}}={\mathit{V}}_{\mathrm{LE}}-\widehat{f}\left(\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)\right).$$

(8)

This article proposes a bias-corrected linking error estimate for SL linking.

The total error (TE) comprises both sources of uncertainty, that is, the standard errorcaused by to randomness due to persons and the linking error caused by to randomness (or random DIF) due to items [9,11,19,21,22]. The conventional estimate of the variance matrix for the total error is provided by

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

(9)

A bias-corrected variance estimate of the total error is provided by

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

(10)

In the next subsections, we discuss estimation of the SE, LE, and TE for SL linking.

2.1. Standard Error Estimation

First, we review the SE estimation of the linking parameter estimate $\widehat{\mathit{\delta }}$ in SL linking. The treatment is based on the delta method, and can be found in the literature [16,17,23,24]. Let ${\mathit{V}}_{\widehat{\mathit{\gamma }}}=\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)$ be the variance matrix of the vector of estimated item parameters $\widehat{\mathit{\gamma }}$. A Taylor approximation of the partial derivative ${\mathit{H}}_{\mathit{\delta }}$ of the linking function H provides

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

(11)

Because it holds ${\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\mathit{\gamma })=\mathbf{0}$, we obtain

$$\widehat{\mathit{\delta }}-\mathit{\delta }={\mathit{A}}_{\mathrm{SE}}(\widehat{\mathit{\gamma }}-\mathit{\gamma }),\mathrm{where}{\mathit{A}}_{\mathrm{SE}}={\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\mathit{\delta },\mathit{\gamma })}^{-1}{\mathit{H}}_{\mathit{\delta }\mathit{\gamma }}(\mathit{\delta },\mathit{\gamma }).$$

(12)

Hence, the variance of $\widehat{\mathit{\delta }}$ due to sampling persons can be written as

$${\mathit{V}}_{\mathrm{SE}}=\mathsf{Var}\left(\widehat{\mathit{\delta }}\right)={\mathit{A}}_{\mathrm{SE}}{\mathit{V}}_{\widehat{\mathit{\gamma }}}{\mathit{A}}_{\mathrm{SE}}^{\top }.$$

(13)

An estimate of the variance is obtained by

$${\mathit{V}}_{\mathrm{SE}}=\mathsf{Var}\left(\widehat{\mathit{\delta }}\right)={\widehat{\mathit{A}}}_{\mathrm{SE}}{\mathit{V}}_{\widehat{\mathit{\gamma }}}{\widehat{\mathit{A}}}_{\mathrm{SE}}^{\top },\mathrm{where}{\widehat{\mathit{A}}}_{\mathrm{SE}}={\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\widehat{\mathit{\delta }},\mathit{\gamma })}^{-1}{\mathit{H}}_{\mathit{\delta }\mathit{\gamma }}(\widehat{\mathit{\delta }},\mathit{\gamma }).$$

(14)

Standard error estimates for $\widehat{\mu }$ or $\widehat{\sigma }$ can be obtained by taking the square root of the corresponding elements in ${\mathit{V}}_{\mathrm{SE}}$.

The variance matrix ${\mathit{V}}_{\mathrm{SE}}$ correctly assesses the variability in $\widehat{\mathit{\delta }}$ if no random DIF occurs. In the presence of random DIF, linking errors must also be considered in order to correctly reflect the uncertainty in the linking parameter estimate. The estimation of these errors is discussed in the following sections.

2.2. Linking Error Estimation Based on Taylor Approximation

The assessment of linking errors for linking methods was investigated using M-estimation theory in [19,20]. This theory has been developed for linking functions that are additive regarding the dependence on items, that is,

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

(15)

However, SL linking does not allow such an additive decomposition. To this end, we derive a variance estimate of the linking error in SL linking that relies on Taylor approximation. We use the notation

$$P_i(\theta ;{\mathit{\gamma }}_{ig})=\Psi \left(a_{ig}(\theta -b_{ig})\right).$$

(16)

The estimating function in SL linking can be written as

$${\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})=\left(\begin{array}{c}{H}_{\mu }(\mathit{\delta },\widehat{\mathit{\gamma }})\\ H_{\sigma }(\mathit{\delta },\widehat{\mathit{\gamma }})\end{array}\right)=\sum _{t=1}^T{\omega }_t\sum _{i=1}^I{Z}_{it}(\mathit{\delta },\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})=\mathbf{0},$$

(17)

where we use the abbreviations

$$Z_{it}(\mathit{\delta },\widehat{\mathit{\gamma }})=P_i(\sigma {\theta }_t+\mu ;{\widehat{\mathit{\gamma }}}_{i1})-P_i({\theta }_t;{\widehat{\mathit{\gamma }}}_{i2})\mathrm{and}$$

(18)

$${\mathit{D}}_{I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})=\left(\begin{array}{c}{D}_{\mu ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})\\ D_{\sigma ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{I}\sum _{i=1}^I{a}_i{P}_i^{\prime }(\sigma {\theta }_t+\mu ;{\widehat{\mathit{\gamma }}}_{i1})}\\ {\displaystyle \frac{1}{I}\sum _{i=1}^I{a}_i{\theta }_t{P}_i^{\prime }(\sigma {\theta }_t+\mu ;{\widehat{\mathit{\gamma }}}_{i1})}\end{array}\right).$$

(19)

Now, we define the stochastic limits

$${\mathit{D}}_t(\mathit{\delta },\widehat{\mathit{\gamma }})=\left(\begin{array}{c}{D}_{\mu ,t}(\mathit{\delta },\widehat{\mathit{\gamma }})\\ D_{\sigma ,t}(\mathit{\delta },\widehat{\mathit{\gamma }})\end{array}\right)=\underset{I\to \infty }{plim}{\mathit{D}}_{I,t}(\mathit{\delta },\widehat{\mathit{\gamma }}),\mathrm{where}$$

(20)

$$D_{\mu ,t}(\mathit{\delta },\widehat{\mathit{\gamma }})=\underset{I\to \infty }{plim}{D}_{\mu ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})\mathrm{and}{D}_{\sigma ,t}(\mathit{\delta },\widehat{\mathit{\gamma }})=\underset{I\to \infty }{plim}{D}_{\sigma ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }}).$$

(21)

We utilize a Taylor expansion of ${\mathit{H}}_{\mathit{\delta }}$ around $\mathit{\delta }$ and replace ${\mathit{D}}_{I,t}$ with ${\mathit{D}}_t$ in (17), yielding

$${\mathit{H}}_{\mathit{\delta }}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})={\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})+{\mathit{H}}_{\mathit{\delta }\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})(\widehat{\mathit{\delta }}-\mathit{\delta })=\mathbf{0},\mathrm{where}$$

(22)

$${\mathit{H}}_{\mathit{\delta }\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})=I\sum _{t=1}^T{\omega }_t\left(\begin{array}{cc}{D}_{\mu ,I,t}^2(\mathit{\delta },\widehat{\mathit{\gamma }})& D_{\mu ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})D_{\sigma ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})\\ D_{\mu ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})D_{\sigma ,I,t}(\mathit{\delta },\widehat{\mathit{\gamma }})& D_{\sigma ,I,t}^2(\mathit{\delta },\widehat{\mathit{\gamma }})\end{array}\right).$$

(23)

From (22), we obtain

$$\widehat{\mathit{\delta }}-\mathit{\delta }={\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\mathit{\delta },\widehat{\mathit{\gamma }})}^{-1}{\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }}).$$

(24)

The linking error can be derived from the variance matrix corresponding to the linear Taylor expansion (24) as

$$\mathsf{Var}\left(\widehat{\mathit{\delta }}\right)={\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\mathit{\delta },\widehat{\mathit{\gamma }})}^{-1}\mathsf{Var}\left({\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})\right){\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\mathit{\delta },\widehat{\mathit{\gamma }})}^{-\top }.$$

(25)

Note that we can rewrite (17) approximately as

$${\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})=\sum _{i=1}^I\sum _{t=1}^T{\omega }_t{Z}_{it}(\mathit{\delta },\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\mathit{\delta },\widehat{\mathit{\gamma }}),$$

(26)

and $\mathsf{Var}\left({\mathit{H}}_{\mathit{\delta }}(\mathit{\delta },\widehat{\mathit{\gamma }})\right)$ can be evaluated as the empirical variance of the random variables

$${\mathit{V}}_i(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=\sum _{t=1}^T{\omega }_t{Z}_{it}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})\mathrm{for}i=1,\dots ,I.$$

(27)

Hence, we finally obtain the variance matrix ${\mathit{V}}_{\mathrm{LE},\mathrm{TAY}}$ due to the linking error as

$${\mathit{V}}_{\mathrm{LE},\mathrm{TAY}}={\displaystyle \frac{I}{I-1}}{\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})}^{-1}\mathsf{Var}\left({\mathit{H}}_{\mathit{\delta }}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})\right){\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})}^{-\top }.$$

(28)

Note that the finite sampling correction factor $I/(I-1)>1$ is included in (28) to improve the parameter coverage performance of confidence interval estimates in finite samples (see [19,20,25]). The linking error estimates can be obtained as the square root of the corresponding diagonal entries in ${\mathit{V}}_{\mathrm{LE},\mathrm{TAY}}$.

2.3. Jackknife Linking Error

The jackknife method as a general method for estimating the linking error has been discussed in [10,11]. The main advantage of the jackknife linking error is that it is generally applicable to any linking method, even in concurrent calibration or fixed item parameter calibration [26]. The idea behind the application of this technique to linking is to repeatedly remove the ith item from the linking procedure and obtain the parameter estimate ${\widehat{\mathit{\delta }}}_{(-i)}$. This estimate differs slightly from the linking parameter estimate $\widehat{\mathit{\delta }}$ that involves all items. The variance estimate of the linking error based on the jackknife technique is defined as

$${\mathit{V}}_{\mathrm{LE},\mathrm{JK}}={\displaystyle \frac{I}{I-1}}\sum _{i=1}^I\left({\widehat{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}\right){\left({\widehat{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}\right)}^{\top }.$$

(29)

Note that the factor $I/(I-1)>1$ in (29) differs from the originally proposed factor $(I-1)/I<1$ in [10,27]. Here, we chose the former one because it yielded more satisfactory coverage rates for linking parameter estimates. Notably, the relevance of this factor diminishes with an increasing number of items I.

The variance estimate ${\mathit{V}}_{\mathrm{LE},\mathrm{JK}}$ contains additional variability due to sampling errors. This variance portion can be removed by the double jackknife method [11,21]. An alternative and less computationally demanding method for removing the bias contribution for the jackknife linking error is presented in Section 2.5.

2.4. Approximate Jackknife Linking Error

The jackknife linking error requires additional computation of the linking method I times (i.e., the number of items). Depending on the complexity of the linking method (e.g., if the estimation needs iterations), the jackknife estimate might be computationally cumbersome. To avoid this, the jackknife parameter estimates can be approximated by closed-form updates that involve only one iteration. We use the notation from Section 2.2 and define

$$Z_{it}\left(\widehat{\mathit{\delta }}\right)=P_i(\widehat{\sigma }{\theta }_t+\widehat{\mu };{\widehat{\mathit{\gamma }}}_{i1})-P_i({\theta }_t;{\widehat{\mathit{\gamma }}}_{i2}).$$

(30)

The linking parameter estimate $\widehat{\mathit{\delta }}$ based on all items is based on the estimating equation

$${\mathit{H}}_{\mathit{\delta }}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=\sum _{t=1}^T{\omega }_t\sum _{i=1}^I{Z}_{it}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=\mathbf{0}.$$

(31)

The parameter estimate ${\tilde{\mathit{\delta }}}_{(-i)}$ by removing item i approximately solves the equation

$$\sum _{t=1}^T{\omega }_t\sum _{{\displaystyle \genfrac{}{}{0pt}{}{j=1}{j\ne i}}}^I{Z}_{jt}({\tilde{\mathit{\delta }}}_{(-i)},\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=\mathbf{0}.$$

(32)

In (32), we assume that ${\mathit{D}}_{I,t}$ is a constant that does not depend on $\mathit{\delta }$. Next, we apply a Taylor expansion of (32) around $\widehat{\mathit{\delta }}$ and use the identity due to (31)

$$\sum _{t=1}^T{\omega }_t\sum _{{\displaystyle \genfrac{}{}{0pt}{}{j=1}{j\ne i}}}^I{Z}_{jt}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=-\sum _{t=1}^T{\omega }_t{Z}_{it}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }}).$$

(33)

The application of the Taylor formula applied to (32) yields

$$\mathbf{0}=-{\mathit{c}}_i+{\mathit{B}}_i({\tilde{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}),\mathrm{where}$$

(34)

$${\mathit{c}}_i=\sum _{t=1}^T{\omega }_t{Z}_{it}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }}){\mathit{D}}_{I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})\mathrm{and}$$

(35)

$${\mathit{B}}_i=\sum _{t=1}^T{\omega }_t\left[\sum _{{\displaystyle \genfrac{}{}{0pt}{}{j=1}{j\ne i}}}^I{a}_{j1}{P}_j^{\prime }(\widehat{\sigma }{\theta }_t+\widehat{\mu };{\widehat{\mathit{\gamma }}}_{j1})\right]\left(\begin{array}{c}{D}_{\mu ,I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})\\ D_{\sigma ,I,t}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})\end{array}\right)\left(\begin{array}{cc}1& {\theta }_t\end{array}\right).$$

(36)

Finally, we obtain

$${\tilde{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}={\mathit{B}}_i^{-1}{\mathit{c}}_i.$$

(37)

In this way, we obtain the variance estimate of the approximate jackknife procedure as

$${\mathit{V}}_{\mathrm{LE},\mathrm{AJK}}={\displaystyle \frac{I}{I-1}}\sum _{i=1}^I\left({\tilde{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}\right){\left({\tilde{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}\right)}^{\top }={\displaystyle \frac{I}{I-1}}\sum _{i=1}^I{\mathit{B}}_i^{-1}{\mathit{c}}_i{\mathit{c}}_i^{\top }{\mathit{B}}_i^{-\top }.$$

(38)

Because the derivatives of the linking function of H are often easy to obtain, the approximate jackknife method for estimating the linking error is computationally much cheaper than exact jackknife linking error estimation.

2.5. Bias-Corrected Approximate Jackknife Linking Error

The jackknife and approximate jackknife approaches to linking error estimation are also affected by sampling errors due to the sampling of persons, which positively biases these variance estimates. In the (approximate) jackknife technique, the expected value of the square of the deviations ${\widehat{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}$ seeks to estimate the square of the deviation ${\mathit{\delta }}_{(-i)}-\mathit{\delta }$. We observe that $\widehat{\mathit{\delta }}$ and ${\widehat{\mathit{\delta }}}_{(-i)}$ solve the nonlinear equations

$${\mathit{H}}_{\mathit{\delta }}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})=\mathbf{0}\mathrm{and}{\mathit{H}}_{\mathit{\delta },(-i)}({\widehat{\mathit{\delta }}}_{(-i)},\widehat{\mathit{\gamma }})=\mathbf{0},$$

(39)

where ${\mathit{H}}_{\mathit{\delta },(-i)}$ denotes the vector of partial derivatives from those linking function in which the ith item is removed. We can use Taylor expansion of ${\mathit{H}}_{\mathit{\delta }}$ and ${\mathit{H}}_{\mathit{\delta },(-i)}$, yielding

$$\widehat{\mathit{\delta }}-\mathit{\delta }=\mathit{U}(\widehat{\mathit{\gamma }}-\mathit{\gamma })\mathrm{with}\mathit{U}={\mathit{H}}_{\mathit{\delta }\mathit{\delta }}{(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})}^{-1}{\mathit{H}}_{\mathit{\delta }\mathit{\gamma }}(\widehat{\mathit{\delta }},\widehat{\mathit{\gamma }})\mathrm{and}$$

(40)

$${\widehat{\mathit{\delta }}}_{(-i)}-{\mathit{\delta }}_{(-i)}={\mathit{U}}_{(-i)}(\widehat{\mathit{\gamma }}-\mathit{\gamma })\mathrm{with}{\mathit{U}}_{(-i)}={\mathit{H}}_{\mathit{\delta }\mathit{\delta },(-i)}{({\widehat{\mathit{\delta }}}_{(-i)},\widehat{\mathit{\gamma }})}^{-1}{\mathit{H}}_{\mathit{\delta }\mathit{\gamma },(-i)}({\widehat{\mathit{\delta }}}_{(-i)},\widehat{\mathit{\gamma }}).$$

(41)

Then, from Equations (40) and (41) we can derive

$$\left[{\widehat{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}\right]-\left[{\mathit{\delta }}_{(-i)}-\mathit{\delta }\right]=\left({\mathit{U}}_{(-i)}-\mathit{U}\right)(\widehat{\mathit{\gamma }}-\mathit{\gamma }).$$

(42)

Equation (42) provides a hint for bias correction in the form of the variance of the linking errors induced by the sampling variability associated with $\widehat{\mathit{\gamma }}$. Consequently, the variance bias can be expressed as

$${\mathcal{V}}_{\mathrm{Bias}}={\displaystyle \frac{I}{I-1}}\sum _{i=1}^I\mathsf{Var}\left({\widehat{\mathit{\delta }}}_{(-i)}-\widehat{\mathit{\delta }}\right)={\displaystyle \frac{I}{I-1}}\sum _{i=1}^I\left({\mathit{U}}_{(-i)}-\mathit{U}\right){\mathit{V}}_{\widehat{\mathit{\gamma }}}{\left({\mathit{U}}_{(-i)}-\mathit{U}\right)}^{\top },$$

(43)

where ${\mathit{V}}_{\widehat{\mathit{\gamma }}}=\mathsf{Var}\left(\widehat{\mathit{\gamma }}\right)$. Equation (43) can be utilized to compute the bias-corrected jackknife estimates

$${\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}},\mathrm{JK}}={\mathit{V}}_{\mathrm{LE},\mathrm{JK}}-{\mathcal{V}}_{\mathrm{Bias}}\mathrm{and}{\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}},\mathrm{AJK}}={\mathit{V}}_{\mathrm{LE},\mathrm{AJK}}-{\mathcal{V}}_{\mathrm{Bias}}.$$

(44)

The corresponding linking error estimates can also be obtained by taking the square root of the corresponding diagonal elements in the variance matrices ${\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}},\mathrm{JK}}$ or ${\mathit{V}}_{{\mathrm{LE}}_{\mathrm{bc}},\mathrm{AJK}}$. In the case of negative variance estimates, the linking error estimate is defined to be zero. Notably, Equation (44) allows an accurate estimation of the bias-corrected linking error to be obtained.

3. Simulation Study 1: Infinite Sample Size

3.1. Method

In Simulation Study 1, only item parameters and no item responses were generated for each replication, as the focus was on a scenario with an infinite sample size. Hence, uncertainty was only caused by DIF in item parameters, not by sampling errors due to the sampling (or choice) of persons. The 2PL model for two groups served as the IRF in the data-generating model (DGM). To ensure identification, the mean and SD of the factor variable $\theta$ for the first group were fixed at 0 and 1, respectively. For the second group, the mean $\mu$ and the SD $\sigma$ were used to represent group differences, with the fixed values of $\mu =0.3$ and $\sigma =1.2$ applied throughout all simulation conditions in the DGM.

We carried out the simulation study for $I=10$, 20, 40, and 80 items. Group-specific parameters $a_{ig}$ and $b_{ig}$ for each item for groups $g=1,2$ were based on fixed base item parameters and random DIF effects that were newly simulated in each replication of the simulation. The item parameters used in this simulation study align with those in [19,20]. For the case of $I=10$ items, base item discrimination values $a_{i0}$ were set at 0.73, 1.25, 1.20, 1.47, 0.97, 1.38, 1.05, 1.14, 1.15, and 0.67, resulting in a mean $M=1.101$ and an SD of 0.257 for the item discriminations. The base item difficulties $b_{i0}$ were set to −1.31, 1.44, −1.20, 0.10, 0.10, −0.74, 1.48, −0.61, 0.82, and −0.07, yielding a mean item difficulty of $M=0.001$ and an SD of 1.002. For item numbers that were multiples of 10, the parameters for the ten items were repeated accordingly. Group-specific item difficulties $b_{ig}$$(g=1,2)$ were then simulated as follows:

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

(45)

where $e_i$ represents a random uniform DIF effect. Note that $e_i=b_{i2}-b_{i1}$ defines the uniform DIF effect as the difference in item difficulties between the two groups. Group-specific item discriminations $a_{ig}$ $(g=1,2)$ were generated as

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

(46)

where $f_i$ introduces a second random DIF effect. This nonuniform DIF effect $f_i$ is defined as the difference in the logarithms of item discriminations, that is, $f_i=log{a}_{i2}-log{a}_{i1}$. In the simulation, the DIF effects $e_i$ and $f_i$ were modeled as uncorrelated, with both having a mean of zero. Both DIF effects were normally distributed. The SD of the DIF effects was set to $\tau$ for $e_i$ and $0.3·\tau$ for $f_i$. In this simulation study, the DIF SDs $\tau$ for $e_i$ (DIF in item difficulties) were simulated with values of 0.20, 0.40, and 0.60, resulting in respective DIF SDs for $f_i$ (DIF in log-transformed item discriminations) of 0.06, 0.12, and 0.18.

The identified item parameters in the infinite-size sample were computed as

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

(47)

SL linking was applied to these obtained item parameters. We compared three methods to assess the linking error. First, we applied the jackknife (JK) method as described in Section 2.3. Second, the approximate jackknife (AJK) method (see Section 2.4) was utilized. Third, the linking error was computed based on Taylor approximation (TAY), as described in Section 2.2.

In each of the 3 (DIF SD $\tau$) × 4 (number of items I) = 12 cells of the simulation, 4000 replications were performed. The adequacy of the alternative linking error estimation methods was assessed for the estimated mean $\widehat{\mu }$ and estimated SD $\widehat{\sigma }$. We evaluated the coverage rate for $\widehat{\mu }$ and $\widehat{\sigma }$ at a 95% confidence level, using the normal distribution to calculate the percentage of instances in which the confidence interval included the true values of $\mu =0.3$ or $\sigma =1.2$, respectively.

All analyses for this simulation study were conducted using R statistical software (Version 4.4.1; [28]). Dedicated R functions were written to carry out this simulation study. These functions, along with replication materials for the study, are available at https://osf.io/4dr2c (accessed on 28 October 2024).

3.2. Results

Table 1. Simulation Study 1: Coverage rates at the 95% confidence level for the mean $\mu$ and the standard deviation $\sigma$ in an infinite sample size as a function of the DIF SD $\tau$ and the number of items I.

		$\mathit{\tau }=\mathbf{0.2}$			$\mathit{\tau }=\mathbf{0.4}$			$\mathit{\tau }=\mathbf{0.6}$
Par	$\mathit{I}$	JK	AJK	TAY	JK	AJK	TAY	JK	AJK	TAY
$\mu$	10	93.1	93.1	90.1	94.1	94.2	91.5	93.5	93.5	91.0
	20	94.6	94.6	93.1	94.2	94.2	92.9	94.4	94.4	93.1
	40	94.6	94.6	93.8	94.9	94.9	94.3	94.8	94.8	94.0
	80	94.8	94.8	94.4	95.4	95.4	94.9	94.7	94.7	94.5
$\sigma$	10	95.6	95.5	91.6	95.0	94.8	91.2	94.6	94.6	89.9
	20	95.1	95.0	92.9	95.1	95.0	93.1	95.1	95.0	92.9
	40	94.6	94.5	93.3	94.9	94.9	94.0	95.0	95.0	93.9
	80	95.5	95.5	94.9	95.0	95.0	94.2	94.5	94.5	94.0

Note. Par = parameter; JK = linking error based on jackknife; AJK = linking error based on approximate jackknife; TAY = linking error based on Taylor approximation; Coverage rates smaller than 92.5 and larger than 97.5 are printed in bold font.

displays the coverage rates at the 95% confidence level for the mean $\mu$ and SD $\sigma$ in an infinite sample size across different SDs $\tau$ of the random DIF and different number of items I. For both $\mu$ and $\sigma$, the JK and AJK methods yielded nearly identical coverage rates, consistently around the target level of 95%; the largest difference between JK and AJK in the coverage rates was 0.2% for $\sigma$ in the condition $\tau =0.4$ and $I=10$. The TAY method showed lower coverage rates in conditions of higher DIF (e.g., $\tau =0.4$ or 0.6), indicating inferior performance of this method in these simulation conditions. Coverage rates for the TAY method decreased with increasing DIF SD $\tau$ for both the $\mu$ and $\sigma$ parameters, which was especially evident at $\tau =0.6$, where the coverage fell below 92.5%. This contrasts with the JK and AJK methods, which maintained closer adherence to the 95% target rate. As the number of items increased from $I=10$ items to $I=80$ items, the differences in the coverage rates between the three methods became smaller. These findings suggest that a larger number of items improved the quality of linking error estimation based on the Taylor method.

Figure 1 visually represents the results from Table 1 for the AJK and TAY methods. The JK method is not displayed in this figure, as its results closely mirrored those of the AJK method. Overall, Figure 1 underscores that the AJK method provides more reliable coverage across various DIF SD and item number conditions. At the same time, TAY is less precise, with a substantial undercoverage, especially in the scenarios with a large DIF SD and small number of items.

4. Simulation Study 2: Finite Sample Sizes

4.1. Method

Simulation Study 2 extends Simulation Study 1 to finite sample sizes. In this case, the variability of the estimated group mean $\widehat{\mu }$ and estimated SD $\widehat{\sigma }$ arises from the sampling of persons (i.e., standard error) and random DIF (i.e., linking error), resulting in a total error that quantifies the overall uncertainty.

Item responses were generated based on the 2PL model for the two groups, with the tests containing $I=10$, 20, or 40 items. The item parameters matched those used in Simulation Study 1 (see Section 3). The factor variable $\theta$ followed a normal distribution in both groups. As in Simulation Study 1, the $\theta$ distribution of the first group had a mean of 0 and an SD of 1. For the second group, the random variable $\theta$ was fixed at a mean $\mu =0.3$ and an SD $\sigma =1.2$. We modeled random DIF effects using a normal distribution, setting the DIF SD $\tau$ of item difficulties to 0, 0.3, and 0.6. We employed the same DGM for item parameters as in Simulation Study 1 (see Section 3, Equations (45) and (46)). The DIF SD in log-transformed item discriminations was chosen as $0.3·\tau$, resulting in SDs of 0, 0.1, and 0.2. Note that the condition $\tau =0$ corresponds to a scenario with no DIF. Sample sizes of $N=500$, 1000, 2000, and 4000 were selected to reflect typical sample sizes used in 2PL model medium- to large-scale testing applications.

Unlike in Simulation Study 1, the item parameters for the 2PL model were estimated separately for each group using MML estimation in the first step. In the second step, SL linking was performed. Because of the findings in Simulation Study 1, we only evaluated the approximate jackknife (AJK) method for linking error estimation, yielding the estimation method “LE” for the linking error (see Section 2.4). Moreover, we probed the bias-correction method for estimating the linking error based on the approximate jackknife technique, denoted by “${\mathrm{LE}}_{\mathrm{bc}}$” (see Section 2.5). Standard errors (SE) were computed based on the delta method as described in Section 2.1. Based on these quantities, the overall uncertainty was quantified by the total error $\mathrm{TE}$ and the bias-corrected total error ${\mathrm{TE}}_{\mathrm{bc}}$ defined by

$$\mathrm{TE}=\sqrt{{\mathrm{SE}}^2+{\mathrm{LE}}^2}\mathrm{and}{\mathrm{TE}}_{\mathrm{bc}}=\sqrt{{\mathrm{SE}}^2+{\mathrm{LE}}_{\mathrm{bc}}^2}.$$

(48)

In all, 4000 replications were conducted in each of the 3 (DIF SD $\tau$) × 3 (number of items I) × 4 (sample size N) = 36 cells of the simulation. In our analysis of this simulation study, we calculated the median of the linking error estimates $\mathrm{LE}$ and ${\mathrm{LE}}_{\mathrm{bc}}$ for $\widehat{\mu }$ and $\widehat{\sigma }$. We opted for the median instead of the mean in order to carefully investigate the performance of the estimate in the scenario with no DIF. As in Simulation Study 1, we mainly focused on assessment of the adequacy coverage rates for the estimates $\widehat{\mu }$ and $\widehat{\sigma }$, which were based on the standard error $\mathrm{SE}$, the uncorrected total error $\mathrm{TE}$, and the bias-corrected total error ${\mathrm{TE}}_{\mathrm{bc}}$.

All analyses in this simulation study were conducted using R software (Version 4.4.1; [28]). The 2PL model was estimated with the sirt::xxirt() function from the R package sirt (Version 4.2-73; [29]). Dedicated R functions were written for SL linking estimation and the computation of the standard, linking, and total errors. Replication materials for Simulation Study 2 are available at https://osf.io/4dr2c (accessed on 28 October 2024).

4.2. Results

Table 2. Simulation Study 2: Median of linking error estimates LE and ${\mathrm{LE}}_{\mathrm{bc}}$ for the mean $\mu$ and the standard deviation $\sigma$ as a function of the DIF SD $\tau$, the number of items I and sample size N.

		$\mathit{\mu }$						$\mathit{\sigma }$
		$\mathit{\tau }=\mathbf{0}$		$\mathit{\tau }=\mathbf{0.3}$		$\mathit{\tau }=\mathbf{0.6}$		$\mathit{\tau }=\mathbf{0}$		$\mathit{\tau }=\mathbf{0.3}$		$\mathit{\tau }=\mathbf{0.6}$
$\mathit{I}$	$\mathit{N}$	LE	${\mathbf{LE}}_{\mathbf{bc}}$	LE	${\mathbf{LE}}_{\mathbf{bc}}$	LE	${\mathbf{LE}}_{\mathbf{bc}}$	LE	${\mathbf{LE}}_{\mathbf{bc}}$	LE	${\mathbf{LE}}_{\mathbf{bc}}$	LE	${\mathbf{LE}}_{\mathbf{bc}}$
10	500	0.053	0.000	0.117	0.102	0.215	0.207	0.085	0.000	0.101	0.044	0.133	0.095
	1000	0.038	0.000	0.109	0.102	0.211	0.207	0.060	0.000	0.079	0.048	0.117	0.098
	2000	0.027	0.000	0.110	0.107	0.213	0.211	0.043	0.000	0.068	0.051	0.111	0.101
	4000	0.019	0.000	0.106	0.104	0.208	0.207	0.030	0.000	0.059	0.050	0.104	0.099
20	500	0.035	0.000	0.080	0.071	0.146	0.141	0.051	0.000	0.062	0.032	0.086	0.067
	1000	0.025	0.000	0.076	0.071	0.145	0.143	0.036	0.000	0.050	0.034	0.077	0.067
	2000	0.018	0.000	0.073	0.071	0.143	0.141	0.025	0.000	0.043	0.034	0.071	0.066
	4000	0.013	0.000	0.072	0.071	0.143	0.142	0.018	0.000	0.039	0.034	0.070	0.067
40	500	0.024	0.000	0.055	0.049	0.102	0.099	0.033	0.000	0.041	0.023	0.057	0.046
	1000	0.017	0.000	0.053	0.050	0.100	0.099	0.023	0.000	0.033	0.024	0.052	0.047
	2000	0.012	0.000	0.051	0.050	0.099	0.099	0.016	0.000	0.029	0.023	0.049	0.046
	4000	0.009	0.000	0.050	0.049	0.099	0.099	0.012	0.000	0.026	0.024	0.048	0.046

Note. LE = linking error; ${\mathrm{LE}}_{\mathrm{bc}}$ = bias-corrected linking error.

shows median estimates of the linking error LE and bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$ for the group mean $\mu$ and group SD $\sigma$ in finite sample settings. The results are analyzed across various conditions of DIF SD $\tau$, number of items I, and sample sizes N. It turned out that both median LE and ${\mathrm{LE}}_{\mathrm{bc}}$ values increased with a higher amount of DIF (i.e., with increasing DIF SD $\tau$), which was particularly noticeable with smaller sample sizes and fewer items. For example, the LE for $\mu$ at $\tau =0.6$ was markedly higher than for $\tau =0$, suggesting that the DIF SD significantly impacted the linking error, as is known from the literature. Across all conditions, the bias-corrected linking error ${\mathrm{LE}}_{\mathrm{bc}}$ was generally lower than the corresponding LE values, illustrating that the bias correction technique was effective in reducing the bias during estimation of the linking error LE under the no-DIF condition ($\tau =0$). The difference between LE and ${\mathrm{LE}}_{\mathrm{bc}}$ was more pronounced with smaller numbers of items and smaller sample sizes. As the sample size or the number of items increased, the median LE and ${\mathrm{LE}}_{\mathrm{bc}}$ estimates decreased, demonstrating that larger samples and more items result in smaller linking errors due to random DIF.

Table 3. Simulation Study 2: Coverage rates at the 95% confidence level for the mean $\mu$ and the standard deviation $\sigma$ as a function of the DIF SD $\tau$, the number of items I and sample size N.

			$\mathit{\tau }=\mathbf{0}$			$\mathit{\tau }=\mathbf{0.3}$			$\mathit{\tau }=\mathbf{0.6}$
Par	$\mathit{I}$	$\mathit{N}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$
$\mu$	10	500	95.4	98.1	96.1	80.8	97.3	95.5	57.1	94.8	93.8
		1000	95.0	97.8	95.5	70.4	96.3	95.0	44.0	94.8	94.3
		2000	91.6	97.7	95.6	55.5	94.9	94.3	32.2	94.5	94.1
		4000	94.9	98.1	95.6	44.0	94.6	94.0	23.8	94.1	94.0
	20	500	95.4	97.0	95.6	85.1	96.1	95.1	67.9	96.0	95.2
		1000	95.3	97.1	95.6	78.1	95.8	95.1	53.7	95.2	94.8
		2000	95.7	97.3	96.0	66.3	95.8	95.2	41.3	94.8	94.6
		4000	95.2	97.0	95.4	51.3	94.8	94.4	29.9	94.8	94.7
	40	500	93.6	94.8	93.6	88.4	94.7	93.8	73.7	94.1	93.8
		1000	94.7	95.6	94.9	84.2	95.2	94.6	64.6	95.3	95.0
		2000	94.0	95.3	94.2	76.7	95.2	94.9	50.4	95.0	94.8
		4000	92.9	94.6	93.1	61.8	94.7	94.4	38.4	95.2	95.1
$\sigma$	10	500	94.9	99.0	95.9	92.0	98.8	95.1	81.8	97.7	92.8
		1000	94.9	99.1	96.2	88.1	98.0	93.9	73.7	97.0	94.1
		2000	92.1	99.1	96.2	77.8	97.3	93.4	57.8	95.4	93.5
		4000	95.1	99.3	96.3	71.4	97.2	93.9	48.6	96.3	95.0
	20	500	95.0	98.0	95.5	93.3	98.0	95.4	86.5	97.7	95.2
		1000	95.1	98.5	95.9	89.9	97.9	95.5	77.7	97.1	95.3
		2000	95.6	98.4	96.3	84.8	97.5	94.9	65.5	96.2	95.0
		4000	94.9	98.2	95.6	77.3	96.8	95.4	53.2	95.4	94.8
	40	500	94.8	97.3	95.0	93.4	97.3	95.0	89.0	96.6	94.9
		1000	95.0	97.1	95.2	91.9	96.9	95.2	83.4	96.6	95.3
		2000	94.4	97.0	94.7	88.1	96.2	94.6	75.0	96.1	95.5
		4000	94.8	97.2	95.1	82.5	96.2	94.9	60.6	95.6	95.2

Note. Par = parameter; SE = standard error; TE = total error; ${\mathrm{TE}}_{\mathrm{bc}}$ = bias-corrected total error; Coverage rates smaller than 92.5 and larger than 97.5 are printed in bold font.

presents coverage rates at the 95% confidence level for the mean $\mu$ and SD $\sigma$ across different DIF SDs $\tau$, item numbers I, and sample sizes N. We compared three error estimation approaches in all conditions: the standard error (SE), total error (TE), and bias-corrected total error (${\mathrm{TE}}_{\mathrm{bc}}$). As the DIF SD $\tau$ increased from 0 to 0.6, coverage rates for SE consistently decreased, particularly for large sample sizes. For instance, when $I=10$ and $N=500$, the SE coverage rate for $\mu$ dropped from 95.4% at $\tau =0$ to 57.1% at $\tau =0.6$, indicating that coverage rates based on SE are inadequate under conditions with large random DIF. In contrast, coverage rates based on TE and ${\mathrm{TE}}_{\mathrm{bc}}$ were more adequate, with the bias-corrected linking error ${\mathrm{TE}}_{\mathrm{bc}}$ typically showing the best performance. For example, with $\tau =0.6$, $I=10$, and $N=500$, the coverage rate of ${\mathrm{TE}}_{\mathrm{bc}}$ for $\sigma$ remained above 93%, while SE coverage was significantly lower at about 82% and TE had an inflated coverage rate of about 98%. This suggests that ${\mathrm{TE}}_{\mathrm{bc}}$ is more resilient to a larger extent of random DIF. Overall, Table 3 highlights that while SE and TE coverage rates respectively resulted in substantial undercoverage or overcoverage rates, ${\mathrm{TE}}_{\mathrm{bc}}$ maintained adequate coverage rates across various simulation conditions.

Finally, Figure 2 summarizes the coverage rates for the SE, TE, and ${\mathrm{TE}}_{\mathrm{bc}}$ error estimation methods across all simulation conditions in the form of a boxplot. The boxplot shows how each error estimation method performed in maintaining accurate 95% coverage rates under varying simulation conditions. The ${\mathrm{TE}}_{\mathrm{bc}}$ method generally showed the most consistent coverage close to the desired 95% coverage rate across all conditions. This finding indicates the effectiveness of the ${\mathrm{TE}}_{\mathrm{bc}}$ estimation method for adequately assessing the impact of random DIF in linking parameter estimates $\widehat{\mu }$ and $\widehat{\sigma }$. Both the SE and TE estimation methods displayed more variability in their coverage rates, with SE often undershooting the target, particularly under higher DIF conditions. This suggests that SE and TE are less reliable in accurately estimating the total amount of uncertainty in linking parameter estimates.

5. Discussion

This article has examined error estimation for Stocking–Lord linking in dichotomous responses, focusing on sampling errors from person sampling and linking errors from random differential item functioning. By introducing alternatives to the computationally demanding jackknife method, we offer two methods—the approximate jackknife approach and a Taylor approximation-based approach—that present efficient substitutes without substantially compromising accuracy. Our results reveal that the approximate jackknife approach performs almost identically to the traditional method, providing a practical tool for scenarios requiring reduced computational overhead.

Moreover, we have developed a bias-corrected linking error estimate that improves reliability, particularly with smaller sample sizes or fewer items. This approach significantly mitigates the influence of sampling error, yielding more robust linking error estimates. Coverage rates based on the bias-corrected total error surpassed those obtained with the conventional total error, highlighting its value for enhancing accuracy in educational testing and other psychometric applications.

However, this study has certain limitations. The error estimation methods were only examined for dichotomous items, suggesting a need for future extensions to polytomous items. Moreover, the linking error estimation technique for Stocking–Lord linking was applied to the 2PL, but can also be used for group comparisons involving the one-parameter, three-parameter, or four-parameter logistic IRT models. Additionally, the studied methods have yet to account for dependent item responses, such as those clustered in testlets. Exploring linking error estimation in multidimensional IRT models presents another promising area for future work. Addressing these limitations could expand the applicability and effectiveness of linking and error estimation methods, especially for complex assessment designs.

As noted by an anonymous reviewer, the group comparison focused solely on two parameters: the mean difference $\mu$, and the relative change in standard deviation $\sigma$. However, other distributional changes, such as quantiles or the proportion of individuals exceeding a specific cutoff value, may also be of interest. When the $\theta$ variable follows a normal distribution, these additional parameters are nonlinear functions of $\mu$ and $\sigma$, and their associated error estimates can be derived using the delta method.

Driven by the premise that items exhibiting DIF should not be excluded from group comparisons, this article has examined linking errors under random DIF. Excluding such items risks altering the meaning of the construct, as DIF items are considered construct-relevant. In contrast, much of the existing literature assumes a subset of items with fixed DIF, with most items remaining unaffected. Under this assumption, failing to exclude DIF items could lead to biased group comparisons. However, this approach treats DIF items as construct-irrelevant, an assumption that we argue is rarely justifiable in empirical contexts.

This paper has focused on comparing two groups, without specifying the precise meaning of the grouping variable. This variable could represent demographic categories such as gender, country, or socioeconomic background; alternatively, it could define groups based on experimental treatments. The proposed methodology is also applicable to comparisons across different grades or time points.

In summary, our findings suggest that the approximate jackknife and bias-corrected linking error estimates offer valuable, efficient alternatives for Stocking–Lord linking. These contributions could enhance both the precision and practicality of linking error estimation, with implications for various psychometric and educational measurement contexts.

6. Conclusions

This article presents an efficient approach for estimating bias-corrected linking error and total error in Stocking–Lord linking. Two simulation studies demonstrate the necessity of bias correction to avoid inflated coverage rates in linking parameter estimates. The jackknife linking error estimate for Stocking–Lord linking is implemented in a computationally efficient manner, eliminating the need for repeated minimization of the linking function.