tabular

# The generalized Hausman test for detecting non-normality in the latent
variable distribution of the two-parameter IRT model

Lucia Guastadisegni<EMAIL_ADDRESS>University of Bologna
Silvia Cagnone  University of Bologna  Irini Moustaki  London School of
Economics and Political Science  Vassilis Vasdekis  Athens University of
Economics and Business

###### Abstract

This paper introduces the generalized Hausman test as a novel method for
detecting non-normality of the latent variable distribution of unidimensional
Item Response Theory (IRT) models for binary data. The test utilizes the
pairwise maximum likelihood estimator obtained for the parameters of the
classical two-parameter IRT model, which assumes normality of the latent
variable, and the quasi-maximum likelihood estimator obtained under a semi-
nonparametric framework, allowing for a more flexible distribution of the
latent variable. The performance of the generalized Hausman test is evaluated
through a simulation study and it is compared with the likelihood-ratio and
the $M_{2}$ test statistics. Additionally, various information criteria are
computed. The simulation results show that the generalized Hausman test
outperforms the other tests under most conditions. However, the results
obtained from the information criteria are somewhat contradictory under
certain conditions, suggesting a need for further investigation and
interpretation.

Keywords— Semi-non-parametric-IRT model, misspecification test, correlated
binary data

##  1 Introduction

The aim of the paper is to propose a generalized Hausman type test for
detecting deviation of the latent variable distribution from normality. The
paper focuses on mixture of normals and skewed normal distributions.

In many fields of research, there is a need to study theoretical constructs,
such as abilities, quality of life, or business confidence, which cannot be
directly observed and measured. To address this, latent variable models are
employed, which measure the constructs of interest also known as latent
variables or factors by analyzing the associations or correlations among
observed variables also known as manifest variables or items (Bartholomew et
al., , 2011). Latent variable models can accommodate both continuous and
discrete variables, including binary and polytomous outcomes. In the case of
binary or polytomous outcomes, Item Response Theory (IRT) models are commonly
used. These models are a type of latent variable model where the observed
outcomes are binary or polytomous, and the latent variables are assumed to be
continuous (van der Linden and Hambleton, , 2013). The use of IRT models is
particularly prevalent in social, psychological, and educational research,
where they are employed to measure various constructs such as attitudes or
abilities.

One of the standard assumptions of IRT models is that the latent variable(s)
follow a normal distribution. However, assuming normality of the latent
variable(s) when the true distribution has a different shape can lead to
biased parameter estimates, especially with binary outcomes (Ma and Genton, ,
2010). Furthermore, assuming an incorrect distribution of the latent variable
can lead to erroneous conclusions when conducting hypothesis testing
(Guastadisegni et al., , 2022). In the literature of the generalized linear
latent variable models (GLLVM) (Bartholomew et al., , 2011, Skrondal and Rabe-
Hesketh, , 2004) and IRT models, several methods that assume a different form
for the distribution of the latent variable have been proposed. Montanari and
Viroli,  (2010) introduce a skew-normal latent variable in the factor model,
while Cagnone and Viroli,  (2012) present a latent trait model where the
factors are distributed as a finite mixture of multivariate gaussians. Within
the GLLVM framework, Ma and Genton,  (2010) propose a semi-parametric method,
consistent for various types of manifest variables under different
distributions of the latent variables and Irincheeva et al.,  (2012) consider
the semi-non-parametric (SNP) approach, introduced by Gallant and Nychka,
(1987). This approach allows for more flexible smooth densities of the latent
variables. The SNP method has been used also in unidimensional IRT model by
Woods and Lin,  (2009) and in multidimensional IRT model by Monroe,  (2014).
For binary responses, Knott and Tzamourani,  (2007) estimate the distribution
of the latent variable from the data, using the empirical histogram method
combined with the bootstrap. Woods,  (2006) proposes the so-called Ramsey
curve IRT model, where the latent variables are splines based densities, that
are linear combination of polynomial functions joint together at knots. This
method implies a modification of the standard E-M algorithm.

In the majority of cases, information criteria, such as the Akaike information
criterion (AIC) (Akaike, , 1974) or the Bayesian information criterion (BIC)
(Schwarz, , 1978), are used to choose between a model where the latent
variables are normally distributed and a model where the latent variables have
a more complex shape (Woods and Lin, , 2009, Irincheeva et al., , 2012,
Monroe, , 2014). With continuous manifest variables, Ma and Genton,  (2010)
perform the Kolmogorov–Smirnov test on the distribution of the continuous
responses to evaluate the normality of the latent variable. However, when the
responses are categorical, detecting non-normality of the latent variables
through a statistical test remains an open issue.

Hausman,  (1978) proposes a specification test to detect failure of the
orthogonality assumption in regression analysis. Due to its simplicity, the
Hausman test can be applied in various contexts to detect different types of
model misspecification. This test is based on the comparison between two
different estimators that are consistent when the model is correctly specified
and one of them is also efficient. In presence of model misspecification, only
the inefficient estimator is consistent. The efficiency assumption simplifies
the computation of the covariance matrix of the difference between the two
estimators. However, this matrix can fail to be positive definite under model
misspecification or in presence of small sample sizes. Moreover, in some cases
none of the two estimators considered are fully efficient.

The generalized version of the Hausman (GH) test, proposed by White,  (1982),
is a more flexible and robust alternative to the original Hausman test.
Indeed, the generalized version allows for both estimators to be inefficient
and to result from two different models. Moreover, the covariance matrix of
the difference of the two estimators is robust and always positive definite.

In the IRT context, as far as we know, the Hausman test has been used only by
Ranger and Much,  (2020) to detect misspecification of the item characteristic
functions and local dependencies among items. They highlight that this test
has good performance in terms of Type I error rates for large sample sizes and
power under most conditions. In generalized linear mixed models (GLMM) for
clustered data, a robust version of the Hausman test, similar to the one by
White,  (1982), has been proposed by Bartolucci et al.,  (2017) when a
discrete distribution for the random effects is assumed. The test can be also
used to detect the possible correlation between random effects and cluster-
specific covariates. With respect to the information criteria, they found that
the robust Hausman test prefers more parsimonious models and it can detect the
presence of endogeneity.

The objective of this work is to extend the GH test to detect non-normality of
the latent variable distribution in unidimensional IRT models for binary data.
To build the test, the estimators resulting from two different models are
considered. The first model is a two-parameter logistic (2PL) unidimensional
IRT model that assumes normality of the latent variable, while the second
model is the unidimensional SNP-IRT model that assumes a more flexible
distribution for the latent variable. The 2PL model is estimated using a
pairwise maximum likelihood (PL) method (Katsikatsou et al., , 2012,Vasdekis
et al., , 2012), which is a composite likelihood method that uses information
from bivariate-order margins (Lindsay, , 1988, Varin, , 2008). The SNP-IRT
model is estimated using a quasi-maximum likelihood (quasi-ML) method (White,
, 1982). The choice of these estimators and models is motivated by the
following reasons. First, both methods are consistent when the latent variable
is normally distributed. Moreover, the quasi-ML method for the SNP-IRT model
is consistent also under different distribution assumptions of the latent
variable (Gallant and Tauchen, , 1989, Irincheeva et al., , 2012). These
conditions on the consistency of the parameter estimators are required to
correctly apply the GH test (White, , 1982). Second, the PL estimator and the
quasi-ML estimator yield different values for the parameter variance. This
implies that, also under normality of the latent variable distribution, the
covariance matrix of the difference of the two estimators involved in the
computation of the GH test is different from zero. The non-zero covariance
matrix avoids numerical instability in the computation of the test.

The theoretical aspects of the models, estimators and matrices involved in the
computation of the GH test are described. Moreover, we carry out an extensive
simulation study to evaluate the performance of the GH test both in terms of
Type I error rates and empirical power. For the latter, we consider both
mixture of normals and skew-normal distributions for the latent variable, with
varying degrees of departure from the normal distribution. Additionally, we
evaluate the asymptotic behavior of the test, in terms of both Type I error
rates and power, using a very large sample size. The performance of the GH
test is compared with the $M_{2}$ test, a statistic based on a quadratic form
in marginal residuals of order one and two, not affected by the problem of
sparse data (Maydeu-Olivares and Joe, , 2005) and the likelihood-ratio (LR)
test for nested models. In addition, some information criteria are computed.
An application to real data is also presented.

The article is organized as follows. In Section 2, we describe the 2PL model
and SNP-IRT models for binary data and the PL and quasi-ML estimators,
respectively. In Section 3, we present the GH test to detect non-normality of
the latent variable distribution. In Section 4 we review the $M_{2}$ test and
the LR test and in Section 5 the information criteria. In Section 6, we
present a Monte Carlo simulation study and in Section 7 the results from a
real data analysis. Finally, in Section 8, some concluding remarks are
presented and discussed.

##  2 The 2PL and SNP-IRT model for binary data

Let us denote by $y_{1},...,y_{p}$ a set of observed binary variables/items,
by $n$ the number of individuals and by $z$ the latent variable with density
function $h(z)$.

According to the 2PL model, the response category probability for the $i$-th
individual to the $j$-th item is modelled using a logistic model (measurement
model)

$\begin{split}&P(y_{ij}=1|z_{i})={\pi_{ij}}(z_{i})=\frac{\exp{(\alpha_{0j}+\alpha_{1j}z_{i})}}{1+\exp{(\alpha_{0j}+\alpha_{1j}z_{i})}},\hskip
22.76228ptj=1,\ldots,p\end{split}$ (1)

where $\alpha_{0j}$ is the item intercept and $\alpha_{1j}$ the item slope
(factor loadings) and $h(z)=\phi(z)$, where $\phi(z)$ is the density of a
standard normal. An extension of the IRT model is given by the SNP-IRT model
that assumes the same response probability as (1) but a semi-non-parametric
parametrization of the latent variable as follows

$\begin{split}h(z_{i})=P_{L}^{2}(z_{i})\phi(z_{i})\hskip
8.5359pt\text{and}\hskip 8.5359ptP_{L}(z_{i})=\sum_{0\leq l\leq
L}a_{i}z_{i}^{l},\end{split}$ (2)

where $a_{0},...,a_{L}$ are the real coefficients of the polynomial
$P_{L}(z_{i})$ and $L$ is the polynomial degree.

In order for $h(z)$ to be a density, the coefficients $a_{0},...,a_{L}$ of
$P_{L}(z)$ should be chosen such that $\int h(z)dz=1$. For this purpose,
Gallant and Tauchen,  (1989) use a proportionality constant $1/\int
P_{L}(z)^{2}\phi(z)dz$ and fix the constant term of the polynomial equal to 1.
Alternatively, Woods and Lin,  (2009) and Irincheeva et al.,  (2012) use the
parametrization proposed by Zhang and Davidian,  (2001), that imposes

$1=\int_{R}P_{L}^{2}(z)\phi(z)dz=E\\{P_{L}^{2}(w)\\}=\textbf{a}^{\prime}E(\boldsymbol{\tilde{w}}\boldsymbol{\tilde{w}}^{\prime})\textbf{a}=\textbf{a}^{\prime}A\textbf{a}$
(3)

with $w\sim N(0,1)$, $P_{L}(w)=\textbf{a}^{\prime}\boldsymbol{\tilde{w}}$, and
$\boldsymbol{\tilde{w}}^{\prime}=(1,w,w^{2},...,w^{L})$. The matrix $A$ is
positive definite by definition and $A=B^{\prime}B$, where $B$ is a positive
definite matrix.

If $\textbf{c}=B\textbf{a}$, equation (3) becomes
$\textbf{c}^{\prime}\textbf{c}=1$ and
$\textbf{c}^{\prime}=(c_{1},...,c_{L+1})$. The elements of c can be
represented using a polar coordinate transformation as
$c_{1}=\sin{\varphi_{1}}$,
$c_{2}=\cos{\varphi_{1}}\sin{\varphi_{2}},\ldots,c_{L}=\cos{\varphi_{1}}\ldots
cos{\varphi_{L-1}}\sin{\varphi_{L}},c_{L+1}=\cos{\varphi_{1}}\cos{\varphi_{2}}\ldots\cos{\varphi_{L-1}}\cos{\varphi_{L}}$,
with angles $-\pi/2\leq\varphi_{l}\leq\pi/2$, $l=1,...,L$. The density of the
latent variable in (2) can be expressed as

$h(z|\boldsymbol{\varphi},L)=(\textbf{a}^{\prime}\tilde{\textbf{z}})^{2}\phi(z),$
(4)

where a can be obtained from c as $\textbf{a}=B^{-1}\textbf{c}$,
$\tilde{\textbf{z}}^{\prime}=(1,z,z^{2},...,z^{L})$ and
$\boldsymbol{\varphi}^{\prime}=(\varphi_{1},...,\varphi_{L})$. When $L=0$ the
distribution of the latent variable reduces to the normal one. When $L=1$,
$P_{L}(z)=a_{0}+a_{1}z$, $a_{0}=\sin{\varphi_{1}}$, $a_{1}=\cos{\varphi_{1}}$.
The SNP parametrization with $L=1$ includes unimodal and bimodal
distributions.

Figure 1 illustrate the SNP densities of $z$ when $L=1$, for different values
of the $\varphi_{1}$ parameter.

Figure 1: SNP densities of $z$ when $L=1$, for different values of the
$\varphi_{1}$ parameter

When $\varphi_{1}$ is negative and close to -1, the distribution is slightly
right-skewed, whereas when it is close to 1, it is slightly left-skewed. When
the values of $\varphi_{1}$ are between -1 and 1, the distributions are
bimodal. Even if not reported in the graph, when
$\varphi_{1}=\pm\frac{\pi}{2}$, the SNP distribution reduces to the normal
one.

