2 The BCSM-MI method

The BCSM-MI is introduced by considering the random item parameter model for binary item responses (Fox, Reference Fox2010, Verhagen & Fox, Reference Verhagen and Fox2013). The random item parameters, denoted as $a_{kj}$ and $b_{kj}$ , are assumed to be normally distributed with mean $a_k$ and $b_k$ and variance $\sigma _{a_k}$ and $\sigma _{b_k}$ , respectively. A latent response variable is used to specify the random item parameter model with $Z_{ijk}$ a normally distributed latent response of person i in group j for item k. Then, the random item effects model is represented by

(1) $$ \begin{align} Z_{ijk} &= a_{kj}\theta_{ij} - b_{kj} + e_{ijk}, \\ & = a_k\theta_{ij} - b_k + \epsilon_{a_{k_j}}\theta_{ij} + \epsilon_{b_{k_j}} + e_{ijk}, \nonumber \end{align} $$

where $\epsilon _{a_{k_j}}\theta _{ij}$ and $\epsilon _{b_{k_j}}$ represents an item’s person–group interaction (DIF in discrimination) and an item–group interaction (DIF in difficulty) in group j, respectively, and assumed to be normally distributed with mean zero and variance $\sigma _{a_k}$ and $\sigma _{b_k}$ , respectively. The random measurement error $e_{ijk}$ is assumed to be standard normally distributed. Following Fox et al. (Reference Fox, Koops, Feskens and Beinhauer2020), the BCSM-MI is derived by integrating out the random item parameters. Without conditioning on the group-specific item parameters, they show that the covariance among two latent responses to item k in group j is given by

(2) $$ \begin{align} Cov\left(Z_{ijk},Z_{ljk} \mid \theta_{ij}, a_k, b_k \right) &= \sigma_{a_k}\theta_{ij} \theta_{lj} + \sigma_{b_k}. \end{align} $$

Given this result, the BCSM-MI for responses to item k in group j, including the dependence structure of the DIF effects, can be represented by

(3) $$ \begin{align} \begin{aligned} \mathbf{Z}_{jk} & \sim \mathcal{N}\left(a_k \boldsymbol{\theta}_{j} - b_k, \boldsymbol{\Sigma}_{jk} \right) \\ \boldsymbol{\Sigma}_{jk} &= \mathbf{I}_{n_j} + \sigma_{b_k}\mathbf{J}_{n_j} + \sigma_{a_k}\boldsymbol{\theta}_j\boldsymbol{\theta}^t_j, \end{aligned} \end{align} $$

where $\boldsymbol {\theta }_j \sim \mathcal {N}\left (\mu _{\theta _j}, \sigma _{\theta _j}\mathbf {I}_{n_j} \right )$ and $\mathbf {Z}_{jk}\in \boldsymbol {\Omega }(\mathbf {z}_{jk})$ with $\boldsymbol {\Omega }(\mathbf {z}_{jk}) =\{\mathbf {z}_{jk}, z_{ijk}\le 0 \text { if } y_{ijk}=0, z_{ijk}>0 \text { if } y_{ijk}=1 \}$ representing the allowable region for the latent responses. The $\mathbf {J}_{n_{j}}$ is a matrix of dimension $n_j$ with all elements equal to one and $\mathbf {I}_{n_j}$ is the identity matrix.

To illustrate the modeling of DIF through the structured covariance matrix, the conditional expectation of $Z_{ijk}$ is derived under the BCSM-MI given item and person parameters and remaining latent responses in group j, denoted as $\mathbf {Z}_{(-i)jk}$ . When $\sigma _{b_k}=0$ representing no DIF in difficulty parameter $b_k$ , it is shown that this conditional expectation equals the mean term $\mu _{ijk}=a_k\theta _{ij} -b_k$ plus an item–person–group interaction, $\theta _{ij}\epsilon _{a_{k_j}}$ , according to Equation (1). From Equation (1), with $\epsilon _{b_{k_j}}=0,$ it follows that the ordinary least squares estimator of $\epsilon _{a_{k_j}}$ can be expressed as

(4) $$ \begin{align} \hat{\epsilon}_{a_{k_j}} & = \frac{\sum_{s\ne i} \theta_{sj} \left({Z}_{sjk} - \mu_{sjk} \right)}{\sum_{s \ne i} \theta^2_{sj}}, \end{align} $$

where observation $Z_{ijk}$ is excluded in the estimator for the $\epsilon _{a_{k_j}}$ .

According to the BCSM-MI, the conditional expectation of $Z_{ijk}$ given $\mathbf {Z}_{(-i)jk}$ can be expressed as (see Appendix A)

(5) $$ \begin{align} E\left(Z_{ijk} \mid \mathbf{Z}_{(-i)jk}\right) & = a_k\theta_{ij} - b_k + \theta_{ij} \frac{\sum_{s \ne i} \theta_{sj} \left(Z_{sjk} - \mu_{sjk} \right)}{1/\sigma_{a_{k}} + \sum_{s \ne i} \theta^2_{sj}} \nonumber \\ & = a_k\theta_{ij} - b_k + \theta_{ij} \left(\frac{\sum_{s \ne i} \theta^2_{sj}}{1/\sigma_{a_{k}} + \sum_{s \ne i} \theta^2_{sj}}\right) \hat{\epsilon}_{a_{k_j}}, \end{align} $$

where, in the last term, a Bayes estimator of the $\epsilon _{a_{k_j}}$ can be recognized, in which a positive $\sigma _{a_k}$ arranges shrinkage of the effect toward zero. As a conclusion, under the BCSM-MI, the expected latent response $Z_{ijk}$ depends on the item–person–group interaction, which represents the contribution of DIF in discrimination parameter $a_k$ . When $\sigma _{a_k}=0$ , the covariance support for this interaction effect is not present in the structured covariance matrix in Equation (3). Subsequently, the expression in Equation (5) is only defined for $\sigma _{a_k} \ne 0$ .

The structured covariance matrix needs to be positive definite, which restricts the parameter space of the covariance parameters. When $\sigma _{b_k}=0$ , the determinant of the covariance matrix needs to be positive and implies a restriction on the $\sigma _{a_k}$ . It follows that

