university of maryland, college park

††footnotetext: Correspondence should be sent to Weimeng Wang
E-Mail<EMAIL_ADDRESS>
Affiliation: University of Maryland, College Park

###### Abstract

Recent advancements in testing differential item functioning (DIF) have
greatly relaxed restrictions made by the conventional multiple group item
response theory (IRT) model with respect to the number of grouping variables
and the assumption of predefined DIF-free anchor items. The application of the
$L_{1}$ penalty in DIF detection has shown promising results in identifying a
DIF item without a priori knowledge on anchor items while allowing the
simultaneous investigation of multiple grouping variables. The least absolute
shrinkage and selection operator (LASSO) is added directly to the loss
function to encourage variable sparsity such that DIF parameters of anchor
items are penalized to be zero. Therefore, no predefined anchor items are
needed. However, DIF detection using LASSO requires a non-trivial model
selection consistency assumption and is difficult to draw statistical
inference. Given the importance of identifying DIF items in test development,
this study aims to apply the decorrelated score test Ning Liu (2017) to test
DIF once the penalized method is used. Unlike the existing regularized DIF
method which is unable to test the statistical significance of a DIF item
selected by LASSO, the decorrelated score test requires weaker assumptions and
is able to provide asymptotically valid inference to test DIF. Additionally,
the deccorrelated score function can be used to construct asymptotically
unbiased normal and efficient DIF parameter estimates via a one-step
correction. The performance of the proposed decorrelated score test and the
one-step estimator are evaluated via a Monte Carlo simulation study.

  * Key words: Differential item functioning, Moderated nonlinear factor analysis, LASSO, $L_{1}$ penalty, LASSO, Score test, Penalized EM algorithm

Abstract

##  1 Introduction

A major goal of educational and psychological assessment is to create
interpretable scores on a relevant construct being measured. A well-
constructed test score representing a latent trait (e.g., proficiency level),
should be reliable and valid. One major threat to the validity of the intended
use of scores from an instrument is the lack of measurement invariance or
existence of differential item functioning (DIF). Sometimes, characteristics
of a test may introduce unintended systematic score differences between
individual test takers or subgroups of examinees of the same ability level,
which can result in invalid score interpretation. Therefore, DIF detection
should be conducted for routine operations of psychological and educational
assessment to ensure valid interpretations of test scores.

Item response theory (IRT), a collection of mathematical and statistical
models, was—and still remains—the foundational modeling system in educational
assessment Embretson Reise (2000). One advantage of an IRT model in conducting
a DIF analysis is that it offers a formal statistical framework for assessing
DIF. Despite the popularity and advantages of DIF detection using IRT models,
the inherent assumptions therein make practical applications of DIF analyses
challenging. For example, IRT DIF detection methods usually explicitly or
implicitly assume the existence of predefined DIF-free items111See
wangliuliu2022 for an exception., also known as anchor items. Typically, DIF-
free items are needed to anchor the latent scale so that DIF can be
distinguished from the between group difference in latent trait distributions.
Correctly specifying anchor items is crucial in correctly identifying a DIF
item. Violations to the DIF-free anchor assumption may lead to inflated false
detection rates in finding DIF items WC. Wang (2004); WC. Wang Yeh (2003);
Woods (2009); Woods . (2013). A practical obstacle is locating DIF-free anchor
items without a priori information. Another challenge is detecting DIF items
associated with multiple grouping variables. Conventional IRT DIF detection
methods typically investigate DIF with respect to one categorical grouping
variable (e.g., gender, ethnicity, SES) at a time. However, it is impossible
to investigate the complex nature of DIF due to interconnected background
characteristics M. Liu (2017) and thus it is unlikely to truly uncover the
source of DIF Shea (2013).

As a solution to these practical difficulties encountered by conventional IRT
DIF methods, regularization has been applied to the latent regression modeling
framework to detect DIF items. In particular, the use of the $L_{1}$ penalty
or the LASSO has shown promising results in DIF detection due to its
efficiency in variable selection <e.g.,
¿bauer2017more,belzak2020improving,magis2015detection,tutz2015penalty. For
example, the $L_{1}$ penalty has been successfully applied to moderated
nonlinear factor analysis (MNLFA) model to identify DIF items Bauer . (2020);
Belzak Bauer (2020). Such a modeling framework offers greater flexibility to
detect DIF for multiple grouping variables simultaneously. Item parameters and
population distribution parameters, upon applying suitable link functions, can
be expressed as linear or nonlinear functions of multiple person
characteristics. Consequently, DIF detection can be treated as a variable
selection problem—non-zero coefficients are only assigned to covariates that
cause DIF. The goal of the penalty procedure is to encourage sparseness as
much as possible with respect to DIF parameters such that DIF parameters of
anchor items are penalized to be zero while those of DIF items are non-zero.
Hence, no predefined anchor items are needed as items with DIF parameters
penalized to be zero are used as anchors.

Despite the flexibility of the $L_{1}$ penalty for DIF detection, several
concerns need to be addressed carefully. First, the accuracy of DIF detection
based on the $L_{1}$ penalty is contingent on non-trivial assumptions. In
order for LASSO to select the right covariate for each item, and thus, items
exhibiting true DIF, a variable selection consistency condition is needed Zhao
Yu (2006). In the current case, failing to meet this condition could result in
mistakingly flagging non-DIF items or miss real DIF items even in large
samples. Second, with a finite sample size especially those encountered in the
social sciences, inferential statistics such as confidence interval estimates
and p-values are extremely important to differentiate a true DIF item from a
sampling error. For example, a common practice in DIF detection in educational
and psychological assessments is effect size reporting after DIF is detected.
In this case, an uncertainty measure may prove to be helpful in deciding
whether a DIF effect size is truly different from zero. Nevertheless, it is
difficult to draw statistical inference based on the LASSO-type estimator due
to shrinkage resulting from penalization Belzak Bauer (2020); Huang (2018);
Lindstrøm Dahl (2020); Tutz Schauberger (2015). Lastly, LASSO creates biased
parameter estimates and non-normal limiting distribution of the parameter
estimatesFu Knight (2000), which makes subsequent analyses difficult Chen .
(2021); Fan Li (2001). For example, DIF effect sizes, item parameter
estimates, and standard error estimates cannot be directly used to quantify
uncertainty in latent score estimates Y. Liu Yang (2018).

##  2 Purpose of the Study

The current study aims to fill the gap in the literature to make inferential
claims once LASSO is used. Specifically, a decorrelated score test Ning Liu
(2017) is proposed to detect DIF for binary response data with multiple
covariates after the $L_{1}$ penalty is used. The decorrelated score test
could potentially be extended to accommodate continuous and discrete response
data and multidimensional latent variables. In the current study, the simplest
and the most commonly used case was considered to better demonstrate the
extension of the decorrelated score test to detect DIF. Unlike the existing
regularized DIF method, the decorrelated score test does not require the
variable selection consistency assumption and is able to provide valid
inference on DIF effects. Specifically, a sparse score vector with respect to
the focal parameter is estimated consistently so that the resulting score test
statistic has an asymptotically normal reference distribution. Additionally,
an asymptotic unbiased estimator can be constructed using a one-step bias
correction using the decorrelated score function.

A Monte Carlo simulation study is conducted to examinee the finite sample
behavior of the decorrelated scodre test. Its performance in controlling the
Type I error rate, establishing sufficient power, controlling the false
detection rate in identifying a DIF item under different DIF-related
conditions are compared with three methods: (1) regularization method based on
LASSO selection only, (2) a naive model refitting method Belzak Bauer (2020),
and (3) the oracle solution assuming known anchors.

##  3 Decorrelated Score Test to Test DIF

Mathematically, DIF can be defined as follows. Considering $n$ independent and
identically distributed multivariate vectors,
$\mathfrak{Y}=(\mathbf{Y}_{1},\cdots,\mathbf{Y}_{n})^{\top}$ and
$\mathbf{Y}_{i}=(Y_{i1},\dots,Y_{iJ})^{\top}$, following a statistical model
$\mathcal{P}=\\{P:\bfxi\in\Xi\\}$, where is a $d-$dimensional vector of
unknown parameters and $\Xi\subset\mathcal{R}^{d}$ is the parameter space. The
categorical random variable $Y_{ij}$ representing the item response from
person $i$, $i=1,\cdots,n$ to item $j\in\mathcal{J}=\\{1,\dots,J\\}$, define
the conditional item response function, denoted as
$f_{j}(y|\theta,\mathbf{x})=P(Y_{ij}=y|\theta_{i}=\theta,\mathbf{x}_{i}=\mathbf{x})$,
as the probability of endorsing a particular answer conditional on the
person’s latent ability level $\theta_{i}=\theta\in\mathcal{R}$ and person
covariate $\mathbf{x}_{i}\in\mathcal{R}^{K}$. More broadly, item $j$ exhibits
DIF if

$\displaystyle f_{j}(y|\theta,\mathbf{x})\neq
f_{j}(y|\theta,\tilde{\mathbf{x}}),$ (1)

for some $\mathbf{x}\neq\tilde{\mathbf{x}}\in\mathcal{R}^{K}$. Equation 1
indicates that DIF exists when the conditional item response function $f_{j}$
differs for any $\mathbf{x}\neq\tilde{\mathbf{x}}$ after controlling for the
same level of latent ability. In other words, measurement invariance holds for
item $j$, only when
$f_{j}(y|\theta,\mathbf{x})=f_{j}(y|\theta,\tilde{\mathbf{x}})$ for all
$\mathbf{x}\neq\tilde{\mathbf{x}}$.

###  3.1 Moderated Nonlinear Factor Analysis Model

The moderated nonlinear factor analysis model (MNLFA) assumes that the item
response for item $j$ from person $i$, (i.e., $Y_{ij}\in\\{0,1\\}$) follows a
Bernoulli distribution
($Y_{ij}|\theta_{i}\sim\hbox{Bern}(P_{ij}(\theta_{i}))$), where $P_{ij}$
denotes the probability of endorsing item $j$ for person $i$. Incorporating
vector of person covariates $\mathbf{x}_{i}$ into the conditional item
response function (cIRF) to affect both the item response and the latent
distribution offers greater flexibility in testing DIF across different levels
of both continuous and categorical variables. It also permits testing DIF for
multiple grouping variables as well as their interaction effects (e.g., gender
$\times$ age) at the same time. Using the same notation, the cIRF $f_{j}$ is
written as