When $L=2$, $P_{L}(z)=a_{0}+a_{1}z+a_{2}z^{2}$,
$a_{0}=\sin{\varphi_{1}}-\frac{1}{\sqrt{2}}\cos{\varphi_{1}}\cos{\varphi_{2}}$,
$a_{1}=\cos{\varphi_{1}}\sin{\varphi_{2}}$ and
$a_{2}=\frac{1}{\sqrt{2}}\cos{\varphi_{1}}\cos{\varphi_{2}}$. The SNP
parametrization with $L=2$ is more flexible than $L=1$ and encompasses
unimodal, multimodal (including up to three modes), and skewed distributions.

Figures 2 displays the SNP densities of $z$ when $L=2$, for different values
of the $\varphi_{1}$ and $\varphi_{2}$ parameters.

Figure 2: SNP densities of $z$ when $L=2$, for different values of the
$\varphi_{1}$ and $\varphi_{2}$ parameters

When the value of $\varphi_{1}$ is negative and the value of $\varphi_{2}$ is
close to -1 or 1, the distributions are bimodal. When both parameters are
close to 0, the distributions are trimodal. When the value of $\varphi_{1}$ is
positive and greater than 0.5, and the value of $\varphi_{2}$ is close to -1
or 1, the distributions are highly right-skewed and left-skewed, respectively.
Even though it is not reported in Figure 2, when
$\varphi_{2}=\pm\frac{\pi}{2}$, the SNP densities reduce to those observed
when $L=1$. Moreover, when both parameters are set to $\pm\frac{\pi}{2}$, the
SNP densities return to the normal case.

We indicate with $SNP_{2}$ the 2PL model for $L=2$, with $SNP_{1}$ the 2PL
model for $L=1$ and with $SNP_{0}$ the 2PL model for $L=0$.

###  2.1 Pairwise estimator for the $\boldsymbol{SNP_{0}}$ model

To implement the GH test, the parameters of the $SNP_{0}$ model are estimated
with the pairwise method. The pairwise log-likelihood, based on the bivariate
marginal densities $f({y}_{ij},y_{ik},\boldsymbol{\theta})$, $j,k=1,..p$ and
$k>j$, is

$\begin{split}&pl_{SNP_{0}}(\mathbf{y},\boldsymbol{\theta})=\sum_{i=1}^{n}\sum_{j=1}^{p}\sum_{k>j}\ln
f({y}_{ij},y_{ik},\boldsymbol{\theta})=\\\
&=\sum_{i=1}^{n}\sum_{j=1}^{p}\sum_{k>j}\ln\int\bigg{[}\pi_{ij}(z_{i})^{y_{ij}}(1-\pi_{ij}(z_{i}))^{1-y_{ij}}\bigg{]}\bigg{[}\pi_{ik}(z_{i})^{y_{ik}}(1-\pi_{ik}(z_{i}))^{1-y_{ik}}\bigg{]}\phi(z_{i})dz_{i}.\end{split}$
(5)

$pl_{SNP_{0}}(\mathbf{y},\boldsymbol{\theta})$ is maximized with respect to
$\boldsymbol{\theta}$, where $\boldsymbol{\theta}$ includes the item
intercepts and slopes. Under correct model specification, the maximum PL
estimator $\tilde{\boldsymbol{\theta}}$ converges in probability to the true
parameter vector
$\boldsymbol{\theta}_{0}^{\prime}=(\boldsymbol{\alpha}_{00}^{\prime},\boldsymbol{\alpha}_{01}^{\prime})$
and

$\tilde{\boldsymbol{\theta}}\xrightarrow[]{p}N(\boldsymbol{\theta}_{0},A^{-1}(\boldsymbol{\theta}_{0})B(\boldsymbol{\theta}_{0})A^{-1}(\boldsymbol{\theta}_{0})),$
(6)

where
$A(\boldsymbol{\theta})=E_{y}\bigg{[}-\frac{\partial^{2}pl_{SNP_{0}}(\mathbf{y},\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\bigg{]}$,
$B=var\bigg{[}\frac{\partial
pl_{SNP_{0}}(\mathbf{y},\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\bigg{]}$
and ${A}(\boldsymbol{\theta})\neq{B}(\boldsymbol{\theta})$ (Lindsay, , 1988,
Varin, , 2008). These matrices can be estimated by their observed versions
evaluated at $\tilde{\boldsymbol{\theta}}$ as

$\hat{A}({\tilde{\boldsymbol{\theta}}})=\left.\sum_{i=1}^{n}\frac{\partial^{2}pl_{SNP_{0}}(\mathbf{y}_{i},\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\right|_{\boldsymbol{\theta}=\tilde{\boldsymbol{\theta}}}$
(7)

and

$\hat{B}(\tilde{\boldsymbol{\theta}})=\left.\sum_{i=1}^{n}\frac{\partial
pl_{SNP_{0}}(\boldsymbol{y}_{i},\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\frac{\partial
pl_{SNP_{0}}(\boldsymbol{y}_{i},\boldsymbol{\theta})}{\partial\boldsymbol{\theta}^{\prime}}\right|_{\boldsymbol{\theta}=\tilde{\boldsymbol{\theta}}}.$
(8)

###  2.2 Quasi-ML estimator for the $\boldsymbol{SNP_{L}}$ model

The parameter vector
$\boldsymbol{\theta}^{{(1)}^{\prime}}=(\boldsymbol{\alpha}_{0}^{\prime},\boldsymbol{\alpha}_{1}^{\prime},\boldsymbol{\varphi}^{\prime})=(\boldsymbol{\theta}^{\prime},\boldsymbol{\varphi}^{\prime})$
of the $SNP_{L}$ model, where $L>0$, is estimated using the quasi-ML method.
The quasi-log-likelihood of the data is

$\begin{split}l_{SNP_{L}}(\mathbf{y},\boldsymbol{\theta}^{{(1)}})&=\sum_{i=1}^{n}\ln
f(\mathbf{y}_{i},\boldsymbol{\theta}^{{(1)}})=\\\
&=\sum_{i=1}^{n}\ln\int\prod_{j=1}^{p}\pi_{ij}(z_{i})^{y_{ij}}(1-\pi_{ij}(z_{i}))^{1-y_{ij}}P_{L}^{2}(z_{i})\exp{\bigg{(}-\frac{1}{2}z_{i}^{\prime}z_{i}\bigg{)}}dz_{i}.\end{split}$
(9)

The integral in $l_{SNP_{L}}(\mathbf{y},\boldsymbol{\theta}^{{(1)}})$ is
approximated by using the Gauss-Hermite quadrature, as in Woods and Lin,
(2009). The degree of the polynomial $L$ is fixed and is not estimated by
maximum likelihood. The quasi-log-likelihood function is maximized with
respect to the unknown vector of parameter $\boldsymbol{\theta}^{{(1)}}$ as
follows

$(\accentset{\mbox{\large.}}{\boldsymbol{{\alpha}}}_{0}^{\prime},\accentset{\mbox{\large.}}{\boldsymbol{{\alpha}}}_{1}^{\prime},\hat{\boldsymbol{\varphi}}^{\prime})=argmax_{\boldsymbol{\theta}}l_{SNP_{L}}(\mathbf{y},\boldsymbol{\theta}^{{(1)}}).$
(10)

For identifiability reasons, the item intercepts and slopes should be rescaled
as described in the next section.

####  2.2.1 The identifiability problem and the final estimators

The standard argument for obtaining the final form of the estimators
$\hat{\boldsymbol{\alpha}}_{0}$ and $\hat{\boldsymbol{\alpha}}_{1}$ relies on
the concept of indeterminacy of the scale of the latent variable (Lord, ,
1980, van der Linden and Barrett, , 2016). Indeed, in latent variable models,
it is possible to change the scale of the latent variable $z$ changing the
parameterization accordingly and obtain the same joint density function of the
observed variables.

After maximizing the log-likelihood in (9), the latent variable $z$ has a
density $h(z|\hat{\boldsymbol{\varphi}},L)$ with an estimated mean
$\tilde{E}(Z)$ and variance $\tilde{V}(Z)$. The estimators
$\accentset{\mbox{\large.}}{\boldsymbol{\alpha}}_{0}$ and
$\accentset{\mbox{\large.}}{\boldsymbol{\alpha}}_{1}$ are on the same scale
defined by $h(z|\hat{\boldsymbol{\varphi}},L)$. To compare the 2PL model,
which assumes $z\sim N(0,1)$, with the SNP-IRT model, we rescale the mean and
the variance of the latent variable of the SNP-IRT model to 0 and 1,
respectively, along with adjusting the parametrization accordingly. After the
optimization process, we can express

$logit(\pi_{j}(z))=\accentset{\mbox{\large.}}{\alpha}_{0j}+\accentset{\mbox{\large.}}{\alpha}_{1j}z\hskip
28.45274ptj=1,...,p$ (11)

and

$z=\sqrt{\tilde{V}(Z)}z_{1}+{\tilde{E}}({Z})$ (12)

where $\tilde{E}({Z})$ and $\tilde{V}({Z})$ are found given
$\hat{\boldsymbol{\varphi}}$ and the SNP density of $z$ and $z_{1}$ has the
same distribution of $z$, but with mean 0 and variance 1.

If we substitute (12) in (11) we get

$logit(\pi_{j}(z))=\accentset{\mbox{\large.}}{\alpha}_{0j}+\accentset{\mbox{\large.}}{\alpha}_{1j}\sqrt{\tilde{V}(Z)}z_{1}+\accentset{\mbox{\large.}}{\alpha}_{1j}\tilde{E}(Z)\hskip
28.45274ptj=1,...,p$ (13)

From equation (13) we get the form of the final estimators, which corresponds
to a latent variable with mean 0 and variance 1 (Irincheeva et al., , 2012):

$\hat{\alpha}_{0j}=\accentset{\mbox{\large.}}{\alpha}_{0j}+\accentset{\mbox{\large.}}{\alpha}_{1j}\tilde{E}(Z)\hskip
28.45274ptj=1,...,p$ (14)
$\hat{\alpha}_{1j}=\accentset{\mbox{\large.}}{\alpha}_{1j}\sqrt{\tilde{V}(Z)}\hskip
28.45274ptj=1,...,p,$ (15)

$\tilde{E}({Z})$ and $\tilde{V}({Z})$ can be computed analytically, given the
values of $\hat{\boldsymbol{\varphi}}$, as shown in the Appendix A. Finally,
$\hat{\boldsymbol{\theta}}^{{(1)}^{\prime}}=(\hat{\boldsymbol{\alpha}}^{{}^{\prime}}_{0},\hat{\boldsymbol{\alpha}}^{{}^{\prime}}_{1},\hat{\boldsymbol{\varphi}}^{{}^{\prime}})$.

Under normal, multi-modal and asymmetric distributions of the latent variables
and if the regularity conditions A2-A6 of White,  (1982) are satisfied,

$\hat{\boldsymbol{\theta}}^{{(1)}}\xrightarrow[]{p}N({\boldsymbol{\theta}^{{(1)}}_{0}},A^{-1}({\boldsymbol{\theta}^{{(1)}}_{0}})B({\boldsymbol{\theta}^{{(1)}}_{0}})A^{-1}({\boldsymbol{\theta}^{{(1)}}_{0}})),$
(16)

where
$\boldsymbol{\theta}^{{(1)}^{\prime}}_{0}=(\boldsymbol{\alpha}_{00}^{\prime},\boldsymbol{\alpha}_{01}^{\prime},\boldsymbol{\varphi_{*}}^{\prime})=(\boldsymbol{\theta}_{0}^{\prime},\boldsymbol{\varphi_{*}}^{\prime})$.
$\boldsymbol{\varphi_{*}}$ is the value of $\boldsymbol{\varphi}$ that
minimizes the Kullback-Leibler information criterion (White, , 1982, Gallant
and Tauchen, , 1989, Irincheeva et al., , 2012). If the true latent variable
follows an $SNP_{L}$ density, the vector $\boldsymbol{\varphi_{*}}$ coincides
with the true parameter value $\boldsymbol{\varphi}_{0}$ and the quasi-ML
method reduces to the classic full-ML method. $A(\boldsymbol{\theta})$ and
$B(\boldsymbol{\theta})$ are the expected Hessian and cross-product matrices,
respectively. Their observed versions can be computed with the Delta method
(Cramér, , 1946) and are defined similarly to (7) and (8), where
$pl_{SNP_{0}}(\mathbf{y}_{i},\boldsymbol{\theta})$ is replaced by
$l_{SNP_{L}}(\mathbf{y}_{i},\boldsymbol{\theta}^{(1)})$.

##  3 The Generalized Hausman Test

In this section, we present the GH test, derived by White,  (1982), applied
here to detect non-normality of the latent variable using the SNP-IRT model.

As in the previous sections, let us denote by $\boldsymbol{\theta}$ the sub-
vector of
$\boldsymbol{\theta}^{{(1)}^{\prime}}=(\boldsymbol{\alpha}_{0}^{\prime},\boldsymbol{\alpha}_{1}^{\prime},\boldsymbol{\varphi}^{\prime})$
that includes the item intercepts $\boldsymbol{\alpha}_{0}$ and slopes
$\boldsymbol{\alpha}_{1}$. $\boldsymbol{\theta}$ has dimension $2p\times 1$,
where $p$ is the number of items.

For the 2PL IRT model ($SNP_{0}$), consider the maximum PL estimator
$\tilde{\boldsymbol{\theta}}_{SNP_{0}}$.

Consider the quasi-ML estimator
$\hat{\boldsymbol{\theta}}^{{(1)}^{\prime}}_{SNP_{L}}=({\hat{\boldsymbol{\theta}}^{{}^{\prime}}_{SNP_{L}}},\hat{\boldsymbol{\varphi}}^{{}^{\prime}})$
of a SNP-IRT model with $L>0$, where the sub-vector of parameter
$\hat{\boldsymbol{\varphi}}$ has dimension $L\times 1$ and so
$\hat{\boldsymbol{\theta}}_{SNP_{L}}^{{(1)}}$ has dimension $(2p+L)\times 1$.
Following White,  (1982), under normality of the latent variable

$\sqrt{n}(\hat{\boldsymbol{\theta}}_{SNP_{L}}-\tilde{\boldsymbol{\theta}}_{SNP_{0}})\xrightarrow[]{d}N(0,S(\boldsymbol{\theta}_{0},{\boldsymbol{\theta}_{0}}^{{(1)}})).$
(17)

An estimator of $S(\boldsymbol{\theta}_{0},{\boldsymbol{\theta}_{0}}^{{(1)}})$
is given by

$\small\begin{split}&\hat{S}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})=\hat{A}^{\boldsymbol{\theta}\boldsymbol{\varphi}}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1}\hat{B}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})\hat{A}^{\boldsymbol{\theta}\boldsymbol{\varphi}}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1^{\prime}}+\hat{A}(\tilde{\boldsymbol{\theta}}_{SNP_{0}})^{-1}\hat{B}(\tilde{\boldsymbol{\theta}}_{SNP_{0}})\hat{A}({\tilde{\boldsymbol{\theta}}_{SNP_{0}}})^{-1^{\prime}}-\\\
&-\hat{A}^{\boldsymbol{\theta}\boldsymbol{\varphi}}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1}\hat{R}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{\prime}\hat{A}(\tilde{\boldsymbol{\theta}}_{SNP_{0}})^{-1^{\prime}}-\hat{A}(\tilde{\boldsymbol{\theta}}_{SNP_{0}})^{-1}\hat{R}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})\hat{A}^{\boldsymbol{\theta}\boldsymbol{\varphi}}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1^{\prime}},\end{split}$
(18)