(6) $$ \begin{align} det\left( \mathbf{I}_{n_j} + \sigma_{a_k} \boldsymbol{\theta}_j\boldsymbol{\theta}_j\right) & = det\left(\mathbf{I}_{n_j}\right)\left(1+\sigma_{a_k} \boldsymbol{\theta}^t_j\boldsymbol{\theta}_j\right)> 0 \nonumber \\ & = 1+ \sigma_{a_k} \sum_{i}\theta^2_{ij} > 0. \end{align} $$

Subsequently, the covariance matrix is positive definite if $\sigma _{a_k}> -1/\sum _{i}\theta ^2_{ij}$ . The parameter space of the covariance parameter $\sigma _{a_k}$ represents three different types of non-uniform DIF. The common non-uniform DIF in discrimination is represented by $\sigma _{a_k}>0$ , which states a non-uniform difference in item characteristic curves across groups. There is no DIF in discrimination if $\sigma _{a_k}=0$ . When $\sigma _{a_k} \in \left (-1/\sum _{i}\theta ^2_{ij},0\right )$ , there is DIF in discrimination within each group j. This means that for persons in group j, the item discriminates differently and this is referred to as within-group non-uniform DIF. Within each group, the item discriminates more highly to some than other group members, and this is manifested as a negative correlation among the (latent) item responses. Within-group DIF is beyond the scope of this research. Nielsen et al. (Reference Nielsen, Smink and Fox2021) and Fox and Smink (Reference Fox and Smink2023) discuss the modeling of negative associations among clustered measurements and the interpretation of negative clustering effects.

3 Related methods on detecting DIF

The psychometric literature on DIF is huge and shows intensive research in assessing the fairness of measurement instruments. The IRT-based methods condition on the latent variable and are originally designed for a two-group comparison (Camilli & Shepard, Reference Camilli and Shepard1994; Millsap, Reference Millsap2011). A likelihood-ratio (LR) test can be applied to test for differences in item parameters across groups (Thissen et al., Reference Thissen, Steinberg, Wainer, Wainer and Braun1988). These methods require prior knowledge about anchor items assumed to be free of DIF. Results from a stepwise procedure to identify anchor items can be sensitive to the initial choice of anchors. Uncertainty in the model selection process can lead to an incorrect selection of items free of DIF. Due to the stepwise procedure of model selection, the validity of the statistical inference on the DIF parameters is not guaranteed. Disadvantages of DIF methods based on pre-knowledge of anchoring items have been discussed by Chen et al. (Reference Chen, Li, Ouyang and Xu2023).

Recent extensions in IRT-based DIF modeling consider additional group-specific parameters for differences across groups in item discrimination and difficulty for several items. Although an IRT-based LR test can be used to identify variations in item discriminations across groups, due to the increase in the number of model parameters maximum likelihood estimation is more likely to fail. Furthermore, the LR method seriously increases the required sample size in each group for accurate parameter estimation.

The performance of DIF detection methods have not been extensively studied, mainly due to the computational burden. Finch (Reference Finch2005) provides an overview of the performance of several popular methods (e.g., MIMIC, Mantel–Haenszel, and LR test), which require anchor items and they are meant for a two-group comparison. It is shown that the performances depend on the number of anchor items, where more anchors increased the power of DIF detection. Anchors containing DIF negatively affect the performances. Furthermore, the Type-1 error was beyond the nominal level especially for relatively short tests (around 20 items). For longer tests, the MIMIC model showed only a small inflation in the Type-I error. Shih and Wang (Reference Shih and Wang2009) examined the performance of the MIMIC model for DIF detection and showed acceptable (average) Type-I error rates. However, for a group size of 500 and 1,000, the mean power (true positive rate) was around 50%–60% and 80%–90%, respectively, where the difference in difficulty parameters between the reference and focal groups was set at 0.4. It is well known that conditional DIF methods, in which DIF effects are estimated and assessed for significance, require moderate to large sample sizes to perform well (Fox et al., Reference Fox, Koops, Feskens and Beinhauer2020).

Due to page limit, a further discussion of more recently developed DIF methods (e.g., regularized DIF and model-based recursive partitioning) is given in Section S1 of the Supplementary Material. Subsequently, in Section S2 of the Supplementary Material, a comparison in performance between regularized DIF and the BCSM-MI method is given through a simulation study.

The BCSM-MI has several advantages over other DIF methods. For each item parameter, DIF can be assessed for any number of groups with a single covariance parameter. The BCSM-MI method avoids the model selection problem to detect anchor items by modeling simultaneously DIF in all items. The item DIF is simultaneously assessed for significance, which leads to a valid inference procedure. The proposed BF test for a covariance parameter is an omnibus test for (non-)uniform DIF across all groups thereby avoiding multiple hypothesis testing problem. The method does not rely on any anchor items to measure all persons on the same scale. The DIF in each item is modeled through the covariance matrix, which adjusts the person and item parameter estimates but this does not lead to changes in the scale of the person parameters. Therefore, the BCSM-MI can be identified by restricting the mean and variance of the scale of the person parameters (possibly of a reference group), independent of the number of items showing DIF.

4 MCMC estimation of BCSM-MI

The MCMC algorithm for the estimation of the BCSM-MI parameters is described in Appendix C. The sampling of the covariance parameters is done through slice sampling and discussed in more detail in this section. The other model parameters are sampled directly from their conditional distributions given augmented item response data. A Metropolis–Hastings step is defined for the sampling of person parameters, when negative covariance values are sampled (see Appendix D).

The BCSM-MI is identified in the same way as an MG-IRT model. The mean and variance of the person parameters of the reference group are restricted to 0 and 1, respectively. It is also possible to place restrictions on the item difficulty and discrimination parameters to identify the latent scale.

4.1 Posterior distribution (DIF) covariance parameters

