Linking Error Estimation in Haberman Linking

Abstract

Haberman linking is a widely used method for comparing groups using the two-parameter logistic item response model. However, the traditional Haberman linking approach relies on joint item parameter estimation, which prevents the application of standard M-estimation theory for linking error calculation in the presence of differential item functioning. To address this limitation, a novel pairwise Haberman linking method is introduced. Pairwise Haberman linking aligns with Haberman linking when no items are missing but eliminates the need for joint item parameters, allowing for the use of M-estimation theory in linking error computation. Theoretical derivations and simulation studies show that pairwise Haberman linking delivers reliable statistical inferences for items and persons, particularly in terms of coverage rates. Furthermore, using a bias-corrected linking error is recommended to reduce the influence of sample size on error estimates.

1. Introduction

Item response theory (IRT) models [1,2] are statistical models to analyze multivariate discrete random variables. In this article, we focus on dichotomous (i.e., binary) random variables and the comparison of multiple groups by means of linking procedures. Let $\mathit{X}=(X_1,\dots ,X_I)$ denotes a vector of I random variables $X_i\in \{0,1\}$, in the literature referred to as items or item responses. The set of items is also referred to as a (e.g., cognitive or personality) test in the social sciences. A unidimensional IRT model [3] is a statistical model for the probability distribution $\mathsf{P}(\mathit{X}=\mathit{x})$ for $\mathit{x}=(x_1,\dots ,x_I)\in {\{0,1\}}^I$

$${\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)

where $\varphi$ is the density of the normal distribution, parameterized by the mean $\mu$ and the standard deviation (SD) $\sigma$. The distribution parameters of the latent ability variable are contained in $\mathit{\delta }=(\mu ,\sigma )$. Note that the normal distribution of the $\theta$ variable is not required [4,5], but is often convenient in estimating the IRT model. Since the normal distribution serves as a prior distribution for the individual posterior distributions $\mathsf{P}\left(\theta \right|\mathit{X}=\mathit{x})$, its influence diminishes as the number of items increases. The vector $\mathit{\gamma }=({\mathit{\gamma }}_1,\dots ,{\mathit{\gamma }}_I)$ includes the item parameters associated with 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 [6] is expressed as

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

(2)

where $a_i$ and $b_i$ denote the item discrimination and the item difficulty $b_i$, respectively, and $\Psi \left(x\right)={(1+exp(-x))}^{-1}$ is the logistic distribution function. Here, the item parameter vector is given by ${\mathit{\gamma }}_i=(a_i,b_i)$.

It is important to note that the item parameters $(a_i,b_i$) for $i=1,\dots ,I$ and distribution parameters $(\mu ,\sigma )$ cannot be disentangled from each other without additional identification constraints [7]. If the parameters $(a_i,b_i)$ and $(\mu ,\sigma )$ represent the multivariate distribution $\mathsf{P}(\mathit{X}=\mathit{x})$ in the IRT model (1), the model could alternatively be represented by item parameters $({\tilde{a}}_i,{\tilde{b}}_i)$ and $\mu =0$ and $\sigma =1$, where

$${\tilde{a}}_i=a_i\sigma \mathrm{and}{\tilde{b}}_i={\sigma }^{-1}(b_i-\mu ).$$

(3)

For a sample of N individuals with independently and identically distributed observations ${\mathit{x}}_1,\dots ,{\mathit{x}}_N$ drawn from the distribution of the random variable $\mathit{X}$, the unknown parameters of the IRT model in (1) can be consistently estimated using marginal maximum likelihood estimation (MML; [8,9]).

IRT models are frequently used to assess and compare the performance of two groups on a test by analyzing differences in the latent variable $\theta$, as defined in the IRT model in Equation (1). This paper focuses on linking methods [10] based on the 2PL model. In the initial step of the linking process, the 2PL model is estimated separately for each group, allowing for the possibility of differential item functioning (DIF), where items may exhibit varying behaviors across groups [11,12,13,14]. The absence of DIF is labeled as measurement invariance [15]. In the second step, differences in item parameters are utilized to determine group differences in $\theta$ through a linking method [10,16,17,18].

The presence of DIF introduces additional variability in the estimated group means and group SDs [19,20,21,22,23]. As a result, the estimated distribution parameters are influenced by the selection of items, even when the sample size of persons is infinite. This variability is captured in the linking error [24,25,26,27,28,29].

The total uncertainty consists of both the ordinary standard error and the linking error, collectively referred to as the total error. Previous work has derived and evaluated estimates for linking errors in methods based on two groups [30,31] or fixed item parameter calibration [32,33]. For comparisons across multiple groups, Haberman linking (HL; [34]) is widely used (e.g., [35,36,37,38,39]). An efficient statistical method for calculating the linking error and total error for HL is still lacking in the literature. An exception is the method proposed in [24], which estimates total error using the computationally intensive double jackknife approach for subjects and items. This article addresses that gap by applying M-estimation theory, as outlined in [30]. The M-estimation theory for the linking method has been described in more detail and refined in [31], but the performance of the linking error and total error estimates was only evaluated for linking methods involving two groups. This paper introduces an adapted version of HL, called pairwise HL (PHL), which avoids the need to estimate joint item parameters. When all items are observed in every group, the HL and PHL methods yield equivalent parameter estimates. This article derives linking errors and total errors for the PHL method based on M-estimation theory and evaluates its performance through a simulation study.

The structure of the article is as follows: Section 2 introduces the proposed PHL approach. Section 3 outlines the estimation of the linking error and total error associated with the PHL method. Section 4 reports the results of a simulation study evaluating the accuracy of the linking error and total error estimates for the PHL approach. Finally, the article concludes with a discussion in Section 5.

2. Pairwise Haberman Linking Approach