where the matrices $\hat{A}(\tilde{\boldsymbol{\theta}}_{SNP_{0}})$ and
$\hat{B}(\tilde{\boldsymbol{\theta}}_{SNP_{0}})$, defined in formulas (7) and
(8), have dimension $2p\times 2p$ and are evaluated at
$\tilde{\boldsymbol{\theta}}_{SNP_{0}}$.
$\hat{A}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$ and
$\hat{B}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$ are the observed
Hessian and cross-product matrix of dimension $(2p+L)\times(2p+L)$ for the
$SNP_{L}$ model, evaluated at $\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}}$.
The matrix
$\hat{A}^{\boldsymbol{\theta}\boldsymbol{\varphi}}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1}$
is obtained by deleting the last $L$ rows from the matrix
$\hat{A}(\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1}$ and has dimension
$2p\times(2p+L)$. The matrix
$\hat{R}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$
has dimension $2p\times(2p+L)$ and can be computed as

$\hat{R}(\boldsymbol{\theta}_{SNP_{0}},{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})=\sum_{i=1}^{n}\frac{\partial
pl_{{SNP_{0}}}(\boldsymbol{y}_{i},\boldsymbol{\theta})}{\partial{{\boldsymbol{\theta}}}}\frac{\partial
l_{{SNP_{L}}}(\boldsymbol{y}_{i},\boldsymbol{\theta}^{{(1)}})}{\partial{\boldsymbol{\theta}}^{{(1)^{\prime}}}},\hskip
8.5359pt$ (19)

where $pl_{{SNP_{0}}}(\boldsymbol{y}_{i},\boldsymbol{\theta})$ is the pairwise
log-likelihood for the individual $i$ under the model $SNP_{0}$ and
$l_{{SNP_{L}}}(\boldsymbol{y}_{i},\boldsymbol{\theta}^{{(1)}})$ is the log-
likelihood for the individual $i$ under the model $SNP_{L}$. The matrix in
(19) is evaluated at
$(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$.
We choose the maximum PL and the quasi-ML estimator for the two models to
avoid that, under correct model specification,
$\tilde{\boldsymbol{\theta}}_{SNP_{0}}$ and
$\hat{\boldsymbol{\theta}}_{SNP_{L}}$ converge to the same covariance matrix,
producing a
$\hat{S}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$
matrix in (18) with all entries close to 0.

Given the theoretical result in (17), the GH test is given by

$GH=(\hat{\boldsymbol{\theta}}_{SNP_{L}}-\tilde{\boldsymbol{\theta}}_{SNP_{0}})^{\prime}\hat{S}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})^{-1}(\hat{\boldsymbol{\theta}}_{SNP_{L}}-\tilde{\boldsymbol{\theta}}_{SNP_{0}}).$
(20)

Under normality of the latent variable, the GH test is asymptotically
distributed as a $\chi^{2}_{2p}$, with $2p$ degrees of freedom, i.e. the
number of parameters in $\boldsymbol{\theta}$.

We consider a simpler version of (20), that does not involve the inversion of
the matrix
$\hat{S}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$.
We consider the following quadratic form

$GH_{T}=(\hat{\boldsymbol{\theta}}_{SNP_{L}}-\tilde{\boldsymbol{\theta}}_{SNP_{0}})^{\prime}I_{2p}^{-1}(\hat{\boldsymbol{\theta}}_{SNP_{L}}-\tilde{\boldsymbol{\theta}}_{SNP_{0}}),$
(21)

where $I_{2p}$ can be omitted from the above formula.

Following Yuan and Bentler,  (2010)

$GH_{T}=\sum_{l=1}^{d}\lambda_{l}\delta_{l}^{2},\hskip
28.45274pt\delta_{l}\sim N(0,1),$ (22)

where $d$ is the rank of
$S({\boldsymbol{\theta}_{0}},{\boldsymbol{\theta}_{0}}^{{(1)}})$, and
$\lambda_{1},...,\lambda_{d}$ are its non-zero eigenvalues.

It is possible to approximate the distribution of $GH_{T}$ using the moment
matching method (Welch, , 1938,Yuan and Bentler, , 2010) as follows

$GH_{T}\sim a\chi^{2}_{b}.$ (23)

The quantity $a$ and $b$ are defined as

$a=\frac{\sum_{l=1}^{d}\lambda_{l}^{2}}{\sum_{l=1}^{d}\lambda_{l}}$ (24)

and

$b=\frac{(\sum_{l=1}^{d}\lambda_{l})^{2}}{\sum_{l=1}^{d}\lambda_{l}^{2}}.$
(25)

Since $S({\boldsymbol{\theta}_{0}},{\boldsymbol{\theta}_{0}}^{{(1)}})$ can be
consistently estimated by
$\hat{S}(\tilde{\boldsymbol{\theta}}_{SNP_{0}},\hat{\boldsymbol{\theta}}^{{(1)}}_{SNP_{L}})$
defined in (18), $a$ and $b$ can be consistently estimated substituting
$\hat{\lambda}_{1},...,\hat{\lambda}_{d}$ in (24) and (25), where $d$ is rank
of $\hat{S}$ and $\hat{\lambda}_{1},...,\hat{\lambda}_{d}$ are its non-zero
eigenvalues. The approximation in (23) matches the first two moments of
$GH_{T}$ with those of $a\chi^{2}_{b}$ (23).

##  4 Other goodness-of-fit tests

In the case of normality of the latent variable, the overall goodness-of-fit
of an IRT model for binary data is usually assessed through the Pearson’s chi-
square test and LR test (Bartholomew et al., , 2011). Since the $SNP_{L}$ and
$SNP_{0}$ model are nested, we can consider the LR test for nested models
(Wilks, , 1938). The $SNP_{L}$ model, if
$\varphi_{1}=\cdots={\varphi_{L}}=\pm\frac{\pi}{2}$, reduces to the $SNP_{0}$
model (Irincheeva et al., , 2012). For the computation of the LR test, the
$SNP_{0}$ model needs to be estimated using full-maximum likelihood (full-ML)
instead of PL, in order to obtain a comparable value of the log-likelihood
function with that of the $SNP_{1}$ model.

Let us denote by
$\boldsymbol{\varphi}^{\prime}=(\varphi_{1},\cdots,{\varphi_{L}})$ and by
$\textbf{c}^{\prime}=(\pm\frac{\pi}{2},\cdots,\pm\frac{\pi}{2})$. The null and
alternative hypotheses can be formulated as follows:

$H_{0}:\boldsymbol{\varphi}=\textbf{c}\hskip 8.5359ptvs\hskip
8.5359ptH_{1}:\boldsymbol{\varphi}\neq\textbf{c}.$ (26)

The test statistic is defined as

$LR=2\\{l_{SNP_{L}}(\mathbf{y},\hat{\boldsymbol{\theta}}^{{(1)}})-l_{SNP_{0}}(\mathbf{y},\hat{\boldsymbol{\theta}})\\},$
(27)

where $l_{SNP_{L}}(\mathbf{y},\hat{\boldsymbol{\theta}}^{{(1)}})$ and
$l_{SNP_{0}}(\mathbf{y},\hat{\boldsymbol{\theta}})$ are the quasi-log-
likelihood and full-log-likelihood functions of the $SNP_{L}$ and $SNP_{0}$
models, respectively, evaluated at their maximum values. Under $H_{0}$, the LR
test is asymptotically distributed as a $\chi^{2}_{L}$.

However, the LR test is affected by the problem of sparseness. Indeed, for a
fixed sample size, the number of empty cells in a frequency table increases
with the number of binary items. In this case, the distribution of the LR test
statistic is not well approximated by the chi-square distribution. To overcome
this problem, Maydeu-Olivares and Joe,  (2005) propose a family of test
statistics $M_{r}$, based on the residuals up to order $r$. The most popular
statistic is $M_{2}$, that uses the univariate and bivariate marginal
information. As data sparseness increases, the empirical Type I error rates of
the $M_{2}$ test remain accurate (Maydeu-Olivares and Joe, , 2005, 2006). As
before, let’s consider $p$ items and a sample size $n$. Under the null
hypothesis we test that the $SNP_{0}$ model holds. The hypotheses $H_{0}$ and
$H_{1}$ can be formalized as follows:

$H_{0}:\boldsymbol{\pi}=\boldsymbol{\pi}(\boldsymbol{\theta})\hskip
14.22636ptvs\hskip
14.22636ptH_{1}:\boldsymbol{\pi}\neq\boldsymbol{\pi}(\boldsymbol{\theta}),$
(28)

where $\boldsymbol{\theta}$, as usual, includes the item intercepts and slopes
and $\boldsymbol{\pi}(\boldsymbol{\theta})$ indicates the response patterns
probabilities.

The statistic $M_{2}$ is (Maydeu-Olivares and Joe, , 2005):

$M_{2}=n\hat{\textbf{e}}_{2}^{\prime}\hat{\boldsymbol{U}}_{2}\hat{\textbf{e}}_{2}.$
(29)

The vector $\hat{\textbf{e}}$ includes the univariate and bivariate residuals
while the matrix $\hat{\boldsymbol{U}}_{2}$ depends on a transformation matrix
and on the Jacobian matrix of the cell probabilities with respect to the items
intercept and slope parameter (more details can be found in Maydeu-Olivares
and Joe, , 2005, 2006). Under $H_{0}$, the statistic $M_{2}$ is asymptotically
distributed as a $\chi^{2}_{m}$, with degrees of freedom
$m=\frac{p(p+1)}{2}-2p$, that is the number of univariate and bivariate
residuals minus the number of estimated parameters of the $SNP_{0}$ model. In
the simulations, we evaluate the performance of the LR and $M_{2}$ tests under
non-normality of the latent variable.

##  5 Model selection criteria

The Akaike information criterion (AIC), the Bayesian information criterion
(BIC) and the Hannan–Quinn criterion (HQ) can be used to choose the degree of
the polynomial $L$ of the SNP-IRT model (Davidian and Gallant, , 1993, Woods
and Lin, , 2009, Irincheeva et al., , 2012, Monroe, , 2014).

The AIC is (Akaike, , 1974):

$AIC=-2l(\mathbf{y},\hat{\boldsymbol{\theta}}^{(1)})+2k,$ (30)

where $l(\mathbf{y},\hat{\boldsymbol{\theta}}^{(1)})$ is the quasi-log-
likelihood function of the $SNP_{L}$ model evaluated at the maximum value and
$k$ is the number of parameter of the model.

The BIC is (Schwarz, , 1978):

$BIC=-2l(\mathbf{y},\hat{\boldsymbol{\theta}}^{(1)})+k\ln{n},$ (31)

where $n$ is the sample size.

The HQ is (Hannan and Quinn, , 1979):

$HQ=-2l(\mathbf{y},\hat{\boldsymbol{\theta}}^{(1)})+2k\ln{\ln{n}}$ (32)

Usually $L=1$ or $L=2$ are enough to detect a departure from normality and to
approximate different shapes of latent variable distributions. Selecting
higher order of the polynomial could result in overfitting (Irincheeva, ,
2011).

When $L=0$, $\hat{\boldsymbol{\theta}}^{(1)}=\hat{\boldsymbol{\theta}}$ in the
above formulas and, as for the LR test, the $SNP_{0}$ model is estimated using
full-ML.

##  6 Simulation study

###  6.1 Design

In this section, we study the performance of the $GH_{T}$ test to assess non-
normality of the latent variable distribution and we compare its performance
with the $M_{2}$ and $LR$ tests. Moreover, for all simulation scenarios, the
BIC, AIC and HQ criteria have been computed.

To evaluate the performance of the $GH_{T}$, $M_{2}$ and $LR$ tests, we have
considered five scenarios (SC), corresponding to five different distribution
assumptions for the latent variable $z$ in the data generating models. The
general model is

$\begin{split}&logit(\pi_{ij})=\alpha_{0j}+\alpha_{1j}z_{i}\hskip
28.45274pti=1,...,n\hskip 28.45274ptj=1,2,...,p\\\ &z\sim f(z)\end{split}$
(33)

Item intercepts have been randomly chosen from the interval [-0.8; 1.12] while
the item slopes from the interval [0.5; 1.5]. To study the Type I error rates
of the $GH_{T}$, $M_{2}$ and LR tests we have considered the following
scenario:

  1. A