To obtain an analytical expression for the posterior distribution of each covariance parameter, a novel method is applied based on eigenvector transformations of the latent item response data. This method generalizes the procedure of Fox (Reference Fox2024), in which one covariance parameter is considered. The method is based on recognizing that the covariance structure of the BCSM-MI in Equation (3) consists of two rank-one components to describe dependence among item responses due to DIF. The two non-zero eigenvalues can be expressed as a function of the covariance parameters. The transformation of the (latent) item response data with the eigenvectors is considered, which has a bivariate normal distribution with a covariance matrix containing the two non-zero eigenvalues. Another set of eigenvectors can be computed which diagonalize this covariance matrix with respect to one of the covariance parameters. Then, an analytical expression of each covariance parameter is obtained by considering the distribution of the latent item response data transformed by both eigenvectors.

By definition, the transformation of the (latent) response data by eigenvectors is sufficient for the covariance parameters. Furthermore, any correlation between the DIF in difficulty and discrimination will not bias the estimation results. The second eigenvector decomposition is an orthogonal transformation with respect to the other covariance parameter thereby properly avoiding the influence of any correlation between the two covariance parameters.

4.1.1 Eigenvector transformations

This procedure starts with the computation of the eigenvectors with non-zero eigenvalues for the two rank one components in the covariance matrix. The following definitions follow directly from theory of a rank-1 matrix. Let $\mathbf {v}_{a_{jk}}=\boldsymbol {\theta }_j/\|\boldsymbol {\theta }_j \|$ be the eigenvector which diagonalize the covariance matrix component $\sigma _{a_k}\boldsymbol {\theta }_j\boldsymbol {\theta }^t_j$ with eigenvalue $\sigma _{a_k}\boldsymbol {\theta }^t_j\boldsymbol {\theta }_j$ , where $\|\boldsymbol {\theta }_j \|$ represents the Frobenius norm of $\boldsymbol {\theta }_j$ . The eigenvector $\mathbf {v}_{b_{jk}}=\mathbf {1}_{n_j}/\|\mathbf {1}_{n_j} \|$ diagonalizes the covariance matrix component $\sigma _{b_k}\mathbf {J}_{n_j}$ with eigenvalue $\sigma _{b_k}n_j$ . The eigenvectors are stored in matrix $\mathbf {v}_{jk}=(\mathbf {v}_{a_{jk}},\mathbf {v}_{b_{jk}})$ . Then, $\tilde {\mathbf {z}}_{jk} = \mathbf {v}^t_{jk}(\mathbf {z}_{jk}-\boldsymbol {\mu }_{jk})$ has a bivariate normal distribution with mean zero and covariance matrix $\mathbf {v}^t_{jk}\boldsymbol {\Sigma }_{jk} \mathbf {v}_{jk}$ .

Consider the projection matrix $\mathcal {P}_a$ and $\mathcal {P}_b$ orthogonal to the column space of $\mathbf {v}^t_{jk}\mathbf {1}_{n_j}$ and $\mathbf {v}^t_{jk}\boldsymbol {\theta }_j$ , respectively, which are given by

$$ \begin{align*} \mathcal{P}_a & = \mathbf{I}_2 - \left(\mathbf{v}^t_{jk}\mathbf{1}_{n_j}\right)\left(\left(\mathbf{v}^t_{jk}\mathbf{1}_{n_j}\right) \left(\mathbf{v}^t_{jk}\mathbf{1}_{n_j}\right)^t\right)^{-1}\left(\mathbf{v}^t_{jk}\mathbf{1}_{n_j}\right)^t, \\ \mathcal{P}_b & = \mathbf{I}_2 - \left(\mathbf{v}^t_{jk}\boldsymbol{\theta}_j\right)\left(\left(\mathbf{v}^t_{jk}\boldsymbol{\theta}_j\right) \left(\mathbf{v}^t_{jk}\boldsymbol{\theta}_j\right)^t\right)^{-1}\left(\mathbf{v}^t_{jk}\boldsymbol{\theta}_j\right)^t, \end{align*} $$

respectively. The projection matrices $\mathcal {P}_a$ and $\mathcal {P}_b$ are each of rank one. Let $\mathbf {w}_{a_{jk}}$ and $\mathbf {w}_{b_{jk}}$ be the eigenvector of projection matrices $\mathcal {P}_a$ and $\mathcal {P}_b$ , respectively, each corresponding to the one non-zero eigenvalue. Then, by definition, $\mathbf {w}_{a_{jk}}$ and $\mathbf {w}_{b_{jk}}$ diagonalize the covariance matrix $\mathbf {v}^t_{jk}\boldsymbol {\Sigma }_{jk}\mathbf {v}_{jk}$ with respect to $\sigma _{a_k}$ and $\sigma _{b_k}$ , respectively. For instance, eigenvectors $\mathbf {w}_{a_{jk}}$ and $\mathbf {v}_{jk}$ diagonalize covariance matrix $\boldsymbol {\Sigma }_{jk}$ with respect to $\sigma _{a_k}$ :

$$ \begin{align*} \mathbf{w}^t_{a_{jk}}\left(\mathbf{v}^t_{jk}\boldsymbol{\Sigma}_{jk} \mathbf{v}_{jk}\right)\mathbf{w}_{a_{jk}} & = \mathbf{w}^t_{a_{jk}}\left(\mathbf{v}^t_{jk}\mathbf{I}_{n_j} \mathbf{v}_{jk} + \sigma_{b_k}\mathbf{v}^t_{jk}\mathbf{J}_{n_j} \mathbf{v}_{jk} + \sigma_{a_k}\mathbf{v}^t_{jk}\boldsymbol{\theta}_j\boldsymbol{\theta}_j^t \mathbf{v}_{jk}\right)\mathbf{w}_{a_{jk}} \\ & = \mathbf{w}^t_{a_{jk}}\mathbf{v}^t_{jk}\mathbf{I}_{n_j} \mathbf{v}_{jk}\mathbf{w}_{a_{jk}} + \sigma_{a_k}\mathbf{w}^t_{a_{jk}}\mathbf{v}^t_{jk}\boldsymbol{\theta}_j\boldsymbol{\theta}_j^t \mathbf{v}_{jk}\mathbf{w}_{a_{jk}}, \\ & = \tilde{\mathbf{w}}^t_{a_{jk}}\tilde{\mathbf{w}}_{a_{jk}} + \sigma_{a_k}\left(\tilde{\mathbf{w}}^t_{a_{jk}}\boldsymbol{\theta}_j\boldsymbol{\theta}_j^t \tilde{\mathbf{w}}_{a_{jk}}\right), \end{align*} $$