In this section, we revisit the estimation of Haberman linking (HL) as originally proposed by [34]. This method estimates group means ($\mathit{\mu }$) and standard deviations ($\mathit{\sigma }$), along with joint item discriminations ($\mathit{a}=(a_1,\dots ,a_I)$) and item difficulties ($\mathit{b}=(b_1,\dots ,b_I)$). However, this formulation has a notable limitation: M-estimation theory cannot be directly applied because the number of estimated parameters increases with the number of items. In the context of statistical inference based on subject sampling, this phenomenon is referred to as the incidental parameter problem in statistics [40,41].

To address this issue, we propose a pairwise Haberman linking (PHL) approach. This method estimates only the distribution parameters ($\mathit{\mu }$ and $\mathit{\sigma }$) and defines an optimization function based on differences in item parameters. When all items are available across all groups, PHL is equivalent to HL. However, discrepancies arise when certain items are missing in some groups.

The HL and PHL methods use estimated item parameters ${\widehat{a}}_{ig}$ and ${\widehat{b}}_{ig}$ ($i=1,\dots ,I$; $g=1,\dots ,G$) as inputs, which are obtained from the separate estimation of the 2PL model in each group with a fixed mean $\mu =0$ and a fixed SD $\sigma =1$. Note that these parameters also depend on the unknown group means ${\mu }_g$ and group SDs ${\sigma }_g$ (see (3)).

The traditional HL method typically involves two steps [34,42]. First, the log-transformed group standard deviations, $s_g=log{\sigma }_g$, are estimated ($g=1,\dots ,G$). In the second step, group means ${\mu }_g$ are estimated. For identification purposes, we set $s_1=0$ (implying ${\sigma }_1=1$) and ${\mu }_1=0$.

Let $\mathit{s}=(s_2,\dots ,s_G)$ and $\mathit{\mu }=({\mu }_2,\dots ,{\mu }_G)$. For each item i and group g, define dummy variables $d_{ig}$, which take the value 1 if item i is observed in group g and 0 if the item is missing. In a fully crossed design where all items are observed in all groups, all $d_{ig}$ values are equal to 1.

Log-transformed SDs $\mathit{s}$ and log-transformed joint item discriminations $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_I)$, where $a_i=exp\left({\alpha }_i\right)$ for $i=1,\dots ,I$, are obtained by minimizing the following optimization function:

$$H_1(\mathit{s},\mathit{\alpha })=\sum _{i=1}^I\sum _{g=1}^G{d}_{ig}\rho \left(log{\widehat{a}}_{ig}-{\alpha }_i-s_g\right),$$

(4)

where $\rho$ denotes the square loss function $\rho \left(x\right)=x^2/2$, also referred to as the $L_2$ loss function. If estimates $\widehat{\mathit{s}}=({\widehat{s}}_2,\dots ,{\widehat{s}}_G)$ were obtained by minimizing $H_1$ in (4), the group SD estimates are calculated as ${\widehat{\sigma }}_g=exp\left(s_g\right)$. Subsequently, group means $\mathit{\mu }$ and joint item difficulties $\mathit{b}$ are determined by minimizing

$$H_2(\mathit{\mu },\mathit{b})=\sum _{i=1}^I\sum _{g=1}^G{d}_{ig}\rho \left(exp\left({\widehat{s}}_g\right){\widehat{b}}_{ig}-b_i+{\mu }_g\right),$$

(5)

where, as before, $\rho$ denotes the square loss function. The distribution parameter estimates $\widehat{\mathit{s}}$ and $\widehat{\mathit{\mu }}$, along with joint item parameter estimates $\widehat{\mathit{\alpha }}$ and $\widehat{\mathit{b}}$, fulfill the following estimating equations, derived by taking the required partial derivatives in (4) and (5):

$$\sum _{g=1}^G{d}_{ig}\left(log{\widehat{a}}_{ig}-{\alpha }_i-s_g\right)=0\mathrm{and}\sum _{g=1}^G{d}_{ig}\left(exp\left(s_g\right){\widehat{b}}_{ig}-b_i+{\mu }_g\right)=0\mathrm{for}i=1,\dots ,I\mathrm{and}$$

(6)

$$\sum _{i=1}^I{d}_{ig}\left(log{\widehat{a}}_{ig}-{\alpha }_i-s_g\right)=0\mathrm{and}\sum _{i=1}^I{d}_{ig}\left(exp\left(s_g\right){\widehat{b}}_{ig}-b_i+{\mu }_g\right)=0\mathrm{for}g=2,\dots ,G.$$

(7)

Haberman discussed the conditions under which a unique solution is obtained in the minimization of $H_1$ and $H_2$ [34].

For motivating the optimization function for the proposed PHL method, a complete crossed design is assumed, where all items are available in all groups, i.e., $d_{ig}=1$ for all $i=1,\dots ,I$ and $g=1,\dots ,G$. From (7), the following identity holds for all $g\ne h$:

$$\sum _{i=1}^I\left(log{\widehat{a}}_{ig}-{\alpha }_i-s_g\right)-\sum _{i=1}^I\left(log{\widehat{a}}_{ih}-{\alpha }_i-s_h\right)=\sum _{i=1}^I\left(log{\widehat{a}}_{ig}-log{\widehat{a}}_{ih}-s_g+s_h\right)=0.$$

(8)

Consequently, it holds for any group g

$$\sum _{h\ne g}\sum _{i=1}^I\left(log{\widehat{a}}_{ig}-log{\widehat{a}}_{ih}-s_g+s_h\right)=0\mathrm{for}g=2,\dots ,G.$$

(9)

Hence, the log-transformed standard deviations $\mathit{s}$ can also be obtained by minimizing the pairwise differences

$$H_1\left(\mathit{s}\right)=\sum _{i=1}^I{\omega }_i\sum _{g\ne h}{d}_{ig}{d}_{ih}\rho \left(log{\widehat{a}}_{ig}-log{\widehat{a}}_{ig}-s_g+s_h\right).$$