$z\sim N(0,1)$

To study the power of the $GH_{T}$, $M_{2}$ and LR tests we have considered
the following cases of two bimodal distributions and two skewed distributions:

  1. B

$z\sim 0.7N(-1,0.7)+0.3N(1,0.8),$

where $z$ has an overall mean equal to -0.40 and variance equal to 1.38.

  1. C

$z\sim 0.1N(-2,0.25)+0.9N(2,1),$

where $z$ has an overall mean equal to 1.6 and variance equal to 2.37.

  1. D

$Z\sim SN(\mu=-2.5,\sigma=2,\lambda=5)$,

where $z$ has mean -0.93 and variance 1.55.

  1. E

$Z\sim SN(\mu=-2.5,\sigma=2,\lambda=10)$,

where $z$ has mean -0.91 and variance 1.47.

The mixture of normals in scenario C wider departs from the normal
distribution than the one in scenario B. Similarly, the skew-normal
distribution in scenario E exhibits a greater departure from the normal
distribution than the one in D.

Two versions of the $GH_{T}$ test have been considered in the simulations. The
first version, denoted as $GH_{T1}$, is based on the $SNP_{0}$ and $SNP_{1}$
models. The choice of the $SNP_{1}$ model has been motivated by the fact that
it can approximate bimodal distributions well, such as the ones in scenarios B
and C.

The optimization of the $SNP_{1}$ model has been achieved in R with direct
maximization via the function “nlminb”, that uses the analytically computed
gradient and Hessian matrix. For the $SNP_{1}$ model, the initial values of
the parameters $\boldsymbol{\alpha}_{0}$ and $\boldsymbol{\alpha}_{1}$ used in
the optimization process are the full-ML parameter estimates obtained with the
$SNP_{0}$ model. In each data replication, for the $\varphi_{1}$ parameter, we
have sampled 10 initial values from a sequence of values, equally spaced by
0.1 in the interval $[-\frac{\pi}{2};\frac{\pi}{2}]$, i.e. the domain of
$\varphi_{1}$. Among the estimated $SNP_{1}$ models in each data replication,
we have selected the one that corresponds to the maximum value of the quasi-
log-likelihood function. The performance of the $GH_{T1}$ test has been
evaluated under all scenarios.

In scenarios D and E, a second version of $GH_{T}$ called $GH_{T2}$ has been
considered in addition to $GH_{T1}$. $GH_{T2}$ is based on the $SNP_{0}$ model
and the $SNP_{2}$ model. This choice is motivated by the fact that the SNP
model provides a representation that captures the highly skewed case
particularly well when $L=2$, as showed in Figure 2. As for the $SNP_{1}$
model, also the optimization of the $SNP_{2}$ model has been conducted in R
using the direct maximization function "nlminb". However, the $SNP_{2}$ model
is highly sensitive to the initial values chosen. To mitigate this, the true
parameter values have been used as initial values for the
$\boldsymbol{\alpha}_{0}$ and $\boldsymbol{\alpha}_{1}$ parameters. The
initial values for the $\varphi_{1}$ and $\varphi_{2}$ parameters have been
set to 0.7 and 1, respectively, which correspond to a skew-normal case. To
compute the $GH_{T1}$ and $GH_{T2}$ tests, the $SNP_{0}$ model has been
estimated using the maximum PL method through the R function "optim". The
Hessian, cross-product matrices and the matrix in formula (19) involved in the
$GH_{T1}$ and $GH_{T2}$ tests have been computed numerically with the
“NumDeriv” R package.

To compute the values of the information criteria and the $LR$ test, the
$SNP_{0}$ model has been estimated using the full-ML method. This estimation
has been achieved in R via the function "optim". The $M_{2}$ test in each data
replication has been computed using the "M2" function of the R package "mirt".

The Type I error rates and power of the $GH_{T1},GH_{T2}$, $M_{2}$ and LR
tests have been computed as follows:
$\hat{p}=\sum_{l=1}^{N_{v}}\frac{I(T_{l}\geq c)}{N_{v}}$, where $N_{v}$
represents the number of valid statistics out of the total replications. Here,
$I$ denotes an indicator function, and $T_{l}$ is the value of the test
statistic evaluated in the $l$-th replication. The critical value $c$
corresponds to the theoretical asymptotic value, specifically the
$(1-\alpha)$th percentile of the $\hat{a}\chi^{2}_{\hat{b}}$ distribution for
the $GH_{T}$ test. The values $\hat{a}$ and $\hat{b}$ are computed as in (24)
and (25). For the $M_{2}$ test, the critical value is associated with the
$\chi^{2}_{{m}}$ distribution, where $m$ is equal to $\frac{p(p+1)}{2}-2p$.

Considering that the $SNP_{1}$ model has one more parameter than the $SNP_{0}$
model, the LR test uses $c$ as the theoretical asymptotic critical value
corresponding to the $(1-\alpha)$th percentile of the $\chi^{2}_{{1}}$
distribution. This specific test is referred to as $LR_{1}$. For the
comparison of $SNP_{0}$ with $SNP_{2}$, the LR test with 2 degrees of freedom
is used, denoted as $LR_{2}$.

To establish the confidence interval (CI) of each rate $\hat{p}$, we calculate
it as $\hat{p}\pm
z_{(1-\frac{\alpha}{2})}\sqrt{\frac{\alpha(1-\alpha)}{N_{v}}}$.

Next, the percentages of times the AIC, BIC, and HQ criteria select the
$SNP_{1}$ over $SNP_{0}$ have been computed as
$\hat{P}=\sum_{l=1}^{N_{v}}\frac{I(IC_{SNP_{1}}<IC_{SNP_{0}})}{N_{v}}100$,
where IC indicates the AIC, BIC, and HQ criteria. Similarly, we computed
percentages for the selection of the $SNP_{2}$ model over the $SNP_{0}$ model
using information criteria.

For scenarios A, B and C, we have considered the following simulation
conditions: number of items $(p=10,20)$ $\times$ sample size
($n=500,1000,5000$)$\times$ test statistics ($GH_{T1}$, $M_{2}$,
$LR_{1}$)$\times$ information criterion (AIC, BIC, HQ). For scenarios D and E,
the simulation conditions are: number of items $(p=10,20)$ $\times$ sample
size ($n=500,1000$)$\times$ test statistics ($GH_{T1}$,$GH_{T2}$, $M_{2}$,
$LR_{1}$, $LR_{2}$)$\times$ information criterion (AIC, BIC, HQ). In all the
simulation scenarios, $R=500$ replications and two nominal levels of $\alpha$
have been considered, that is $\alpha=0.05,0.01$.

Preliminary results on the performance of the $GH_{T1}$ test under scenario A
and C have been obtained in Guastadisegni et al.,  (2023). The test showed
good performance both in terms of Type I error rates and power, especially for
large sample sizes.

Some additional results, that include the bias for the parameter estimates of
the $SNP_{0}$, $SNP_{1}$ and $SNP_{2}$ model, have been reported in the
Appendix B. If the mean of the true latent variable is different from 0 and
the variance different from 1, to compute the parameter bias, the true
$\boldsymbol{\alpha}_{00}$ and $\boldsymbol{\alpha}_{01}$ parameters have been
rescaled accordingly to formula (14) and (15), replacing
$\accentset{\mbox{\large.}}{\alpha}_{0j}$ with the true intercept,
$\accentset{\mbox{\large.}}{\alpha}_{1j}$ with the true slope, $\tilde{E}(Z)$
and $\tilde{V}(Z)$ with the mean and the variance of the true latent variable.
In this way, the parameter bias is due only to the misspecification of the
shape of the latent variable, and not to the misspecification of the moments.
A similar procedure to compute parameter bias has been adopted by Irincheeva,
(2011).

###  6.2 Results

Table 1 reports the Type I error rates of the $GH_{T1}$, $M_{2}$ and $LR_{1}$
tests for scenario A.

Table 1: Scenario A: Type I error rates of the $GH_{T1}$, $M_{2}$ and $LR_{1}$
tests, $p=10,20$, $n=500,1000,5000$

|  | $\alpha=0.05$ | $\alpha=0.01$
---|---|---|---
$p$ | $n$ | $GH_{T1}$ | $M_{2}$ | $LR_{1}$ | $GH_{T1}$ | $M_{2}$ | $LR_{1}$
10 | 500 | 0.018 | 0.066 | 0 | 0.002 | 0.012 | 0
| 1000 | 0.044 | 0.046 | 0.002 | 0.006 | 0.01 | 0
| 5000 | 0.056 | 0.04 | 0.042 | 0.02 | 0.002 | 0
20 | 500 | 0.026 | 0.036 | 0.006 | 0.008 | 0 | 0
| 1000 | 0.042 | 0.054 | 0.054 | 0.008 | 0.006 | 0
| 5000 | 0.042 | 0.05 | 0.264 | 0.01 | 0.016 | 0.094

11footnotetext:

Note 1: Values in boldface indicate that the nominal level $\alpha$ is not
included in their confidence interval

Overall, the $GH_{T1}$ test has good performance in terms of Type I error
rates under most conditions. Indeed, $GH_{T1}$ is more conservative than
expected only for $\alpha=0.05$, 10 and 20 items and small sample size.
$M_{2}$ has good performance in terms of Type I error rates for all values of
$\alpha$, number of items and sample sizes considered. Among the three tests
considered, the $LR_{1}$ test has the worst performance. For $\alpha=0.05$, it
rejects less than it should with 10 items, $n=500,1000$ and 20 items, $n=500$.
Moreover, for all level of $\alpha$ considered, it has seriously inflated Type
I error rates with 20 items and $n=5000$.

Table 2 shows the percentages of times AIC, BIC and HQ select $SNP_{0}$
instead of $SNP_{1}$ for scenario A.

Table 2: Scenario A: percentages of times AIC, BIC and HQ select $SNP_{0}$ instead of $SNP_{1}$, $p=10,20$, $n=500,1000,5000$ $p$ | $n$ | AIC | BIC | HQ
---|---|---|---|---
10 | 500 | 97% | 100% | 99.8%
| 1000 | 93.2 % | 100% | 98.8%
| 5000 | 90.6 % | 100% | 96.6 %
20 | 500 | 87.8% | 99.8% | 99.2%
| 1000 | 78.4% | 100% | 94.6 %
| 5000 | 54.4% | 97 % | 78.4%

Among the three information criteria considered, the AIC has the worst
performance. This is evident especially with 20 items and all sample sizes,
where it selects the $SNP_{0}$ model more or less from 54% to 88% of times.
The BIC has the best performance and it selects the $SNP_{0}$ model almost in
the totality of cases under all conditions. The performance of the HQ to
select the $SNP_{0}$ model is in between the performance of the AIC and BIC
criteria.

Table 3 presents the empirical power of the $GH_{T1}$, $M_{2}$ and $LR_{1}$
tests for scenarios B and C.

Table 3: Scenarios B and C: empirical power of the of the $GH_{T1}$, $M_{2}$ and $LR_{1}$ tests, $p=10,20$, $n=500,1000,5000$ |  |  | $\alpha=0.05$ | $\alpha=0.01$
---|---|---|---|---
SC | $p$ | $n$ | $GH_{T1}$ | $M_{2}$ | $LR_{1}$ | $GH_{T1}$ | $M_{2}$ | $LR_{1}$
B | 10 | 500 | 0.612 | 0.03 | 0.572 | 0.536 | 0.004 | 0.462
|  | 1000 | 0.868 | 0.028 | 0.792 | 0.796 | 0.006 | 0.558
|  | 5000 | 1 | 0.052 | 0.936 | 1 | 0.014 | 0.93
| 20 | 500 | 0.886 | 0.022 | 0.596 | 0.76 | 0.006 | 0.428
|  | 1000 | 0.984 | 0.036 | 0.63 | 0.96 | 0.002 | 0.60
|  | 5000 | 0.996 | 0.034 | 0.678 | 0.996 | 0.004 | 0.666
C | 10 | 500 | 0.924 | 0.01 | 0.912 | 0.852 | 0.002 | 0.8
|  | 1000 | 0.998 | 0.006 | 0.968 | 0.994 | 0 | 0.96
|  | 5000 | 1 | 0.028 | 0.972 | 1 | 0.004 | 0.972
| 20 | 500 | 0.988 | 0 | 0.782 | 0.95 | 0 | 0.72
|  | 1000 | 0.996 | 0 | 0.822 | 0.994 | 0 | 0.756
|  | 5000 | 0.998 | 0.004 | 0.922 | 0.998 | 0 | 0.892

Among the three tests and under all levels of $\alpha$ considered, $GH_{T1}$
has the highest power when the true latent variable is generated from a
mixture of normals, that is under scenarios B and C. In general, the power of
the $GH_{T1}$ and $LR_{1}$ tests increases as the sample size and the number
of items increase and as the true latent variable distribution largely departs
from the normal one. For what concerns the $M_{2}$ test, it has very low or no
power to detect non-normality of the latent variable distribution under both
scenarios considered.

Table 4 shows the percentages of times AIC, BIC and HQ select $SNP_{1}$
instead of $SNP_{0}$ for scenarios B and C.

Table 4: Scenarios B and C: percentages of times AIC, BIC and HQ select $SNP_{1}$ instead of $SNP_{0}$, $10,20$, $n=500,1000,5000$ SC | $p$ | $n$ | AIC | BIC | HQ
---|---|---|---|---|---
B | 10 | 500 | 75.6 % | 47.6% | 58.8 %
|  | 1000 | 85.2% | 55.4 % | 78.8 %
|  | 5000 | 94% | 93% | 93.4%
| 20 | 500 | 66.2 % | 47% | 60.2 %
|  | 1000 | 65.2% | 60% | 63%
|  | 5000 | 70% | 65.4 % | 67.6 %
C | 10 | 500 | 95.6% | 80.6% | 91.8%
|  | 1000 | 97% | 95.6% | 96.8%
|  | 5000 | 97.6% | 97.2% | 97.2%
| 20 | 500 | 85.8% | 73% | 79.4%
|  | 1000 | 87.6% | 75.4% | 82%
|  | 5000 | 95% | 87.2% | 92 %