where $\tilde {\mathbf {w}}^t_{a_{jk}}=\mathbf {w}^t_{a_{jk}}\mathbf {v}^t_{jk}$ and $\tilde {\mathbf {w}}^t_{a_{jk}}\mathbf {1}_{n_j}=\mathbf {0}$ , since $\mathcal {P}_a$ is orthogonal to the column space of $\mathbf {v}^t_{jk}\mathbf {1}_{n_j}$ .

4.1.2 Posterior distribution of $\sigma _{a_k}$

Consider the distribution of $\tilde {\tilde {z}}_{a_{jk}}=\tilde {\mathbf {w}}^t_{a_{jk}}\tilde {\mathbf {z}}_{jk}$ . Let $s_{a_{jk}}=\tilde {\mathbf {w}}^t_{a_{jk}}\tilde {\mathbf {w}}_{a_{jk}}$ and $\lambda _{a_{jk}}=\tilde {\mathbf {w}}^t_{a_{jk}}\boldsymbol {\theta }_j\boldsymbol {\theta }_j^t\tilde {\mathbf {w}}_{a_{jk}}$ . The distribution of $\tilde {\tilde {\mathbf {z}}}_{a_k}$ given $\sigma _{a_k}$ , $s_{a_{jk}}$ , and $\lambda _{a_{jk}}$ is the product of independently normally distributed variables and is represented by

(7) $$ \begin{align} p\left(\tilde{\tilde{\mathbf{z}}}_{a_k} \mid \sigma_{a_k}, \boldsymbol{\lambda}_{a_k},\mathbf{s}_{a_k} \right) & \propto \prod_j \left(\sigma_{a_k}\lambda_{a_{jk}} + s_{a_{jk}} \right)^{-1/2} \exp\left(\frac{-SS_{a_{jk}}/2}{\sigma_{a_k}\lambda_{a_{jk}} + s_{a_{jk}}}\right), \end{align} $$

where $SS_{a_{jk}}= \left (\tilde {\tilde {z}}_{a_{jk}}\right )^2$ .

There are two restrictions on the parameter space of $\sigma _{a_k}$ . The first restriction states that $\sigma _{a_k}\lambda _{a_{jk}} +s_{a_{jk}}$ is greater than zero for all j, since it represents an eigenvalue of the covariance matrix $\boldsymbol {\Sigma }_{jk}$ . A lower bound for $\sigma _{a_{k}}$ can be derived, and it follows that $LB_1(\sigma _{a_k}) = -\min _j s_{a_{jk}}/\lambda _{a_{jk}}$ .

The second restriction follows from the positive-definite restriction on the covariance matrix $\boldsymbol {\Sigma }_{jk}$ for all j. The determinant of the covariance matrix can be explicitly derived though the matrix determinant lemma (Appendix B). It follows that

(8) $$ \begin{align} LB_2\left(\sigma_{a_k}\right) = -\min_j \left(\boldsymbol{\theta}^t_j\boldsymbol{\theta}_j - \frac{(\sum_{i\mid j} \theta_{ij})^2}{1/\sigma_{b_k}+n_j} \right)^{-1}. \end{align} $$

Then, a shifted-inverse gamma prior can be defined for $\sigma _{a_k}$ , which integrates the lower bound restrictions. Therefore, let $s_{\sigma _{a_k}}=\min \left (-LB_1\left (\sigma _{a_k}\right ), -LB_2\left (\sigma _{a_k}\right ) \right )$ , and

(9) $$ \begin{align} p\left(\sigma_{a_k} \mid s_{\sigma_{a_k}} \right) \propto \left(\sigma_{a_k} + s_{\sigma_{a_k}} \right)^{-g_1-1} \exp\left(-\frac{g_2}{\sigma_{a_k} + s_{\sigma_{a_k}}} \right). \end{align} $$

The posterior distribution of $\sigma _{a_k}$ is proportional to the product of the conditional distribution in Equation (7) and the prior in Equation (9), and cannot be obtained in closed form.

4.1.3 Posterior distribution of $\sigma _{b_k}$

The restrictions on parameter $\sigma _{b_k}$ follow in a similar way. The eigenvalue $\sigma _{b_k}\lambda _{b_{jk}} +s_{b_{jk}}$ is greater than zero for all j, which leads to a lower bound for $\sigma _{b_{k}}$ , $LB_1(\sigma _{b_k}) = -\min _j \frac {s_{b_{jk}}}{\lambda _{b_{jk}}}$ . The determinant of the covariance matrix needs to be greater than zero, which also restricts the parameter space for $\sigma _{b_{k}}$ , and it follows that

(10) $$ \begin{align} LB_2\left(\sigma_{b_k}\right) & = -\min_j \frac{\left(1+\sigma_{a_k}\boldsymbol{\theta}^t_j\boldsymbol{\theta}_j \right)}{n_j\left(1+\sigma_{a_k}\sum_{i\mid j}\left(\theta_{ij}-\overline{\theta}_{j}\right)^2\right)}. \end{align} $$

Then, let $s_{\sigma _{b_k}}=\min \left (-LB_1\left (\sigma _{b_k}\right ), -LB_2\left (\sigma _{b_k}\right ) \right )$ , and a prior for $\sigma _{b_k}$ is a shifted-inverse gamma distribution represented by

(11) $$ \begin{align} p\left(\sigma_{b_k} \mid s_{\sigma_{b_k}} \right) \propto \left(\sigma_{b_k} + s_{\sigma_{b_k}} \right)^{-g_1-1} \exp\left(-\frac{g_2}{\sigma_{b_k} + s_{\sigma_{b_k}}} \right). \end{align} $$