$\displaystyle f_{j}(Y_{ij}|\theta,\mathbf{x})$
$\displaystyle=P_{ij}^{Y_{ij}}(1-P_{ij})^{1-Y_{ij}},\,\hbox{where}$
$\displaystyle P_{ij}=P(Y_{ij}$
$\displaystyle=1|\theta_{i}=\theta,\mathbf{x}_{i}=\mathbf{x})=\frac{1}{1+\exp[-\alpha_{j}(\mathbf{x})-\beta_{j}(\mathbf{x})\theta]},$
(2)

where $\alpha_{j}(\mathbf{x})$ and $\beta_{j}(\mathbf{x})$ are the item
intercept and item slope functions, respectively, which can be further
expressed as functions of the vector of person covariate
$\mathbf{x}\in\mathcal{R}^{K}$. Latent ability $\theta_{i}$ is assumed to be
normally distributed:
$\theta_{i}\sim\mathcal{N}(\mu(\mathbf{x}),\sigma^{2}(\mathbf{x}))$, in which
the latent population parameters: mean ($\mu$) and variance ($\sigma^{2}$) can
also be expressed as functions of the person covariate vector $\mathbf{x}$
defined as follows:

$\displaystyle\mu(\mathbf{x})={}^{\top}\mathbf{x}\quad\sigma^{2}(\mathbf{x})=\exp({}^{\top}\mathbf{x}),$
(3)

where $K$-dimensional vectors and represent effects of the person covariate
vector $\mathbf{x}$ on the latent mean and variance on a logarithmic scale.
These two vectors are further referred to as latent population parameters.
Also, item intercept and item slope functions are expressed as

$\displaystyle\alpha_{j}(\mathbf{x})$
$\displaystyle=d_{j}+{}_{0j}^{\top}\mathbf{x},$ (4)
$\displaystyle\beta_{j}(\mathbf{x})$
$\displaystyle=a_{j}+{}_{1j}^{\top}\mathbf{x},$ (5)

where $d_{j}$ and $a_{j}$ are the item location and item slope parameters. The
$K$-dimensional vectors 0j and 1j stand for person covariate effects on the
item intercept and item slope parameter, respectively. To detect DIF for item
$j$, a composite hypothesis test is typically conducted to test against the
null hypothesis

$\displaystyle\hbox{H}_{0}:\,{}_{0j}={}_{1j}=\bf{0}.$ (6)

Typically, a minimal constraint is needed to identify the MNLFA model. This
requires that there exists at least one item for each column of 0j and 1j is
zero/anchor. Without predefined anchor items or constrained parameters, the
MNLFA model can be estimated using the Penalized Expectation Maximization (EM)
Algorithm and DIF items are those items with remaining DIF parameters.

###  3.2 Regularized Differential Item Functioning

The marginal likelihood function $f(\mathbf{y}_{i}|\mathbf{x}_{i})$ of an
observed individual response vector
$\mathbf{y}_{i}=(y_{i1},\dots,y_{1J})^{\top}$ assuming conditional
independence among item responses given the latent variable and person
covariates $\mathbf{x}_{i}=(x_{i1},\dots,x_{iK})^{\top}$ is expressed as

$\displaystyle
f(\mathbf{y}_{i}|\mathbf{x}_{i})=\int\prod_{j=1}^{J}f_{j}(y_{ij}|\theta_{i},\mathbf{x}_{i})\phi(\theta_{i}|\mathbf{x}_{i})d\theta_{i},$
(7)

in which $\phi(\cdot)$ is the probably density function of a normal
distribution governed by the population parameters and . Assuming that the
item response vectors $\mathfrak{Y}$ is independent and identically
distributed (i.i.d), the sample likelihood function is written as

$\displaystyle
f_{n}(\mathbf{Y}|\mathbf{X})=\prod_{i=1}^{n}f(\mathbf{y}_{i}|\mathbf{x}_{i}),$
(8)

where $\mathbf{Y}=(\mathbf{y}_{1},\cdots,\mathbf{y}_{n})^{\top}$ and
$\mathbf{X}=(\mathbf{x}_{1},\cdots,\mathbf{x}_{n})^{\top}$ represent observed
data matrices. Subsequently, the corresponding marginal log-likelihood is
expressed as

$\displaystyle\ell_{n}(\bfxi;\mathbf{Y}|\mathbf{X})=\frac{1}{n}\log
f_{n}(\mathbf{Y}|\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}\log
f(\mathbf{y}_{i}|\mathbf{x}_{i}),$ (9)

where
$\bfxi=(a_{1},d_{1},{}_{01}^{\top},{}_{11}^{\top},\dots,a_{J},d_{J},{}_{0J}^{\top},{}_{1J}^{\top},{}^{\top},{}^{\top})^{\top}$
is a collection of model parameters. To estimate the model using the penalized
EM algorithm, the $L_{1}$ penalty function is added directly to the loss
function

$\displaystyle
p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})=-\ell_{n}(\bfxi;\mathbf{Y}|\mathbf{X})+\lambda\|{}_{j}\|_{1},$
(10)

where ${}_{j}=({}_{0j}^{\top},{}_{1j}^{\top})^{\top}$ denotes all DIF
parameters, $\ell_{n}(\bfxi;\mathbf{Y}|\mathbf{X})$ is the sample log-
likelihood expressed in Equation 9 and $\|\|_{1}$ indicates the $L_{1}$ norm
that sums absolute values of all elements weighted by the penalty weight
$\lambda$. Then the local minimizer of
$p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})$ with respect to

$\displaystyle\hat{\bfxi}$
$\displaystyle=\operatorname*{\arg\min}p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})$
(11)

can be estimated by the penalized EM algorithm, which is documented in
Appendix A.

In the upcoming section, the decorrelated score test will be introduced. The
intuition about the decorrelated score test resembles the efficient score test
when the model is identified. Thus, the efficient score test for an identified
model will first be introduced followed by the decorrelated score test when
the model is not identified with redundant parameters.

###  3.3 Decorrelated Score Test

Partition model parameters into $\bfxi=(\bfpsi,{}^{\top})^{\top}$, where is
the $d_{0}-$dimensional vector of focal parameters and ⊤ stands for
$d_{1}=(d-d_{0})-$dimensional vector of nuisance parameters. Given the loss
function $\ell(\bfxi,\mathbf{Y})=-\ell_{n}(\bfxi,\mathbf{Y})$, which is the
average sample log-likelihood, define
$\boldsymbol{\mathcal{I}}=\operatorname{\mathbb{E}}(\nabla^{2}\ell(\bfxi,\mathbf{Y}))$.
Let ∗ be the true value of . Similarly, denote
$\boldsymbol{\mathcal{I}}^{*}=\operatorname{\mathbb{E}}_{{}^{*}}(\nabla^{2}\ell({}^{*},\mathbf{Y}))$.
Note that for the rest of the study, $d$ is fixed. Let’s first consider the
case when the model is identified. For example, an MNLFA model with
sufficiently large subset of known anchors. The asymptotic normality of the ML
estimator [$\hat{\bfxi}=\hbox{argmin}_{\bfxi\in\Xi}\ell(\bfxi,\mathbf{Y})$]
follows from a first-order Taylor series expansion of the score function at ∗,

$\displaystyle\mathbb{0}$
$\displaystyle=\nabla\ell(\hat{\bfxi},\mathbf{Y})=\nabla\ell({{}^{*}},\mathbf{Y})+\nabla^{2}_{\bfxi,\bfxi}\ell({{}^{*}},\mathbf{Y})(\hat{\bfxi}-{}^{*})+R_{n}$

where the remainder term is of the form
$R_{n}=\frac{1}{2}(\hat{\bfxi}-\bar{\bfxi})^{T}\nabla_{\bfxi,\bfxi}^{3}\ell(\bar{\bfxi},\mathbf{Y})(\hat{\bfxi}-\bar{\bfxi})$
and $\bar{\bfxi}$ is in between $\hat{\bfxi}$ and ∗. By rearranging this
equation,

$\displaystyle\sqrt{n}\,(\hat{\bfxi}-{}^{*})$
$\displaystyle=-\sqrt{n}\,\boldsymbol{\mathcal{I}}^{*-1}\nabla\ell({{}^{*}},\mathbf{Y})+\sqrt{n}R_{n}$
$\displaystyle\Leftrightarrow\sqrt{n}\,\begin{pmatrix}\hat{\bfpsi}-{}^{*}\\\
\hat{\bfeta}-{}^{*}\end{pmatrix}$
$\displaystyle=-\sqrt{n}\,\begin{pmatrix}\boldsymbol{\mathcal{I}}^{*}_{\bfpsi,\bfpsi}&\boldsymbol{\mathcal{I}}^{*}_{\bfpsi,\bfeta}\\\
\boldsymbol{\mathcal{I}}^{*}_{\bfeta,\bfpsi}&\boldsymbol{\mathcal{I}}^{*}_{\bfeta,\bfeta}\end{pmatrix}^{-1}\begin{pmatrix}\nabla\ell({{}^{*}},\mathbf{Y})\\\
\nabla\ell({{}^{*}},\mathbf{Y})\end{pmatrix}-\sqrt{n}R_{n},$

where subscripts of $\boldsymbol{\mathcal{I}}^{*}$ are corresponding
partitions of the matrix. Drop $\mathbf{Y}$ from the rest of the notations. As
$\boldsymbol{\mathcal{I}}_{\bfeta,\bfeta}$ is invertible, using the block
inverse formula, it follows that

$\displaystyle\sqrt{n}\,(\hat{\bfpsi}-{}^{*})=-\sqrt{n}\,\boldsymbol{\mathcal{I}}^{*-1}_{\bfpsi|\bfeta}s({{}^{*}})+\sqrt{n}R_{n},\quad
R_{n}=o_{p}(1/\sqrt{n})$