(10)

Note that item-specific weights ${\omega }_i$ are introduced in (10), where they equal 1 in complete crossed designs. In designs with missing items, however, the weights will be chosen differently, as will be clarified later.

Similarly, the identity for group means $\mathit{\mu }$ can be derived as

$$\sum _{i=1}^I\left(exp\left(s_g\right){\widehat{b}}_{ig}-b_i-{\mu }_g\right)-\sum _{i=1}^I\left(exp\left(s_h\right){\widehat{b}}_{ih}-b_i-{\mu }_h\right)=\sum _{i=1}^I\left(exp\left(s_g\right){\widehat{b}}_{ig}-exp\left(s_h\right){\widehat{b}}_{ih}+{\mu }_g-{\mu }_h\right)=0.$$

(11)

Therefore, we arrive at

$$\sum _{h\ne g}\sum _{i=1}^I\left(exp\left(s_g\right){\widehat{b}}_{ig}-exp\left(s_h\right){\widehat{b}}_{ih}+{\mu }_g-{\mu }_h\right)=0\mathrm{for}g=2,\dots ,G.$$

(12)

The vector of group means $\mathit{\mu }$ can be estimated in PHL by minimizing

$$H_2\left(\mathit{\mu }\right)=\sum _{i=1}^I{\omega }_i\sum _{g\ne h}{d}_{ig}{d}_{ih}\rho \left(exp\left({\widehat{s}}_g\right){\widehat{b}}_{ig}-exp\left({\widehat{s}}_h\right){\widehat{b}}_{ih}+{\mu }_g-{\mu }_h\right).$$

(13)

The linking method based on the optimization functions $H_1$ in (10) and $H_2$ in (13) forms the PHL method, which closely resembles the optimization functions used in invariance alignment (IA; [43,44]). Note that PHL will differ from HL in designs with missing items. In the remainder of the article, we consider two versions of PHL: either ${\omega }_i=1$ (method PHL1) or ${\omega }_i=I/G_i$ (method PHL2), where $G_i={\sum }_{g=1}^G{d}_{ig}$. The latter method is motivated by the fact that if item i is observed in $G_i$ groups, there are $G_i(G_i-1)/2$ pairwise differences in the PHL optimization function, whereas HL involves only $G_i$ terms for item i. To make PHL more similar to HL, the weights in PHL2 are chosen accordingly.

The original HL method is implemented in the function equateIRT::multiec() in the R (Version 4.4.1) [45] package equateIRT [46,47], as well as in the functions sirt::linking.haberman() and sirt::linking.haberman.lq() (with the argument method=’joint’) in the R package sirt [45]. Both PHL methods are available in the function sirt::linking.haberman.lq() (with method=’pw1’ or method=’pw2’ for PHL1 and PHL2, respectively) in the same package. SAS code for HL is also available [48].

3. Linking Error Estimation in the Pairwise Haberman Linking Approach

This section presents the derivation of standard errors, linking errors, and total errors for the newly introduced PHL method in Section 2. The theory for deriving the estimates was proposed in earlier work [31] and carved out in more detail for the two-group case in [30].

This article focuses on estimating the uncertainty associated with the group means $\widehat{\mathit{\mu }}$ and log-transformed group standard deviations $\widehat{\mathit{s}}$ within the vector $\mathit{\delta }=(\mathit{\mu },\mathit{s})$. The uncertainty in the linking parameter estimate $\widehat{\mathit{\delta }}$ arises from two sources: the sampling (or selection) of individuals, which corresponds to the scheme $N\to \infty$, and the selection (or sampling) of items, as well as the modeling of variability in group comparisons, which corresponds to the scheme $I\to \infty$.

Much of the linking literature addresses the issue of uncertainty in item parameter estimates $\widehat{\mathit{\gamma }}$ by calculating the standard error (SE) of $\widehat{\mathit{\delta }}$ for a fixed number of items [46,49,50,51,52]. In this context, variability in population-level item parameters $\widehat{\mathit{\gamma }}$ (i.e., for $N=\infty$) arises due to sampling variability in the estimated item parameters ${\widehat{\mathit{\gamma }}}_i$, which results from the sampling of individuals. If DIF occurs, the estimated linking parameter $\widehat{\mathit{\delta }}$ will depend on the selected set of items, even with an infinite sample size. This variability, known as linking error (LE) [22,25,27], is handled under the scheme $I\to \infty$ for variance estimation.

The total error (TE) encompasses both sources of uncertainty: the standard error arising from randomness in the sampling of persons and the linking error resulting from variability in item selection [22,24,27,28,53]. However, it has been suggested that the traditional LE estimate may be influenced by sampling error [30,31]. To address this, a bias-corrected version of the LE is additionally proposed as an alternative to the traditional LE estimate.

We now present the theoretical framework for a general linking function based on M-estimation [54,55,56,57]. The vector of input parameters for the linking method is given by ${\widehat{\mathit{\gamma }}}_i=({\widehat{a}}_{i1},{\widehat{b}}_{i1},d_{i1},\dots ,{\widehat{a}}_{iG},{\widehat{b}}_{iG},d_{iG})$. As the sample size N increases (assuming that the sample sizes per group also increase) ${\widehat{\mathit{\gamma }}}_i$ stochastically converges to the true parameter vector ${\mathit{\gamma }}_i=(a_{i1},b_{i1},d_{i1},\dots ,a_{iG},b_{iG},d_{iG})$. Let $\widehat{\mathit{\gamma }}=({\widehat{\mathit{\gamma }}}_1,\dots ,{\widehat{\mathit{\gamma }}}_I)$ be the vector of all item parameters, with an estimated variance matrix ${\mathit{V}}_{\widehat{\mathit{\gamma }}}$, typically obtained from software that estimates IRT models. The aim of a linking method is to estimate the vector of distribution parameters $\mathit{\delta }$ using the set of all estimated item parameters ${\widehat{\gamma }}_i$ for $i=1,\dots ,I$.