Under all scenarios, AIC has the best performance when $n=500,1000$ and it has
the highest percentages of selections as the $SNP_{1}$ model. Under scenarios
B and C, AIC has the best performance also for $n=5000$. In the majority of
cases, BIC has the worst performance when $n=500,1000$. This is evident
especially under scenario B and $n=500$, where it selects the $SNP_{1}$ models
only around 47% of times. The performance of HQ to select the $SNP_{1}$ model
under all scenarios for small sample size is in between the performance of AIC
and BIC criteria.

Table 5 reports the empirical power of the $GH_{T1}$, $LR_{1}$, $GH_{T2}$,
$LR_{2}$ and $M_{2}$ tests for scenarios D and E.

Table 5: Scenarios D and E: empirical power of the $GH_{T1}$, $LR_{1}$, $GH_{T2}$, $LR_{2}$ and $M_{2}$ tests, $p=10,20$, $n=500,1000$ |  |  | $\alpha=0.05$ | $\alpha=0.01$
---|---|---|---|---
SC | $p$ | $n$ | $GH_{T1}$ | $LR_{1}$ | $GH_{T2}$ | $LR_{2}$ | $M_{2}$ | $GH_{T1}$ | $LR_{1}$ | $GH_{T2}$ | $LR_{2}$ | $M_{2}$
D | 10 | 500 | 0.554 | 0.374 | 0.492 | 0.7 | 0.024 | 0.416 | 0.27 | 0.338 | 0.462 | 0.004
|  | 1000 | 0.764 | 0.588 | 0.722 | 0.956 | 0.028 | 0.636 | 0.264 | 0.608 | 0.85 | 0.008
| 20 | 500 | 0.786 | 0.358 | 0.89 | 0.972 | 0.026 | 0.576 | 0.196 | 0.762 | 0.876 | 0.008
|  | 1000 | 0.962 | 0.38 | 0.992 | 1 | 0.02 | 0.924 | 0.348 | 0.972 | 1 | 0.002
E | 10 | 500 | 0.588 | 0.43 | 0.562 | 0.818 | 0.028 | 0.458 | 0.334 | 0.404 | 0.62 | 0.01
|  | 1000 | 0.894 | 0.736 | 0.814 | 0.98 | 0.036 | 0.834 | 0.448 | 0.688 | 0.92 | 0.014
| 20 | 500 | 0.874 | 0.45 | 0.916 | 0.99 | 0.022 | 0.71 | 0.282 | 0.836 | 0.946 | 0.006
|  | 1000 | 0.986 | 0.474 | 1 | 1 | 0.026 | 0.948 | 0.458 | 0.998 | 1 | 0.026

Under scenarios D and E, with 10 items, $GH_{T1}$ exhibits a slightly higher
power compared to $GH_{T2}$. As reported in the Appendix B, the $SNP_{1}$ and
$SNP_{2}$ models result in a similar reduction in parameter bias compared to
$SNP_{0}$. However, $SNP_{2}$ better approximates the shape of the true latent
variable distribution (see Appendix B) and with 20 items, the power of the
$GH_{T2}$ test is slightly higher compared to $GH_{T1}$. In general, both
tests achieve a high power with large sample sizes and as the skew-normal
becomes more extreme. $LR_{1}$ has low power under all conditions, while
$LR_{2}$ exhibits the highest power among all scenarios. Although increasing
the degree of the polynomial can enhance the $LR$ and $GH_{T}$ tests
performance, the $SNP_{2}$ method requires accurate initial values for the
parameters to yield good results, making it challenging to use in practice. As
observed in scenarios B and C, the $M_{2}$ test has very low or no power in
detecting misspecifications of the latent variable distribution.

Table 6 shows the percentages of times AIC, BIC and HQ select $SNP_{1}$ over
$SNP_{0}$ and $SNP_{2}$ over $SNP_{0}$ for scenarios D and E.

Table 6: Scenarios D and E: percentages of times AIC, BIC and HQ select $SNP_{1}$ over $SNP_{0}$ and $SNP_{2}$ over $SNP_{0}$, $p=10,20$, $n=500,1000$ |  |  |  $SNP_{1}$ over $SNP_{0}$ |  $SNP_{2}$ over $SNP_{0}$ |
---|---|---|---|---|---
SC | $p$ | $n$ | AIC | BIC | HQ | AIC | BIC | HQ |
D | 10 | 500 | 58.4% | 28.4% | 39% | 87% | 26.2% | 60.4% |
|  | 1000 | 70.4% | 24.4% | 58.4% | 99% | 61.8% | 90% |
| 20 | 500 | 44.8% | 22% | 36.6% | 98.8% | 75.4% | 94% |
|  | 1000 | 41.4% | 34.6% | 38% | 100% | 98.2% | 100% |
E | 10 | 500 | 64.6% | 34.4% | 44.2% | 91.6% | 42.2% | 72.6% |
|  | 1000 | 81.2% | 44% | 73.2% | 99% | 76% | 95% |
| 20 | 500 | 51.2% | 31% | 45.2% | 99.6% | 84.4% | 97.8% |
|  | 1000 | 49.6% | 45.6% | 47.4% | 100% | 99.6% | 100% |

None of the criteria have good performance in selecting $SNP_{1}$ over
$SNP_{0}$ for both scenarios and considered sample sizes. However, since the
true latent variables are skewed and $SNP_{2}$ better approximates this type
of distributions, the performance of all the information criteria improves
when selecting between $SNP_{0}$ and $SNP_{2}$, with AIC showing the best
performance and BIC performing the worst. Overall, from the results in Table
2, Table 4 and Table 6, it has not been possible to identify the criterion
that has the best performance under normality and non-normality of the latent
variable.

##  7 Real data application: American students exposure to school and
neighbourhood violence