where
$\mathbf{s}({{}^{*}})=\nabla\ell({}^{*})-\boldsymbol{\mathcal{I}}^{*}\boldsymbol{\mathcal{I}}^{*-1}\nabla\ell({}^{*})$
and
$\boldsymbol{\mathcal{I}}^{*}_{\bfpsi|\bfeta}=\boldsymbol{\mathcal{I}}^{*}-\boldsymbol{\mathcal{I}}^{*}\boldsymbol{\mathcal{I}}^{*-1}\boldsymbol{\mathcal{I}}^{*}$
are the efficient score and efficient information, respectively. Note that
efficient score can be interpreted as the projection of the
$\nabla\ell(\bfxi)$ to the orthogonal complement of the score function with
respect to the nuisance parameters Vaart (1998). Under
$H_{0}:{}^{*}=\mathbf{0}$, it holds that
$n\nabla(\mathbf{0},\hat{\bfeta})^{\top}[\hat{\boldsymbol{\mathcal{I}}}_{\bfpsi|\bfeta}]^{-1}\nabla(\mathbf{0},\hat{\bfeta})\overset{\mathcal{D}}{\rightarrow}\chi^{2}_{d_{0}}$,
where $\hat{\bfeta}=\hbox{argmin}\ell(\mathbf{0},\bfeta)$ is constrained
parameter estimates. It is straightforward to see that the estimated efficient
score ($\mathbf{s}(\mathbf{0},\hat{\bfeta})$) equals
$\nabla(\mathbf{0},\hat{\bfeta})$ due to the fact that
$\nabla\ell(\mathbf{0},\hat{\bfeta})=\mathbf{0}$.

However, when the model is not identified,
$\boldsymbol{\mathcal{I}}_{\bfeta,\bfeta}$, in general, is not invertible.
Intuitively, there exists redundant parameters in the model. A natural
extension of the above case when the model is not identified is estimating an
“efficient score” so that the influence of entries of the redundant parameters
in the nuisance score is minimized. Following the same logic as the efficient
score, a sparse score can be estimated by projecting $\nabla\ell(\bfxi)$ to
the orthogonal complement of a low-dimensional subspace spanned by the non-
redundant part of the nuisance score vector. Therefore, a sparse vector/matrix
is needed to find the best sparse linear combination of $\nabla\ell(\bfxi)$ to
approximate $\nabla\ell(\bfxi)$. Mathematically, the decorrelated score, or
the extension of the efficient score when the model is not identified, has no
correlation with the score function with respect to the nuisance score (i.e.,
$\operatorname{\mathbb{E}}(\mathbf{s}(\bfxi)^{\top}\nabla\ell(\bfxi))=\mathbf{0}$).
Geometrically, the decorrelated score has the same interpretation which is the
projection of $\nabla\ell(\bfxi)$ to the orthogonal complement of the linear
space spanned by the nuisance score function $\nabla\ell(\bfxi)$. Following
this intuition, define the decorrelated score function as

$\displaystyle\mathbf{s}(\bfxi)$
$\displaystyle=\nabla\ell(\bfxi)-\mathbf{W}^{\top}\nabla\ell(\bfxi),\hbox{
where
}\mathbf{W}^{\top}=\boldsymbol{\mathcal{I}}\boldsymbol{\mathcal{I}}^{-1}\in\mathcal{R}^{d_{0}\times
d_{1}}.$ (12)

To find the projection of $\nabla\ell(\bfxi)$ to the orthogonal complement of
the linear space span by the nuisance score function $\nabla\ell(\bfxi)$, the
sparse matrix $\mathbf{W}$ can be estimated using Algorithm 1 (see Table 1).
As can be seen, the key is to estimate a sparse matrix
$\hat{\mathbf{W}}=(\hat{\mathbf{W}}_{*1},\cdots,\hat{\mathbf{W}}_{*d_{0}})$
column by column to construct the decorrelated score
($\mathbf{s}(\bfpsi,\bfeta)$) so that the additional and redundant parameters
do not influence the decorrelated score.

=========================

Insert Table 1 about here

=========================

Further the Dscore test statistic can be constructed

$\displaystyle\hat{T}_{Dscore}=n\hat{\mathbf{s}}(\mathbf{0},\bfeta)^{\top}[\hat{\boldsymbol{\mathcal{I}}}_{\bfpsi|\bfeta}]^{-1}\hat{\mathbf{s}}(\mathbf{0},\bfeta),$
(13)

where
$\hat{\boldsymbol{\mathcal{I}}}_{\bfpsi|\bfeta}=\nabla^{2}_{\bfpsi,\bfpsi}\ell(\hat{\bfxi})-\hat{\mathbf{W}}^{\top}\nabla^{2}_{\bfeta,\bfpsi}\ell(\hat{\bfxi})$.
Under some technical assumptions222These assumptions are documented in the
Appendix Appendix B and are not yet verified in the latent variable models.
The intuition of the assumption verification is discussed in the Appendix
Appendix B, it can be proved that

$\displaystyle\hat{T}_{Dscore}\xrightarrow{\cal{D}}\chi^{2}_{d_{0}}.$ (14)

Moreover, a one-step asymptotical unbiased estimator $\tilde{\bfpsi}$ can be
constructed using a single Newton step as

$\displaystyle\tilde{\bfpsi}=\hat{\bfpsi}-\hat{\boldsymbol{\mathcal{I}}}^{-1}_{\bfpsi|\bfeta}\hat{\mathbf{s}}(\hat{\bfxi})$
(15)

and with assumptions 1 to 4 documented in Appendix Appendix B, ning2017general
showed that

$\displaystyle\sqrt{n}(\tilde{\bfpsi}-{}^{*})\xrightarrow{\cal{D}}{{\cal
N}(0,\boldsymbol{\mathcal{I}}_{\bfpsi|\bfeta}^{*-1})}.$ (16)

Subsequently, the $(1-\alpha)\times 100\%$ confidence interval can be
constructed for a linear combination of ∗ (i.e., $\mathbf{c}^{\top}{}^{*}$
where $\mathbf{c}^{\top}$ is a $d_{0}-$dimensional constant vector) as
$[\tilde{\bfpsi}-n^{-1/2}\Phi^{-1}(1-\alpha/2)(\mathbf{c}^{\top}\hat{\boldsymbol{\mathcal{I}}}_{\bfpsi|\bfeta}\mathbf{c})^{-1/2},\tilde{\bfpsi}+n^{-1/2}\Phi^{-1}(1-\alpha/2)(\mathbf{c}^{\top}\hat{\boldsymbol{\mathcal{I}}}_{\bfpsi|\bfeta}\mathbf{c})^{-1/2}]$.

##  4 Monte Carlo Simulation

In order to investigate the finite sample behavior of the proposed Dscore test
in testing DIF under different conditions, a Monte Carlo study is conducted.
Performance of the proposed Dscore test is evaluated in comparison to three
methods: (1) the Wald test assuming known and correctly specified anchors, (2)
the Reg-DIF method (i.e., results based on the LASSO selection only), and (3)
and the naive model refitting approach Belzak Bauer (2020). The simulation
study design and evaluation criteria are discussed in Section 4.1 and Section
4.4, respectively. All computations are performed in R R Core Team (2020). The
source code will be made publicly available through the Open Science Framework
(OSF) upon acceptance of the manuscript.

###  4.1 Study Design

Binary response data are generated from a unidimensional MNLFA model under two
conditions—with or without DIF items. Two factors are manipulated including
(1) the total sample size ($n=500,1,000,\&\ 2,500$) and (2) the percentage of
DIF items ($0\%,\,25\%,\,\&\ 50\%$). As a result, there are $3\times 3=9$
fully crossed conditions. These manipulating factors are selected due to their
relevance to the DIF detection mechanism of the Dscore test. Moreover, actual
values of the manipulated factors were chosen to conform to real-world data
analytic scenarios and align with previous methodological studies
<e.g.,¿bauer2017more,belzak2020improving. Due to the limited DIF research on
the topic, model parameters were generated based on a real data analysis using
the UK normative sample data of the Revised Eysenck Personality Questionnaire
<EPQ-R, ¿eysenck1985revised333We are grateful to Dr. Paul Barrett for granting
us access to the data. True data generating values of these DIF items are
described in details in the following section.

###  4.2 The True Data Generation Model

The binary response for examinee $i$ for item $j$, denoted
$Y_{ij}\in\\{0,1\\}$, $i=1,\dots,n$, $j=1,\dots,J$, was generated from a MNLFA
model (see Equation 3.1). The total number of items was fixed at $J=12$. Model
parameters of the data generating model are tabulated in Table 2.
Specifically, item discrimination and intercept parameters are the ML
estimates of item parameters using the EPQ-R data with five items identified
as anchors. Three person covariates are considered in the study which mimic
gender, age, and their product. Specially, the dichotomous grouping variable
is generated from a Bernoulli distribution with success probability of $.5$
(i.e., $x\sim\hbox{Bern}(0.5)$). Then, for those in category 1, age is
generated from $\mathcal{N}(0.2,1)$. Otherwise, age is generated from
$\mathcal{N}(0,1)$, which creates a correlation between age and gender of 0.1.
An interaction effect is then created by multiplying age and gender. By
including an interaction effect, the correlations between the interaction
effect and the two variables is large ($>.6$). Furthermore, three types of DIF
items were generated. Under the $25\%$ DIF condition, items 1 to 3 are DIF
items representing items with large, medium, and small DIF effect sizes. When
the proportion of DIF items is increased to $50\%$, items 4 to 6 are added to
the DIF set and their DIF effect sizes replicate those of items 1-3.

=========================

Insert Table 2 about here

=========================

To visualize the generated items, Figure 1 displays the true probability of
each person endorsing items with the condition of $50\%$ DIF items and
$n=2,500$. Plots on the left correspond to DIF items whereas those on the
right visualize non-DIF items. Gender influence on the item response can be
seen from different colors and shapes.

=========================

Insert Figure 1 about here

=========================

Lastly, the latent variable follows a normal distribution ${\cal
N}({}^{\top}\mathbf{x},{}^{\top}\mathbf{x})$ conditional on $\mathbf{x}$, and
similarly true population parameter values are shown in Table 2, which are
also enlightened by real data. A total of 500 replications for each condition
are implemented ensure that the 95% normal approximation confidence band at
the nominal level $0.05$ to be $[0.031,0.069]$.

###  4.3 Estimation

For each replication under each condition, the following four methods are
fitted to the binary response data to test DIF at the item level and the
parameter level. MNLFA models are estimated using penalized maximum likelihood
for the Reg-DIF method and the Dscore method or ML for the refit method and
oracle solution. For both estimation methods, integration with respect to the
latent variable is approximated by a 49-point Gauss–Hermite quadrature. The
four methods are discussed in detail here.