The linking parameter estimate $\widehat{\mathit{\delta }}$ is obtained by solving a multivariate estimating equation with respect to $\mathit{\delta }$

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

(14)

where $H_{\mathit{\delta }}$ and $h_{\mathit{\delta }}$ denote the partial derivatives of univariate or multivariate functions H and h with respect to the parameter $\mathit{\delta }$. In HL and PHL, the vector $H_{\mathit{\delta }}$ consists of the derivatives of the optimization functions $H_1$ and $H_2$.

In the PHL method, the item-wise contributions $h_{\mathit{\delta }}(\mathit{\delta };{\widehat{\mathit{\gamma }}}_i)$ in (14) are given by (see Section 2)

$$h_{\mathit{\delta }}(\mathit{\delta };{\widehat{\mathit{\gamma }}}_i)={\omega }_i\left(\begin{array}{c}{\displaystyle \sum _{h\ne 2}{d}_{i2}{d}_{ih}\left(exp\left(s_2\right){\widehat{b}}_{i2}-exp\left(s_h\right){\widehat{b}}_{ih}+{\mu }_2-{\mu }_h\right)}\\ ⋮\\ {\displaystyle \sum _{h\ne G}{d}_{iG}{d}_{ih}\left(exp\left(s_G\right){\widehat{b}}_{iG}-exp\left(s_h\right){\widehat{b}}_{ih}+{\mu }_G-{\mu }_h\right)}\\ {\displaystyle \sum _{h\ne 2}{d}_{i2}{d}_{ih}\left(log{\widehat{a}}_{i2}-log{\widehat{a}}_{ih}-s_2+s_h\right)}\\ ⋮\\ {\displaystyle \sum _{h\ne G}{d}_{iG}{d}_{ih}\left(log{\widehat{a}}_{iG}-log{\widehat{a}}_{ih}-s_G+s_h\right)}\end{array}\right).$$

(15)

Note that the dummy variables $d_{ig}$ in (15) indicate which terms in the estimating equations are used when an item is observed in a particular combination of groups g and h.

3.1. Standard Error

We now compute the standard error of $\widehat{\mathit{\delta }}$ due to sampling of persons (see [51,58,59]). A Taylor approximation of $h_{\mathit{\delta }}$ around $(\mathit{\delta },{\mathit{\gamma }}_i)$ is performed, yielding

$$h_{\mathit{\delta }}(\widehat{\mathit{\delta }};{\widehat{\mathit{\gamma }}}_i)\begin{array}{c}\hfill \simeq \end{array}{h}_{\mathit{\delta }}(\mathit{\delta },{\mathit{\gamma }}_i)+h_{\delta \gamma }(\mathit{\delta },{\mathit{\gamma }}_i)({\widehat{\mathit{\gamma }}}_i-{\mathit{\gamma }}_i)+h_{\delta \delta }(\mathit{\delta },{\mathit{\gamma }}_i)(\widehat{\mathit{\delta }}-\mathit{\delta }).$$

(16)

We can now use

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

(17)

$$\widehat{\mathit{\delta }}-\mathit{\delta }\begin{array}{c}\hfill \simeq \end{array}{\left(\sum _{i=1}^I{h}_{\delta \delta }(\mathit{\delta },{\mathit{\gamma }}_i)\right)}^{-1}\sum _{i=1}^I{h}_{\delta \gamma }(\mathit{\delta },{\mathit{\gamma }}_i)({\widehat{\mathit{\gamma }}}_i-{\mathit{\gamma }}_i)={\mathit{A}}^{-1}\mathit{C}(\widehat{\mathit{\gamma }}-\mathit{\gamma }),\mathrm{where}$$

(18)

$$\mathit{A}=\sum _{i=1}^I{h}_{\delta \delta }(\mathit{\delta },{\mathit{\gamma }}_i)\mathrm{and}\mathit{C}=\left(\begin{array}{ccc}{h}_{\delta \gamma }(\mathit{\delta };{\mathit{\gamma }}_1)& \cdots & h_{\delta \gamma }(\mathit{\delta };{\mathit{\gamma }}_I)\end{array}\right).$$

(19)

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

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

(20)

The unknown quantities in (20) can be estimated using

$$\widehat{\mathit{C}}=\left(\begin{array}{ccc}{h}_{\delta \gamma }(\widehat{\mathit{\delta }},{\widehat{\mathit{\gamma }}}_1)& \dots & h_{\delta \gamma }(\widehat{\mathit{\delta }},{\widehat{\mathit{\gamma }}}_I)\end{array}\right)\mathrm{and}\widehat{\mathit{A}}=\sum _{i=1}^I{h}_{\delta \delta }(\widehat{\mathit{\delta }},{\widehat{\mathit{\gamma }}}_i).$$

(21)

Standard errors for entries in $\widehat{\mathit{\delta }}$ are obtained by taking the square root of the entries in the diagonal of the estimate of ${\mathit{V}}_{\mathrm{SE}}$.

The entries in $h_{\delta \delta }$ and $h_{\delta \gamma }$ for the PHL method can be computed in a straightforward manner. For example, the entries in $h_{\delta \delta }$ are given by

$$h_{s_g{s}_g}=-\sum _{i=1}^I{\omega }_i{d}_{ig}(G_i-1),h_{s_g{s}_h}=\sum _{i=1}^I{\omega }_i{d}_{ig}{d}_{ih}\mathrm{for}g\ne h,$$

(22)

$$h_{{\mu }_g{\mu }_g}=\sum _{i=1}^I{\omega }_i{d}_{ig}(G_i-1),h_{{\mu }_g{\mu }_h}=-\sum _{i=1}^I{\omega }_i{d}_{ig}{d}_{ih}\mathrm{for}g\ne h,$$