This research is based on data from the National Longitudinal Survey of
Freshmen (NLSF), a project designed by Douglas S. Massey and Camille Z.
Charles and funded by the Mellon Foundation and the Atlantic Philanthropies
(available at
http://oprdata.princeton.edu/archive/restricted). The aim of the NLSF is to
collect data to explain the minority underachievement in higher education.
Data have been collected from 1999 to 2003, in four waves, to capture emergent
psychological processes, measuring the degree of social integration and
intellectual engagement. The survey included equal-sized samples of white,
black, Asian, and Latino freshmen entering selective colleges and
universities. We have analyzed only a part of the questionnaire that refers to
the year 1999. In particular we have selected 9 binary items that measure
violence in the neighbourhood. The items description is reported in Table 7.

Table 7: NLSF data: item description Item | Question
---|---
1 | In your neighborhood, before you were ten
| do you remember seeing homeless people on the street?
2 | Prostitutes on street?
3 | Gang members hanging out on the street?
4 | Drug paraphernalia on the street?
5 | People selling illegal drugs in public?
6 | People using illegal drugs in public?
7 | People drinking or drunk in public?
8 | Physical violence in public?
9 | Hearing the sound of gunshots?

The original sample of observations is composed by 3924 observations. Possible
responses are “no”, “yes”, “don’t know” and “refused”. Individuals that have
responded “don’t know” or “refused” have been excluded from the analysis,
while responses “no” have been coded as 0 and “yes” as 1. The data set that
has been analyzed is finally composed by 3891 observations. A sub-sample of
400 observations and a larger number of items of the questionnaire that refers
to the year 1999 has been analyzed also by Cagnone and Viroli,  (2012), that
fitted a latent trait models with two factors distributed as a mixture of
normals. They found that the item reported in Table 7 are highly loaded only
on one factor, distributed as a mixture of normals.

Some preliminary descriptive analysis on this dataset have been performed with
the “ltm” R package. Even if not reported in the tables, all items are
statistically associated and in all the items the proportion of “0” responses
is higher than “1”.

The first step of the analysis has been to fit the $SNP_{0}$ and $SNP_{1}$
models to the data. Table 8 reports the full-ML and quasi-ML parameter
estimates of the $SNP_{0}$ and $SNP_{1}$ models, respectively, and related
standard errors, based on the sandwich covariance matrix.

Table 8: NLSF data: $SNP_{0}$ and $SNP_{1}$ parameter estimates and related standard errors  | Estimates
---|---
Parameter | $SNP_{0}$ | $SNP_{1}$
$\alpha_{01}$ | -1.85 (0.07) | -2.61 (0.13)
$\alpha_{02}$ | -5.21 (0.22) | -6.56 (0.34)
$\alpha_{03}$ | -3.43 (0.14) | -4.88 (0.25)
$\alpha_{04}$ | -5.05 (0.27) | -7.49 (0.52)
$\alpha_{05}$ | -6.85 (0.48) | -9.43 (0.68)
$\alpha_{06}$ | -5.75 (0.31) | -7.76 (0.46)
$\alpha_{07}$ | -1.61 (0.09) | -2.78 (0.17)
$\alpha_{08}$ | -1.99 (0.09) | -2.99 (0.16)
$\alpha_{09}$ | -2.76 (0.10) | -3.68 (0.18)
$\alpha_{11}$ | 1.94 (0.09) | 2.93 (0.18)
$\alpha_{12}$ | 2.40 (0.15) | 4.09 (0.28)
$\alpha_{13}$ | 2.84 (0.15) | 4.64 (0.26)
$\alpha_{14}$ | 3.77 (0.23) | 6.71 (0.47)
$\alpha_{15}$ | 4.56 (0.37) | 7.75 (0.57)
$\alpha_{16}$ | 3.38 (0.22) | 5.83 (0.38)
$\alpha_{17}$ | 2.77 (0.15) | 4.16 (0.26)
$\alpha_{18}$ | 2.32 (0.11) | 3.58 (0.21)
$\alpha_{19}$ | 2.00 (0.10) | 3.19 (0.19)
$\varphi_{1}$ | - | 0.23(0.04)

Despite both methods being on the same scale, we notice dissimilar parameter
estimates. Irincheeva,  (2011) highlights that when the $SNP_{0}$ and
$SNP_{L}$ methods are employed on real binary data with true latent variables
that may deviate from a normal distribution, they yield significantly
different parameter estimates. The $SNP_{0}$ model might not be sufficient to
effectively capture the data, resulting in notable bias in both the item
intercepts and slopes (Ma and Genton, , 2010).

To choose the best model for the data, first we have evaluated the fit of the
$SNP_{0}$ model.

Since the data are sparse and 174 observed response patterns out of the total
231 observed response patterns have expected frequencies less than 5,
residuals calculated from marginal frequencies have been inspected. We have
considered the rule of thumb that residuals greater than 4 are indicators of
bad fit of the correspondent pair or triplets of items (Bartholomew et al., ,
2011). Even if not reported in the tables, the $SNP_{0}$ model does not have a
good fit for some pairs and triplets of items. To evaluate if the $SNP_{1}$
model has a better fit than the $SNP_{0}$ model, information criteria have
been computed.

Table 9 reports the values of the AIC, BIC and HQ criteria for the $SNP_{0}$
and $SNP_{1}$ models.

Table 9: NLSF data: Information criteria for the $SNP_{0}$ and $SNP_{1}$ models | AIC | BIC | HQ
---|---|---|---
$SNP_{0}$ | 21309.32 | 21422.12 | 21349.36
$SNP_{1}$ | 21303.67 | 21422.73 | 21345.93

The information criteria give conflicting results. Indeed AIC and HQ select
the $SNP_{1}$ model, while BIC the $SNP_{0}$ model.

In coherence with the simulation study, the $M_{2}$, $LR_{1}$ and $GH_{T1}$
tests have been computed.

Table 10 reports the value of the $M_{2}$, $LR_{1}$ and $GH_{T1}$ statistics,
the degrees of freedom (d.o.f) and the associated $p$-values.

Table 10: NLSF data: $M_{2}$, $LR_{1}$ and $GH_{T1}$ test statistics and associated degrees of freedom and $p$-values Test | Value | d.o.f |  $p$-value
---|---|---|---
$M_{2}$ | 182.33 | 27 | 0
$LR_{1}$ | 7.66 | 1 | 0.006
$GH_{T1}$ | 222.92 | 2.65 | 0

According to the $M_{2}$ test, the $SNP_{0}$ model does not have a good fit to
the data. However, the $M_{2}$ test does not reveal the source of misfit. It
could be the non-normality of the latent variable even if, as showed in the
simulation study, the $M_{2}$ test has a very low power to detect this type of
model misspecification, or other types of model misspecification. According to
the $LR_{1}$ and $GH_{T1}$ tests, the null hypothesis that the latent variable
is normally distributed has been rejected. This result is coherent with the
one of Cagnone and Viroli,  (2012) and with the simulation study, where both
the $GH_{T1}$ and $LR_{1}$ tests show good performance in terms of power to
detect non-normality of the latent variable distribution with many items and
large sample sizes. However, the result of $GH_{T1}$ is more reliable than
$LR_{1}$ because $GH_{T1}$ never has inflated Type I error rates with large
sample sizes. In this data analysis, we did not consider $GH_{T2}$ and
$LR_{2}$ tests. This decision was made because both tests already reject the
null hypothesis when considering only $L=1$, and there is also the issue of
the initial values used in the optimization process.

##  8 Discussion

In this work, we have extended the use of the GH test to detect non-normality
of the latent variable distribution in unidimensional IRT model for binary
data. The GH test has been obtained by comparing the PL estimator of the 2PL
IRT model for binary data with the quasi-ML estimator of the SNP-IRT model,
which allows for a more flexible shape of the latent variable distribution. A
simpler version than the GH test, referred to as the $GH_{T}$ test, has been
employed because it does not require the inversion of the covariance matrix.
Its approximated distribution has been derived using the moment matching
method. Two versions of $GH_{T}$ have been considered. First, we evaluated the
performance of the $GH_{T1}$ test based on the $SNP_{1}$ model through a
simulation study and real data analysis. We have compared the performance of
the $GH_{T1}$ test with the $M_{2}$ and the $LR_{1}$ tests and computed
information criteria, such as AIC, BIC, and HQ. The simulation study has shown
that the $GH_{T1}$ test has good performance in terms of Type I error rates
under most conditions. It also exhibits the highest power to detect non-
normality of the latent variable distribution when the true latent variable
followed a mixture of normals. The $LR_{1}$ test has inflated or deflated Type
I error rates under certain conditions and consistently lower power compared
to the $GH_{T1}$ test. Considering that the SNP density when $L=1$ can only
capture bimodal and slightly skew-normal distributions, we have adopted the
SNP density with $L=2$ in the context of skew-normal distributions. This
particular parameterization enables the representation of highly skew-normal
distributions for specific combinations of parameter values. For this reason,
we have employed the $GH_{T2}$ test based on the $SNP_{2}$ model, along with
the $LR_{2}$ test. The $GH_{T1}$ and $GH_{T2}$ tests perform similarly under
the skew-normal scenarios and exhibit high power with large sample sizes. The
$LR_{2}$ test has the highest power. However, both the $GH_{T2}$ and $LR_{2}$
tests required the estimation of the $SNP_{2}$ model. This poses a challenge
since accurate initial values for the true parameter values are required in
the optimization process to attain a reliable approximation of the true latent
variable and minimize bias in parameter estimates compared to $SNP_{0}$. As a
result, the application of the $SNP_{2}$ model in the analysis of real data
becomes impractical due to these demanding requirements. The $M_{2}$ test has
good performance in terms of Type I error rates but very low or no power to
detect the non-normality of the latent variable. Similar results on the low
power of the $M_{2}$ test to detect non-normality of the latent variable
distribution have been found by Paek et al.,  (2019). Also Ranger and Much,
(2020) have found that the $M_{2}$ test has very low or no power when the
misspecification was in the form of an upper boundary item characteristic
function. From the simulations it has not been possible to identify the
information criterion that has the best performance both under normality and
non-normality of the latent variable. AIC tends to select the more complex
model in the highest percentages of cases, regardless of the underlying
distribution of the latent variable. This is because it tends to favor models
with more parameters, which may lead to overfitting. On the other hand, BIC
has the best performance under normality of the latent variable, as it
penalizes the number of parameters more strongly than AIC, which makes it more
suitable for selecting parsimonious models. However, BIC has the worst
performance under non-normality of the latent variable for small sample sizes,
as it selects overly simple models that do not capture the complexity of the
data. The performance of the HQ criterion is in between the one of AIC and
BIC. The results from the real data analysis confirmed the findings from the
simulation study, where the $GH_{T1}$ test proved to be very useful in
detecting non-normality of the latent variable when the information criteria
yielded contradictory results. Overall, the $GH_{T1}$ test emerged as the most
powerful tool for detecting non-normality of the latent variable distribution.

Further studies could focus on addressing the issue of initial values in the
SNP estimation when $L>1$, making it more applicable in practical contexts.
The performance of the GH test implemented with higher-order polynomials could
also be evaluated by means of simulations and in real data analysis.
Increasing the degree of the polynomial allows for more flexibility in
modeling the shape of the latent variable distribution and can enhance both
information criteria and test performance.

In addition, the performance of the GH test in the IRT context could include
other types of model violations, as local dependence or violation of the item
characteristic function. In these cases, other types of estimation methods
consistent under model misspecification should be considered in order to apply
the test.

## References

  * Akaike, (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE transactions on automatic control, 19(6):716–723.
  * Bartholomew et al., (2011) Bartholomew, D. J., Knott, M., and Moustaki, I. (2011).  Latent Variable Models and Factor Analysis: A Unified Approach.  Wiley, 3 edition.
  * Bartolucci et al., (2017) Bartolucci, F., Bacci, S., and Pigini, C. (2017).  Misspecification test for random effects in generalized linear finite-mixture models for clustered binary and ordered data.  Econometrics and Statistics, 3:112–131.
  * Cagnone and Viroli, (2012) Cagnone, S. and Viroli, C. (2012).  A factor mixture analysis model for multivariate binary data.  Statistical Modelling, 12(3):257–277.
  * Cramér, (1946) Cramér, H. (1946).  Mathematical Methods of Statistics.  Princeton University Press.
  * Davidian and Gallant, (1993) Davidian, M. and Gallant, A. R. (1993).  The nonlinear mixed effects model with a smooth random effects density.  Biometrika, 80(3):475–488.
  * Gallant and Nychka, (1987) Gallant, A. R. and Nychka, D. W. (1987).  Semi-nonparametric maximum likelihood estimation.  Econometrica, 55(2):363–390.
  * Gallant and Tauchen, (1989) Gallant, A. R. and Tauchen, G. (1989).  Seminonparametric estimation of conditionally constrained heterogeneous processes: Asset pricing applications.  Econometrica, 57(5):1091–1120.
  * Guastadisegni et al., (2022) Guastadisegni, L., Cagnone, S., Moustaki, I., and Vasdekis, V. (2022).  Use of the Lagrange multiplier test for assessing measurement invariance under model misspecification.  Educational and Psychological Measurement, 82(2):254–280.
  * Guastadisegni et al., (2023) Guastadisegni, L., Moustaki, I., Vasdekis, V., and Cagnone, S. (2023).  Detecting latent variable non-normality through the generalized hausman test.  In Wiberg, M., Molenaar, D., González, J., Kim, J.-S., and Hwang, H., editors, Quantitative Psychology, pages 107–118, Cham. Springer Nature Switzerland.
  * Hannan and Quinn, (1979) Hannan, E. J. and Quinn, B. G. (1979).  The determination of the order of an autoregression.  Journal of the Royal Statistical Society: Series B (Methodological), 41(2):190–195.
  * Hausman, (1978) Hausman, J. A. (1978).  Specification tests in econometrics.  Econometrica, 46(6):1251–1271.
  * Irincheeva, (2011) Irincheeva, I. (2011).  Generalized linear latent variable models with flexible distributions.  PhD thesis, University of Geneva.
  * Irincheeva et al., (2012) Irincheeva, I., Cantoni, E., and Genton, M. G. (2012).  Generalized linear latent variable models with flexible distribution of latent variables.  Scandinavian Journal of Statistics, 39(4):663–680.
  * Katsikatsou et al., (2012) Katsikatsou, M., Moustaki, I., Yang-Wallentin, F., and Jöreskog, K. G. (2012).  Pairwise likelihood estimation for factor analysis models with ordinal data.  Computational Statistics & Data Analysis, 56(12):4243–4258.
  * Knott and Tzamourani, (2007) Knott, M. and Tzamourani, P. (2007).  Bootstrapping the estimated latent distribution of the two-parameter latent trait model.  British Journal of Mathematical and Statistical Psychology, 60(1):175–191.
  * Lindsay, (1988) Lindsay, B. G. (1988).  Composite likelihood methods.  Contemporary mathematics, 80(1):221–239.
  * Lord, (1980) Lord, F. M. (1980).  Applications of Item Response Theory to Practical Testing Problems.  Routledge, 1 edition.
  * Ma and Genton, (2010) Ma, Y. and Genton, M. G. (2010).  Explicit estimating equations for semiparametric generalized linear latent variable models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(4):475–495.
  * Maydeu-Olivares and Joe, (2005) Maydeu-Olivares, A. and Joe, H. (2005).  Limited-and full-information estimation and goodness-of-fit testing in $2^{n}$ contingency tables.  Journal of the American Statistical Association, 100(471):1009–1020.
  * Maydeu-Olivares and Joe, (2006) Maydeu-Olivares, A. and Joe, H. (2006).  Limited information goodness-of-fit testing in multidimensional contingency tables.  Psychometrika, 71(4):713–732.
  * Monroe, (2014) Monroe, S. L. (2014).  Multidimensional item factor analysis with semi-nonparametric latent densities.  PhD thesis, UCLA.
  * Montanari and Viroli, (2010) Montanari, A. and Viroli, C. (2010).  A skew-normal factor model for the analysis of student satisfaction towards university courses.  Journal of Applied Statistics, 37(3):473–487.
  * Paek et al., (2019) Paek, I., Xu, J., and Lin, Z. (2019).  Detection rates of the m2 test for nonzero lower asymptotes under normal and nonnormal ability distributions in the applications of irt.  Applied psychological measurement, 43(1):84–88.
  * Ranger and Much, (2020) Ranger, J. and Much, S. (2020).  Analyzing the fit of IRT models with the Hausman test.  Frontiers in Psychology, 11.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the Dimension of a Model.  The Annals of Statistics, 6(2):461 – 464.
  * Skrondal and Rabe-Hesketh, (2004) Skrondal, A. and Rabe-Hesketh, S. (2004).  Generalized latent variable modeling: Multilevel, longitudinal, and structural equation models.  CRC Press.
  * van der Linden and Barrett, (2016) van der Linden, W. J. and Barrett, M. D. (2016).  Linking item response model parameters.  Psychometrika, 81(3):650–673.
  * van der Linden and Hambleton, (2013) van der Linden, W. J. and Hambleton, R. K. (2013).  Handbook of modern item response theory.  Springer Science & Business Media.
  * Varin, (2008) Varin, C. (2008).  On composite marginal likelihoods.  Advances in Statistical Analysis, 92(1):1–28.
  * Vasdekis et al., (2012) Vasdekis, V. G., Cagnone, S., and Moustaki, I. (2012).  A composite likelihood inference in latent variable models for ordinal longitudinal responses.  Psychometrika, 77:425–441.
  * Welch, (1938) Welch, B. L. (1938).  The significance of the difference between two means when the population variances are unequal.  Biometrika, 29(3-4):350–362.
  * White, (1982) White, H. (1982).  Maximum likelihood estimation of misspecified models.  Econometrica, 50(1):1–25.
  * Wilks, (1938) Wilks, S. S. (1938).  The large-sample distribution of the likelihood ratio for testing composite hypotheses.  The annals of mathematical statistics, 9(1):60–62.
  * Woods, (2006) Woods, C. M. (2006).  Ramsay-curve item response theory (RC-IRT) to detect and correct for nonnormal latent variables.  Psychological Methods, 11(3):253–270.
  * Woods and Lin, (2009) Woods, C. M. and Lin, N. (2009).  Item response theory with estimation of the latent density using Davidian curves.  Applied Psychological Measurement, 33(2):102–117.
  * Yuan and Bentler, (2010) Yuan, K.-H. and Bentler, P. M. (2010).  Two simple approximations to the distributions of quadratic forms.  British Journal of Mathematical and Statistical Psychology, 63(2):273–291.
  * Zhang and Davidian, (2001) Zhang, D. and Davidian, M. (2001).  Linear mixed models with flexible distributions of random effects for longitudinal data.  Biometrics, 57(3):795–802.

##  Appendix A The mean and variance of the SNP latent variable

To compute the final estimator $\hat{\boldsymbol{\alpha}}_{0}$ in formula (14)
and $\hat{\boldsymbol{\alpha}}_{1}$ in (15), it is necessary to compute
$\tilde{E}(Z)$ and $\tilde{V}(Z)$ for the latent variable with density in (2).
These quantities can be derived analytically.

These quantity are derived for $L=2$. After the optimization process,
$P_{L}(z)=a_{0}+a_{1}z+a_{2}z^{2}$ and
$h(z|\hat{\boldsymbol{\varphi}})=P_{L}^{2}(z)\phi(z)$, where
$a_{0}=sin\hat{\varphi}_{1}-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}$,
$a_{1}=cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}$,
$a_{2}=\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}$.

From Zhang and Davidian,  (2001)

$\tilde{E}(Z)=a^{\prime}M^{*}a$ (34)

where the element in the $i$-th row and $j$-th column of $M^{*}$ is
$E(z^{i+j-1})$ and $z\sim N(0,1)$. The matrix $M^{*}$ includes the moment of a
standard normal distribution.

When $L=2$

$M^{*}=E\begin{pmatrix}z&z^{2}&z^{3}\\\ z^{2}&z^{3}&z^{4}\\\
z^{3}&z^{4}&z^{5}\\\ \end{pmatrix}=\begin{pmatrix}0&1&0\\\ 1&0&3\\\ 0&3&0\\\
\end{pmatrix}$ (35)

and

$\small\begin{split}\tilde{E}(Z)&=\begin{pmatrix}sin\hat{\varphi}_{1}-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}&cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}&\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}\begin{pmatrix}0&1&0\\\
1&0&3\\\ 0&3&0\\\
\end{pmatrix}\begin{pmatrix}sin\hat{\varphi}_{1}-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\\\
cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}\\\
\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}=\\\
&=\begin{pmatrix}cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}&sin\hat{\varphi}_{1}+\frac{2}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}&3cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}\end{pmatrix}\begin{pmatrix}sin\hat{\varphi}_{1}-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\\\
cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}\\\
\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}=\\\
&=2sin\hat{\varphi}_{1}cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}+\frac{4}{\sqrt{2}}{cos\hat{\varphi}_{1}}^{2}cos\hat{\varphi}_{2}sin\hat{\varphi}_{2}\end{split}$
(36)

To compute $\tilde{V}(Z)$ we need also $\tilde{E}(Z^{2})$. It can be computed
as as $a^{\prime}M^{**}a$, where the element in the $i$-th row and $j$-th
column of $M^{**}$ is $E(z^{i+j})$, and $z\sim N(0,1)$ ( Zhang and Davidian, ,
2001). When $L=2$