First, each item is tested for DIF based on LASSO selection only (i.e., the
Reg-DIF method). The MNLFA model fitted with all three grouping covariates is
estimated using the penalized EM algorithm. To select the optimal tuning
parameter under each replication, typically BIC is used in latent variable
models. However, the current study fixed $\lambda$ at an estimated optimal
value (i.e., $c\sqrt{1/N}$ where the constant $c$ is estimated by
$\hat{c}=\sqrt{n}\sum^{R}_{r=1}(\hat{\lambda}_{r})/R$)444A pilot study found
that BIC calculated based on the penalized EM parameter estimates always
select a relatively smaller $\lambda$ value when the DIF effect size is large,
which results in too many FDR. This phenomena persists even when sample size
is large $n=2,500$ or $n=5,000$. We suspect that the less optimal $\lambda$
selection is due to parameter bias due to the regularization. This can be
verified that when the $\lambda$ value is selected by the BIC value calculated
based on the maximum likelihood estimates (i.e., the model refitting method),
the FDR rate is more controlled especially under the large sample size
condition.. The conditions with $n=500$ were used to estimate the constant.
Results show that $\hat{c}=0.8291,0.6883,0.5727$ when the number of DIF items
is 0, 3, and 6, respectively. Therefore, $\lambda$ is fixed at the estimated
optimal value for each condition. For the reg-DIF method, an item does not
exhibit DIF if both 0j and 1j are zero vectors. Otherwise, the item is
considered as a DIF item.

Next, the naive model refitting method can be applied by refitting the same
MNLFA model with the anchors selected using method 1. This model is estimated
using the marginal ML estimation method with the EM algorithm. The marginal
likelihood function can be approximated using the same configuration as in the
penalized EM algorithm. The convergence tolerance for the log-likelihood
change for the penalized EM algorithm and the M-step are set to be $10^{-4}$
and $10^{-6}$, respectively. The maximum number of iterations is set to be
$500$. Additionally, the Dscore test is conducted after the initial run
(method 1) using the penalized EM algorithm. Details of the Dscore test was
described in Section 3.3 and thus are not repeated here. A critical step to
conduct the Dscore test is to estimate $\hat{\mathbf{W}}$ as displayed in
Equation A16. As the theory only requires that $\lambda$ and
$\lambda^{\prime}$ are of the same rate, $\lambda^{\prime}$ will be set to the
selected $\lambda$ value in the initial parameter estimation step to speed up
the computation and save time. ning2017general and fang2017testing have found
that the decorrelated score test is not sensitive to $\lambda^{\prime}$ and
both fixed $\lambda^{\prime}=0.5\sqrt{\log d/n}$ in their simulation studies.
Furthermore, asymptotic unbiased parameter estimates can be estimated using
the one-step bias correction (see Equation 15). Note that although the focal
parameters can be multidimensional or unidimensional (i.e., the debias step
can be conducted at the item level or parameter level), the current study
investigates the one-step bias correction by parameter type. For instance, for
each item, the bias correction treats $(a_{j},d_{j},{}_{0j},{}_{1j})^{\top}$
as the focal parameters and everything else as nuisance parameters. The
population parameter estimates are corrected at the same time by treating all
item parameters as nuisance parameters.

Lastly, to compare the performance of the Dscore test with the oracle solution
assuming anchors are known, Wald tests are performed to test DIF at the item-
level. For each replication and each condition, item 11 and item 12 are
treated as anchors. In addition, final model parameter estimates and their
corresponding standard errors are estimated using anchors selected by the Wald
test.

###  4.4 Evaluation Criteria

The comparative inferential performance of the four methods is evaluated in
terms of (1) hypothesis testing in testing a DIF item measured by rejection
rates at the $0.05$ level, (2) parameter recovery including bias and variance
of model parameter estimates including a-DIF and d-DIF, and (3) recovery of
standard errors. Standard error estimates are evaluated by comparing the
square root of the mean of the variance of parameter estimates against the
empirical standard errors (i.e., Monte Carlo standard deviation of the
parameter estimates).

##  5 Results

###  5.1 Results: Hypothesis Testing

In this section, the Type I error rate, false positive rate (FDR), and power
of detecting a DIF item are reported. To visualize the comparative
performance, the empirical rejection rate at the nominal level $0.05$ was
plotted under (1) the null condition when there were no DIF items (i.e., Type
I error), (2) the alternative condition when there was a mix of DIF and DIF-
free items for DIF items only (i.e., Power), and (3) the alternative condition
when there was a mix of DIF and DIF-free items for DIF-free items only (i.e.,
FDR).

####  5.1.1 Type I Error Rate

Figure 2 shows the Type I error rate of detecting a DIF item under the null
condition when there were no DIF items. The horizontal dashed lines show the
$0.05$ nominal level and the dotted lines represent the 95% normal
approximation confidence band at the nominal level $0.05$ across 500
replications. If the rejection rate falls within the 95% confidence band
(i.e., 95% confidence interval =$[0.031,0.069]$), the specific method is
considered to have well-controlled Type I error rate of incorrectly detecting
a DIF item. Larger values than the upper bound of the confidence band at the
nominal level 0.05 represent over-rejecting the null hypothesis whereas
smaller values than the lower bound of the confidence band indicate under-
rejecting the null hypothesis.

Overall, the Reg-DIF, model refit method, and decorrelated score test showed
controlled Type I error under all manipulated sample size conditions as the
rejection rate of all these methods fall within the 95% confidence band.

=========================

Insert Figure 2 about here

=========================

####  5.1.2 Power

Figure 3 shows the power of detecting DIF items under different conditions.
The power can be negatively impacted by the sample size while more robust to
the number of true DIF items. Overall, the decorrelated score test was more
powerful than the Reg-DIF and the model refit methods across all conditions.
The difference was more obvious under more challenging conditions where the
sample size is small. The performance of the Reg-DIF method and the model
refit method was nearly identical when the total number of DIF items was
small. The difference between the two methods was more obvious when the number
of DIF items was large. Specifically, the model refit method was slightly less
powerful than Reg-DIF by design. The reason is that once an item is identified
as an anchor item, its DIF parameters are excluded in the refitted model. In
other words, this item can no longer be identified as a DIF item.

=========================

Insert Figure 3 about here

=========================

####  5.1.3  False Detection Rate

Finally, results for incorrectly detecting an anchor item as a DIF item can be
visualized using the empirical rejection rate for the true anchor item under
alternative conditions. As displayed in Figure 4, the performance of the
decorrelated score test is outstanding as the FDR of the decorrelated score
test falls within the 95% normal-approximation confidence band, indicating
controlled FDR. Conversely, the Reg-DIF and model refit method under-rejected
the null hypothesis when the number of DIF items is small. This is
understandable, as was previously shown that the power of identifying the true
DIF item was in general smaller than the decorrelated score test and the
oracle solution. However, when the number of DIF items was large, the Reg-DIF
method over-rejected the null hypothesis. As was previously illustrated, the
model refit method is more conservative in rejecting a null hypothesis as
compared to the Reg-DIF method. Intuitively, it makes sense that the FDR is
more controlled as compared to the Reg-DIF method. Nevertheless, it is still
overly rejecting the null hypothesis when the sample size is not sufficiently
large. As the sample size increases, the FDR using the model refit method
falls within the 95% confidence band at the nominal 0.05 level.

=========================

Insert Figure 4 about here

=========================

###  5.2 Results: Parameter Recovery

To investigate parameter recovery, the bias and variance of the parameter
estimates using each method were computed. Discussion is focussed on the DIF
parameters, the parameters that inferences will be drawn on. Bias and variance
of other model parameter estimates including item slope ($a_{j}$), item
intercept ($d_{j}$), and population parameters (,) using different methods are
included in the supplementary materials. For all methods other than the one-
step debiased estimator using the decorrelated score function, DIF parameters
penalized to be 0 or not included in the final model are treated as 0 to
compute bias and variance of an estimator.

####  5.2.1 Bias

When there was no DIF, all methods recover model parameters well with bias
generally less than $0.05$. This is understandable because when there was no
DIF, the Type I error rate, as was previously shown, was well-controlled. The
penalized EM performed similarly as the EM since in most cases because
$\lambda$ was large enough to penalize all DIF parameters to be 0. Even if
this does not happen, the model refit method is always fitting the correct
model under the null condition.

When there was a mix of DIF and non-DIF items, the performance of different
methods varied. The Reg-DIF method often resulted in biased DIF parameters for
the true DIF items (i.e., items 1 to 3 in the 3 DIF item condition and item 1
to 6 in the 6 DIF item condition) due to the shrinkage. The bias for the DIF
parameters for the anchor items were relatively small. Compared with the Reg-
DIF method, the model refit method reduced bias for the true DIF items by
$3\%$ to $92\%$ depending on the item. Note that bias remained relatively
large for the a-DIF parameter even sample size was large (see last row of the
Figure 5) due to failing to select the right DIF-effect. It was more clear
when the sampling distribution of non-zero DIF parameters was investigated
(see Figure 7). As can be seen, the sampling distribution of the non-zero DIF
effect (blue line) using the model refit method was bi-modal with one mode at
0 indicating that the effect was not selected by the LASSO and another non-
zero mode.

As a comparison, the one-step debiased parameter estimate using the
decorrelated score function performed remarkably well. DIF parameters were
recovered well as the bias ranged from -0.03 to 0.06 across conditions. Its
advantages in reducing bias due to shrinkage were more obvious as compared
with the model refit method under smaller sample size or more DIF items
conditions. Under the most difficult condition (i.e., small sample size and
large number of DIF items), it performed equally well as compared to the
oracle solution. More importantly, as the sample size increased, the bias
decreased.

=========================

Insert Figure 5 about here

=========================

=========================

Insert Figure 6 about here

=========================

=========================

Insert Figure 7 about here

=========================

####  5.2.2 Variance of DIF Parameter Estimates