(23)

$$h_{{\mu }_g{s}_g}=exp\left(s_g\right)\sum _{i=1}^I{\omega }_i{d}_{ig}(G_i-1){\widehat{b}}_{ig},h_{{\mu }_g{s}_h}=-exp\left(s_h\right)\sum _{i=1}^I{\omega }_i{d}_{ig}{d}_{ih}{\widehat{b}}_{ih}\mathrm{for}g\ne h.$$

(24)

Moreover, we have $h_{s_g{\mu }_h}=0$ for all $g,h=2,\dots ,G$.

3.2. Linking Error

The linear Taylor approximation (18) can be used to derive the variance to item selection. The variance matrix based on M-estimation is given by

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

(25)

where $\mathit{A}$ is given in (19) and $\mathit{B}$ is obtained by

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

(26)

In (26), the independence assumption of items is applied. In M-estimation, the matrix $\mathit{A}$ is referred to as the bread matrix, and $\mathit{B}$ as the meat matrix. The unknown quantities in (25) can be estimated by $\widehat{\mathit{A}}$ (see (21)) and

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

(27)

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

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

(28)

The factor $I/(I-1)$ in (28) is included to improve the statistical properties of the linking error estimate for a small number of items [30,60,61,62].

The corresponding linking errors are computed by taking the square roots of the diagonal elements in ${\widehat{\mathit{V}}}_{\mathrm{LE}}$. In a practical implementation of the linking error in PHL, the partial derivatives ${\mathit{h}}_{\mathit{\delta }{\mathit{\gamma }}_i}$ must be computed. These can also be straightforwardly obtained, similar to the derivations in (22)–(24).

3.3. Bias-Corrected Linking Error

The bias-corrected estimate of the linking error variance matrix ${\mathit{V}}_{\mathrm{LE}}$ is derived as follows. The estimated meat matrix $\widehat{\mathit{B}}$ is expressed as

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

(29)

However, the population-level linking error should ideally be computed using the true item parameters ${\mathit{\gamma }}_i$ rather than the estimated item parameters ${\widehat{\mathit{\gamma }}}_i$ that appear in (29). A linear Taylor approximation provides

$$h_{\mathit{\delta }}(\widehat{\mathit{\delta }};{\widehat{\mathit{\gamma }}}_i)\begin{array}{c}\hfill \simeq \end{array}{h}_{\mathit{\delta }}(\widehat{\mathit{\delta }};{\mathit{\gamma }}_i)+h_{\delta \gamma }(\widehat{\mathit{\delta }},{\mathit{\gamma }}_i)({\widehat{\mathit{\gamma }}}_i-{\mathit{\gamma }}_i).$$

(30)

Thus, the inflated variance contribution in $\widehat{\mathit{B}}$ caused by sampling error is given by

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

(31)

where the approximate independence of item parameters across items was assumed. Based on this, a bias-corrected meat matrix is computed as

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

(32)

Finally, the bias-corrected variance matrix for the linking error is given as

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

(33)

3.4. Total Error

The total uncertainty in $\widehat{\mathit{\delta }}$ (i.e., the total error) can now be quantified. The variance as the total error, defined as the sum of the variances due to sampling error and linking error, is expressed as follows (see [22,27,30]):

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

(34)

The bias-corrected 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}}.$$

(35)

If negative variances are encountered for the bias-corrected linking error or total error estimates, the corresponding error estimate is set to zero.

3.5. Nonlinear Transformation

The error estimates for $\widehat{\mathit{\delta }}$ presented in the previous subsections are computed for group means $\widehat{\mathit{\mu }}$ and log-transformed group SDs $\widehat{\mathit{s}}$. To obtain statistical inference for the SDs $\widehat{\mathit{\sigma }}$, the nonlinear transformation ${\sigma }_g=exp\left(s_g\right)$ must be applied in the delta method. This allows the derivation of an error estimate (i.e., standard error, linking error, or total error) for ${\widehat{\sigma }}_g$ based on the error estimate $\mathrm{Err}\left({\widehat{s}}_g\right)$ for ${\widehat{s}}_g$. The delta formula provides the corresponding estimate for the SD ${\sigma }_g$ as

$$\mathrm{Err}\left({\widehat{\sigma }}_g\right)=exp\left({\widehat{s}}_g\right)\mathrm{Err}\left({\widehat{s}}_g\right).$$

(36)

The error estimates ${\mu }_g$ remain unaffected by the transformation of $\widehat{\mathit{s}}$.

4. Simulation Study

4.1. Method

The 2PL model for $G=4$ groups was used as the IRF in the data-generating model. The $\theta$ variable was assumed to be normally distributed in the four groups with group means ${\mu }_1=0$, ${\mu }_2=0.3$, ${\mu }_3=0.6$, and ${\mu }_4=-0.3$, respectively. The SDs were chosen as ${\sigma }_1=1$, ${\sigma }_2=1.2$, ${\sigma }_3=0.8$, and ${\sigma }_4=1$.

The simulation study was conducted for $I=20$ and $I=40$ items. Group-specific parameters $a_{ig}$ and $b_{ig}$ for each item for groups $g=1,\dots ,4$ were based on fixed base item parameters and DIF effects that were newly simulated in each replication of the simulation. The base item discriminations $a_i$ for the first 10 items were chosen as 0.89, 0.67, 1.44, 1.28, 1.16, 0.84, 1.76, 0.49, 0.80, and 1.19, which resulted in a mean item discrimination of $M=1.052$ and $SD=0.385$. The base item difficulties $b_i$ for the first 10 items were chosen as 0.72, 1.49, −0.88, 0.94, 0.95, −0.66, −1.23, 0.85, 0.53, and 0.29, yielding a mean item difficulty of $M=0.300$ and $SD=0.909$. For the $I=20$ condition, the parameters were doubled, while for $I=40$, they were replicated four times. The item parameters are also available at https://osf.io/bp6em (accessed on 19 November 2024). The group-specific item difficulties $b_{ig}$ ($i=1,\dots ,I;g=1,\cdots ,4$) were simulated as