To define the posterior distribution of $\sigma _{b_k}$ , consider the distribution of $\tilde {\tilde {\mathbf {z}}}_{b_k} = \tilde {\mathbf {w}}^t_{b_{jk}}\tilde {\mathbf {z}}_{jk}$ . Let $s_{b_{jk}}=\tilde {\mathbf {w}}^t_{b_{jk}}\tilde {\mathbf {w}}_{b_{jk}}$ and $\lambda _{b_{jk}}=\tilde {\mathbf {w}}^t_{b_{jk}}\boldsymbol {\theta }_j\boldsymbol {\theta }_j^t\tilde {\mathbf {w}}_{b_{jk}}$ . The distribution of $\tilde {\tilde {\mathbf {z}}}_{b_k}$ given $\sigma _{b_k}$ , $s_{b_{jk}}$ , and $\lambda _{b_{jk}}$ is represented by

(12) $$ \begin{align} p\left(\tilde{\tilde{\mathbf{z}}}_{b_k} \mid \sigma_{b_k}, \boldsymbol{\lambda}_{b_k},\mathbf{s}_{b_k} \right) & \propto \prod_j \left(\sigma_{b_k}\lambda_{b_{jk}} + s_{b_{jk}} \right)^{-1/2} \exp\left(\frac{-SS_{b_{jk}}/2}{\sigma_{b_k}\lambda_{b_{jk}} + s_{b_{jk}}}\right), \end{align} $$

where $SS_{b_{jk}}= \left (\tilde {\tilde {z}}_{b_{jk}}\right )^2$ . It follows that the posterior distribution of $\sigma _{b_k}$ is proportional to the product of the prior in Equation (11) and the conditional distribution in Equation (12).

5 Bayes factor tests for non-uniform DIF

The BF is used to assess hypotheses concerning the covariance parameters $\sigma _{a_k}$ and $\sigma _{b_k}$ . For each covariance parameter, there are three hypotheses of interest. The covariance parameter can be equal to zero, for instance, $H_0: \sigma _{a_k}=0$ , which represents no DIF in discrimination for item k. The parameter can be greater than zero, $H_2: \sigma _{a_k}>0$ , which represents DIF in discrimination. A special case can be considered, $H_1: \sigma _{a_k} < 0$ , where DIF is present among group members (instead of DIF between groups). The practical application of this case is not explored in this study. The main purpose of the BF tests is to measure the data support in favor of no (non-uniform) DIF ( $H_0$ ) and DIF ( $H_2$ ). An unrestricted hypothesis $H_u$ is used as a reference to examine the preference of each hypothesis in comparison to the unrestricted hypothesis.

The normalizing constant of the marginal likelihood of the (latent) item response data given the covariance parameter in Equations (7) and (12) is unknown. A parameter transformation is used to avoid the computation of a normalizing constant under each hypothesis. It is shown that the marginal likelihood under each hypothesis for the transformed parameter is computationally much more efficient.

5.1 Change of variables

The distribution of $\tilde {\tilde {\mathbf {z}}}_{a_k}$ given $\sigma _{a_k}$ is a product of independently normally distributed variables with different variance components, but each including $\sigma _{a_k}$ (Equation (7)). A parameter transformation is considered for which the posterior distribution has a known analytical form. The transformation function $\omega _{a_k} = \frac {\sigma _{a_k}\lambda _{a_{jk}}}{\sigma _{a_k}\lambda _{a_{jk}}+s_{a_{jk}}}$ is defined with inverse function and derivative

$$ \begin{align*} \sigma_{a_k}(\omega_{a_k}) & = \frac{\omega_{a_k}}{1-\omega_{a_k}}\left(s_{a_{jk}}/\lambda_{a_{jk}} \right) \\ \frac{d}{d\omega_{a_k}}\sigma_{a_k}(\omega_{a_k}) & = \frac{s_{a_{jk}}}{\lambda_{a_{jk}}(1-\omega_{a_k})^2}, \end{align*} $$

respectively. It can be shown that the distribution of the $\tilde {\tilde {\mathbf {z}}}_{a_k}$ given $\omega _{a_k}$ is a gamma distribution

(13) $$ \begin{align} p\left(\tilde{\tilde{\mathbf{z}}}_{a_k} \mid \omega_{a_k}, \boldsymbol{\lambda}_{a_k},\mathbf{s}_{a_k} \right) & \propto \prod_j \frac{s_{a_{jk}}}{\lambda_{a_{jk}}(1-\omega_{a_k})^2}\left(\frac{\omega_{a_k}}{1-\omega_{a_k}}s_{a_{jk}} + s_{a_{jk}} \right)^{-1/2} \nonumber \\ & \exp\left(\frac{-SS_{a_{jk}}/2}{\frac{\omega_{a_k}}{1-\omega_{a_k}}s_{a_{jk}} + s_{a_{jk}}}\right)\nonumber \\ & \propto \prod_j \left(1 - \omega_{a_k}\right)^{1/2} \exp\left(-\left(1-\omega_{a_k}\right)SS_{a_{jk}}/(2s_{a_{jk}}) \right) \nonumber \\ & = \left(1 - \omega_{a_k}\right)^{J/2} \exp\left(-\left(1-\omega_{a_k}\right)\sum_{j}SS_{a_{jk}}/(2s_{a_{jk}}) \right) \nonumber \\ & = \tilde{\omega}_{a_k}^{J/2} \exp\left(-\tilde{\omega}_{a_k}\sum_{j}SS_{a_{jk}}/(2s_{a_{jk}}) \right), \end{align} $$

where $\tilde {\omega }_{a_k}=1-\omega _{a_k}$ , and assume a gamma prior for $\tilde {\omega }_{a_k}$ . Then, $\tilde {\omega }_{a_k}$ has a gamma posterior distribution given $\tilde {\tilde {\mathbf {z}}}_{a_k}$ with shape parameter $g_1+J/2$ and scale parameter $g_2+\sum _j SS_{a_{jk}}/(2s_{a_{jk}})$ . The lower bound for $\sigma _{a_k}$ is translated to a boundary value for $\omega _{a_k}$ , denoted as $UB(\tilde {\omega }_{a_k})$ .