Figure 8 and Figure 9 present the variance of the DIF parameter estimates.
Large values indicate more uncertainty of the parameter estimates while
smaller values indicate less uncertainty. As the sample size increases,
variance should, in general, decrease. Among all methods, the Reg-DIF method
produced the least variable parameter estimates. The model refit method
produced similarly variable a-DIF parameter estimates but notably more
variable d-DIF parameter estimates as compared to the Reg-DIF especially for
the DIF items. Lastly, the one-step debiased parameter estimate using the
decorrelated score function was the most variable estimator for the a-DIF
parameter as compared to all other methods especially when the sample size was
small. However, it was less variable for the d-DIF parameter for the DIF item
as compared with the model refit method.

=========================

Insert Figure 8 about here

=========================

=========================

Insert Figure 9 about here

=========================

###  5.3 Results: Standard Error Estimates

The recovery of SE was calculated as the ratio between the square root of the
mean of the variance of parameter estimates over the empirical standard
deviation of the Monte Carlo parameter estimates. If the ratio is closer to 1,
the method has a more valid uncertainty measure. One caveat when calculating
the SE using the model refit method is that there is no valid SE estimates for
DIF parameters penalized to be 0. If the Reg-DIF method does not select a
specific DIF effect for a specific replication, the corresponding DIF
covariate is excluded from the model refit method. In such cases, SE estimates
of the unselected DIF effect is treated as 0. These arbitrary 0 SEs would
artificially decrease the average of the estimated SEs across replications and
thus should be interpreted with caution. Alternatively, omitting DIF effects
penalized to be 0 will often exclude too many DIF effects across replications
resulting in instability of the empirical SEs and inaccurate recovery SE
measure555Recovery of SEs measure for the model refit method based on non-zero
DIF effects only has strange large values. Similar findings are reported from
chen2021advantages.

The SE recovery for the a-DIF and d-DIF parameters can be visualized in Figure
10 to 11. As there are no valid SE estimates from the penalized EM alogrithm,
the Reg-DIF method currently presented in the following figures are based on
the SE estimates using the model refit method. Note that as some DIF
parameters of certain items were never selected by the Reg-DIF method leaving
the recovery of SEs unavailable for the model refit method. This is actually
one of the disadvantages of using the model refit method that inference can
only be drawn on the selected DIF parameters. Again, the model refit method
tended to underestimate SEs of DIF parameters. The underestimation of the SEs
of DIF effects could contribute to incorrectly detecting an anchor item as a
DIF item as was shown previously shown in Figure 4. As a comparison, the
decorrelated score method not only can recover SEs of non-zero DIF effects but
also zero DIF effects. For a-DIF and d-DIF, SEs ranged from 0.95 to 1.08 and
0.90 to 1.02, respectively. Note that it appears that SEs of a-DIF and d-DIF
effects for item 7 were underestimated. However, given that variances of a-DIF
and d-DIF parameters for item 7 were strangely large due to the skewness of
the sampling distribution, it makes sense that SEs of DIF parameters of item 7
seem to be underestimated.

=========================

Insert Figure 10 about here

=========================

=========================

Insert Figure 11 about here

=========================

##  6 Discussion and Conclusion

###  6.1 Summary

One of the fundamental issues and remaining challenges inherent in DIF
detection is finding the correct anchor items. Although recent development in
the DIF detection literature (such as using regularization in the MNLFA
modeling framework) has shown promising results in detecting DIF items without
predefined anchor items, issues such as inflated false detection rate and
inability to draw valid inference still remain.

The goal of the current study is to apply the decorrelated score test to test
DIF items based on the $L_{1}-$ penalized maximum likelihood solution under
the MNLFA modeling framework. Unlike DIF detection based on regularization
only, the decorrelated score test is valid for all DIF and DIF-free items
across all covariates. Additionally, it has been shown that the decorrelated
score function can be further used to construct an asymptotically unbiased
estimator. Furthermore, the simulation study has shown the comparative
performance of the decorrelated score test, DIF detection using regularization
only, the model refit method, and the Wald test assuming the correct anchor
items in hypothesis testing and parameter recovery.

Overall, the decorrelated score is advantageous than the other two anchor-free
DIF detection methods in three aspects. First, the decorrelated score test
shows more statistical power in detecting a DIF item while maintaining
controlled false detection rate even compared with the oracle solution with
two anchor items available. As summarized and highlighted by previous studies
Belzak Bauer (2020); Jacobucci . (2016), the Reg-DIF method and the model
refit method often result in inflated false detection rate in identifying a
true DIF item especially when the sample size is small. Similar to these
studies, the current simulation study also shows inflated false detection rate
especially under the most difficult yet realistic conditions (i.e., when the
sample size is small and the number of DIF items is large). As was previously
discussed, the current simulation study also shows that the model refit method
has similar performance as the decorrelated score test in correctly
identifying a true DIF item while controlling for the Type I error when the
sample size is large.

Another important finding from the current simulation study is that the one-
step debiased estimator using the decorrelated score function yields less
biased item and DIF parameter estimates. Specifically, as compared with the
model refit method, the one-step debiased estimator recovers the item slopes,
item intercepts, a-DIF, and d-DIF parameters well and sometimes even
comparable to the oracle solution. As expected, the model refit method reduces
the bias for model parameters estimated from the penalized EM algorithm.
However, bias still remains. This is expected as LASSO seems to be more
sensitive to detecting d-DIF rather than the a-DIF, a point that we return to
at the end of this section, it is likely that a-DIF parameters are likely to
be mistakenly unselected. Accordingly, the model refit method fits an
incorrect model and, therefore, bias remains. Given the less biased item and
DIF parameter estimates, it is recommended to use the one-step debiased
estimator to produce reliable point estimates without fitting an additional
MNFLA model using EM. However, less bias often indicates more parameter
variability. This is evidenced by the fact that there are more uncertainty in
the item parameter and DIF parameter estimates using the one-step debiased
estimator based on the decorrelated score function. Another interesting
observation is that consistent with the model refit method, the one-step
debiased estimator always recover item intercept parameters better as compared
to item slope parameters (see the supplementary materials). This finding is
also aligned other studies <see¿[for examples]bauer2020simplifying,
belzak2021using that d-DIF parameters are in general recovered better than
a-DIF parameters, which could potentially explain the superiority of the item
intercept parameter recovery. Consequently, a-DIF and item slope parameters
have better initial parameter estimates to be plugged into the bias correction
function.

Lastly, the decorrelated score test is able to draw reliable statistical
inference on a-DIF and d-DIF parameters for zero and non-zero effects,
respectively, where the model refit method may fail or is not able to generate
one when the DIF effect is zero. As discussed earlier, SEs cannot be computed
if DIF effects are penalized to be zero. To this end, the one-step debiased
estimator is the only method considered in the present study that can produce
reliable statistical inference for the zero DIF effects. As for the remaining
DIF effects that are not penalized to be zero, the model refit method often
underestimate the SEs of the DIF effects Chen . (2021); Huang (2018) leading
to incorrect inferences. This may be of concern, if substantive researchers
are interested in testing a-DIF and d-DIF separately. If so, the decorrelated
score test offers promising and valid statistical test to identify DIF at both
the item and the parameter levels.

##  7 Future Studies

Although the proposed decorrelated score test and the one-step debiased
estimator have shown promising results in detecting DIF items without
predefined anchor items and providing valid uncertainty measures for DIF
effects, limitations and issues still remain to be addressed by future
studies.

First, the simulation study has suggested that the one-step debiased estimator
does not seem to reduce the bias of the penalized ML estimator and that the
associated SEs are unreliably estimated using the efficient information. Given
that the ill recovery with the uncertainty measure also happens with
unpenalized item parameters, it is conjectured that additional penalty weights
might be needed for the population parameter or a more careful selection of
penalty weights is needed. Some preliminary investigation suggests that if the
population parameter is penalized in the initial $L_{1}$ penalization stage,
the bias correction using the decorrelated score function shows better
performance. Future studies are needed to examine the behavior of the
decorrelated score function when population parameters are penalized.
Alternatively, finding the optimal $\lambda^{\prime}$ for the population
parameter for each iteration might also improve the unpredicted behavior of
SEs. In the current simulation study, $\lambda^{\prime}$ is set to the fixed
$\lambda$ based on the estimated rate as the theory suggests
$\lambda\asymp\sqrt{1/n}$ and $\lambda^{\prime}\asymp\sqrt{1/n}$ are
approximately the same rate to ensure the $l_{2}$ error bound for initial
parameter estimates. It might worth to explore the influence of
$\lambda^{\prime}$ on the SE recovery of the unpenalized parameter. Although
ning2017general and fang2017testing have reported that the performance of the
decorrelated score test is insensitive to the actual $\lambda^{\prime}$ value
in the high-dimensional regression model and proportional hazard model, it is
not clear whether $\lambda^{\prime}$ will impact the one-step debiased
estimator especially when there are unpenalized parameters in the model.

Second, as $\lambda$ is critically important for the initial parameter
estimates and the accuracy of DIF detection using the Reg-DIF method and the
model refit method, future research is encouraged to investigate different
model selection critera in addition to BIC. As mentioned in Section 4, a pilot
study found that BIC calculated based on the penalized EM parameter estimates
always selects a relatively smaller $\lambda$ value when the DIF effect size
is relatively large and thus results in too many false positives. The current
study used a fixed $\lambda$ value estimated from the rate calculated from the
$n=500$ condition based on the ML estimator which potentially avoids the issue
of poor performance of BIC. This may also explains the relatively smaller
inflated false detection rate and power as compared with other Reg-DIF studies
<e.g.,¿bauer2020simplifying,belzak2020improving. Given the sparse literature
on the model selection accuracy in the regularized DIF framework, future
research is needed to find viable and efficient penalty selection approaches.

Third, future research is needed to investigate the impact of different
penalty functions on the performance of the decorrelated score test in DIF
detection and the asymptotic behavior of the debiased parameter estimator. The
general theory provided in ning2017general can be directly applied to many
penalty functions <e.g., adaptive LASSO,¿[]zou2006adaptive. Our simulation
study along with the other studies <e.g.,¿belzak2020improving showed similar
differential performance in recovering different parameters and SEs. In
addition, applying different penalty weights for a-DIF and d-DIF may yield
better initial parameter estimates and thus improve the performance of the
one-step debias estimator.