$$b_{ig}=b_i+e_i\mathrm{for}g=1,\dots ,4,$$

(37)

where $e_i$ is a random uniform DIF effect. Group-specific item discriminations $a_{ig}(g=1,\dots ,4)$ were simulated as

$$a_{ig}=a_{i}exp\left(f_i\right)\mathrm{for}g=1,\dots ,G,$$

(38)

where $f_i$ is a nonuniform random DIF effect. In the simulation, we assumed that the random DIF effects $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$ and followed a normal distribution. The DIF standard deviation for DIF effects $e_i$ in item difficulties was chosen as 0, 0.2, and 0.4, referring to no DIF, small DIF, and large DIF. The corresponding DIF SDs for DIF effects $f_i$ in log-transformed item discriminations were 0, 0.06, and 0.12, respectively.

We also varied the proportion of missing items per group to represent situations where certain items were not administered or deleted due to technical reasons (e.g., translation errors) to some groups. Missing item rates were set at 0%, 10%, and 30% per group with missing items randomly selected in each group for each replication.

Sample sizes per group were chosen as $N=250$, 500, 1000, and 2000 to reflect typical sample sizes in small- to large-scale testing applications of the 2PL model [63,64].

In each of the 4 (sample size N) × 3 (DIF SD $\tau$) × 2 (number of items I) × 3 (proportion of missing items) = 96 cells of the simulation, 3000 replications were conducted. The 2PL model was first estimated separately for each of the four groups. Linking was performed with Haberman linking (HL) and the two pairwise Haberman linking methods, PHL1 and PHL2, based on the item parameter estimates from the 2PL model. For identification reasons, the distribution parameters in the first group were fixed at ${\mu }_1=0$ and ${\sigma }_1=1$.

We evaluated the bias and the root mean square error (RMSE) for the ${\widehat{\mu }}_g$ and ${\widehat{\sigma }}_g$ estimates. Additionally, a  relative RMSE was calculated with HL serving as the reference method (set to a relative RMSE of 100). For the PHL1 and PHL2 methods, standard errors (SE), linking errors (LE), bias-corrected linking errors (${\mathrm{LE}}_{\mathrm{bc}}$), total errors (TE), and bias-corrected total errors (${\mathrm{TE}}_{\mathrm{bc}}$) were computed. Finally, the coverage rates for the ${\widehat{\mu }}_g$ and ${\widehat{\sigma }}_g$ estimates at the 95% confidence level were assessed, defined as the percentage of replications in which the estimated confidence intervals contained the true values of ${\mu }_g$ or ${\sigma }_g$ (for $g=2,3,4$).

All analyses for this simulation study were performed using the statistical software R (Version 4.4.1; [45]). The 2PL model, along with the corresponding standard errors of item parameters, was estimated with the sirt::xxirt() function from the R package sirt (Version 4.2-89; [65]). The HL, PHL1, and PHL2 methods were estimated using the sirt::linking.haberman.lq() function from the same R package, which also provides standard errors, linking errors, and total errors for the PHL1 and PHL2 methods in the output. Replication material for this simulation study is available at https://osf.io/bp6em (accessed on 19 November 2024).

4.2. Results

To assess the similarity between the PHL1 and PHL2 methods with HL, we computed the average absolute differences for all distribution parameters (i.e., ${\widehat{\mu }}_g$ and ${\widehat{\sigma }}_g$ for $g=2,3,4$) in each replication. Figure 1 illustrates the average absolute differences in the estimates of $\widehat{\mu }$ and $\widehat{\sigma }$ between Haberman Linking (HL) and the pairwise Haberman linking methods (PHL1 and PHL2) across selected subsets of simulation conditions. In the 0% missing items condition, the PHL methods coincided with HL by definition of the linking methods. Across all conditions, the absolute differences between PHL1 and HL, and between PHL2 and HL, were very small, indicating that both PHL1 and PHL2 produced results highly comparable to HL. The average absolute differences tended to be slightly smaller as the number of items increased from $I=20$ to $I=40$. Additionally, the absolute differences between the PHL methods and HL increased as the percentage of missing items rose. Overall, the absolute differences for PHL2 relative to HL were generally slightly smaller than those for PHL1, suggesting that PHL2 aligned marginally better with HL, though the differences were minimal.

Figure 2 illustrates the bias in the estimates of the group mean ${\widehat{\mu }}_g$ and SD ${\widehat{\sigma }}_g$ ($g=2,3,4$) for the three linking methods HL, PHL1, and PHL2, pooling all bias values across subsets of simulation conditions. Across all methods, the bias in the ${\widehat{\mu }}_g$ and ${\widehat{\sigma }}_g$ estimates was small. There was no significant difference between HL, PHL1, and PHL2 in terms of bias in $\widehat{\mu }$, indicating that all three linking methods performed comparably in estimating the group means and group SDs. A higher number of items (i.e., $I=40$) slightly decreased the absolute values of the bias.

Figure 3 presents the relative RMSE of ${\widehat{\mu }}_g$ and ${\widehat{\sigma }}_g$ estimates across subsets of simulation conditions. Overall, the RMSE values for the PHL1 method were lower than those for PHL2, though both were generally close to 100, suggesting that the PHL methods did not cause a significant efficiency loss compared with the traditional HL linking method. Notably, a few RMSE values deviated notably from 100, but these occurrences were more frequent when the number of items was small (i.e., $I=20$).

Below, we report only the results for the PHL2 method, as the PHL1 method showed larger mean absolute differences compared to the HL method and led to more variable linking parameter estimates.