In the same way, it can be shown with transformation function $\omega _{b_k} = \frac {\sigma _{b_k}\lambda _{b_{jk}}}{\sigma _{b_k}\lambda _{b_{jk}}+s_{b_{jk}}}$ that $\tilde {\omega }_{b_k}=1-\omega _{b_k}$ has a gamma distribution. Then, with a gamma prior for the $\tilde {\omega }_{b_k}$ , the posterior distribution of $\tilde {\omega }_{b_k}$ is a gamma distribution with shape parameter $g_1+J/2$ and scale parameter $g_2+\sum _j SS_{b_{jk}}/(2s_{b_{jk}})$ , where $SS_{b_{jk}} = (\tilde {\tilde {z}}_{b_{jk}})^2$ . An upper bound for $\tilde {\omega }_{b_k}$ , $UB(\tilde {\omega }_{b_k})$ , follows from the lower bound restriction on $\sigma _{b_k}$

5.2 Bayes factor testing

The hypotheses for $\sigma _{a_k}$ can be translated to those for $\tilde {\omega }_{a_k}$ , taking into account that the upper bound $UB\left (\tilde {\omega }_{a_k}\right )$ is greater than one ( $LB\left (\sigma _{a_k}\right )$ is less than zero). For hypothesis $H_2 (\sigma _{a_k}>0)$ , $\tilde {\omega }_{a_k} \in (-\infty ,1)$ , for $H_1 (\sigma _{a_k}<0)$ , $\tilde {\omega }_{a_k} \in (1,UB\left (\tilde {\omega }_{a_k}\right )$ , and for $H_0 (\sigma _{a_k}=0), \tilde {\omega }_{a_k}=1$ . For the unrestricted hypothesis, let $\pi _u\left (\tilde {\omega }_{a_k}\right )$ be the (unrestricted) gamma prior for $\tilde {\omega }_{a_k}$ with shape $g_1$ and scale parameter $g_2$ . Following the encompassing prior approach (Böing-Messing & Mulder, Reference Böing-Messing and Mulder2020; Klugkist & Hoijtink, Reference Klugkist and Hoijtink2007), for each hypothesis $H_t$ for $t=0,1,2$ , the prior under $H_t$ is a restricted version of the encompassing gamma prior, $\pi _t\left (\tilde {\omega }_{a_k}\right ) = \pi _u\left (\tilde {\omega }_{a_k}\right )\chi \left (\tilde {\omega }_{a_k} \mid H_t\right )$ , where

$$ \begin{align*} \chi\left(\tilde{\omega}_{a_k} \mid H_t\right)=\frac{\mathbb{I}\left(\tilde{\omega}_{a_k} \in H_t\right)}{P\left(\tilde{\omega}_{a_k} \mid H_t\right)}. \end{align*} $$

The BF for the restricted hypothesis $H_t$ for item k concerning $\tilde {\omega }_{a_k}$ against the unrestricted hypothesis $H_u$ can be expressed as

(14) $$ \begin{align} \begin{aligned} BF_{tu}(\tilde{\omega}_{a_k}) & = \int_{\mathbb{Z}_k} \left[\int_{-\infty}^{UB(\tilde{\omega}_{a_k})} p\left(\tilde{\omega}_{a_k} \mid \tilde{\tilde{\mathbf{z}}}_{k} \right) \chi\left(\tilde{\omega}_{a_k} \mid H_t\right) d\tilde{\omega}_{a_k} \right] p\left(\tilde{\tilde{\mathbf{z}}}_{k} \mid \mathbf{y}_k \right) d\tilde{\tilde{\mathbf{z}}}_{k} \\ & = E\left[\frac{P\left(\tilde{\omega}_{a_k} \in H_t, \tilde{\omega}_{a_k} < UB(\sigma_{a_k}) \mid \tilde{\tilde{\mathbf{z}}}_{k} \right)/P\left(\tilde{\omega}_{a_k} < UB(\sigma_{a_k}) \mid \tilde{\tilde{\mathbf{z}}}_{k} \right)}{P\left(\tilde{\omega}_{a_k} \in H_t, \tilde{\omega}_{a_k} < UB(\sigma_{a_k}) \right)/P\left(\tilde{\omega}_{a_k} < UB(\sigma_{a_k}) \right)}\mid \mathbf{y}_k \right], \end{aligned} \end{align} $$

where $\mathbb {Z}_k$ represents the support for $\tilde {\tilde {\mathbf {z}}}_{k}$ . The BF in Equation (14) is expressed as the expected value of a ratio of gamma cumulative distribution functions, which leads to an easy computation of the BF. Furthermore, the BF for two constrained hypotheses also follows easily, since $BF_{tt'}=BF_{tu}/BF_{t'u}$ for $t,t'=0,1,2$ . The logarithm of the BF (log-BF) can be computed, where a log- $BF_{tt'}$ of 3–5 is considered strong evidence and above 5 very strong (decisive) evidence in favor of hypothesis t (Kass & Raftery, Reference Kass and Raftery1995).

Let $\pi _u\left (\tilde {\omega }_{b_k}\right )$ denote an unrestricted gamma prior for $\tilde {\omega }_{b_k}$ with shape $g_1$ and scale parameter $g_2$ . Then, for each hypothesis, the prior under $H_t$ is a restriction on the encompassing gamma prior. Subsequently, the BFs for examining hypotheses about $\sigma _{b_k}$ follow in a similar way and can also be expressed as a ratio of gamma cumulative distribution functions, while using the appropriate upper bound for $\tilde {\omega }_{b_k}$ . The corresponding hypotheses for $\tilde {\omega }_{b_k}$ are given by $H_2$ ( $\sigma _{b_k}>0$ ), $\tilde {\omega }_{b_k} \in (-\infty ,1)$ , $H_1 (\sigma _{b_k}<0), \tilde {\omega }_{b_k} \in (1,UB(\tilde {\omega }_{b_k}))$ , and $H_0 (\sigma _{b_k}=0)$ , $\tilde {\omega }_{b_k}=1$ .

The BF for covariance parameters can be applied to any number of groups and for any number of items. It provides information about which hypothesis is more likely without requiring that any of the hypotheses is true. The BF for testing a precise null hypothesis ( $H_0$ ) against an unconstrained alternative hypothesis ( $H_u$ ) is not favoring the null regardless of the information in the data (known as Bartlett’s paradox), since the prior probability of $\tilde {\omega }_{b_k}=1$ is not sensitive to changes in the variance of the gamma prior.