Lastly, DIF effect size measures can be critically important when deciding to
remove, modify, or keep a flagged DIF item. In practice it might be more
desirable or more realistic to flag items with large DIF impact, which makes
effect size reporting more essential in DIF detection. Latent regression
models such as MNLFA models provide a natural DIF effect size measure (i.e.,
the person covariate effect on the item slope 1j and item intercept 0j).
However, more thorough and careful study on the meaningful cut-off values are
needed to facilitate decision making. Alternatively, item level or scale level
effects can also be helpful and can be extended to the MNLFA modeling
framework. For example, average unsigned difference <see¿woods2011dif which
calculates weighted differences in the expected response functions between the
focal and reference groups can be extended to the MNLFA model for selected
values or levels of person covariate to evaluate the magnitude of the DIF
effect at the item level. In addition to the item level DIF impact, the
differential test function (DTF) index Roju . (1995) or the expected total
test score difference due to DTF Stark . (2004) can be computed to inform the
overall DIF effect on the scale level.

In sum, the proposed decorrelated score test and its one-step debiased
estimator based on the regularized moderated nonlinear factor analysis model
offers a promising solution to the practical obstacles encountered by
conventional IRT methods. Valid inference at the DIF parameter level opens
doors for more complicated substantive research questions such as
investigating the complex nature of DIF due to interconnection of background
characteristics.

## Appendix A

###  A.1 Penalized Expectation-Maximization Algorithm

The Expectation-Maximization (EM) algorithm Bock Aitkin (1982) is used to find
the local minimizer of $p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})$ with
respect to . The basic idea of the penalized EM algorithm is to repeatedly
approximate the upper bound of the penalized negative sample log-likelihood
$p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})\leq
q_{n}(\bfxi;{}^{(r)}),r=0,1,2,\cdots$ defined in Equation A1 at each iteration
(E-step) and obtain the optimizer (M-step). As a result, the alternation
between the E-step and the M-step produces a sequence of parameter updates
(r). The final parameter estimates $\hat{\bfxi}$ is referred to as the
penalized ML estimator. A critical property of the EM algorithm is that the
parameter estimate updates the penalized negative sample log-likelihood in a
non-increasing fashion (i.e.,
$p_{n}({}^{(r+1)};\mathbf{Y},\lambda|\mathbf{X})-p_{n}({}^{(r)};\mathbf{Y},\lambda|\mathbf{X})\leq
q_{n}({}^{(r+1)},{}^{(r)})-q_{n}({}^{(r)},{}^{(r)})\leq 0$ where
$q_{n}(\bfxi,{}^{r})$ can be viewed as the upper bound of
$p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})$). Specifically, at iteration $r$,
$q_{n}(\bfxi;{}^{(r)})$ is defined as

$\displaystyle q_{n}(\bfxi;{}^{(r)})$
$\displaystyle=\mathbb{E}_{(\theta_{i}|{}^{(r)},\mathbf{Y},\mathbf{X})}\left[p_{n}(\bfxi;\mathbf{Y};\theta_{i}|\mathbf{X})\right]$
(A1)
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{(\theta_{i}|{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})}\left[\log
f(\bfxi;\mathbf{y}_{i};\theta_{i}|\mathbf{x}_{i})\right]+\sum_{j=1}^{J}p_{\lambda}({}_{j})$
(A2) $\displaystyle=-\frac{1}{n}\sum_{i}^{n}\int
f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})\biggl{(}\log\left(f(\mathbf{y}_{i}|\theta_{i},\bfxi,\mathbf{x}_{i})f(\theta_{i}|\bfxi,\mathbf{x}_{i})\right)\biggl{)}d\theta_{i}$
$\displaystyle+\sum_{j=1}^{J}p_{\lambda}({}_{j})$ (A3)

Dropping $\sum_{j=1}^{J}p_{\lambda}({}_{j})$ for simplicity, the first term of
the right hand side (RHS) equals

First term of the RHS $\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\int
f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})\biggl{(}\sum_{j=1}^{J}\log
f_{j}(\mathbf{y}_{ij}|\theta_{i},{}_{j},\mathbf{x}_{i})+\log
f(\theta_{i}|\bfgamma,\bfdelta,\mathbf{x}_{i})\biggl{)}d\theta_{i}$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\int
f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})\biggl{(}\sum_{j=1}^{J}\log
f_{j}(\mathbf{y}_{ij}|\theta_{i},{}_{j},\mathbf{x}_{i})\biggl{)}d\theta_{i}$
$\displaystyle\quad-\frac{1}{n}\sum_{i=1}^{n}\int
f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})\biggl{(}\log
f(\theta_{i}|\bfgamma,\bfdelta,\mathbf{x}_{i})\biggl{)}d\theta_{i}$
$\displaystyle=-\frac{1}{n}\sum_{j=1}^{J}\sum_{i=1}^{n}\int
f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})\biggl{(}\log
f_{j}(\mathbf{y}_{ij}|\theta_{i},{}_{j},\mathbf{x}_{i})\biggl{)}d\theta_{i}$
$\displaystyle\quad-\frac{1}{n}\sum_{i=1}^{n}\int
f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})\biggl{(}\log
f(\theta_{i}|\bfgamma,\bfdelta,\mathbf{x}_{i})\biggl{)}d\theta_{i}.$

Further, denote $\theta_{iq},q=1,\cdots,Q$ and $w_{iq}$ as quadrature nodes
and weights for person $i$, respectively. The intractable integral in the
$q_{n}(\bfxi;{}^{(r)})$ can be approximated by summations on this quadrature
grid as

$\displaystyle q_{n}(\bfxi;{}^{(r)})$
$\displaystyle\approx-\frac{1}{n}\sum_{j=1}^{J}\sum_{i=1}^{n}\sum_{q=1}^{Q}e_{iq}^{(r)}\biggl{(}\log
f_{j}(\mathbf{y}_{ij}|\theta_{iq},{}_{j},\mathbf{x}_{i})\biggl{)}+\sum_{j}^{J}p_{\lambda}({}_{j})$
(A4)
$\displaystyle\quad-\frac{1}{n}\sum_{i=1}^{n}\sum_{q=1}^{Q}e_{iq}^{(r)}\biggl{(}\log
f(\theta_{iq}|\bfgamma,\bfdelta,\mathbf{x}_{i})\biggl{)},$ (A5)

where
$e_{iq}^{(r)}=f(\theta_{i}|{}^{(r)},{}^{(r)},\mathbf{y}_{i},\mathbf{x}_{i})$
is the posterior probability of $\theta_{i}$ at iteration $r$, which is
approximated by

$\displaystyle
e_{iq}^{(r)}=\frac{\prod_{j=1}^{J}f_{j}(y_{ij}|\theta_{iq},{}_{j}^{(r)},\mathbf{x}_{i})w_{iq}}{\sum_{q^{\prime}=1}^{Q}\prod_{j=1}^{J}f_{j}(y_{ij}|\theta_{iq^{\prime}},{}_{j}^{(r)},\mathbf{x}_{i})w_{iq^{\prime}}}.$
(A6)

Finally, the Bock-Aitkin EM algorithm with the $L_{1}$ penalty can be achieve
using the following steps until convergence is reached:

  * E-step. Compute the posterior weights $e_{iq}^{(r)}$, $i=1,\dots,n$, $q=1,\dots,Q$;

  * M-step.

    * For latent population parameters, compute

$\displaystyle({}^{(r+1)},{}^{(r+1)})^{\prime}=\arg\min_{\bfgamma,\bfdelta}\Biggl{\\{}-\frac{1}{n}\sum_{i=1}^{n}\sum_{q=1}^{Q}e_{iq}^{(r)}\biggl{(}\log
f(\theta_{iq}|\bfgamma,\bfdelta,\mathbf{x}_{i})\biggl{)}\Biggr{\\}}$ (A7)

    * For parameters for item $j$, compute

$\displaystyle{}_{j}^{(r+1)}=\arg\min_{{}_{j}}\Biggl{\\{}-\frac{1}{n}\sum_{i=1}^{n}\sum_{q=1}^{Q}e_{iq}^{(r)}\biggl{(}\log
f_{j}(y_{ij}|\theta_{iq},{}_{j},\mathbf{x}_{i})\biggl{)}+p_{\lambda}({}_{j})\Biggr{\\}}.$
(A8)

###  B.2 M-step Optimization

As was previously shown, at each iteration $r+1$, the M-step optimizes
${}^{(r+1)}=\arg\min p_{n}(\bfxi;\mathbf{Y},\lambda|\mathbf{X})$, which can be
split into two optimization problems shown in Equations A7 and A8. These two
equations update the population parameters (i.e., and ) and item parameters,
respectively. Equation A7 can be obtained by a Newton-type optimizer. The
Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is used in the current
study. In contrast, Equation A8 for each item $j$ needs to be handled
separately due to the non-differentiability of the $L_{1}$ penalty. Finding
${}_{j}^{(r+1)}$, $j=1,\dots,J$ in the M-step amounts to a conditional density
estimation problem in a sample size $nQ$ with weights $e_{iq}^{(r)}$.
Specifically, maximizing the unpenalized conditional density can be achieved
by solving the reweighted least square problem. With the $L_{1}$ penalty, the
coordinate descent algorithm Friedman . (2010) can be used to solve the
penalized weighted least square problem with pseudo data of a sample size
$nQ$.

## Appendix B

As mentioned before, the general theory of the decorrelated score test and its
one-step estimator are established under some assumptions Ning Liu (2017).
These assumptions need to be verified to guarantee the asymptotic properties
of the Dscore test and the asymptotic unbiased estimator as shown in Equations
14 and 16. The following four assumptions are briefly summarized and the
verification will be provided in the final dissertation.

######  Assumption 1 (Consistency conditions for initial parameter estimates)

For some sequence $a_{1}(n)$ and $a_{2}(n)$ converge to $0$ as
$n\rightarrow\infty$, it holds

$\displaystyle\lim_{n\to\infty}P_{{}^{*}}(\|\hat{\bfxi}-{}^{*}\|_{1}$
$\displaystyle\lesssim a_{1}(n))=1$ (A9)
$\displaystyle\lim_{n\to\infty}P_{{}^{*}}(\|\hat{{\mathbf{W}}}-{\mathbf{W}}^{*}\|_{1}$
$\displaystyle\lesssim a_{2}(n))=1,$ (A10)