Figure 4 illustrates the median values of the linking error estimate (LE) and the bias-corrected linking error (${\mathrm{LE}}_{\mathrm{bc}}$) for the estimated group mean (${\widehat{\mu }}_2$) in the second group, across selected subsets of simulation conditions for the PHL2 linking method. Overall, LE decreased with larger sample sizes, while ${\mathrm{LE}}_{\mathrm{bc}}$ showed the opposite trend. Notably, the bias-corrected linking error exhibited less dependency on sample size than LE. In the absence of DIF (i.e., $\tau =0$), the population-level linking error is zero. The true value of 0 was obtained for ${\mathrm{LE}}_{\mathrm{bc}}$ at the median, while LE was substantially larger and only decreased with larger sample sizes. In the presence of DIF (i.e., $\tau >0$), LE and ${\mathrm{LE}}_{\mathrm{bc}}$ converged as sample size increased.

Table 1. Simulation Study: Coverage rates of of the estimated mean ${\widehat{\mu }}_2$ for the pairwise Haberman linking method PHL2 as a function of the DIF SD $\tau$, number of items I and sample size N for different percentages of missing items.

			0% Missing Items			10% Missing Items			30% Missing Items
$\mathit{\tau }$	$\mathit{I}$	$\mathit{N}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$
0	20	250	95.5	97.8	95.8	95.6	98.1	96.1	95.8	98.4	96.4
		500	94.6	97.4	94.9	95.6	98.0	95.9	95.8	98.7	96.1
		1000	95.2	97.8	95.6	94.6	97.6	94.9	95.0	98.3	95.4
		2000	94.6	97.3	94.9	94.3	97.4	94.8	93.5	98.1	95.5
	40	250	94.5	96.8	94.8	94.6	96.8	94.9	95.2	97.3	95.7
		500	94.3	96.3	94.5	94.9	96.7	95.0	94.9	97.4	95.2
		1000	95.3	96.9	95.3	95.8	97.6	96.0	94.9	97.2	95.1
		2000	94.6	96.6	94.6	94.2	96.5	94.3	95.4	97.5	95.6
0.2	20	250	91.6	97.3	93.9	91.7	97.6	94.2	92.2	98.0	94.6
		500	89.1	96.9	93.8	89.0	96.9	94.0	87.1	97.1	93.6
		1000	82.1	96.0	93.6	82.4	96.9	94.4	78.7	95.9	92.3
		2000	73.6	95.9	94.3	72.6	95.8	94.3	65.1	94.2	91.6
	40	250	92.8	96.5	94.5	92.9	97.3	94.7	92.0	97.2	94.5
		500	90.8	95.9	93.8	90.6	96.5	94.2	89.7	96.8	94.2
		1000	86.7	96.4	94.8	85.4	95.3	93.6	82.6	96.4	94.1
		2000	80.4	96.0	94.8	78.8	95.4	94.0	75.7	95.8	94.4
0.4	20	250	83.5	96.5	94.0	83.2	96.9	94.0	81.9	95.8	92.9
		500	75.0	95.9	94.2	72.9	95.4	93.3	71.9	95.8	93.3
		1000	61.5	94.9	93.4	59.8	95.1	94.0	57.1	94.8	93.2
		2000	48.4	94.5	93.4	45.1	93.8	93.1	41.6	93.9	93.0
	40	250	85.7	95.2	93.8	86.5	96.3	95.0	84.4	96.3	94.3
		500	80.5	95.4	93.8	79.3	95.8	94.7	74.5	95.1	93.6
		1000	70.4	95.0	94.3	68.3	95.7	94.8	64.3	94.3	93.2
		2000	57.7	94.5	93.9	54.1	94.9	94.4	50.5	94.3	93.7

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

provides an overview of the coverage rates for the estimated group mean ${\widehat{\mu }}_2$ for PHL2. Across all conditions, higher percentages of missing items generally led to slightly lower coverage rates. For 0% missing items, ${\mathrm{TE}}_{\mathrm{bc}}$ consistently provided coverage rates within the desired range (close to 95%), outperforming both SE and TE in terms of accuracy. At 30% missing items, ${\mathrm{TE}}_{\mathrm{bc}}$ still maintained coverage near the nominal level of 95%, while TE and SE often fell below this threshold. Overall, statistical inference based on SE resulted in only 33% of the cells meeting the acceptable coverage range, with a systematic undercoverage that produced an average rate of 82.6%. TE performed satisfactorily in 68% of the cells, exhibiting overcoverage with an average coverage rate of 96.4%. In contrast, the ${\mathrm{TE}}_{\mathrm{bc}}$ method excelled in 96% of the cells, with an average coverage rate of 94.4%. Coverage rates based on the bias-corrected total error ${\mathrm{TE}}_{\mathrm{bc}}$ remained robust across all levels of DIF SD $\tau$, consistently yielding better coverage rates compared with TE and SE. When $I=40$, ${\mathrm{TE}}_{\mathrm{bc}}$ achieved better coverage rates than $I=20$ under most conditions, especially for higher levels of missing items and a larger DIF SD $\tau$.

Table 2. Simulation Study: Coverage rates of of the estimated standard deviation ${\widehat{\sigma }}_2$ for the pairwise Haberman linking method PHL2 as a function of the DIF SD $\tau$, number of items I and sample size N for different percentages of missing items.

			0% Missing Items			10% Missing Items			30% Missing Items