$M^{**}=E\begin{pmatrix}z^{2}&z^{3}&z^{4}\\\ z^{3}&z^{4}&z^{5}\\\
z^{4}&z^{5}&z^{6}\\\ \end{pmatrix}=\begin{pmatrix}1&0&3\\\ 0&3&0\\\ 3&0&15\\\
\end{pmatrix}$ (37)

and

$\small\begin{split}\tilde{E}(Z^{2})&=\begin{pmatrix}sin\hat{\varphi}1-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}&cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}&\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}\begin{pmatrix}1&0&3\\\
0&3&0\\\ 3&0&15\\\
\end{pmatrix}\begin{pmatrix}sin\hat{\varphi}_{1}-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\\\
cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}\\\
\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}=\\\
&=\begin{pmatrix}sin\hat{\varphi}_{1}+\frac{2}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}&3cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}&3sin\hat{\varphi}_{1}+\frac{12}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}\begin{pmatrix}sin\hat{\varphi}_{1}-\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\\\
cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}\\\
\frac{1}{\sqrt{2}}cos\hat{\varphi}_{1}cos\hat{\varphi}_{2}\end{pmatrix}=\\\
&=sin{\hat{\varphi}_{1}}^{2}+\frac{4}{\sqrt{2}}cos\hat{\varphi}_{1}sin\hat{\varphi}_{1}cos\hat{\varphi}_{2}+3{cos\hat{\varphi}_{1}}^{2}sin{\hat{\varphi}_{2}}^{2}+5{cos\hat{\varphi}_{1}}^{2}cos{\hat{\varphi}_{2}}^{2}\end{split}$
(38)

The variance of the latent variable with a SNP density is computed as
$\tilde{V}(Z)=\tilde{E}(Z^{2})-\tilde{E}(Z)^{2}$, and we obtain

$\begin{split}\tilde{V}(Z)&=sin{\hat{\varphi}_{1}}^{2}+\frac{4}{\sqrt{2}}cos\hat{\varphi}_{1}sin\hat{\varphi}1cos\hat{\varphi}_{2}+3cos{\hat{\varphi}_{1}}^{2}sin{\hat{\varphi}_{2}}^{2}+5cos{\hat{\varphi}_{1}}^{2}cos{\hat{\varphi}_{2}}^{2}-\\\
&-(2sin\hat{\varphi}_{1}cos\hat{\varphi}_{1}sin\hat{\varphi}_{2}+\frac{4}{\sqrt{2}}cos{\hat{\varphi}_{1}}^{2}cos\hat{\varphi}_{2}sin\hat{\varphi}_{2})^{2}\end{split}$
(39)

To get the mean and the variance of the latent variable for the $SNP_{1}$
model, $\hat{\varphi}_{2}$ should be set equal to $\frac{\pi}{2}$in equations
(36) and (39).

##  Appendix B Additional simulation results

###  B.1 Parameter bias for the $\boldsymbol{SNP_{0}}$, $\boldsymbol{SNP_{1}}$
and $\boldsymbol{SNP_{2}}$ models

In this section, we report the bias observed in the quasi-ML estimates of the
$SNP_{1}$ model and both full-ML and PL estimates of the $SNP_{0}$ under
scenario C. Additionally, for scenario E, we also report the bias of the
quasi-ML estimates of the $SNP_{2}$ model. In each scenario considered, the
mean bias of each model parameter, indicated with $\theta$, is computed as:

$Bias_{\hat{\theta}}=\frac{\sum_{l=1}^{N}|\hat{\theta}_{l}-{\theta}_{0}|}{N},$

where $\theta_{0}$ is the true parameter, $\hat{\theta}_{l}$ is the estimate
of the parameter $\theta$ in the $l$-th replication and $N$ is the number of
replications, equal to 500.

Table S1 presents the bias of the parameter estimates for the $SNP_{1}$ and
$SNP_{0}$ models under scenario $C$, $p=10$, $n=5000$.

Table 11: Bias of the parameter estimates for the $SNP_{1}$ and $SNP_{0}$ models under scenario $C$, $p=10$, $n=5000$ $\theta$ |  |  $SNP_{1}$(quasi-ML) |  $SNP_{0}$(full-ML) |  $SNP_{0}$(PL)
---|---|---|---|---
$\alpha_{0_{1}}$ |  | 0.13 | 0.45 | 0.71
$\alpha_{0_{2}}$ |  | 0.06 | 0.16 | 0.30
$\alpha_{0_{3}}$ |  | 0.05 | 0.16 | 0.24
$\alpha_{0_{4}}$ |  | 0.21 | 0.83 | 1.52
$\alpha_{0_{5}}$ |  | 0.06 | 0.18 | 0.28
$\alpha_{0_{6}}$ |  | 0.15 | 0.54 | 1.08
$\alpha_{0_{7}}$ |  | 0.06 | 0.19 | 0.28
$\alpha_{0_{8}}$ |  | 0.15 | 0.55 | 0.87
$\alpha_{0_{9}}$ |  | 0.04 | 0.10 | 0.18
$\alpha_{0_{10}}$ |  | 0.03 | 0.04 | 0.06
$\alpha_{1_{1}}$ |  | 0.15 | 0.52 | 0.74
$\alpha_{1_{2}}$ |  | 0.08 | 0.28 | 0.43
$\alpha_{1_{3}}$ |  | 0.07 | 0.15 | 0.20
$\alpha_{1_{4}}$ |  | 0.24 | 0.96 | 1.57
$\alpha_{1_{5}}$ |  | 0.08 | 0.25 | 0.36
$\alpha_{1_{6}}$ |  | 0.18 | 0.70 | 1.21
$\alpha_{1_{7}}$ |  | 0.08 | 0.20 | 0.28
$\alpha_{1_{8}}$ |  | 0.17 | 0.61 | 0.88
$\alpha_{1_{9}}$ |  | 0.06 | 0.19 | 0.30
$\alpha_{1_{10}}$ |  | 0.04 | 0.07 | 0.13

Table S2 presents the bias of the parameter estimates for the $SNP_{2}$,
$SNP_{1}$ and $SNP_{0}$ models under scenario $E$, $p=10$, $n=1000$.

Table 12: Bias of the parameter estimates for the $SNP_{2}$, $SNP_{1}$ and $SNP_{0}$ models under scenario $E$, $p=10$, $n=1000$ $\theta$ |  $SNP_{2}$(quasi-ML) |  $SNP_{1}$(quasi-ML) |  $SNP_{0}$(full-ML) |  $SNP_{0}$(PL)
---|---|---|---|---
$\alpha_{0_{1}}$ | 0.09 | 0.11 | 0.14 | 0.20
$\alpha_{0_{2}}$ | 0.06 | 0.06 | 0.06 | 0.07
$\alpha_{0_{3}}$ | 0.08 | 0.10 | 0.11 | 0.14
$\alpha_{0_{4}}$ | 0.09 | 0.14 | 0.16 | 0.23
$\alpha_{0_{5}}$ | 0.06 | 0.08 | 0.07 | 0.09
$\alpha_{0_{6}}$ | 0.07 | 0.09 | 0.09 | 0.12
$\alpha_{0_{7}}$ | 0.08 | 0.09 | 0.10 | 0.14
$\alpha_{0_{8}}$ | 0.09 | 0.13 | 0.15 | 0.23
$\alpha_{0_{9}}$ | 0.05 | 0.06 | 0.05 | 0.06
$\alpha_{0_{10}}$ | 0.05 | 0.05 | 0.05 | 0.05
$\alpha_{1_{1}}$ | 0.13 | 0.19 | 0.17 | 0.25
$\alpha_{1_{2}}$ | 0.08 | 0.09 | 0.08 | 0.07
$\alpha_{1_{3}}$ | 0.10 | 0.13 | 0.15 | 0.21
$\alpha_{1_{4}}$ | 0.15 | 0.23 | 0.18 | 0.24
$\alpha_{1_{5}}$ | 0.08 | 0.10 | 0.09 | 0.11
$\alpha_{1_{6}}$ | 0.12 | 0.13 | 0.11 | 0.11
$\alpha_{1_{7}}$ | 0.10 | 0.12 | 0.14 | 0.19
$\alpha_{1_{8}}$ | 0.14 | 0.23 | 0.20 | 0.29
$\alpha_{1_{9}}$ | 0.08 | 0.08 | 0.08 | 0.08
$\alpha_{1_{10}}$ | 0.07 | 0.07 | 0.07 | 0.07

##  Appendix C Graphs of SNP densities

In this section, we present some graphs that include the true non-normal
latent variable densities used for the simulations, as well as the SNP
densities, for scenario C and E. The $\varphi_{1}$ and $\varphi_{2}$
parameters in the SNP densities correspond to the median parameter estimates
obtained across simulations.

Figure 3 displays the true density of the latent variable and the estimated
SNP density with $L=1$ for scenario C, $p=10$ and $n=5000$.

Figure 3: True density and estimated SNP density for scenario C

Figure 4 illustrates the true density of the latent variable and the estimated
SNP densities with $L=1$ and $L=2$ for scenario E, $p=10$ and $n=1000$.

Figure 4: True density and estimated SNP densities for scenario E


then $j\not\in\operatorname{BR}_{\delta}(\boldsymbol{x},u_{f})$. That is,
given $j\not\in\operatorname{BR}_{\delta+2\epsilon}(x,\widetilde{u}_{f})$, we
know there exists $j^{\prime}$ such that
$\widetilde{u}_{f}(\boldsymbol{x},j)-\widetilde{u}_{f}(\boldsymbol{x},j^{\prime})\leq-\delta-2\epsilon$.
Reusing the fact that for any $||\widetilde{u}-u||_{\infty}\leq\epsilon$,
$\boldsymbol{x}\in\Delta^{m}$,
$|\widetilde{u}(\boldsymbol{x},j)-u(\boldsymbol{x},j)|\leq\epsilon$, we have
${u}_{f}(\boldsymbol{x},j)-{u}_{f}(\boldsymbol{x},j^{\prime})\leq\widetilde{u}_{f}(\boldsymbol{x},j)-\widetilde{u}_{f}(\boldsymbol{x},j^{\prime})+2\epsilon\leq-\delta$
and thus $j\not\in\operatorname{BR}_{\delta}(\boldsymbol{x},u_{f})$.

Since
$\operatorname{BR}_{\delta+2\epsilon}(\boldsymbol{x},\widetilde{u}_{f})\supseteq\operatorname{BR}_{\delta}(\boldsymbol{x},u_{f})$,
we have by definition,

$V(\boldsymbol{x};u_{l},u_{f},\delta)=\mathop{min}_{j\in\operatorname{BR}_{\delta}(\boldsymbol{x},{u}_{f})}u_{l}(\boldsymbol{x},j)\geq\mathop{min}_{j\in\operatorname{BR}_{\delta+2\epsilon}(\boldsymbol{x},\widetilde{u}_{f})}u_{l}(\boldsymbol{x},j)=V(\boldsymbol{x};u_{l},\widetilde{u}_{f},\delta+2\epsilon).$

This leads to the following inequalities,

$V^{*}(u_{l},\widetilde{u}_{f},\delta+2\epsilon)=V(\boldsymbol{x}_{1};u_{l},\widetilde{u}_{f},\delta+2\epsilon)\leq
V(\boldsymbol{x}_{1};u_{l},u_{f},\delta)\leq
V(\boldsymbol{x}_{2};u_{l},u_{f},\delta)=V^{*}(u_{l},u_{f},\delta),$

where $\boldsymbol{x}_{1}$ is $(\delta+2\epsilon)$-RSE strategy of Stackelberg
game $u_{l},\widetilde{u}_{f}$, and $\boldsymbol{x}_{2}$ is $\delta$-RSE
strategy of Stackelberg game $u_{l},u_{f}$.

###  C.2 Proof of Proposition 5

Proof. We divide the proof into two parts: $\delta<\Delta$ or
$\delta\geq\Delta$. The proofs for two parts are similar but require different
instance construction.

$\Omega(1/\sqrt{T})$ lower bound instance when $\delta<\Delta$:

Table 7 illustrates a game in which the inducibility gap is $\Delta$. Suppose
$\Delta>\delta$. Consider the following two Stackelberg games $G_{1}$, $G_{2}$
where $r_{l}(i,j)\sim\text{Bern}(u_{l}(i,j))$ and
$r_{f}(i,j)\sim\text{Bern}(u_{f}(i,j))$, with mean values shown in Table 7 and
$\epsilon\in[0,1]$ is a parameter to be determined.

$u_{l},u_{f}$ | $j_{1}$ | $j_{2}$
---|---|---
$i_{1}$ | $1,\frac{1+\epsilon}{2}+\delta$ | $0,\frac{1-\epsilon}{2}$
$i_{2}$ | $0,0$ | $0,\Delta$
$i_{3}$ | $0,\Delta$ | $0,0$

$u_{l},u_{f}$ | $j_{1}$ | $j_{2}$
---|---|---
$i_{1}$ | $1,\frac{1-\epsilon}{2}+\delta$ | $0,\frac{1+\epsilon}{2}$
$i_{2}$ | $0,0$ | $0,\Delta$
$i_{3}$ | $0,\Delta$ | $0,0$

Table 7: Example instances $G_{1}=\\{u_{l},u_{f}\\}$ (left) and
$G_{2}=\\{u_{l},u_{f}\\}$ (right). Each table represents a Stackelberg game
$u_{l},u_{f}\in\mathbb{R}^{3\times 2}$ with inducibility gap $\Delta$.

In Stackelberg game $G_{1}$, the leader strategy in $\delta$-RSE can be
computed as

$\boldsymbol{x}^{*}_{1}=\mathop{argmax}_{\boldsymbol{x}\in\Delta^{m}}\mathop{min}_{j\in\operatorname{BR}_{\delta}(\boldsymbol{x})}u_{l}(\boldsymbol{x},j)=(1,0,0).$