where $\|\|_{1}$ stands for the $L_{1}$ operator norm of a matrix (e.g.,
$\|\mathbf{A}\|_{1}=\max_{1\leq j\leq m}\sum_{i=1}^{n}|a_{ij}|$) and
$\lesssim$ denotes that the left side is less than or equal to the right hand
side times some constant $C>0$. Although estimation consistency of LASSO has
been studied in linear and generalized linear models Knight Fu (2000); Ning
Liu (2017), it has never been studied in the latent variable modeling
framework. As the loss function or the negative sample log-likelihood of MNLFA
model is not strictly convex, the proof of parameter estimation consistency
can be non-trivial. Moreover, as mentioned by ning2017general, there can be
extra difficulty in bounding $\|\hat{{\mathbf{W}}}-{\mathbf{W}}^{*}\|$, as
$\hat{{\mathbf{W}}}$ depends on $\hat{\bfxi}$.

######  Assumption 2 (Concentration of the gradient and Hessian)

Let ${\bf V}^{*}=({\bf I}_{d_{0}\times d_{0}},-{\mathbf{W}^{*\top}})^{\top}$,
then assume

$\displaystyle\|\nabla\ell({}^{*})\|_{\infty}$
$\displaystyle=\text{O}_{p}(\sqrt{\log d/n})$ (A11)
$\displaystyle\|{\mathbf{V}}^{*\top}\nabla^{2}\ell({}^{*})-\operatorname{\mathbb{E}}_{{}^{*}}({\mathbf{V}}^{*\top}\nabla^{2}\ell({}^{*}))\|_{\infty}$
$\displaystyle=\text{O}_{p}(\sqrt{\log d/n}),$ (A12)

where $\|\|_{\infty}$ of a vector $\mathrm{vec}{A}$ for example indicates
$\|\mathrm{vec}{A}\|_{\infty}=\max_{1\leq i\leq d}|a_{i}|$ and of a matrix is
the maximum absolute row sum of the matrix (i.e.,
$\|{\mathbf{A}}\|_{\infty}=\max_{1\leq i\leq n}\sum_{j=1}^{m}|a_{ij}|$ ). In
the low dimensional setting, it might be sufficient to prove under some finite
moment assumptions on the gradient and Hessian matrix.

######  Assumption 3 (Local smoothness on the loss function)

Let $\hat{\bfxi}_{0}=({\bf 0},\hat{\bfeta}^{\top})^{\top}$,
$\hat{{\mathbf{V}}}=({\bf I}_{d_{0}\times
d_{0}},-\hat{\mathbf{W}}^{\top})^{\top}$, and ${{\mathbf{V}}}^{*}=({\bf
I}_{d_{0}\times d_{0}},-{\mathbf{W}}^{*\top})^{\top}$. For both
$\check{\bfxi}=\hat{\bfxi}_{0}$ and $\check{\bfxi}=\hat{\bfxi}$, it holds that

$\displaystyle\|{\mathbf{V}}^{*\top}\\{\nabla\ell(\check{\bfxi})-\nabla\ell({}^{*})-\nabla^{2}\ell({}^{*})(\check{\bfxi}-{}^{*})\\}\|_{\infty}$
$\displaystyle=o_{p}(n^{-1/2})$ (A13)
$\displaystyle\|(\hat{{\mathbf{V}}}-{\mathbf{V}}^{*})^{\top}(\nabla\ell(\check{\bfxi})-\nabla\ell({}^{*}))\|_{\infty}$
$\displaystyle=o_{p}(n^{-1/2}).$ (A14)

This assumption implicitly assumes that the $\ell(\bfxi)$ is second-order
differentiable. The verification of the assumption should be straightforward
if the loss function is a quadratic form of .

######  Assumption 4 (Convergence of the score function)

Let ${}^{*}=\lim_{n\to\infty}\hbox{Var}(n^{1/2}\nabla\ell({}^{*}))$. Then the
score function holds that

$\displaystyle\sqrt{n}\nabla\ell({}^{*})^{\top}{\mathbf{V}}^{*}{({\mathbf{V}}^{*\top}{}^{*}{\mathbf{V}}^{*})}^{-1}{\mathbf{V}}^{*\top}\nabla\ell({}^{*})\xrightarrow{\mathcal{D}}{\chi}^{2}_{d_{0}}.$
(A15)

This assumption can be established by verifying the Lindeberg’s condition,
which is a sufficient condition for a sequence of independent random
variables.

## References

  * Bauer (2017) bauer2017moreBauer, DJ.  2017\.  A more general model for testing measurement invariance and differential item functioning. A more general model for testing measurement invariance and differential item functioning.  Psychological Methods223507–526.
  * Bauer . (2020) bauer2020simplifyingBauer, DJ., Belzak, WC. Cole, VT.  2020\.  Simplifying the assessment of measurement invariance over multiple background variables: Using regularized moderated nonlinear factor analysis to detect differential item functioning Simplifying the assessment of measurement invariance over multiple background variables: Using regularized moderated nonlinear factor analysis to detect differential item functioning.  Structural Equation Modeling: A Multidisciplinary Journal27143–55.
  * Belzak (2021) belzak2021usingBelzak, W.  2021\.  Using Regularization to Evaluate Differential Item Functioning among Multiple Covariates: A Penalized Expectation-Maximization Algorithm via Coordinate Descent and Soft-Thresholding Using regularization to evaluate differential item functioning among multiple covariates: A penalized expectation-maximization algorithm via coordinate descent and soft-thresholding .  The University of North Carolina at Chapel Hill.
  * Belzak Bauer (2020) belzak2020improvingBelzak, W. Bauer, DJ.  2020\.  Improving the assessment of measurement invariance: Using regularization to select anchor items and identify differential item functioning. Improving the assessment of measurement invariance: Using regularization to select anchor items and identify differential item functioning.  Psychological Methods256673–690.
  * Bock Aitkin (1982) bock1982marginalBock, R. Aitkin, M.  1982\.  Marginal maximum likelihood estimation of item parameters Marginal maximum likelihood estimation of item parameters.  Psychometrika473369–369.
  * Chen . (2021) chen2021advantagesChen, SM., Bauer, DJ., Belzak, WM. Brandt, H.  2021\.  Advantages of Spike and Slab Priors for Detecting Differential Item Functioning Relative to Other Bayesian Regularizing Priors and Frequentist Lasso Advantages of spike and slab priors for detecting differential item functioning relative to other bayesian regularizing priors and frequentist lasso.  Structural Equation Modeling: A Multidisciplinary Journal1–18.
  * Embretson Reise (2000) EmbretsonEmbretson, SE. Reise, SP.  2000\.  Item response theory for psychologists Item response theory for psychologists.  Mahwah, New JerseyLawrence Erlbaum Associates, Inc.
  * Eysenck . (1985) eysenck1985revisedEysenck, SB., Eysenck, HJ. Barrett, P.  1985\.  A revised version of the psychoticism scale A revised version of the psychoticism scale.  Personality and Individual Differences6121–29.
  * Fan Li (2001) fan2001variableFan, J. Li, R.  2001\.  Variable selection via nonconcave penalized likelihood and its oracle properties Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American statistical Association964561348–1360.
  * Fang . (2017) fang2017testingFang, EX., Ning, Y. Liu, H.  2017\.  Testing and confidence intervals for high dimensional proportional hazards models Testing and confidence intervals for high dimensional proportional hazards models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)7951415–1437.
  * Friedman . (2010) friedman2010Friedman, J., Hastie, T. Tibshirani, R.  2010\.  Regularization Paths for Generalized Linear Models via Coordinate Descent Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software3311–22.  https://www.jstatsoft.org/v033/i01  10.18637/jss.v033.i01
  * Fu Knight (2000) fu2000asymptoticsFu, W. Knight, K.  2000\.  Asymptotics for lasso-type estimators Asymptotics for lasso-type estimators.  The Annals of Statistics2851356–1378.
  * Huang (2018) huang2018penalizedHuang, PH.  2018\.  A penalized likelihood method for multi-group structural equation modelling A penalized likelihood method for multi-group structural equation modelling.  British Journal of Mathematical and Statistical Psychology713499–522.
  * Jacobucci . (2016) jacobucci2016regularizedJacobucci, R., Grimm, KJ. McArdle, JJ.  2016\.  Regularized structural equation modeling Regularized structural equation modeling.  Structural Equation Modeling: A Multidisciplinary Journal234555–566.
  * Knight Fu (2000) knight2000asymptoticsKnight, K. Fu, W.  2000\.  Asymptotics for lasso-type estimators Asymptotics for lasso-type estimators.  The Annals of Statistics2851356–1378.
  * Lindstrøm Dahl (2020) lindstrom2020modelLindstrøm, JC. Dahl, FA.  2020\.  Model selection with lasso in multi-group structural equation models Model selection with lasso in multi-group structural equation models.  Structural Equation Modeling: A Multidisciplinary Journal27133–42.
  * M. Liu (2017) liu2017differentialLiu, M.  2017\.  Differential item functioning in large-scale mathematics assessments: Comparing the capabilities of the Rasch trees model to traditional approaches Differential item functioning in large-scale mathematics assessments: Comparing the capabilities of the rasch trees model to traditional approaches .  University of Toledo.
  * Y. Liu Yang (2018) liu2018bootstrapLiu, Y. Yang, JS.  2018\.  Bootstrap-calibrated interval estimates for latent variable scores in item response theory Bootstrap-calibrated interval estimates for latent variable scores in item response theory.  Psychometrika832333–354.
  * Magis . (2015) magis2015detectionMagis, D., Tuerlinckx, F. De Boeck, P.  2015\.  Detection of differential item functioning using the lasso approach Detection of differential item functioning using the lasso approach.  Journal of Educational and Behavioral Statistics402111–135.
  * Ning Liu (2017) ning2017generalNing, Y. Liu, H.  2017\.  A general theory of hypothesis tests and confidence regions for sparse high dimensional models A general theory of hypothesis tests and confidence regions for sparse high dimensional models.  The Annals of Statistics451158–195.
  * R Core Team (2020) RR Core Team.  2020\.  R: A Language and Environment for Statistical Computing R: A language and environment for statistical computing [].  Vienna, Austria.  https://www.R-project.org/
  * Roju . (1995) roju1995irtRoju, NS., Van der Linden, WJ. Fleer, PF.  1995\.  IRT-based internal measures of differential functioning of items and tests Irt-based internal measures of differential functioning of items and tests.  Applied Psychological Measurement194353–368.
  * Shea (2013) shea2013usingShea, CA.  2013\.  Using a mixture IRT model to understand English learner performance on large-scale assessments Using a mixture IRT model to understand english learner performance on large-scale assessments.  University of Massachusetts Amherst.
  * Stark . (2004) stark2004examiningStark, S., Chernyshenko, OS. Drasgow, F.  2004\.  Examining the effects of differential item (functioning and differential) test functioning on selection decisions: When are statistically significant effects practically important? Examining the effects of differential item (functioning and differential) test functioning on selection decisions: When are statistically significant effects practically important?  Journal of Applied Psychology893497.
  * Tutz Schauberger (2015) tutz2015penaltyTutz, G. Schauberger, G.  2015\.  A penalty approach to differential item functioning in Rasch models A penalty approach to differential item functioning in rasch models.  Psychometrika80121–43.
  * Vaart (1998) vaart_1998Vaart, AWvd.  1998\.  Semiparametric Models Semiparametric models.  Asymptotic Statistics Asymptotic statistics ( 358–432).  Cambridge University Press.  10.1017/CBO9780511802256.026
  * W. Wang . (2022) wangliuliu2022Wang, W., Liu, Y. Liu, H.  2022\.  Testing Differential Item Functioning Without Predefined Anchor Items Using Robust Regression Testing differential item functioning without predefined anchor items using robust regression.  Journal of Educational and Behavioral Statistics0010769986221109208.  https://doi.org/10.3102/10769986221109208  10.3102/10769986221109208
  * WC. Wang (2004) wang2004effectsWang, WC.  2004\.  Effects of anchor item methods on the detection of differential item functioning within the family of Rasch models Effects of anchor item methods on the detection of differential item functioning within the family of rasch models.  The Journal of Experimental Education723221–261.
  * WC. Wang Yeh (2003) wang2003effectsWang, WC. Yeh, YL.  2003\.  Effects of anchor item methods on differential item functioning detection with the likelihood ratio test Effects of anchor item methods on differential item functioning detection with the likelihood ratio test.  Applied Psychological Measurement276479–498.
  * Woods (2009) woods2009empiricaldifWoods, CM.  2009\.  Empirical selection of anchors for tests of differential item functioning Empirical selection of anchors for tests of differential item functioning.  Applied Psychological Measurement33142–57.
  * Woods (2011) woods2011difWoods, CM.  2011\.  DIF testing for ordinal items with poly-SIBTEST, the Mantel and GMH tests, and IRT-LR-DIF when the latent distribution is nonnormal for both groups Dif testing for ordinal items with poly-sibtest, the mantel and gmh tests, and irt-lr-dif when the latent distribution is nonnormal for both groups.  Applied Psychological Measurement352145–164.
  * Woods . (2013) woodscaiWangWoods, CM., Cai, L. Wang, M.  2013\.  The Langer-Improved Wald Test for DIF Testing With Multiple Groups: Evaluation and Comparison to Two-Group IRT The Langer-improved wald test for DIF testing with multiple groups: Evaluation and comparison to two-group IRT.  Educational and Psychological Measurement733532-547.  https://doi.org/10.1177/0013164412464875
  * Zhao Yu (2006) zhao2006modelZhao, P. Yu, B.  2006\.  On model selection consistency of Lasso On model selection consistency of lasso.  The Journal of Machine Learning Research72541–2563.
  * Zou (2006) zou2006adaptiveZou, H.  2006\.  The adaptive lasso and its oracle properties The adaptive lasso and its oracle properties.  Journal of the American statistical association1014761418–1429.