$\mathit{\tau }$	$\mathit{I}$	$\mathit{N}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$	SE	TE	${\mathbf{TE}}_{\mathbf{bc}}$
0	20	250	94.8	98.2	95.4	94.8	98.6	95.5	94.6	98.4	95.7
		500	95.1	98.4	95.7	94.9	98.4	95.3	94.9	99.0	95.5
		1000	95.2	98.3	95.6	94.2	98.4	94.8	94.5	98.7	95.0
		2000	94.7	98.1	95.3	95.4	98.9	96.1	94.3	98.9	95.8
	40	250	94.4	97.6	94.9	94.7	97.7	95.1	94.3	97.9	95.1
		500	94.3	97.3	94.5	95.4	97.7	95.6	94.2	98.0	94.6
		1000	95.1	97.3	95.2	94.4	97.4	94.6	95.1	98.4	95.4
		2000	95.0	97.3	95.0	94.8	97.8	95.0	95.4	98.5	95.8
0.2	20	250	94.4	98.4	95.3	93.9	98.2	95.1	94.2	98.3	95.4
		500	94.1	98.1	95.2	93.9	98.7	94.9	93.6	98.8	95.2
		1000	92.8	98.3	94.9	91.9	98.2	94.3	91.4	98.2	93.6
		2000	88.9	97.4	93.8	88.4	97.8	93.6	86.9	97.5	93.2
	40	250	93.5	96.8	94.3	93.7	97.2	94.7	93.5	97.5	94.6
		500	93.3	97.1	94.1	93.1	97.3	94.1	93.3	97.7	94.0
		1000	93.0	96.8	94.7	91.9	97.0	93.5	91.9	97.8	94.0
		2000	91.1	96.9	94.3	90.0	96.7	94.2	90.0	97.8	94.6
0.4	20	250	92.1	97.7	94.7	90.5	97.5	93.9	91.3	98.1	94.7
		500	88.8	97.8	93.7	89.5	97.3	93.5	89.2	98.1	93.5
		1000	83.5	97.5	93.5	83.7	96.9	93.5	83.0	96.7	92.2
		2000	74.6	96.2	93.3	74.3	96.1	93.3	71.9	95.4	91.7
	40	250	90.5	96.6	93.4	91.0	96.2	93.7	91.4	97.3	94.2
		500	92.1	97.4	95.1	91.0	97.1	94.6	89.1	97.0	93.4
		1000	87.0	96.4	94.2	86.5	97.2	94.8	85.0	96.8	93.8
		2000	79.8	96.6	94.3	78.9	96.0	94.2	76.0	96.0	93.6

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

presents the coverage rates for the estimated standard deviation (${\widehat{\sigma }}_2$) using the pairwise Haberman linking method (PHL2). Coverage rates based on TE slightly declined as the percentage of missing items increased. In contrast, coverage rates for ${\mathrm{TE}}_{\mathrm{bc}}$ remained resilient to missing items, staying close to 95% even with 30% missing items. As $\tau$ increased from 0 to 0.4, coverage rates for TE significantly decreased, particularly in scenarios with smaller sample sizes or fewer items, resulting in overcoverage with an average coverage rate of 97.6%. Meanwhile, ${\mathrm{TE}}_{\mathrm{bc}}$ maintained coverage rates near the nominal 95% threshold across all conditions, demonstrating its reliability even under substantial DIF effects, with an average coverage rate of 94.5%. For ${\widehat{\mu }}_2$, coverage rates based on SE were inadequate, showing substantial undercoverage, especially in large samples and with large DIF SD ($\tau$).

5. Discussion

HL is a widely used linking method in educational assessment, particularly in complex designs involving multiple groups. When items exhibit DIF, linking errors can be computed to address the uncertainty in linking parameter estimates resulting from item selection. However, until now, there has been no straightforward solution for computing linking errors for HL. Since HL estimates joint item parameters, standard M-estimation theory for linking error estimation is not applicable. To address this, an alternative pairwise HL approach (denoted PHL) is proposed in this article, which aligns with the traditional HL method when there are no missing items across groups. When some items are missing, PHL provides linking parameter estimates close to those from the original HL approach. A key advantage of the PHL method is that it does not rely on joint item parameters, allowing for the application of M-estimation theory to compute linking errors. This article demonstrates, through theory and a simulation study, that the PHL method can yield adequate statistical inferences for persons and items. Furthermore, it is shown that a bias-corrected linking error should replace the ordinary linking error to prevent the estimate from being influenced by sample size.

As with any simulation study, our research has limitations. First, we focused solely on dichotomous items. Future studies could explore linking error estimation for the PHL method applied to polytomous items [66]. Second, it was assumed that DIF effects were independent across items. In cases where items are grouped within subsets that share a common stimulus (e.g., in reading comprehension tests with a common text stimulus; [67]), the dependence structure of DIF effects should be incorporated into the computation of linking errors. Third, linking error estimates could be extended to more complex data-generating IRT models beyond the 2PL analysis model (e.g., [68]), which may lead to misspecification of the IRT model. Finally, future research could compare the linking errors from PHL using M-estimation theory with those obtained from jackknife methods [24].

As noted by an anonymous reviewer, the assumption of independent DIF effects across items is a strong assumption. This assumption is certainly violated when items are arranged into testlets (i.e., groups of items) that refer to a common item stimulus. In an example of a reading comprehension test, several items could be administered that refer to the same stimulus (i.e., reading text). In the case of cross-country comparisons in an educational assessment that involve the translation of reading texts into several languages, it is likely that DIF is not solely an item but also a testlet (i.e., a reading text) property [69,70].

DIF effects are assumed to be independent across items. This is a strong assumption. It would be useful to provide situations or realistic applications where (approximate) independence is a reasonable condition.

Ultimately, researchers must assess whether DIF is relevant to the construct being measured or not [71,72,73]. When the original HL method or the PHL method is used with a square loss function, all items are included in the linking process, meaning that group differences are based on all items. This approach aligns with the notion that no item should be excluded from group comparisons if DIF is considered construct-relevant. In contrast, when a robust loss function is applied in Haberman linking [42], items with large DIF effects within a group are effectively excluded from comparisons [74]. Such an approach is only appropriate if DIF is deemed construct-irrelevant. Future research could explore linking error estimation for robust HL methods.