6 Simulation study

For each data replication, cluster sizes were sampled from a uniform distribution between ( $n-10$ , $n+10$ ) with n the average cluster size and $n_j$ the size of cluster j. The discrimination and difficulty parameters were sampled from a log-normal and normal distribution with mean zero and standard deviation 0.2 and 0.5, respectively. Subject parameters were sampled from a standard normal distribution. Interest was in the covariance parameters, although (latent) cluster means were also estimated. True covariance values for $\sigma _{a_{k}}$ and $\sigma _{b_{k}}$ were fixed across data replications and are given in the first column of Table 1. The covariances values were set at $-$ 0.002 and increased with a step size of 0.02 from item 1 to 10. A positive covariance value represents variance across group-specific item parameters. The grid of covariance values reflects no DIF, low, to a high level of DIF across item parameters. Thus, for instance, item 3, the $\sigma _{a_{k}}=0.038$ represents a range of $a_k \pm 0.40$ for the group-specific discriminations. For item 10, it is approximately $a_k \pm 0.84$ . The grid of covariance values, specifying the level of DIF, was equal for the discrimination and difficulty parameters.

Table 1

Simulation study: Covariance estimation results for different sample sizes averaged across 1,000 data replications

Table 1 Long description

The table is divided into three horizontal sections based on sample size parameters. The columns are Item (1 to 10), True sigma a = sigma b, and two sets of results for sigma a and sigma b, each containing Median, S E, and B F sub 0 2.

* Section 1: J = 10, n = 200. True values range from -0.002 to 0.178. For sigma a, B F sub 0 2 starts at 6.939 for Item 1 and decreases to -114.539 for Item 10. For sigma b, B F sub 0 2 starts at 8.314 and decreases to -142.483.

* Section 2: J = 5, n = 200. True values range from -0.002 to 0.178. For sigma a, B F sub 0 2 starts at 6.917 and decreases to -40.667. For sigma b, B F sub 0 2 starts at 5.594 and decreases to -71.162.

* Section 3: J = 2, n = 250. True values range from -0.002 to 0.178. For sigma a, B F sub 0 2 starts at 7.230 and decreases to -1.440. For sigma b, B F sub 0 2 starts at 7.130 and decreases to -30.487.

S E represents the Monte Carlo standard error, and B F sub 0 2 represents the log-B F of no D I F versus between-cluster D I F.

Note: SE represents the Monte Carlo standard error, and $BF_{02}$ represents the log-BF of no DIF versus (between-cluster) DIF.

The BCSM-MI was identified by restricting the scale of the person parameters to have a mean of zero and a variance of one. The prior for each covariance parameter is a shifted-inverse gamma, in Equations (9) and (11), with shape and rate parameter equal to 0.001. The hyper prior parameter settings were equal to those specified in Fox (Reference Fox2024, Appendix C).

Three conditions were considered with 10 items, and J = 10, 5, and 2 clusters, and with average cluster sizes n = 200, 200, and 250, respectively. A total of 1,000 data replications were made in each condition. Average results are reported in Table 1. The MCMC algorithm was ran for 5,000 iterations with a burn-in period of 1,000 iterations. The effective sample size was at least 400 for each parameter for each data set. The log-BFs were computed of the null hypothesis $H_0$ : $\sigma _{a_k}=0$ (no DIF) versus the alternative hypothesis $H_2$ : $\sigma _{a_k}> 0$ (DIF), represented as $BF_{02}$ , for all items k. A similar log-BF was computed for $\sigma _{b_k}$ for all items k.

7 Real data study

Data from the Programme for International Student Assessment or PISA was analyzed for DIF between countries and mode of administration. A total of 14 math items of booklet B_2015_52 was administered through a computer-based assessment (CBA) in PISA 2015 (OECD, 2016). The same set of items was administered within “standard” paper-based assessment (PBA) in 2012 in booklet B_2012_4. The item responses from Belgium (BEL), Germany (DEU), and The Netherlands (NED) were considered. The response patterns of 249 (BEL), 148 (DEU), and 183 (NED) respondents were collected from a PBA in 2012, and of 677 (BEL), 338 (DEU), and 351 (NED) respondents from a CBA in 2015.

An overview of assessment and score differences between the 2012 and 2015 cycle is given by Jerrim (Reference Jerrim2016). It was concluded that pupils who did the computer-based test under performed to those who did the paper-based test. The mode of test administration had an impact on the pupils’ scores. For several items, there was a mode effect in which, for instance, items appeared to be less difficult in a PBA than in a CBA. The BCSM-MI method was used to examine country-specific item–mode interactions and country–item-specific person–mode interactions, while accounting for latent group means and variances. The BCSM-MI did not require any anchor items, which enabled a simultaneous DIF analysis of the 14 items across the three selected countries, each with two modes of administration.

In Table 2, for each item parameter, the covariance value, SE, and the level of $BF_{02}$ (no DIF versus between-cluster DIF) are given. For item discrimination 10, the $BF_{02}$ shows decisive evidence in favor of DIF. Although the $BF_{02}$ is also negative for item discrimination 7, it shows only weak evidence in favor of DIF. For item difficulties 1, 2, and 12, there is also decisive evidence in favor of DIF. For the item parameters associated with DIF, the size of the $BF_{02}$ is associated with the level of covariance, where more negative $BF_{02}$ values correspond to higher covariance values. According to the BF results, those items (rows) shown in light gray have a negative $BF_{12}$ for at least one item parameter, and item parameters in dark gray show DIF.

Table 2

PISA 2012–2015: DIF examination of paper-based and computed-based assessed math items across three countries

Table 2 Long description

The table is divided into two main sections: Discrimination (sigma sub a) and Difficulty (sigma sub b), each listing items 1 through 14.

Columns include:

* Item number

* Median

* S E (Monte Carlo standard error)

* B F sub 0 2

* L R (Likelihood-ratio test)

* P (L R) (P-value L R test)

* L R sub seq (sequential final L R test)