That is, in the $\delta$-RSE of Stackelberg game $G_{1}$, the leader plays
pure strategy $i_{1}$ and $\operatorname{BR}_{\delta}(i_{1})=\\{j_{1}\\}$. As
a result, we have $u_{l}^{*}(\delta)=1$ for $G_{1}$.

On the other hand, for Stackelberg game $G_{2}$, note that
$u_{f}(i_{1},j_{1})-u_{f}(i_{1},j_{2})=\delta-\epsilon<\delta$. Thus, $j_{2}$
is included in the follower’s $\delta$-optimal response set if the leader
plays pure strategy $i_{1}$, i.e. the $\boldsymbol{x}^{*}_{1}$ for the other
game $G_{1}$. Consequently, the leader receives $0$ utility if the leader
plays pure strategy $i_{1}$ in the Stackelberg game $G_{2}$. However, the
leader can play a mixed strategy $\boldsymbol{x}^{*}_{2}$ that includes both
pure actions $i_{1}$ and $i_{3}$ to increase the follower’s utility difference
between responding with $j_{1}$ and $j_{2}$. When the probability of playing
$i_{3}$ is high enough, $\boldsymbol{x}^{*}_{2}$ can make
$u_{f}(\boldsymbol{x}^{*}_{2},j_{1})-u_{f}(\boldsymbol{x}^{*}_{2},j_{2})=\delta$.
Specifically, let $\boldsymbol{x}^{*}_{2}=(p,0,1-p)$, then we have the
following constraint on $p$ in order to exclude $j_{2}$ from the
$\delta$-optimal response set:

$p\cdot\left(\frac{1-\epsilon}{2}+\delta\right)+(1-p)\cdot\Delta-p\cdot\frac{1+\epsilon}{2}=\delta$
(23)

which makes $p=\frac{\Delta-\delta}{\Delta-\delta-\epsilon}$. As a result, for
the Stackelberg game $G_{2}$, we have
$\operatorname{BR}_{\delta}(\boldsymbol{x}^{*}_{2})=\\{j_{1}\\}$ and
$u_{\textup{RSE}}(\delta)=u_{l}(\boldsymbol{x}^{*}_{2},j_{1})=\frac{\Delta-\delta}{\Delta-\delta-\epsilon}$.

Therefore, if the learning algorithm cannot distinguish the above two
Stackelberg games with finite samples and makes a wrong estimation, the
outputted leader strategy will suffer an error of at least
$\frac{\epsilon}{\Delta-\delta+\epsilon}$. Specifically, the following
incorrect estimations can happen:

  1. 1.

If the true game is $G_{1}$, but the learning algorithm estimated the game as
$G_{2}$ and outputs leader strategy as
$\widehat{\boldsymbol{x}}=\boldsymbol{x}^{*}_{2}$. Then the leader utility of
playing $\boldsymbol{x}^{*}_{2}$ will be
$\frac{\Delta-\delta}{\Delta-\delta+\epsilon}$ while
$u_{\textup{RSE}}(\delta)=1$, this has the learning error of
$\frac{\epsilon}{\Delta-\delta+\epsilon}$.

  2. 2.

If the true game is $G_{2}$, but the learning algorithm estimated the game as
$G_{1}$ and outputs leader strategy as
$\widehat{\boldsymbol{x}}=\boldsymbol{x}^{*}_{1}$. Then the leader utility of
playing $\boldsymbol{x}^{*}_{1}$ will be $0$ while
$u_{\textup{RSE}}(\delta)=\frac{\Delta-\delta}{\Delta-\delta+\epsilon}$, and
this has the learning error of $\frac{\Delta-\delta}{\Delta-\delta+\epsilon}$
which is even larger than $\frac{\epsilon}{\Delta-\delta+\epsilon}$.

Consequently, if the learning algorithm can output a leader strategy
$\widehat{\boldsymbol{x}}$ that violates equation (16), i.e.
$\mathop{min}_{j\in\operatorname{BR}_{\delta}(\widehat{\boldsymbol{x}})}u_{l}(\widehat{\boldsymbol{x}},j)>u_{\textup{RSE}}(\delta)-\frac{\epsilon}{\Delta-\delta+\epsilon}$,
then for the game $G_{1}$ we have
$u_{l}(\widehat{\boldsymbol{x}},j)>\frac{\Delta-\delta}{\Delta-\delta+\epsilon}$,
while for the game $G_{2}$ we have
$u_{\textup{RSE}}(\delta)=\frac{\Delta-\delta}{\Delta-\delta+\epsilon}$. In
other words, the learning algorithm can distinguish $G_{1}$ from $G_{2}$ with
$T$ samples.

Therefore, consider the learning algorithm that samples $r_{f}(i,j)$ for $T$
times, and the goal is to identify if $u_{f}\in G_{1}$ or $u_{f}\in G_{2}$. We
prove the proposition by contradiction. Suppose Proposition 5 is not correct,
then with $T$ samples, we can identify if $u_{f}\in G_{1}$ or $u_{f}\in G_{2}$
with probability more than $\frac{2}{3}$. Since $u_{f}(i_{2},j_{1})$,
$u_{f}(i_{2},j_{2})$, $u_{f}(i_{3},j_{1})$, and $u_{f}(i_{3},j_{2})$ are the
same for $G_{1}$ and $G_{2}$, the algorithm can only identify $G_{1}$ and
$G_{2}$ by sampling $u_{f}(i_{1},j_{1})$ or $u_{f}(i_{1},j_{2})$. Formally, we
define the problem as follows. Let $\Omega=[0,1]^{T}$ to be the sample space
for the outcome of $T$ samples of $r_{f}(i_{1},j_{1})$, our goal is to have
the following decision rule

$\text{Rule: }\Omega\rightarrow\\{G_{1},G_{2}\\},$ (24)

which satisfies the following two properties for any $\omega\in\Omega$:

$\begin{split}\mathop{\mathbf{Pr}}[u_{f}\in
G_{1}|\text{Rule}(\omega)=G_{1}]>\frac{2}{3}\text{ and
}\mathop{\mathbf{Pr}}[u_{f}\in
G_{2}|\text{Rule}(\omega)=G_{2}]>\frac{2}{3}.\end{split}$ (25)

As a result, let $\omega_{o}\in\Omega$ be the event this Rule returns $G_{1}$
(i.e., $\text{Rule}(\omega_{o})=G_{1}$). Then we have:

$\mathop{\mathbf{Pr}}[u_{f}\in G_{1}|\omega_{o}]-\mathop{\mathbf{Pr}}[u_{f}\in
G_{2}|\omega_{o}]>\frac{1}{3}$ (26)

Next, for any event $\omega\in\Omega$, let
$P_{k}(\omega)=\mathop{\mathbf{Pr}}[u_{f}\in G_{k}|\omega]$ where $k=1,2$.
Then we have $P_{k}=P_{k,1}\times\cdots\times P_{k,T}$ where $P_{k,t}$ is the
distribution of $t$’th sample of $r_{f}(i_{1},j_{1})$. As a result, by
applying a basic KL-divergence argument to distributions $P_{1}$ and $P_{2}$
[41, Lemma 2.5], for any event $\omega$ we have
$|P_{1}(\omega)-P_{2}(\omega)|\leq\epsilon\sqrt{T}$. Plugging
$\omega=\omega_{0}$ and $\epsilon=\frac{1}{3\sqrt{T}}$ we have
$|P_{1}(\omega)-P_{2}(\omega)|\leq\frac{1}{3}$, contradicting with equation
(26). Therefore, the Stackelberg game instances $G_{1}$ and $G_{2}$ in table 7
with $\epsilon=\frac{1}{3\sqrt{T}}$ proves the proposition when
$\delta<\Delta$.

$\Omega(1)$ lower bound instance when $\delta\geq\Delta$:

Consider the following two Stackelberg games $G_{1}$ and $G_{2}$, where
$r_{l}(i,j)\sim\text{Bern}(u_{l}(i,j))$ and
$r_{f}(i,j)\sim\text{Bern}(u_{f}(i,j))$, with mean values shown in Table 8 and
$\epsilon\in[0,1]$ is a parameter to be determined. Note that in $G_{1}$
follower action $j_{2}$ is dominated by action $j_{1}$. Thus, we have
$\Delta=0$.

$u_{l},u_{f}$ | $j_{1}$ | $j_{2}$
---|---|---
$i_{1}$ | $1,\frac{1+\epsilon}{2}+\delta$ | $0,\frac{1-\epsilon}{2}$
$i_{2}$ | $\frac{1}{2},1$ | $\frac{1}{2},1$

$u_{l},u_{f}$ | $j_{1}$ | $j_{2}$
---|---|---
$i_{1}$ | $1,\frac{1-\epsilon}{2}+\delta$ | $0,\frac{1+\epsilon}{2}$
$i_{2}$ | $\frac{1}{2},1$ | $\frac{1}{2},1$

Table 8: Example instances $G_{1}=\\{u_{l},u_{f}\\}$ (left) and
$G_{2}=\\{u_{l},u_{f}\\}$ (right). Each table represents a Stackelberg game
$u_{l},u_{f}\in\mathbb{R}^{2\times 2}$.

According to Table 8, it is straightforward that the $\delta$-RSE of
Stackelberg game $G_{1}$ is $\langle i_{1},j_{1}\rangle$ where the leader
plays pure strategy $i_{1}$ and follower responds with pure strategy $j_{1}$.
Leader’s utility at $\delta$-RSE is $u_{\textup{RSE}}(\delta)=1$ for $G_{1}$.
On the other hand, for Stackelberg game $G_{2}$, we have
$u_{f}(\boldsymbol{x},j_{1})-u_{f}(\boldsymbol{x},j_{2})<\delta$ for any
leader strategy $\boldsymbol{x}$. Thus, $j_{2}$ is always in the follower’s
$\delta$-optimal response set. To achieve the best utility, the leader plays
pure strategy $i_{2}$ while the follower is indifferent between playing
$j_{1}$ or $j_{2}$. As a result, $u_{\textup{RSE}}(\delta)=\frac{1}{2}$ for
$G_{2}$.

Suppose the learning algorithm can output leader strategy
$\widehat{\boldsymbol{x}}$ that violates equation (16), i.e.,

$\mathop{min}_{j\in\operatorname{BR}_{\delta}(\widehat{\boldsymbol{x}})}u_{l}(\widehat{\boldsymbol{x}},j)>u_{\textup{RSE}}(\delta)-1/2.$

Then for game $G_{1}$, we have
$\mathop{min}_{j\in\operatorname{BR}_{\delta}(\widehat{\boldsymbol{x}})}u_{l}(\widehat{\boldsymbol{x}},j)>\frac{1}{2}$
while
$u_{\textup{RSE}}(\delta)=\mathop{max}_{\boldsymbol{x}\in\Delta^{m}}\mathop{min}_{j\in\operatorname{BR}_{\delta}(\boldsymbol{x})}u_{l}(\boldsymbol{x},j)=\frac{1}{2}$
for game $G_{2}$. Then this means the learning algorithm can identify $G_{1}$
from $G_{2}$.

Therefore, consider the learning algorithm that samples $r_{f}(i,j)$ for $T$
times, and the goal is to identify if $u_{f}\in G_{1}$ or $u_{f}\in G_{2}$. We
prove the proposition by contradiction. Suppose Proposition 5 is not correct,
then with $T$ samples, we can identify if $u_{f}\in G_{1}$ or $u_{f}\in G_{2}$
with probability more than $\frac{2}{3}$. Since $u_{f}(i_{2},j_{1})$ and
$u_{f}(i_{2},j_{2})$ are the same for $G_{1}$ and $G_{2}$, the algorithm can
only identify $G_{1}$ and $G_{2}$ by sampling $u_{f}(i_{1},j_{1})$ or
$u_{f}(i_{1},j_{2})$. Formally, we define the problem as follows. Let
$\Omega=[0,1]^{T}$ to be the sample space for the outcome of $T$ samples of
$r_{f}(i_{1},j_{1})$ and our goal is to have the following decision rule

$\text{Rule: }\Omega\rightarrow\\{G_{1},G_{2}\\}$ (27)

which satisfies the following two properties:

$\begin{split}\mathop{\mathbf{Pr}}[u_{f}\in
G_{1}|\text{Rule(observations)}=G_{1}]>\frac{2}{3}\\\
\mathop{\mathbf{Pr}}[u_{f}\in
G_{2}|\text{Rule(observations)}=G_{2}]>\frac{2}{3}\end{split}$ (28)

Therefore, let $\omega_{o}\in\Omega$ be the event this Rule returns $G_{1}$.
Then we have:

$\mathop{\mathbf{Pr}}[u_{f}\in G_{1}|\omega_{o}]-\mathop{\mathbf{Pr}}[u_{f}\in
G_{2}|\omega_{o}]>\frac{1}{3}$ (29)

Next, for any event $\omega\in\Omega$, let
$P_{k}(\omega)=\mathop{\mathbf{Pr}}[u_{f}\in G_{k}|\omega]$ where $k=1,2$.
Then we have $P_{k}=P_{k,1}\times\cdots\times P_{k,T}$ where $P_{k,t}$ is the
distribution of $t$’th sample of $r_{f}(i_{1},j_{1})$. As a result, by
applying a basic KL-divergence argument to distributions $P_{1}$ and $P_{2}$
[41, Lemma 2.5], for any event $\omega$ we have
$|P_{1}(\omega)-P_{2}(\omega)|\leq\epsilon\sqrt{T}$. Plugging
$\omega=\omega_{0}$ and $\epsilon=\frac{1}{3\sqrt{T}}$ we have
$|P_{1}(\omega)-P_{2}(\omega)|\leq\frac{1}{3}$, contradicting with equation
(29).

Therefore, choosing the problem class to be $G_{1}$ and $G_{2}$ in table 8
with $\epsilon=\frac{1}{3\sqrt{T}}$, finishes the proof.

∎