## Figures

Figure 1: Conditional probability of endorsing an item. $2,500$ item responses
to all 12 items were generated using the moderated nonlinear factor analysis
model with three covariates. Items 1 to 6 are DIF items. Figure 2: Type I
error results of incorrectly detecting a DIF item under the null condition
when there is no DIF items. Different methods are displayed in different
colors and shapes. Two dotted horizontal line shows the 95% normal-
approximation confidence band at the nominal level $0.05$ (horizontal dashed
line). The oracle solution only performed on item 1 to item 10 as the last two
items are treated as anchors. Figure 3: Power results of correctly detecting a
DIF item under the alternative condition when there is a mix of DIF and DIF-
free items. Different methods are displayed in different colors and shapes.
For each method, the empirical rejection rate under the nominal level $0.05$
are calculated for each DIF item under each condition. The column shows the
condition when the number of DIF items is 3 or 6 out of 12 items. Each row
represents a different sample size condition. The reference dashed line shows
the empirical rejection rate = 0.05 Figure 4: False detection rate results of
incorrectly detecting a DIF item under the alternative condition when there is
a mix of DIF and DIF-free items. Different methods are displayed in different
colors and shapes. For each method, the empirical rejection rate at the
nominal level $0.05$ are calculated for each DIF item under each condition.
Two dotted horizontal line shows the 95% normal approximation confidence band
for the nominal level $0.05$ (horizontal dashed lines). The column shows the
condition when the number of DIF items is 3 or 6 items out of 12 items. Each
row represents a sample size condition. Figure 5: Average bias of a-DIF
parameters. Different methods are displayed in different colors and shapes.
The column shows the condition when the number of DIF items is 0, 3 or 6 out
of 12 items. Each row represents a sample size condition. DIF items are shown
in the grey shaded area. Figure 6: Average bias of d-DIF parameters. Different
methods are displayed in different colors and shapes. The column shows the
condition when the number of DIF items is 0, 3 or 6 out of 12 items. Each row
represents a sample size condition. DIF items are shown in the grey shaded
area. Figure 7: Density plot of a non-zero DIF parameter by different method.
The vertical reference lines indicate the means of the parameter estimates of
each effect for each method. The example is plotted using the sample size
$2,500$ and first 3-item DIF condition. The non-zero DIF effect is the
continuous DIF effect of the first item. Figure 8: Average variance of a-DIF
parameters. Different methods are displayed in different colors and shapes.
The column shows conditions when the number of DIF items is 0, 3 or 6 out of
12 items. Each row represents a specific sample size condition. The grey
dashed reference line displays the mean variance across all items of the
oracle solution for each condition to be used as a benchmark. DIF items are
shown in the grey shaded area. Figure 9: Average variance of d-DIF parameters.
Different methods are displayed in different colors and shapes. The column
shows conditions when the number of DIF items is 0, 3 or 6 out of 12 items.
Each row represents a specific sample size condition. The grey dashed
reference line displays the mean variance across all items of the oracle
solution for each condition to be used as a benchmark. DIF items are shown in
the grey shaded area. Figure 10: Average standard error recovery of a-DIF
parameters. Different methods are displayed in different colors and shapes.
The column shows conditions when the number of DIF items is 0, 3 or 6 out of
12 items. Each row represents a specific sample size condition. The grey
dashed reference line displays ratio of 1 indicating perfect recovery. Figure
11: Average standard error recovery of d-DIF parameters. Different methods are
displayed in different colors and shapes. The column shows conditions when the
number of DIF items is 0, 3 or 6 out of 12 items. Each row represents a
specific sample size condition. The grey dashed reference line displays ratio
of 1 indicating perfect recovery.

## Tables

Table 1: Algorithm 1: Decorrelated Score Function Estimation Algorithm 1
Estimated the decorrelated score function

1:Negative sample log-likelihood $\ell(\bfpsi,\bfeta)$, penalty function
$p_{\lambda}()$, and tuning parameters $\lambda$ and $\lambda^{\prime}$.

2:Estimate $\hat{\bfxi}$ using penalized ML as in Equation 10 and partition
$\hat{\bfxi}$ into
$\hat{\bfxi}=(\hat{\bfpsi}^{\top},\hat{\bfeta}^{\top})^{\top}$

3:Estimate $\mathbf{W}$ column by column
$\displaystyle\hat{\mathbf{W}}_{*j}=\operatorname*{\arg\min}_{{\bf
w}_{j}}\frac{1}{2n}\sum_{i}^{n}\\{\nabla_{{}_{j}}\ell_{i}(\hat{\bfxi})-\mathbf{W}_{*j}^{\top}\nabla\ell_{i}(\hat{\bfxi})\\}^{2}+p_{\lambda^{\prime}}(\mathbf{W}_{*j})$
(A16)

4:Calculate the estimated descorrelated score function using
$\displaystyle\hat{\mathbf{s}}(\bfpsi,\hat{\bfeta})$
$\displaystyle=\nabla\ell(\bfpsi,\hat{\bfeta})-\hat{\mathbf{W}}^{\top}\nabla\ell(\bfpsi,\hat{\bfeta})$
(A17) Return $\hat{\mathbf{s}}(\bfpsi,\hat{\bfeta})$

Table 2: Model Parameters of the True Data Generating Model

|  |  | a-DIF | d-DIF
---|---|---|---|---
|  |  |  Age |  Gender |  Product |  Age |  Gender |  Product
Item |  $d_{j}$ |  $a_{j}$ |  $\beta_{11}$ |  $\beta_{12}$ |  $\beta_{13}$ |  $\beta_{01}$ |  $\beta_{02}$ |  $\beta_{03}$
1 |  0.00 |  2.00 |  0.20 |  0.50 |  0.20 |  0.20 |  -0.50 |  -0.20
2 |  1.20 |  1.20 |  -0.20 |  -0.50 |  0.00 |  -0.20 |  0.25 |  0.00
3 |  -0.20 |  2.00 |  -0.25 |  0.25 |  0.10 |  -0.15 |  -0.25 |  -0.15
$4^{*}$ |  1.50 |  1.50 |  0.20 |  -0.50 |  -0.20 |  0.20 |  0.50 |  0.20
$5^{*}$ |  1.20 |  1.20 |  -0.20 |  -0.50 |  0.00 |  -0.20 |  0.25 |  0.00
$6^{*}$ |  1.10 |  1.90 |  -0.25 |  0.25 |  -0.10 |  -0.15 |  -0.25 |  0.15
7 |  -1.80 |  2.40 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00
8 |  0.50 |  1.50 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00
9 |  0.60 |  1.40 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00
10 |  -2.00 |  1.80 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00
11 |  0.60 |  2.30 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00
12 |  1.60 |  1.80 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00 |  0.00
| $(-0.2,-0.2,-0.2)^{\top}$
| $(-0.1,0.3,0.1)^{\top}$

Note. * indicates that effect sizes of these items under the 25% DIF condition
are 0.