* P (L R sub seq) (P-value sequential final L R test)

* M D (mean deviation)

* R M S D (root mean square deviation)

Key data points in Discrimination section:

* Item 1: Median 0.011, S E 0.002, B F sub 0 2 8.104, L R 5.576, P (L R) 0.350.

* Item 10: Median 0.115, S E 0.007, B F sub 0 2 minus 33.694, L R 34.445, P (L R) 0.000, L R sub seq 34.221, P (L R sub seq) 0.000.

Key data points in Difficulty section:

* Item 1: Median 0.111, S E 0.007, B F sub 0 2 minus 34.109, L R 67.065, P (L R) 0.000, L R sub seq 67.065, P (L R sub seq) 0.000, M D 0.141, R M S D 0.120.

* Item 12: Median 0.060, S E 0.003, B F sub 0 2 minus 11.973, L R 31.068, P (L R) 0.000, L R sub seq 31.068, P (L R sub seq) 0.000, M D 0.086, R M S D 0.045.

Note: SE represents the Monte Carlo standard error. $LR$ : Likelihood-ratio test. $P(LR)$ : P-value LR test. $LR_{seq}$ (sequential) final LR test. $P(LR_{seq})$ : P-value (sequential) final LR test, MD: mean deviation, and RMSD: root mean square deviation.

The items for which the $BF_{12}$ was positive for discrimination and difficulty were considered invariant and used as anchors (i.e., free of DIF) in an MG-IRT analysis. A baseline MG-IRT model was defined using the anchors for an LR test procedure (Millsap, Reference Millsap2011). The latent mean and variance parameters were considered free parameters, given a reference group for which the latent mean and variance were restricted to 0 and 1, respectively. Items 1, 2, 7, 10, and 12 were examined for DIF in discrimination and difficulty, which are colored light gray in Table 2. The R-package mirt (Chalmers, Reference Chalmers2012) was used to perform LR testing for DIF by comparing the baseline model to various constrained MG-IRT models, in which item parameters were constrained to be equal across groups.

For each item examined for DIF, an MG-IRT model with the discrimination parameter constrained was compared to the baseline model using the LR test. In a similar way, MG-IRT models were defined for which each of the item difficulty parameters was restricted. Thus, a total of 10 models were needed to examine each item parameter of the five examined items. In Table 2, the LR values and corresponding P-values for the discriminations and difficulties are given under the header LR and P(LR), respectively. High LR values correspond to negative $BF_{02}$ values and small P-values. There appears to be an agreement between the methods about which item parameters show DIF. However, the difficulty parameter of item 7 appears to show DIF according to the LR test, while the BF test shows weak support for no DIF. For the discrimination parameter of item 7, the $BF_{02}$ shows weak evidence for DIF, while the LR test does not. However, note that the log-BF test is a simultaneous DIF analysis of all items and is not based on a comparison(s) with a baseline model.

Subsequently, DIF in discriminations and difficulties was each sequentially tested under the MG-IRT model. In this sequential procedure, items were restricted to be invariant when satisfying the AIC criterion, and the remaining set of items was examined until no more invariant items were found. In Table 2, for the final set of item parameters that displayed DIF the LR statistic and P-value are reported, which are denoted as $LR_{seq}$ and $P(LR_{seq})$ , respectively. For the discriminations, items 10 and 12 are marked as DIF, and for the difficulties, items 1, 2, 7, and 12. However, when using the sample-adjusted BIC as a selection criterion to decide on which item is restricted to be invariant, the LR sequential test results changed. In that case, only the discrimination parameter of item 10 is still marked to show DIF, while the discrimination of item 12 is considered to be free of DIF. This is in line with the BCSM-MI analysis. Furthermore, only the difficulties of items 1 and 12 show DIF, and difficulties of items 2 and 7 do not show DIF. This is not in line with the BCSM-MI results, where the difficulty of item 2 also showed DIF.

The mean deviation (MD) and the root mean square deviation (RMSD) were computed for each item (Köhler et al., Reference Köhler, Hartig, Khorramdel and Pokropek2024, Equations (2) and (3)). In a baseline MG-IRT model, all item parameters were restricted across groups to estimate the common set of item parameters. Subsequently, group-specific item parameters were estimated using the anchors from the BCSM-MI analysis. In Table 2, for the examined set of items, the maximum MD and RMSD across groups are represented in the rows of the difficulty parameters. The MD is considered to assess the difference in item difficulty between the group specific item response function and the common item response function, where the RMSD measures this (squared) difference with respect to difficulty and discrimination. The estimated levels of MD are in correspondence with the BF. Items 1, 2, and 12 have relatively high MD values corresponding with high BF values for the difficulty parameters. Item 7 has a relatively high MD value for which the RMSD value is also relatively high, which corresponds to the weak support of the BF for DIF in discrimination.

In the PISA study of Köhler et al. (Reference Köhler, Hartig, Khorramdel and Pokropek2024), under a two-parameter IRT model, an item’s response function is considered to differ across groups, when the RMSD and the absolute value of MD are greater than 0.12. According to that criterion, only item 1 is considered to show DIF. It could be that for large sample sizes comparable to PISA studies, this threshold leads to acceptable Type-1 errors. Obviously, a high threshold will reduce the number of false positives. However, for item 1, an estimated difficulty covariance of 0.111 indicates that the range of difficulty values across groups is around $\pm 0.67$ . An RMSD value greater than 0.15 to identify items behaving differently across groups was used in von Davier et al. (Reference von Davier, Yamamoto, Shin, Chen, Khorramdel, Weeks, Davis, Kong and Kandathil2019). When using this threshold, none of the items would be classified to show DIF. Note that the simulation study showed that for five groups with a group size of around 200, a discrimination covariance of 0.06, and a difficulty covariance of 0.06, already led to strong evidence in favor of DIF. Therefore, the item response function of item 10 differs across groups with respect to the discrimination. The item response functions of items 1, 2, and 12 also differ across groups with respect to the item difficulty. For these item response functions, the log-BF shows very strong evidence in favor of DIF, and the very high LR values also lead to rejecting the null hypothesis of a constrained (invariant) item parameter.